Rapidity shift calculated in AliGenMC::Init()
[u/mrichter/AliRoot.git] / THijing / AliGenHijing.cxx
CommitLineData
36b81802 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$
71ea527c 18Revision 1.2 2003/04/04 08:14:02 morsch
19Boost() moved to AliGenMC
20
0c85382f 21Revision 1.1 2003/03/15 14:45:57 morsch
22Classes imported from EVGEN
23
36b81802 24Revision 1.47 2003/01/14 10:50:18 alibrary
25Cleanup of STEER coding conventions
26
27Revision 1.46 2003/01/07 14:12:33 morsch
28Provides collision geometry.
29
30Revision 1.45 2002/12/16 09:44:49 morsch
31Default for fRadiation is 3.
32
33Revision 1.44 2002/10/14 14:55:35 hristov
34Merging the VirtualMC branch to the main development branch (HEAD)
35
36Revision 1.42.4.1 2002/08/28 15:06:50 alibrary
37Updating to v3-09-01
38
39Revision 1.43 2002/08/09 12:09:52 morsch
40Direct gamma trigger correctly included.
41
42Revision 1.42 2002/03/12 11:07:08 morsch
43Add status code of particle to SetTrack call.
44
45Revision 1.41 2002/03/05 11:25:33 morsch
46- New quenching options
47- Correction in CheckTrigger()
48
49Revision 1.40 2002/02/12 11:05:53 morsch
50Get daughter indices right.
51
52Revision 1.39 2002/02/12 09:16:39 morsch
53Correction in SelectFlavor()
54
55Revision 1.38 2002/02/12 08:53:21 morsch
56SetNoGammas can be used to inhibit writing of gammas and pi0.
57
58Revision 1.37 2002/02/08 16:50:50 morsch
59Add name and title in constructor.
60
61Revision 1.36 2002/01/31 20:17:55 morsch
62Allow for triggered jets with simplified topology: Exact pT, back-to-back
63
64Revision 1.35 2001/12/13 07:56:25 hristov
65Set pointers to zero in the default constructor
66
67Revision 1.34 2001/12/11 16:55:42 morsch
68Correct initialization for jet phi-range.
69
70Revision 1.33 2001/12/05 10:18:51 morsch
71Possibility of kinematic biasing of jet phi range. (J. Klay)
72
73Revision 1.32 2001/11/28 13:51:11 morsch
74Introduce kinematic biasing (etamin, etamax) of jet trigger. Bookkeeping
75(number of trials) done in AliGenHijingEventHeader.
76
77Revision 1.31 2001/11/06 12:30:34 morsch
78Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric collision systems.
79
80Revision 1.30 2001/10/21 18:35:56 hristov
81Several pointers were set to zero in the default constructors to avoid memory management problems
82
83Revision 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
87Revision 1.28 2001/10/08 11:55:25 morsch
88Store 4-momenta of trigegred jets in event header.
89Possibility to switch of initial and final state radiation.
90
91Revision 1.27 2001/10/08 07:13:14 morsch
92Add setter for minimum transverse momentum of triggered jet.
93
94Revision 1.26 2001/10/04 08:12:24 morsch
95Redefinition of stable condition.
96
97Revision 1.25 2001/07/27 17:09:36 morsch
98Use local SetTrack, KeepTrack and SetHighWaterMark methods
99to delegate either to local stack or to stack owned by AliRun.
100(Piotr Skowronski, A.M.)
101
102Revision 1.24 2001/07/20 09:34:56 morsch
103Count the number of spectator neutrons and protons and add information
104to the event header. (Chiara Oppedisano)
105
106Revision 1.23 2001/07/13 17:30:22 morsch
107Derive from AliGenMC.
108
109Revision 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
113Revision 1.21 2001/02/28 17:35:24 morsch
114Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
115
116Revision 1.20 2001/02/14 15:50:40 hristov
117The last particle in event marked using SetHighWaterMark
118
119Revision 1.19 2000/12/21 16:24:06 morsch
120Coding convention clean-up
121
122Revision 1.18 2000/12/06 17:46:30 morsch
123Avoid random numbers 1 and 0.
124
125Revision 1.17 2000/12/04 11:22:03 morsch
126Init of sRandom as in 1.15
127
128Revision 1.16 2000/12/02 11:41:39 morsch
129Use SetRandom() to initialize random number generator in constructor.
130
131Revision 1.15 2000/11/30 20:29:02 morsch
132Initialise static variable sRandom in constructor: sRandom = fRandom;
133
134Revision 1.14 2000/11/30 07:12:50 alibrary
135Introducing new Rndm and QA classes
136
137Revision 1.13 2000/11/09 17:40:27 morsch
138Possibility to select/unselect spectator protons and neutrons.
139Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
140
141Revision 1.12 2000/10/20 13:38:38 morsch
142Debug printouts commented.
143
144Revision 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
149Revision 1.10 2000/10/17 15:10:20 morsch
150Write first all the parent particles to the stack and then the final state particles.
151
152Revision 1.9 2000/10/17 13:38:59 morsch
153Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..) (FCA)
154
155Revision 1.8 2000/10/17 12:46:31 morsch
156Protect EvaluateCrossSections() against division by zero.
157
158Revision 1.7 2000/10/02 21:28:06 fca
159Removal of useless dependecies via forward declarations
160
161Revision 1.6 2000/09/11 13:23:37 morsch
162Write last seed to file (fortran lun 50) and reed back from same lun using calls to
163luget_hijing and luset_hijing.
164
165Revision 1.5 2000/09/07 16:55:40 morsch
166fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
167
168Revision 1.4 2000/07/11 18:24:56 fca
169Coding convention corrections + few minor bug fixes
170
171Revision 1.3 2000/06/30 12:08:36 morsch
172In member data: char* replaced by TString, Init takes care of resizing the strings to
1738 characters required by Hijing.
174
175Revision 1.2 2000/06/15 14:15:05 morsch
176Add possibility for heavy flavor selection: charm and beauty.
177
178Revision 1.1 2000/06/09 20:47:27 morsch
179AliGenerator 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
205AliGenHijing::AliGenHijing()
206 :AliGenMC()
207{
208// Constructor
209 fParticles = 0;
210 fHijing = 0;
211 fDsigmaDb = 0;
212 fDnDb = 0;
213}
214
215AliGenHijing::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();
36b81802 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
255AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
256{
257// copy constructor
258}
259
260
261AliGenHijing::~AliGenHijing()
262{
263// Destructor
264 if ( fDsigmaDb) delete fDsigmaDb;
265 if ( fDnDb) delete fDnDb;
266 delete fParticles;
267}
268
269void 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
331void 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 }
71ea527c 383 if (fLHC) Boost();
36b81802 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
519void AliGenHijing::KeepFullEvent()
520{
521 fKeep=1;
522}
523
524void 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
590Bool_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
619Bool_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
642Bool_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
36b81802 656
657void 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
703Bool_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
763AliGenHijing& 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
782extern "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}