Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / EVGEN / AliGenPythia.cxx
index 3336bc4..0820dcb 100644 (file)
 
 /*
 $Log$
+Revision 1.52.4.4  2002/10/10 16:40:08  hristov
+Updating VirtualMC to v3-09-02
+
+Revision 1.62  2002/09/24 10:00:01  morsch
+CheckTrigger() corrected.
+
+Revision 1.61  2002/07/30 09:52:38  morsch
+Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
+
+Revision 1.60  2002/07/19 14:49:03  morsch
+Typo corrected.
+
+Revision 1.59  2002/07/19 14:35:36  morsch
+Count total number of trials. Print mean Q, x1, x2.
+
+Revision 1.58  2002/07/17 10:04:09  morsch
+SetYHard method added.
+
+Revision 1.57  2002/05/22 13:22:53  morsch
+Process kPyMbNonDiffr added.
+
+Revision 1.56  2002/04/26 10:30:01  morsch
+Option kPyBeautyPbMNR added. (N. Carrer).
+
+Revision 1.55  2002/04/17 10:23:56  morsch
+Coding Rule violations corrected.
+
+Revision 1.54  2002/03/28 11:49:10  morsch
+Pass status code in SetTrack.
+
+Revision 1.53  2002/03/25 14:51:13  morsch
+New stack-fill and count options introduced (N. Carrer).
+
+Revision 1.51  2002/03/06 08:46:57  morsch
+- Loop until np-1
+- delete dyn. alloc. arrays (N. Carrer)
+
+Revision 1.50  2002/03/03 13:48:50  morsch
+Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
+NLO calculations (Nicola Carrer).
+
+Revision 1.49  2002/02/08 16:50:50  morsch
+Add name and title in constructor.
+
+Revision 1.48  2001/12/20 11:44:28  morsch
+Add kinematic bias for direct gamma production.
+
+Revision 1.47  2001/12/19 14:45:00  morsch
+Store number of trials in header.
+
+Revision 1.46  2001/12/19 10:36:19  morsch
+Add possibility if jet kinematic biasing.
+
+Revision 1.45  2001/11/28 08:06:52  morsch
+Use fMaxLifeTime parameter.
+
+Revision 1.44  2001/11/27 13:13:07  morsch
+Maximum lifetime for long-lived particles to be put on the stack is parameter.
+It can be set via SetMaximumLifetime(..).
+
+Revision 1.43  2001/10/21 18:35:56  hristov
+Several pointers were set to zero in the default constructors to avoid memory management problems
+
+Revision 1.42  2001/10/15 08:21:55  morsch
+Vertex truncation settings moved to AliGenMC.
+
+Revision 1.41  2001/10/08 08:45:42  morsch
+Possibility of vertex cut added.
+
 Revision 1.40  2001/09/25 11:30:23  morsch
 Pass event vertex to header.
 
@@ -117,6 +186,17 @@ Revision 1.11  1999/09/29 09:24:14  fca
 Introduction of the Copyright and cvs Log
 */
 
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+
 #include "AliGenPythia.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliDecayerPythia.h"
