[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 the cocktail 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(),
  fPtSelect(0),
  fCentrality(0),
  fV2Systematic(0),
  fForceConv(kFALSE),
  fHeaviestParticle(kGENs)
{
  // 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));

  //Initialize user selection for Pt Parameterization and centrality:
  AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);

  // Create and add electron sources to the generator

  // pizero
  AliGenParam *genpizero=0;
  Char_t namePizero[10];    
  snprintf(namePizero,10,"Pizero");    
  genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
  AddSource2Generator(namePizero,genpizero);
  TF1 *fPtPizero = genpizero->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

  // eta  
  if(fHeaviestParticle<kGenEta)return;
  AliGenParam *geneta=0;
  Char_t nameEta[10];    
  snprintf(nameEta,10,"Eta");    
  geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
  AddSource2Generator(nameEta,geneta);
  TF1 *fPtEta = geneta->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

  // rho  
  if(fHeaviestParticle<kGenRho)return;
  AliGenParam *genrho=0;
  Char_t nameRho[10];    
  snprintf(nameRho,10,"Rho");    
  genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
  AddSource2Generator(nameRho,genrho);
  TF1 *fPtRho = genrho->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif
  
  // omega
  if(fHeaviestParticle<kGenOmega)return;
  AliGenParam *genomega=0;
  Char_t nameOmega[10];    
  snprintf(nameOmega,10,"Omega");    
  genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
  AddSource2Generator(nameOmega,genomega);
  TF1 *fPtOmega = genomega->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

  // etaprime
  if(fHeaviestParticle<kGenEtaprime)return;
  AliGenParam *genetaprime=0;
  Char_t nameEtaprime[10];    
  snprintf(nameEtaprime,10,"Etaprime");    
  genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
  AddSource2Generator(nameEtaprime,genetaprime);
  TF1 *fPtEtaprime = genetaprime->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

  // phi  
  if(fHeaviestParticle<kGenPhi)return;
  AliGenParam *genphi=0;
  Char_t namePhi[10];    
  snprintf(namePhi,10,"Phi");    
  genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
  AddSource2Generator(namePhi,genphi);
  TF1 *fPtPhi = genphi->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

  // jpsi  
  if(fHeaviestParticle<kGenJpsi)return;
  AliGenParam *genjpsi=0;
  Char_t nameJpsi[10];    
  snprintf(nameJpsi,10,"Jpsi");    
  genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
  AddSource2Generator(nameJpsi,genjpsi);
  TF1 *fPtJpsi = genjpsi->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
  fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
#else
  fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif

}

//-------------------------------------------------------------------
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);
  genSource->SetForceGammaConversion(fForceConv);
  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;
  Float_t evPlane;
  Rndm(&evPlane,1);
  evPlane*=TMath::Pi()*2;
  while((entry = (AliGenCocktailEntry*)next())) {
    gen = entry->Generator();
    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->SetEventPlane(evPlane);
      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;
    case 443:
      dNdy = fYieldArray[kGenJpsi];
      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
4ae1c9f0	18	// Class to create the cocktail for physics with electrons, di-electrons,
e40b9538	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
4ae1c9f0	28
e40b9538	29	#include <TObjArray.h>
	30	#include <TParticle.h>
	31	#include <TF1.h>
	32	#include <TVirtualMC.h>
	33	#include <TPDGCode.h>
	34	#include "AliGenCocktailEventHeader.h"
	35
	36	#include "AliGenCocktailEntry.h"
	37	#include "AliGenEMCocktail.h"
	38	#include "AliGenEMlib.h"
	39	#include "AliGenBox.h"
	40	#include "AliGenParam.h"
	41	#include "AliMC.h"
	42	#include "AliRun.h"
	43	#include "AliStack.h"
	44	#include "AliDecayer.h"
	45	#include "AliDecayerPythia.h"
	46	#include "AliLog.h"
	47	#include "AliGenCorrHF.h"
	48
	49	ClassImp(AliGenEMCocktail)
	50
	51	//________________________________________________________________________
	52	AliGenEMCocktail::AliGenEMCocktail()
4ae1c9f0	53	:AliGenCocktail(),
e40b9538	54	fDecayer(0),
	55	fDecayMode(kAll),
	56	fWeightingMode(kNonAnalog),
	57	fNPart(1000),
4ae1c9f0	58	fYieldArray(),
	59	fPtSelect(0),
	60	fCentrality(0),
	61	fV2Systematic(0),
	62	fForceConv(kFALSE),
	63	fHeaviestParticle(kGENs)
e40b9538	64	{
4ae1c9f0	65	// Constructor
e40b9538	66
	67	}
	68
	69	//_________________________________________________________________________
	70	AliGenEMCocktail::~AliGenEMCocktail()
	71	{
4ae1c9f0	72	// Destructor
e40b9538	73
	74	}
	75
	76	//_________________________________________________________________________
	77	void AliGenEMCocktail::CreateCocktail()
	78	{
4ae1c9f0	79	// create and add sources to the cocktail
e40b9538	80
	81	fDecayer->SetForceDecay(fDecayMode);
	82	fDecayer->ForceDecay();
	83
4ae1c9f0	84	// Set kinematic limits
e40b9538	85	Double_t ptMin = fPtMin;
	86	Double_t ptMax = fPtMax;
	87	Double_t yMin = fYMin;;
	88	Double_t yMax = fYMax;;
	89	Double_t phiMin = fPhiMin*180./TMath::Pi();
	90	Double_t phiMax = fPhiMax*180./TMath::Pi();
	91	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));
	92	AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
	93
