5710647f5ab0f3d7853aba39983e0bbddffce768
[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 <TObjArray.h>
36 #include <TParticle.h>
37
38 #include "AliGenCocktailEntry.h"
39 #include "AliGenMUONCocktail.h"
40 #include "AliGenMUONlib.h"
41 #include "AliGenParam.h"
42 #include "AliMC.h"
43 #include "AliRun.h"
44 #include "AliStack.h"
45
46 ClassImp(AliGenMUONCocktail)
47   
48   
49 //________________________________________________________________________
50 AliGenMUONCocktail::AliGenMUONCocktail()
51   :AliGenCocktail()
52 {
53 // Constructor
54   fTotalRate =0;  
55   fNSucceded=0; 
56   fNGenerated=0; 
57   fMuonMultiplicity=1;
58   fMuonPtCut= 1.;
59   fMuonThetaMinCut=171.; 
60   fMuonThetaMaxCut=178.;
61   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
62     //
63 }
64 //_________________________________________________________________________
65 AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
66     AliGenCocktail(cocktail)
67 {
68 // Copy constructor
69   fTotalRate =0; 
70   fNSucceded=0; 
71   fNGenerated=0; 
72   fMuonMultiplicity=1;
73   fMuonPtCut= 1.;
74   fMuonThetaMinCut=171.; 
75   fMuonThetaMaxCut=178.;
76   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
77     //
78 }
79 //_________________________________________________________________________
80 AliGenMUONCocktail::~AliGenMUONCocktail()
81 {
82 // Destructor
83 }
84
85 //_________________________________________________________________________
86 void AliGenMUONCocktail::Init()
87 {
88   // Defining MUON physics cocktail
89
90   // Kinematical limits for particle generation
91   Float_t ptMin  = fPtMin;
92   Float_t ptMax  = fPtMax;
93   Float_t yMin   = fYMin;;
94   Float_t yMax   = fYMax;;
95   Float_t phiMin = fPhiMin*180./TMath::Pi();
96   Float_t phiMax = fPhiMax*180./TMath::Pi();
97   printf(">>> Kinematical range pT:%f:%f  y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
98
99   Float_t sigmaReaction =   0.072;   // MINB pp at LHC energies 72 mb
100
101   // Generating J/Psi Physics
102   AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
103   // 4pi Generation 
104   genjpsi->SetPtRange(0,100.);
105   genjpsi->SetYRange(-8.,8);
106   genjpsi->SetPhiRange(0.,360.);
107   genjpsi->SetForceDecay(kDiMuon);
108   genjpsi->SetTrackingFlag(1);
109   // Calculation of the paritcle multiplicity per event in the muonic channel
110   Float_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
111   Float_t sigmajpsi     = 31.0e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
112   Float_t brjpsi        = 0.0588;           // Branching Ratio for JPsi
113   genjpsi->Init();  // Generating pT and Y parametrsation for the 4pi
114   ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
115   printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratiojpsi,genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmajpsi );
116   // Generation in the kinematical limits of AliGenMUONCocktail
117   genjpsi->SetPtRange(ptMin, ptMax);  
118   genjpsi->SetYRange(yMin, yMax);
119   genjpsi->SetPhiRange(phiMin, phiMax);
120   genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
121   // Adding Generator 
122   AddGenerator(genjpsi, "Jpsi", ratiojpsi);
123   fTotalRate+=ratiojpsi;
124
125  // Generating Psi prime Physics
126   AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "Vogt", "PsiP");
127   // 4pi Generation 
128   genpsiP->SetPtRange(0,100.);
129   genpsiP->SetYRange(-8.,8);
130   genpsiP->SetPhiRange(0.,360.);
131   genpsiP->SetForceDecay(kDiMuon);
132   genpsiP->SetTrackingFlag(1);
133   // Calculation of the paritcle multiplicity per event in the muonic channel
134   Float_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
135   Float_t sigmapsiP     = 4.68e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
136   Float_t brpsiP        = 0.0103;           // Branching Ratio for PsiP
137   genpsiP->Init();  // Generating pT and Y parametrsation for the 4pi
138   ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
139   printf(">>> ratio psiP %g et %g Ncol %g sigma %g\n",ratiopsiP,genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmapsiP );
140   // Generation in the kinematical limits of AliGenMUONCocktail
141   genpsiP->SetPtRange(ptMin, ptMax);  
142   genpsiP->SetYRange(yMin, yMax);
143   genpsiP->SetPhiRange(phiMin, phiMax);
144   genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
145   // Adding Generator 
146   AddGenerator(genpsiP, "PsiP", ratiopsiP);
147   fTotalRate+=ratiopsiP;
148
149
150 // Generating Upsilon Physics
151   AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");  
152   genupsilon->SetPtRange(0,100.);  
153   genupsilon->SetYRange(-8.,8);
154   genupsilon->SetPhiRange(0.,360.);
155   genupsilon->SetForceDecay(kDiMuon);
156   genupsilon->SetTrackingFlag(1);
157   Float_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
158   Float_t sigmaupsilon     = 0.501e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
159   Float_t brupsilon        = 0.0248;  // Branching Ratio for Upsilon
160   genupsilon->Init();  // Generating pT and Y parametrsation for the 4pi
161   ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
162   printf(">>> ratio upsilon %g et %g\n",ratioupsilon, genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
163   genupsilon->SetPtRange(ptMin, ptMax);  
164   genupsilon->SetYRange(yMin, yMax);
165   genupsilon->SetPhiRange(phiMin, phiMax);
166   genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
167   AddGenerator(genupsilon,"Upsilon", ratioupsilon);
168   fTotalRate+=ratioupsilon;
169
170 // Generating UpsilonP Physics
171   AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "Vogt", "UpsilonP");  
172   genupsilonP->SetPtRange(0,100.);  
173   genupsilonP->SetYRange(-8.,8);
174   genupsilonP->SetPhiRange(0.,360.);
175   genupsilonP->SetForceDecay(kDiMuon);
176   genupsilonP->SetTrackingFlag(1);
177   Float_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
178   Float_t sigmaupsilonP     = 0.246e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
179   Float_t brupsilonP        = 0.0131;  // Branching Ratio for UpsilonP
180   genupsilonP->Init();  // Generating pT and Y parametrsation for the 4pi
181   ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
182   printf(">>> ratio upsilonP %g et %g\n",ratioupsilonP, genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
183   genupsilonP->SetPtRange(ptMin, ptMax);  
184   genupsilonP->SetYRange(yMin, yMax);
185   genupsilonP->SetPhiRange(phiMin, phiMax);
186   genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
187   AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
188   fTotalRate+=ratioupsilonP;
189
190
191 // Generating UpsilonPP Physics
192   AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "Vogt", "UpsilonPP");  
193   genupsilonPP->SetPtRange(0,100.);  
194   genupsilonPP->SetYRange(-8.,8);
195   genupsilonPP->SetPhiRange(0.,360.);
196   genupsilonPP->SetForceDecay(kDiMuon);
197   genupsilonPP->SetTrackingFlag(1);
198   Float_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
199   Float_t sigmaupsilonPP     = 0.100e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
200   Float_t brupsilonPP        = 0.0181;  // Branching Ratio for UpsilonPP
201   genupsilonPP->Init();  // Generating pT and Y parametrsation for the 4pi
202   ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
203   printf(">>> ratio upsilonPP %g et %g\n",ratioupsilonPP, genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
204   genupsilonPP->SetPtRange(ptMin, ptMax);  
205   genupsilonPP->SetYRange(yMin, yMax);
206   genupsilonPP->SetPhiRange(phiMin, phiMax);
207   genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
208   AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
209   fTotalRate+=ratioupsilonPP;
210
211
212 // Generating Charm Physics 
213   AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");  
214   gencharm->SetPtRange(0,100.);  
215   gencharm->SetYRange(-8.,8);
216   gencharm->SetPhiRange(0.,360.);
217   gencharm->SetForceDecay(kSemiMuonic);
218   gencharm->SetTrackingFlag(1);
219   Float_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
220   Float_t sigmacharm     = 2. * 6.64e-3 * 0.65 ;   // 
221   Float_t brcharm        = 0.12;  // Branching Ratio for Charm
222   gencharm->Init();  // Generating pT and Y parametrsation for the 4pi
223   ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * 
224     gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
225   gencharm->SetPtRange(ptMin, ptMax);  
226   gencharm->SetYRange(yMin, yMax);
227   gencharm->SetPhiRange(phiMin, phiMax);
228   gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
229
230   printf(">>> ratio charm %f\n",ratiocharm);
231   AddGenerator(gencharm,"Charm", ratiocharm);
232   fTotalRate+=ratiocharm;
233
234 // Generating Beauty Physics "Correlated Pairs"
235   AliGenParam * genbeauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");  
236   genbeauty->SetPtRange(0,100.);  
237   genbeauty->SetYRange(-8.,8);
238   genbeauty->SetPhiRange(0.,360.);
239   genbeauty->SetForceDecay(kSemiMuonic);
240   genbeauty->SetTrackingFlag(1);
241   Float_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
242   Float_t sigmabeauty     = 2. * 0.210e-3 * 0.84;   // 
243   Float_t brbeauty        = 0.15;  // Branching Ratio for Beauty
244   genbeauty->Init();  // Generating pT and Y parametrsation for the 4pi
245   ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction * 
246     genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
247   genbeauty->SetPtRange(ptMin, ptMax);  
248   genbeauty->SetYRange(yMin, yMax);
249   genbeauty->SetPhiRange(phiMin, phiMax);
250   genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
251
252   printf(">>> ratio beauty %f\n",ratiobeauty);
253   AddGenerator(genbeauty,"Beauty", ratiobeauty);
254   fTotalRate+=ratiobeauty;
255
256 // Generating Pion Physics
257   AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");  
258   genpion->SetPtRange(0,100.);  
259   genpion->SetYRange(-8.,8);
260   genpion->SetPhiRange(0.,360.);
261   genpion->SetForceDecay(kPiToMu);
262   genpion->SetTrackingFlag(1);
263   Float_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
264   Float_t sigmapion     = 0.93e-2; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
265   Float_t brpion        = 0.9999;  // Branching Ratio for Pion 
266   genpion->Init();  // Generating pT and Y parametrsation for the 4pi
267   ratiopion = sigmapion * brpion *  (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
268   genpion->SetPtRange(ptMin, ptMax);  
269   genpion->SetYRange(yMin, yMax);
270   genpion->SetPhiRange(phiMin, phiMax);
271   genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
272   printf(">>> ratio pion %f\n",ratiopion);
273   AddGenerator(genpion,"Pion", ratiopion);
274   fTotalRate+=ratiopion;
275
276 // Generating Kaon Physics
277   AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");  
278   genkaon->SetPtRange(0,100.);  
279   genkaon->SetYRange(-8.,8);
280   genkaon->SetPhiRange(0.,360.);
281   genkaon->SetForceDecay(kKaToMu);
282   genkaon->SetTrackingFlag(1);
283   Float_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
284   Float_t sigmakaon     = 1.23e-4;   // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06 
285   Float_t brkaon        = 0.6351 ;  // Branching Ratio for Kaon 
286   genkaon->Init();  // Generating pT and Y parametrsation for the 4pi
287   ratiokaon = sigmakaon * brkaon * (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
288   genkaon->SetPtRange(ptMin, ptMax);  
289   genkaon->SetYRange(yMin, yMax);
290   genkaon->SetPhiRange(phiMin, phiMax);
291   genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
292   printf(">>> ratio kaon %f\n",ratiokaon);
293   AddGenerator(genkaon,"Kaon", ratiokaon);
294   fTotalRate+=ratiokaon;
295 }
296
297 //_________________________________________________________________________
298 void AliGenMUONCocktail::Generate()
299 {
300   //
301 // Generate event 
302     TIter next(fEntries);
303     AliGenCocktailEntry *entry = 0;
304     AliGenCocktailEntry *preventry = 0;
305     AliGenerator* gen = 0;
306
307     TObjArray *partArray = gAlice->GetMCApp()->Particles();
308
309 //
310 //  Generate the vertex position used by all generators
311 //    
312     if(fVertexSmear == kPerEvent) Vertex();
313     Bool_t primordialTrigger = kFALSE;
314
315     while(!primordialTrigger) {
316
317       //Reseting stack
318       AliRunLoader * runloader = gAlice->GetRunLoader();
319       if (runloader)
320         if (runloader->Stack())
321           runloader->Stack()->Reset();
322       //
323       // Loop over generators and generate events
324       Int_t igen=0;
325       Int_t npart =0;
326       
327       while((entry = (AliGenCocktailEntry*)next())) {
328         gen = entry->Generator();
329         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
330         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
331           igen++;       
332           if (igen ==1) entry->SetFirst(0);
333           else  entry->SetFirst((partArray->GetEntriesFast())+1);
334           gen->SetNumberParticles(npart);
335           gen->Generate();
336           entry->SetLast(partArray->GetEntriesFast());
337           preventry = entry;
338         }
339       }  
340       next.Reset();
341
342       // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance 171,178 and pTCut
343       Int_t iPart;
344       fNGenerated++;
345       Int_t numberOfMuons=0;
346       //      printf(">>>fNGenerated is %d\n",fNGenerated);
347       for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){      
348         //      gAlice->GetMCApp()->Particle(iPart)->Print();
349         if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
350              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
351              (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
352              (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) { 
353           gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
354           numberOfMuons++;
355         }
356       }
357       //  printf(">>> Number of Muons is %d \n", numberOfMuons);
358       if (numberOfMuons >= fMuonMultiplicity ) primordialTrigger = kTRUE;
359     }
360     //printf(">>> Trigger Accepted!!! \n");
361     fNSucceded++;
362     //   Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
363     //    printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
364 }
365
366
367
368
369
370