[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(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);
  const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
  const Float_t kconvL=1./10; // mm to cm conversion
  x=x4.get(1)*kconvL;//[cm]
  y=x4.get(2)*kconvL;//[cm]
  z=x4.get(3)*kconvL;//[cm]
  t=x4.get(0)*kconvT;//[s]

  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 kHadronicD:
      SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOD.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 kHadronicDWithout4Bodies:
     case kPhiKK:
     case kOmega:
     case kLambda:	 
     case kNoDecay:
     case kNoDecayHeavy:
     case kNeutralPion:
     case kBPsiPrimeDiElectron:
     case kBeautyUpgrade:
     case kBJpsiUndecayed: 
     case kDiElectronEM:
     case kElectronEM:
     case kGammaEM:
     case kJpsiDiMuon:
     case kNoDecayBeauty:
     case kPsiPrimeJpsiDiElectron:
      AliWarning(Form("Warning: case %d not implemented for this class!",(int)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);
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	//____________________________________________________________
	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);
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
	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"));
941c36b4	210	break;
	211	case kHadronicD:
	212	SetDecayTablePath(gSystem->ExpandPathName("$ALICE_ROOT/TEvtGen/EvtGen/DecayTable/BTOD.DEC"));
	213	break;
da0e9ce3	214	case kHardMuons:
	215	case kChiToJpsiGammaToMuonMuon:
	216	case kChiToJpsiGammaToElectronElectron:
	217	case kBSemiMuonic:
	218	case kSemiMuonic:
	219	case kDiMuon:
	220	case kSemiElectronic:
	221	case kDiElectron:
	222	case kBPsiPrimeDiMuon:
	223	case kPiToMu:
	224	case kKaToMu:
	225	case kAllMuonic:
	226	case kWToMuon:
	227	case kWToCharm:
	228	case kWToCharmToMuon:
	229	case kZDiMuon:
	230	case kZDiElectron:
da0e9ce3	231	case kHadronicDWithout4Bodies:
	232	case kPhiKK:
	233	case kOmega:
	234	case kLambda:
	235	case kNoDecay:
	236	case kNoDecayHeavy:
	237	case kNeutralPion:
	238	case kBPsiPrimeDiElectron:
d77a8584	239	case kBeautyUpgrade:
f2f8bbb9	240	case kBJpsiUndecayed:
49686315	241	case kDiElectronEM:
	242	case kElectronEM:
	243	case kGammaEM:
	244	case kJpsiDiMuon:
	245	case kNoDecayBeauty:
	246	case kPsiPrimeJpsiDiElectron:
f2f8bbb9	247	AliWarning(Form("Warning: case %d not implemented for this class!",(int)decay));
da0e9ce3	248	break;
	249	}
	250	ReadDecayTable();
	251	}
	252
	253	Float_t AliDecayerEvtGen::GetPartialBranchingRatio(Int_t)
	254	{
	255	// This method is dummy
	256	return 1.;
	257	}
	258
	259	Float_t AliDecayerEvtGen::GetLifetime(Int_t kf)
	260	{
	261	//
	262	//Input: pdg code of a particle
	263	//return lifetime in sec for a particle with particle code kf
	264	//
	265	EvtId IdPart=EvtPDL::evtIdFromStdHep(kf);
	266	Double_t lifetime = EvtPDL::getctau(IdPart); //c*tau (mm)
	267	AliDebug(1,Form("lifetime is %f (mum) , particle id= %d",lifetime*1000,kf));
	268	return lifetime*kconv; //tau (sec)
	269	}
	270
	271	void AliDecayerEvtGen::ReadDecayTable()
	272	{
	273	//Input none - Output none
	274	//Read the decay table that correspond to the path
	275	//fDecayTablePath
	276	//
	277	TString temp = fDecayTablePath;
	278	if(!temp.EndsWith("DECAY.DEC"))
	279	fGenerator->readUDecay(fDecayTablePath);
	280	}
	281	////////////////////////////////////////////////////////////
	282	Bool_t AliDecayerEvtGen::SetDecayTablePath(Char_t *path)
	283	{
	284	//
	285	//Set the path of the decay table read to force particle decays
	286	//
	287	if(gSystem->AccessPathName(path))
	288	{
	289	AliWarning("Attention: This path not exist!\n");
	290	return kFALSE;
	291	}
	292	fDecayTablePath = path;
	293	return kTRUE;
	294	}
	295
	296
	297