4ae1c9f0	94	//Initialize user selection for Pt Parameterization and centrality:
	95	AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
	96
	97	// Create and add electron sources to the generator
e40b9538	98
4ae1c9f0	99	// pizero
4ae1c9f0	100	AliGenParam *genpizero=0;
e40b9538	101	Char_t namePizero[10];
4ae1c9f0	102	snprintf(namePizero,10,"Pizero");
e40b9538	103	genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
	104	AddSource2Generator(namePizero,genpizero);
	105	TF1 *fPtPizero = genpizero->GetPt();
faf00b0d	106	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	107	fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
	108	#else
4ae1c9f0	109	fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	110	#endif
4ae1c9f0	111
	112	// eta
	113	if(fHeaviestParticle<kGenEta)return;
	114	AliGenParam *geneta=0;
e40b9538	115	Char_t nameEta[10];
4ae1c9f0	116	snprintf(nameEta,10,"Eta");
e40b9538	117	geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
	118	AddSource2Generator(nameEta,geneta);
	119	TF1 *fPtEta = geneta->GetPt();
4ae1c9f0	120	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
faf00b0d	121	fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
4ae1c9f0	122	#else
4ae1c9f0	123	fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	124	#endif
4ae1c9f0	125
	126	// rho
	127	if(fHeaviestParticle<kGenRho)return;
	128	AliGenParam *genrho=0;
e40b9538	129	Char_t nameRho[10];
4ae1c9f0	130	snprintf(nameRho,10,"Rho");
e40b9538	131	genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
	132	AddSource2Generator(nameRho,genrho);
	133	TF1 *fPtRho = genrho->GetPt();
faf00b0d	134	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	135	fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
	136	#else
4ae1c9f0	137	fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	138	#endif
4ae1c9f0	139
	140	// omega
	141	if(fHeaviestParticle<kGenOmega)return;
	142	AliGenParam *genomega=0;
e40b9538	143	Char_t nameOmega[10];
4ae1c9f0	144	snprintf(nameOmega,10,"Omega");
e40b9538	145	genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
	146	AddSource2Generator(nameOmega,genomega);
	147	TF1 *fPtOmega = genomega->GetPt();
faf00b0d	148	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	149	fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
	150	#else
4ae1c9f0	151	fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	152	#endif
4ae1c9f0	153
	154	// etaprime
	155	if(fHeaviestParticle<kGenEtaprime)return;
	156	AliGenParam *genetaprime=0;
e40b9538	157	Char_t nameEtaprime[10];
4ae1c9f0	158	snprintf(nameEtaprime,10,"Etaprime");
e40b9538	159	genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
	160	AddSource2Generator(nameEtaprime,genetaprime);
	161	TF1 *fPtEtaprime = genetaprime->GetPt();
faf00b0d	162	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	163	fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
	164	#else
4ae1c9f0	165	fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	166	#endif
4ae1c9f0	167
	168	// phi
	169	if(fHeaviestParticle<kGenPhi)return;
	170	AliGenParam *genphi=0;
e40b9538	171	Char_t namePhi[10];
4ae1c9f0	172	snprintf(namePhi,10,"Phi");
e40b9538	173	genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
	174	AddSource2Generator(namePhi,genphi);
	175	TF1 *fPtPhi = genphi->GetPt();
faf00b0d	176	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	177	fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
	178	#else
4ae1c9f0	179	fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d	180	#endif
4ae1c9f0	181
	182	// jpsi
	183	if(fHeaviestParticle<kGenJpsi)return;
	184	AliGenParam *genjpsi=0;
	185	Char_t nameJpsi[10];
	186	snprintf(nameJpsi,10,"Jpsi");
	187	genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
	188	AddSource2Generator(nameJpsi,genjpsi);
	189	TF1 *fPtJpsi = genjpsi->GetPt();
	190	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
	191	fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
	192	#else
	193	fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
	194	#endif
	195
