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