[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.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.                  *
 **************************************************************************/

/* $Id: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */

// Class to create cocktails for physics with electrons, di-electrons,
// and photons from the decay of the following sources:
// pizero, eta, rho, omega, etaprime, phi
// Kinematic distributions of the sources are taken from AliGenEMlib.
// Decay channels can be selected via the method SetDecayMode.
// Particles can be generated flat in pT with weights according to the
// chosen pT distributions from AliGenEMlib (weighting mode: kNonAnalog),
// or they are generated according to the pT distributions themselves
// (weighting mode: kAnalog)  
 
#include <TObjArray.h>
#include <TParticle.h>
#include <TF1.h>
#include <TVirtualMC.h>
#include <TPDGCode.h>
#include "AliGenCocktailEventHeader.h"

#include "AliGenCocktailEntry.h"
#include "AliGenEMCocktail.h"
#include "AliGenEMlib.h"
#include "AliGenBox.h"
#include "AliGenParam.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliStack.h"
#include "AliDecayer.h"
#include "AliDecayerPythia.h"
#include "AliLog.h"
#include "AliGenCorrHF.h"

ClassImp(AliGenEMCocktail)  
  
//________________________________________________________________________
AliGenEMCocktail::AliGenEMCocktail()
  :AliGenCocktail(),
   fDecayer(0),
   fDecayMode(kAll),
   fWeightingMode(kNonAnalog),
   fNPart(1000),
   fYieldArray()
{
// Constructor

}

//_________________________________________________________________________
AliGenEMCocktail::~AliGenEMCocktail()
{
// Destructor

}

//_________________________________________________________________________
void AliGenEMCocktail::CreateCocktail()
{
// create and add sources to the cocktail

  fDecayer->SetForceDecay(fDecayMode);
  fDecayer->ForceDecay();

// Set kinematic limits
  Double_t ptMin  = fPtMin;
  Double_t ptMax  = fPtMax;
  Double_t yMin   = fYMin;;
  Double_t yMax   = fYMax;;
  Double_t phiMin = fPhiMin*180./TMath::Pi();
  Double_t phiMax = fPhiMax*180./TMath::Pi();
  AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
  AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));

// Create and add electron sources to the generator

// pizero
  AliGenParam * genpizero=0;
  Char_t namePizero[10];    
  sprintf(namePizero,"Pizero");    
  genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
  AddSource2Generator(namePizero,genpizero);
  TF1 *fPtPizero = genpizero->GetPt();
  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
// eta  
  AliGenParam * geneta=0;
  Char_t nameEta[10];    
  sprintf(nameEta,"Eta");    
  geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
  AddSource2Generator(nameEta,geneta);
  TF1 *fPtEta = geneta->GetPt();
  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
// rho  
  AliGenParam * genrho=0;
  Char_t nameRho[10];    
  sprintf(nameRho,"Rho");    
  genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
  AddSource2Generator(nameRho,genrho);
  TF1 *fPtRho = genrho->GetPt();
  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
// omega
  AliGenParam * genomega=0;
  Char_t nameOmega[10];    
  sprintf(nameOmega,"Omega");    
  genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
  AddSource2Generator(nameOmega,genomega);
  TF1 *fPtOmega = genomega->GetPt();
  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
// etaprime
  AliGenParam * genetaprime=0;
  Char_t nameEtaprime[10];    
  sprintf(nameEtaprime,"Etaprime");    
  genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
  AddSource2Generator(nameEtaprime,genetaprime);
  TF1 *fPtEtaprime = genetaprime->GetPt();
  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
// phi  
  AliGenParam * genphi=0;
  Char_t namePhi[10];    
  sprintf(namePhi,"Phi");    
  genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
  AddSource2Generator(namePhi,genphi);
  TF1 *fPtPhi = genphi->GetPt();
  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
}

