[u/mrichter/AliRoot.git] / TEvtGen / AliGenEvtGen.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.                  *
***************************************************************************/
///////////////////////////////////////////////////////////////////////////
//  Class to generate decays of particles generated by a                 //
//  previous generator. It works as a generator, but pratically it       //
//  performs only decays. It works with this scheme: first loops over    // 
//  particles on the stack, selects those to be decayed, decays them     //
//  and then pushes the decay products on the stack.                     //
//                                                                       //
//  Giuseppe E. Bruno            &    Fiorella Fionda                    //
//  (Giuseppe.Bruno@ba.infn.it)       (Fiorella.Fionda@ba.infn.it)       //
///////////////////////////////////////////////////////////////////////////

#include "AliStack.h"
#include "AliGenEvtGen.h"
#include "AliRun.h"
#include "AliLog.h"
#include <TParticle.h>

ClassImp(AliGenEvtGen)
///////////////////////////////////////////////////////////////////////////
AliGenEvtGen::AliGenEvtGen(): 
  fStack(0x0),    
  fDecayer(0x0),
  fForceDecay(kAll),
  fSwitchOff(kBeautyPart),
  fUserDecay(kFALSE),
  fUserDecayTablePath(0x0)
  {
  //
  // Default Construction
  //
  }
///////////////////////////////////////////////////////////////////////////////////////////
AliGenEvtGen::~AliGenEvtGen() 
  {
  //
  // Standard Destructor
  //
  if(fStack) {delete fStack;}
  fStack = 0; 
  if(fDecayer) {delete fDecayer;}
  fDecayer = 0;
  if(fUserDecayTablePath) {delete fUserDecayTablePath;}
  fUserDecayTablePath = 0;
  }
///////////////////////////////////////////////////////////////////////////////////////////

void AliGenEvtGen::Init() 
  {
  //
  // Standard AliGenerator Initializer - no input 
  // 1) initialize EvtGen with default decay and particle table
  // 2) set the decay mode to force particle 
  // 3) set a user decay table if defined
  //
  if(fDecayer)
  {
  AliWarning("AliGenEvtGen already initialized!!!");
  return;
  }
  fDecayer = new AliDecayerEvtGen();
  fDecayer->Init(); //read the default decay table DECAY.DEC and particle table

  //if is set a decay mode: default decay mode is kAll  
  fDecayer->SetForceDecay(fForceDecay);
  fDecayer->ForceDecay();

  //if is defined a user decay table
  if(fUserDecay) 
    {
    fDecayer->SetDecayTablePath(fUserDecayTablePath);
    fDecayer->ReadDecayTable();
    }
  }

