f4bfd574d93113c51d65561d98fb055e5b614d55
[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   if(AliLog::GetDebugLevel("TEvtGen","AliDecayerEvtGen") > 0) 
122     froot_part->printTree(); //to print the decay chain 
123   froot_part->deleteTree();
124   }
125
126 //____________________________________________________________
127 Int_t AliDecayerEvtGen::ImportParticles(TClonesArray *particles)
128   {
129   //
130   //Input: pointer to a TClonesArray - Output(Int_t): number of decay products    
131   //Put all the informations about the decay products in the 
132   //TClonesArray particles
133   //
134   if (particles == 0) return 0;
135   TClonesArray &clonesParticles = *particles;
136   clonesParticles.Clear();
137
138   int j;
139   int istat;
140   int partnum;
141   double px,py,pz,e;
142   double x,y,z,t;
143   EvtVector4R p4,x4;
144
145   Int_t npart=fEvtstdhep->getNPart();
146   for(int i=0;i<fEvtstdhep->getNPart();i++){
147   j=i+1;
148   int jmotherfirst=fEvtstdhep->getFirstMother(i)+1;
149   int jmotherlast=fEvtstdhep->getLastMother(i)+1;
150   int jdaugfirst=fEvtstdhep->getFirstDaughter(i)+1;
151   int jdauglast=fEvtstdhep->getLastDaughter(i)+1;
152
153   partnum=fEvtstdhep->getStdHepID(i);
154   
155   //verify if all particles of decay chain are in the TDatabasePDG
156   TParticlePDG *partPDG = TDatabasePDG::Instance()->GetParticle(partnum);
157   if(!partPDG)
158     {
159     AliWarning("Particle code non known in TDatabasePDG - set pdg = 89");
160     partnum=89; //internal use for unspecified resonance data
161     }
162
163   istat=fEvtstdhep->getIStat(i);
164
165   if(istat!=1 && istat!=2) Info("ImportParticles","Attention: unknown status code!");
166   if(istat == 2) istat = 11; //status decayed
167
168   p4=fEvtstdhep->getP4(i);
169   x4=fEvtstdhep->getX4(i);
170   px=p4.get(1);
171   py=p4.get(2);
172   pz=p4.get(3);
173   e=p4.get(0);
174   const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
175   const Float_t kconvL=1./10; // mm to cm conversion
176   x=x4.get(1)*kconvL;//[cm]
177   y=x4.get(2)*kconvL;//[cm]
178   z=x4.get(3)*kconvL;//[cm]
179   t=x4.get(0)*kconvT;//[s]
180
181   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));
182
183   new(clonesParticles[i]) TParticle(partnum,istat,jmotherfirst,-1,jdaugfirst,jdauglast,px,py,pz,e,x,y,z,t);
184
185   //set polarization!!!
186   }
187
188   return npart; 
189
190   }
191
192 void  AliDecayerEvtGen::Copy(TObject &) const
193   {
194   //
195   // Copy *this onto AliDecayerEvtGen -- not implemented
196   //
197    Fatal("Copy","Not implemented!\n");
198   }
199
200 void AliDecayerEvtGen::ForceDecay()
201   {
202   //
203   // Intupt: none - Output: none
204   // Set the decay mode to decay particles: for each case is read a 
205   // different decay table. case kAll read the default decay table only   
206   //
207   Decay_t decay = fDecay;
208   switch(decay)
209     {
210      case kAll: // particles decayed "naturally" according to $ALICE_ROOT/TEvtGen/EvtGen/DECAY.DEC
211       break;
212      case kBJpsiDiElectron:
213       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOELE.DEC"));
214       break;
215      case kBJpsi:
216       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSI.DEC"));
217       break;
218      case kBJpsiDiMuon:
219       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOMU.DEC"));
220       break;
221      case kBSemiElectronic:
222       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOELE.DEC"));
223       break;
224      case kHadronicD:
225       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HADRONICD.DEC"));
226       break;
227      case kHadronicDWithout4Bodies:
228       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HADRONICDWITHOUT4BODIES.DEC"));
229       break;
230      case kChiToJpsiGammaToElectronElectron:
231       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/CHICTOJPSITOELE.DEC"));
232       break;
233      case kChiToJpsiGammaToMuonMuon:
234       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/CHICTOJPSITOMUON.DEC"));
235       break;
236      case kSemiElectronic:
237       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BANDCTOELE.DEC"));
238       break;
239      case kBSemiMuonic:
240       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOMU.DEC"));
241       break;
242      case kSemiMuonic:
243       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BANDCTOMU.DEC"));
244       break;
245      case kDiElectron:
246       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIELECTRON.DEC"));
247       break;
248      case kDiMuon:
249       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIMUON.DEC"));
250       break;
251      case kBPsiPrimeDiMuon:
252       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOPSIPRIMETODIMUON.DEC"));
253       break;
254      case kBPsiPrimeDiElectron:
255       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOPSIPRIMETODIELECTRON.DEC"));
256       break;
257      case kJpsiDiMuon:
258       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/JPSIDIMUON.DEC"));
259       break;
260      case kPsiPrimeJpsiDiElectron:
261       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/PSIPRIMETOJPSITOMU.DEC"));
262       break;
263      case kPhiKK:
264       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/PHITOK.DEC"));
265       break;
266      case kOmega:
267       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/OMEGATOLAMBDAK.DEC"));
268       break;
269      case kLambda:
270       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/LAMBDATOPROTPI.DEC"));
271       break;
272      case kHardMuons:
273       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HARDMUONS.DEC"));
274       break;
275      case kElectronEM:
276       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/ELECTRONEM.DEC"));
277       break;
278      case kDiElectronEM:
279       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/DIELECTRONEM.DEC"));
280       break;
281      case kGammaEM:
282       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/GAMMAEM.DEC"));
283       break;
284      case kBeautyUpgrade:
285       SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BEAUTYUPGRADE.DEC"));
286       break;     
287      case kPiToMu:
288      case kKaToMu:
289      case kAllMuonic:
290      case kWToMuon:
291      case kWToCharm:
292      case kWToCharmToMuon:
293      case kZDiMuon:
294      case kZDiElectron:
295      case kNoDecay:
296      case kNoDecayHeavy:
297      case kNeutralPion:
298      case kBJpsiUndecayed: 
299      case kNoDecayBeauty:
300       AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
301      break;
302      }
303      ReadDecayTable();
304   }
305
306 Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
307   {
308    // This method is dummy
309    return  1.;
310   }
311
312 Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
313   {    
314   //
315   //Input: pdg code of a particle   
316   //return lifetime in sec for a particle with particle code kf
317   //
318   EvtId IdPart=EvtPDL::evtIdFromStdHep(kf); 
319   Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
320   AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
321   return lifetime*kconv; //tau (sec)
322   }
323
324 void AliDecayerEvtGen::ReadDecayTable()
325     {
326      //Input none - Output none 
327      //Read the decay table that correspond to the path 
328      //fDecayTablePath
329      // 
330      TString temp = fDecayTablePath;
331      if(!temp.EndsWith("DECAY.DEC"))
332      fGenerator->readUDecay(fDecayTablePath);
333     }
334 ////////////////////////////////////////////////////////////
335 Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
336   {
337   // 
338   //Set the path of the decay table read to force particle decays 
339   //
340   if(gSystem->AccessPathName(path)) 
341   {
342    AliWarning("Attention: This path not exist!\n");
343    return kFALSE;
344   }
345   fDecayTablePath = path;
346   return kTRUE;
347   } 
348
349
350