//-------------------------------------------------------------------
void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
					 AliGenParam* const genSource)
{
// add sources to the cocktail
  Double_t phiMin = fPhiMin*180./TMath::Pi();
  Double_t phiMax = fPhiMax*180./TMath::Pi();

  genSource->SetPtRange(fPtMin, fPtMax);  
  genSource->SetYRange(fYMin, fYMax);
  genSource->SetPhiRange(phiMin, phiMax);
  genSource->SetWeighting(fWeightingMode);
  if (!gMC) genSource->SetDecayer(fDecayer);  
  genSource->Init();
    
  AddGenerator(genSource,nameSource,1.); // Adding Generator    
}

//-------------------------------------------------------------------
void AliGenEMCocktail::Init()
{
// Initialisation
  TIter next(fEntries);
  AliGenCocktailEntry *entry;
  if (fStack) {
    while((entry = (AliGenCocktailEntry*)next())) {
      entry->Generator()->SetStack(fStack);
    }
  }
}

//_________________________________________________________________________
void AliGenEMCocktail::Generate()
{
// Generate event 
  TIter next(fEntries);
  AliGenCocktailEntry *entry = 0;
  AliGenCocktailEntry *preventry = 0;
  AliGenerator* gen = 0;

  if (fHeader) delete fHeader;
  fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");

  const TObjArray *partArray = gAlice->GetMCApp()->Particles();
    
// Generate the vertex position used by all generators    
  if(fVertexSmear == kPerEvent) Vertex();

//Reseting stack
  AliRunLoader * runloader = AliRunLoader::Instance();
  if (runloader)
    if (runloader->Stack())
      runloader->Stack()->Clean();
  
// Loop over generators and generate events
  Int_t igen = 0;
  const char* genName = 0;
  while((entry = (AliGenCocktailEntry*)next())) {
    gen = entry->Generator();
    genName = entry->GetName();
    gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
    
    if (fNPart > 0) {
      igen++;	
      if (igen == 1) entry->SetFirst(0);		
      else  entry->SetFirst((partArray->GetEntriesFast())+1);
      gen->SetNumberParticles(fNPart);		
      gen->Generate();
      entry->SetLast(partArray->GetEntriesFast());
      preventry = entry;
    }
  }  
  next.Reset();

// Setting weights for proper absolute normalization
  Int_t iPart, iMother;
  Int_t pdgMother = 0;
  Double_t weight = 0.;
  Double_t dNdy = 0.;
  Int_t maxPart = partArray->GetEntriesFast();
  for(iPart=0; iPart<maxPart; iPart++){      
    TParticle *part = gAlice->GetMCApp()->Particle(iPart);
    iMother = part->GetFirstMother();
    TParticle *mother = 0;
    if (iMother>=0){
      mother = gAlice->GetMCApp()->Particle(iMother);
      pdgMother = mother->GetPdgCode();
    }
    else
      pdgMother = part->GetPdgCode();
    switch (pdgMother){
    case 111:
      dNdy = fYieldArray[kGenPizero];
      break;
    case 221:
      dNdy = fYieldArray[kGenEta];
      break;
    case 113:
      dNdy = fYieldArray[kGenRho];
      break;
    case 223:
      dNdy = fYieldArray[kGenOmega];
      break;
    case 331:
      dNdy = fYieldArray[kGenEtaprime];
      break;
    case 333:
      dNdy = fYieldArray[kGenPhi];
      break;
      
    default:
      dNdy = 0.;
    }
    
    weight = dNdy*part->GetWeight();
    part->SetWeight(weight);
  }	
  
  fHeader->SetNProduced(maxPart);


  TArrayF eventVertex;
  eventVertex.Set(3);
  for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
  
  fHeader->SetPrimaryVertex(eventVertex);

  gAlice->SetGenEventHeader(fHeader);
}
Commit	Line	Data
e40b9538	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	/* $Id: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */
	17
	18	// Class to create cocktails for physics with electrons, di-electrons,
	19	// and photons from the decay of the following sources:
	20	// pizero, eta, rho, omega, etaprime, phi
	21	// Kinematic distributions of the sources are taken from AliGenEMlib.
	22	// Decay channels can be selected via the method SetDecayMode.
	23	// Particles can be generated flat in pT with weights according to the
	24	// chosen pT distributions from AliGenEMlib (weighting mode: kNonAnalog),
	25	// or they are generated according to the pT distributions themselves
	26	// (weighting mode: kAnalog)
	27
	28	#include <TObjArray.h>
	29	#include <TParticle.h>
	30	#include <TF1.h>
	31	#include <TVirtualMC.h>
	32	#include <TPDGCode.h>
	33	#include "AliGenCocktailEventHeader.h"
	34
	35	#include "AliGenCocktailEntry.h"
	36	#include "AliGenEMCocktail.h"
	37	#include "AliGenEMlib.h"
	38	#include "AliGenBox.h"
	39	#include "AliGenParam.h"
	40	#include "AliMC.h"
	41	#include "AliRun.h"
	42	#include "AliStack.h"
	43	#include "AliDecayer.h"
	44	#include "AliDecayerPythia.h"
	45	#include "AliLog.h"
	46	#include "AliGenCorrHF.h"
	47
	48	ClassImp(AliGenEMCocktail)
	49
	50	//________________________________________________________________________
	51	AliGenEMCocktail::AliGenEMCocktail()
	52	:AliGenCocktail(),
	53	fDecayer(0),
	54	fDecayMode(kAll),
	55	fWeightingMode(kNonAnalog),
	56	fNPart(1000),
	57	fYieldArray()
	58	{
	59	// Constructor
	60
	61	}
	62
	63	//_________________________________________________________________________
	64	AliGenEMCocktail::~AliGenEMCocktail()