/////////////////////////////////////////////////////////////////////////////////////////////
void AliGenEvtGen::Generate() 
  {
  //
  //Generate method - Input none - Output none
  //For each event:
  //1)return the stack of the previous generator and select particles to be decayed by EvtGen
  //2)decay particles selected and put the decay products on the stack
  //
  //
  Float_t polar[3]= {0,0,0};  // Polarisation of daughter particles 
  Float_t origin0[3];         // Origin of the parent particle 
  Float_t pc[3], och[3]; // Momentum and origin of the children particles from EvtGen
  Int_t nt;
  Float_t tof;
  Int_t nPrimsPythia;
  TLorentzVector *mom=new TLorentzVector();
  static TClonesArray *particles;
  if(!particles) particles = new TClonesArray("TParticle",1000);
  fStack = AliRunLoader::Instance()->Stack();
  if(!fStack) {Info("Generate","Error: No stack found!"); return;}
  nPrimsPythia = fStack->GetNprimary();  
  AliDebug(1,Form("nPrimsPythia = %d \n",nPrimsPythia));
  for (Int_t iTrack = 0; iTrack < nPrimsPythia; ++iTrack) {
  TParticle *part = fStack->Particle(iTrack);
  Int_t pdg=part->GetPdgCode();

  AliDebug(1,Form("GetFlavour = %d e pdg = %d \n",GetFlavour(pdg),pdg));
  
  switch(fSwitchOff)
  {
  case kAllPart: 
  break;
  case kBeautyPart:
  if(GetFlavour(pdg)!=5) continue;
  break; 
  case kCharmPart: 
  if(GetFlavour(pdg)!=4) continue;
  break;
  } 

  //check if particle is already decayed by Pythia  
  if(part->GetStatusCode() != 1 || part->GetNDaughters()>0) 
    {
    Info("AliGenEvtGen","Attention: particle %d is already decayed by Pythia!",pdg); 
    continue;
    }
  
  part->SetStatusCode(11); //Set particle as decayed : change the status code
  
  mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
  Int_t np;

   do{
     fDecayer->Decay(part->GetPdgCode(),mom);
     np = fDecayer->ImportParticles(particles);
     }while(np<0);

  Int_t* trackIt    = new Int_t[np];
  Int_t* pParent   = new Int_t[np];
  AliDebug(1,Form("np = %d \n",np)); 
  
  for (int i = 0; i < np; i++) {
             pParent[i]   = -1;
             trackIt[i]   =  0;
         }
  //select trackable particle  
  if (np >1) {
  TParticle* iparticle =  (TParticle *) particles->At(0);//parent particle
   for (int i = 1; i<np ; i++) { 
     iparticle = (TParticle*) particles->At(i);
     Int_t ks = iparticle->GetStatusCode();

     //track last decay products
     if(ks==1) trackIt[i]=1; 
   
     }//decay particles loop

  }// if decay products   
 
  origin0[0]=part->Vx(); //[cm]
  origin0[1]=part->Vy(); //[cm]
  origin0[2]=part->Vz(); //[cm]
  //
  // Put decay products on the stack
  //
  for (int i = 1; i < np; i++) {
         TParticle* iparticle = (TParticle *) particles->At(i);
         Int_t kf   = iparticle->GetPdgCode();
         Int_t ksc  = iparticle->GetStatusCode();
         Int_t jpa  = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
         Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;      

         och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
         och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
         och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
         pc[0]  = iparticle->Px(); //[GeV/c]
         pc[1]  = iparticle->Py(); //[GeV/c]
         pc[2]  = iparticle->Pz(); //[GeV/c]
         tof = part->T()+kconv*iparticle->T(); 

         AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));         
          
         PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc); 
         if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
         pParent[i] = nt;
         KeepTrack(nt);
         SetHighWaterMark(nt);
         }// Particle loop
         particles->Clear();
     	 if (trackIt)    delete[] trackIt;
         if (pParent)    delete[] pParent;
       }
  Info("Generate","AliGenEvtGen DONE");  
  }

//////////////////////////////////////////////////////////////////////////////////////////
Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)  
  { 
  //
  // return the flavour of a particle
  // input: pdg code of the particle
  // output: Int_t 
  //         3 in case of strange (open and hidden)
  //         4 in case of charm (")
  //         5 in case of beauty (")
  //
  Int_t pdg = TMath::Abs(pdgCode);
  //Resonance
  if (pdg > 100000) pdg %= 100000;
  if(pdg > 10000)  pdg %= 10000;
  // meson ?
  if(pdg > 10) pdg/=100;
  // baryon ?
  if(pdg > 10) pdg/=10;
  return pdg;
  }

