b99b54dcf37832861d03cf1b902a1a38877affaf
[u/mrichter/AliRoot.git] / TEvtGen / 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/HADRONICD.DEC"));
225       break;
226      case kHadronicDWithout4Bodies:
227       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HADRONICDWITHOUT4BODIES.DEC"));
228       break;
229      case kChiToJpsiGammaToElectronElectron:
230       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/CHICTOJPSITOELE.DEC"));
231       break;
232      case kChiToJpsiGammaToMuonMuon:
233       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/CHICTOJPSITOMUON.DEC"));
234       break;
235      case kSemiElectronic:
236       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BANDCTOELE.DEC"));
237       break;
238      case kBSemiMuonic:
239       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOMU.DEC"));
240       break;
241      case kSemiMuonic:
242       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BANDCTOMU.DEC"));
243       break;
244      case kDiElectron:
245       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIELECTRON.DEC"));
246       break;
247      case kDiMuon:
248       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIMUON.DEC"));
249       break;
250      case kBPsiPrimeDiMuon:
251       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOPSIPRIMETODIMUON.DEC"));
252       break;
253      case kBPsiPrimeDiElectron:
254       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOPSIPRIMETODIELECTRON.DEC"));
255       break;
256      case kJpsiDiMuon:
257       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/JPSIDIMUON.DEC"));
258       break;
259      case kPsiPrimeJpsiDiElectron:
260       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/PSIPRIMETOJPSITOMU.DEC"));
261       break;
262      case kPhiKK:
263       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/PHITOK.DEC"));
264       break;
265      case kOmega:
266       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/OMEGATOLAMBDAK.DEC"));
267       break;
268      case kLambda:
269       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/LAMBDATOPROTPI.DEC"));
270       break;
271      case kHardMuons:
272       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HARDMUONS.DEC"));
273       break;
274      case kElectronEM:
275       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/ELECTRONEM.DEC"));
276       break;
277      case kDiElectronEM:
278       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIELECTRONEM.DEC"));
279       break;
280      case kGammaEM:
281       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/GAMMAEM.DEC"));
282       break;
283      case kBeautyUpgrade:
284       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BEAUTYUPGRADE.DEC"));
285       break;     
286      case kPiToMu:
287      case kKaToMu:
288      case kAllMuonic:
289      case kWToMuon:
290      case kWToCharm:
291      case kWToCharmToMuon:
292      case kZDiMuon:
293      case kZDiElectron:
294      case kNoDecay:
295      case kNoDecayHeavy:
296      case kNeutralPion:
297      case kBJpsiUndecayed: 
298      case kNoDecayBeauty:
299       AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
300      break;
301      }
302      ReadDecayTable();
303   }
304
305 Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
306   {
307    // This method is dummy
308    return  1.;
309   }
310
311 Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
312   {    
313   //
314   //Input: pdg code of a particle   
315   //return lifetime in sec for a particle with particle code kf
316   //
317   EvtId IdPart=EvtPDL::evtIdFromStdHep(kf); 
318   Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
319   AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
320   return lifetime*kconv; //tau (sec)
321   }
322
323 void AliDecayerEvtGen::ReadDecayTable()
324     {
325      //Input none - Output none 
326      //Read the decay table that correspond to the path 
327      //fDecayTablePath
328      // 
329      TString temp = fDecayTablePath;
330      if(!temp.EndsWith("DECAY.DEC"))
331      fGenerator->readUDecay(fDecayTablePath);
332     }
333 ////////////////////////////////////////////////////////////
334 Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
335   {
336   // 
337   //Set the path of the decay table read to force particle decays 
338   //
339   if(gSystem->AccessPathName(path)) 
340   {
341    AliWarning("Attention: This path not exist!\n");
342    return kFALSE;
343   }
344   fDecayTablePath = path;
345   return kTRUE;
346   } 
347
348
349