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