e40b9538	196	}
	197
	198	//-------------------------------------------------------------------
	199	void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
	200	AliGenParam* const genSource)
	201	{
	202	// add sources to the cocktail
	203	Double_t phiMin = fPhiMin*180./TMath::Pi();
	204	Double_t phiMax = fPhiMax*180./TMath::Pi();
	205
	206	genSource->SetPtRange(fPtMin, fPtMax);
	207	genSource->SetYRange(fYMin, fYMax);
	208	genSource->SetPhiRange(phiMin, phiMax);
	209	genSource->SetWeighting(fWeightingMode);
4ae1c9f0	210	genSource->SetForceGammaConversion(fForceConv);
e40b9538	211	if (!gMC) genSource->SetDecayer(fDecayer);
	212	genSource->Init();
	213
	214	AddGenerator(genSource,nameSource,1.); // Adding Generator
	215	}
	216
	217	//-------------------------------------------------------------------
	218	void AliGenEMCocktail::Init()
	219	{
4ae1c9f0	220	// Initialisation
e40b9538	221	TIter next(fEntries);
	222	AliGenCocktailEntry *entry;
	223	if (fStack) {
	224	while((entry = (AliGenCocktailEntry*)next())) {
	225	entry->Generator()->SetStack(fStack);
	226	}
	227	}
	228	}
	229
	230	//_________________________________________________________________________
	231	void AliGenEMCocktail::Generate()
	232	{
4ae1c9f0	233	// Generate event
e40b9538	234	TIter next(fEntries);
	235	AliGenCocktailEntry *entry = 0;
	236	AliGenCocktailEntry *preventry = 0;
	237	AliGenerator* gen = 0;
	238
	239	if (fHeader) delete fHeader;
	240	fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
	241
	242	const TObjArray *partArray = gAlice->GetMCApp()->Particles();
	243
4ae1c9f0	244	// Generate the vertex position used by all generators
e40b9538	245	if(fVertexSmear == kPerEvent) Vertex();
e40b9538	246
4ae1c9f0	247	//Reseting stack
e40b9538	248	AliRunLoader * runloader = AliRunLoader::Instance();
	249	if (runloader)
	250	if (runloader->Stack())
	251	runloader->Stack()->Clean();
	252
4ae1c9f0	253	// Loop over generators and generate events
e40b9538	254	Int_t igen = 0;
6078e216	255	Float_t evPlane;
	256	Rndm(&evPlane,1);
	257	evPlane=TMath::Pi()2;
e40b9538	258	while((entry = (AliGenCocktailEntry*)next())) {
e40b9538	259	gen = entry->Generator();
e40b9538	260	gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
	261
	262	if (fNPart > 0) {
	263	igen++;
	264	if (igen == 1) entry->SetFirst(0);
	265	else entry->SetFirst((partArray->GetEntriesFast())+1);
6078e216	266	gen->SetEventPlane(evPlane);
e40b9538	267	gen->SetNumberParticles(fNPart);
	268	gen->Generate();
	269	entry->SetLast(partArray->GetEntriesFast());
	270	preventry = entry;
	271	}
	272	}
	273	next.Reset();
	274
4ae1c9f0	275	// Setting weights for proper absolute normalization
e40b9538	276	Int_t iPart, iMother;
	277	Int_t pdgMother = 0;
	278	Double_t weight = 0.;
	279	Double_t dNdy = 0.;
	280	Int_t maxPart = partArray->GetEntriesFast();
	281	for(iPart=0; iPart<maxPart; iPart++){
	282	TParticle *part = gAlice->GetMCApp()->Particle(iPart);
	283	iMother = part->GetFirstMother();
	284	TParticle *mother = 0;
	285	if (iMother>=0){
	286	mother = gAlice->GetMCApp()->Particle(iMother);
	287	pdgMother = mother->GetPdgCode();
	288	}
	289	else
	290	pdgMother = part->GetPdgCode();
	291	switch (pdgMother){
	292	case 111:
	293	dNdy = fYieldArray[kGenPizero];
	294	break;
	295	case 221:
	296	dNdy = fYieldArray[kGenEta];
	297	break;
	298	case 113:
	299	dNdy = fYieldArray[kGenRho];
	300	break;
	301	case 223:
	302	dNdy = fYieldArray[kGenOmega];
	303	break;
	304	case 331:
	305	dNdy = fYieldArray[kGenEtaprime];
	306	break;
	307	case 333:
	308	dNdy = fYieldArray[kGenPhi];
	309	break;
4ae1c9f0	310	case 443:
	311	dNdy = fYieldArray[kGenJpsi];
	312	break;
e40b9538	313
	314	default:
	315	dNdy = 0.;
	316	}
e40b9538	317	weight = dNdy*part->GetWeight();
	318	part->SetWeight(weight);
	319	}
	320
	321	fHeader->SetNProduced(maxPart);
	322
	323
	324	TArrayF eventVertex;
	325	eventVertex.Set(3);
	326	for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
	327
	328	fHeader->SetPrimaryVertex(eventVertex);
	329
	330	gAlice->SetGenEventHeader(fHeader);
	331	}