639c895a7d4c50f148d1690e2d2c85d2fa9bd7d9
[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.29  2001/10/15 08:12:24  morsch
19 - Vertex smearing with truncated gaussian.
20 - Store triggered jet info before and after final state radiation into mc-heade
21
22 Revision 1.28  2001/10/08 11:55:25  morsch
23 Store 4-momenta of trigegred jets in event header.
24 Possibility to switch of initial and final state radiation.
25
26 Revision 1.27  2001/10/08 07:13:14  morsch
27 Add setter for minimum transverse momentum of triggered jet.
28
29 Revision 1.26  2001/10/04 08:12:24  morsch
30 Redefinition of stable condition.
31
32 Revision 1.25  2001/07/27 17:09:36  morsch
33 Use local SetTrack, KeepTrack and SetHighWaterMark methods
34 to delegate either to local stack or to stack owned by AliRun.
35 (Piotr Skowronski, A.M.)
36
37 Revision 1.24  2001/07/20 09:34:56  morsch
38 Count the number of spectator neutrons and protons and add information
39 to the event header. (Chiara Oppedisano)
40
41 Revision 1.23  2001/07/13 17:30:22  morsch
42 Derive from AliGenMC.
43
44 Revision 1.22  2001/06/11 13:09:23  morsch
45 - Store cross-Section and number of binary collisions as a function of impact parameter
46 - Pass AliGenHijingEventHeader to gAlice.
47
48 Revision 1.21  2001/02/28 17:35:24  morsch
49 Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
50
51 Revision 1.20  2001/02/14 15:50:40  hristov
52 The last particle in event marked using SetHighWaterMark
53
54 Revision 1.19  2000/12/21 16:24:06  morsch
55 Coding convention clean-up
56
57 Revision 1.18  2000/12/06 17:46:30  morsch
58 Avoid random numbers 1 and 0.
59
60 Revision 1.17  2000/12/04 11:22:03  morsch
61 Init of sRandom as in 1.15
62
63 Revision 1.16  2000/12/02 11:41:39  morsch
64 Use SetRandom() to initialize random number generator in constructor.
65
66 Revision 1.15  2000/11/30 20:29:02  morsch
67 Initialise static variable sRandom in constructor: sRandom = fRandom;
68
69 Revision 1.14  2000/11/30 07:12:50  alibrary
70 Introducing new Rndm and QA classes
71
72 Revision 1.13  2000/11/09 17:40:27  morsch
73 Possibility to select/unselect spectator protons and neutrons.
74 Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
75
76 Revision 1.12  2000/10/20 13:38:38  morsch
77 Debug printouts commented.
78
79 Revision 1.11  2000/10/20 13:22:26  morsch
80 - skip particle type 92 (string)
81 - Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
82   mesons.
83
84 Revision 1.10  2000/10/17 15:10:20  morsch
85 Write first all the parent particles to the stack and then the final state particles.
86
87 Revision 1.9  2000/10/17 13:38:59  morsch
88 Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
89
90 Revision 1.8  2000/10/17 12:46:31  morsch
91 Protect EvaluateCrossSections() against division by zero.
92
93 Revision 1.7  2000/10/02 21:28:06  fca
94 Removal of useless dependecies via forward declarations
95
96 Revision 1.6  2000/09/11 13:23:37  morsch
97 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
98 luget_hijing and luset_hijing.
99
100 Revision 1.5  2000/09/07 16:55:40  morsch
101 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
102
103 Revision 1.4  2000/07/11 18:24:56  fca
104 Coding convention corrections + few minor bug fixes
105
106 Revision 1.3  2000/06/30 12:08:36  morsch
107 In member data: char* replaced by TString, Init takes care of resizing the strings to
108 8 characters required by Hijing.
109
110 Revision 1.2  2000/06/15 14:15:05  morsch
111 Add possibility for heavy flavor selection: charm and beauty.
112
113 Revision 1.1  2000/06/09 20:47:27  morsch
114 AliGenerator interface class to HIJING using THijing (test version)
115
116 */
117
118
119
120 // Generator using HIJING as an external generator
121 // The main HIJING options are accessable for the user through this interface.
122 // Uses the THijing implementation of TGenerator.
123 //
124 // andreas.morsch@cern.ch
125
126 #include "AliGenHijing.h"
127 #include "AliGenHijingEventHeader.h"
128 #include "AliRun.h"
129 #include "AliPDG.h"
130
131 #include <TArrayI.h>
132 #include <TParticle.h>
133 #include <THijing.h>
134 #include <TGraph.h>
135 #include <TLorentzVector.h>
136
137
138  ClassImp(AliGenHijing)
139
140 AliGenHijing::AliGenHijing()
141                  :AliGenMC()
142 {
143 // Constructor
144 }
145
146 AliGenHijing::AliGenHijing(Int_t npart)
147     :AliGenMC(npart)
148 {
149 // Default PbPb collisions at 5. 5 TeV
150 //
151     SetEnergyCMS();
152     SetImpactParameterRange();
153     SetTarget();
154     SetProjectile();
155     fKeep       =  0;
156     fQuench     =  1;
157     fShadowing  =  1;
158     fTrigger    =  0;
159     fDecaysOff  =  1;
160     fEvaluate   =  0;
161     fSelectAll  =  0;
162     fFlavor     =  0;
163     fSpectators =  1;
164     fDsigmaDb   =  0;
165     fDnDb       =  0;
166     fPtMinJet   = -2.5;         
167     fRadiation  =  1;
168     fEventVertex.Set(3);
169 //
170 // Set random number generator   
171     sRandom = fRandom;
172     fHijing = 0;
173
174 }
175
176 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
177 {
178 // copy constructor
179 }
180
181
182 AliGenHijing::~AliGenHijing()
183 {
184 // Destructor
185     if ( fDsigmaDb) delete  fDsigmaDb;  
186     if ( fDnDb)     delete  fDnDb;  
187 }
188
189 void AliGenHijing::Init()
190 {
191 // Initialisation
192     fFrame.Resize(8);
193     fTarget.Resize(8);
194     fProjectile.Resize(8);
195     
196     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
197                       fAProjectile, fZProjectile, fATarget, fZTarget, 
198                       fMinImpactParam, fMaxImpactParam));
199
200     fHijing=(THijing*) fgMCEvGen;
201     fHijing->SetIHPR2(2,  fRadiation);
202     fHijing->SetIHPR2(3,  fTrigger);
203     fHijing->SetIHPR2(4,  fQuench);
204     fHijing->SetIHPR2(6,  fShadowing);
205     fHijing->SetIHPR2(12, fDecaysOff);    
206     fHijing->SetIHPR2(21, fKeep);
207     fHijing->SetHIPR1(10, fPtMinJet);   
208     fHijing->Initialize();
209
210     
211 //
212     if (fEvaluate) EvaluateCrossSections();
213 //
214 //
215 //  Initialize random generator
216 }
217
218 void AliGenHijing::Generate()
219 {
220 // Generate one event
221
222   Float_t polar[3]    =   {0,0,0};
223   Float_t origin[3]   =   {0,0,0};
224   Float_t origin0[3]  =   {0,0,0};
225   Float_t p[3], random[6];
226   Float_t tof;
227
228   static TClonesArray *particles;
229 //  converts from mm/c to s
230   const Float_t kconv = 0.001/2.999792458e8;
231 //
232   Int_t nt  = 0;
233   Int_t jev = 0;
234   Int_t j, kf, ks, imo;
235   kf = 0;
236     
237   if(!particles) particles = new TClonesArray("TParticle",10000);
238     
239   fTrials = 0;
240   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
241   if(fVertexSmear == kPerEvent) {
242       Float_t dv[3];
243       dv[2] = 1.e10;
244       while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
245           Rndm(random,6);
246           for (j=0; j < 3; j++) {
247               dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
248                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
249           }
250       }
251       for (j=0; j < 3; j++) origin0[j] += dv[j];
252   } else if (fVertexSmear == kPerTrack) {
253 //          fHijing->SetMSTP(151,0);
254       for (j = 0; j < 3; j++) {
255 //            fHijing->SetPARP(151+j, fOsigma[j]*10.);
256       }
257   }
258   while(1)
259   {
260 //    Generate one event
261 // --------------------------------------------------------------------------
262       fSpecn    = 0;  
263       fSpecp    = 0;
264 // --------------------------------------------------------------------------
265       fHijing->GenerateEvent();
266       fTrials++;
267       fHijing->ImportParticles(particles,"All");
268       Int_t np = particles->GetEntriesFast();
269       printf("\n **************************************************%d\n",np);
270       Int_t nc = 0;
271       if (np == 0 ) continue;
272       Int_t i;
273       Int_t * newPos = new Int_t[np];
274
275       for (i = 0; i < np; i++) *(newPos+i) = i;
276 //      Get event vertex
277 //
278       TParticle *  iparticle = (TParticle *) particles->At(0);
279       fEventVertex[0] = origin0[0];
280       fEventVertex[1] = origin0[1];     
281       fEventVertex[2] = origin0[2];
282       
283 //
284 //      First write parent particles
285 //
286
287       for (i = 0; i < np; i++) {
288           iparticle       = (TParticle *) particles->At(i);
289 // Is this a parent particle ?
290           if (Stable(iparticle)) continue;
291 //
292         Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
293         Bool_t  selected             =  kTRUE;
294         Bool_t  hasSelectedDaughters =  kFALSE;
295             
296             
297         kf        = iparticle->GetPdgCode();
298         ks        = iparticle->GetStatusCode();
299         if (kf == 92) continue;
300             
301         if (!fSelectAll) selected = KinematicSelection(iparticle, 0)&&SelectFlavor(kf);
302         hasSelectedDaughters = DaughtersSelection(iparticle, particles);
303 //
304 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
305 //
306         if (selected || hasSelectedDaughters) {
307             nc++;
308             p[0] = iparticle->Px();
309             p[1] = iparticle->Py();
310             p[2] = iparticle->Pz();
311             origin[0] = origin0[0]+iparticle->Vx()/10;
312             origin[1] = origin0[1]+iparticle->Vy()/10;
313             origin[2] = origin0[2]+iparticle->Vz()/10;
314             tof = kconv*iparticle->T();
315             imo = -1;
316             TParticle* mother = 0;
317             if (hasMother) {
318                 imo = iparticle->GetFirstMother();
319                 mother = (TParticle *) particles->At(imo);
320                 imo = (mother->GetPdgCode() != 92) ? imo =* (newPos+imo) : -1;
321                 
322             }
323 // Put particle on the stack ... 
324 //              printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
325             
326             SetTrack(0,imo,kf,p,origin,polar, tof,kPPrimary,nt);
327 // ... and keep it there
328             KeepTrack(nt);
329 //
330             *(newPos+i)=nt;
331         } // selected
332       } // particle loop parents
333 //
334 // Now write the final state particles
335 //
336
337       for (i = 0; i<np; i++) {
338         TParticle *  iparticle       = (TParticle *) particles->At(i);
339 // Is this a final state particle ?
340         if (!Stable(iparticle)) continue;
341
342         Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
343         Bool_t  selected             =  kTRUE;
344         kf        = iparticle->GetPdgCode();
345         ks        = iparticle->GetStatusCode();
346 // --------------------------------------------------------------------------
347 // Count spectator neutrons and protons
348         if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
349           if(kf == kNeutron) fSpecn += 1;
350           if(kf == kProton)  fSpecp += 1;
351         }
352 // --------------------------------------------------------------------------
353 //          
354         if (!fSelectAll) {
355             selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
356             if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
357                                                       && ks != 11);
358         }
359 //
360 // Put particle on the stack if selected
361 //
362         if (selected) {
363             nc++;
364             p[0] = iparticle->Px();
365             p[1] = iparticle->Py();
366             p[2] = iparticle->Pz();
367             origin[0] = origin0[0]+iparticle->Vx()/10;
368             origin[1] = origin0[1]+iparticle->Vy()/10;
369             origin[2] = origin0[2]+iparticle->Vz()/10;
370             tof = kconv*iparticle->T();
371             imo = -1;
372             TParticle* mother = 0;
373             if (hasMother) {
374                 imo = iparticle->GetFirstMother();
375                 mother = (TParticle *) particles->At(imo);
376                 imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
377             }   
378 // Put particle on the stack
379             SetTrack(fTrackIt,imo,kf,p,origin,polar,
380                                                      tof,kPNoProcess,nt);
381             KeepTrack(nt);
382             *(newPos+i)=nt;
383         } // selected
384       } // particle loop final state
385  
386       delete[] newPos;
387
388       printf("\n I've put %i particles on the stack \n",nc);
389       if (nc > 0) {
390           jev += nc;
391           if (jev >= fNpart || fNpart == -1) {
392               fKineBias = Float_t(fNpart)/Float_t(fTrials);
393               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
394               break;
395           }
396       }
397   } // event loop
398   MakeHeader();
399   
400   SetHighWaterMark(nt);
401 }
402
403 void AliGenHijing::KeepFullEvent()
404 {
405     fKeep=1;
406 }
407
408 void AliGenHijing::EvaluateCrossSections()
409 {
410 //     Glauber Calculation of geometrical x-section
411 //
412     Float_t xTot       = 0.;          // barn
413     Float_t xTotHard   = 0.;          // barn 
414     Float_t xPart      = 0.;          // barn
415     Float_t xPartHard  = 0.;          // barn 
416     Float_t sigmaHard  = 0.1;         // mbarn
417     Float_t bMin       = 0.;
418     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
419     const Float_t kdib = 0.2;
420     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
421
422
423     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
424     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
425     Int_t i;
426     Float_t oldvalue= 0.;
427
428     Float_t* b   = new Float_t[kMax];
429     Float_t* si1 = new Float_t[kMax];    
430     Float_t* si2 = new Float_t[kMax];    
431     
432     for (i = 0; i < kMax; i++)
433     {
434         Float_t xb  = bMin+i*kdib;
435         Float_t ov;
436         ov=fHijing->Profile(xb);
437         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
438         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
439         xTot+=gb;
440         xTotHard += gbh;
441         if (xb > fMinImpactParam && xb < fMaxImpactParam)
442         {
443             xPart += gb;
444             xPartHard += gbh;
445         }
446         
447         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
448         oldvalue = xTot;
449         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
450         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
451         if (i>0) {
452             si1[i] = gb/kdib;
453             si2[i] = gbh/gb;
454             b[i]  = xb;
455         }
456     }
457
458     printf("\n Total cross section (barn): %f \n",xTot);
459     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
460     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
461     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
462
463 //  Store result as a graph
464     b[0] = 0;
465     si1[0] = 0;
466     si2[0]=si2[1];
467     
468     fDsigmaDb  = new TGraph(i, b, si1);
469     fDnDb      = new TGraph(i, b, si2);
470 }
471
472 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
473 {
474 //
475 // Looks recursively if one of the daughters has been selected
476 //
477 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
478     Int_t imin = -1;
479     Int_t imax = -1;
480     Int_t i;
481     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
482     Bool_t selected = kFALSE;
483     if (hasDaughters) {
484         imin = iparticle->GetFirstDaughter();
485         imax = iparticle->GetLastDaughter();       
486         for (i = imin; i <= imax; i++){
487             TParticle *  jparticle = (TParticle *) particles->At(i);    
488             Int_t ip = jparticle->GetPdgCode();
489             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
490                 selected=kTRUE; break;
491             }
492             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
493         }
494     } else {
495         return kFALSE;
496     }
497     return selected;
498 }
499
500
501 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
502 {
503 // Select flavor of particle
504 // 0: all
505 // 4: charm and beauty
506 // 5: beauty
507     if (fFlavor == 0) return kTRUE;
508     
509     Int_t ifl = TMath::Abs(pid/100);
510     if (ifl > 10) ifl/=10;
511     return (fFlavor == ifl);
512 }
513
514 Bool_t AliGenHijing::Stable(TParticle*  particle)
515 {
516 // Return true for a stable particle
517 //
518     if (particle->GetFirstDaughter() < 0 )
519     {
520         return kTRUE;
521     } else {
522         return kFALSE;
523     }
524 }
525
526 void AliGenHijing::MakeHeader()
527 {
528 // Builds the event header, to be called after each event
529     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
530     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
531     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
532     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
533     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
534     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
535     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
536                                                        fHijing->GetN01(),
537                                                        fHijing->GetN10(),
538                                                        fHijing->GetN11());
539     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
540
541 // 4-momentum vectors of the triggered jets.
542 //
543 // Before final state gluon radiation.
544     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
545                                               fHijing->GetHINT1(22),
546                                               fHijing->GetHINT1(23),
547                                               fHijing->GetHINT1(24));
548
549     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
550                                               fHijing->GetHINT1(32),
551                                               fHijing->GetHINT1(33),
552                                               fHijing->GetHINT1(34));
553 // After final state gluon radiation.
554     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
555                                               fHijing->GetHINT1(27),
556                                               fHijing->GetHINT1(28),
557                                               fHijing->GetHINT1(29));
558
559     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
560                                               fHijing->GetHINT1(37),
561                                               fHijing->GetHINT1(38),
562                                               fHijing->GetHINT1(39));
563     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
564 // Event Vertex
565     header->SetPrimaryVertex(fEventVertex);
566     gAlice->SetGenEventHeader(header);    
567 }
568
569 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
570 {
571 // Assignment operator
572     return *this;
573 }
574
575 #ifndef WIN32
576 # define rluget_hijing rluget_hijing_
577 # define rluset_hijing rluset_hijing_
578 # define rlu_hijing    rlu_hijing_
579 # define type_of_call
580 #else
581 # define rluget_hijing RLUGET_HIJING
582 # define rluset_hijing RLUSET_HIJING
583 # define rlu_hijing    RLU_HIJING
584 # define type_of_call _stdcall
585 #endif
586
587
588 extern "C" {
589   void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
590   {printf("Dummy version of rluget_hijing reached\n");}
591
592   void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
593   {printf("Dummy version of rluset_hijing reached\n");}
594
595   Double_t type_of_call rlu_hijing(Int_t & /*idum*/) 
596   {
597       Float_t r;
598       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
599       return r;
600   }
601 }