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