Adding STEER/STEER to EINCLUDE
[u/mrichter/AliRoot.git] / 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 kHardMuons:
227      case kChiToJpsiGammaToMuonMuon:
228      case kChiToJpsiGammaToElectronElectron:
229      case kBSemiMuonic:
230      case kSemiMuonic:
231      case kDiMuon:
232      case kSemiElectronic:
233      case kDiElectron:
234      case kBPsiPrimeDiMuon:
235      case kPiToMu:
236      case kKaToMu:
237      case kAllMuonic:
238      case kWToMuon:
239      case kWToCharm:
240      case kWToCharmToMuon:
241      case kZDiMuon:
242      case kZDiElectron:
243      case kHadronicDWithout4Bodies:
244      case kPhiKK:
245      case kOmega:
246      case kLambda:
247      case kNoDecay:
248      case kNoDecayHeavy:
249      case kNeutralPion:
250      case kBPsiPrimeDiElectron:
251      case kBeautyUpgrade:
252      case kBJpsiUndecayed: 
253      case kDiElectronEM:
254      case kElectronEM:
255      case kGammaEM:
256      case kJpsiDiMuon:
257      case kNoDecayBeauty:
258      case kPsiPrimeJpsiDiElectron:
259       AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
260      break;
261      }
262      ReadDecayTable();
263   }
264
265 Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
266   {
267    // This method is dummy
268    return  1.;
269   }
270
271 Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
272   {    
273   //
274   //Input: pdg code of a particle   
275   //return lifetime in sec for a particle with particle code kf
276   //
277   EvtId IdPart=EvtPDL::evtIdFromStdHep(kf); 
278   Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
279   AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
280   return lifetime*kconv; //tau (sec)
281   }
282
283 void AliDecayerEvtGen::ReadDecayTable()
284     {
285      //Input none - Output none 
286      //Read the decay table that correspond to the path 
287      //fDecayTablePath
288      // 
289      TString temp = fDecayTablePath;
290      if(!temp.EndsWith("DECAY.DEC"))
291      fGenerator->readUDecay(fDecayTablePath);
292     }
293 ////////////////////////////////////////////////////////////
294 Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
295   {
296   // 
297   //Set the path of the decay table read to force particle decays 
298   //
299   if(gSystem->AccessPathName(path)) 
300   {
301    AliWarning("Attention: This path not exist!\n");
302    return kFALSE;
303   }
304   fDecayTablePath = path;
305   return kTRUE;
306   } 
307
308
309