Changes needed for MC production of J/psi from B with radiative decays
[u/mrichter/AliRoot.git] / 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"
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
37ClassImp(AliDecayerEvtGen)
38//____________________________________________________________
39AliDecayerEvtGen::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//_________________________________________________________________
53AliDecayerEvtGen::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//____________________________________________________________
66AliDecayerEvtGen::~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//___________________________________________________________
82void 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//____________________________________________________________
96void 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);
34dbffe9 104 EvtVector4R p_init(p->E(),p->Px(),p->Py(),p->Pz());
da0e9ce3 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//____________________________________________________________
114Int_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);
774ceaaf 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]
da0e9ce3 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
179void AliDecayerEvtGen::Copy(TObject &) const
180 {
181 //
182 // Copy *this onto AliDecayerEvtGen -- not implemented
183 //
184 Fatal("Copy","Not implemented!\n");
185 }
186
187void 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:
d77a8584 237 case kBeautyUpgrade:
f2f8bbb9 238 case kBJpsiUndecayed:
239 AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
da0e9ce3 240 break;
241 }
242 ReadDecayTable();
243 }
244
245Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
246 {
247 // This method is dummy
248 return 1.;
249 }
250
251Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
252 {
253 //
254 //Input: pdg code of a particle
255 //return lifetime in sec for a particle with particle code kf
256 //
257 EvtId IdPart=EvtPDL::evtIdFromStdHep(kf);
258 Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
259 AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
260 return lifetime*kconv; //tau (sec)
261 }
262
263void AliDecayerEvtGen::ReadDecayTable()
264 {
265 //Input none - Output none
266 //Read the decay table that correspond to the path
267 //fDecayTablePath
268 //
269 TString temp = fDecayTablePath;
270 if(!temp.EndsWith("DECAY.DEC"))
271 fGenerator->readUDecay(fDecayTablePath);
272 }
273////////////////////////////////////////////////////////////
274Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
275 {
276 //
277 //Set the path of the decay table read to force particle decays
278 //
279 if(gSystem->AccessPathName(path))
280 {
281 AliWarning("Attention: This path not exist!\n");
282 return kFALSE;
283 }
284 fDecayTablePath = path;
285 return kTRUE;
286 }
287
288
289