Add possibility for heavy flavor selection: charm and beauty.
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.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 /*
17 $Log$
18 Revision 1.1  2000/06/09 20:47:27  morsch
19 AliGenerator interface class to HIJING using THijing (test version)
20
21 */
22
23 #include "AliGenHijing.h"
24 #include "AliRun.h"
25
26 #include <TArrayI.h>
27 #include <TParticle.h>
28 #include <THijing.h>
29
30  ClassImp(AliGenHijing)
31
32 AliGenHijing::AliGenHijing()
33                  :AliGenerator()
34 {
35 // Constructor
36 }
37
38 AliGenHijing::AliGenHijing(Int_t npart)
39                  :AliGenerator(npart)
40 {
41 // Default PbPb collisions at 5. 5 TeV
42 //
43     SetEnergyCMS();
44     SetImpactParameterRange();
45     SetTarget();
46     SetProjectile();
47     fKeep=0;
48     fQuench=1;
49     fShadowing=1;
50     fTrigger=0;
51     fDecaysOff=1;
52     fEvaluate=0;
53     fSelectAll=0;
54     fFlavor=0;
55 }
56
57 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
58 {
59 // copy constructor
60 }
61
62
63 AliGenHijing::~AliGenHijing()
64 {
65 // Destructor
66 }
67
68 void AliGenHijing::Init()
69 {
70 // Initialisation
71     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
72                       fAProjectile, fZProjectile, fATarget, fZTarget, 
73                       fMinImpactParam, fMaxImpactParam));
74
75     fHijing=(THijing*) fgMCEvGen;
76     fHijing->Initialize();
77     fHijing->SetIHPR2(3,  fTrigger);
78     fHijing->SetIHPR2(4,  fQuench);
79     fHijing->SetIHPR2(6,  fShadowing);
80     fHijing->SetIHPR2(12, fDecaysOff);    
81     fHijing->SetIHPR2(21, fKeep);
82 //
83     if (fEvaluate) EvaluateCrossSections();
84 }
85
86 void AliGenHijing::Generate()
87 {
88 // Generate one event
89
90     Float_t polar[3] =   {0,0,0};
91     Float_t origin[3]=   {0,0,0};
92     Float_t origin0[3]=  {0,0,0};
93     Float_t p[3], random[6];
94     Float_t tof;
95
96     static TClonesArray *particles;
97 //  converts from mm/c to s
98     const Float_t kconv=0.001/2.999792458e8;
99 //
100     Int_t nt=0;
101     Int_t jev=0;
102     Int_t j, kf, ks, imo;
103     kf=0;
104     
105     if(!particles) particles=new TClonesArray("TParticle",10000);
106     
107     fTrials=0;
108     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
109     if(fVertexSmear==perEvent) {
110         gMC->Rndm(random,6);
111         for (j=0;j<3;j++) {
112             origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
113                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
114 //          fHijing->SetMSTP(151,0);
115         }
116     } else if (fVertexSmear==perTrack) {
117 //      fHijing->SetMSTP(151,0);
118         for (j=0;j<3;j++) {
119 //          fHijing->SetPARP(151+j, fOsigma[j]*10.);
120         }
121     }
122     while(1)
123     {
124
125         fHijing->GenerateEvent();
126         fTrials++;
127         fHijing->ImportParticles(particles,"All");
128         Int_t np = particles->GetEntriesFast();
129         printf("\n **************************************************%d\n",np);
130         Int_t nc=0;
131         if (np == 0 ) continue;
132         Int_t i;
133         Int_t * newPos = new Int_t[np];
134         for (i = 0; i<np; i++) *(newPos+i)=i;
135         
136         for (i = 0; i<np; i++) {
137             TParticle *  iparticle       = (TParticle *) particles->At(i);
138
139             Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
140             Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
141             Bool_t  selected             =  kTRUE;
142             Bool_t  hasSelectedDaughters =  kFALSE;
143
144
145             kf        = iparticle->GetPdgCode();
146             if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
147             if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
148 //
149 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
150 //
151
152             if (selected || hasSelectedDaughters) {
153                 nc++;
154                 ks        = iparticle->GetStatusCode();
155                 p[0]=iparticle->Px();
156                 p[1]=iparticle->Py();
157                 p[2]=iparticle->Pz();
158                 origin[0]=origin0[0]+iparticle->Vx()/10;
159                 origin[1]=origin0[1]+iparticle->Vy()/10;
160                 origin[2]=origin0[2]+iparticle->Vz()/10;
161                 tof=kconv*iparticle->T();
162                 imo=-1;
163                 if (hasMother) {
164                     imo=iparticle->GetFirstMother();
165                     imo=*(newPos+imo);
166                 }
167                 
168 //              printf("\n selected iparent %d %d %d \n",i, kf, imo);
169                 if (hasDaughter) {
170                     gAlice->SetTrack(0,imo,kf,p,origin,polar,
171                                      tof,"Primary",nt);
172                 } else {
173                     gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
174                                      tof,"Secondary",nt);
175                 }
176                 *(newPos+i)=nt;
177             } // selected
178         } // particle loop 
179         delete newPos;
180         printf("\n I've put %i particles on the stack \n",nc);
181         if (nc > 0) {
182             jev+=nc;
183             if (jev >= fNpart || fNpart == -1) {
184                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
185                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
186                 break;
187             }
188         }
189     } // event loop
190 }
191
192 Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
193 {
194 // Perform kinematic selection
195     Float_t px=particle->Px();
196     Float_t py=particle->Py();
197     Float_t pz=particle->Pz();
198     Float_t  e=particle->Energy();
199
200 //
201 //  transverse momentum cut    
202     Float_t pt=TMath::Sqrt(px*px+py*py);
203     if (pt > fPtMax || pt < fPtMin) 
204     {
205 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
206         return kFALSE;
207     }
208 //
209 // momentum cut
210     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
211     if (p > fPMax || p < fPMin) 
212     {
213 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
214         return kFALSE;
215     }
216     
217 //
218 // theta cut
219     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
220     if (theta > fThetaMax || theta < fThetaMin) 
221     {
222         
223 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
224         return kFALSE;
225     }
226
227 //
228 // rapidity cut
229     Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
230     if (y > fYMax || y < fYMin)
231     {
232 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
233         return kFALSE;
234     }
235
236 //
237 // phi cut
238     Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
239     if (phi > fPhiMax || phi < fPhiMin)
240     {
241 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
242         return kFALSE;
243     }
244
245     return kTRUE;
246 }
247
248 void AliGenHijing::KeepFullEvent()
249 {
250     fKeep=1;
251 }
252
253 void AliGenHijing::EvaluateCrossSections()
254 {
255 //     Glauber Calculation of geometrical x-section
256 //
257     Float_t xTot=0.;          // barn
258     Float_t xTotHard=0.;      // barn 
259     Float_t xPart=0.;         // barn
260     Float_t xPartHard=0.;     // barn 
261     Float_t sigmaHard=0.1;    // mbarn
262     Float_t bMin=0.;
263     Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
264     const Float_t kdib=0.2;
265     Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;
266
267
268     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
269     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
270     Int_t i;
271     Float_t oldvalue=0.;
272     
273     for (i=0; i<kMax; i++)
274     {
275         Float_t xb=bMin+i*kdib;
276         Float_t ov;
277         ov=fHijing->Profile(xb);
278         Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
279         Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
280         xTot+=gb;
281         xTotHard+=gbh;
282         if (xb > fMinImpactParam && xb < fMaxImpactParam)
283         {
284             xPart+=gb;
285             xPartHard+=gbh;
286         }
287         
288         if ((xTot-oldvalue)/oldvalue<0.0001) break;
289         oldvalue=xTot;
290         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
291         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
292     }
293     printf("\n Total cross section (barn): %f \n",xTot);
294     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
295     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
296     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
297 }
298
299 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
300 {
301 //
302 // Looks recursively if one of the daughters has been selected
303 //
304 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
305     Int_t imin=-1;
306     Int_t imax=-1;
307     Int_t i;
308     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
309     Bool_t selected=kFALSE;
310     if (hasDaughters) {
311         imin=iparticle->GetFirstDaughter();
312         imax=iparticle->GetLastDaughter();       
313         for (i=imin; i<= imax; i++){
314             TParticle *  jparticle       = (TParticle *) particles->At(i);      
315             Int_t ip=jparticle->GetPdgCode();
316             if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
317             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
318         }
319     } else {
320         return kFALSE;
321     }
322
323     return selected;
324 }
325
326
327 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
328 {
329 // Select flavor of particle
330 // 0: all
331 // 4: charm and beauty
332 // 5: beauty
333     if (fFlavor == 0) return kTRUE;
334     
335     Int_t ifl=TMath::Abs(pid/100);
336     if (ifl > 10) ifl/=10;
337     return ((fFlavor==4 && (ifl==4 || ifl==5))  || 
338             (fFlavor==5 &&  ifl==5));
339
340 }
341
342 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
343 {
344 // Assignment operator
345     return *this;
346 }