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 // Class to create the coktail for physics with muons for pp collisions
20 // using the The followoing sources:
21 // jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
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)
32 // 2009: added polarization (L. Bianchi)
34 // 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid).
35 // Cross-sections for charmonia are now directly taken from the Yellow Report
36 // (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy.
37 // usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
38 //------------------------
40 // - CMS energy passed via parameter
41 // i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C
42 // - resonances now added to the cocktail via AddReso2Generator
44 // B.Vulpescu & P.Crochet
46 #include <TObjArray.h>
47 #include <TParticle.h>
49 #include <TVirtualMC.h>
51 #include "AliGenCocktailEventHeader.h"
53 #include "AliGenCocktailEntry.h"
54 #include "AliGenMUONCocktailpp.h"
55 #include "AliGenMUONlib.h"
56 #include "AliGenParam.h"
60 #include "AliDecayer.h"
62 #include "AliGenCorrHF.h"
63 #include "AliDecayerPolarized.h"
65 ClassImp(AliGenMUONCocktailpp)
67 //________________________________________________________________________
68 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
71 fDecayModeResonance(kAll),
72 fDecayModePythia(kAll),
77 fMuonThetaMaxCut(180.),
78 fMuonOriginCut(-999.),
94 fSigmaReactionArray(),
104 fSigmaUpsilonArray(),
106 fSigmaUpsilonPArray(),
108 fSigmaUpsilonPPArray(),
118 // x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
119 // bottomnium as for 10 TeV
120 fCMSEnergyTeVArray[0] = 7.00;
121 fSigmaReactionArray[0] = 0.0695;
122 fSigmaJPsiArray[0] = 21.8e-6;
123 fSigmaChic1Array[0] = 21.1e-6;
124 fSigmaChic2Array[0] = 34.9e-6;
125 fSigmaPsiPArray[0] = 4.93e-6;
126 fSigmaUpsilonArray[0] = 0.463e-6;
127 fSigmaUpsilonPArray[0] = 0.154e-6;
128 fSigmaUpsilonPPArray[0] = 0.0886e-6;
129 fSigmaCCbarArray[0] = 6.91e-3;
130 fSigmaBBbarArray[0] = 0.232e-3;
132 //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
133 // scaled down according to ccbar and bbar cross-sections
134 fCMSEnergyTeVArray[1] = 10.00;
135 fSigmaReactionArray[1] = 0.0695;
136 fSigmaJPsiArray[1] = 26.06e-6;
137 fSigmaChic1Array[1] = 25.18e-6;
138 fSigmaChic2Array[1] = 41.58e-6;
139 fSigmaPsiPArray[1] = 5.88e-6;
140 fSigmaUpsilonArray[1] = 0.658e-6;
141 fSigmaUpsilonPArray[1] = 0.218e-6;
142 fSigmaUpsilonPPArray[1] = 0.122e-6;
143 fSigmaCCbarArray[1] = 8.9e-3;
144 fSigmaBBbarArray[1] = 0.33e-3;
146 //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
147 // bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that
148 // feed-down from chib is included
149 fCMSEnergyTeVArray[2] = 14.00;
150 fSigmaReactionArray[2] = 0.070;
151 fSigmaJPsiArray[2] = 32.9e-6;
152 fSigmaChic1Array[2] = 31.8e-6;
153 fSigmaChic2Array[2] = 52.5e-6;
154 fSigmaPsiPArray[2] = 7.43e-6;
155 fSigmaUpsilonArray[2] = 0.989e-6;
156 fSigmaUpsilonPArray[2] = 0.502e-6;
157 fSigmaUpsilonPPArray[2] = 0.228e-6;
158 fSigmaCCbarArray[2] = 11.2e-3;
159 fSigmaBBbarArray[2] = 0.51e-3;
163 //_________________________________________________________________________
164 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
170 //_________________________________________________________________________
171 void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
173 // setter for CMSEnergy and corresponding cross-sections
174 fCMSEnergyTeV = fCMSEnergyTeVArray[cmsEnergy];
175 fSigmaReaction = fSigmaReactionArray[cmsEnergy];
176 fSigmaJPsi = fSigmaJPsiArray[cmsEnergy];
177 fSigmaChic1 = fSigmaChic1Array[cmsEnergy];
178 fSigmaChic2 = fSigmaChic2Array[cmsEnergy];
179 fSigmaPsiP = fSigmaPsiPArray[cmsEnergy];
180 fSigmaUpsilon = fSigmaUpsilonArray[cmsEnergy];
181 fSigmaUpsilonP = fSigmaUpsilonPArray[cmsEnergy];
182 fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy];
183 fSigmaCCbar = fSigmaCCbarArray[cmsEnergy];
184 fSigmaBBbar = fSigmaBBbarArray[cmsEnergy];
187 //_________________________________________________________________________
188 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
189 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
190 // setter for resonances polarization
191 if (strcmp(PolFrame,"kColSop")==0){
192 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
193 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
194 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
195 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
196 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
198 } else if (strcmp(PolFrame,"kHelicity")==0){
199 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
200 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
201 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
202 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
203 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
207 AliInfo(Form("The polarization frame is not valid"));
208 AliInfo(Form("No polarization will be set"));
217 //_________________________________________________________________________
218 void AliGenMUONCocktailpp::CreateCocktail()
220 // create and add resonances and open HF to the coctail
221 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
223 // These limits are only used for renormalization of quarkonia cross section
224 // Pythia events are generated in 4pi
225 Double_t ptMin = fPtMin;
226 Double_t ptMax = fPtMax;
227 Double_t yMin = fYMin;;
228 Double_t yMax = fYMax;;
229 Double_t phiMin = fPhiMin*180./TMath::Pi();
230 Double_t phiMax = fPhiMax*180./TMath::Pi();
231 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));
233 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
234 // corrected from feed down of higher resonances
236 Double_t sigmajpsi = fSigmaJPsi;
237 Double_t sigmachic1 = fSigmaChic1;
238 Double_t sigmachic2 = fSigmaChic2;
239 Double_t sigmapsiP = fSigmaPsiP;
240 Double_t sigmaupsilon = fSigmaUpsilon;
241 Double_t sigmaupsilonP = fSigmaUpsilonP;
242 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
243 Double_t sigmaccbar = fSigmaCCbar;
244 Double_t sigmabbbar = fSigmaBBbar;
246 // Cross sections corrected with the BR in mu+mu-
247 // (only in case of use of AliDecayerPolarized)
249 if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;}
250 if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent
251 if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent
252 if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;}
253 if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;}
254 if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
255 if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
257 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
259 // Create and add resonances to the generator
260 AliGenParam * genjpsi=0;
261 AliGenParam * genchic1=0;
262 AliGenParam * genchic2=0;
263 AliGenParam * genpsiP=0;
264 AliGenParam * genupsilon=0;
265 AliGenParam * genupsilonP=0;
266 AliGenParam * genupsilonPP=0;
269 Char_t nameChic1[10];
270 Char_t nameChic2[10];
274 Char_t nameUpsPP[10];
276 sprintf(nameJpsi,"Jpsi");
277 sprintf(nameChic1,"Chic1");
278 sprintf(nameChic2,"Chic2");
279 sprintf(namePsiP,"PsiP");
280 sprintf(nameUps,"Ups");
281 sprintf(nameUpsP,"UpsP");
282 sprintf(nameUpsPP,"UpsPP");
285 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
286 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
287 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
288 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
289 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
291 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
292 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
293 } else if (cmsEnergy == 7){
294 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi");
295 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 7", "Chic1");
296 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 7", "Chic2");
297 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 7", "PsiP");
299 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 7", "Upsilon");
300 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 7", "UpsilonP");
301 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 7", "UpsilonPP");
302 } else if (cmsEnergy == 14){
303 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
304 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
305 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
306 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
308 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
309 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
311 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
314 AddReso2Generator(nameJpsi,genjpsi,sigmajpsi,fJpsiPol);
315 AddReso2Generator(nameChic1,genchic2,sigmachic1,fChic2Pol);
316 AddReso2Generator(nameChic2,genpsiP,sigmapsiP,fPsiPPol);
317 AddReso2Generator(namePsiP,genchic1,sigmachic1,fChic1Pol);
318 AddReso2Generator(nameUps,genupsilon,sigmaupsilon,fUpsPol);
319 AddReso2Generator(nameUpsP,genupsilonP,sigmaupsilonP,fUpsPPol);
320 AddReso2Generator(nameUpsPP,genupsilonPP,sigmaupsilonPP,fUpsPPPol);
322 //------------------------------------------------------------------
323 // Generator of charm
324 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
325 gencharm->SetMomentumRange(0,9999);
326 gencharm->SetForceDecay(kAll);
327 Double_t ratioccbar = sigmaccbar/fSigmaReaction;
328 if (!gMC) gencharm->SetDecayer(fDecayer);
331 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
332 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
334 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
335 //------------------------------------------------------------------
336 // Generator of beauty
337 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
338 genbeauty->SetMomentumRange(0,9999);
339 genbeauty->SetForceDecay(kAll);
340 Double_t ratiobbbar = sigmabbbar/fSigmaReaction;
341 if (!gMC) genbeauty->SetDecayer(fDecayer);
344 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
345 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
347 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
349 //-------------------------------------------------------------------
352 // This has to go into the Config.C
354 // AliGenPythia *pythia = new AliGenPythia(1);
355 // pythia->SetProcess(kPyMbMSEL1);
356 // pythia->SetStrucFunc(kCTEQ5L);
357 // pythia->SetEnergyCMS(14000.);
358 // AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
359 // Decay_t dt = gener->GetDecayModePythia();
360 // pythia->SetForceDecay(dt);
361 // pythia->SetPtRange(0.,100.);
362 // pythia->SetYRange(-8.,8.);
363 // pythia->SetPhiRange(0.,360.);
364 // pythia->SetPtHard(2.76,-1.0);
365 // pythia->SwitchHFOff();
367 // AddGenerator(pythia,"Pythia",1);
371 //-------------------------------------------------------------------
372 void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
373 AliGenParam* const genReso,
377 // add resonances to the cocktail
378 Double_t phiMin = fPhiMin*180./TMath::Pi();
379 Double_t phiMax = fPhiMax*180./TMath::Pi();
381 // first step: generation in 4pi
382 genReso->SetPtRange(0.,100.);
383 genReso->SetYRange(-8.,8.);
384 genReso->SetPhiRange(0.,360.);
385 genReso->SetForceDecay(fDecayModeResonance);
386 if (!gMC) genReso->SetDecayer(fDecayer);
387 genReso->Init(); // generation in 4pi
388 // Ratios with respect to the reaction cross-section in the
389 // kinematics limit of the MUONCocktail
390 Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
392 AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
393 AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
395 // second step: generation in selected kinematical range
396 genReso->SetPtRange(fPtMin, fPtMax);
397 genReso->SetYRange(fYMin, fYMax);
398 genReso->SetPhiRange(phiMin, phiMax);
399 genReso->Init(); // generation in selected kinematical range
401 TVirtualMCDecayer *decReso = 0;
402 if(TMath::Abs(polReso) > 1.e-30){
403 AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
405 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
406 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
409 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
410 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
412 decReso->SetForceDecay(kAll);
414 genReso->SetDecayer(decReso);
417 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
420 //-------------------------------------------------------------------
421 void AliGenMUONCocktailpp::Init()
424 TIter next(fEntries);
425 AliGenCocktailEntry *entry;
427 while((entry = (AliGenCocktailEntry*)next())) {
428 entry->Generator()->SetStack(fStack);
433 //_________________________________________________________________________
434 void AliGenMUONCocktailpp::Generate()
437 TIter next(fEntries);
438 AliGenCocktailEntry *entry = 0;
439 AliGenCocktailEntry *preventry = 0;
440 AliGenerator* gen = 0;
442 if (fHeader) delete fHeader;
443 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
445 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
447 // Generate the vertex position used by all generators
448 if(fVertexSmear == kPerEvent) Vertex();
450 // Loop on primordialTrigger:
451 // minimum muon multiplicity above a pt cut in a theta acceptance region
453 Bool_t primordialTrigger = kFALSE;
454 while(!primordialTrigger) {
456 AliRunLoader * runloader = AliRunLoader::Instance();
458 if (runloader->Stack())
459 runloader->Stack()->Clean();
460 // Loop over generators and generate events
463 const char* genName = 0;
464 while((entry = (AliGenCocktailEntry*)next())) {
465 gen = entry->Generator();
466 genName = entry->GetName();
467 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
469 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
470 gRandom->Poisson(entry->Rate());
474 if (igen == 1) entry->SetFirst(0);
475 else entry->SetFirst((partArray->GetEntriesFast())+1);
477 gen->SetNumberParticles(npart);
479 entry->SetLast(partArray->GetEntriesFast());
486 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
487 // in the muon spectrometer acceptance
490 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
491 for(iPart=0; iPart<maxPart; iPart++){
493 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
494 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
495 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
496 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
497 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
498 (part->Pt()>fMuonPtCut) &&
499 (part->P()>fMuonPCut)) {
504 if (numberOfMuons >= fMuonMultiplicity) {
505 primordialTrigger = kTRUE;
506 fHeader->SetNProduced(maxPart);
514 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
516 fHeader->SetPrimaryVertex(eventVertex);
518 gAlice->SetGenEventHeader(fHeader);
520 // AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
521 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));