[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 <TDatabasePDG.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(AliGenEMlib::kPizero7TeVpp),
  fCentrality(AliGenEMlib::kpp),
  fV2Systematic(AliGenEMlib::kNoV2Sys),
  fForceConv(kFALSE),
  fSelectedParticles(kGenPizero)
{
  // Constructor
  SetHeaviestHadron(kGenPhi);
}

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

}

//_________________________________________________________________________
void AliGenEMCocktail::SetHeaviestHadron(ParticeGenerator_t part)
{
  Int_t val=kGenPizero;
  while(val<part) val|=val<<1;

  fSelectedParticles=val;
}

//_________________________________________________________________________
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
  if(fSelectedParticles&kGenPizero){
  AliGenParam *genpizero=0;
  Char_t namePizero[10];    
  snprintf(namePizero,10,"Pizero");    
    //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
  genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
    //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
  genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
  AddSource2Generator(namePizero,genpizero);
  TF1 *fPtPizero = genpizero->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
    fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
#else
    fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif
  }

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

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

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

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

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

  // direct gamma
  if(fDecayMode!=kGammaEM) return;
  if(fSelectedParticles&kGenDirectRealGamma){
    AliGenParam *genDirectRealG=0;
    Char_t nameDirectRealG[10];    
    snprintf(nameDirectRealG,10,"DirectRealGamma");    
    genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
    genDirectRealG->SetYRange(fYMin, fYMax);
    AddSource2Generator(nameDirectRealG,genDirectRealG);
    TF1 *fPtDirectRealG = genDirectRealG->GetPt();
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
    fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
#else
    fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
#endif
  }

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