Option to write only the heavy flavor chain.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jan 2012 16:59:11 +0000 (16:59 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jan 2012 16:59:11 +0000 (16:59 +0000)
PYTHIA6/AliGenPythia.cxx
PYTHIA6/AliGenPythia.h

index b6cc0ef..10843ce 100644 (file)
@@ -757,6 +757,7 @@ void AliGenPythia::Generate()
                    if (kfl == fFlavorSelect) flavorOK = kTRUE;
                }
                switch (fStackFillOpt) {
+               case kHeavyFlavor:
                case kFlavorSelection:
                    selectOK = kTRUE;
                    break;
@@ -1165,10 +1166,12 @@ Int_t  AliGenPythia::GenerateMB()
        kf = CheckPDGCode(iparticle->GetPdgCode());
        Int_t ks = iparticle->GetStatusCode();
        Int_t km = iparticle->GetFirstMother();
-       if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
-           (ks != 1) ||
-           ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
-           nc++;
+       if (
+           (((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
+           ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
+           ) 
+         {
+            nc++;
            if (ks == 1) trackIt = 1;
            Int_t ipa = iparticle->GetFirstMother()-1;
            
@@ -2359,3 +2362,27 @@ wSD = -1;
 
   return kFALSE;
 }
+
+Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
+{
+// Check if this is a heavy flavor decay product
+  TParticle *  part = (TParticle *) fParticles.At(ipart);
+  Int_t mpdg = TMath::Abs(part->GetPdgCode());
+  Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
+  //
+  // Light hadron
+  if (mfl >= 4 && mfl < 6) return kTRUE;
+  
+  Int_t imo = part->GetFirstMother()-1;
+  TParticle* pm = part;
+  //
+  // Heavy flavor hadron produced by generator
+  while (imo >  -1) {
+    pm  =  (TParticle*)fParticles.At(imo);
+    mpdg = TMath::Abs(pm->GetPdgCode());
+    mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
+    if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
+    imo = pm->GetFirstMother()-1;
+  }
+  return kFALSE;
+}
index 4b5a253..38dbf36 100644 (file)
@@ -31,7 +31,7 @@ class AliGenPythia : public AliGenMC
 {
  public:
 
-    typedef enum {kFlavorSelection, kParentSelection} StackFillOpt_t;
+    typedef enum {kFlavorSelection, kParentSelection, kHeavyFlavor} StackFillOpt_t;
     typedef enum {kCountAll, kCountParents, kCountTrackables} CountMode_t;
     typedef enum {kCluster, kCell} JetRecMode_t;
          
@@ -187,6 +187,7 @@ class AliGenPythia : public AliGenMC
     //
     Bool_t IsInEMCAL(Float_t phi, Float_t eta) const;
     Bool_t IsInPHOS(Float_t phi, Float_t eta) const;
+    Bool_t IsFromHeavyFlavor(Int_t ipart);
     //
     virtual void FinishRun();
     Bool_t CheckTrigger(const TParticle* jet1, const TParticle* jet2);