1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 // 2009: added possibility to hide x-sections (B. Vulpescu)
29 // 2009: added possibility to have the cocktail (fast generator and param.)
30 // for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
31 // 2009: added polarization (L. Bianchi)
33 #include <TObjArray.h>
34 #include <TParticle.h>
36 #include <TVirtualMC.h>
38 #include "AliGenCocktailEventHeader.h"
40 #include "AliGenCocktailEntry.h"
41 #include "AliGenMUONCocktailpp.h"
42 #include "AliGenMUONlib.h"
43 #include "AliGenParam.h"
47 #include "AliDecayer.h"
49 #include "AliGenCorrHF.h"
50 #include "AliDecayerPolarized.h"
52 ClassImp(AliGenMUONCocktailpp)
54 //________________________________________________________________________
55 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
58 fDecayModeResonance(kAll),
59 fDecayModePythia(kAll),
64 fMuonThetaMaxCut(180.),
65 fMuonOriginCut(-999.),
76 //x-sections for pp @ 10 TeV
78 fSigmaReaction(0.0695),
81 fSigmaUpsilon(0.658e-6),
82 fSigmaUpsilonP(0.218e-6),
83 fSigmaUpsilonPP(0.122e-6),
87 //x-sections for pp @ 14 TeV
89 fSigmaReaction(0.070),
92 fSigmaUpsilon(0.989e-6),
93 fSigmaUpsilonP(0.502e-6),
94 fSigmaUpsilonPP(0.228e-6),
96 fSigmaBBbar(0.51e-3),*/
103 //_________________________________________________________________________
104 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
110 //_________________________________________________________________________
111 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
112 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
114 if (strcmp(PolFrame,"kColSop")==0){
115 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
116 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
117 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
118 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
119 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
121 } else if (strcmp(PolFrame,"kHelicity")==0){
122 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
123 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
124 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
125 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
126 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
130 AliInfo(Form("The polarization frame is not valid"));
131 AliInfo(Form("No polarization will be set"));
140 //_________________________________________________________________________
141 void AliGenMUONCocktailpp::CreateCocktail()
144 // These limits are only used for renormalization of quarkonia cross section
145 // Pythia events are generated in 4pi
146 Double_t ptMin = fPtMin;
147 Double_t ptMax = fPtMax;
148 Double_t yMin = fYMin;;
149 Double_t yMax = fYMax;;
150 Double_t phiMin = fPhiMin*180./TMath::Pi();
151 Double_t phiMax = fPhiMax*180./TMath::Pi();
152 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));
154 // Ratios with respect to the reaction cross-section in the
155 // kinematics limit of the MUONCocktail
158 Double_t ratioupsilon;
159 Double_t ratioupsilonP;
160 Double_t ratioupsilonPP;
165 Int_t cmsEnergy = fCMSEnergy;
167 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
168 // corrected from feed down of higher resonances
170 Double_t sigmajpsi = fSigmaJPsi;
171 Double_t sigmapsiP = fSigmaPsiP;
172 Double_t sigmaupsilon = fSigmaUpsilon;
173 Double_t sigmaupsilonP = fSigmaUpsilonP;
174 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
175 Double_t sigmaccbar = fSigmaCCbar;
176 Double_t sigmabbbar = fSigmaBBbar;
178 // Cross sections corrected with the BR in mu+mu-
179 // (only in case of use of AliDecayerPolarized)
181 if(fJpsiPol != 0) {sigmajpsi = fSigmaJPsi*0.0593;}
182 if(fPsiPPol != 0) {sigmapsiP = fSigmaPsiP*0.0075;}
183 if(fUpsPol != 0) {sigmaupsilon = fSigmaUpsilon*0.0248;}
184 if(fUpsPPol != 0) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
185 if(fUpsPPPol != 0) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
187 // MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
188 Double_t sigmaReaction = fSigmaReaction;
190 Int_t eincStart = 10;
192 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
194 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
195 //----------------------------------------------------------------------
197 AliGenParam * genjpsi;
198 if(cmsEnergy == eincStart){
199 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
201 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
203 // first step: generation in 4pi
204 genjpsi->SetPtRange(0.,100.);
205 genjpsi->SetYRange(-8.,8.);
206 genjpsi->SetPhiRange(0.,360.);
207 genjpsi->SetForceDecay(fDecayModeResonance);
208 //if (!gMC) genjpsi->SetDecayer(fDecayer);
210 genjpsi->Init(); // generation in 4pi
211 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
213 AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
214 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
216 // second step: generation in selected kinematical range
217 genjpsi->SetPtRange(ptMin, ptMax);
218 genjpsi->SetYRange(yMin, yMax);
219 genjpsi->SetPhiRange(phiMin, phiMax);
220 genjpsi->Init(); // generation in selected kinematic range
222 TVirtualMCDecayer *Jpsidec = 0;
224 AliInfo(Form("******Setting polarized decayer for J/psi"));
226 Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
227 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
230 Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
231 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
233 Jpsidec->SetForceDecay(kAll);
235 genjpsi->SetDecayer(Jpsidec);
238 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
239 //------------------------------------------------------------------
240 // Psi prime generator
241 AliGenParam * genpsiP;
242 if(cmsEnergy == eincStart){
243 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
245 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
247 // first step: generation in 4pi
248 genpsiP->SetPtRange(0.,100.);
249 genpsiP->SetYRange(-8.,8.);
250 genpsiP->SetPhiRange(0.,360.);
251 genpsiP->SetForceDecay(fDecayModeResonance);
252 //if (!gMC) genpsiP->SetDecayer(fDecayer);
254 genpsiP->Init(); // generation in 4pi
255 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
257 AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
258 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
260 // second step: generation in selected kinematical range
261 genpsiP->SetPtRange(ptMin, ptMax);
262 genpsiP->SetYRange(yMin, yMax);
263 genpsiP->SetPhiRange(phiMin, phiMax);
264 genpsiP->Init(); // generation in selected kinematic range
266 TVirtualMCDecayer *PsiPdec = 0;
268 AliInfo(Form("******Setting polarized decayer for psi'"));
270 PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
271 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
274 PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
275 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
277 PsiPdec->SetForceDecay(kAll);
279 genpsiP->SetDecayer(PsiPdec);
282 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
283 //------------------------------------------------------------------
284 // Upsilon 1S generator
285 AliGenParam * genupsilon;
286 if(cmsEnergy == eincStart) {
287 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
289 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
291 // first step: generation in 4pi
292 genupsilon->SetPtRange(0.,100.);
293 genupsilon->SetYRange(-8.,8.);
294 genupsilon->SetPhiRange(0.,360.);
295 genupsilon->SetForceDecay(fDecayModeResonance);
296 //if (!gMC) genupsilon->SetDecayer(fDecayer);
297 genupsilon->Init(); // generation in 4pi
298 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
300 AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
301 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
303 // second step: generation in selected kinematical range
304 genupsilon->SetPtRange(ptMin, ptMax);
305 genupsilon->SetYRange(yMin, yMax);
306 genupsilon->SetPhiRange(phiMin, phiMax);
307 genupsilon->Init(); // generation in selected kinematical range
309 TVirtualMCDecayer *Upsdec = 0;
311 AliInfo(Form("******Setting polarized decayer for Upsilon"));
313 Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
314 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
317 Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
318 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
320 Upsdec->SetForceDecay(kAll);
322 genupsilon->SetDecayer(Upsdec);
325 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
326 //------------------------------------------------------------------
327 // Upsilon 2S generator
328 AliGenParam * genupsilonP;
329 if(cmsEnergy == eincStart){
330 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
332 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
335 // first step: generation in 4pi
336 genupsilonP->SetPtRange(0.,100.);
337 genupsilonP->SetYRange(-8.,8.);
338 genupsilonP->SetPhiRange(0.,360.);
339 genupsilonP->SetForceDecay(fDecayModeResonance);
340 //if (!gMC) genupsilonP->SetDecayer(fDecayer);
341 genupsilonP->Init(); // generation in 4pi
342 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
344 AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
345 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
347 // second step: generation in the kinematical range
348 genupsilonP->SetPtRange(ptMin, ptMax);
349 genupsilonP->SetYRange(yMin, yMax);
350 genupsilonP->SetPhiRange(phiMin, phiMax);
351 genupsilonP->Init(); // generation in selected kinematical range
353 TVirtualMCDecayer *UpsPdec = 0;
355 AliInfo(Form("******Setting polarized decayer for Upsilon'"));
357 UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
358 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
361 UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
362 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
364 UpsPdec->SetForceDecay(kAll);
366 genupsilonP->SetDecayer(UpsPdec);
369 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
371 //------------------------------------------------------------------
372 // Upsilon 3S generator
373 AliGenParam * genupsilonPP;
374 if(cmsEnergy == eincStart){
375 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
378 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
381 // first step: generation in 4pi
382 genupsilonPP->SetPtRange(0.,100.);
383 genupsilonPP->SetYRange(-8.,8.);
384 genupsilonPP->SetPhiRange(0.,360.);
385 genupsilonPP->SetForceDecay(fDecayModeResonance);
386 //if (!gMC) genupsilonPP->SetDecayer(fDecayer);
387 genupsilonPP->Init(); // generation in 4pi
388 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
390 AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
391 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
393 // second step: generation in selected kinematical range
394 genupsilonPP->SetPtRange(ptMin, ptMax);
395 genupsilonPP->SetYRange(yMin, yMax);
396 genupsilonPP->SetPhiRange(phiMin, phiMax);
397 genupsilonPP->Init(); // generation in selected kinematical range
399 TVirtualMCDecayer *UpsPPdec = 0;
401 AliInfo(Form("******Setting polarized decayer for Upsilon''"));
403 UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
404 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
407 UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
408 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
410 UpsPPdec->SetForceDecay(kAll);
412 genupsilonPP->SetDecayer(UpsPPdec);
415 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
417 //------------------------------------------------------------------
418 // Generator of charm
420 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
421 gencharm->SetMomentumRange(0,9999);
422 gencharm->SetForceDecay(kAll);
423 ratioccbar = sigmaccbar/sigmaReaction;
424 //if (!gMC) gencharm->SetDecayer(fDecayer);
427 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
428 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
430 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
432 //------------------------------------------------------------------
433 // Generator of beauty
435 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
436 genbeauty->SetMomentumRange(0,9999);
437 genbeauty->SetForceDecay(kAll);
438 ratiobbbar = sigmabbbar/sigmaReaction;
439 //if (!gMC) genbeauty->SetDecayer(fDecayer);
442 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
443 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
445 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
447 //-------------------------------------------------------------------
450 // This has to go into the Config.C
452 // AliGenPythia *pythia = new AliGenPythia(1);
453 // pythia->SetProcess(kPyMbMSEL1);
454 // pythia->SetStrucFunc(kCTEQ5L);
455 // pythia->SetEnergyCMS(14000.);
456 // AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
457 // Decay_t dt = gener->GetDecayModePythia();
458 // pythia->SetForceDecay(dt);
459 // pythia->SetPtRange(0.,100.);
460 // pythia->SetYRange(-8.,8.);
461 // pythia->SetPhiRange(0.,360.);
462 // pythia->SetPtHard(2.76,-1.0);
463 // pythia->SwitchHFOff();
465 // AddGenerator(pythia,"Pythia",1);
468 void AliGenMUONCocktailpp::Init()
472 TIter next(fEntries);
473 AliGenCocktailEntry *entry;
475 while((entry = (AliGenCocktailEntry*)next())) {
476 entry->Generator()->SetStack(fStack);
481 //_________________________________________________________________________
482 void AliGenMUONCocktailpp::Generate()
485 TIter next(fEntries);
486 AliGenCocktailEntry *entry = 0;
487 AliGenCocktailEntry *preventry = 0;
488 AliGenerator* gen = 0;
490 if (fHeader) delete fHeader;
491 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
493 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
495 // Generate the vertex position used by all generators
496 if(fVertexSmear == kPerEvent) Vertex();
498 // Loop on primordialTrigger:
499 // minimum muon multiplicity above a pt cut in a theta acceptance region
501 Bool_t primordialTrigger = kFALSE;
502 while(!primordialTrigger) {
504 AliRunLoader * runloader = AliRunLoader::Instance();
506 if (runloader->Stack())
507 runloader->Stack()->Clean();
508 // Loop over generators and generate events
511 const char* genName = 0;
512 while((entry = (AliGenCocktailEntry*)next())) {
513 gen = entry->Generator();
514 genName = entry->GetName();
515 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
517 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
518 gRandom->Poisson(entry->Rate());
522 if (igen == 1) entry->SetFirst(0);
523 else entry->SetFirst((partArray->GetEntriesFast())+1);
525 gen->SetNumberParticles(npart);
527 entry->SetLast(partArray->GetEntriesFast());
533 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
534 // in the muon spectrometer acceptance
537 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
538 for(iPart=0; iPart<maxPart; iPart++){
540 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
541 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
542 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
543 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
544 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
545 (part->Pt()>fMuonPtCut) &&
546 (part->P()>fMuonPCut)) {
551 if (numberOfMuons >= fMuonMultiplicity) {
552 primordialTrigger = kTRUE;
553 fHeader->SetNProduced(maxPart);
561 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
563 fHeader->SetPrimaryVertex(eventVertex);
565 gAlice->SetGenEventHeader(fHeader);
567 // AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
568 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));