Provides collision geometry.
[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.45  2002/12/16 09:44:49  morsch
19 Default for fRadiation is 3.
20
21 Revision 1.44  2002/10/14 14:55:35  hristov
22 Merging the VirtualMC branch to the main development branch (HEAD)
23
24 Revision 1.42.4.1  2002/08/28 15:06:50  alibrary
25 Updating to v3-09-01
26
27 Revision 1.43  2002/08/09 12:09:52  morsch
28 Direct gamma trigger correctly included.
29
30 Revision 1.42  2002/03/12 11:07:08  morsch
31 Add status code of particle to SetTrack call.
32
33 Revision 1.41  2002/03/05 11:25:33  morsch
34 - New quenching options
35 - Correction in CheckTrigger()
36
37 Revision 1.40  2002/02/12 11:05:53  morsch
38 Get daughter indices right.
39
40 Revision 1.39  2002/02/12 09:16:39  morsch
41 Correction in SelectFlavor()
42
43 Revision 1.38  2002/02/12 08:53:21  morsch
44 SetNoGammas can be used to inhibit writing of gammas and pi0.
45
46 Revision 1.37  2002/02/08 16:50:50  morsch
47 Add name and title in constructor.
48
49 Revision 1.36  2002/01/31 20:17:55  morsch
50 Allow for triggered jets with simplified topology: Exact pT, back-to-back
51
52 Revision 1.35  2001/12/13 07:56:25  hristov
53 Set pointers to zero in the default constructor
54
55 Revision 1.34  2001/12/11 16:55:42  morsch
56 Correct initialization for jet phi-range.
57
58 Revision 1.33  2001/12/05 10:18:51  morsch
59 Possibility of kinematic biasing of jet phi range. (J. Klay)
60
61 Revision 1.32  2001/11/28 13:51:11  morsch
62 Introduce kinematic biasing (etamin, etamax) of jet trigger. Bookkeeping
63 (number of trials) done in AliGenHijingEventHeader.
64
65 Revision 1.31  2001/11/06 12:30:34  morsch
66 Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric collision systems.
67
68 Revision 1.30  2001/10/21 18:35:56  hristov
69 Several pointers were set to zero in the default constructors to avoid memory management problems
70
71 Revision 1.29  2001/10/15 08:12:24  morsch
72 - Vertex smearing with truncated gaussian.
73 - Store triggered jet info before and after final state radiation into mc-heade
74
75 Revision 1.28  2001/10/08 11:55:25  morsch
76 Store 4-momenta of trigegred jets in event header.
77 Possibility to switch of initial and final state radiation.
78
79 Revision 1.27  2001/10/08 07:13:14  morsch
80 Add setter for minimum transverse momentum of triggered jet.
81
82 Revision 1.26  2001/10/04 08:12:24  morsch
83 Redefinition of stable condition.
84
85 Revision 1.25  2001/07/27 17:09:36  morsch
86 Use local SetTrack, KeepTrack and SetHighWaterMark methods
87 to delegate either to local stack or to stack owned by AliRun.
88 (Piotr Skowronski, A.M.)
89
90 Revision 1.24  2001/07/20 09:34:56  morsch
91 Count the number of spectator neutrons and protons and add information
92 to the event header. (Chiara Oppedisano)
93
94 Revision 1.23  2001/07/13 17:30:22  morsch
95 Derive from AliGenMC.
96
97 Revision 1.22  2001/06/11 13:09:23  morsch
98 - Store cross-Section and number of binary collisions as a function of impact parameter
99 - Pass AliGenHijingEventHeader to gAlice.
100
101 Revision 1.21  2001/02/28 17:35:24  morsch
102 Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
103
104 Revision 1.20  2001/02/14 15:50:40  hristov
105 The last particle in event marked using SetHighWaterMark
106
107 Revision 1.19  2000/12/21 16:24:06  morsch
108 Coding convention clean-up
109
110 Revision 1.18  2000/12/06 17:46:30  morsch
111 Avoid random numbers 1 and 0.
112
113 Revision 1.17  2000/12/04 11:22:03  morsch
114 Init of sRandom as in 1.15
115
116 Revision 1.16  2000/12/02 11:41:39  morsch
117 Use SetRandom() to initialize random number generator in constructor.
118
119 Revision 1.15  2000/11/30 20:29:02  morsch
120 Initialise static variable sRandom in constructor: sRandom = fRandom;
121
122 Revision 1.14  2000/11/30 07:12:50  alibrary
123 Introducing new Rndm and QA classes
124
125 Revision 1.13  2000/11/09 17:40:27  morsch
126 Possibility to select/unselect spectator protons and neutrons.
127 Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
128
129 Revision 1.12  2000/10/20 13:38:38  morsch
130 Debug printouts commented.
131
132 Revision 1.11  2000/10/20 13:22:26  morsch
133 - skip particle type 92 (string)
134 - Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
135   mesons.
136
137 Revision 1.10  2000/10/17 15:10:20  morsch
138 Write first all the parent particles to the stack and then the final state particles.
139
140 Revision 1.9  2000/10/17 13:38:59  morsch
141 Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
142
143 Revision 1.8  2000/10/17 12:46:31  morsch
144 Protect EvaluateCrossSections() against division by zero.
145
146 Revision 1.7  2000/10/02 21:28:06  fca
147 Removal of useless dependecies via forward declarations
148
149 Revision 1.6  2000/09/11 13:23:37  morsch
150 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
151 luget_hijing and luset_hijing.
152
153 Revision 1.5  2000/09/07 16:55:40  morsch
154 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
155
156 Revision 1.4  2000/07/11 18:24:56  fca
157 Coding convention corrections + few minor bug fixes
158
159 Revision 1.3  2000/06/30 12:08:36  morsch
160 In member data: char* replaced by TString, Init takes care of resizing the strings to
161 8 characters required by Hijing.
162
163 Revision 1.2  2000/06/15 14:15:05  morsch
164 Add possibility for heavy flavor selection: charm and beauty.
165
166 Revision 1.1  2000/06/09 20:47:27  morsch
167 AliGenerator interface class to HIJING using THijing (test version)
168
169 */
170
171
172
173 // Generator using HIJING as an external generator
174 // The main HIJING options are accessable for the user through this interface.
175 // Uses the THijing implementation of TGenerator.
176 //
177 // andreas.morsch@cern.ch
178
179 #include "AliGenHijing.h"
180 #include "AliGenHijingEventHeader.h"
181 #include "AliRun.h"
182 #include "AliPDG.h"
183
184 #include <TArrayI.h>
185 #include <TParticle.h>
186 #include <THijing.h>
187 #include <TGraph.h>
188 #include <TLorentzVector.h>
189
190
191  ClassImp(AliGenHijing)
192
193 AliGenHijing::AliGenHijing()
194                  :AliGenMC()
195 {
196 // Constructor
197     fParticles = 0;
198     fHijing    = 0;
199     fDsigmaDb  = 0;
200     fDnDb      = 0;
201 }
202
203 AliGenHijing::AliGenHijing(Int_t npart)
204     :AliGenMC(npart)
205 {
206 // Default PbPb collisions at 5. 5 TeV
207 //
208     fName = "Hijing";
209     fTitle= "Particle Generator using HIJING";
210
211     SetEnergyCMS();
212     SetImpactParameterRange();
213     SetTarget();
214     SetProjectile();
215     SetBoostLHC();
216     SetJetEtaRange();
217     SetJetPhiRange();
218     
219     fKeep       =  0;
220     fQuench     =  1;
221     fShadowing  =  1;
222     fTrigger    =  0;
223     fDecaysOff  =  1;
224     fEvaluate   =  0;
225     fSelectAll  =  0;
226     fFlavor     =  0;
227     fSpectators =  1;
228     fDsigmaDb   =  0;
229     fDnDb       =  0;
230     fPtMinJet   = -2.5;         
231     fRadiation  =  3;
232     fEventVertex.Set(3);
233 //
234     SetSimpleJets();
235     SetNoGammas();
236 //
237     fParticles = new TClonesArray("TParticle",10000);    
238 //
239 // Set random number generator   
240     sRandom = fRandom;
241     fHijing = 0;
242
243 }
244
245 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
246 {
247 // copy constructor
248 }
249
250
251 AliGenHijing::~AliGenHijing()
252 {
253 // Destructor
254     if ( fDsigmaDb) delete  fDsigmaDb;  
255     if ( fDnDb)     delete  fDnDb;  
256     delete fParticles;
257 }
258
259 void AliGenHijing::Init()
260 {
261 // Initialisation
262     fFrame.Resize(8);
263     fTarget.Resize(8);
264     fProjectile.Resize(8);
265     
266     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
267                       fAProjectile, fZProjectile, fATarget, fZTarget, 
268                       fMinImpactParam, fMaxImpactParam));
269
270     fHijing=(THijing*) fgMCEvGen;
271     fHijing->SetIHPR2(2,  fRadiation);
272     fHijing->SetIHPR2(3,  fTrigger);
273     fHijing->SetIHPR2(6,  fShadowing);
274     fHijing->SetIHPR2(12, fDecaysOff);    
275     fHijing->SetIHPR2(21, fKeep);
276     fHijing->SetHIPR1(10, fPtMinJet);   
277     fHijing->SetHIPR1(50, fSimpleJet);
278 //
279 //  Quenching
280 //
281 //
282 //  fQuench = 0:  no quenching
283 //  fQuench = 1:  hijing default
284 //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
285 //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
286 //  fQuench = 4:  new LHC  parameters with log(e) dependence
287 //  fQuench = 5:  new RHIC parameters with log(e) dependence
288     fHijing->SetIHPR2(50, 0);
289     if (fQuench > 0) 
290         fHijing->SetIHPR2(4,  1);
291     else
292         fHijing->SetIHPR2(4,  0);
293 // New LHC parameters from Xin-Nian Wang
294     if (fQuench == 2) {
295         fHijing->SetHIPR1(14, 1.1);
296         fHijing->SetHIPR1(11, 3.7);
297     } else if (fQuench == 3) {
298         fHijing->SetHIPR1(14, 0.20);
299         fHijing->SetHIPR1(11, 2.5);
300     } else if (fQuench == 4) {
301         fHijing->SetIHPR2(50, 1);
302         fHijing->SetHIPR1(14, 4.*0.34);
303         fHijing->SetHIPR1(11, 3.7);
304     } else if (fQuench == 5) {
305         fHijing->SetIHPR2(50, 1);
306         fHijing->SetHIPR1(14, 0.34);
307         fHijing->SetHIPR1(11, 2.5);
308     }
309     
310     
311     
312 //
313 //  Initialize Hijing  
314 //    
315     fHijing->Initialize();
316 //
317     if (fEvaluate) EvaluateCrossSections();
318 //
319 }
320
321 void AliGenHijing::Generate()
322 {
323 // Generate one event
324
325   Float_t polar[3]    =   {0,0,0};
326   Float_t origin[3]   =   {0,0,0};
327   Float_t origin0[3]  =   {0,0,0};
328   Float_t p[3], random[6];
329   Float_t tof;
330
331 //  converts from mm/c to s
332   const Float_t kconv = 0.001/2.999792458e8;
333 //
334   Int_t nt  = 0;
335   Int_t jev = 0;
336   Int_t j, kf, ks, imo;
337   kf = 0;
338     
339
340     
341   fTrials = 0;
342   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
343   if(fVertexSmear == kPerEvent) {
344       Float_t dv[3];
345       dv[2] = 1.e10;
346       while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
347           Rndm(random,6);
348           for (j=0; j < 3; j++) {
349               dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
350                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
351           }
352       }
353       for (j=0; j < 3; j++) origin0[j] += dv[j];
354   } else if (fVertexSmear == kPerTrack) {
355 //          fHijing->SetMSTP(151,0);
356       for (j = 0; j < 3; j++) {
357 //            fHijing->SetPARP(151+j, fOsigma[j]*10.);
358       }
359   }
360   while(1)
361   {
362 //    Generate one event
363 // --------------------------------------------------------------------------
364       fSpecn    = 0;  
365       fSpecp    = 0;
366 // --------------------------------------------------------------------------
367       fHijing->GenerateEvent();
368       fTrials++;
369       fHijing->ImportParticles(fParticles,"All");
370       if (fTrigger != kNoTrigger) {
371           if (!CheckTrigger()) continue;
372       }
373       if (fLHC) Boost();
374       
375       
376       Int_t np = fParticles->GetEntriesFast();
377       printf("\n **************************************************%d\n",np);
378       Int_t nc = 0;
379       if (np == 0 ) continue;
380       Int_t i;
381       Int_t* newPos     = new Int_t[np];
382       Int_t* pSelected  = new Int_t[np];
383
384       for (i = 0; i < np; i++) {
385           newPos[i]    = i;
386           pSelected[i] = 0;
387       }
388       
389 //      Get event vertex
390 //
391       TParticle *  iparticle = (TParticle *) fParticles->At(0);
392       fEventVertex[0] = origin0[0];
393       fEventVertex[1] = origin0[1];     
394       fEventVertex[2] = origin0[2];
395       
396 //
397 //      First select parent particles
398 //
399
400       for (i = 0; i < np; i++) {
401           iparticle = (TParticle *) fParticles->At(i);
402
403 // Is this a parent particle ?
404           if (Stable(iparticle)) continue;
405 //
406           Bool_t  selected             =  kTRUE;
407           Bool_t  hasSelectedDaughters =  kFALSE;
408           
409           
410           kf        = iparticle->GetPdgCode();
411           ks        = iparticle->GetStatusCode();
412           if (kf == 92) continue;
413             
414           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
415                                SelectFlavor(kf);
416           hasSelectedDaughters = DaughtersSelection(iparticle);
417 //
418 // Put particle on the stack if it is either selected or 
419 // it is the mother of at least one seleted particle
420 //
421           if (selected || hasSelectedDaughters) {
422               nc++;
423               pSelected[i] = 1;
424           } // selected
425       } // particle loop parents
426 //
427 // Now select the final state particles
428 //
429
430       for (i = 0; i<np; i++) {
431           TParticle *  iparticle = (TParticle *) fParticles->At(i);
432 // Is this a final state particle ?
433           if (!Stable(iparticle)) continue;
434       
435           Bool_t  selected             =  kTRUE;
436           kf        = iparticle->GetPdgCode();
437           ks        = iparticle->GetStatusCode();
438           
439 // --------------------------------------------------------------------------
440 // Count spectator neutrons and protons
441           if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
442               if(kf == kNeutron) fSpecn += 1;
443               if(kf == kProton)  fSpecp += 1;
444           }
445 // --------------------------------------------------------------------------
446 //          
447           if (!fSelectAll) {
448               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
449               if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
450                                                         && ks != 11);
451           }
452 //
453 // Put particle on the stack if selected
454 //
455           if (selected) {
456               nc++;
457               pSelected[i] = 1;
458           } // selected
459       } // particle loop final state
460 //
461 // Write particles to stack
462 //
463       for (i = 0; i<np; i++) {
464           TParticle *  iparticle = (TParticle *) fParticles->At(i);
465           Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
466           Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
467
468           if (pSelected[i]) {
469               kf   = iparticle->GetPdgCode();
470               ks   = iparticle->GetStatusCode();
471               p[0] = iparticle->Px();
472               p[1] = iparticle->Py();
473               p[2] = iparticle->Pz();
474               origin[0] = origin0[0]+iparticle->Vx()/10;
475               origin[1] = origin0[1]+iparticle->Vy()/10;
476               origin[2] = origin0[2]+iparticle->Vz()/10;
477               tof = kconv*iparticle->T();
478               imo = -1;
479               TParticle* mother = 0;
480               if (hasMother) {
481                   imo = iparticle->GetFirstMother();
482                   mother = (TParticle *) fParticles->At(imo);
483                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
484               } // if has mother   
485               Bool_t tFlag = (fTrackIt && !hasDaughter);
486               SetTrack(tFlag,imo,kf,p,origin,polar,
487                        tof,kPNoProcess,nt, 1., ks);
488               KeepTrack(nt);
489               newPos[i] = nt;
490           } // if selected
491       } // particle loop
492       delete[] newPos;
493       delete[] pSelected;
494       
495       printf("\n I've put %i particles on the stack \n",nc);
496       if (nc > 0) {
497           jev += nc;
498           if (jev >= fNpart || fNpart == -1) {
499               fKineBias = Float_t(fNpart)/Float_t(fTrials);
500               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
501               break;
502           }
503       }
504   } // event loop
505   MakeHeader();
506   SetHighWaterMark(nt);
507 }
508
509 void AliGenHijing::KeepFullEvent()
510 {
511     fKeep=1;
512 }
513
514 void AliGenHijing::EvaluateCrossSections()
515 {
516 //     Glauber Calculation of geometrical x-section
517 //
518     Float_t xTot       = 0.;          // barn
519     Float_t xTotHard   = 0.;          // barn 
520     Float_t xPart      = 0.;          // barn
521     Float_t xPartHard  = 0.;          // barn 
522     Float_t sigmaHard  = 0.1;         // mbarn
523     Float_t bMin       = 0.;
524     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
525     const Float_t kdib = 0.2;
526     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
527
528
529     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
530     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
531     Int_t i;
532     Float_t oldvalue= 0.;
533
534     Float_t* b   = new Float_t[kMax];
535     Float_t* si1 = new Float_t[kMax];    
536     Float_t* si2 = new Float_t[kMax];    
537     
538     for (i = 0; i < kMax; i++)
539     {
540         Float_t xb  = bMin+i*kdib;
541         Float_t ov;
542         ov=fHijing->Profile(xb);
543         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
544         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
545         xTot+=gb;
546         xTotHard += gbh;
547         if (xb > fMinImpactParam && xb < fMaxImpactParam)
548         {
549             xPart += gb;
550             xPartHard += gbh;
551         }
552         
553         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
554         oldvalue = xTot;
555         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
556         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
557         if (i>0) {
558             si1[i] = gb/kdib;
559             si2[i] = gbh/gb;
560             b[i]  = xb;
561         }
562     }
563
564     printf("\n Total cross section (barn): %f \n",xTot);
565     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
566     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
567     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
568
569 //  Store result as a graph
570     b[0] = 0;
571     si1[0] = 0;
572     si2[0]=si2[1];
573     
574     fDsigmaDb  = new TGraph(i, b, si1);
575     fDnDb      = new TGraph(i, b, si2);
576 }
577
578 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
579 {
580 //
581 // Looks recursively if one of the daughters has been selected
582 //
583 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
584     Int_t imin = -1;
585     Int_t imax = -1;
586     Int_t i;
587     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
588     Bool_t selected = kFALSE;
589     if (hasDaughters) {
590         imin = iparticle->GetFirstDaughter();
591         imax = iparticle->GetLastDaughter();       
592         for (i = imin; i <= imax; i++){
593             TParticle *  jparticle = (TParticle *) fParticles->At(i);   
594             Int_t ip = jparticle->GetPdgCode();
595             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
596                 selected=kTRUE; break;
597             }
598             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
599         }
600     } else {
601         return kFALSE;
602     }
603     return selected;
604 }
605
606
607 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
608 {
609 // Select flavor of particle
610 // 0: all
611 // 4: charm and beauty
612 // 5: beauty
613     Bool_t res = 0;
614     
615     if (fFlavor == 0) {
616         res = kTRUE;
617     } else {
618         Int_t ifl = TMath::Abs(pid/100);
619         if (ifl > 10) ifl/=10;
620         res = (fFlavor == ifl);
621     }
622 //
623 //  This part if gamma writing is inhibited
624     if (fNoGammas) 
625         res = res && (pid != kGamma && pid != kPi0);
626 //
627     return res;
628 }
629
630 Bool_t AliGenHijing::Stable(TParticle*  particle)
631 {
632 // Return true for a stable particle
633 //
634     
635     if (particle->GetFirstDaughter() < 0 )
636     {
637         return kTRUE;
638     } else {
639         return kFALSE;
640     }
641 }
642
643
644 void AliGenHijing::Boost()
645 {
646 //
647 // Boost cms into LHC lab frame
648 //
649     Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
650                                       (Double_t(fZTarget)    * Double_t(fAProjectile)));
651     Double_t beta  = TMath::TanH(dy);
652     Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
653     Double_t gb    = gamma * beta;
654
655     printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
656     
657     Int_t i;
658     Int_t np = fParticles->GetEntriesFast();
659     for (i = 0; i < np; i++) 
660     {
661         TParticle* iparticle = (TParticle*) fParticles->At(i);
662
663         Double_t e   = iparticle->Energy();
664         Double_t px  = iparticle->Px();
665         Double_t py  = iparticle->Py();
666         Double_t pz  = iparticle->Pz();
667
668         Double_t eb  = gamma * e -      gb * pz;
669         Double_t pzb =   -gb * e +   gamma * pz;
670
671         iparticle->SetMomentum(px, py, pzb, eb);
672     }
673 }
674
675
676 void AliGenHijing::MakeHeader()
677 {
678 // Builds the event header, to be called after each event
679     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
680     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
681     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
682     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
683     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
684     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
685     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
686                                                        fHijing->GetN01(),
687                                                        fHijing->GetN10(),
688                                                        fHijing->GetN11());
689     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
690
691 // 4-momentum vectors of the triggered jets.
692 //
693 // Before final state gluon radiation.
694     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
695                                               fHijing->GetHINT1(22),
696                                               fHijing->GetHINT1(23),
697                                               fHijing->GetHINT1(24));
698
699     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
700                                               fHijing->GetHINT1(32),
701                                               fHijing->GetHINT1(33),
702                                               fHijing->GetHINT1(34));
703 // After final state gluon radiation.
704     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
705                                               fHijing->GetHINT1(27),
706                                               fHijing->GetHINT1(28),
707                                               fHijing->GetHINT1(29));
708
709     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
710                                               fHijing->GetHINT1(37),
711                                               fHijing->GetHINT1(38),
712                                               fHijing->GetHINT1(39));
713     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
714 // Bookkeeping for kinematic bias
715     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
716 // Event Vertex
717     header->SetPrimaryVertex(fEventVertex);
718     gAlice->SetGenEventHeader(header);   
719     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
720 }
721
722 Bool_t AliGenHijing::CheckTrigger()
723 {
724 // Check the kinematic trigger condition
725 //
726     Bool_t   triggered = kFALSE;
727  
728     if (fTrigger == 1) {
729 //
730 //  jet-jet Trigger     
731         
732         TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), 
733                                                   fHijing->GetHINT1(27),
734                                                   fHijing->GetHINT1(28),
735                                                   fHijing->GetHINT1(29));
736         
737         TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), 
738                                                   fHijing->GetHINT1(37),
739                                                   fHijing->GetHINT1(38),
740                                                   fHijing->GetHINT1(39));
741         Double_t eta1      = jet1->Eta();
742         Double_t eta2      = jet2->Eta();
743         Double_t phi1      = jet1->Phi();
744         Double_t phi2      = jet2->Phi();
745 //    printf("\n Trigger: %f %f %f %f",
746 //         fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
747         if (
748             (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
749              phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
750             ||
751             (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
752              phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
753             ) 
754             triggered = kTRUE;
755     } else if (fTrigger == 2) {
756 //  Gamma Jet
757 //
758         Int_t np = fParticles->GetEntriesFast();
759         for (Int_t i = 0; i < np; i++) {
760             TParticle* part = (TParticle*) fParticles->At(i);
761             Int_t kf = part->GetPdgCode();
762             Int_t ks = part->GetStatusCode();
763             if (kf == 22 && ks == 40) {
764                 Float_t phi = part->Phi();
765                 Float_t eta = part->Eta();
766                 if  (eta < fEtaMaxJet && 
767                      eta > fEtaMinJet &&
768                      phi < fPhiMaxJet && 
769                      phi > fPhiMinJet) {
770                     triggered = 1;
771                     break;
772                 } // check phi,eta within limits
773             } // direct gamma ? 
774         } // particle loop
775     } // fTrigger == 2
776     return triggered;
777 }
778
779
780
781
782 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
783 {
784 // Assignment operator
785     return *this;
786 }
787
788 #ifndef WIN32
789 # define rluget_hijing rluget_hijing_
790 # define rluset_hijing rluset_hijing_
791 # define rlu_hijing    rlu_hijing_
792 # define type_of_call
793 #else
794 # define rluget_hijing RLUGET_HIJING
795 # define rluset_hijing RLUSET_HIJING
796 # define rlu_hijing    RLU_HIJING
797 # define type_of_call _stdcall
798 #endif
799
800
801 extern "C" {
802   void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
803   {printf("Dummy version of rluget_hijing reached\n");}
804
805   void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
806   {printf("Dummy version of rluset_hijing reached\n");}
807
808   Double_t type_of_call rlu_hijing(Int_t & /*idum*/) 
809   {
810       Float_t r;
811       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
812       return r;
813   }
814 }