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