65	{
66	// Destructor
67
68	}
69
70	//_________________________________________________________________________
71	void AliGenEMCocktail::CreateCocktail()
72	{
73	// create and add sources to the cocktail
74
75	fDecayer->SetForceDecay(fDecayMode);
76	fDecayer->ForceDecay();
77
78	// Set kinematic limits
79	Double_t ptMin = fPtMin;
80	Double_t ptMax = fPtMax;
81	Double_t yMin = fYMin;;
82	Double_t yMax = fYMax;;
83	Double_t phiMin = fPhiMin*180./TMath::Pi();
84	Double_t phiMax = fPhiMax*180./TMath::Pi();
85	AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
86	AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
87
88	// Create and add electron sources to the generator
89
90	// pizero
91	AliGenParam * genpizero=0;
92	Char_t namePizero[10];
93	sprintf(namePizero,"Pizero");
94	genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
95	AddSource2Generator(namePizero,genpizero);
96	TF1 *fPtPizero = genpizero->GetPt();
97	fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
98	// eta
99	AliGenParam * geneta=0;
100	Char_t nameEta[10];
101	sprintf(nameEta,"Eta");
102	geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
103	AddSource2Generator(nameEta,geneta);
104	TF1 *fPtEta = geneta->GetPt();
105	fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
106	// rho
107	AliGenParam * genrho=0;
108	Char_t nameRho[10];
109	sprintf(nameRho,"Rho");
110	genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
111	AddSource2Generator(nameRho,genrho);
112	TF1 *fPtRho = genrho->GetPt();
113	fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
114	// omega
115	AliGenParam * genomega=0;
116	Char_t nameOmega[10];
117	sprintf(nameOmega,"Omega");
118	genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
119	AddSource2Generator(nameOmega,genomega);
120	TF1 *fPtOmega = genomega->GetPt();
121	fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
122	// etaprime
123	AliGenParam * genetaprime=0;
124	Char_t nameEtaprime[10];
125	sprintf(nameEtaprime,"Etaprime");
126	genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
127	AddSource2Generator(nameEtaprime,genetaprime);
128	TF1 *fPtEtaprime = genetaprime->GetPt();
129	fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
130	// phi
131	AliGenParam * genphi=0;
132	Char_t namePhi[10];
133	sprintf(namePhi,"Phi");
134	genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
135	AddSource2Generator(namePhi,genphi);
136	TF1 *fPtPhi = genphi->GetPt();
137	fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
138	}
139
140	//-------------------------------------------------------------------
141	void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
142	AliGenParam* const genSource)
143	{
144	// add sources to the cocktail
145	Double_t phiMin = fPhiMin*180./TMath::Pi();
146	Double_t phiMax = fPhiMax*180./TMath::Pi();
147
148	genSource->SetPtRange(fPtMin, fPtMax);
149	genSource->SetYRange(fYMin, fYMax);
150	genSource->SetPhiRange(phiMin, phiMax);
151	genSource->SetWeighting(fWeightingMode);
152	if (!gMC) genSource->SetDecayer(fDecayer);
153	genSource->Init();
154
155	AddGenerator(genSource,nameSource,1.); // Adding Generator
156	}
157
158	//-------------------------------------------------------------------
159	void AliGenEMCocktail::Init()
160	{
161	// Initialisation
162	TIter next(fEntries);
163	AliGenCocktailEntry *entry;
164	if (fStack) {
165	while((entry = (AliGenCocktailEntry*)next())) {
166	entry->Generator()->SetStack(fStack);
167	}
168	}
169	}
170
171	//_________________________________________________________________________
172	void AliGenEMCocktail::Generate()
173	{
174	// Generate event
175	TIter next(fEntries);
176	AliGenCocktailEntry *entry = 0;
177	AliGenCocktailEntry *preventry = 0;
178	AliGenerator* gen = 0;
179
180	if (fHeader) delete fHeader;
181	fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
182
183	const TObjArray *partArray = gAlice->GetMCApp()->Particles();
184
185	// Generate the vertex position used by all generators
186	if(fVertexSmear == kPerEvent) Vertex();
187
188	//Reseting stack
189	AliRunLoader * runloader = AliRunLoader::Instance();
190	if (runloader)
191	if (runloader->Stack())
192	runloader->Stack()->Clean();
193
194	// Loop over generators and generate events
195	Int_t igen = 0;
196	const char* genName = 0;
197	while((entry = (AliGenCocktailEntry*)next())) {
198	gen = entry->Generator();
199	genName = entry->GetName();
200	gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
201
202	if (fNPart > 0) {
203	igen++;
204	if (igen == 1) entry->SetFirst(0);
205	else entry->SetFirst((partArray->GetEntriesFast())+1);
206	gen->SetNumberParticles(fNPart);
207	gen->Generate();
208	entry->SetLast(partArray->GetEntriesFast());
209	preventry = entry;
210	}
211	}
212	next.Reset();
213
214	// Setting weights for proper absolute normalization
215	Int_t iPart, iMother;
216	Int_t pdgMother = 0;
217	Double_t weight = 0.;
218	Double_t dNdy = 0.;
219	Int_t maxPart = partArray->GetEntriesFast();
220	for(iPart=0; iPart<maxPart; iPart++){
221	TParticle *part = gAlice->GetMCApp()->Particle(iPart);
222	iMother = part->GetFirstMother();
223	TParticle *mother = 0;
224	if (iMother>=0){
225	mother = gAlice->GetMCApp()->Particle(iMother);
226	pdgMother = mother->GetPdgCode();
227	}
228	else
229	pdgMother = part->GetPdgCode();
230	switch (pdgMother){
231	case 111:
232	dNdy = fYieldArray[kGenPizero];
233	break;
234	case 221:
235	dNdy = fYieldArray[kGenEta];
236	break;
237	case 113:
238	dNdy = fYieldArray[kGenRho];
239	break;
240	case 223:
241	dNdy = fYieldArray[kGenOmega];
242	break;
243	case 331:
244	dNdy = fYieldArray[kGenEtaprime];
245	break;
246	case 333:
247	dNdy = fYieldArray[kGenPhi];
248	break;
249
250	default:
251	dNdy = 0.;
252	}
253
254	weight = dNdy*part->GetWeight();
255	part->SetWeight(weight);
256	}
257
258	fHeader->SetNProduced(maxPart);
259
260
261	TArrayF eventVertex;
262	eventVertex.Set(3);
263	for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
264
265	fHeader->SetPrimaryVertex(eventVertex);
266
267	gAlice->SetGenEventHeader(fHeader);
268	}
269
270