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