STARLIGHT update:

[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / gammagammaleptonpair.cpp

diff --git a/STARLIGHT/starlight/src/gammagammaleptonpair.cpp b/STARLIGHT/starlight/src/gammagammaleptonpair.cpp
index cdcb458130a36a58b0d8e720bfd2af35045da134..e4a751b17c68ff821711d68a136d79e704decb96 100644 (file)
--- a/STARLIGHT/starlight/src/gammagammaleptonpair.cpp
+++ b/STARLIGHT/starlight/src/gammagammaleptonpair.cpp
@@ -385,17 +385,17 @@ double Gammagammaleptonpair::pp_2(double E)
 
 //______________________________________________________________________________
 void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
-                                    double  ,  // E (unused)
-                                    double  W,
-                                    double  px0, double  py0, double  pz0,
-                                    double& px1, double& py1, double& pz1,
-                                    double& px2, double& py2, double& pz2,
-                                    int&    iFbadevent)
+                                       double  ,  // E (unused)
+                                       double  W,
+                                       double  px0, double  py0, double  pz0,
+                                       double& px1, double& py1, double& pz1, double& E1,
+                                       double& px2, double& py2, double& pz2, double& E2,
+                                       double& mdec,
+                                       int&    iFbadevent)
 {
     //     This routine decays a particle into two particles of mass mdec,
     //     taking spin into account
 
-    double mdec=0.,E1=0.,E2=0.;
     double pmag, anglelep[20001];
     // double ytest=0.,dndtheta;
     double phi,theta,xtest,Ecm;
@@ -536,7 +536,8 @@ starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
     int iFbadevent=0;
     starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
        
-    double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
+    double px2=0.,px1=0.,py2=0.,pt2=0.,py1=0.,pz2=0.,pz1=0.,pt1=0.;
+    double E1=0.,E2=0.,mdec=0.;
 //this function decays particles and writes events to a file
     //zero out the event structure
     leptonpair._numberOfTracks=0;
@@ -554,32 +555,25 @@ starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
     picky(rapidity);
 
     pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
-    twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+    twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px1,px2,py2,E2,mdec,iFbadevent);
     if (iFbadevent==0){
-       int q1=0,q2=0; 
-
-       double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-       if (xtest<0.5)
-       {
-           q1=1;
-           q2=-1;
-       }
-       else{
-           q1=-1;
-           q2=1;
-       }       
+      const bool xtest(randyInstance.Rndom() < 0.5);
+      const bool charges[2] = {
+       ( xtest) ? +1 : -1,
+       (!xtest) ? -1 : +1
+      };
        leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
        leptonpair.px[0]=px1;
        leptonpair.py[0]=py1;
        leptonpair.pz[0]=pz1;
        leptonpair._fsParticle[0]=ipid; 
-       leptonpair._charge[0]=q1;
+       leptonpair._charge[0]=charges[0];
 
        leptonpair.px[1]=px2;
        leptonpair.py[1]=py2;
        leptonpair.pz[1]=pz2;
        leptonpair._fsParticle[1]=ipid;
-       leptonpair._charge[1]=q2;
+       leptonpair._charge[1]=charges[1];
 
        ievent=ievent+1;
     }
@@ -602,7 +596,7 @@ upcEvent Gammagammaleptonpair::produceEvent()
    int iFbadevent=0;
    starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
    
-   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.;
+   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.,mdec=0.;
    bool accepted = false;
    do{ 
      //this function decays particles and writes events to a file
@@ -615,7 +609,7 @@ upcEvent Gammagammaleptonpair::produceEvent()
    
   
      _nmbAttempts++;
-     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px2,py2,pz2,E2,mdec,iFbadevent);
      double pt1chk = sqrt(px1*px1+py1*py1);
      double pt2chk = sqrt(px2*px2+py2*py2);
     
@@ -649,30 +643,18 @@ upcEvent Gammagammaleptonpair::produceEvent()
    //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
    
    if (iFbadevent==0){
-     int q1=0,q2=0; 
-     
-     double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-     if (xtest<0.5)
-       {
-        q1=1;
-        q2=-1;
-       }
-     else{
-       q1=-1;
-       q2=1;
-     }
+     const bool xtest(randyInstance.Rndom() < 0.5);
+     const int charges[2] = {
+       ( xtest) ? +1 : -1,
+       (!xtest) ? -1 : +1
+     };
 
      // The new stuff
-     double mlepton = getMass(); 
-     E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 ); 
-     E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 ); 
-
-     starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1);
+     starlightParticle particle1(px1, py1, pz1, E1, mdec, -charges[0]*ipid, charges[0]);
      event.addParticle(particle1);
      
-     starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2);
-     event.addParticle(particle2);
-     
+     starlightParticle particle2(px2, py2, pz2, E2, mdec, -charges[1]*ipid, charges[1]);
+     event.addParticle(particle2);     
     }
    return event;
 }