* EVGEN/AliGenEMCocktail.cxx (& .h) is the class that generates the cocktail follow...
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 25 Sep 2010 08:58:05 +0000 (08:58 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 25 Sep 2010 08:58:05 +0000 (08:58 +0000)
 *  EVGEN/AliGenEMlib.cxx (& .h) is the class that contains the kinematic distributions of the sources, also following the muon example. At the moment, it contains
    only one pT distribution plus mT scaling. In the near future I'll extend this class with further parameterizations. Also, I'll add the functionality to set
    parameters by proper set methods (instead of the hard-wired parameters as used currently in all of the similar lib classes).
 *  EVGEN/libEVGEN.pkg and EVGEN/EVGENLinkDef.h have been modified to include the two new classes.
 *  EVGEN/AliDecayer.h has been modified by adding three new decay modes (kElectronEM, kDiElectronEM, and kGammaEM) between which one can select in the cocktail
    class.
 *  EVGEN/AliGenMC.cxx has been modified such that it "knows about" the three new decay modes in EVGEN/AliDecayer.h
 *  PYTHIA6/AliDecayerPythia.cxx has been modified accordingly. Here, the corresponding particle decays are forced depending on the selected decay mode.
 *  FASTSIM/fastGenEMCocktail.C is an example macro to generate a single electron cocktail.
(Ralf Averbeck)

EVGEN/AliDecayer.h
EVGEN/AliGenEMCocktail.cxx [new file with mode: 0644]
EVGEN/AliGenEMCocktail.h [new file with mode: 0644]
EVGEN/AliGenEMlib.cxx [new file with mode: 0644]
EVGEN/AliGenEMlib.h [new file with mode: 0644]
EVGEN/AliGenMC.cxx
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg
FASTSIM/fastGenEMCocktail.C [new file with mode: 0644]
PYTHIA6/AliDecayerPythia.cxx

index 3fe4f37..8f7d752 100644 (file)
@@ -19,7 +19,8 @@ typedef enum
     kNoDecay, kHadronicD, kHadronicDWithout4Bodies, kOmega, kLambda, kPhiKK, 
     kAll, kNoDecayHeavy, kHardMuons, kBJpsi,
     kWToMuon,kWToCharm, kWToCharmToMuon, kZDiMuon, kZDiElectron, kNeutralPion, kAllMuonic,
-    kChiToJpsiGammaToMuonMuon, kChiToJpsiGammaToElectronElectron, kNoDecayBeauty, kPsiPrimeJpsiDiElectron
+    kChiToJpsiGammaToMuonMuon, kChiToJpsiGammaToElectronElectron, kNoDecayBeauty, kPsiPrimeJpsiDiElectron,
+    kElectronEM, kGammaEM, kDiElectronEM
 } Decay_t;
 #endif
 
