[u/mrichter/AliRoot.git] / TEvtGen / AliDecayerEvtGen.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
///////////////////////////////////////////////////////////////////////////
//  Implementation of AliDecayer using EvtGen package.                   //
//                                                                       //
//  Giuseppe E. Bruno            &    Fiorella Fionda                    //
//  (Giuseppe.Bruno@ba.infn.it)       (Fiorella.Fionda@ba.infn.it)       //
///////////////////////////////////////////////////////////////////////////

#include <TSystem.h>
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TLorentzVector.h>

#include "EvtGenBase/EvtStdHep.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "AliDecayerEvtGen.h"
#include "AliLog.h"

ClassImp(AliDecayerEvtGen)
//____________________________________________________________
AliDecayerEvtGen::AliDecayerEvtGen():
  fRandomEngine(0x0),
  fGenerator(0x0),
  fEvtstdhep(0x0),
  fDecayTablePath(0x0),
  fParticleTablePath(0x0),
  fDecay(kAll)
  {
  // Default constructor
  fEvtstdhep = new EvtStdHep();
  fDecayTablePath = gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DECAY.DEC"); //default decay table
  fParticleTablePath = gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/evt.pdl"); //particle table
  }
//_________________________________________________________________
AliDecayerEvtGen::AliDecayerEvtGen(const AliDecayerEvtGen &decayer):
  AliDecayer(decayer),
  fRandomEngine(decayer.fRandomEngine),
  fGenerator(decayer.fGenerator),
  fEvtstdhep(decayer.fEvtstdhep),
  fDecayTablePath(decayer.fDecayTablePath),
  fParticleTablePath(decayer.fParticleTablePath),
  fDecay(decayer.fDecay)
  {
  // Copy Constructor
  decayer.Copy(*this);
  }
//____________________________________________________________
AliDecayerEvtGen::~AliDecayerEvtGen()
  {
  // Destructor
  if(fRandomEngine) {delete fRandomEngine;}
  fRandomEngine = 0;
  if(fGenerator) {delete fGenerator;}
  fGenerator = 0;
  if(fEvtstdhep) {delete fEvtstdhep;}
  fEvtstdhep = 0;
  if(fDecayTablePath) {delete fDecayTablePath;}
  fDecayTablePath = 0;
  if(fParticleTablePath) {delete fParticleTablePath;}
  fParticleTablePath = 0;
  }

//___________________________________________________________
void AliDecayerEvtGen::Init()
  {
  //Standard AliDecayerEvtGen initializer:
  //initialize EvtGen with default decay table (DECAY.DEC), particle table (evt.pdl)
  //and fRandomEngine for generation of random numbers
  //
  if(fGenerator){
  AliWarning(" AliDecayerEvtGen already initialized!!!!\n");
  return;
  }
  fRandomEngine=new EvtNUMRandomEngine();
  fGenerator=new EvtGen(fDecayTablePath,fParticleTablePath,fRandomEngine);
  }
//____________________________________________________________
void AliDecayerEvtGen::Decay(Int_t ipart, TLorentzVector *p)
  {
  //
  //Decay a particle
  //input: pdg code and momentum of the particle to be decayed  
  //all informations about decay products are stored in fEvtstdhep 
  //
  EvtId IPART=EvtPDL::evtIdFromStdHep(ipart);
  EvtVector4R p_init=EvtVector4R::EvtVector4R(p->E(),p->Px(),p->Py(),p->Pz());
  EvtParticle *froot_part=EvtParticleFactory::particleFactory(IPART,p_init);
  fGenerator->generateDecay(froot_part);
  fEvtstdhep->init();
  froot_part->makeStdHep(*fEvtstdhep);
  //froot_part->printTree(); //to print the decay chain 
  froot_part->deleteTree();
  }