@@ -132,8 +212,14 @@ AliGenPythia::AliGenPythia()
                  :AliGenMC()
 {
 // Default Constructor
-  fDecayer = new AliDecayerPythia();
+  fParticles = 0;
+  fPythia    = 0;
+  fDecayer   = new AliDecayerPythia();
   SetEventListRange();
+  SetJetPhiRange();
+  SetJetEtaRange();
+  SetGammaPhiRange();
+  SetGammaEtaRange();
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
@@ -143,6 +229,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
 // semimuonic decay
 // structure function GRVHO
 //
+    fName = "Pythia";
+    fTitle= "Particle Generator using PYTHIA";
     fXsection  = 0.;
     fNucA1=0;
     fNucA2=0;
@@ -150,20 +238,33 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetStrucFunc();
     SetForceDecay();
     SetPtHard();
+    SetYHard();
     SetEnergyCMS();
     fDecayer = new AliDecayerPythia();
     // Set random number generator 
     sRandom=fRandom;
-    SetEventListRange();
     fFlavorSelect   = 0;
     // Produced particles  
     fParticles = new TClonesArray("TParticle",1000);
     fEventVertex.Set(3);
+    SetEventListRange();
+    SetJetPhiRange();
+    SetJetEtaRange();
+    SetGammaPhiRange();
+    SetGammaEtaRange();
+    // Options determining what to keep in the stack (Heavy flavour generation)
+    fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
+    fFeedDownOpt = kTRUE;             // allow feed down from higher family
+    // Fragmentation on/off
+    fFragmentation = kTRUE;
+    // Default counting mode
+    fCountMode = kCountAll;
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 {
 // copy constructor
+    Pythia.Copy(*this);
 }
 
 AliGenPythia::~AliGenPythia()
@@ -194,8 +295,18 @@ void AliGenPythia::Init()
 
 
     fPythia->SetCKIN(3,fPtHardMin);
-    fPythia->SetCKIN(4,fPtHardMax);    
+    fPythia->SetCKIN(4,fPtHardMax);
+    fPythia->SetCKIN(7,fYHardMin);
+    fPythia->SetCKIN(8,fYHardMax);
+    
     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
+    // Fragmentation?
+    if (fFragmentation) {
+      fPythia->SetMSTP(111,1);
+    } else {
+      fPythia->SetMSTP(111,0);
+    }
+
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 
     //    fPythia->Pylist(0);
@@ -204,20 +315,20 @@ void AliGenPythia::Init()
     switch (fProcess) 
     {
     case kPyCharm:
-       fParentSelect[0] =  411;
-       fParentSelect[1] =  421;
-       fParentSelect[2] =  431;
-       fParentSelect[3] = 4122;        
-       fFlavorSelect    = 4;
-       break;
     case kPyCharmUnforced:
-       fParentSelect[0] =  411;
-       fParentSelect[1] =  421;
-       fParentSelect[2] =  431;
-       fParentSelect[3]=  4122;
-       fFlavorSelect    = 4;   
+    case kPyCharmPbMNR:
+       fParentSelect[0] =   411;
+       fParentSelect[1] =   421;
+       fParentSelect[2] =   431;
+       fParentSelect[3] =  4122;
+       fFlavorSelect    =  4;  
+       break;
+    case kPyD0PbMNR:
+       fParentSelect[0] =   421;
+       fFlavorSelect    =   4; 
        break;
     case kPyBeauty:
+    case kPyBeautyPbMNR:
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;
@@ -242,10 +353,19 @@ void AliGenPythia::Init()
        fParentSelect[0] = 443;
        break;
     case kPyMb:
+    case kPyMbNonDiffr:
     case kPyJets:
     case kPyDirectGamma:
        break;
     }
+//
+//  This counts the total number of calls to Pyevnt() per run.
+    fTrialsRun = 0;
+    fQ         = 0.;
+    fX1        = 0.;
+    fX2        = 0.;    
+    fNev       = 0 ;
+//    
     AliGenMC::Init();
 }
 
@@ -265,7 +385,7 @@ void AliGenPythia::Generate()
     Int_t j, kf;
     fTrials=0;
 
