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