New stack-fill and count options introduced (N. Carrer).

diff --git a/EVGEN/AliGenPythia.cxx b/EVGEN/AliGenPythia.cxx
index 3d9a6bd..8886114 100644 (file)
--- a/EVGEN/AliGenPythia.cxx
+++ b/EVGEN/AliGenPythia.cxx
@@ -203,6 +203,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetEventListRange();
     SetJetPhiRange();
     SetJetEtaRange();
+    // 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)
@@ -240,6 +247,13 @@ void AliGenPythia::Init()
     fPythia->SetCKIN(3,fPtHardMin);
     fPythia->SetCKIN(4,fPtHardMax);    
     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);
@@ -256,6 +270,10 @@ void AliGenPythia::Init()
        fParentSelect[3] =  4122;
        fFlavorSelect    =  4;  
        break;
+    case kPyD0PbMNR:
+       fParentSelect[0] =   421;
+       fFlavorSelect    =  4;  
+       break;
     case kPyBeauty:
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
@@ -345,9 +363,12 @@ void AliGenPythia::Generate()
        for (i=0; i< np; i++) {
            pParent[i]   = -1;
            pSelected[i] =  0;
+           trackIt[i]   =  0;
        }
-       printf("\n **************************************************%d\n",np);
-       Int_t nc = 0;
+       // printf("\n **************************************************%d\n",np);
+       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) {
            
            for (i = 0; i<np; i++) {
@@ -376,25 +397,44 @@ void AliGenPythia::Generate()
                    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 
@@ -422,6 +462,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;
@@ -436,16 +478,16 @@ void AliGenPythia::Generate()
                } else {
                    if (ParentSelected(kf)) trackIt[i] = 0;
                }
+               if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 //
 //
 
            } // particle selection loop
-           if (nc > -1) {
+           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();
@@ -456,7 +498,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, ks);
+                                    iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
                    pParent[i] = nt;
                    KeepTrack(nt); 
                } //  SetTrack loop
@@ -470,7 +512,20 @@ void AliGenPythia::Generate()
        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);
@@ -532,7 +587,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, 1., ks);
+                        tof, kPPrimary, nt);
            KeepTrack(nt);
            pParent[i] = nt;
        } // select particle

diff --git a/EVGEN/AliGenPythia.h b/EVGEN/AliGenPythia.h
index b7d2408..095f80c 100644 (file)
--- a/EVGEN/AliGenPythia.h
+++ b/EVGEN/AliGenPythia.h
@@ -16,6 +16,10 @@ class TParticle;
 class AliGenPythia : public AliGenMC
 {
  public:
+
+  typedef enum {kFlavorSelection, kParentSelection} StackFillOpt_t;
+  typedef enum {kCountAll, kCountParents, kCountTrackables} CountMode_t;
+
     AliGenPythia();
     AliGenPythia(Int_t npart);
     AliGenPythia(const AliGenPythia &Pythia);
@@ -44,6 +48,23 @@ class AliGenPythia : public AliGenMC
        {fEtaMinGamma = etamin; fEtaMaxGamma = etamax;}
     virtual void    SetGammaPhiRange(Float_t phimin = -180., Float_t phimax = 180.)
        {fPhiMinGamma = TMath::Pi()*phimin/180.; fPhiMaxGamma = TMath::Pi()*phimax/180.;}
+    // Set option for feed down from higher family
+    virtual void SetFeedDownHigherFamily(Bool_t opt) {
+      fFeedDownOpt = opt;
+    }
+    // Set option for selecting particles kept in stack according to flavor
+    // or to parent selection
+    virtual void SetStackFillOpt(StackFillOpt_t opt) {
+      fStackFillOpt = opt;
+    }
+    // Set fragmentation option
+    virtual void SetFragmentation(const Bool_t opt) {
+      fFragmentation = opt;
+    }
+    // Set counting mode
+    virtual void SetCountMode(const CountMode_t mode) {
+      fCountMode = mode;
+    }
     
     // get cross section of process
     virtual Float_t GetXsection() {return fXsection;}      
@@ -84,6 +105,21 @@ class AliGenPythia : public AliGenMC
     Float_t     fPhiMinGamma;    // Minimum phi of triggered gamma
     Float_t     fPhiMaxGamma;    // Maximum phi of triggered gamma
 
+    StackFillOpt_t fStackFillOpt; // Stack filling with all particles with
+                                  // that flavour or only with selected
+                                  // parents and their decays
+    Bool_t fFeedDownOpt;          // Option to set feed down from higher
+                                  // quark families (e.g. b->c)
+    Bool_t fFragmentation;        // Option to activate fragmentation by Pythia
+    //
+    // Options for counting when the event will be finished.
+    // fCountMode = kCountAll         --> All particles that end up in the
+    //                                    stack are counted
+    // fCountMode = kCountParents     --> Only selected parents are counted
+    // fCountMode = kCountTrackabless --> Only particles flagged for tracking
+    //                                     are counted
+    CountMode_t fCountMode;
+
  private:
     // adjust the weight from kinematic cuts
     void   AdjustWeights();