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