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