-//  Set collision vertex position 
+    //  Set collision vertex position 
     if(fVertexSmear==kPerEvent) {
        fPythia->SetMSTP(151,1);
        for (j=0;j<3;j++) {
@@ -281,6 +401,7 @@ void AliGenPythia::Generate()
        if (gAlice->GetEvNumber()>=fDebugEventFirst &&
            gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
        fTrials++;
+       
        fPythia->ImportParticles(fParticles,"All");
 
 //
@@ -291,25 +412,31 @@ void AliGenPythia::Generate()
        Int_t np = fParticles->GetEntriesFast();
        if (np == 0 ) continue;
 // Get event vertex and discard the event if the Z coord. is too big   
-       TParticle *  iparticle = (TParticle *) fParticles->At(0);
-       Float_t distz = iparticle->Vz()/10. - fOrigin.At(2);
-       if(TMath::Abs(distz)>fCutVertexZ)continue;
+       TParticle *iparticle = (TParticle *) fParticles->At(0);
+       Float_t distz = iparticle->Vz()/10.;
+       if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 //
-       fEventVertex[0] = iparticle->Vx()/10.;
-       fEventVertex[1] = iparticle->Vy()/10.;  
-       fEventVertex[2] = iparticle->Vz()/10.;
+       fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
+       fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
+       fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
 //
        Int_t* pParent   = new Int_t[np];
        Int_t* pSelected = new Int_t[np];
        Int_t* trackIt   = new Int_t[np];
-       for (i=0; i< np-1; i++) {
+       for (i=0; i< np; i++) {
            pParent[i]   = -1;
            pSelected[i] =  0;
+           trackIt[i]   =  0;
        }
-       printf("\n **************************************************%d\n",np);
-       Int_t nc = 0;
-       if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
-           for (i = 0; i<np-1; i++) {
+
+       Int_t nc = 0;        // Total n. of selected particles
+       Int_t nParents = 0;  // Selected parents
+       Int_t nTkbles = 0;   // Trackable particles
+       if (fProcess != kPyMb && fProcess != kPyJets && 
+           fProcess != kPyDirectGamma &&
+           fProcess != kPyMbNonDiffr) {
+           
+           for (i = 0; i<np; i++) {
                iparticle = (TParticle *) fParticles->At(i);
                Int_t ks = iparticle->GetStatusCode();
                kf = CheckPDGCode(iparticle->GetPdgCode());
@@ -334,26 +461,44 @@ void AliGenPythia::Generate()
                    TParticle *  mother = (TParticle *) fParticles->At(ipa);
                    kfMo = TMath::Abs(mother->GetPdgCode());
                }
-//             printf("\n particle (all)  %d %d %d", i, pSelected[i], kf);
-               if (kfl >= fFlavorSelect) { 
+               // What to keep in Stack?
+               Bool_t flavorOK = kFALSE;
+               Bool_t selectOK = kFALSE;
+               if (fFeedDownOpt) {
+                 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
+               } else {
+                 if (kfl > fFlavorSelect) {
+                   nc = -1;
+                   break;
+                 }
+                 if (kfl == fFlavorSelect) flavorOK = kTRUE;
+               }
+               switch (fStackFillOpt) {
+               case kFlavorSelection:
+                 selectOK = kTRUE;
+                 break;
+               case kParentSelection:
+                 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
+                 break;
+               }
+               if (flavorOK && selectOK) { 
 //
 // Heavy flavor hadron or quark
 //
 // Kinematic seletion on final state heavy flavor mesons
                    if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
                    {
-                       nc = -1;
-                       break;
+                     continue;
                    }
                    pSelected[i] = 1;
+                   if (ParentSelected(kf)) ++nParents; // Update parent count
 //                 printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
                } else {
 // Kinematic seletion on decay products
                    if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
                        && !KinematicSelection(iparticle, 1))
                    {
-                       nc = -1;
-                       break;
+                     continue;
                    }
 //
 // Decay products 
@@ -369,7 +514,7 @@ void AliGenPythia::Generate()
 // Semi-stable or de-selected: diselect decay products:
 // 
 //
-                       if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > 10e-15)
+                       if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
                        {
                            Int_t ipF = iparticle->GetFirstDaughter();
                            Int_t ipL = iparticle->GetLastDaughter();   
@@ -381,6 +526,8 @@ void AliGenPythia::Generate()
                }
                if (pSelected[i] == -1) pSelected[i] = 0;
                if (!pSelected[i]) continue;
+               // Count quarks only if you did not include fragmentation
+               if (fFragmentation && kf <= 10) continue;
                nc++;
 // Decision on tracking
                trackIt[i] = 0;
@@ -388,22 +535,24 @@ void AliGenPythia::Generate()
 // Track final state particle
                if (ks == 1) trackIt[i] = 1;
 // Track semi-stable particles
-               if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15))  trackIt[i] = 1;
+               if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 // Track particles selected by process if undecayed. 
                if (fForceDecay == kNoDecay) {
                    if (ParentSelected(kf)) trackIt[i] = 1;
                } else {
                    if (ParentSelected(kf)) trackIt[i] = 0;
                }
+               if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 //
 //
 
            } // particle selection loop
-           if (nc > -1) {
-               for (i = 0; i<np-1; i++) {
+           if (nc > 0) {
+               for (i = 0; i<np; i++) {
                    if (!pSelected[i]) continue;
                    TParticle *  iparticle = (TParticle *) fParticles->At(i);
                    kf = CheckPDGCode(iparticle->GetPdgCode());
+                   Int_t ks = iparticle->GetStatusCode();  
                    p[0] = iparticle->Px();
                    p[1] = iparticle->Py();
                    p[2] = iparticle->Pz();
@@ -414,7 +563,7 @@ void AliGenPythia::Generate()
                    Int_t ipa     = iparticle->GetFirstMother()-1;
                    Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
                    SetTrack(fTrackIt*trackIt[i] ,
-                                    iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
+                                    iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
                    pParent[i] = nt;
                    KeepTrack(nt); 
                } //  SetTrack loop
@@ -423,11 +572,35 @@ void AliGenPythia::Generate()
            nc = GenerateMB();
        } // mb ?
 
+       if (pParent)   delete[] pParent;
+       if (pSelected) delete[] pSelected;
+       if (trackIt)   delete[] trackIt;
+
        if (nc > 0) {
-           jev+=nc;
+         switch (fCountMode) {
+         case kCountAll:
+           // printf(" Count all \n");
+           jev += nc;
+           break;
+         case kCountParents:
+           // printf(" Count parents \n");
+           jev += nParents;
+           break;
+         case kCountTrackables:
+           // printf(" Count trackable \n");
+           jev += nTkbles;
+           break;
+         }
            if (jev >= fNpart || fNpart == -1) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
                printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+
+               fQ  += fPythia->GetVINT(51);
+               fX1 += fPythia->GetVINT(41);
+               fX2 += fPythia->GetVINT(42);
+               fTrialsRun += fTrials;
+               fNev++;
+               MakeHeader();
                break;
            }
        }
@@ -441,6 +614,9 @@ void AliGenPythia::Generate()
 
 Int_t  AliGenPythia::GenerateMB()
 {
+//
+// Min Bias selection and other global selections
+//
     Int_t i, kf, nt, iparent;
     Int_t nc = 0;
     Float_t p[3];
@@ -451,25 +627,28 @@ Int_t  AliGenPythia::GenerateMB()
     
     Int_t np = fParticles->GetEntriesFast();
     Int_t* pParent = new Int_t[np];
-    for (i=0; i< np-1; i++) pParent[i] = -1;
-
-    for (i = 0; i<np-1; i++) {
+    for (i=0; i< np; i++) pParent[i] = -1;
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+       TParticle* jet1 = (TParticle *) fParticles->At(6);
+       TParticle* jet2 = (TParticle *) fParticles->At(7);
+       if (!CheckTrigger(jet1, jet2)) return 0;
+    }
+    
+    for (i = 0; i<np; i++) {
        Int_t trackIt = 0;
        TParticle *  iparticle = (TParticle *) fParticles->At(i);
        kf = CheckPDGCode(iparticle->GetPdgCode());
        Int_t ks = iparticle->GetStatusCode();
        Int_t km = iparticle->GetFirstMother();
-//     printf("\n Particle: %d %d %d", i, kf, ks);
-       
        if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
            (ks != 1) ||
            (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
            nc++;
            if (ks == 1) trackIt = 1;
            Int_t ipa = iparticle->GetFirstMother()-1;
-
+           
            iparent = (ipa > -1) ? pParent[ipa] : -1;
-
+           
 //
 // store track information
            p[0] = iparticle->Px();
@@ -480,7 +659,7 @@ Int_t  AliGenPythia::GenerateMB()
            origin[2] = fOrigin[2]+iparticle->Vz()/10.;
            Float_t tof=kconv*iparticle->T();
            SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
-                        tof, kPPrimary, nt);
+                    tof, kPPrimary, nt, 1., ks);
            KeepTrack(nt);
            pParent[i] = nt;
        } // select particle
@@ -489,7 +668,6 @@ Int_t  AliGenPythia::GenerateMB()
     if (pParent) delete[] pParent;
     
     printf("\n I've put %i particles on the stack \n",nc);
-    MakeHeader();
     return nc;
 }
 
@@ -498,6 +676,13 @@ void AliGenPythia::FinishRun()
 {
 // Print x-section summary
     fPythia->Pystat(1);
+    fQ  /= fNev;
+    fX1 /= fNev;
+    fX2 /= fNev;    
+    printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
+    printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
+    
+
 }
 
 void AliGenPythia::AdjustWeights()
@@ -520,15 +705,67 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 }
 
 
-void AliGenPythia::MakeHeader()
+void AliGenPythia::MakeHeader() const
 {
 // Builds the event header, to be called after each event
     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
     ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+    ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
     header->SetPrimaryVertex(fEventVertex);
     gAlice->SetGenEventHeader(header);
 }
        
+
+Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const
+{
+// Check the kinematic trigger condition
+//
+    Double_t eta[2];
+    eta[0] = jet1->Eta();
+    eta[1] = jet2->Eta();
+    Double_t phi[2];
+    phi[0] = jet1->Phi();
+    phi[1] = jet2->Phi();
+    Int_t    pdg[2]; 
+    pdg[0] = jet1->GetPdgCode();
+    pdg[1] = jet2->GetPdgCode();    
+    Bool_t   triggered = kFALSE;
+
+    if (fProcess == kPyJets) {
+       //Check eta range first...
+       if (
+           ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet)   &&
+            (phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet)) 
+           
+           ||
+           
+           ((eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet)   &&
+            (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet))
+           )
+       {
+           triggered = kTRUE;
+       }
+    } else {
+       Int_t ij = 0;
+       Int_t ig = 1;
+       if (pdg[0] == kGamma) {
+           ij = 1;
+           ig = 0;
+       }
+       //Check eta range first...
+       if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
+           (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
+       {
+           //Eta is okay, now check phi range
+           if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
+               (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
+           {
+               triggered = kTRUE;
+           }
+       }
+    }
+    return triggered;
+}
          
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {