Possibility to select/unselect spectator protons and neutrons.
[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.12  2000/10/20 13:38:38  morsch
19 Debug printouts commented.
20
21 Revision 1.11  2000/10/20 13:22:26  morsch
22 - skip particle type 92 (string)
23 - Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
24   mesons.
25
26 Revision 1.10  2000/10/17 15:10:20  morsch
27 Write first all the parent particles to the stack and then the final state particles.
28
29 Revision 1.9  2000/10/17 13:38:59  morsch
30 Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
31
32 Revision 1.8  2000/10/17 12:46:31  morsch
33 Protect EvaluateCrossSections() against division by zero.
34
35 Revision 1.7  2000/10/02 21:28:06  fca
36 Removal of useless dependecies via forward declarations
37
38 Revision 1.6  2000/09/11 13:23:37  morsch
39 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
40 luget_hijing and luset_hijing.
41
42 Revision 1.5  2000/09/07 16:55:40  morsch
43 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
44
45 Revision 1.4  2000/07/11 18:24:56  fca
46 Coding convention corrections + few minor bug fixes
47
48 Revision 1.3  2000/06/30 12:08:36  morsch
49 In member data: char* replaced by TString, Init takes care of resizing the strings to
50 8 characters required by Hijing.
51
52 Revision 1.2  2000/06/15 14:15:05  morsch
53 Add possibility for heavy flavor selection: charm and beauty.
54
55 Revision 1.1  2000/06/09 20:47:27  morsch
56 AliGenerator interface class to HIJING using THijing (test version)
57
58 */
59
60 #include "AliGenHijing.h"
61 #include "AliGenHijingEventHeader.h"
62 #include "AliRun.h"
63 #include "AliMC.h"
64
65 #include <TArrayI.h>
66 #include <TParticle.h>
67 #include <THijing.h>
68
69  ClassImp(AliGenHijing)
70
71 AliGenHijing::AliGenHijing()
72                  :AliGenerator()
73 {
74 // Constructor
75 }
76
77 AliGenHijing::AliGenHijing(Int_t npart)
78                  :AliGenerator(npart)
79 {
80 // Default PbPb collisions at 5. 5 TeV
81 //
82     SetEnergyCMS();
83     SetImpactParameterRange();
84     SetTarget();
85     SetProjectile();
86     fKeep=0;
87     fQuench=1;
88     fShadowing=1;
89     fTrigger=0;
90     fDecaysOff=1;
91     fEvaluate=0;
92     fSelectAll=0;
93     fFlavor=0;
94     fSpectators=1;
95 }
96
97 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
98 {
99 // copy constructor
100 }
101
102
103 AliGenHijing::~AliGenHijing()
104 {
105 // Destructor
106 }
107
108 void AliGenHijing::Init()
109 {
110 // Initialisation
111     fFrame.Resize(8);
112     fTarget.Resize(8);
113     fProjectile.Resize(8);
114     
115     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
116                       fAProjectile, fZProjectile, fATarget, fZTarget, 
117                       fMinImpactParam, fMaxImpactParam));
118
119     fHijing=(THijing*) fgMCEvGen;
120
121     fHijing->SetIHPR2(3,  fTrigger);
122     fHijing->SetIHPR2(4,  fQuench);
123     fHijing->SetIHPR2(6,  fShadowing);
124     fHijing->SetIHPR2(12, fDecaysOff);    
125     fHijing->SetIHPR2(21, fKeep);
126     fHijing->Rluset(50,0);
127     fHijing->Initialize();
128
129     
130 //
131     if (fEvaluate) EvaluateCrossSections();
132 }
133
134 void AliGenHijing::Generate()
135 {
136 // Generate one event
137
138     Float_t polar[3] =   {0,0,0};
139     Float_t origin[3]=   {0,0,0};
140     Float_t origin0[3]=  {0,0,0};
141     Float_t p[3], random[6];
142     Float_t tof;
143
144     static TClonesArray *particles;
145 //  converts from mm/c to s
146     const Float_t kconv=0.001/2.999792458e8;
147 //
148     Int_t nt=0;
149     Int_t jev=0;
150     Int_t j, kf, ks, imo;
151     kf=0;
152     
153     if(!particles) particles=new TClonesArray("TParticle",10000);
154     
155     fTrials=0;
156     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
157     if(fVertexSmear==kPerEvent) {
158         gMC->Rndm(random,6);
159         for (j=0;j<3;j++) {
160             origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
161                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
162 //          fHijing->SetMSTP(151,0);
163         }
164     } else if (fVertexSmear==kPerTrack) {
165 //      fHijing->SetMSTP(151,0);
166         for (j=0;j<3;j++) {
167 //          fHijing->SetPARP(151+j, fOsigma[j]*10.);
168         }
169     }
170     while(1)
171     {
172
173         fHijing->GenerateEvent();
174         fTrials++;
175         fHijing->ImportParticles(particles,"All");
176         Int_t np = particles->GetEntriesFast();
177         printf("\n **************************************************%d\n",np);
178         Int_t nc=0;
179         if (np == 0 ) continue;
180         Int_t i;
181         Int_t * newPos = new Int_t[np];
182
183         for (i = 0; i<np; i++) *(newPos+i)=i;
184 //
185 //      First write parent particles
186 //
187
188         for (i = 0; i<np; i++) {
189             TParticle *  iparticle       = (TParticle *) particles->At(i);
190 // Is this a parent particle ?
191             if (Stable(iparticle)) continue;
192 //
193             Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
194             Bool_t  selected             =  kTRUE;
195             Bool_t  hasSelectedDaughters =  kFALSE;
196             
197             
198             kf        = iparticle->GetPdgCode();
199             ks        = iparticle->GetStatusCode();
200             if (kf == 92) continue;
201             
202             if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
203             hasSelectedDaughters = DaughtersSelection(iparticle, particles);
204 //
205 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
206 //
207             if (selected || hasSelectedDaughters) {
208                 nc++;
209                 p[0]=iparticle->Px();
210                 p[1]=iparticle->Py();
211                 p[2]=iparticle->Pz();
212                 origin[0]=origin0[0]+iparticle->Vx()/10;
213                 origin[1]=origin0[1]+iparticle->Vy()/10;
214                 origin[2]=origin0[2]+iparticle->Vz()/10;
215                 tof=kconv*iparticle->T();
216                 imo=-1;
217                 if (hasMother) {
218                     imo=iparticle->GetFirstMother();
219                     TParticle* mother= (TParticle *) particles->At(imo);
220                     imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
221                 }
222 // Put particle on the stack ... 
223 //              printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
224
225                 gAlice->SetTrack(0,imo,kf,p,origin,polar,
226                                  tof,"Primary",nt);
227 // ... and keep it there
228                 gAlice->KeepTrack(nt);
229 //
230                 *(newPos+i)=nt;
231             } // selected
232         } // particle loop parents
233 //
234 // Now write the final state particles
235 //
236
237         for (i = 0; i<np; i++) {
238             TParticle *  iparticle       = (TParticle *) particles->At(i);
239 // Is this a final state particle ?
240             if (!Stable(iparticle)) continue;
241 //          
242             Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
243             Bool_t  selected             =  kTRUE;
244             kf        = iparticle->GetPdgCode();
245             ks        = iparticle->GetStatusCode();
246             if (!fSelectAll) {
247               selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
248               if (!fSpectators && selected) selected = (ks != 0 && ks != 10);
249             }
250 //
251 // Put particle on the stack if selected
252 //
253             if (selected) {
254                 nc++;
255                 p[0]=iparticle->Px();
256                 p[1]=iparticle->Py();
257                 p[2]=iparticle->Pz();
258                 origin[0]=origin0[0]+iparticle->Vx()/10;
259                 origin[1]=origin0[1]+iparticle->Vy()/10;
260                 origin[2]=origin0[2]+iparticle->Vz()/10;
261                 tof=kconv*iparticle->T();
262                 imo=-1;
263
264                 if (hasMother) {
265                     imo=iparticle->GetFirstMother();
266                     TParticle* mother= (TParticle *) particles->At(imo);
267                     imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
268                 }
269 // Put particle on the stack
270                 gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
271                                  tof,"Secondary",nt);
272
273 //              printf("\n set track final: %d %d %d",imo, kf, nt);
274                 gAlice->KeepTrack(nt);
275                 *(newPos+i)=nt;
276             } // selected
277         } // particle loop final state
278  
279         delete newPos;
280
281         printf("\n I've put %i particles on the stack \n",nc);
282         if (nc > 0) {
283             jev+=nc;
284             if (jev >= fNpart || fNpart == -1) {
285                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
286                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
287                 break;
288             }
289         }
290     } // event loop
291     fHijing->Rluget(50,-1);
292 }
293
294 Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
295 {
296 // Perform kinematic selection
297     Double_t px=particle->Px();
298     Double_t py=particle->Py();
299     Double_t pz=particle->Pz();
300     Double_t  e=particle->Energy();
301
302 //
303 //  transverse momentum cut    
304     Double_t pt=TMath::Sqrt(px*px+py*py);
305     if (pt > fPtMax || pt < fPtMin) 
306     {
307 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
308         return kFALSE;
309     }
310 //
311 // momentum cut
312     Double_t p=TMath::Sqrt(px*px+py*py+pz*pz);
313     if (p > fPMax || p < fPMin) 
314     {
315 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
316         return kFALSE;
317     }
318     
319 //
320 // theta cut
321     Double_t  theta = Double_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
322     if (theta > fThetaMax || theta < fThetaMin) 
323     {
324         
325 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
326         return kFALSE;
327     }
328
329 //
330 // rapidity cut
331     Double_t y;
332     if(e<=pz) y = 99;
333     else if (e<=-pz)  y = -99;
334     else y = 0.5*TMath::Log((e+pz)/(e-pz));
335     if (y > fYMax || y < fYMin)
336     {
337 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
338         return kFALSE;
339     }
340
341 //
342 // phi cut
343     Double_t phi=Double_t(TMath::ATan2(Double_t(py),Double_t(px)));
344     if (phi > fPhiMax || phi < fPhiMin)
345     {
346 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
347         return kFALSE;
348     }
349
350     return kTRUE;
351 }
352
353 void AliGenHijing::KeepFullEvent()
354 {
355     fKeep=1;
356 }
357
358 void AliGenHijing::EvaluateCrossSections()
359 {
360 //     Glauber Calculation of geometrical x-section
361 //
362     Float_t xTot=0.;          // barn
363     Float_t xTotHard=0.;      // barn 
364     Float_t xPart=0.;         // barn
365     Float_t xPartHard=0.;     // barn 
366     Float_t sigmaHard=0.1;    // mbarn
367     Float_t bMin=0.;
368     Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
369     const Float_t kdib=0.2;
370     Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;
371
372
373     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
374     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
375     Int_t i;
376     Float_t oldvalue=0.;
377     
378     for (i=0; i<kMax; i++)
379     {
380         Float_t xb=bMin+i*kdib;
381         Float_t ov;
382         ov=fHijing->Profile(xb);
383         Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
384         Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
385         xTot+=gb;
386         xTotHard+=gbh;
387         if (xb > fMinImpactParam && xb < fMaxImpactParam)
388         {
389             xPart+=gb;
390             xPartHard+=gbh;
391         }
392         
393         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
394         oldvalue=xTot;
395         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
396         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
397     }
398     printf("\n Total cross section (barn): %f \n",xTot);
399     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
400     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
401     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
402 }
403
404 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
405 {
406 //
407 // Looks recursively if one of the daughters has been selected
408 //
409 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
410     Int_t imin=-1;
411     Int_t imax=-1;
412     Int_t i;
413     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
414     Bool_t selected=kFALSE;
415     if (hasDaughters) {
416         imin=iparticle->GetFirstDaughter();
417         imax=iparticle->GetLastDaughter();       
418         for (i=imin; i<= imax; i++){
419             TParticle *  jparticle       = (TParticle *) particles->At(i);      
420             Int_t ip=jparticle->GetPdgCode();
421             if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {
422                 selected=kTRUE; break;
423             }
424             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
425         }
426     } else {
427         return kFALSE;
428     }
429
430     return selected;
431 }
432
433
434 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
435 {
436 // Select flavor of particle
437 // 0: all
438 // 4: charm and beauty
439 // 5: beauty
440     if (fFlavor == 0) return kTRUE;
441     
442     Int_t ifl=TMath::Abs(pid/100);
443     if (ifl > 10) ifl/=10;
444     return (fFlavor == ifl);
445 }
446
447 Bool_t AliGenHijing::Stable(TParticle*  particle)
448 {
449     Int_t kf = TMath::Abs(particle->GetPdgCode());
450     
451     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
452          
453     {
454         return kTRUE;
455     } else {
456         return kFALSE;
457     }
458 }
459
460 void AliGenHijing::MakeHeader()
461 {
462 // Builds the event header, to be called after each event
463     AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
464 //    header->SetDate(date);
465 //    header->SetRunNumber(run);
466 //    header->SetEventNumber(event);
467     header->SetNProduced(fHijing->GetNATT());
468     header->SetImpactParameter(fHijing->GetHINT1(19));
469     header->SetTotalEnergy(fHijing->GetEATT());
470     header->SetHardScatters(fHijing->GetJATT());
471     header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
472     header->SetCollisions(fHijing->GetN0(),
473                           fHijing->GetN01(),
474                           fHijing->GetN10(),
475                           fHijing->GetN11());
476 }
477
478 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
479 {
480 // Assignment operator
481     return *this;
482 }
483
484
485
486
487
488
489
490
491
492
493
494
495