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