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