First commit.
[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 // Classe to create the MUON coktail for physics in the Alice muon spectrometer
19 // Gines Martinez, jan 2004, Nantes martinez@in2p3.fr
20
21
22 //
23
24 #include <TList.h>
25 #include <TObjArray.h>
26 #include <TF1.h>
27 #include <TParticle.h>
28
29 #include "AliGenParam.h"
30 #include "AliGenMUONlib.h"
31 #include "AliGenMUONCocktail.h"
32 #include "AliGenCocktailEntry.h"
33 #include "AliCollisionGeometry.h"
34 #include "AliRun.h"
35 #include "AliMC.h"
36 #include "AliStack.h"
37
38 ClassImp(AliGenMUONCocktail)
39   
40   
41 //________________________________________________________________________
42 AliGenMUONCocktail::AliGenMUONCocktail()
43   :AliGenCocktail()
44 {
45 // Constructor
46   fTotalRate =0;  
47   fNSucceded=0; 
48   fNGenerated=0; 
49   fMuonMultiplicity=1;
50   fMuonPtCut= 1.;
51   fMuonThetaMinCut=171.; 
52   fMuonThetaMaxCut=178.;
53   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
54     //
55 }
56 //_________________________________________________________________________
57 AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
58     AliGenCocktail(cocktail)
59 {
60 // Copy constructor
61   fTotalRate =0; 
62   fNSucceded=0; 
63   fNGenerated=0; 
64   fMuonMultiplicity=1;
65   fMuonPtCut= 1.;
66   fMuonThetaMinCut=171.; 
67   fMuonThetaMaxCut=178.;
68   fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
69     //
70 }
71 //_________________________________________________________________________
72 AliGenMUONCocktail::~AliGenMUONCocktail()
73 {
74 // Destructor
75     delete fEntries;
76 }
77
78 //_________________________________________________________________________
79 void AliGenMUONCocktail::Init()
80 {
81   // Defining MUON physics cocktail
82
83   // Kinematical limits for particle generation
84   Float_t ptMin  = fPtMin;
85   Float_t ptMax  = fPtMax;
86   Float_t yMin   = fYMin;;
87   Float_t yMax   = fYMax;;
88   Float_t phiMin = fPhiMin*180./TMath::Pi();
89   Float_t phiMax = fPhiMax*180./TMath::Pi();
90   printf(">>> Kinematical range pT:%f:%f  y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
91
92   Float_t sigma_reaction =   0.072;   // MINB pp at LHC energies 72 mb
93   
94   // Generating J/Psi Physics
95   AliGenParam * gen_jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
96   // 4pi Generation 
97   gen_jpsi->SetPtRange(0,100.);
98   gen_jpsi->SetYRange(-8.,8);
99   gen_jpsi->SetPhiRange(0.,360.);
100   gen_jpsi->SetForceDecay(kDiMuon);
101   gen_jpsi->SetTrackingFlag(1);
102   // Calculation of the paritcle multiplicity per event in the muonic channel
103   Float_t ratio_jpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
104   Float_t sigma_jpsi     = 31.0e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
105   Float_t br_jpsi        = 0.0588;           // Branching Ratio for JPsi
106   gen_jpsi->Init();  // Generating pT and Y parametrsation for the 4pi
107   ratio_jpsi = sigma_jpsi * br_jpsi * fNumberOfCollisions / sigma_reaction * gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
108   printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratio_jpsi,gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigma_jpsi );
109   // Generation in the kinematical limits of AliGenMUONCocktail
110   gen_jpsi->SetPtRange(ptMin, ptMax);  
111   gen_jpsi->SetYRange(yMin, yMax);
112   gen_jpsi->SetPhiRange(phiMin, phiMax);
113   gen_jpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
114   // Adding Generator 
115   AddGenerator(gen_jpsi, "Jpsi", ratio_jpsi);
116   fTotalRate+=ratio_jpsi;
117
118 // Generating Upsilon Physics
119   AliGenParam * gen_upsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");  
120   gen_upsilon->SetPtRange(0,100.);  
121   gen_upsilon->SetYRange(-8.,8);
122   gen_upsilon->SetPhiRange(0.,360.);
123   gen_upsilon->SetForceDecay(kDiMuon);
124   gen_upsilon->SetTrackingFlag(1);
125   Float_t ratio_upsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
126   Float_t sigma_upsilon     = 0.501e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
127   Float_t br_upsilon        = 0.0248;  // Branching Ratio for Upsilon
128   gen_upsilon->Init();  // Generating pT and Y parametrsation for the 4pi
129   ratio_upsilon = sigma_upsilon * br_upsilon * fNumberOfCollisions / sigma_reaction * gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
130   printf(">>> ratio upsilon %g et %g\n",ratio_upsilon, gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
131   gen_upsilon->SetPtRange(ptMin, ptMax);  
132   gen_upsilon->SetYRange(yMin, yMax);
133   gen_upsilon->SetPhiRange(phiMin, phiMax);
134   gen_upsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
135   AddGenerator(gen_upsilon,"Upsilon", ratio_upsilon);
136   fTotalRate+=ratio_upsilon;
137
138 // Generating Charm Physics 
139   AliGenParam * gen_charm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");  
140   gen_charm->SetPtRange(0,100.);  
141   gen_charm->SetYRange(-8.,8);
142   gen_charm->SetPhiRange(0.,360.);
143   gen_charm->SetForceDecay(kSemiMuonic);
144   gen_charm->SetTrackingFlag(1);
145   Float_t ratio_charm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
146   Float_t sigma_charm     = 2. * 6.64e-3 * 0.65 ;   // 
147   Float_t br_charm        = 0.12;  // Branching Ratio for Charm
148   gen_charm->Init();  // Generating pT and Y parametrsation for the 4pi
149   ratio_charm = sigma_charm * br_charm * fNumberOfCollisions / sigma_reaction * 
150     gen_charm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
151   gen_charm->SetPtRange(ptMin, ptMax);  
152   gen_charm->SetYRange(yMin, yMax);
153   gen_charm->SetPhiRange(phiMin, phiMax);
154   gen_charm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
155
156   printf(">>> ratio charm %f\n",ratio_charm);
157   AddGenerator(gen_charm,"Charm", ratio_charm);
158   fTotalRate+=ratio_charm;
159
160 // Generating Beauty Physics "Correlated Pairs"
161   AliGenParam * gen_beauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");  
162   gen_beauty->SetPtRange(0,100.);  
163   gen_beauty->SetYRange(-8.,8);
164   gen_beauty->SetPhiRange(0.,360.);
165   gen_beauty->SetForceDecay(kSemiMuonic);
166   gen_beauty->SetTrackingFlag(1);
167   Float_t ratio_beauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
168   Float_t sigma_beauty     = 2. * 0.210e-3 * 0.84;   // 
169   Float_t br_beauty        = 0.15;  // Branching Ratio for Beauty
170   gen_beauty->Init();  // Generating pT and Y parametrsation for the 4pi
171   ratio_beauty = sigma_beauty * br_beauty * fNumberOfCollisions / sigma_reaction * 
172     gen_beauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
173   gen_beauty->SetPtRange(ptMin, ptMax);  
174   gen_beauty->SetYRange(yMin, yMax);
175   gen_beauty->SetPhiRange(phiMin, phiMax);
176   gen_beauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
177
178   printf(">>> ratio beauty %f\n",ratio_beauty);
179   AddGenerator(gen_beauty,"Beauty", ratio_beauty);
180   fTotalRate+=ratio_beauty;
181
182 // Generating Pion Physics
183   AliGenParam * gen_pion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");  
184   gen_pion->SetPtRange(0,100.);  
185   gen_pion->SetYRange(-8.,8);
186   gen_pion->SetPhiRange(0.,360.);
187   gen_pion->SetForceDecay(kPiToMu);
188   gen_pion->SetTrackingFlag(1);
189   Float_t ratio_pion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
190   Float_t sigma_pion     = 1.6e-2;   // A ojo TO be studied in detail.  
191   Float_t br_pion        = 0.9999;  // Branching Ratio for Pion 
192   gen_pion->Init();  // Generating pT and Y parametrsation for the 4pi
193   ratio_pion = sigma_pion * br_pion * fNumberOfParticipants / sigma_reaction * gen_pion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
194   gen_pion->SetPtRange(ptMin, ptMax);  
195   gen_pion->SetYRange(yMin, yMax);
196   gen_pion->SetPhiRange(phiMin, phiMax);
197   gen_pion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
198   printf(">>> ratio pion %f\n",ratio_pion);
199   AddGenerator(gen_pion,"Pion", ratio_pion);
200   fTotalRate+=ratio_pion;
201
202 // Generating Kaon Physics
203   AliGenParam * gen_kaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");  
204   gen_kaon->SetPtRange(0,100.);  
205   gen_kaon->SetYRange(-8.,8);
206   gen_kaon->SetPhiRange(0.,360.);
207   gen_kaon->SetForceDecay(kKaToMu);
208   gen_kaon->SetTrackingFlag(1);
209   Float_t ratio_kaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
210   Float_t sigma_kaon     = 1.8e-4;   // OJO 
211   Float_t br_kaon        = 0.6351 ;  // Branching Ratio for Kaon 
212   gen_kaon->Init();  // Generating pT and Y parametrsation for the 4pi
213   ratio_kaon = sigma_kaon * br_kaon * fNumberOfParticipants/ sigma_reaction * gen_kaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
214   gen_kaon->SetPtRange(ptMin, ptMax);  
215   gen_kaon->SetYRange(yMin, yMax);
216   gen_kaon->SetPhiRange(phiMin, phiMax);
217   gen_kaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
218   printf(">>> ratio kaon %f\n",ratio_kaon);
219   AddGenerator(gen_kaon,"Kaon", ratio_kaon);
220   fTotalRate+=ratio_kaon;
221 }
222
223 //_________________________________________________________________________
224 void AliGenMUONCocktail::Generate()
225 {
226   //
227 // Generate event 
228     TIter next(fEntries);
229     AliGenCocktailEntry *entry = 0;
230     AliGenCocktailEntry *preventry = 0;
231     AliGenerator* gen = 0;
232
233     TObjArray *partArray = gAlice->GetMCApp()->Particles();
234
235 //
236 //  Generate the vertex position used by all generators
237 //    
238     if(fVertexSmear == kPerEvent) Vertex();
239     Bool_t PrimordialTrigger = kFALSE;
240
241     while(!PrimordialTrigger) {
242
243       //Reseting stack
244       AliRunLoader * runloader = gAlice->GetRunLoader();
245       if (runloader)
246         if (runloader->Stack())
247           runloader->Stack()->Reset();
248       //
249       // Loop over generators and generate events
250       Int_t igen=0;
251       Int_t npart =0;
252       
253       while((entry = (AliGenCocktailEntry*)next())) {
254         gen = entry->Generator();
255         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
256         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
257           igen++;       
258           if (igen ==1) entry->SetFirst(0);
259           else  entry->SetFirst((partArray->GetEntriesFast())+1);
260           gen->SetNumberParticles(npart);
261           gen->Generate();
262           entry->SetLast(partArray->GetEntriesFast());
263           preventry = entry;
264         }
265       }  
266       next.Reset();
267
268       // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance 171,178 and pTCut
269       Int_t iPart;
270       fNGenerated++;
271       Int_t numberOfMuons=0;
272       //      printf(">>>fNGenerated is %d\n",fNGenerated);
273       for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){      
274         //      gAlice->GetMCApp()->Particle(iPart)->Print();
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           gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
280           numberOfMuons++;
281         }
282       }
283       //  printf(">>> Number of Muons is %d \n", numberOfMuons);
284       if (numberOfMuons >= fMuonMultiplicity ) PrimordialTrigger = kTRUE;
285     }
286     //printf(">>> Trigger Accepted!!! \n");
287     fNSucceded++;
288     //   Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
289     //    printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
290 }
291
292
293
294
295
296