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