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