Changes related to the initialization of random numbers generators. Now one can use...
[u/mrichter/AliRoot.git] / THijing / 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 /* $Id$ */
17
18 // Generator using HIJING as an external generator
19 // The main HIJING options are accessable for the user through this interface.
20 // Uses the THijing implementation of TGenerator.
21 //
22 // andreas.morsch@cern.ch
23
24 #include <TArrayI.h>
25 #include <TGraph.h>
26 #include <THijing.h>
27 #include <TLorentzVector.h>
28 #include <TPDGCode.h>
29 #include <TParticle.h>
30
31 #include "AliGenHijing.h"
32 #include "AliGenHijingEventHeader.h"
33 #include "AliRun.h"
34 #include "AliHijingRndm.h"
35
36 ClassImp(AliGenHijing)
37
38 AliGenHijing::AliGenHijing()
39                  :AliGenMC()
40 {
41 // Constructor
42     fParticles = 0;
43     fHijing    = 0;
44     fDsigmaDb  = 0;
45     fDnDb      = 0;
46     AliHijingRndm::SetHijingRandom(GetRandom());
47 }
48
49 AliGenHijing::AliGenHijing(Int_t npart)
50     :AliGenMC(npart)
51 {
52 // Default PbPb collisions at 5. 5 TeV
53 //
54     fName = "Hijing";
55     fTitle= "Particle Generator using HIJING";
56
57     SetEnergyCMS();
58     SetImpactParameterRange();
59     SetBoostLHC();
60     SetJetEtaRange();
61     SetJetPhiRange();
62     
63     fKeep       =  0;
64     fQuench     =  1;
65     fShadowing  =  1;
66     fTrigger    =  0;
67     fDecaysOff  =  1;
68     fEvaluate   =  0;
69     fSelectAll  =  0;
70     fFlavor     =  0;
71     fSpectators =  1;
72     fDsigmaDb   =  0;
73     fDnDb       =  0;
74     fPtMinJet   = -2.5;         
75     fRadiation  =  3;
76     fEventVertex.Set(3);
77 //
78     SetSimpleJets();
79     SetNoGammas();
80 //
81     fParticles = new TClonesArray("TParticle",10000);    
82 //
83 // Set random number generator   
84     AliHijingRndm::SetHijingRandom(GetRandom());
85     fHijing = 0;
86
87 }
88
89 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
90 {
91 // copy constructor
92 }
93
94
95 AliGenHijing::~AliGenHijing()
96 {
97 // Destructor
98     if ( fDsigmaDb) delete  fDsigmaDb;  
99     if ( fDnDb)     delete  fDnDb;  
100     delete fParticles;
101 }
102
103 void AliGenHijing::Init()
104 {
105 // Initialisation
106     fFrame.Resize(8);
107     fTarget.Resize(8);
108     fProjectile.Resize(8);
109     
110     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
111                       fAProjectile, fZProjectile, fATarget, fZTarget, 
112                       fMinImpactParam, fMaxImpactParam));
113
114     fHijing=(THijing*) fgMCEvGen;
115     fHijing->SetIHPR2(2,  fRadiation);
116     fHijing->SetIHPR2(3,  fTrigger);
117     fHijing->SetIHPR2(6,  fShadowing);
118     fHijing->SetIHPR2(12, fDecaysOff);    
119     fHijing->SetIHPR2(21, fKeep);
120     fHijing->SetHIPR1(10, fPtMinJet);   
121     fHijing->SetHIPR1(50, fSimpleJet);
122 //
123 //  Quenching
124 //
125 //
126 //  fQuench = 0:  no quenching
127 //  fQuench = 1:  hijing default
128 //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
129 //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
130 //  fQuench = 4:  new LHC  parameters with log(e) dependence
131 //  fQuench = 5:  new RHIC parameters with log(e) dependence
132     fHijing->SetIHPR2(50, 0);
133     if (fQuench > 0) 
134         fHijing->SetIHPR2(4,  1);
135     else
136         fHijing->SetIHPR2(4,  0);
137 // New LHC parameters from Xin-Nian Wang
138     if (fQuench == 2) {
139         fHijing->SetHIPR1(14, 1.1);
140         fHijing->SetHIPR1(11, 3.7);
141     } else if (fQuench == 3) {
142         fHijing->SetHIPR1(14, 0.20);
143         fHijing->SetHIPR1(11, 2.5);
144     } else if (fQuench == 4) {
145         fHijing->SetIHPR2(50, 1);
146         fHijing->SetHIPR1(14, 4.*0.34);
147         fHijing->SetHIPR1(11, 3.7);
148     } else if (fQuench == 5) {
149         fHijing->SetIHPR2(50, 1);
150         fHijing->SetHIPR1(14, 0.34);
151         fHijing->SetHIPR1(11, 2.5);
152     }
153     
154     
155     AliGenMC::Init();
156     
157 //
158 //  Initialize Hijing  
159 //    
160     fHijing->Initialize();
161 //
162     if (fEvaluate) EvaluateCrossSections();
163 //
164 }
165
166 void AliGenHijing::Generate()
167 {
168 // Generate one event
169
170   Float_t polar[3]    =   {0,0,0};
171   Float_t origin[3]   =   {0,0,0};
172   Float_t origin0[3]  =   {0,0,0};
173   Float_t p[3], random[6];
174   Float_t tof;
175
176 //  converts from mm/c to s
177   const Float_t kconv = 0.001/2.999792458e8;
178 //
179   Int_t nt  = 0;
180   Int_t jev = 0;
181   Int_t j, kf, ks, imo;
182   kf = 0;
183     
184
185     
186   fTrials = 0;
187   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
188   if(fVertexSmear == kPerEvent) {
189       Float_t dv[3];
190       dv[2] = 1.e10;
191       while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
192           Rndm(random,6);
193           for (j=0; j < 3; j++) {
194               dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
195                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
196           }
197       }
198       for (j=0; j < 3; j++) origin0[j] += dv[j];
199   } else if (fVertexSmear == kPerTrack) {
200 //          fHijing->SetMSTP(151,0);
201       for (j = 0; j < 3; j++) {
202 //            fHijing->SetPARP(151+j, fOsigma[j]*10.);
203       }
204   }
205   while(1)
206   {
207 //    Generate one event
208 // --------------------------------------------------------------------------
209       fSpecn    = 0;  
210       fSpecp    = 0;
211 // --------------------------------------------------------------------------
212       fHijing->GenerateEvent();
213       fTrials++;
214       fHijing->ImportParticles(fParticles,"All");
215       if (fTrigger != kNoTrigger) {
216           if (!CheckTrigger()) continue;
217       }
218       if (fLHC) Boost();
219       
220       
221       Int_t np = fParticles->GetEntriesFast();
222       printf("\n **************************************************%d\n",np);
223       Int_t nc = 0;
224       if (np == 0 ) continue;
225       Int_t i;
226       Int_t* newPos     = new Int_t[np];
227       Int_t* pSelected  = new Int_t[np];
228
229       for (i = 0; i < np; i++) {
230           newPos[i]    = i;
231           pSelected[i] = 0;
232       }
233       
234 //      Get event vertex
235 //
236       TParticle *  iparticle = (TParticle *) fParticles->At(0);
237       fEventVertex[0] = origin0[0];
238       fEventVertex[1] = origin0[1];     
239       fEventVertex[2] = origin0[2];
240       
241 //
242 //      First select parent particles
243 //
244
245       for (i = 0; i < np; i++) {
246           iparticle = (TParticle *) fParticles->At(i);
247
248 // Is this a parent particle ?
249           if (Stable(iparticle)) continue;
250 //
251           Bool_t  selected             =  kTRUE;
252           Bool_t  hasSelectedDaughters =  kFALSE;
253           
254           
255           kf        = iparticle->GetPdgCode();
256           ks        = iparticle->GetStatusCode();
257           if (kf == 92) continue;
258             
259           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
260                                SelectFlavor(kf);
261           hasSelectedDaughters = DaughtersSelection(iparticle);
262 //
263 // Put particle on the stack if it is either selected or 
264 // it is the mother of at least one seleted particle
265 //
266           if (selected || hasSelectedDaughters) {
267               nc++;
268               pSelected[i] = 1;
269           } // selected
270       } // particle loop parents
271 //
272 // Now select the final state particles
273 //
274
275       for (i = 0; i<np; i++) {
276           TParticle *  iparticle = (TParticle *) fParticles->At(i);
277 // Is this a final state particle ?
278           if (!Stable(iparticle)) continue;
279       
280           Bool_t  selected             =  kTRUE;
281           kf        = iparticle->GetPdgCode();
282           ks        = iparticle->GetStatusCode();
283           
284 // --------------------------------------------------------------------------
285 // Count spectator neutrons and protons
286           if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
287               if(kf == kNeutron) fSpecn += 1;
288               if(kf == kProton)  fSpecp += 1;
289           }
290 // --------------------------------------------------------------------------
291 //          
292           if (!fSelectAll) {
293               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
294               if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
295                                                         && ks != 11);
296           }
297 //
298 // Put particle on the stack if selected
299 //
300           if (selected) {
301               nc++;
302               pSelected[i] = 1;
303           } // selected
304       } // particle loop final state
305 //
306 // Write particles to stack
307 //
308       for (i = 0; i<np; i++) {
309           TParticle *  iparticle = (TParticle *) fParticles->At(i);
310           Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
311           Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
312
313           if (pSelected[i]) {
314               kf   = iparticle->GetPdgCode();
315               ks   = iparticle->GetStatusCode();
316               p[0] = iparticle->Px();
317               p[1] = iparticle->Py();
318               p[2] = iparticle->Pz();
319               origin[0] = origin0[0]+iparticle->Vx()/10;
320               origin[1] = origin0[1]+iparticle->Vy()/10;
321               origin[2] = origin0[2]+iparticle->Vz()/10;
322               tof = kconv*iparticle->T();
323               imo = -1;
324               TParticle* mother = 0;
325               if (hasMother) {
326                   imo = iparticle->GetFirstMother();
327                   mother = (TParticle *) fParticles->At(imo);
328                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
329               } // if has mother   
330               Bool_t tFlag = (fTrackIt && !hasDaughter);
331               SetTrack(tFlag,imo,kf,p,origin,polar,
332                        tof,kPNoProcess,nt, 1., ks);
333               KeepTrack(nt);
334               newPos[i] = nt;
335           } // if selected
336       } // particle loop
337       delete[] newPos;
338       delete[] pSelected;
339       
340       printf("\n I've put %i particles on the stack \n",nc);
341       if (nc > 0) {
342           jev += nc;
343           if (jev >= fNpart || fNpart == -1) {
344               fKineBias = Float_t(fNpart)/Float_t(fTrials);
345               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
346               break;
347           }
348       }
349   } // event loop
350   MakeHeader();
351   SetHighWaterMark(nt);
352 }
353
354 void AliGenHijing::KeepFullEvent()
355 {
356     fKeep=1;
357 }
358
359 void AliGenHijing::EvaluateCrossSections()
360 {
361 //     Glauber Calculation of geometrical x-section
362 //
363     Float_t xTot       = 0.;          // barn
364     Float_t xTotHard   = 0.;          // barn 
365     Float_t xPart      = 0.;          // barn
366     Float_t xPartHard  = 0.;          // barn 
367     Float_t sigmaHard  = 0.1;         // mbarn
368     Float_t bMin       = 0.;
369     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
370     const Float_t kdib = 0.2;
371     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
372
373
374     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
375     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
376     Int_t i;
377     Float_t oldvalue= 0.;
378
379     Float_t* b   = new Float_t[kMax];
380     Float_t* si1 = new Float_t[kMax];    
381     Float_t* si2 = new Float_t[kMax];    
382     
383     for (i = 0; i < kMax; i++)
384     {
385         Float_t xb  = bMin+i*kdib;
386         Float_t ov;
387         ov=fHijing->Profile(xb);
388         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
389         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
390         xTot+=gb;
391         xTotHard += gbh;
392         printf("profile %f %f %f\n", xb, ov, fHijing->GetHINT1(12));
393         
394         if (xb > fMinImpactParam && xb < fMaxImpactParam)
395         {
396             xPart += gb;
397             xPartHard += gbh;
398         }
399         
400         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
401         oldvalue = xTot;
402         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
403         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
404         if (i>0) {
405             si1[i] = gb/kdib;
406             si2[i] = gbh/gb;
407             b[i]  = xb;
408         }
409     }
410
411     printf("\n Total cross section (barn): %f \n",xTot);
412     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
413     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
414     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
415
416 //  Store result as a graph
417     b[0] = 0;
418     si1[0] = 0;
419     si2[0]=si2[1];
420     
421     fDsigmaDb  = new TGraph(i, b, si1);
422     fDnDb      = new TGraph(i, b, si2);
423 }
424
425 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
426 {
427 //
428 // Looks recursively if one of the daughters has been selected
429 //
430 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
431     Int_t imin = -1;
432     Int_t imax = -1;
433     Int_t i;
434     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
435     Bool_t selected = kFALSE;
436     if (hasDaughters) {
437         imin = iparticle->GetFirstDaughter();
438         imax = iparticle->GetLastDaughter();       
439         for (i = imin; i <= imax; i++){
440             TParticle *  jparticle = (TParticle *) fParticles->At(i);   
441             Int_t ip = jparticle->GetPdgCode();
442             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
443                 selected=kTRUE; break;
444             }
445             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
446         }
447     } else {
448         return kFALSE;
449     }
450     return selected;
451 }
452
453
454 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
455 {
456 // Select flavor of particle
457 // 0: all
458 // 4: charm and beauty
459 // 5: beauty
460     Bool_t res = 0;
461     
462     if (fFlavor == 0) {
463         res = kTRUE;
464     } else {
465         Int_t ifl = TMath::Abs(pid/100);
466         if (ifl > 10) ifl/=10;
467         res = (fFlavor == ifl);
468     }
469 //
470 //  This part if gamma writing is inhibited
471     if (fNoGammas) 
472         res = res && (pid != kGamma && pid != kPi0);
473 //
474     return res;
475 }
476
477 Bool_t AliGenHijing::Stable(TParticle*  particle)
478 {
479 // Return true for a stable particle
480 //
481     
482     if (particle->GetFirstDaughter() < 0 )
483     {
484         return kTRUE;
485     } else {
486         return kFALSE;
487     }
488 }
489
490
491
492 void AliGenHijing::MakeHeader()
493 {
494 // Builds the event header, to be called after each event
495     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
496     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
497     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
498     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
499     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
500     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
501     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
502                                                        fHijing->GetN01(),
503                                                        fHijing->GetN10(),
504                                                        fHijing->GetN11());
505     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
506
507 // 4-momentum vectors of the triggered jets.
508 //
509 // Before final state gluon radiation.
510     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
511                                               fHijing->GetHINT1(22),
512                                               fHijing->GetHINT1(23),
513                                               fHijing->GetHINT1(24));
514
515     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
516                                               fHijing->GetHINT1(32),
517                                               fHijing->GetHINT1(33),
518                                               fHijing->GetHINT1(34));
519 // After final state gluon radiation.
520     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
521                                               fHijing->GetHINT1(27),
522                                               fHijing->GetHINT1(28),
523                                               fHijing->GetHINT1(29));
524
525     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
526                                               fHijing->GetHINT1(37),
527                                               fHijing->GetHINT1(38),
528                                               fHijing->GetHINT1(39));
529     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
530 // Bookkeeping for kinematic bias
531     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
532 // Event Vertex
533     header->SetPrimaryVertex(fEventVertex);
534     gAlice->SetGenEventHeader(header);   
535     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
536 }
537
538 Bool_t AliGenHijing::CheckTrigger()
539 {
540 // Check the kinematic trigger condition
541 //
542     Bool_t   triggered = kFALSE;
543  
544     if (fTrigger == 1) {
545 //
546 //  jet-jet Trigger     
547         
548         TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), 
549                                                   fHijing->GetHINT1(27),
550                                                   fHijing->GetHINT1(28),
551                                                   fHijing->GetHINT1(29));
552         
553         TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), 
554                                                   fHijing->GetHINT1(37),
555                                                   fHijing->GetHINT1(38),
556                                                   fHijing->GetHINT1(39));
557         Double_t eta1      = jet1->Eta();
558         Double_t eta2      = jet2->Eta();
559         Double_t phi1      = jet1->Phi();
560         Double_t phi2      = jet2->Phi();
561 //    printf("\n Trigger: %f %f %f %f",
562 //         fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
563         if (
564             (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
565              phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
566             ||
567             (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
568              phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
569             ) 
570             triggered = kTRUE;
571     } else if (fTrigger == 2) {
572 //  Gamma Jet
573 //
574         Int_t np = fParticles->GetEntriesFast();
575         for (Int_t i = 0; i < np; i++) {
576             TParticle* part = (TParticle*) fParticles->At(i);
577             Int_t kf = part->GetPdgCode();
578             Int_t ks = part->GetStatusCode();
579             if (kf == 22 && ks == 40) {
580                 Float_t phi = part->Phi();
581                 Float_t eta = part->Eta();
582                 if  (eta < fEtaMaxJet && 
583                      eta > fEtaMinJet &&
584                      phi < fPhiMaxJet && 
585                      phi > fPhiMinJet) {
586                     triggered = 1;
587                     break;
588                 } // check phi,eta within limits
589             } // direct gamma ? 
590         } // particle loop
591     } // fTrigger == 2
592     return triggered;
593 }
594
595
596
597
598 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
599 {
600 // Assignment operator
601     return *this;
602 }
603