EVGEN/AliGenMUONCocktail.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 MUON coktail for physics in the Alice muon spectrometer
  20 // The followoing muons sources are included in this cocktail:
  21 //     jpsi, upsilon, non-correlated open and beauty, and muons from pion and kaons.
  22 // The free parameeters are :
  23 //      pp reaction cross-section
  24 //      production cross-sections in pp collisions and
  25 //      branching ratios in the muon channel
  26 // Hard probes are supposed to scale with Ncoll and hadronic production with (0.8Ncoll+0.2*Npart)
  27 // There is a primordial trigger wiche requires :
  28 //      a minimum muon multiplicity above a pT cut in a theta acceptance cone
  29 //
  30 // Gines Martinez, jan 2004, Nantes  martinez@in2p3.fr
  31
  32
  33 //
  34
  35 //#include <TList.h>
  36 #include <TObjArray.h>
  37 #include <TF1.h>
  38 #include <TParticle.h>
  39
  40 #include "AliGenParam.h"
  41 #include "AliGenMUONlib.h"
  42 #include "AliGenMUONCocktail.h"
  43 #include "AliGenCocktailEntry.h"
  44 #include "AliCollisionGeometry.h"
  45 #include "AliRun.h"
  46 #include "AliMC.h"
  47 #include "AliStack.h"
  48
  49 ClassImp(AliGenMUONCocktail)
  50
  51
  52 //________________________________________________________________________
  53 AliGenMUONCocktail::AliGenMUONCocktail()
  54   :AliGenCocktail()
  55 {
  56 // Constructor
  57   fTotalRate =0;
  58   fNSucceded=0;
  59   fNGenerated=0;
  60   fMuonMultiplicity=1;
  61   fMuonPtCut= 1.;
  62   fMuonThetaMinCut=171.;
  63   fMuonThetaMaxCut=178.;
  64   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
  65     //
  66 }
  67 //_________________________________________________________________________
  68 AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
  69     AliGenCocktail(cocktail)
  70 {
  71 // Copy constructor
  72   fTotalRate =0;
  73   fNSucceded=0;
  74   fNGenerated=0;
  75   fMuonMultiplicity=1;
  76   fMuonPtCut= 1.;
  77   fMuonThetaMinCut=171.;
  78   fMuonThetaMaxCut=178.;
  79   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
  80     //
  81 }
  82 //_________________________________________________________________________
  83 AliGenMUONCocktail::~AliGenMUONCocktail()
  84 {
  85 // Destructor
  86     delete fEntries;
  87 }
  88
  89 //_________________________________________________________________________
  90 void AliGenMUONCocktail::Init()
  91 {
  92   // Defining MUON physics cocktail
  93
  94   // Kinematical limits for particle generation
  95   Float_t ptMin  = fPtMin;
  96   Float_t ptMax  = fPtMax;
  97   Float_t yMin   = fYMin;;
  98   Float_t yMax   = fYMax;;
  99   Float_t phiMin = fPhiMin*180./TMath::Pi();
 100   Float_t phiMax = fPhiMax*180./TMath::Pi();
 101   printf(">>> Kinematical range pT:%f:%f  y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 102
 103   Float_t sigmaReaction =   0.072;   // MINB pp at LHC energies 72 mb
 104
 105   // Generating J/Psi Physics
 106   AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
 107   // 4pi Generation
 108   genjpsi->SetPtRange(0,100.);
 109   genjpsi->SetYRange(-8.,8);
 110   genjpsi->SetPhiRange(0.,360.);
 111   genjpsi->SetForceDecay(kDiMuon);
 112   genjpsi->SetTrackingFlag(1);
 113   // Calculation of the paritcle multiplicity per event in the muonic channel
 114   Float_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 115   Float_t sigmajpsi     = 31.0e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 116   Float_t brjpsi        = 0.0588;           // Branching Ratio for JPsi
 117   genjpsi->Init();  // Generating pT and Y parametrsation for the 4pi
 118   ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 119   printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratiojpsi,genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmajpsi );
 120   // Generation in the kinematical limits of AliGenMUONCocktail
 121   genjpsi->SetPtRange(ptMin, ptMax);
 122   genjpsi->SetYRange(yMin, yMax);
 123   genjpsi->SetPhiRange(phiMin, phiMax);
 124   genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 125   // Adding Generator
 126   AddGenerator(genjpsi, "Jpsi", ratiojpsi);
 127   fTotalRate+=ratiojpsi;
 128
 129  // Generating Psi prime Physics
 130   AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "Vogt", "PsiP");
 131   // 4pi Generation
 132   genpsiP->SetPtRange(0,100.);
 133   genpsiP->SetYRange(-8.,8);
 134   genpsiP->SetPhiRange(0.,360.);
 135   genpsiP->SetForceDecay(kDiMuon);
 136   genpsiP->SetTrackingFlag(1);
 137   // Calculation of the paritcle multiplicity per event in the muonic channel
 138   Float_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 139   Float_t sigmapsiP     = 4.68e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 140   Float_t brpsiP        = 0.0103;           // Branching Ratio for PsiP
 141   genpsiP->Init();  // Generating pT and Y parametrsation for the 4pi
 142   ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 143   printf(">>> ratio psiP %g et %g Ncol %g sigma %g\n",ratiopsiP,genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmapsiP );
 144   // Generation in the kinematical limits of AliGenMUONCocktail
 145   genpsiP->SetPtRange(ptMin, ptMax);
 146   genpsiP->SetYRange(yMin, yMax);
 147   genpsiP->SetPhiRange(phiMin, phiMax);
 148   genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 149   // Adding Generator
 150   AddGenerator(genpsiP, "PsiP", ratiopsiP);
 151   fTotalRate+=ratiopsiP;
 152
 153
 154 // Generating Upsilon Physics
 155   AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");
 156   genupsilon->SetPtRange(0,100.);
 157   genupsilon->SetYRange(-8.,8);
 158   genupsilon->SetPhiRange(0.,360.);
 159   genupsilon->SetForceDecay(kDiMuon);
 160   genupsilon->SetTrackingFlag(1);
 161   Float_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 162   Float_t sigmaupsilon     = 0.501e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 163   Float_t brupsilon        = 0.0248;  // Branching Ratio for Upsilon
 164   genupsilon->Init();  // Generating pT and Y parametrsation for the 4pi
 165   ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 166   printf(">>> ratio upsilon %g et %g\n",ratioupsilon, genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
 167   genupsilon->SetPtRange(ptMin, ptMax);
 168   genupsilon->SetYRange(yMin, yMax);
 169   genupsilon->SetPhiRange(phiMin, phiMax);
 170   genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 171   AddGenerator(genupsilon,"Upsilon", ratioupsilon);
 172   fTotalRate+=ratioupsilon;
 173
 174 // Generating UpsilonP Physics
 175   AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "Vogt", "UpsilonP");
 176   genupsilonP->SetPtRange(0,100.);
 177   genupsilonP->SetYRange(-8.,8);
 178   genupsilonP->SetPhiRange(0.,360.);
 179   genupsilonP->SetForceDecay(kDiMuon);
 180   genupsilonP->SetTrackingFlag(1);
 181   Float_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 182   Float_t sigmaupsilonP     = 0.246e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 183   Float_t brupsilonP        = 0.0131;  // Branching Ratio for UpsilonP
 184   genupsilonP->Init();  // Generating pT and Y parametrsation for the 4pi
 185   ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 186   printf(">>> ratio upsilonP %g et %g\n",ratioupsilonP, genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
 187   genupsilonP->SetPtRange(ptMin, ptMax);
 188   genupsilonP->SetYRange(yMin, yMax);
 189   genupsilonP->SetPhiRange(phiMin, phiMax);
 190   genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 191   AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
 192   fTotalRate+=ratioupsilonP;
 193
 194
 195 // Generating UpsilonPP Physics
 196   AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "Vogt", "UpsilonPP");
 197   genupsilonPP->SetPtRange(0,100.);
 198   genupsilonPP->SetYRange(-8.,8);
 199   genupsilonPP->SetPhiRange(0.,360.);
 200   genupsilonPP->SetForceDecay(kDiMuon);
 201   genupsilonPP->SetTrackingFlag(1);
 202   Float_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 203   Float_t sigmaupsilonPP     = 0.100e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 204   Float_t brupsilonPP        = 0.0181;  // Branching Ratio for UpsilonPP
 205   genupsilonPP->Init();  // Generating pT and Y parametrsation for the 4pi
 206   ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 207   printf(">>> ratio upsilonPP %g et %g\n",ratioupsilonPP, genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
 208   genupsilonPP->SetPtRange(ptMin, ptMax);
 209   genupsilonPP->SetYRange(yMin, yMax);
 210   genupsilonPP->SetPhiRange(phiMin, phiMax);
 211   genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 212   AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
 213   fTotalRate+=ratioupsilonPP;
 214
 215
 216 // Generating Charm Physics
 217   AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");
 218   gencharm->SetPtRange(0,100.);
 219   gencharm->SetYRange(-8.,8);
 220   gencharm->SetPhiRange(0.,360.);
 221   gencharm->SetForceDecay(kSemiMuonic);
 222   gencharm->SetTrackingFlag(1);
 223   Float_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 224   Float_t sigmacharm     = 2. * 6.64e-3 * 0.65 ;   //
 225   Float_t brcharm        = 0.12;  // Branching Ratio for Charm
 226   gencharm->Init();  // Generating pT and Y parametrsation for the 4pi
 227   ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction *
 228     gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 229   gencharm->SetPtRange(ptMin, ptMax);
 230   gencharm->SetYRange(yMin, yMax);
 231   gencharm->SetPhiRange(phiMin, phiMax);
 232   gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 233
 234   printf(">>> ratio charm %f\n",ratiocharm);
 235   AddGenerator(gencharm,"Charm", ratiocharm);
 236   fTotalRate+=ratiocharm;
 237
 238 // Generating Beauty Physics "Correlated Pairs"
 239   AliGenParam * genbeauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");
 240   genbeauty->SetPtRange(0,100.);
 241   genbeauty->SetYRange(-8.,8);
 242   genbeauty->SetPhiRange(0.,360.);
 243   genbeauty->SetForceDecay(kSemiMuonic);
 244   genbeauty->SetTrackingFlag(1);
 245   Float_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 246   Float_t sigmabeauty     = 2. * 0.210e-3 * 0.84;   //
 247   Float_t brbeauty        = 0.15;  // Branching Ratio for Beauty
 248   genbeauty->Init();  // Generating pT and Y parametrsation for the 4pi
 249   ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction *
 250     genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 251   genbeauty->SetPtRange(ptMin, ptMax);
 252   genbeauty->SetYRange(yMin, yMax);
 253   genbeauty->SetPhiRange(phiMin, phiMax);
 254   genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
 255
 256   printf(">>> ratio beauty %f\n",ratiobeauty);
 257   AddGenerator(genbeauty,"Beauty", ratiobeauty);
 258   fTotalRate+=ratiobeauty;
 259
 260 // Generating Pion Physics
 261   AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");
 262   genpion->SetPtRange(0,100.);
 263   genpion->SetYRange(-8.,8);
 264   genpion->SetPhiRange(0.,360.);
 265   genpion->SetForceDecay(kPiToMu);
 266   genpion->SetTrackingFlag(1);
 267   Float_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 268   Float_t sigmapion     = 0.93e-2; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
 269   Float_t brpion        = 0.9999;  // Branching Ratio for Pion
 270   genpion->Init();  // Generating pT and Y parametrsation for the 4pi
 271   ratiopion = sigmapion * brpion *  (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 272   genpion->SetPtRange(ptMin, ptMax);
 273   genpion->SetYRange(yMin, yMax);
 274   genpion->SetPhiRange(phiMin, phiMax);
 275   genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 276   printf(">>> ratio pion %f\n",ratiopion);
 277   AddGenerator(genpion,"Pion", ratiopion);
 278   fTotalRate+=ratiopion;
 279
 280 // Generating Kaon Physics
 281   AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");
 282   genkaon->SetPtRange(0,100.);
 283   genkaon->SetYRange(-8.,8);
 284   genkaon->SetPhiRange(0.,360.);
 285   genkaon->SetForceDecay(kKaToMu);
 286   genkaon->SetTrackingFlag(1);
 287   Float_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 288   Float_t sigmakaon     = 1.23e-4;   // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
 289   Float_t brkaon        = 0.6351 ;  // Branching Ratio for Kaon
 290   genkaon->Init();  // Generating pT and Y parametrsation for the 4pi
 291   ratiokaon = sigmakaon * brkaon * (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 292   genkaon->SetPtRange(ptMin, ptMax);
 293   genkaon->SetYRange(yMin, yMax);
 294   genkaon->SetPhiRange(phiMin, phiMax);
 295   genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 296   printf(">>> ratio kaon %f\n",ratiokaon);
 297   AddGenerator(genkaon,"Kaon", ratiokaon);
 298   fTotalRate+=ratiokaon;
 299 }
 300
 301 //_________________________________________________________________________
 302 void AliGenMUONCocktail::Generate()
 303 {
 304   //
 305 // Generate event
 306     TIter next(fEntries);
 307     AliGenCocktailEntry *entry = 0;
 308     AliGenCocktailEntry *preventry = 0;
 309     AliGenerator* gen = 0;
 310
 311     TObjArray *partArray = gAlice->GetMCApp()->Particles();
 312
 313 //
 314 //  Generate the vertex position used by all generators
 315 //
 316     if(fVertexSmear == kPerEvent) Vertex();
 317     Bool_t primordialTrigger = kFALSE;
 318
 319     while(!primordialTrigger) {
 320
 321       //Reseting stack
 322       AliRunLoader * runloader = gAlice->GetRunLoader();
 323       if (runloader)
 324         if (runloader->Stack())
 325           runloader->Stack()->Reset();
 326       //
 327       // Loop over generators and generate events
 328       Int_t igen=0;
 329       Int_t npart =0;
 330
 331       while((entry = (AliGenCocktailEntry*)next())) {
 332         gen = entry->Generator();
 333         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 334         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
 335           igen++;
 336           if (igen ==1) entry->SetFirst(0);
 337           else  entry->SetFirst((partArray->GetEntriesFast())+1);
 338           gen->SetNumberParticles(npart);
 339           gen->Generate();
 340           entry->SetLast(partArray->GetEntriesFast());
 341           preventry = entry;
 342         }
 343       }
 344       next.Reset();
 345
 346       // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance 171,178 and pTCut
 347       Int_t iPart;
 348       fNGenerated++;
 349       Int_t numberOfMuons=0;
 350       //      printf(">>>fNGenerated is %d\n",fNGenerated);
 351       for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
 352         //      gAlice->GetMCApp()->Particle(iPart)->Print();
 353         if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
 354              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 355              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 356              (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) {
 357           gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);
 358           numberOfMuons++;
 359         }
 360       }
 361       //  printf(">>> Number of Muons is %d \n", numberOfMuons);
 362       if (numberOfMuons >= fMuonMultiplicity ) primordialTrigger = kTRUE;
 363     }
 364     //printf(">>> Trigger Accepted!!! \n");
 365     fNSucceded++;
 366     //   Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
 367     //    printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
 368 }
 369
 370
 371
 372
 373
 374