Coding rule violations corrected (G. Martinez)
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.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 MUON coktail for physics in the Alice muon spectrometer
20 // The followoing muons sources are included in this cocktail: 
21 //     jpsi, upsilon, non-correlated open and beauty, and muons from pion and kaons.
22 // The free parameeters are :
23 //      pp reaction cross-section
24 //      production cross-sections in pp collisions and 
25 //      branching ratios in the muon channel
26 // Hard probes are supposed to scale with Ncoll and hadronic production with (0.8Ncoll+0.2*Npart)
27 // There is a primordial trigger wiche requires :
28 //      a minimum muon multiplicity above a pT cut in a theta acceptance cone
29 //
30 // Gines Martinez, jan 2004, Nantes  martinez@in2p3.fr
31
32
33 //
34
35 //#include <TList.h>
36 #include <TObjArray.h>
37 #include <TF1.h>
38 #include <TParticle.h>
39
40 #include "AliGenParam.h"
41 #include "AliGenMUONlib.h"
42 #include "AliGenMUONCocktail.h"
43 #include "AliGenCocktailEntry.h"
44 #include "AliCollisionGeometry.h"
45 #include "AliRun.h"
46 #include "AliMC.h"
47 #include "AliStack.h"
48
49 ClassImp(AliGenMUONCocktail)
50   
51   
52 //________________________________________________________________________
53 AliGenMUONCocktail::AliGenMUONCocktail()
54   :AliGenCocktail()
55 {
56 // Constructor
57   fTotalRate =0;  
58   fNSucceded=0; 
59   fNGenerated=0; 
60   fMuonMultiplicity=1;
61   fMuonPtCut= 1.;
62   fMuonThetaMinCut=171.; 
63   fMuonThetaMaxCut=178.;
64   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
65     //
66 }
67 //_________________________________________________________________________
68 AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
69     AliGenCocktail(cocktail)
70 {
71 // Copy constructor
72   fTotalRate =0; 
73   fNSucceded=0; 
74   fNGenerated=0; 
75   fMuonMultiplicity=1;
76   fMuonPtCut= 1.;
77   fMuonThetaMinCut=171.; 
78   fMuonThetaMaxCut=178.;
79   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
80     //
81 }
82 //_________________________________________________________________________
83 AliGenMUONCocktail::~AliGenMUONCocktail()
84 {
85 // Destructor
86     delete fEntries;
87 }
88
89 //_________________________________________________________________________
90 void AliGenMUONCocktail::Init()
91 {
92   // Defining MUON physics cocktail
93
94   // Kinematical limits for particle generation
95   Float_t ptMin  = fPtMin;
96   Float_t ptMax  = fPtMax;
97   Float_t yMin   = fYMin;;
98   Float_t yMax   = fYMax;;
99   Float_t phiMin = fPhiMin*180./TMath::Pi();
100   Float_t phiMax = fPhiMax*180./TMath::Pi();
101   printf(">>> Kinematical range pT:%f:%f  y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
102
103   Float_t sigmaReaction =   0.072;   // MINB pp at LHC energies 72 mb
104   
105   // Generating J/Psi Physics
106   AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
107   // 4pi Generation 
108   genjpsi->SetPtRange(0,100.);
109   genjpsi->SetYRange(-8.,8);
110   genjpsi->SetPhiRange(0.,360.);
111   genjpsi->SetForceDecay(kDiMuon);
112   genjpsi->SetTrackingFlag(1);
113   // Calculation of the paritcle multiplicity per event in the muonic channel
114   Float_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
115   Float_t sigmajpsi     = 31.0e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
116   Float_t brjpsi        = 0.0588;           // Branching Ratio for JPsi
117   genjpsi->Init();  // Generating pT and Y parametrsation for the 4pi
118   ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
119   printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratiojpsi,genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmajpsi );
120   // Generation in the kinematical limits of AliGenMUONCocktail
121   genjpsi->SetPtRange(ptMin, ptMax);  
122   genjpsi->SetYRange(yMin, yMax);
123   genjpsi->SetPhiRange(phiMin, phiMax);
124   genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
125   // Adding Generator 
126   AddGenerator(genjpsi, "Jpsi", ratiojpsi);
127   fTotalRate+=ratiojpsi;
128
129 // Generating Upsilon Physics
130   AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");  
131   genupsilon->SetPtRange(0,100.);  
132   genupsilon->SetYRange(-8.,8);
133   genupsilon->SetPhiRange(0.,360.);
134   genupsilon->SetForceDecay(kDiMuon);
135   genupsilon->SetTrackingFlag(1);
136   Float_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
137   Float_t sigmaupsilon     = 0.501e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
138   Float_t brupsilon        = 0.0248;  // Branching Ratio for Upsilon
139   genupsilon->Init();  // Generating pT and Y parametrsation for the 4pi
140   ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
141   printf(">>> ratio upsilon %g et %g\n",ratioupsilon, genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
142   genupsilon->SetPtRange(ptMin, ptMax);  
143   genupsilon->SetYRange(yMin, yMax);
144   genupsilon->SetPhiRange(phiMin, phiMax);
145   genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
146   AddGenerator(genupsilon,"Upsilon", ratioupsilon);
147   fTotalRate+=ratioupsilon;
148
149 // Generating Charm Physics 
150   AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");  
151   gencharm->SetPtRange(0,100.);  
152   gencharm->SetYRange(-8.,8);
153   gencharm->SetPhiRange(0.,360.);
154   gencharm->SetForceDecay(kSemiMuonic);
155   gencharm->SetTrackingFlag(1);
156   Float_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
157   Float_t sigmacharm     = 2. * 6.64e-3 * 0.65 ;   // 
158   Float_t brcharm        = 0.12;  // Branching Ratio for Charm
159   gencharm->Init();  // Generating pT and Y parametrsation for the 4pi
160   ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * 
161     gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
162   gencharm->SetPtRange(ptMin, ptMax);  
163   gencharm->SetYRange(yMin, yMax);
164   gencharm->SetPhiRange(phiMin, phiMax);
165   gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
166
167   printf(">>> ratio charm %f\n",ratiocharm);
168   AddGenerator(gencharm,"Charm", ratiocharm);
169   fTotalRate+=ratiocharm;
170
171 // Generating Beauty Physics "Correlated Pairs"
172   AliGenParam * genbeauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");  
173   genbeauty->SetPtRange(0,100.);  
174   genbeauty->SetYRange(-8.,8);
175   genbeauty->SetPhiRange(0.,360.);
176   genbeauty->SetForceDecay(kSemiMuonic);
177   genbeauty->SetTrackingFlag(1);
178   Float_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
179   Float_t sigmabeauty     = 2. * 0.210e-3 * 0.84;   // 
180   Float_t brbeauty        = 0.15;  // Branching Ratio for Beauty
181   genbeauty->Init();  // Generating pT and Y parametrsation for the 4pi
182   ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction * 
183     genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
184   genbeauty->SetPtRange(ptMin, ptMax);  
185   genbeauty->SetYRange(yMin, yMax);
186   genbeauty->SetPhiRange(phiMin, phiMax);
187   genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
188
189   printf(">>> ratio beauty %f\n",ratiobeauty);
190   AddGenerator(genbeauty,"Beauty", ratiobeauty);
191   fTotalRate+=ratiobeauty;
192
193 // Generating Pion Physics
194   AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");  
195   genpion->SetPtRange(0,100.);  
196   genpion->SetYRange(-8.,8);
197   genpion->SetPhiRange(0.,360.);
198   genpion->SetForceDecay(kPiToMu);
199   genpion->SetTrackingFlag(1);
200   Float_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
201   Float_t sigmapion     = 0.93e-2; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
202   Float_t brpion        = 0.9999;  // Branching Ratio for Pion 
203   genpion->Init();  // Generating pT and Y parametrsation for the 4pi
204   ratiopion = sigmapion * brpion *  (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
205   genpion->SetPtRange(ptMin, ptMax);  
206   genpion->SetYRange(yMin, yMax);
207   genpion->SetPhiRange(phiMin, phiMax);
208   genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
209   printf(">>> ratio pion %f\n",ratiopion);
210   AddGenerator(genpion,"Pion", ratiopion);
211   fTotalRate+=ratiopion;
212
213 // Generating Kaon Physics
214   AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");  
215   genkaon->SetPtRange(0,100.);  
216   genkaon->SetYRange(-8.,8);
217   genkaon->SetPhiRange(0.,360.);
218   genkaon->SetForceDecay(kKaToMu);
219   genkaon->SetTrackingFlag(1);
220   Float_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
221   Float_t sigmakaon     = 1.23e-4;   // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06 
222   Float_t brkaon        = 0.6351 ;  // Branching Ratio for Kaon 
223   genkaon->Init();  // Generating pT and Y parametrsation for the 4pi
224   ratiokaon = sigmakaon * brkaon * (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
225   genkaon->SetPtRange(ptMin, ptMax);  
226   genkaon->SetYRange(yMin, yMax);
227   genkaon->SetPhiRange(phiMin, phiMax);
228   genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
229   printf(">>> ratio kaon %f\n",ratiokaon);
230   AddGenerator(genkaon,"Kaon", ratiokaon);
231   fTotalRate+=ratiokaon;
232 }
233
234 //_________________________________________________________________________
235 void AliGenMUONCocktail::Generate()
236 {
237   //
238 // Generate event 
239     TIter next(fEntries);
240     AliGenCocktailEntry *entry = 0;
241     AliGenCocktailEntry *preventry = 0;
242     AliGenerator* gen = 0;
243
244     TObjArray *partArray = gAlice->GetMCApp()->Particles();
245
246 //
247 //  Generate the vertex position used by all generators
248 //    
249     if(fVertexSmear == kPerEvent) Vertex();
250     Bool_t primordialTrigger = kFALSE;
251
252     while(!primordialTrigger) {
253
254       //Reseting stack
255       AliRunLoader * runloader = gAlice->GetRunLoader();
256       if (runloader)
257         if (runloader->Stack())
258           runloader->Stack()->Reset();
259       //
260       // Loop over generators and generate events
261       Int_t igen=0;
262       Int_t npart =0;
263       
264       while((entry = (AliGenCocktailEntry*)next())) {
265         gen = entry->Generator();
266         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
267         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
268           igen++;       
269           if (igen ==1) entry->SetFirst(0);
270           else  entry->SetFirst((partArray->GetEntriesFast())+1);
271           gen->SetNumberParticles(npart);
272           gen->Generate();
273           entry->SetLast(partArray->GetEntriesFast());
274           preventry = entry;
275         }
276       }  
277       next.Reset();
278
279       // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance 171,178 and pTCut
280       Int_t iPart;
281       fNGenerated++;
282       Int_t numberOfMuons=0;
283       //      printf(">>>fNGenerated is %d\n",fNGenerated);
284       for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){      
285         //      gAlice->GetMCApp()->Particle(iPart)->Print();
286         if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
287              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
288              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
289              (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) { 
290           gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
291           numberOfMuons++;
292         }
293       }
294       //  printf(">>> Number of Muons is %d \n", numberOfMuons);
295       if (numberOfMuons >= fMuonMultiplicity ) primordialTrigger = kTRUE;
296     }
297     //printf(">>> Trigger Accepted!!! \n");
298     fNSucceded++;
299     //   Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
300     //    printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
301 }
302
303
304
305
306
307