- Some coded moved to AliGenMC
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 13 Jul 2001 10:58:54 +0000 (10:58 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 13 Jul 2001 10:58:54 +0000 (10:58 +0000)
- Improved handling of secondary vertices.

EVGEN/AliGenParam.cxx
EVGEN/AliGenParam.h
EVGEN/AliGenPythia.cxx
EVGEN/AliGenPythia.h

index 80dc91c..07dcca9 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.30  2001/06/15 07:55:04  morsch
+Put only first generation decay products on the stack.
+
 Revision 1.29  2001/03/27 10:58:41  morsch
 Initialize decayer before generation. Important if run inside cocktail.
 
@@ -111,18 +114,13 @@ AliGenParam::AliGenParam()
     fYPara  = 0;
     fParam  = 0;
     fAnalog = kAnalog;
-    SetCutOnChild();
-    SetChildMomentumRange();
-    SetChildPtRange();
-    SetChildPhiRange();
-    SetChildThetaRange();  
     SetDeltaPt();
 //
 // Set random number generator   
     sRandom = fRandom;
 }
 
-AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenerator(npart)
+AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenMC(npart)
 {
 // Constructor using number of particles parameterisation id and library
     
@@ -134,14 +132,7 @@ AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* t
     fYPara  = 0;
     fParam  = param;
     fAnalog = kAnalog;
-    fChildSelect.Set(5);
-    for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
     SetForceDecay();
-    SetCutOnChild();
-    SetChildMomentumRange();
-    SetChildPtRange();
-    SetChildPhiRange();
-    SetChildThetaRange(); 
     SetDeltaPt(); 
 //
 // Set random number generator   
@@ -150,7 +141,7 @@ AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* t
 
 //____________________________________________________________
 
-AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenerator(npart)
+AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenMC(npart)
 {
 // Constructor using parameterisation id and number of particles
 //      
@@ -179,7 +170,7 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param,
                          Double_t (*PtPara) (Double_t*, Double_t*),
                          Double_t (*YPara ) (Double_t* ,Double_t*),
                         Int_t    (*IpPara) (TRandom *))                 
-    :AliGenerator(npart)
+    :AliGenMC(npart)
 {
 // Constructor
 // Gines Martinez 1/10/99 
@@ -265,34 +256,7 @@ void AliGenParam::Init()
     fDecayer->Init();
 
 //
-    switch (fForceDecay) 
-    {
-    case kSemiElectronic:
-    case kDiElectron:
-    case kBJpsiDiElectron:
-    case kBPsiPrimeDiElectron:
-       fChildSelect[0]=11;     
-       break;
-    case kSemiMuonic:
-    case kDiMuon:
-    case kBJpsiDiMuon:
-    case kBPsiPrimeDiMuon:
-       fChildSelect[0]=13;
-       break;
-    case kPiToMu:
-       fChildSelect[0]=13;
-       break;
-    case kKaToMu:
-       fChildSelect[0]=13;
-       break;
-    case kHadronicD:
-// Implement me !!
-       break;
-    case kNoDecay:
-       break;
-    case kAll:
-       break;
-    }
+    AliGenMC::Init();
 }
 
 //____________________________________________________________
@@ -317,14 +281,14 @@ void AliGenParam::Generate()
   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
   Float_t p[3], pc[3], 
-          och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
+          och[3];             // Momentum, polarisation and origin of the children particles from lujet
   Float_t ty, xmt;
-  Int_t nt, i, j, kfch[10];
+  Int_t nt, i, j;
   Float_t  wgtp, wgtch;
   Double_t dummy;
   static TClonesArray *particles;
   //
-  if(!particles) particles=new TClonesArray("TParticle",1000);
+  if(!particles) particles = new TClonesArray("TParticle",1000);
 
   static TDatabasePDG *pDataBase = new TDatabasePDG();
   if(!pDataBase) pDataBase = new TDatabasePDG();
@@ -334,18 +298,12 @@ void AliGenParam::Generate()
 // Calculating vertex position per event
   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
   if(fVertexSmear==kPerEvent) {
-//      Rndm(random,6);
-//      for (j=0;j<3;j++) {
-//       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
-//           TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-//           TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-//      }
-//    }
       Vertex();
       for (j=0;j<3;j++) origin0[j]=fVertex[j];
   }
   
   Int_t ipa=0;
+  
 // Generating fNpart particles
   while (ipa<fNpart) {
       while(1) {
@@ -384,6 +342,7 @@ void AliGenParam::Generate()
          ptot=TMath::Sqrt(pt*pt+pl*pl);
 // Cut on momentum
          if(ptot<fPMin || ptot>fPMax) continue;
+//
          p[0]=pt*TMath::Cos(phi);
          p[1]=pt*TMath::Sin(phi);
          p[2]=pl;
@@ -403,6 +362,7 @@ void AliGenParam::Generate()
 // if fForceDecay == none Primary particle is tracked by GEANT 
 // (In the latest, make sure that GEANT actually does all the decays you want)   
 //
+
          if (fForceDecay != kNoDecay) {
 // Using lujet to decay particle
              Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
@@ -412,50 +372,71 @@ void AliGenParam::Generate()
 // select decay particles
              Int_t np=fDecayer->ImportParticles(particles);
              Int_t ncsel=0;
-
+             Int_t* pFlag      = new Int_t[np];
+             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; i++) {
+                 pFlag[i]     =  0;
+                 pSelected[i] =  0;
+                 pParent[i]   = -1;
+             }
+             
              if (np >1) {
                  TParticle* iparticle =  (TParticle *) particles->At(0);
-                 Int_t ipF = iparticle->GetFirstDaughter();
-                 Int_t ipL = iparticle->GetLastDaughter();           
-                 for (i = ipF-1; i<ipL; i++) {
+                 Int_t ipF, ipL;
+                 for (i = 1; i<np ; i++) {
+                     trackIt[i] = 1;
                      iparticle = (TParticle *) particles->At(i);
                      Int_t kf = iparticle->GetPdgCode();
+                     Int_t ks = iparticle->GetStatusCode();
+// flagged particle
+
+                     if (pFlag[i] == 1) {
+                         printf("\n deselected %d", i);
+                         ipF = iparticle->GetFirstDaughter();
+                         ipL = iparticle->GetLastDaughter();   
+                         if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+                         continue;
+                     }
+
+// flag decay products of particles with long life-time (c tau > .3 mum)                     
+                     
+                     if (ks != 1) { 
+                         TParticlePDG *particle = pDataBase->GetParticle(kf);
+                         
+                         Double_t lifeTime = fDecayer->GetLifetime(kf);
+                         Double_t mass     = particle->Mass();
+                         Double_t width    = particle->Width();
+                         printf("\n lifetime %d %e %e %e ", i, lifeTime, mass, width);
+                         if (lifeTime > (Double_t) 1.e-15) {
+                             ipF = iparticle->GetFirstDaughter();
+                             ipL = iparticle->GetLastDaughter();       
+                             if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+                         } else{
+                             trackIt[i]     = 0;
+                             pSelected[i]   = 1;
+                         }
+                     } // ks==1 ?
 //
 // children
-                     if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll)
+                     if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
                      {
-                         pc[0]=iparticle->Px();
-                         pc[1]=iparticle->Py();
-                         pc[2]=iparticle->Pz();
-                         och[0]=origin0[0]+iparticle->Vx()/10;
-                         och[1]=origin0[1]+iparticle->Vy()/10;
-                         och[2]=origin0[2]+iparticle->Vz()/10;
                          if (fCutOnChild) {
-                             Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
-                             Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]);
-                             Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
-                             Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
-                             Bool_t childok = 
-                                 ((ptChild    > fChildPtMin    && ptChild    <fChildPtMax)      &&
-                                  (pChild     > fChildPMin     && pChild     <fChildPMax)       &&
-                                  (thetaChild > fChildThetaMin && thetaChild <fChildThetaMax)   &&
-                                  (phiChild   > fChildPhiMin   && phiChild   <fChildPhiMax));
-                             if(childok)
-                             {
-                                 pch[ncsel][0]=pc[0];
-                                 pch[ncsel][1]=pc[1];
-                                 pch[ncsel][2]=pc[2];
-                                 kfch[ncsel]=kf;
+                             pc[0]=iparticle->Px();
+                             pc[1]=iparticle->Py();
+                             pc[2]=iparticle->Pz();
+                             Bool_t  childok = KinematicSelection(iparticle, 1);
+                             if(childok) {
+                                 pSelected[i]  = 1;
                                  ncsel++;
                              } else {
                                  ncsel=-1;
                                  break;
                              } // child kine cuts
                          } else {
-                             pch[ncsel][0]=pc[0];
-                             pch[ncsel][1]=pc[1];
-                             pch[ncsel][2]=pc[2];
-                             kfch[ncsel]=kf;
+                             pSelected[i]  = 1;
                              ncsel++;
                          } // if child selection
                      } // select muon
@@ -466,18 +447,45 @@ void AliGenParam::Generate()
              if ((fCutOnChild && ncsel >0) || !fCutOnChild){
                  ipa++;
 //
-// parent
-                 gAlice->
-                     SetTrack(0,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
-                 iparent=nt;
+// Parent
+                 gAlice->SetTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+                 pParent[0] = nt;
                  gAlice->KeepTrack(nt); 
-                 for (i=0; i< ncsel; i++) {
-                     gAlice->SetTrack(fTrackIt,iparent,kfch[i],
-                                      &pch[i][0],och,polar,
-                                      0,kPDecay,nt,wgtch);
-                     gAlice->KeepTrack(nt); 
+//
+// Decay Products
+//               
+                 for (i = 1; i < np; i++) {
+                     if (pSelected[i]) {
+                         TParticle* iparticle = (TParticle *) particles->At(i);
+                         Int_t kf  = iparticle->GetPdgCode();
+                         Int_t ipa = iparticle->GetFirstMother()-1;
+                         
+                         och[0] = origin0[0]+iparticle->Vx()/10;
+                         och[1] = origin0[1]+iparticle->Vy()/10;
+                         och[2] = origin0[2]+iparticle->Vz()/10;
+                         pc[0]  = iparticle->Px();
+                         pc[1]  = iparticle->Py();
+                         pc[2]  = iparticle->Pz();
+                         
+                         if (ipa > -1) {
+                             iparent = pParent[ipa];
+                         } else {
+                             iparent = -1;
+                         }
+                         printf("\n SetTrack %d %d %d %d", i, kf, iparent, trackIt[i]);
+                         gAlice->SetTrack(fTrackIt*trackIt[i], iparent, kf,
+                                          pc, och, polar,
+                                          0, kPDecay, nt, wgtch);
+                         pParent[i] = nt;
+                         gAlice->KeepTrack(nt); 
+                     }
                  }
              }  // Decays by Lujet
+
+             if (pFlag)      delete[] pFlag;
+             if (pParent)    delete[] pParent;
+             if (pSelected)  delete[] pSelected;          
+             if (trackIt)    delete[] trackIt;     
          } // kinematic selection
          else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
          {
@@ -489,48 +497,8 @@ void AliGenParam::Generate()
     } // while
   } // event loop
   gAlice->SetHighWaterMark(nt);
-  
 }
 
-Bool_t AliGenParam::ChildSelected(Int_t ip)
-{
-// True if particle is in list of selected children
-    for (Int_t i=0; i<5; i++)
-    {
-       if (fChildSelect[i]==ip) return kTRUE;
-    }
-    return kFALSE;
-}
-
-Bool_t AliGenParam::KinematicSelection(TParticle *particle)
-{
-// Perform kinematic cuts
-    Float_t px=particle->Px();
-    Float_t py=particle->Py();
-    Float_t pz=particle->Pz();
-//
-// momentum cut
-    Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
-    if (p > fPMax || p < fPMin) 
-    {
-//     printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
-       return kFALSE;
-    }
-    Float_t pt=TMath::Sqrt(px*px+py*py);
-    
-//
-// theta cut
-    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
-    if (theta > fThetaMax || theta < fThetaMin) 
-    {
-//     printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
-       return kFALSE;
-    }
-
-    return kTRUE;
-}
-
-
 AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
 {
 // Assignment operator
index ec57eec..dfd4dc8 100644 (file)
@@ -5,8 +5,7 @@
 
 /* $Id$ */
 
-#include "AliGenerator.h"
-#include "AliDecayer.h"
+#include "AliGenMC.h"
 #include <TArrayI.h>
 
 class AliPythia;
@@ -16,7 +15,7 @@ class TF1;
 
 typedef enum { kAnalog, kNonAnalog} Weighting_t;
 //-------------------------------------------------------------
-class AliGenParam : public AliGenerator
+class AliGenParam : public AliGenMC
 {
  public:
     AliGenParam();
@@ -34,19 +33,7 @@ class AliGenParam : public AliGenerator
     // select particle type
     virtual void SetParam(Int_t param) {fParam = param;}
     // force decay type
-    virtual void SetForceDecay(Decay_t decay = kDiMuon) {fForceDecay = decay;}
     virtual void SetWeighting(Weighting_t flag = kAnalog) {fAnalog = flag;}    
-    virtual void SetCutOnChild(Int_t flag = 0) {fCutOnChild = flag;}
-    virtual void SetChildMomentumRange(Float_t pmin = 0, Float_t pmax = 1.e10)
-       {fChildPMin = pmin; fChildPMax = pmax;}
-    virtual void SetChildPtRange(Float_t ptmin = 0, Float_t ptmax = 20.)
-       {fChildPtMin = ptmin; fChildPtMax = ptmax;}
-    virtual void SetChildPhiRange(Float_t phimin = -180., Float_t phimax = 180)
-       {fChildPhiMin = TMath::Pi()*phimin/180;
-       fChildPhiMax  = TMath::Pi()*phimax/180;}
-    virtual void SetChildThetaRange(Float_t thetamin = 0, Float_t thetamax = 180)
-       {fChildThetaMin = TMath::Pi()*thetamin/180;
-       fChildThetaMax  = TMath::Pi()*thetamax/180;}
     virtual void SetDeltaPt(Float_t delta=0.01) {fDeltaPt = delta;}
     
     AliGenParam & operator=(const AliGenParam & rhs);
@@ -62,25 +49,8 @@ class AliGenParam : public AliGenerator
     Float_t     fPtWgt;        // Pt-weight
     Float_t     fBias;         // Biasing factor
     Int_t       fTrials;       // Number of trials
-    Decay_t     fForceDecay;   // Decay channel forced
-    Int_t       fCutOnChild;   // Cuts on decay products (children)  are enabled/disabled
-    Float_t     fChildPtMin;   // Children minimum pT
-    Float_t     fChildPtMax;   // Children maximum pT
-    Float_t     fChildPMin;    // Children minimum p
-    Float_t     fChildPMax;    // Children maximum p
-    Float_t     fChildPhiMin;  // Children minimum phi
-    Float_t     fChildPhiMax;  // Children maximum phi
-    Float_t     fChildThetaMin;// Children minimum theta
-    Float_t     fChildThetaMax;// Children maximum theta
     Float_t     fDeltaPt;      // pT sampling in steps of fDeltaPt
-    TArrayI     fChildSelect;  // Children to be selected from decay products
     AliDecayer  *fDecayer;     // ! Pointer to pythia object for decays
- private:
-    // check if particle is selected as child
-    Bool_t ChildSelected(Int_t ip);
-    // all kinematic selection goes here
-    Bool_t KinematicSelection(TParticle *particle);
-
   ClassDef(AliGenParam,1) // Generator using parameterised pt- and y-distribution
 };
 #endif
index 418910f..d4cf4fc 100644 (file)
 
 /*
 $Log$
+Revision 1.37  2001/06/28 11:17:28  morsch
+SetEventListRange setter added. Events in specified range are listed for
+debugging. (Yuri Kharlov)
+
 Revision 1.36  2001/03/30 07:05:49  morsch
 Final print-out in finish run.
 Write parton system for jet-production (preliminary solution).
@@ -102,6 +106,7 @@ Introduction of the Copyright and cvs Log
 */
 
 #include "AliGenPythia.h"
+#include "AliGenPythiaEventHeader.h"
 #include "AliDecayerPythia.h"
 #include "AliRun.h"
 #include "AliPythia.h"
@@ -112,7 +117,7 @@ Introduction of the Copyright and cvs Log
  ClassImp(AliGenPythia)
 
 AliGenPythia::AliGenPythia()
-                 :AliGenerator()
+                 :AliGenMC()
 {
 // Default Constructor
   fDecayer = new AliDecayerPythia();
@@ -120,7 +125,7 @@ AliGenPythia::AliGenPythia()
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
-                 :AliGenerator(npart)
+                 :AliGenMC(npart)
 {
 // default charm production at 5. 5 TeV
 // semimuonic decay
@@ -129,9 +134,6 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fXsection  = 0.;
     fNucA1=0;
     fNucA2=0;
-    fParentSelect.Set(5);
-    fChildSelect.Set(5);
-    for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
     SetProcess();
     SetStrucFunc();
     SetForceDecay();
@@ -141,6 +143,9 @@ AliGenPythia::AliGenPythia(Int_t npart)
     // Set random number generator 
     sRandom=fRandom;
     SetEventListRange();
+    fFlavorSelect   = 0;
+    // Produced particles  
+    fParticles = new TClonesArray("TParticle",1000);
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
@@ -171,10 +176,6 @@ void AliGenPythia::Init()
     fParentWeight=1./Float_t(fNpart);
 //
 //  Forward Paramters to the AliPythia object
-    //    gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1");
-    //    fPythia->Pyupda(2,1);    
-    //    gSystem->Exec("rm fort.1");
-
     fDecayer->SetForceDecay(fForceDecay);    
     fDecayer->Init();
 
@@ -190,63 +191,49 @@ void AliGenPythia::Init()
     switch (fProcess) 
     {
     case kPyCharm:
-       fParentSelect[0]=411;
-       fParentSelect[1]=421;
-       fParentSelect[2]=431;
-       fParentSelect[3]=4122;  
+       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;  
+       fParentSelect[0] =  411;
+       fParentSelect[1] =  421;
+       fParentSelect[2] =  431;
+       fParentSelect[3]=  4122;
+       fFlavorSelect    = 4;   
        break;
     case kPyBeauty:
-       fParentSelect[0]=511;
-       fParentSelect[1]=521;
-       fParentSelect[2]=531;
-       fParentSelect[3]=5122;  
+       fParentSelect[0]=  511;
+       fParentSelect[1]=  521;
+       fParentSelect[2]=  531;
+       fParentSelect[3]= 5122;
+       fParentSelect[4]= 5132;
+       fParentSelect[5]= 5232;
+       fParentSelect[6]= 5332;
+       fFlavorSelect   = 5;    
        break;
     case kPyBeautyUnforced:
-       fParentSelect[0]=511;
-       fParentSelect[1]=521;
-       fParentSelect[2]=531;
-       fParentSelect[3]=5122;  
+       fParentSelect[0] =  511;
+       fParentSelect[1] =  521;
+       fParentSelect[2] =  531;
+       fParentSelect[3] = 5122;
+       fParentSelect[4] = 5132;
+       fParentSelect[5] = 5232;
+       fParentSelect[6] = 5332;
+       fFlavorSelect    = 5;   
        break;
     case kPyJpsiChi:
     case kPyJpsi:
-       fParentSelect[0]=443;
+       fParentSelect[0] = 443;
        break;
     case kPyMb:
     case kPyJets:
     case kPyDirectGamma:
        break;
     }
-
-    switch (fForceDecay) 
-    {
-    case kSemiElectronic:
-    case kDiElectron:
-    case kBJpsiDiElectron:
-    case kBPsiPrimeDiElectron:
-       fChildSelect[0]=kElectron;      
-       break;
-    case kSemiMuonic:
-    case kDiMuon:
-    case kBJpsiDiMuon:
-    case kBPsiPrimeDiMuon:
-    case kPiToMu:
-    case kKaToMu:
-       fChildSelect[0]=kMuonMinus;
-       break;
-    case kHadronicD:
-      fChildSelect[0]=kPiPlus;
-      fChildSelect[1]=kKPlus;
-      break;
-    case kAll:
-    case kNoDecay:
-      break;
-    }
+    AliGenMC::Init();
 }
 
 void AliGenPythia::Generate()
@@ -256,138 +243,147 @@ void AliGenPythia::Generate()
 
     Float_t polar[3]   =   {0,0,0};
     Float_t origin[3]  =   {0,0,0};
-    Float_t originP[3] =   {0,0,0};
     Float_t origin0[3] =   {0,0,0};
-    Float_t p[3], pP[4];
-//    Float_t random[6];
-    static TClonesArray *particles;
+    Float_t p[3];
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
-    
-    
 //
     Int_t nt=0;
-    Int_t ntP=0;
     Int_t jev=0;
     Int_t j, kf;
-
-    if(!particles) particles=new TClonesArray("TParticle",1000);
-    
     fTrials=0;
+
+//  Set collision vertex position 
     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
     if(fVertexSmear==kPerEvent) {
        fPythia->SetMSTP(151,1);
        for (j=0;j<3;j++) {
-           fPythia->SetPARP(151+j, fOsigma[j]/10.);
+           fPythia->SetPARP(151+j, fOsigma[j]*10.);
        }
     } else if (fVertexSmear==kPerTrack) {
        fPythia->SetMSTP(151,0);
     }
-    
+//  event loop    
     while(1)
     {
        fPythia->Pyevnt();
        if (gAlice->GetEvNumber()>=fDebugEventFirst &&
            gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
        fTrials++;
-       fPythia->ImportParticles(particles,"All");
-       Int_t np = particles->GetEntriesFast();
+       fPythia->ImportParticles(fParticles,"All");
+
+//
+//
+//
+       Int_t i;
+       
+       Int_t np = fParticles->GetEntriesFast();
+       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++) {
+           pParent[i]   = -1;
+           pSelected[i] =  0;
+       }
        printf("\n **************************************************%d\n",np);
-       Int_t nc=0;
+       Int_t nc = 0;
        if (np == 0 ) continue;
        if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
-           for (Int_t i = 0; i<np-1; i++) {
-               TParticle *  iparticle = (TParticle *) particles->At(i);
+           for (i = 0; i<np-1; i++) {
+               TParticle *  iparticle = (TParticle *) fParticles->At(i);
                Int_t ks = iparticle->GetStatusCode();
                kf = CheckPDGCode(iparticle->GetPdgCode());
+// No initial state partons
                if (ks==21) continue;
-               fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight;      
 //
-// Parent
-               if (ParentSelected(TMath::Abs(kf))) {
-                   if (KinematicSelection(iparticle)) {
-                       if (nc==0) {
+// Heavy Flavor Selection
 //
-// Store information concerning the hard scattering process
+               // quark ?
+               kf = TMath::Abs(kf);
+               Int_t kfl = kf;
+               // meson ?
+               if  (kfl > 10) kfl/=100;
+               // baryon
+               if (kfl > 10) kfl/=10;
+               if (kfl > 10) kfl/=10;
+
+               Int_t ipa = iparticle->GetFirstMother()-1;
+               Int_t kfMo = 0;
+               
+               if (ipa > -1) {
+                   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) { 
 //
-                           Float_t massP  = fPythia->GetPARI(13);
-                           Float_t   ptP  = fPythia->GetPARI(17);
-                           Float_t    yP  = fPythia->GetPARI(37);
-                           Float_t  xmtP  = sqrt(ptP*ptP+massP*massP);
-                           Float_t    ty  = Float_t(TMath::TanH(yP));
-                           pP[0] = ptP;
-                           pP[1] = 0;
-                           pP[2] = xmtP*ty/sqrt(1.-ty*ty);
-                           pP[3] = massP;
-                           gAlice->SetTrack(0,-1,-1,
-                                            pP,originP,polar,
-                                            0,kPPrimary,ntP,fParentWeight);
-//                                          0,"Hard Scat.",ntP,fParentWeight);
-                           gAlice->KeepTrack(ntP);
-                       }
-                       nc++;
+// Heavy flavor hadron or quark
 //
-// store parent track information
-                       p[0]=iparticle->Px();
-                       p[1]=iparticle->Py();
-                       p[2]=iparticle->Pz();
-                       origin[0]=origin0[0]+iparticle->Vx()/10.;
-                       origin[1]=origin0[1]+iparticle->Vy()/10.;
-                       origin[2]=origin0[2]+iparticle->Vz()/10.;
-
-                       Int_t ifch=iparticle->GetFirstDaughter();
-                       Int_t ilch=iparticle->GetLastDaughter();        
-
-                       if ((ifch !=0 && ilch !=0) || fForceDecay == kNoDecay) {
-                         Int_t trackit=0;
-                         if (fForceDecay == kNoDecay) trackit = 1;
-                           gAlice->SetTrack(trackit,ntP,kf,
-                                            p,origin,polar,
-                                            0,kPPrimary,nt,fParentWeight);
-                           gAlice->KeepTrack(nt);
-                           Int_t iparent = nt;
+// Kinematic seletion on final state heavy flavor mesons
+                   if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
+                   {
+                       nc = -1;
+                       break;
+                   }
+                   pSelected[i] = 1;
+//                 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;
+                   }
 //
-// Children        
-                           if (fForceDecay != kNoDecay) {
-                             for (j=ifch; j<=ilch; j++)
-                               {
-                                 TParticle *  ichild = 
-                                   (TParticle *) particles->At(j-1);
-                                 kf = CheckPDGCode(ichild->GetPdgCode());
+// Decay products 
+// Select if mother was selected and is not tracked
+
+                   if (pSelected[ipa] && 
+                       !trackIt[ipa]  &&     // mother will be  tracked ?
+                       kfMo !=  5 &&         // mother is b-quark, don't store fragments          
+                       kfMo !=  4 &&         // mother is c-quark, don't store fragments 
+                       kf   != 92)           // don't store string
+                   {
 //
+// Semi-stable or de-selected: diselect decay products:
 // 
-                                 if (ChildSelected(TMath::Abs(kf))) {
-                                   origin[0]=origin0[0]+ichild->Vx()/10.;
-                                   origin[1]=origin0[1]+ichild->Vy()/10.;
-                                   origin[2]=origin0[2]+ichild->Vz()/10.;              
-                                   p[0]=ichild->Px();
-                                   p[1]=ichild->Py();
-                                   p[2]=ichild->Pz();
-                                   Float_t tof=kconv*ichild->T();
-                                   gAlice->SetTrack(fTrackIt, iparent, kf,
-                                                    p,origin,polar,
-                                                    tof,kPDecay,nt,fChildWeight);
-                                   gAlice->KeepTrack(nt);
-                                 } // select child
-                               } // child loop
-                           } 
-                       }
-                   } // kinematic selection
-               } // select particle
-           } // particle loop
-       } else {
-           for (Int_t i = 0; i<np-1; i++) {
-               TParticle *  iparticle = (TParticle *) particles->At(i);
-               kf = CheckPDGCode(iparticle->GetPdgCode());
-               Int_t ks = iparticle->GetStatusCode();
-               Int_t km = iparticle->GetFirstMother();
-               //      printf("\n process %d %d\n", ks,km);
-               
-               if ((ks==1 && kf!=0 && KinematicSelection(iparticle)) ||
-                   (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
-                   nc++;
 //
-// store track information
+                       if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > 10e-15)
+                       {
+                           Int_t ipF = iparticle->GetFirstDaughter();
+                           Int_t ipL = iparticle->GetLastDaughter();   
+                           if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
+                       }
+//                     printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
+                       pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
+                   }
+               }
+               if (pSelected[i] == -1) pSelected[i] = 0;
+               if (!pSelected[i]) continue;
+               nc++;
+// Decision on tracking
+               trackIt[i] = 0;
+//
+// 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;
+// Track particles selected by process if undecayed. 
+               if (fForceDecay == kNoDecay) {
+                   if (ParentSelected(kf)) trackIt[i] = 1;
+               } else {
+                   if (ParentSelected(kf)) trackIt[i] = 0;
+               }
+//
+//
+
+           } // particle selection loop
+           if (nc > -1) {
+               for (i = 0; i<np-1; i++) {
+                   if (!pSelected[i]) continue;
+                   TParticle *  iparticle = (TParticle *) fParticles->At(i);
+                   kf = CheckPDGCode(iparticle->GetPdgCode());
                    p[0]=iparticle->Px();
                    p[1]=iparticle->Py();
                    p[2]=iparticle->Pz();
@@ -395,13 +391,19 @@ void AliGenPythia::Generate()
                    origin[1]=origin0[1]+iparticle->Vy()/10.;
                    origin[2]=origin0[2]+iparticle->Vz()/10.;
                    Float_t tof=kconv*iparticle->T();
-                   gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
-                                    tof,kPPrimary,nt);
-                   gAlice->KeepTrack(nt);
-               } // select particle
-           } // particle loop 
-           printf("\n I've put %i particles on the stack \n",nc);
+                   Int_t ipa = iparticle->GetFirstMother()-1;
+                   Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
+//                 printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]);
+                   gAlice->SetTrack(fTrackIt*trackIt[i] ,
+                                    iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
+                   pParent[i] = nt;
+                   gAlice->KeepTrack(nt); 
+               } //  SetTrack loop
+           }
+       } else {
+           nc = GenerateMB();
        } // mb ?
+
        if (nc > 0) {
            jev+=nc;
            if (jev >= fNpart || fNpart == -1) {
@@ -418,93 +420,68 @@ void AliGenPythia::Generate()
     fXsection=fPythia->GetPARI(1);
 }
 
-void AliGenPythia::FinishRun()
-{
-// Print x-section summary
-    fPythia->Pystat(1);
-}
-
-Bool_t AliGenPythia::ParentSelected(Int_t ip)
-{
-// True if particle is in list of parent particles to be selected
-    for (Int_t i=0; i<5; i++)
-    {
-       if (fParentSelect[i]==ip) return kTRUE;
-    }
-    return kFALSE;
-}
-
-Bool_t AliGenPythia::ChildSelected(Int_t ip)
+Int_t  AliGenPythia::GenerateMB()
 {
-// True if particle is in list of decay products to be selected
-    if (fForceDecay == kAll) return kTRUE;
+    Int_t i, kf, nt, iparent;
+    Int_t nc = 0;
+    Float_t p[3];
+    Float_t polar[3]   =   {0,0,0};
+    Float_t origin[3]  =   {0,0,0};
+    Float_t origin0[3] =   {0,0,0};
+//  converts from mm/c to s
+    const Float_t kconv=0.001/2.999792458e8;
     
-    for (Int_t i=0; i<5; i++)
-    {
-       if (fChildSelect[i]==ip) return kTRUE;
-    }
-    return kFALSE;
-}
+    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++) {
+       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;
 
-Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
-{
-// Perform kinematic selection
-    Float_t px=particle->Px();
-    Float_t py=particle->Py();
-    Float_t pz=particle->Pz();
-    Float_t  e=particle->Energy();
+           iparent = (ipa > -1) ? pParent[ipa] : -1;
 
 //
-//  transverse momentum cut    
-    Float_t pt=TMath::Sqrt(px*px+py*py);
-    if (pt > fPtMax || pt < fPtMin) 
-    {
-//     printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
-       return kFALSE;
-    }
-//
-// momentum cut
-    Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
-    if (p > fPMax || p < fPMin) 
-    {
-//     printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
-       return kFALSE;
-    }
+// store track information
+           p[0]=iparticle->Px();
+           p[1]=iparticle->Py();
+           p[2]=iparticle->Pz();
+           origin[0]=origin0[0]+iparticle->Vx()/10.;
+           origin[1]=origin0[1]+iparticle->Vy()/10.;
+           origin[2]=origin0[2]+iparticle->Vz()/10.;
+           Float_t tof=kconv*iparticle->T();
+           gAlice->SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
+                            tof, kPPrimary, nt);
+           gAlice->KeepTrack(nt);
+           pParent[i] = nt;
+       } // select particle
+    } // particle loop 
+
+    if (pParent) delete[] pParent;
     
-//
-// theta cut
-    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
-    if (theta > fThetaMax || theta < fThetaMin) 
-    {
-//     printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
-       return kFALSE;
-    }
-
-//
-// rapidity cut
-    if ( (e-pz)<=0 || (e+pz)<=0 ) {
-      return kFALSE;
-    }
-    else {
-      Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
-      if (y > fYMax || y < fYMin)
-        {
-//     printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
-          return kFALSE;
-        }
-    }
+    printf("\n I've put %i particles on the stack \n",nc);
+    MakeHeader();
+    return nc;
+}
 
-//
-// phi cut
-    Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
-    if (phi > fPhiMax || phi < fPhiMin)
-    {
-//     printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
-       return kFALSE;
-    }
 
-    return kTRUE;
+void AliGenPythia::FinishRun()
+{
+// Print x-section summary
+    fPythia->Pystat(1);
 }
+
 void AliGenPythia::AdjustWeights()
 {
 // Adjust the weights after generation of all events
@@ -516,40 +493,6 @@ void AliGenPythia::AdjustWeights()
        part->SetWeight(part->GetWeight()*fKineBias);
     }
 }
-
-Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
-{
-//
-//  If the particle is in a diffractive state, then take action accordingly
-  switch (pdgcode) {
-  case 91:
-    return 92;
-  case 110:
-    //rho_diff0 -- difficult to translate, return rho0
-    return 113;
-  case 210:
-    //pi_diffr+ -- change to pi+
-    return 211;
-  case 220:
-    //omega_di0 -- change to omega0
-    return 223;
-  case 330:
-    //phi_diff0 -- return phi0
-    return 333;
-  case 440:
-    //J/psi_di0 -- return J/psi
-    return 443;
-  case 2110:
-    //n_diffr -- return neutron
-    return 2112;
-  case 2210:
-    //p_diffr+ -- return proton
-    return 2212;
-  }
-  //non diffractive state -- return code unchanged
-  return pdgcode;
-}
-
     
 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 {
@@ -557,6 +500,15 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
     fNucA1 = a1;
     fNucA2 = a2;
 }
+
+
+void AliGenPythia::MakeHeader()
+{
+// Builds the event header, to be called after each event
+    AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
+   ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+   gAlice->SetGenEventHeader(header);
+}
        
          
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
index 5057608..daf0e3d 100644 (file)
@@ -6,15 +6,13 @@
 /* $Id$ */
 
 
-#include "AliGenerator.h"
+#include "AliGenMC.h"
 #include "AliPythia.h"
-#include "AliDecayer.h"
-#include <TArrayI.h>    
 
 class AliPythia;
 class TParticle;
 
-class AliGenPythia : public AliGenerator
+class AliGenPythia : public AliGenMC
 {
  public:
     AliGenPythia();
@@ -33,27 +31,26 @@ class AliGenPythia : public AliGenerator
        {fPtHardMin = ptmin; fPtHardMax = ptmax; }
     // set centre of mass energy
     virtual void    SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;}
-    // force decay type
-    virtual void    SetForceDecay(Decay_t decay = kSemiMuonic) {fForceDecay = decay;}
     // treat protons as inside nuclei
     virtual void    SetNuclei(Int_t a1, Int_t a2);
     // get cross section of process
     virtual Float_t GetXsection() {return fXsection;}      
-    // Check PDG code
-    virtual Int_t   CheckPDGCode(Int_t pdgcode);
     virtual void    FinishRun();
     
     // Assignment Operator
     AliGenPythia & operator=(const AliGenPythia & rhs);
+ private:
+    Int_t  GenerateMB();
+    virtual void    MakeHeader();    
  protected:
+    TClonesArray* fParticles;  // Particle  List
+    
     Process_t   fProcess;       // Process type
     StrucFunc_t fStrucFunc;     // Structure Function
-    Decay_t     fForceDecay;    // Decay channel  are forced
     Float_t     fEnergyCMS;     // Centre of mass energy
     Float_t     fKineBias;      // Bias from kinematic selection
     Int_t       fTrials;        // Number of trials
-    TArrayI     fParentSelect;  // Parent particles to be selected 
-    TArrayI     fChildSelect;   // Decay products to be selected
+    Int_t       fFlavorSelect;  // Heavy Flavor Selection
     Float_t     fXsection;      // Cross-section
     AliPythia   *fPythia;       //! Pythia 
     Float_t     fPtHardMin;     // lower pT-hard cut 
@@ -61,16 +58,10 @@ class AliGenPythia : public AliGenerator
     Int_t       fNucA1;         // mass number nucleus side 1
     Int_t       fNucA2;         // mass number nucleus side 2
     Bool_t      fFullEvent;     // Write Full event if true
-    AliDecayer  *fDecayer;      // ! pointer to the decayer instance
+    AliDecayer  *fDecayer;      //! Pointer to the decayer instance
     Int_t       fDebugEventFirst; // First event to debug
     Int_t       fDebugEventLast;  // Last  event to debug
  private:
-    // check if particle is selected as parent particle
-    Bool_t ParentSelected(Int_t ip);
-    // check if particle is selected as child particle
-    Bool_t ChildSelected(Int_t ip);
-    // all kinematic selection cuts go here 
-    Bool_t KinematicSelection(TParticle *particle);
     // adjust the weight from kinematic cuts
     void   AdjustWeights();