TEvtGen/AliDecayerEvtGen.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 //  Implementation of AliDecayer using EvtGen package.                   //
  17 //                                                                       //
  18 //  Giuseppe E. Bruno            &    Fiorella Fionda                    //
  19 //  (Giuseppe.Bruno@ba.infn.it)       (Fiorella.Fionda@ba.infn.it)       //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22 #include <TSystem.h>
  23 #include <TDatabasePDG.h>
  24 #include <TParticle.h>
  25 #include <TClonesArray.h>
  26 #include <TLorentzVector.h>
  27
  28 #include "EvtGenBase/EvtStdHep.hh"
  29 #include "EvtGenBase/EvtRandomEngine.hh"
  30 #include "EvtGen/EvtGen.hh"
  31 #include "EvtGenBase/EvtParticle.hh"
  32 #include "EvtGenBase/EvtPDL.hh"
  33 #include "EvtGenBase/EvtParticleFactory.hh"
  34 #include "AliDecayerEvtGen.h"
  35 #include "AliLog.h"
  36
  37 ClassImp(AliDecayerEvtGen)
  38 //____________________________________________________________
  39 AliDecayerEvtGen::AliDecayerEvtGen():
  40   fRandomEngine(0x0),
  41   fGenerator(0x0),
  42   fEvtstdhep(0x0),
  43   fDecayTablePath(0x0),
  44   fParticleTablePath(0x0),
  45   fDecay(kAll)
  46   {
  47   // Default constructor
  48   fEvtstdhep = new EvtStdHep();
  49   fDecayTablePath = gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DECAY.DEC"); //default decay table
  50   fParticleTablePath = gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/evt.pdl"); //particle table
  51   }
  52 //_________________________________________________________________
  53 AliDecayerEvtGen::AliDecayerEvtGen(const AliDecayerEvtGen &decayer):
  54   AliDecayer(decayer),
  55   fRandomEngine(decayer.fRandomEngine),
  56   fGenerator(decayer.fGenerator),
  57   fEvtstdhep(decayer.fEvtstdhep),
  58   fDecayTablePath(decayer.fDecayTablePath),
  59   fParticleTablePath(decayer.fParticleTablePath),
  60   fDecay(decayer.fDecay)
  61   {
  62   // Copy Constructor
  63   decayer.Copy(*this);
  64   }
  65 //____________________________________________________________
  66 AliDecayerEvtGen::~AliDecayerEvtGen()
  67   {
  68   // Destructor
  69   if(fRandomEngine) {delete fRandomEngine;}
  70   fRandomEngine = 0;
  71   if(fGenerator) {delete fGenerator;}
  72   fGenerator = 0;
  73   if(fEvtstdhep) {delete fEvtstdhep;}
  74   fEvtstdhep = 0;
  75   if(fDecayTablePath) {delete fDecayTablePath;}
  76   fDecayTablePath = 0;
  77   if(fParticleTablePath) {delete fParticleTablePath;}
  78   fParticleTablePath = 0;
  79   }
  80
  81 //___________________________________________________________
  82 void AliDecayerEvtGen::Init()
  83   {
  84   //Standard AliDecayerEvtGen initializer:
  85   //initialize EvtGen with default decay table (DECAY.DEC), particle table (evt.pdl)
  86   //and fRandomEngine for generation of random numbers
  87   //
  88   if(fGenerator){
  89   AliWarning(" AliDecayerEvtGen already initialized!!!!\n");
  90   return;
  91   }
  92   fRandomEngine=new EvtNUMRandomEngine();
  93   fGenerator=new EvtGen(fDecayTablePath,fParticleTablePath,fRandomEngine);
  94   }
  95 //____________________________________________________________
  96 void AliDecayerEvtGen::Decay(Int_t ipart, TLorentzVector *p)
  97   {
  98   //
  99   //Decay a particle
 100   //input: pdg code and momentum of the particle to be decayed
 101   //all informations about decay products are stored in fEvtstdhep
 102   //
 103   EvtId IPART=EvtPDL::evtIdFromStdHep(ipart);
 104   EvtVector4R p_init(p->E(),p->Px(),p->Py(),p->Pz());
 105   EvtParticle *froot_part=EvtParticleFactory::particleFactory(IPART,p_init);
 106   fGenerator->generateDecay(froot_part);
 107   fEvtstdhep->init();
 108   froot_part->makeStdHep(*fEvtstdhep);
 109   //froot_part->printTree(); //to print the decay chain
 110   froot_part->deleteTree();
 111   }
 112
 113 //____________________________________________________________
 114 Int_t AliDecayerEvtGen::ImportParticles(TClonesArray *particles)
 115   {
 116   //
 117   //Input: pointer to a TClonesArray - Output(Int_t): number of decay products
 118   //Put all the informations about the decay products in the
 119   //TClonesArray particles
 120   //
 121   if (particles == 0) return 0;
 122   TClonesArray &clonesParticles = *particles;
 123   clonesParticles.Clear();
 124
 125   int j;
 126   int istat;
 127   int partnum;
 128   double px,py,pz,e;
 129   double x,y,z,t;
 130   EvtVector4R p4,x4;
 131
 132   Int_t npart=fEvtstdhep->getNPart();
 133   for(int i=0;i<fEvtstdhep->getNPart();i++){
 134   j=i+1;
 135   int jmotherfirst=fEvtstdhep->getFirstMother(i)+1;
 136   int jmotherlast=fEvtstdhep->getLastMother(i)+1;
 137   int jdaugfirst=fEvtstdhep->getFirstDaughter(i)+1;
 138   int jdauglast=fEvtstdhep->getLastDaughter(i)+1;
 139
 140   partnum=fEvtstdhep->getStdHepID(i);
 141
 142   //verify if all particles of decay chain are in the TDatabasePDG
 143   TParticlePDG *partPDG = TDatabasePDG::Instance()->GetParticle(partnum);
 144   if(!partPDG)
 145     {
 146     AliWarning("Particle code non known in TDatabasePDG - set pdg = 89");
 147     partnum=89; //internal use for unspecified resonance data
 148     }
 149
 150   istat=fEvtstdhep->getIStat(i);
 151
 152   if(istat!=1 && istat!=2) Info("ImportParticles","Attention: unknown status code!");
 153   if(istat == 2) istat = 11; //status decayed
 154
 155   p4=fEvtstdhep->getP4(i);
 156   x4=fEvtstdhep->getX4(i);
 157   px=p4.get(1);
 158   py=p4.get(2);
 159   pz=p4.get(3);
 160   e=p4.get(0);
 161   const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
 162   const Float_t kconvL=1./10; // mm to cm conversion
 163   x=x4.get(1)*kconvL;//[cm]
 164   y=x4.get(2)*kconvL;//[cm]
 165   z=x4.get(3)*kconvL;//[cm]
 166   t=x4.get(0)*kconvT;//[s]
 167
 168   AliDebug(1,Form("partnum = %d istat = %d primaMadre = %d ultimaMadre = %d primaF = %d ultimaF=%d x=%f y=%f z=%f t=%f e=%f px=%f \n",partnum,istat,jmotherfirst,jmotherlast,jdaugfirst,jdauglast,x,y,z,t,e,px));
 169
 170   new(clonesParticles[i]) TParticle(partnum,istat,jmotherfirst,-1,jdaugfirst,jdauglast,px,py,pz,e,x,y,z,t);
 171
 172   //set polarization!!!
 173   }
 174
 175   return npart;
 176
 177   }
 178
 179 void  AliDecayerEvtGen::Copy(TObject &) const
 180   {
 181   //
 182   // Copy *this onto AliDecayerEvtGen -- not implemented
 183   //
 184    Fatal("Copy","Not implemented!\n");
 185   }
 186
 187 void AliDecayerEvtGen::ForceDecay()
 188   {
 189   //
 190   // Intupt: none - Output: none
 191   // Set the decay mode to decay particles: for each case is read a
 192   // different decay table. case kAll read the default decay table only
 193   //
 194   Decay_t decay = fDecay;
 195   switch(decay)
 196     {
 197      case kAll:
 198      break;
 199      case kBJpsiDiElectron:
 200      SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOELE.DEC"));
 201      break;
 202      case kBJpsi:
 203      SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSI.DEC"));
 204      break;
 205      case kBJpsiDiMuon:
 206      SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOMU.DEC"));
 207      break;
 208      case kBSemiElectronic:
 209      SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOELE.DEC"));
 210      break;
 211      case kHardMuons:
 212      case kChiToJpsiGammaToMuonMuon:
 213      case kChiToJpsiGammaToElectronElectron:
 214      case kBSemiMuonic:
 215      case kSemiMuonic:
 216      case kDiMuon:
 217      case kSemiElectronic:
 218      case kDiElectron:
 219      case kBPsiPrimeDiMuon:
 220      case kPiToMu:
 221      case kKaToMu:
 222      case kAllMuonic:
 223      case kWToMuon:
 224      case kWToCharm:
 225      case kWToCharmToMuon:
 226      case kZDiMuon:
 227      case kZDiElectron:
 228      case kHadronicD:
 229      case kHadronicDWithout4Bodies:
 230      case kPhiKK:
 231      case kOmega:
 232      case kLambda:
 233      case kNoDecay:
 234      case kNoDecayHeavy:
 235      case kNeutralPion:
 236      case kBPsiPrimeDiElectron:
 237      case kBeautyUpgrade:
 238      AliWarning(Form("Warning: case %s not implemented for this class!",(Char_t)decay));
 239      break;
 240      }
 241      ReadDecayTable();
 242   }
 243
 244 Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
 245   {
 246    // This method is dummy
 247    return  1.;
 248   }
 249
 250 Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
 251   {
 252   //
 253   //Input: pdg code of a particle
 254   //return lifetime in sec for a particle with particle code kf
 255   //
 256   EvtId IdPart=EvtPDL::evtIdFromStdHep(kf);
 257   Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
 258   AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
 259   return lifetime*kconv; //tau (sec)
 260   }
 261
 262 void AliDecayerEvtGen::ReadDecayTable()
 263     {
 264      //Input none - Output none
 265      //Read the decay table that correspond to the path
 266      //fDecayTablePath
 267      //
 268      TString temp = fDecayTablePath;
 269      if(!temp.EndsWith("DECAY.DEC"))
 270      fGenerator->readUDecay(fDecayTablePath);
 271     }
 272 ////////////////////////////////////////////////////////////
 273 Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
 274   {
 275   //
 276   //Set the path of the decay table read to force particle decays
 277   //
 278   if(gSystem->AccessPathName(path))
 279   {
 280    AliWarning("Attention: This path not exist!\n");
 281    return kFALSE;
 282   }
 283   fDecayTablePath = path;
 284   return kTRUE;
 285   }
 286
 287
 288