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