New class used for primary vertex finding (AliITSVertexerTracks)

[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.cxx

diff --git a/EVGEN/AliGenHIJINGpara.cxx b/EVGEN/AliGenHIJINGpara.cxx
index 565a1ac48dae38de71036e3f110e9304e5a6db5a..17fae20c87a23625da0f2e929e0f270a4420ee1b 100644 (file)
--- a/EVGEN/AliGenHIJINGpara.cxx
+++ b/EVGEN/AliGenHIJINGpara.cxx
@@ -15,6 +15,33 @@
 
 /*
 $Log$
+Revision 1.19  2003/01/14 10:50:18  alibrary
+Cleanup of STEER coding conventions
+
+Revision 1.18  2002/12/11 11:58:11  morsch
+Bug in formula for pi0 energy for decay corrected.
+
+Revision 1.17  2002/12/10 17:44:57  morsch
+Correct mother child relation for pi0.
+
+Revision 1.16  2002/11/28 11:46:15  morsch
+Don't track pi0 if already decayed.
+
+Revision 1.15  2002/11/28 11:38:53  morsch
+Typo corrected.
+
+Revision 1.14  2002/11/26 17:12:36  morsch
+Decay pi0 if requested.
+
+Revision 1.13  2002/10/14 14:55:35  hristov
+Merging the VirtualMC branch to the main development branch (HEAD)
+
+Revision 1.11.4.1  2002/07/24 08:56:28  alibrary
+Updating EVGEN on TVirtulaMC
+
+Revision 1.12  2002/06/19 06:56:34  hristov
+Memory leak corrected
+
 Revision 1.11  2002/03/20 10:21:13  hristov
 Set fPtMax to 15 GeV in order to avoid some numerical problems
 
@@ -78,14 +105,19 @@ All coding rule violations except RS3 corrected (AM)
 //                                                               //
 ///////////////////////////////////////////////////////////////////
 
-#include "AliGenHIJINGpara.h"
+#include <TArrayF.h>
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TF1.h>
+#include <TParticle.h>
+#include <TPDGCode.h>
+
+#include "AliConst.h"
+#include "AliDecayer.h"
+#include "AliDecayerPythia.h"
 #include "AliGenEventHeader.h"
+#include "AliGenHIJINGpara.h"
 #include "AliRun.h"
-#include "AliConst.h"
-#include "AliPDG.h"
-
-#include <TF1.h>
-#include <TArrayF.h>
 
 ClassImp(AliGenHIJINGpara)
 
@@ -197,12 +229,15 @@ AliGenHIJINGpara::AliGenHIJINGpara()
     //
     // Default constructor
     //
-    fPtpi = 0;
-    fPtka = 0;
-    fETApic = 0;
-    fETAkac = 0;
+    fPtpi    =  0;
+    fPtka    =  0;
+    fETApic  =  0;
+    fETAkac  =  0;
+    fDecayer =  0;
+    fNt      = -1;
     SetCutVertexZ();
     SetPtRange();
+    SetPi0Decays();
 }
 
 //_____________________________________________________________________________
@@ -212,14 +247,17 @@ AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
   // 
   // Standard constructor
   //
-    fName="HIGINGpara";
+    fName="HIJINGpara";
     fTitle="HIJING Parametrisation Particle Generator";
-    fPtpi = 0;
-    fPtka = 0;
-    fETApic = 0;
-    fETAkac = 0;
+    fPtpi    =  0;
+    fPtka    =  0;
+    fETApic  =  0;
+    fETAkac  =  0;
+    fDecayer =  0;
+    fNt      = -1;
     SetCutVertexZ();
     SetPtRange();
+    SetPi0Decays();
 }
 
 //_____________________________________________________________________________
@@ -286,8 +324,12 @@ void AliGenHIJINGpara::Init()
        printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");           
        printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax);
     }
+//
+//
+    if (fPi0Decays) fDecayer = new AliDecayerPythia();
 }
 
+
 //_____________________________________________________________________________
 void AliGenHIJINGpara::Generate()
 {
@@ -307,7 +349,7 @@ void AliGenHIJINGpara::Generate()
     Float_t pt, pl, ptot;
     Float_t phi, theta;
     Float_t p[3];
-    Int_t i, part, nt, j;
+    Int_t i, part, j;
     //
     TF1 *ptf;
     TF1 *etaf;
@@ -340,7 +382,7 @@ void AliGenHIJINGpara::Generate()
            if(random[0]<kBorne) {
                part=kPions[Int_t (random[1]*3)];
                ptf=fPtpi;
-             etaf=fETApic;
+               etaf=fETApic;
            } else {
                part=kKaons[Int_t (random[1]*4)];
                ptf=fPtka;
@@ -363,10 +405,23 @@ void AliGenHIJINGpara::Generate()
                        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
                }
            }
-           SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
+           if (part == kPi0 && fPi0Decays){
+//
+//          Decay pi0 if requested
+               SetTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
+               KeepTrack(fNt);
+               DecayPi0(origin, p);
+           } else {
+               SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
+               KeepTrack(fNt);
+           }
+
            break;
        }
+       SetHighWaterMark(fNt);
     }
+//
+
 // Header
     AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
 // Event Vertex
@@ -381,6 +436,41 @@ AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
 }
 
 void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
-  AliGenerator::SetPtRange(ptmin, ptmax);
+    AliGenerator::SetPtRange(ptmin, ptmax);
 }
 
+void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p) 
+{
+//
+//    Decay the pi0
+//    and put decay products on the stack
+//
+    static TClonesArray *particles;
+    if(!particles) particles = new TClonesArray("TParticle",1000);
+//    
+    const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
+    Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]+ kMass * kMass);
+//
+//  Decay the pi0    
+    TLorentzVector pmom(p[0], p[1], p[2], e);
+    fDecayer->Decay(kPi0, &pmom);
+    
+//
+// Put decay particles on the stack
+//
+    Float_t polar[3] = {0., 0., 0.};
+    Int_t np = fDecayer->ImportParticles(particles);
+    Int_t nt;    
+    for (Int_t i = 1; i < np; i++)
+    {
+       TParticle* iParticle =  (TParticle *) particles->At(i);
+       p[0] = iParticle->Px();
+       p[1] = iParticle->Py();
+       p[2] = iParticle->Pz();
+       Int_t part = iParticle->GetPdgCode();
+
+       SetTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
+       KeepTrack(nt);
+    }
+    fNt = nt;
+}