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