EVGEN/AliGenMUONCocktailpp.cxx

   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$ */
  17
  18 //
  19 // Classe to create the coktail for physics with muons for pp at 14 TeV
  20 // The followoing sources:
  21 // jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
  22 // The free parameeters are :
  23 //      pp reaction cross-section
  24 //      production cross-sections in pp collisions
  25 //
  26
  27 #include <TObjArray.h>
  28 #include <TParticle.h>
  29 #include <TF1.h>
  30
  31 #include "AliGenCocktailEventHeader.h"
  32
  33 #include "AliGenCocktailEntry.h"
  34 #include "AliGenMUONCocktailpp.h"
  35 #include "AliGenMUONlib.h"
  36 #include "AliGenParam.h"
  37 #include "AliMC.h"
  38 #include "AliRun.h"
  39 #include "AliStack.h"
  40 #include "AliLog.h"
  41 #include "AliGenPythia.h"
  42
  43
  44 ClassImp(AliGenMUONCocktailpp)
  45
  46 //________________________________________________________________________
  47 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
  48   :AliGenCocktail()
  49 {
  50 // Constructor
  51   fTotalRate =0;
  52   fNSucceded=0;
  53   fNGenerated=0;
  54   fMuonMultiplicity=0;
  55   fMuonPtCut= 0.;
  56   fMuonThetaMinCut=0.;
  57   fMuonThetaMaxCut=0.;
  58 }
  59 //_________________________________________________________________________
  60 AliGenMUONCocktailpp::AliGenMUONCocktailpp(const AliGenMUONCocktailpp & cocktail):
  61     AliGenCocktail(cocktail)
  62 {
  63 // Copy constructor
  64   fTotalRate =0;
  65   fNSucceded=0;
  66   fNGenerated=0;
  67   fMuonMultiplicity=0;
  68   fMuonPtCut=0.;
  69   fMuonThetaMinCut=0.;
  70   fMuonThetaMaxCut=0.;
  71 }
  72 //_________________________________________________________________________
  73 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
  74 // Destructor
  75 {
  76
  77 }
  78
  79 //_________________________________________________________________________
  80 void AliGenMUONCocktailpp::Init()
  81 {
  82 // MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
  83     Double_t sigmaReaction =   0.070;
  84
  85 // These limits are only used for renormalization of quarkonia cross section
  86 // Pythia events are generated in 4pi
  87     Double_t ptMin  = fPtMin;
  88     Double_t ptMax  = fPtMax;
  89     Double_t yMin   = fYMin;;
  90     Double_t yMax   = fYMax;;
  91     Double_t phiMin = fPhiMin*180./TMath::Pi();
  92     Double_t phiMax = fPhiMax*180./TMath::Pi();
  93     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));
  94
  95 // Ratios with respect to the reaction cross-section in the
  96 // kinematics limit of the MUONCocktail
  97     Double_t ratiojpsi;
  98     Double_t ratiopsiP;
  99     Double_t ratioupsilon;
 100     Double_t ratioupsilonP;
 101     Double_t ratioupsilonPP;
 102
 103 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV
 104
 105     Double_t sigmajpsi = 53.85e-6;
 106     Double_t sigmapsiP = 7.67e-6;
 107     Double_t sigmaupsilon = 1.15e-6;
 108     Double_t sigmaupsilonP = 0.526e-6;
 109     Double_t sigmaupsilonPP = 0.228e-6;
 110
 111 // Generation using CDF scaled distribution
 112 //(See http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 113
 114 //------------------------------------------------------------------
 115 // Jpsi generator
 116     AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
 117 // first step: generation in 4pi
 118     genjpsi->SetPtRange(0.,100.);
 119     genjpsi->SetYRange(-8.,8.);
 120     genjpsi->SetPhiRange(0.,360.);
 121     genjpsi->SetForceDecay(kAll);
 122     genjpsi->Init();  // generation in 4pi
 123     ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
 124     AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
 125     AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
 126 // second step: generation in selected kinematical range
 127     genjpsi->SetPtRange(ptMin, ptMax);
 128     genjpsi->SetYRange(yMin, yMax);
 129     genjpsi->SetPhiRange(phiMin, phiMax);
 130     genjpsi->Init(); // generation in selected kinematic range
 131     AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
 132     fTotalRate+=ratiojpsi;
 133
 134 //------------------------------------------------------------------
 135 // Psi prime generator
 136     AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP");
 137 // first step: generation in 4pi
 138     genpsiP->SetPtRange(0.,100.);
 139     genpsiP->SetYRange(-8.,8.);
 140     genpsiP->SetPhiRange(0.,360.);
 141     genpsiP->SetForceDecay(kAll);
 142     genpsiP->Init();  // generation in 4pi
 143     ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 144     AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
 145     AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
 146 // second step: generation in selected kinematical range
 147     genpsiP->SetPtRange(ptMin, ptMax);
 148     genpsiP->SetYRange(yMin, yMax);
 149     genpsiP->SetPhiRange(phiMin, phiMax);
 150     genpsiP->Init(); // generation in selected kinematic range
 151     AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
 152     fTotalRate+=ratiopsiP;
 153
 154 //------------------------------------------------------------------
 155 // Upsilon 1S generator
 156     AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon");
 157 // first step: generation in 4pi
 158     genupsilon->SetPtRange(0.,100.);
 159     genupsilon->SetYRange(-8.,8.);
 160     genupsilon->SetPhiRange(0.,360.);
 161     genupsilon->SetForceDecay(kAll);
 162     genupsilon->Init();  // generation in 4pi
 163     ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 164     AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
 165     AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
 166 // second step: generation in selected kinematical range
 167     genupsilon->SetPtRange(ptMin, ptMax);
 168     genupsilon->SetYRange(yMin, yMax);
 169     genupsilon->SetPhiRange(phiMin, phiMax);
 170     genupsilon->Init(); // generation in selected kinematical range
 171     AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
 172     fTotalRate+=ratioupsilon;
 173
 174 //------------------------------------------------------------------
 175 // Upsilon 2S generator
 176     AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF scaled", "UpsilonP");
 177 // first step: generation in 4pi
 178     genupsilonP->SetPtRange(0.,100.);
 179     genupsilonP->SetYRange(-8.,8.);
 180     genupsilonP->SetPhiRange(0.,360.);
 181     genupsilonP->SetForceDecay(kAll);
 182     genupsilonP->Init();  // generation in 4pi
 183     ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 184     AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
 185     AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
 186 // second step: generation in the kinematical range
 187     genupsilonP->SetPtRange(ptMin, ptMax);
 188     genupsilonP->SetYRange(yMin, yMax);
 189     genupsilonP->SetPhiRange(phiMin, phiMax);
 190     genupsilonP->Init(); // generation in selected kinematical range
 191     AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
 192     fTotalRate+=ratioupsilonP;
 193
 194 //------------------------------------------------------------------
 195 // Upsilon 3S generator
 196     AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF scaled", "UpsilonPP");
 197 // first step: generation in 4pi
 198     genupsilonPP->SetPtRange(0.,100.);
 199     genupsilonPP->SetYRange(-8.,8.);
 200     genupsilonPP->SetPhiRange(0.,360.);
 201     genupsilonPP->SetForceDecay(kAll);
 202     genupsilonPP->Init();  // generation in 4pi
 203     ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 204     AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
 205     AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
 206 // second step: generation in selected kinematical range
 207     genupsilonPP->SetPtRange(ptMin, ptMax);
 208     genupsilonPP->SetYRange(yMin, yMax);
 209     genupsilonPP->SetPhiRange(phiMin, phiMax);
 210     genupsilonPP->Init(); // generation in selected kinematical range
 211     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
 212     fTotalRate+=ratioupsilonPP;
 213
 214 //------------------------------------------------------------------
 215 // Pythia generator
 216     AliGenPythia *pythia = new AliGenPythia(1);
 217     pythia->SetProcess(kPyMbMSEL1);
 218     pythia->SetStrucFunc(kCTEQ5L);
 219     pythia->SetEnergyCMS(14000.);
 220     pythia->SetForceDecay(kAll);
 221     pythia->SetPtRange(0.,100.);
 222     pythia->SetYRange(-8.,8.);
 223     pythia->SetPhiRange(0.,360.);
 224     pythia->SetPtHard(2.76,-1.0);
 225     pythia->Init();
 226     AddGenerator(pythia,"Pythia",1);
 227     fTotalRate+=1.;
 228
 229 }
 230
 231 //_________________________________________________________________________
 232 void AliGenMUONCocktailpp::Generate()
 233 {
 234 // Generate event
 235     TIter next(fEntries);
 236     AliGenCocktailEntry *entry = 0;
 237     AliGenCocktailEntry *preventry = 0;
 238     AliGenerator* gen = 0;
 239
 240     TObjArray *partArray = gAlice->GetMCApp()->Particles();
 241
 242 // Generate the vertex position used by all generators
 243     if(fVertexSmear == kPerEvent) Vertex();
 244
 245 // Loop on primordialTrigger:
 246 // minimum muon multiplicity above a pt cut in a theta acceptance region
 247
 248     Bool_t primordialTrigger = kFALSE;
 249     while(!primordialTrigger) {
 250         //Reseting stack
 251         AliRunLoader * runloader = gAlice->GetRunLoader();
 252         if (runloader)
 253             if (runloader->Stack())
 254                 runloader->Stack()->Reset();
 255         // Loop over generators and generate events
 256         Int_t igen = 0;
 257         Int_t npart = 0;
 258         const char* genName = 0;
 259         while((entry = (AliGenCocktailEntry*)next())) {
 260             gen = entry->Generator();
 261             genName = entry->GetName();
 262             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 263
 264             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
 265                 gRandom->Poisson(entry->Rate());
 266
 267             if (npart > 0) {
 268                 igen++;
 269                 if (igen == 1) entry->SetFirst(0);
 270                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 271
 272                 gen->SetNumberParticles(npart);
 273                 gen->Generate();
 274                 entry->SetLast(partArray->GetEntriesFast());
 275                 preventry = entry;
 276             }
 277         }
 278         next.Reset();
 279
 280 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
 281 // in the muon spectrometer acceptance
 282         Int_t iPart;
 283         fNGenerated++;
 284         Int_t numberOfMuons=0;
 285
 286         for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
 287
 288             if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
 289                  (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 290                  (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 291                  (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) {
 292
 293                 numberOfMuons++;
 294             }
 295         }
 296
 297         if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
 298     }
 299     fNSucceded++;
 300
 301     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 302     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 303 }
 304
 305