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