TEvtGen/AliGenEvtGen.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14 ***************************************************************************/
  15 ///////////////////////////////////////////////////////////////////////////
  16 //  Class to generate decays of particles generated by a                 //
  17 //  previous generator. It works as a generator, but pratically it       //
  18 //  performs only decays. It works with this scheme: first loops over    //
  19 //  particles on the stack, selects those to be decayed, decays them     //
  20 //  and then pushes the decay products on the stack.                     //
  21 //                                                                       //
  22 //  Giuseppe E. Bruno            &    Fiorella Fionda                    //
  23 //  (Giuseppe.Bruno@ba.infn.it)       (Fiorella.Fionda@ba.infn.it)       //
  24 ///////////////////////////////////////////////////////////////////////////
  25
  26 #include "AliStack.h"
  27 #include "AliGenEvtGen.h"
  28 #include "AliRun.h"
  29 #include "AliLog.h"
  30 #include <TParticle.h>
  31
  32 ClassImp(AliGenEvtGen)
  33 ///////////////////////////////////////////////////////////////////////////
  34 AliGenEvtGen::AliGenEvtGen():
  35   fStack(0x0),
  36   fDecayer(0x0),
  37   fForceDecay(kAll),
  38   fSwitchOff(kBeautyPart),
  39   fUserDecay(kFALSE),
  40   fUserDecayTablePath(0x0)
  41   {
  42   //
  43   // Default Construction
  44   //
  45   }
  46 ///////////////////////////////////////////////////////////////////////////////////////////
  47 AliGenEvtGen::~AliGenEvtGen()
  48   {
  49   //
  50   // Standard Destructor
  51   //
  52   if(fStack) {delete fStack;}
  53   fStack = 0;
  54   if(fDecayer) {delete fDecayer;}
  55   fDecayer = 0;
  56   if(fUserDecayTablePath) {delete fUserDecayTablePath;}
  57   fUserDecayTablePath = 0;
  58   }
  59 ///////////////////////////////////////////////////////////////////////////////////////////
  60
  61 void AliGenEvtGen::Init()
  62   {
  63   //
  64   // Standard AliGenerator Initializer - no input
  65   // 1) initialize EvtGen with default decay and particle table
  66   // 2) set the decay mode to force particle
  67   // 3) set a user decay table if defined
  68   //
  69   if(fDecayer)
  70   {
  71   AliWarning("AliGenEvtGen already initialized!!!");
  72   return;
  73   }
  74   fDecayer = new AliDecayerEvtGen();
  75   fDecayer->Init(); //read the default decay table DECAY.DEC and particle table
  76
  77   //if is set a decay mode: default decay mode is kAll
  78   fDecayer->SetForceDecay(fForceDecay);
  79   fDecayer->ForceDecay();
  80
  81   //if is defined a user decay table
  82   if(fUserDecay)
  83     {
  84     fDecayer->SetDecayTablePath(fUserDecayTablePath);
  85     fDecayer->ReadDecayTable();
  86     }
  87   }
  88
  89 /////////////////////////////////////////////////////////////////////////////////////////////
  90 void AliGenEvtGen::Generate()
  91   {
  92   //
  93   //Generate method - Input none - Output none
  94   //For each event:
  95   //1)return the stack of the previous generator and select particles to be decayed by EvtGen
  96   //2)decay particles selected and put the decay products on the stack
  97   //
  98   //
  99   Float_t polar[3]= {0,0,0};  // Polarisation of daughter particles
 100   Float_t origin0[3];         // Origin of the parent particle
 101   Float_t pc[3], och[3]; // Momentum and origin of the children particles from EvtGen
 102   Int_t nt;
 103   Float_t tof;
 104   Int_t nPrimsPythia;
 105   TLorentzVector *mom=new TLorentzVector();
 106   static TClonesArray *particles;
 107   if(!particles) particles = new TClonesArray("TParticle",1000);
 108   fStack = AliRunLoader::Instance()->Stack();
 109   if(!fStack) {Info("Generate","Error: No stack found!"); return;}
 110   nPrimsPythia = fStack->GetNprimary();
 111   AliDebug(1,Form("nPrimsPythia = %d \n",nPrimsPythia));
 112   for (Int_t iTrack = 0; iTrack < nPrimsPythia; ++iTrack) {
 113   TParticle *part = fStack->Particle(iTrack);
 114   Int_t pdg=part->GetPdgCode();
 115
 116   AliDebug(1,Form("GetFlavour = %d e pdg = %d \n",GetFlavour(pdg),pdg));
 117
 118   switch(fSwitchOff)
 119   {
 120   case kAllPart:
 121   break;
 122   case kBeautyPart:
 123   if(GetFlavour(pdg)!=5) continue;
 124   break;
 125   case kCharmPart:
 126   if(GetFlavour(pdg)!=4) continue;
 127   break;
 128   }
 129
 130   //check if particle is already decayed by Pythia
 131   if(part->GetStatusCode() != 1 || part->GetNDaughters()>0)
 132     {
 133     Info("AliGenEvtGen","Attention: particle %d is already decayed by Pythia!",pdg);
 134     continue;
 135     }
 136
 137   part->SetStatusCode(11); //Set particle as decayed : change the status code
 138   part->SetBit(kDoneBit);
 139   part->ResetBit(kTransportBit);
 140
 141   mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
 142   Int_t np;
 143
 144    do{
 145      fDecayer->Decay(part->GetPdgCode(),mom);
 146      np = fDecayer->ImportParticles(particles);
 147      }while(np<0);
 148
 149   Int_t* trackIt    = new Int_t[np];
 150   Int_t* pParent   = new Int_t[np];
 151   AliDebug(1,Form("np = %d \n",np));
 152
 153   for (int i = 0; i < np; i++) {
 154              pParent[i]   = -1;
 155              trackIt[i]   =  0;
 156          }
 157   //select trackable particle
 158   if (np >1) {
 159   TParticle* iparticle =  (TParticle *) particles->At(0);//parent particle
 160    for (int i = 1; i<np ; i++) {
 161      iparticle = (TParticle*) particles->At(i);
 162      Int_t ks = iparticle->GetStatusCode();
 163
 164      //track last decay products
 165      if(ks==1) trackIt[i]=1;
 166
 167      }//decay particles loop
 168
 169   }// if decay products
 170
 171   origin0[0]=part->Vx(); //[cm]
 172   origin0[1]=part->Vy(); //[cm]
 173   origin0[2]=part->Vz(); //[cm]
 174   //
 175   // Put decay products on the stack
 176   //
 177   for (int i = 1; i < np; i++) {
 178          TParticle* iparticle = (TParticle *) particles->At(i);
 179          Int_t kf   = iparticle->GetPdgCode();
 180          Int_t ksc  = iparticle->GetStatusCode();
 181          Int_t jpa  = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
 182          Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;
 183
 184          och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
 185          och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
 186          och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
 187          pc[0]  = iparticle->Px(); //[GeV/c]
 188          pc[1]  = iparticle->Py(); //[GeV/c]
 189          pc[2]  = iparticle->Pz(); //[GeV/c]
 190          tof = part->T()+kconv*iparticle->T();
 191
 192          AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));
 193
 194          PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc);
 195          if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
 196          pParent[i] = nt;
 197          KeepTrack(nt);
 198          SetHighWaterMark(nt);
 199          }// Particle loop
 200          particles->Clear();
 201          if (trackIt)    delete[] trackIt;
 202          if (pParent)    delete[] pParent;
 203        }
 204   Info("Generate","AliGenEvtGen DONE");
 205   }
 206
 207 //////////////////////////////////////////////////////////////////////////////////////////
 208 Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)
 209   {
 210   //
 211   // return the flavour of a particle
 212   // input: pdg code of the particle
 213   // output: Int_t
 214   //         3 in case of strange (open and hidden)
 215   //         4 in case of charm (")
 216   //         5 in case of beauty (")
 217   //
 218   Int_t pdg = TMath::Abs(pdgCode);
 219   //Resonance
 220   if (pdg > 100000) pdg %= 100000;
 221   if(pdg > 10000)  pdg %= 10000;
 222   // meson ?
 223   if(pdg > 10) pdg/=100;
 224   // baryon ?
 225   if(pdg > 10) pdg/=10;
 226   return pdg;
 227   }
 228
 229 /////////////////////////////////////////////////////////////////////////////////////////
 230 Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
 231   {
 232   //
 233   //Set the path of user decay table if it is defined
 234   //
 235   //put a comment to control if path exists
 236   if(gSystem->AccessPathName(path))
 237     {
 238     AliWarning("Attention: This path not exist!\n");
 239     return kFALSE;
 240     }
 241   fUserDecayTablePath = path;
 242   fUserDecay = kTRUE;
 243   return kTRUE;
 244   }