//____________________________________________________________
Int_t AliDecayerEvtGen::ImportParticles(TClonesArray *particles)
  {
  //
  //Input: pointer to a TClonesArray - Output(Int_t): number of decay products    
  //Put all the informations about the decay products in the 
  //TClonesArray particles
  //
  if (particles == 0) return 0;
  TClonesArray &clonesParticles = *particles;
  clonesParticles.Clear();

  int j;
  int istat;
  int partnum;
  double px,py,pz,e;
  double x,y,z,t;
  EvtVector4R p4,x4;

  Int_t npart=fEvtstdhep->getNPart();
  for(int i=0;i<fEvtstdhep->getNPart();i++){
  j=i+1;
  int jmotherfirst=fEvtstdhep->getFirstMother(i)+1;
  int jmotherlast=fEvtstdhep->getLastMother(i)+1;
  int jdaugfirst=fEvtstdhep->getFirstDaughter(i)+1;
  int jdauglast=fEvtstdhep->getLastDaughter(i)+1;

  partnum=fEvtstdhep->getStdHepID(i);
  
  //verify if all particles of decay chain are in the TDatabasePDG
  TParticlePDG *partPDG = TDatabasePDG::Instance()->GetParticle(partnum);
  if(!partPDG)
    {
    AliWarning("Particle code non known in TDatabasePDG - set pdg = 89");
    partnum=89; //internal use for unspecified resonance data
    }

  istat=fEvtstdhep->getIStat(i);

  if(istat!=1 && istat!=2) Info("ImportParticles","Attention: unknown status code!");
  if(istat == 2) istat = 11; //status decayed

  p4=fEvtstdhep->getP4(i);
  x4=fEvtstdhep->getX4(i);
  px=p4.get(1);
  py=p4.get(2);
  pz=p4.get(3);
  e=p4.get(0);

  x=x4.get(1);//[mm]
  y=x4.get(2);//[mm]
  z=x4.get(3);//[mm]
  t=x4.get(0);//[mm]

  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));

  new(clonesParticles[i]) TParticle(partnum,istat,jmotherfirst,-1,jdaugfirst,jdauglast,px,py,pz,e,x,y,z,t);

  //set polarization!!!
  }

  return npart; 

  }

void  AliDecayerEvtGen::Copy(TObject &) const
  {
  //
  // Copy *this onto AliDecayerEvtGen -- not implemented
  //
   Fatal("Copy","Not implemented!\n");
  }

void AliDecayerEvtGen::ForceDecay()
  {
  //
  // Intupt: none - Output: none
  // Set the decay mode to decay particles: for each case is read a 
  // different decay table. case kAll read the default decay table only   
  //
  Decay_t decay = fDecay;
  switch(decay)
    {
     case kAll: 
     break;
     case kBJpsiDiElectron:
     SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOELE.DEC"));
     break;
     case kBJpsi:
     SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSI.DEC"));
     break;
     case kBJpsiDiMuon:
     SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOMU.DEC"));
     break;
     case kBSemiElectronic:
     SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOELE.DEC"));
     break;  
     case kHardMuons:
     case kChiToJpsiGammaToMuonMuon:
     case kChiToJpsiGammaToElectronElectron:
     case kBSemiMuonic:
     case kSemiMuonic:
     case kDiMuon:
     case kSemiElectronic:
     case kDiElectron:
     case kBPsiPrimeDiMuon:
     case kPiToMu:
     case kKaToMu:
     case kAllMuonic:
     case kWToMuon:
     case kWToCharm:
     case kWToCharmToMuon:
     case kZDiMuon:
     case kZDiElectron:
     case kHadronicD:
     case kHadronicDWithout4Bodies:
     case kPhiKK:
     case kOmega:
     case kLambda:	 
     case kNoDecay:
     case kNoDecayHeavy:
     case kNeutralPion:
     case kBPsiPrimeDiElectron:
     case kLambdaBLambda:
     AliWarning(Form("Warning: case %s not implemented for this class!",(Char_t)decay));
     break;
     }
     ReadDecayTable();
  }

Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
  {
   // This method is dummy
   return  1.;
  }

Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
  {    
  //
  //Input: pdg code of a particle   
  //return lifetime in sec for a particle with particle code kf
  //
  EvtId IdPart=EvtPDL::evtIdFromStdHep(kf); 
  Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
  AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
  return lifetime*kconv; //tau (sec)
  }

void AliDecayerEvtGen::ReadDecayTable()
    {
     //Input none - Output none 
     //Read the decay table that correspond to the path 
     //fDecayTablePath
     // 
     TString temp = fDecayTablePath;
     if(!temp.EndsWith("DECAY.DEC"))
     fGenerator->readUDecay(fDecayTablePath);
    }
