The decay tree is printed only if the debug level is at least one.
[u/mrichter/AliRoot.git] / TEvtGen / TEvtGen / AliDecayerEvtGen.cxx
CommitLineData
da0e9ce3 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"
0ca57c2f 30#include "EvtGenBase/EvtStdlibRandomEngine.hh"
da0e9ce3 31#include "EvtGen/EvtGen.hh"
32#include "EvtGenBase/EvtParticle.hh"
33#include "EvtGenBase/EvtPDL.hh"
34#include "EvtGenBase/EvtParticleFactory.hh"
35#include "AliDecayerEvtGen.h"
0ca57c2f 36#include "EvtGenExternal/EvtExternalGenList.hh"
37#include "EvtGenBase/EvtAbsRadCorr.hh"
da0e9ce3 38#include "AliLog.h"
39
40ClassImp(AliDecayerEvtGen)
41//____________________________________________________________
42AliDecayerEvtGen::AliDecayerEvtGen():
43 fRandomEngine(0x0),
0ca57c2f 44 fRadCorrEngine(0x0),
da0e9ce3 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//_________________________________________________________________
57AliDecayerEvtGen::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//____________________________________________________________
70AliDecayerEvtGen::~AliDecayerEvtGen()
71 {
72 // Destructor
73 if(fRandomEngine) {delete fRandomEngine;}
74 fRandomEngine = 0;
0ca57c2f 75 if(fRadCorrEngine) {delete fRadCorrEngine;}
76 fRadCorrEngine = 0;
da0e9ce3 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//___________________________________________________________
88void 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 }
0ca57c2f 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);
da0e9ce3 106 }
107//____________________________________________________________
108void 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);
34dbffe9 116 EvtVector4R p_init(p->E(),p->Px(),p->Py(),p->Pz());
da0e9ce3 117 EvtParticle *froot_part=EvtParticleFactory::particleFactory(IPART,p_init);
118 fGenerator->generateDecay(froot_part);
119 fEvtstdhep->init();
120 froot_part->makeStdHep(*fEvtstdhep);
38515416 121 if(AliLog::GetDebugLevel("TEvtGen","AliDecayerEvtGen") > 0)
122 froot_part->printTree(); //to print the decay chain
da0e9ce3 123 froot_part->deleteTree();
124 }
125
126//____________________________________________________________
127Int_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);
774ceaaf 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]
da0e9ce3 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
192void AliDecayerEvtGen::Copy(TObject &) const
193 {
194 //
195 // Copy *this onto AliDecayerEvtGen -- not implemented
196 //
197 Fatal("Copy","Not implemented!\n");
198 }
199
200void 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 {
0ca57c2f 210 case kAll: // particles decayed "naturally" according to $ALICE_ROOT/TEvtGen/EvtGen/DECAY.DEC
3a6cdcc8 211 break;
da0e9ce3 212 case kBJpsiDiElectron:
3a6cdcc8 213 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOELE.DEC"));
214 break;
da0e9ce3 215 case kBJpsi:
3a6cdcc8 216 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSI.DEC"));
217 break;
da0e9ce3 218 case kBJpsiDiMuon:
3a6cdcc8 219 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOMU.DEC"));
220 break;
da0e9ce3 221 case kBSemiElectronic:
3a6cdcc8 222 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOELE.DEC"));
223 break;
941c36b4 224 case kHadronicD:
de323cb7 225 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HADRONICD.DEC"));
941c36b4 226 break;
de323cb7 227 case kHadronicDWithout4Bodies:
228 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/HADRONICDWITHOUT4BODIES.DEC"));
229 break;
230 case kChiToJpsiGammaToElectronElectron:
3a6cdcc8 231 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/CHICTOJPSITOELE.DEC"));
232 break;
da0e9ce3 233 case kChiToJpsiGammaToMuonMuon:
3a6cdcc8 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;
da0e9ce3 239 case kBSemiMuonic:
3a6cdcc8 240 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOMU.DEC"));
241 break;
da0e9ce3 242 case kSemiMuonic:
3a6cdcc8 243 SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BANDCTOMU.DEC"));
244 break;
da0e9ce3 245 case kDiElectron:
3a6cdcc8 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;
da0e9ce3 251 case kBPsiPrimeDiMuon:
3a6cdcc8 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;
de323cb7 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;
3a6cdcc8 272 case kHardMuons:
de323cb7 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;
da0e9ce3 287 case kPiToMu:
288 case kKaToMu:
289 case kAllMuonic:
290 case kWToMuon:
291 case kWToCharm:
292 case kWToCharmToMuon:
293 case kZDiMuon:
294 case kZDiElectron:
da0e9ce3 295 case kNoDecay:
296 case kNoDecayHeavy:
297 case kNeutralPion:
f2f8bbb9 298 case kBJpsiUndecayed:
49686315 299 case kNoDecayBeauty:
f2f8bbb9 300 AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
da0e9ce3 301 break;
302 }
303 ReadDecayTable();
304 }
305
306Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
307 {
308 // This method is dummy
309 return 1.;
310 }
311
312Float_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
324void 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////////////////////////////////////////////////////////////
335Bool_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