cut on identical particles
I'm using batch mode to analyse some process with two photons p,p -> A,A in the final state. But the cut
Cut parameter: Z(A)
Cut invert: False
Cut min: 55
Cut max:
seems to only apply for one of the two photons.
I checked this with usrfun.c. How to define a cut requires both photons should have transverse energy larger than 100GeV for example? Thanks.
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- CalcHEP Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Shao-Feng Ge
- Solved:
- Last query:
- Last reply:
Revision history for this message
|
#1 |
Z(A) is applied for each A. You can check it using GUI interface
Best
Alexander Pukhov
On 01/18/2018 10:03 AM, Shao-Feng Ge wrote:
> Question #663205 on CalcHEP changed:
> https:/
>
> Description changed to:
> I'm using batch mode to analyse some process with two photons p,p -> A,A
> in the final state. But the cut
>
> Cut parameter: Z(A)
> Cut invert: False
> Cut min: 55
> Cut max:
>
> seems to only apply for one of the two photons.
>
> I checked this with usrfun.c. How to define a cut requires both photons
> should have transverse energy larger than 100GeV for example? Thanks.
>
Revision history for this message
|
#2 |
I'm trying to use batch mode to verify this explicitly and encounter extra problems. First, I define the usrfun.c as following:
=======
double usrfun(char * name, int nIn, int nOut, double * pvect,char*
{ char p1[10],p2[10]; // for 2 particles in MT(p1,p2)
double tValue = 1.;
double invMass;
if(
{
// momenta
double *q1=pvect+4*2, *q2=pvect+4*3;
/*
double m1 = q1[0]*q1[0] - q1[1]*q1[1] - q1[2]*q1[2] - q1[3]*q1[3];
double m2 = q2[0]*q2[0] - q2[1]*q2[1] - q2[2]*q2[2] - q2[3]*q2[3];
printf("q1 = {%f, %f, %f, %f} [m2 = %f]\n q2 = {%f, %f, %f, %f} [m2 = %f]\n", q1[0], q1[1], q1[2], q1[3], m1, q1[4], q1[5], q1[6], q1[7], m2);
*/
// invariant mass
invMass = sqrt( pow(q1[0] + q2[0], 2.)
- pow(q1[1] + q2[1], 2.)
- pow(q1[2] + q2[2], 2.)
- pow(q1[3] + q2[3], 2.) );
// total momentum [ref momentum]
double qt[4];
for (int i = 0; i < 4; i++)
qt[i] = q1[i] + q2[i];
// gamma
double gamma = qt[0] / invMass;
// beta vec3
double betaX = - qt[1] / qt[0];
double betaY = - qt[2] / qt[0];
double betaZ = - qt[3] / qt[0];
// beta
double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
double beta = sqrt(beta2);
// boosted momentum
double Q1[8]; double *Q2 = Q1+4;
for (int i = 0; i < 2; i++)
{
// parallel norm
double para = (betaX*q1[1+4*i] + betaY*q1[2+4*i] + betaZ*q1[3+4*i]) / beta;
// parallel vec3
double pPara[3] = {para*betaX/beta, para*betaY/beta, para*betaZ/beta};
// original energy
double pE = q1[0+4*i];
// boosted energy
double bE = gamma*(pE + beta*para);
// boosted parallel norm
para = gamma*(para + beta*pE);
// boosted vector
Q1[0+4*i] = bE;
Q1[1+4*i] = q1[1+4*i] + para*betaX/beta - pPara[0];
Q1[2+4*i] = q1[2+4*i] + para*betaY/beta - pPara[1];
Q1[3+4*i] = q1[3+4*i] + para*betaZ/beta - pPara[2];
}
/*
double M1 = Q1[0]*Q1[0] - Q1[1]*Q1[1] - Q1[2]*Q1[2] - Q1[3]*Q1[3];
double M2 = Q2[0]*Q2[0] - Q2[1]*Q2[1] - Q2[2]*Q2[2] - Q2[3]*Q2[3];
printf("Q1 = {%f, %f, %f, %f} [M1 = %f]\n Q2 = {%f, %f, %f, %f} [M2 = %f]\n\n", Q1[0], Q1[1], Q1[2], Q1[3], M1, Q2[0], Q2[1], Q2[2], Q2[3], M2);
*/
double cosTheta = Q1[3] / Q1[0];
if (abs(cosTheta) > 0.9)
tValue = 0.;
}
if(
{
// momenta
double *qa=pvect+4*2;
double pE = qa[0];
double pZ = qa[3];
double cosTheta = pZ/pE;
tValue = - log((1 - cosTheta) / (1 + cosTheta)) / 2.;
}
if(
{
// momenta
double *qa=pvect+4*3;
double pE = qa[0];
double pZ = qa[3];
double cosTheta = pZ/pE;
tValue = - log((1 - cosTheta) / (1 + cosTheta)) / 2.;
}
return tValue;
}
=======
Then, I use batch script
=======
Model: Standard Model
Model changed: False
Gauge: Feynman
Composite: p=u,U,d,
Process: p,p->A,A
pdf1: cteq6l1 (proton)
pdf2: cteq6l1 (proton)
p1: 6500
p2: 6500
Cut parameter: N(A)
Cut invert: False
Cut min:-2.37
Cut max: 2.37
Cut parameter: N(A)
Cut invert: True
Cut min: 1.37
Cut max: 1.52
Cut parameter: N(A)
Cut invert: True
Cut min:-1.52
Cut max:-1.37
Cut parameter: UEtM(A,A)
Cut invert: False
Cut min: 0.5
Cut max: 1.5
nSess_1 : 5
nCalls_1 : 100000
nSess_2 : 5
nCalls_2 : 1000000
Dist parameter: UNa(A,A)
Dist min: -2.4
Dist max: 2.4
Dist n bins: 24
Dist title: UNa(A)
Dist x-title: Na(A)
Dist parameter: UNb(A,A)
Dist min: -2.4
Dist max: 2.4
Dist n bins: 24
Dist title: UNb(A)
Dist x-title: Nb(A)
Dist parameter: M(A,A)
Dist min: 200
Dist max: 2700
Dist n bins: 125
Dist title: M(A,A)
Dist x-title: M(AA) [GeV]
Filename: SM
=======
running batch mode would return:
=======
Error in calculation of cross section for u,U->A,A
Please have a look at cs.o and cs.out in Processes/
Killing all jobs and quitting.
=======
As I tried, the code can successfully run after deleting the UNa(A,A) & UNb(A,A) distributions in the batch script. Any idea?
Revision history for this message
|
#3 |
You need
#include <math.h>
abs => fabs
UEtM return 0 or 1. Thus cuts which you apply to UEtM will not work as
you expect.
In GUI mode it works somehow. But we need cut for small T(A) to avoid
infinities caused zero mass of quark.
Best
Alexander Pukhov
PS: I think it is a bad idea to check implementation of new functions
in batch mode.
On 01/19/2018 02:03 AM, Shao-Feng Ge wrote:
> Question #663205 on CalcHEP changed:
> https:/
>
> Status: Answered => Open
>
> Shao-Feng Ge is still having a problem:
> I'm trying to use batch mode to verify this explicitly and encounter
> extra problems. First, I define the usrfun.c as following:
>
> =======
> double usrfun(char * name, int nIn, int nOut, double * pvect,char*
> { char p1[10],p2[10]; // for 2 particles in MT(p1,p2)
> double tValue = 1.;
> double invMass;
>
> if(name=
> {
> // momenta
> double *q1=pvect+4*2, *q2=pvect+4*3;
>
> /*
> double m1 = q1[0]*q1[0] - q1[1]*q1[1] - q1[2]*q1[2] - q1[3]*q1[3];
> double m2 = q2[0]*q2[0] - q2[1]*q2[1] - q2[2]*q2[2] - q2[3]*q2[3];
>
> printf("q1 = {%f, %f, %f, %f} [m2 = %f]\n q2 = {%f, %f, %f, %f} [m2 = %f]\n", q1[0], q1[1], q1[2], q1[3], m1, q1[4], q1[5], q1[6], q1[7], m2);
> */
>
> // invariant mass
> invMass = sqrt( pow(q1[0] + q2[0], 2.)
> - pow(q1[1] + q2[1], 2.)
> - pow(q1[2] + q2[2], 2.)
> - pow(q1[3] + q2[3], 2.) );
>
> // total momentum [ref momentum]
> double qt[4];
> for (int i = 0; i < 4; i++)
> qt[i] = q1[i] + q2[i];
>
> // gamma
> double gamma = qt[0] / invMass;
> // beta vec3
> double betaX = - qt[1] / qt[0];
> double betaY = - qt[2] / qt[0];
> double betaZ = - qt[3] / qt[0];
> // beta
> double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
> double beta = sqrt(beta2);
>
> // boosted momentum
> double Q1[8]; double *Q2 = Q1+4;
>
> for (int i = 0; i < 2; i++)
> {
> // parallel norm
> double para = (betaX*q1[1+4*i] + betaY*q1[2+4*i] + betaZ*q1[3+4*i]) / beta;
> // parallel vec3
> double pPara[3] = {para*betaX/beta, para*betaY/beta, para*betaZ/beta};
> // original energy
> double pE = q1[0+4*i];
> // boosted energy
> double bE = gamma*(pE + beta*para);
> // boosted parallel norm
> para = gamma*(para + beta*pE);
> // boosted vector
> Q1[0+4*i] = bE;
> Q1[1+4*i] = q1[1+4*i] + para*betaX/beta - pPara[0];
> Q1[2+4*i] = q1[2+4*i] + para*betaY/beta - pPara[1];
> Q1[3+4*i] = q1[3+4*i] + para*betaZ/beta - pPara[2];
> }
>
> /*
> double M1 = Q1[0]*Q1[0] - Q1[1]*Q1[1] - Q1[2]*Q1[2] - Q1[3]*Q1[3];
> double M2 = Q2[0]*Q2[0] - Q2[1]*Q2[1] - Q2[2]*Q2[2] - Q2[3]*Q2[3];
>
> printf("Q1 = {%f, %f, %f, %f} [M1 = %f]\n Q2 = {%f, %f, %f, %f} [M2 = %f]\n\n", Q1[0], Q1[1], Q1[2], Q1[3], M1, Q2[0], Q2[1], Q2[2], Q2[3], M2);
> */
>
> double cosTheta = Q1[3] / Q1[0];
>
> if (abs(cosTheta) > 0.9)
> tValue = 0.;
> }
>
> if(name=
> {
> // momenta
> double *qa=pvect+4*2;
> double pE = qa[0];
> double pZ = qa[3];
>
> double cosTheta = pZ/pE;
> tValue = - log((1 - cosTheta) / (1 + cosTheta)) / 2.;
>
> }
>
> if(name=
> {
> // momenta
> double *qa=pvect+4*3;
> double pE = qa[0];
> double pZ = qa[3];
>
> double cosTheta = pZ/pE;
> tValue = - log((1 - cosTheta) / (1 + cosTheta)) / 2.;
> }
>
> return tValue;
> }
>
> =======
>
> Then, I use batch script
>
> =======
>
> Model: Standard Model
> Model changed: False
> Gauge: Feynman
> Composite: p=u,U,d,
> Process: p,p->A,A
>
> pdf1: cteq6l1 (proton)
> pdf2: cteq6l1 (proton)
>
> p1: 6500
> p2: 6500
>
> Cut parameter: N(A)
> Cut invert: False
> Cut min:-2.37
> Cut max: 2.37
>
> Cut parameter: N(A)
> Cut invert: True
> Cut min: 1.37
> Cut max: 1.52
>
> Cut parameter: N(A)
> Cut invert: True
> Cut min:-1.52
> Cut max:-1.37
>
> Cut parameter: UEtM(A,A)
> Cut invert: False
> Cut min: 0.5
> Cut max: 1.5
>
> nSess_1 : 5
> nCalls_1 : 100000
> nSess_2 : 5
> nCalls_2 : 1000000
>
> Dist parameter: UNa(A,A)
> Dist min: -2.4
> Dist max: 2.4
> Dist n bins: 24
> Dist title: UNa(A)
> Dist x-title: Na(A)
>
> Dist parameter: UNb(A,A)
> Dist min: -2.4
> Dist max: 2.4
> Dist n bins: 24
> Dist title: UNb(A)
> Dist x-title: Nb(A)
>
> Dist parameter: M(A,A)
> Dist min: 200
> Dist max: 2700
> Dist n bins: 125
> Dist title: M(A,A)
> Dist x-title: M(AA) [GeV]
>
> Filename: SM
>
> =======
>
> running batch mode would return:
>
> =======
>
> Error in calculation of cross section for u,U->A,A
> Please have a look at cs.o and cs.out in Processes/
> Killing all jobs and quitting.
>
> =======
>
> As I tried, the code can successfully run after deleting the UNa(A,A) &
> UNb(A,A) distributions in the batch script. Any idea?
>
Revision history for this message
|
#4 |
Did it solve your problem?
Revision history for this message
|
#5 |
Sorry for not responding. Yes, it does solve my problem. Thanks a lot.