diff --git a/EVGEN/AliGenEMCocktail.cxx b/EVGEN/AliGenEMCocktail.cxx
new file mode 100644 (file)
index 0000000..c4b256a
--- /dev/null
@@ -0,0 +1,270 @@
+/**************************************************************************
+ * 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);
+}
+
+    
diff --git a/EVGEN/AliGenEMCocktail.h b/EVGEN/AliGenEMCocktail.h
new file mode 100644 (file)
index 0000000..db823b9
--- /dev/null
@@ -0,0 +1,60 @@
+#ifndef ALIGENEMCOCKTAIL_H
+#define ALIGENEMCOCKTAIL_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//
+// 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
+//
+
+#include "AliGenCocktail.h"
+#include "AliGenParam.h"
+#include "AliDecayer.h"
+#include "AliPythia.h"
+
+class AliGenCocktailEntry;
+
+class AliGenEMCocktail : public AliGenCocktail
+{
+ public:
+
+    AliGenEMCocktail();
+    enum GeneratorCode { kGenPizero, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGENs };    
+
+    virtual ~AliGenEMCocktail();    
+    virtual void Init();
+    virtual void CreateCocktail();
+    virtual void Generate();    
+    Float_t GetDecayMode()         const {return fDecayMode;}
+    Float_t GetWeightingMode()     const {return fWeightingMode;}
+    
+    void    SetDecayer(AliDecayer* const decayer){fDecayer = decayer;}
+    void    SetDecayMode(Decay_t decay){ fDecayMode = decay;}
+    void    SetWeightingMode(Weighting_t weight){ fWeightingMode = weight;}
+    void    SetNPart(Int_t npart){ fNPart = npart;}
+    
+ protected:
+
+    //
+ private:
+    AliGenEMCocktail(const AliGenEMCocktail &cocktail); 
+    AliGenEMCocktail & operator=(const AliGenEMCocktail &cocktail); 
+
+    void AddSource2Generator(Char_t *nameReso, AliGenParam* const genReso);
+    
+    AliDecayer* fDecayer;        // External decayer
+    Decay_t fDecayMode;          // decay mode in which resonances are forced to decay, default: kAll
+    Weighting_t fWeightingMode;  // weighting mode: kAnalog or kNonAnalog
+    
+    Int_t    fNPart;             // multiplicity of each source per event
+    Double_t fYieldArray[kGENs]; // array of dN/dy for each source
+
+    ClassDef(AliGenEMCocktail,1)  //  cocktail for EM physics
+};
+
+#endif
+
+
+
diff --git a/EVGEN/AliGenEMlib.cxx b/EVGEN/AliGenEMlib.cxx
new file mode 100644 (file)
index 0000000..b3ba4fa
--- /dev/null
@@ -0,0 +1,359 @@
+/**************************************************************************
+ * 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: AliGenEMlib.cxx 30052 2008-11-25 14:54:18Z morsch $ */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// Implementation of AliGenEMlib for electron, di-electron, and photon     //
+// cocktail calculations.                                                  //
+// It is based on AliGenGSIlib.                                            //
+//                                                                         //
+// Responsible: R.Averbeck@gsi.de                                          //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+
+#include "TMath.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "AliGenEMlib.h"
+
+
+ClassImp(AliGenEMlib)
+
+//==========================================================================
+//
+//              Definition of Particle Distributions
+//                    
+//==========================================================================
+
+//--------------------------------------------------------------------------
+//
+//                              Pizero
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpPizero(TRandom *)
+{
+// Return pizero pdg code
+  return 111;     
+}
+
+Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Generate pizero pT distribution from modified Hagedorn parameterization
+// taken from fit to unidentified hadrons in pp at 7 TeV
+  const Double_t kc=0.000565;
+  const Double_t kp0=0.2472;
+  const Double_t kp1=4.354;
+  const Double_t kn=7.007;
+  Double_t invYield;
+  Double_t x=*px;
+
+  invYield = kc/TMath::Power(kp0+x/kp1,kn);
+
+  return invYield*(2*TMath::Pi()*x);
+   
+}
+
+Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Eta
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpEta(TRandom *)
+{
+// Return eta pdg code
+  return 221;     
+}
+
+Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Eta pT
+  return MtScal(*px,2);
+}
+
+Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Rho
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpRho(TRandom *)
+{
+// Return rho pdg code
+  return 113;     
+}
+
+Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Rho pT
+  return MtScal(*px,3);
+}
+
+Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Omega
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpOmega(TRandom *)
+{
+// Return omega pdg code
+  return 223;     
+}
+
+Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Omega pT
+  return MtScal(*px,4);
+}
+
+Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Etaprime
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpEtaprime(TRandom *)
+{
+// Return etaprime pdg code
+  return 331;     
+}
+
+Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Eta pT
+  return MtScal(*px,5);
+}
+
+Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Phi
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpPhi(TRandom *)
+{
+// Return phi pdg code
+  return 333;     
+}
+
+Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
+{
+// Phi pT
+  return MtScal(*px,6);
+}
+
+Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::YFlat(Double_t y)
+{
+//--------------------------------------------------------------------------
+//
+//                    flat rapidity distribution 
+//
+//--------------------------------------------------------------------------
+
+  Double_t dNdy = 1.;   
+
+  return dNdy;
+
+}
+
+//============================================================= 
+//
+//                    Mt-scaling  
+//
+//============================================================= 
+//
+ Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
+{
+// Function for the calculation of the Pt distribution for a 
+// given particle np, from the pizero Pt distribution using  
+// mt scaling. 
+
+// MASS   1=>PIZERO, 2=>ETA, 3=>RHO, 4=>OMEGA, 5=>ETAPRIME, 6=>PHI
+
+  const Double_t khm[6] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946};
+  np--;
+
+  Double_t scaledPt = sqrt(pt*pt + khm[np]*khm[np] - khm[0]*khm[0]);
+  Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
+
+  //     VALUE MESON/PI AT 5 GEV
+
+  Double_t normPt = 5.;
+  Double_t scaledNormPt = sqrt(normPt*normPt + khm[np]*khm[np] - khm[0]*khm[0]);
+  const Double_t kfmax[6]={1., 0.48, 1.0, 0.9, 0.25, 0.4};
+
+  Double_t norm = kfmax[np] * (PtPizero(&normPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
+
+  return norm*(pt/scaledPt)*scaledYield;
+}
+
+//==========================================================================
+//
+//                     Set Getters 
+//    
+//==========================================================================
+
+typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
+
+typedef Int_t (*GenFuncIp) (TRandom *);
+
+GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
+{
+// Return pointer to pT parameterisation
+   GenFunc func=0;
+   TString sname(tname);
+
+   switch (param) 
+    {
+    case kPizero:
+      func=PtPizero;
+      break;
+    case kEta:
+      func=PtEta;
+      break;
+    case kRho:
+      func=PtRho;
+      break;
+    case kOmega:
+      func=PtOmega;
+      break;
+    case kEtaprime:
+      func=PtEtaprime;
+      break;
+    case kPhi:
+      func=PtPhi;
+      break;
+
+    default:
+      func=0;
+      printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
+    }
+   return func;
+}
+
+GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
+{
+// Return pointer to y- parameterisation
+   GenFunc func=0;
+   TString sname(tname);
+
+   switch (param) 
+    {
+    case kPizero:
+         func=YPizero;
+         break;
+    case kEta:
+         func=YEta;
+         break;
+    case kRho:
+         func=YRho;
+         break;
+    case kOmega:
+         func=YOmega;
+         break;
+    case kEtaprime:
+         func=YEtaprime;
+         break;
+    case kPhi:
+         func=YPhi;
+         break;
+
+    default:
+        func=0;
+        printf("<AliGenEMlib::GetY> unknown parametrisation\n");
+    }
+    return func;
+}
+
+GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
+{
+// Return pointer to particle type parameterisation
+   GenFuncIp func=0;
+   TString sname(tname);
+
+   switch (param) 
+    {
+    case kPizero:
+         func=IpPizero;
+         break;
+    case kEta:
+         func=IpEta;
+         break;
+    case kRho:
+         func=IpRho;
+         break;
+    case kOmega:
+         func=IpOmega;
+         break;
+    case kEtaprime:
+         func=IpEtaprime;
+         break;
+    case kPhi:
+         func=IpPhi;
+         break;
+
+    default:
+        func=0;
+        printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
+    }
+    return func;
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/EVGEN/AliGenEMlib.h b/EVGEN/AliGenEMlib.h
new file mode 100644 (file)
index 0000000..76ae20b
--- /dev/null
@@ -0,0 +1,70 @@
+#ifndef ALIGENEMLIB_H
+#define ALIGENEMLIB_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliGenEMlib.h 30052 2008-11-25 14:54:18Z morsch $ */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// Implementation of AliGenEMlib for electron, di-electron, and photon     //
+// cocktail calculations.                                                  //
+// It is based on AliGenGSIlib.                                            //
+//                                                                         //
+// Responsible: R.Averbeck@gsi.de                                          //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#include "AliGenLib.h"
+class TRandom;
+
+class AliGenEMlib :public AliGenLib {
+ public:
+    GenFunc   GetPt(Int_t param, const char * tname=0) const;
+    GenFunc   GetY(Int_t param, const char * tname=0) const;
+    GenFuncIp GetIp(Int_t param, const char * tname=0) const;    
+
+    enum constants{kPizero, kEta, kRho, kOmega, kEtaprime, kPhi};
+
+ private:
+
+// Pizero
+    static Int_t    IpPizero(TRandom *ran);
+    static Double_t PtPizero( const Double_t *px, const Double_t *dummy );
+    static Double_t YPizero(const Double_t *py, const Double_t *dummy);
+    static Double_t MtScal(Double_t pt, Int_t np);
+
+// Eta
+    static Int_t    IpEta(TRandom *ran);
+    static Double_t PtEta( const Double_t *px, const Double_t *dummy );
+    static Double_t YEta(const Double_t *py, const Double_t *dummy);
+
+// Rho
+    static Int_t    IpRho(TRandom *ran);
+    static Double_t PtRho( const Double_t *px, const Double_t *dummy );
+    static Double_t YRho(const Double_t *py, const Double_t *dummy);
+
+// Omega
+    static Int_t    IpOmega(TRandom *ran);
+    static Double_t PtOmega( const Double_t *px, const Double_t *dummy );
+    static Double_t YOmega(const Double_t *py, const Double_t *dummy);
+
+// Etaprime
+    static Int_t    IpEtaprime(TRandom *ran);
+    static Double_t PtEtaprime( const Double_t *px, const Double_t *dummy );
+    static Double_t YEtaprime(const Double_t *py, const Double_t *dummy);
+
+// Phi
+    static Int_t    IpPhi(TRandom *ran);
+    static Double_t PtPhi( const Double_t *px, const Double_t *dummy );
+    static Double_t YPhi(const Double_t *py, const Double_t *dummy);
+
+// General
+    static Double_t PtFlat(const Double_t *px, const Double_t *dummy);
+    static Double_t YFlat(Double_t y);
+
+
+  ClassDef(AliGenEMlib,0)
+};
+
+#endif
index 6660454..1ce5e58 100644 (file)
@@ -121,6 +121,8 @@ void AliGenMC::Init()
     case kDiElectron:
     case kBJpsiDiElectron:
     case kBPsiPrimeDiElectron:
+    case kElectronEM:
+    case kDiElectronEM:
        fChildSelect[0] = kElectron;    
        break;
     case kHardMuons:   
@@ -163,9 +165,12 @@ void AliGenMC::Init()
         fChildSelect[1]= 211;
        break;
     case kPsiPrimeJpsiDiElectron:
-      fChildSelect[0]= 211;
-      fChildSelect[1]= 11;
-      break;
+        fChildSelect[0]= 211;
+        fChildSelect[1]= 11;
+       break;
+    case kGammaEM:
+        fChildSelect[0] = kGamma;      
+        break;
     case kOmega:
     case kAll:
     case kAllMuonic:
index 9bf193d..e734358 100644 (file)
@@ -21,6 +21,7 @@
 #pragma link C++ class  AliGenCocktail+;
 #pragma link C++ class  AliGenMUONCocktail+;
 #pragma link C++ class  AliGenMUONCocktailpp+;
+#pragma link C++ class  AliGenEMCocktail+;
 #pragma link C++ class  AliGenCocktailAfterBurner+;
 #pragma link C++ class  AliGenCocktailEntry+;
 #pragma link C++ class  AliGenExtFile+;
@@ -35,6 +36,7 @@
 #pragma link C++ class  AliDimuCombinator+;
 #pragma link C++ class  AliGenPHOSlib+;
 #pragma link C++ class  AliGenGSIlib+;
+#pragma link C++ class  AliGenEMlib+;
 #pragma link C++ class  AliGenPMDlib+;
 #pragma link C++ class  AliGenSTRANGElib+;
 #pragma link C++ class  AliGenMC+;
index 4ce1be3..4c760bc 100644 (file)
@@ -8,8 +8,8 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                 AliGenHaloProtvino.cxx \
                 AliGenExtFile.cxx AliGenScan.cxx AliGenPHOSlib.cxx \
                 AliGenDoubleScan.cxx AliGenCocktailEntry.cxx \
-               AliGenGSIlib.cxx AliGenPMDlib.cxx\
-               AliGenMC.cxx\
+               AliGenGSIlib.cxx AliGenEMlib.cxx \
+               AliGenPMDlib.cxx AliGenMC.cxx\
                AliGenReader.cxx AliGenReaderCwn.cxx AliGenReaderTreeK.cxx \
                 AliGenReaderEcalHijing.cxx AliGenReaderEcalJets.cxx\
                AliGenHIJINGparaBa.cxx AliGeVSimParticle.cxx AliGenGeVSim.cxx\
@@ -17,7 +17,8 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliGenAfterBurnerFlow.cxx \
                AliGenSlowNucleons.cxx \
                AliSlowNucleonModel.cxx AliSlowNucleonModelExp.cxx \
-               AliGenMUONCocktail.cxx AliGenMUONCocktailpp.cxx AliGenHBTosl.cxx \
+               AliGenMUONCocktail.cxx AliGenMUONCocktailpp.cxx \
+               AliGenEMCocktail.cxx AliGenHBTosl.cxx \
                AliGenReaderEMD.cxx AliDecayerPolarized.cxx AliGenCorrHF.cxx AliGenCosmicsParam.cxx \
                AliGenKrypton.cxx AliGenThermalPhotons.cxx AliGenPromptPhotons.cxx\
                AliGenPileup.cxx AliGenFunction.cxx AliGenTHnSparse.cxx\
diff --git a/FASTSIM/fastGenEMCocktail.C b/FASTSIM/fastGenEMCocktail.C
new file mode 100644 (file)
index 0000000..b6e822a
--- /dev/null
@@ -0,0 +1,116 @@
+void fastGenEMCocktail(Int_t nev = 1, char* filename = "galice.root")
+{
+// load libraries
+
+  gSystem->Load("liblhapdf.so");      // Parton density functions
+  gSystem->Load("libpythia6.so");     // Pythia
+  gSystem->Load("libEG");
+  gSystem->Load("libEGPythia6");
+  gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
+
+// Runloader
+    
+  AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
+    
+  rl->SetCompressionLevel(2);
+  rl->SetNumberOfEventsPerFile(nev);
+  rl->LoadKinematics("RECREATE");
+  rl->MakeTree("E");
+  gAlice->SetRunLoader(rl);
+
+// Create stack
+  rl->MakeStack();
+  AliStack* stack      = rl->Stack();
+// Header
+  AliHeader* header = rl->GetHeader();
+
+// Create and Initialize Generator
+
+  AliGenEMCocktail *gener = new AliGenEMCocktail();
+
+//=======================================================================
+// Set External decayer
+  TVirtualMCDecayer *decayer = new AliDecayerPythia();
+
+  gener->SetNPart(1000);               // source multiplicity per event
+  gener->SetPtRange(0.,20.);
+  gener->SetYRange(-1.,1.);
+  gener->SetPhiRange(0., 360.);
+  gener->SetOrigin(0.,0.,0.); 
+  gener->SetSigma(0.,0.,0.);
+  gener->SetVertexSmear(kPerEvent);
+  gener->SetTrackingFlag(0);
+  gener->SetDecayMode(kElectronEM);    // select cocktail:
+                                       // kElectronEM   => single electron
+                                       // kDiElectronEM => electron-positron
+                                       // kGammaEM      => single photon
+  gener->SetDecayer(decayer);
+  gener->SetWeightingMode(kNonAnalog); // select weighting:
+                                       // kNonAnalog => weight ~ dN/dp_T
+                                       // kAnalog    => weight ~ 1
+  gener->CreateCocktail();
+  gener->Init();
+  gener->SetStack(stack);
+
+//
+//                        Event Loop
+//
+  Int_t iev;
+     
+  for (iev = 0; iev < nev; iev++) {
+
+    printf("\n \n Event number %d \n \n", iev);
+       
+// Initialize event
+    header->Reset(0,iev);
+    rl->SetEventNumber(iev);
+    stack->Reset();
+    rl->MakeTree("K");
+//  stack->ConnectTree();
+    
+// Generate event
+    gener->Generate();
+// Analysis
+    Int_t npart = stack->GetNprimary();
+    printf("Analyse %d Particles\n", npart);
+    for (Int_t part=0; part<npart; part++) {
+      TParticle *MPart = stack->Particle(part);
+      Int_t mpart  = MPart->GetPdgCode();
+//      printf("Particle %d\n", mpart);
+    }
+       
+// Finish event
+    header->SetNprimary(stack->GetNprimary());
+    header->SetNtrack(stack->GetNtrack());  
+
+// I/O
+//     
+    stack->FinishEvent();
+    header->SetStack(stack);
+    rl->TreeE()->Fill();
+    rl->WriteKinematics("OVERWRITE");
+
+  } // event loop
+
+//
+//                         Termination
+// Generator
+  gener->FinishRun();
+// Write file
+  rl->WriteHeader("OVERWRITE");
+  gener->Write();
+  rl->Write();
+}
+
+
+
+
+
+
+
+
+
+
+
+
index 19e895b..5cadc49 100644 (file)
@@ -391,6 +391,30 @@ void AliDecayerPythia::ForceDecay()
     case kNoDecayBeauty:                                                    
         SwitchOffBDecay(); 
        break;
+    case kElectronEM:
+        ForceParticleDecay(  111,11,1); // pi^0     
+       ForceParticleDecay(  221,11,1); // eta     
+       ForceParticleDecay(  113,11,1); // rho     
+       ForceParticleDecay(  223,11,1); // omega     
+       ForceParticleDecay(  331,11,1); // etaprime     
+       ForceParticleDecay(  333,11,1); // phi     
+       break;
+    case kDiElectronEM:
+       ForceParticleDecay(  111,11,2); // pi^0     
+       ForceParticleDecay(  221,11,2); // eta     
+       ForceParticleDecay(  113,11,2); // rho     
+       ForceParticleDecay(  223,11,2); // omega     
+       ForceParticleDecay(  331,11,2); // etaprime     
+       ForceParticleDecay(  333,11,2); // phi     
+       break;
+    case kGammaEM:
+       ForceParticleDecay(  111,22,1); // pi^0     
+       ForceParticleDecay(  221,22,1); // eta     
+       ForceParticleDecay(  113,22,1); // rho     
+       ForceParticleDecay(  223,22,1); // omega     
+       ForceParticleDecay(  331,22,1); // etaprime     
+       ForceParticleDecay(  333,22,1); // phi     
+       break;
     }
 }