cut on identical particles

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://answers.launchpad.net/calchep/+question/663205
>
> 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.
>

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**pName,int*pCode)
{ char p1[10],p2[10]; // for 2 particles in MT(p1,p2)
    double tValue = 1.;
  double invMass;

    if(name==strstr(name,"EtM(")) // name is started from "EtM("
    {
   // 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==strstr(name,"Na("))
    {
   // 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==strstr(name,"Nb("))
    {
   // 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,D,c,C,s,S,b,B
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/m1/Feynman/p1/single.
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?

You need

           #include <math.h>

           #include<string.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://answers.launchpad.net/calchep/+question/663205
>
> 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**pName,int*pCode)
> { char p1[10],p2[10]; // for 2 particles in MT(p1,p2)
> double tValue = 1.;
> double invMass;
>
> if(name==strstr(name,"EtM(")) // name is started from "EtM("
> {
> // 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==strstr(name,"Na("))
> {
> // 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==strstr(name,"Nb("))
> {
> // 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,D,c,C,s,S,b,B
> 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/m1/Feynman/p1/single.
> 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?
>

Did it solve your problem?

Sorry for not responding. Yes, it does solve my problem. Thanks a lot.