/////////////////////////////////////////////////////////////////////////////////////////
Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
  {
  //
  //Set the path of user decay table if it is defined
  //
  //put a comment to control if path exists 
  if(gSystem->AccessPathName(path))
    {
    AliWarning("Attention: This path not exist!\n");
    return kFALSE;
    }
  fUserDecayTablePath = path;
  fUserDecay = kTRUE;
  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	// Class to generate decays of particles generated by a //
	17	// previous generator. It works as a generator, but pratically it //
	18	// performs only decays. It works with this scheme: first loops over //
	19	// particles on the stack, selects those to be decayed, decays them //
	20	// and then pushes the decay products on the stack. //
	21	// //
	22	// Giuseppe E. Bruno & Fiorella Fionda //
	23	// (Giuseppe.Bruno@ba.infn.it) (Fiorella.Fionda@ba.infn.it) //
	24	///////////////////////////////////////////////////////////////////////////
	25
	26	#include "AliStack.h"
	27	#include "AliGenEvtGen.h"
	28	#include "AliRun.h"
51998ca8	29	#include "AliLog.h"
da0e9ce3	30	#include <TParticle.h>
	31
	32	ClassImp(AliGenEvtGen)
	33	///////////////////////////////////////////////////////////////////////////
	34	AliGenEvtGen::AliGenEvtGen():
	35	fStack(0x0),
	36	fDecayer(0x0),
	37	fForceDecay(kAll),
	38	fSwitchOff(kBeautyPart),
	39	fUserDecay(kFALSE),
	40	fUserDecayTablePath(0x0)
	41	{
	42	//
	43	// Default Construction
	44	//
	45	}
	46	///////////////////////////////////////////////////////////////////////////////////////////
	47	AliGenEvtGen::~AliGenEvtGen()
	48	{
	49	//
	50	// Standard Destructor
	51	//
	52	if(fStack) {delete fStack;}
	53	fStack = 0;
	54	if(fDecayer) {delete fDecayer;}
	55	fDecayer = 0;
	56	if(fUserDecayTablePath) {delete fUserDecayTablePath;}
	57	fUserDecayTablePath = 0;
	58	}
	59	///////////////////////////////////////////////////////////////////////////////////////////
	60
	61	void AliGenEvtGen::Init()
	62	{
	63	//
	64	// Standard AliGenerator Initializer - no input
	65	// 1) initialize EvtGen with default decay and particle table
	66	// 2) set the decay mode to force particle
	67	// 3) set a user decay table if defined
	68	//
	69	if(fDecayer)
	70	{
	71	AliWarning("AliGenEvtGen already initialized!!!");
	72	return;
	73	}
	74	fDecayer = new AliDecayerEvtGen();
	75	fDecayer->Init(); //read the default decay table DECAY.DEC and particle table
	76
	77	//if is set a decay mode: default decay mode is kAll
	78	fDecayer->SetForceDecay(fForceDecay);
	79	fDecayer->ForceDecay();
	80
	81	//if is defined a user decay table
	82	if(fUserDecay)
	83	{
	84	fDecayer->SetDecayTablePath(fUserDecayTablePath);
	85	fDecayer->ReadDecayTable();
	86	}
	87	}
	88
	89	/////////////////////////////////////////////////////////////////////////////////////////////
	90	void AliGenEvtGen::Generate()
	91	{
	92	//
	93	//Generate method - Input none - Output none
94	//For each event:
95	//1)return the stack of the previous generator and select particles to be decayed by EvtGen
96	//2)decay particles selected and put the decay products on the stack
97	//
98	//
99	Float_t polar[3]= {0,0,0}; // Polarisation of daughter particles
100	Float_t origin0[3]; // Origin of the parent particle
101	Float_t pc[3], och[3]; // Momentum and origin of the children particles from EvtGen
102	Int_t nt;
103	Float_t tof;
104	Int_t nPrimsPythia;
105	TLorentzVector *mom=new TLorentzVector();
106	static TClonesArray *particles;
107	if(!particles) particles = new TClonesArray("TParticle",1000);
108	fStack = AliRunLoader::Instance()->Stack();
109	if(!fStack) {Info("Generate","Error: No stack found!"); return;}
110	nPrimsPythia = fStack->GetNprimary();
111	AliDebug(1,Form("nPrimsPythia = %d \n",nPrimsPythia));
112	for (Int_t iTrack = 0; iTrack < nPrimsPythia; ++iTrack) {
113	TParticle *part = fStack->Particle(iTrack);
114	Int_t pdg=part->GetPdgCode();
115
116	AliDebug(1,Form("GetFlavour = %d e pdg = %d \n",GetFlavour(pdg),pdg));
117
118	switch(fSwitchOff)
119	{
120	case kAllPart:
121	break;
122	case kBeautyPart:
123	if(GetFlavour(pdg)!=5) continue;
124	break;
125	case kCharmPart:
126	if(GetFlavour(pdg)!=4) continue;
127	break;
128	}
129
130	//check if particle is already decayed by Pythia
131	if(part->GetStatusCode() != 1 \|\| part->GetNDaughters()>0)
132	{
133	Info("AliGenEvtGen","Attention: particle %d is already decayed by Pythia!",pdg);
134	continue;
135	}
136
137	part->SetStatusCode(11); //Set particle as decayed : change the status code
138
139	mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
140	Int_t np;
141
142	do{
143	fDecayer->Decay(part->GetPdgCode(),mom);
144	np = fDecayer->ImportParticles(particles);
145	}while(np<0);
146
147	Int_t* trackIt = new Int_t[np];
148	Int_t* pParent = new Int_t[np];
149	AliDebug(1,Form("np = %d \n",np));
150
151	for (int i = 0; i < np; i++) {
152	pParent[i] = -1;
153	trackIt[i] = 0;
154	}
155	//select trackable particle
156	if (np >1) {
157	TParticle* iparticle = (TParticle *) particles->At(0);//parent particle
158	for (int i = 1; i<np ; i++) {
159	iparticle = (TParticle*) particles->At(i);
160	Int_t ks = iparticle->GetStatusCode();
161
162	//track last decay products
163	if(ks==1) trackIt[i]=1;
164
165	}//decay particles loop
166
167	}// if decay products
168
169	origin0[0]=part->Vx(); //[cm]
170	origin0[1]=part->Vy(); //[cm]
171	origin0[2]=part->Vz(); //[cm]
172	//
173	// Put decay products on the stack
174	//
175	for (int i = 1; i < np; i++) {
176	TParticle* iparticle = (TParticle *) particles->At(i);
177	Int_t kf = iparticle->GetPdgCode();
178	Int_t ksc = iparticle->GetStatusCode();
179	Int_t jpa = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
180	Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;
181
182	och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
183	och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
184	och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
185	pc[0] = iparticle->Px(); //[GeV/c]
186	pc[1] = iparticle->Py(); //[GeV/c]
187	pc[2] = iparticle->Pz(); //[GeV/c]
188	tof = part->T()+kconv*iparticle->T();
189
190	AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));
191
192	PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc);
193	if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
194	pParent[i] = nt;
195	KeepTrack(nt);
196	SetHighWaterMark(nt);
197	}// Particle loop
198	particles->Clear();
199	if (trackIt) delete[] trackIt;
200	if (pParent) delete[] pParent;
201	}
202	Info("Generate","AliGenEvtGen DONE");
203	}
204
205	//////////////////////////////////////////////////////////////////////////////////////////
206	Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)
207	{
208	//
209	// return the flavour of a particle
210	// input: pdg code of the particle
211	// output: Int_t
212	// 3 in case of strange (open and hidden)
213	// 4 in case of charm (")
214	// 5 in case of beauty (")
215	//
216	Int_t pdg = TMath::Abs(pdgCode);
217	//Resonance
218	if (pdg > 100000) pdg %= 100000;
219	if(pdg > 10000) pdg %= 10000;
220	// meson ?
221	if(pdg > 10) pdg/=100;
222	// baryon ?
223	if(pdg > 10) pdg/=10;
224	return pdg;
225	}
226
227	/////////////////////////////////////////////////////////////////////////////////////////
228	Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
229	{
230	//
231	//Set the path of user decay table if it is defined
232	//
233	//put a comment to control if path exists
234	if(gSystem->AccessPathName(path))
235	{
236	AliWarning("Attention: This path not exist!\n");
237	return kFALSE;
238	}
239	fUserDecayTablePath = path;
240	fUserDecay = kTRUE;
241	return kTRUE;
242	}