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