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