////////////////////////////////////////////////////////////
Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
  {
  // 
  //Set the path of the decay table read to force particle decays 
  //
  if(gSystem->AccessPathName(path)) 
  {
   AliWarning("Attention: This path not exist!\n");
   return kFALSE;
  }
  fDecayTablePath = path;
  return kTRUE;
  }
Commit	Line	Data
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
	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=EvtVector4R::EvtVector4R(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
162	x=x4.get(1);//[mm]
163	y=x4.get(2);//[mm]
164	z=x4.get(3);//[mm]
165	t=x4.get(0);//[mm]
166
167	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));
168
169	new(clonesParticles[i]) TParticle(partnum,istat,jmotherfirst,-1,jdaugfirst,jdauglast,px,py,pz,e,x,y,z,t);
170
171	//set polarization!!!
172	}
173
174	return npart;
175
176	}
177
178	void AliDecayerEvtGen::Copy(TObject &) const
179	{
180	//
181	// Copy *this onto AliDecayerEvtGen -- not implemented
182	//
183	Fatal("Copy","Not implemented!\n");
184	}
185
186	void AliDecayerEvtGen::ForceDecay()
187	{
188	//
189	// Intupt: none - Output: none
190	// Set the decay mode to decay particles: for each case is read a
191	// different decay table. case kAll read the default decay table only
192	//
193	Decay_t decay = fDecay;
194	switch(decay)
195	{
196	case kAll:
197	break;
198	case kBJpsiDiElectron:
199	SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOELE.DEC"));
200	break;
201	case kBJpsi:
202	SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSI.DEC"));
203	break;
204	case kBJpsiDiMuon:
205	SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOJPSITOMU.DEC"));
206	break;
207	case kBSemiElectronic:
208	SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOELE.DEC"));
209	break;
210	case kHardMuons:
211	case kChiToJpsiGammaToMuonMuon:
212	case kChiToJpsiGammaToElectronElectron:
213	case kBSemiMuonic:
214	case kSemiMuonic:
215	case kDiMuon:
216	case kSemiElectronic:
217	case kDiElectron:
218	case kBPsiPrimeDiMuon:
219	case kPiToMu:
220	case kKaToMu:
221	case kAllMuonic:
222	case kWToMuon:
223	case kWToCharm:
224	case kWToCharmToMuon:
225	case kZDiMuon:
226	case kZDiElectron:
227	case kHadronicD:
228	case kHadronicDWithout4Bodies:
229	case kPhiKK:
230	case kOmega:
231	case kLambda:
232	case kNoDecay:
233	case kNoDecayHeavy:
234	case kNeutralPion:
235	case kBPsiPrimeDiElectron:
e549a63e	236	case kLambdaBLambda:
da0e9ce3	237	AliWarning(Form("Warning: case %s not implemented for this class!",(Char_t)decay));
	238	break;
	239	}
	240	ReadDecayTable();
	241	}
	242
	243	Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
	244	{
	245	// This method is dummy
	246	return 1.;
	247	}
	248
	249	Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
	250	{
	251	//
	252	//Input: pdg code of a particle
	253	//return lifetime in sec for a particle with particle code kf
	254	//
	255	EvtId IdPart=EvtPDL::evtIdFromStdHep(kf);
	256	Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
	257	AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
	258	return lifetime*kconv; //tau (sec)
	259	}
	260
	261	void AliDecayerEvtGen::ReadDecayTable()
	262	{
	263	//Input none - Output none
	264	//Read the decay table that correspond to the path
	265	//fDecayTablePath
	266	//
	267	TString temp = fDecayTablePath;
	268	if(!temp.EndsWith("DECAY.DEC"))
	269	fGenerator->readUDecay(fDecayTablePath);
	270	}
	271	////////////////////////////////////////////////////////////
	272	Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
	273	{
	274	//
	275	//Set the path of the decay table read to force particle decays
	276	//
	277	if(gSystem->AccessPathName(path))
	278	{
	279	AliWarning("Attention: This path not exist!\n");
	280	return kFALSE;
	281	}
	282	fDecayTablePath = path;
	283	return kTRUE;
	284	}
	285
	286
	287