65543c2b6212a6cf35189faf2f930f2df5212db6
[u/mrichter/AliRoot.git] / 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