fTrackit removed (fTrackIt exists instead)
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.cxx
CommitLineData
110410f9 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$
675e9664 18Revision 1.18 2000/12/06 17:46:30 morsch
19Avoid random numbers 1 and 0.
20
7498b3ad 21Revision 1.17 2000/12/04 11:22:03 morsch
22Init of sRandom as in 1.15
23
097ec7c1 24Revision 1.16 2000/12/02 11:41:39 morsch
25Use SetRandom() to initialize random number generator in constructor.
26
1ed25e11 27Revision 1.15 2000/11/30 20:29:02 morsch
28Initialise static variable sRandom in constructor: sRandom = fRandom;
29
ad1df918 30Revision 1.14 2000/11/30 07:12:50 alibrary
31Introducing new Rndm and QA classes
32
65fb704d 33Revision 1.13 2000/11/09 17:40:27 morsch
34Possibility to select/unselect spectator protons and neutrons.
35Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
36
2705454f 37Revision 1.12 2000/10/20 13:38:38 morsch
38Debug printouts commented.
39
1df2195f 40Revision 1.11 2000/10/20 13:22:26 morsch
41- skip particle type 92 (string)
42- Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
43 mesons.
44
608708ba 45Revision 1.10 2000/10/17 15:10:20 morsch
46Write first all the parent particles to the stack and then the final state particles.
47
e06552ba 48Revision 1.9 2000/10/17 13:38:59 morsch
49Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..) (FCA)
50
c007646a 51Revision 1.8 2000/10/17 12:46:31 morsch
52Protect EvaluateCrossSections() against division by zero.
53
73bc6a56 54Revision 1.7 2000/10/02 21:28:06 fca
55Removal of useless dependecies via forward declarations
56
94de3818 57Revision 1.6 2000/09/11 13:23:37 morsch
58Write last seed to file (fortran lun 50) and reed back from same lun using calls to
59luget_hijing and luset_hijing.
60
b78ff279 61Revision 1.5 2000/09/07 16:55:40 morsch
62fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
63
ba244bac 64Revision 1.4 2000/07/11 18:24:56 fca
65Coding convention corrections + few minor bug fixes
66
aee8290b 67Revision 1.3 2000/06/30 12:08:36 morsch
68In member data: char* replaced by TString, Init takes care of resizing the strings to
698 characters required by Hijing.
70
3b8bd63f 71Revision 1.2 2000/06/15 14:15:05 morsch
72Add possibility for heavy flavor selection: charm and beauty.
73
eb342d2d 74Revision 1.1 2000/06/09 20:47:27 morsch
75AliGenerator interface class to HIJING using THijing (test version)
76
110410f9 77*/
78
675e9664 79
80
81// Generator using HIJING as an external generator
82// The main HIJING options are accessable for the user through this interface.
83// Uses the THijing implementation of TGenerator.
84//
85// andreas.morsch@cern.ch
86
110410f9 87#include "AliGenHijing.h"
3b8bd63f 88#include "AliGenHijingEventHeader.h"
110410f9 89#include "AliRun.h"
90
91#include <TArrayI.h>
92#include <TParticle.h>
93#include <THijing.h>
94
ad1df918 95
110410f9 96 ClassImp(AliGenHijing)
97
98AliGenHijing::AliGenHijing()
99 :AliGenerator()
100{
101// Constructor
102}
103
104AliGenHijing::AliGenHijing(Int_t npart)
ad1df918 105 :AliGenerator(npart)
110410f9 106{
107// Default PbPb collisions at 5. 5 TeV
108//
109 SetEnergyCMS();
110 SetImpactParameterRange();
111 SetTarget();
112 SetProjectile();
113 fKeep=0;
114 fQuench=1;
115 fShadowing=1;
116 fTrigger=0;
117 fDecaysOff=1;
118 fEvaluate=0;
119 fSelectAll=0;
eb342d2d 120 fFlavor=0;
2705454f 121 fSpectators=1;
ad1df918 122//
1ed25e11 123// Set random number generator
097ec7c1 124 sRandom = fRandom;
110410f9 125}
126
127AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
128{
129// copy constructor
130}
131
132
133AliGenHijing::~AliGenHijing()
134{
135// Destructor
136}
137
138void AliGenHijing::Init()
139{
140// Initialisation
3b8bd63f 141 fFrame.Resize(8);
142 fTarget.Resize(8);
143 fProjectile.Resize(8);
144
110410f9 145 SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
146 fAProjectile, fZProjectile, fATarget, fZTarget,
147 fMinImpactParam, fMaxImpactParam));
148
149 fHijing=(THijing*) fgMCEvGen;
ba244bac 150
110410f9 151 fHijing->SetIHPR2(3, fTrigger);
152 fHijing->SetIHPR2(4, fQuench);
153 fHijing->SetIHPR2(6, fShadowing);
154 fHijing->SetIHPR2(12, fDecaysOff);
155 fHijing->SetIHPR2(21, fKeep);
b78ff279 156 fHijing->Rluset(50,0);
ba244bac 157 fHijing->Initialize();
b78ff279 158
159
110410f9 160//
161 if (fEvaluate) EvaluateCrossSections();
ad1df918 162//
163//
164// Initialize random generator
110410f9 165}
166
167void AliGenHijing::Generate()
168{
169// Generate one event
170
171 Float_t polar[3] = {0,0,0};
172 Float_t origin[3]= {0,0,0};
173 Float_t origin0[3]= {0,0,0};
174 Float_t p[3], random[6];
175 Float_t tof;
176
177 static TClonesArray *particles;
178// converts from mm/c to s
179 const Float_t kconv=0.001/2.999792458e8;
180//
181 Int_t nt=0;
182 Int_t jev=0;
183 Int_t j, kf, ks, imo;
eb342d2d 184 kf=0;
185
110410f9 186 if(!particles) particles=new TClonesArray("TParticle",10000);
187
188 fTrials=0;
189 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b 190 if(fVertexSmear==kPerEvent) {
65fb704d 191 Rndm(random,6);
110410f9 192 for (j=0;j<3;j++) {
193 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
194 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
195// fHijing->SetMSTP(151,0);
196 }
aee8290b 197 } else if (fVertexSmear==kPerTrack) {
110410f9 198// fHijing->SetMSTP(151,0);
199 for (j=0;j<3;j++) {
200// fHijing->SetPARP(151+j, fOsigma[j]*10.);
201 }
202 }
203 while(1)
204 {
205
206 fHijing->GenerateEvent();
207 fTrials++;
208 fHijing->ImportParticles(particles,"All");
209 Int_t np = particles->GetEntriesFast();
210 printf("\n **************************************************%d\n",np);
211 Int_t nc=0;
212 if (np == 0 ) continue;
213 Int_t i;
214 Int_t * newPos = new Int_t[np];
e06552ba 215
eb342d2d 216 for (i = 0; i<np; i++) *(newPos+i)=i;
e06552ba 217//
218// First write parent particles
219//
220
eb342d2d 221 for (i = 0; i<np; i++) {
110410f9 222 TParticle * iparticle = (TParticle *) particles->At(i);
e06552ba 223// Is this a parent particle ?
608708ba 224 if (Stable(iparticle)) continue;
e06552ba 225//
226 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
110410f9 227 Bool_t selected = kTRUE;
228 Bool_t hasSelectedDaughters = kFALSE;
e06552ba 229
230
eb342d2d 231 kf = iparticle->GetPdgCode();
608708ba 232 ks = iparticle->GetStatusCode();
233 if (kf == 92) continue;
234
eb342d2d 235 if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
e06552ba 236 hasSelectedDaughters = DaughtersSelection(iparticle, particles);
110410f9 237//
238// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
239//
110410f9 240 if (selected || hasSelectedDaughters) {
241 nc++;
110410f9 242 p[0]=iparticle->Px();
243 p[1]=iparticle->Py();
244 p[2]=iparticle->Pz();
245 origin[0]=origin0[0]+iparticle->Vx()/10;
246 origin[1]=origin0[1]+iparticle->Vy()/10;
247 origin[2]=origin0[2]+iparticle->Vz()/10;
248 tof=kconv*iparticle->T();
249 imo=-1;
250 if (hasMother) {
251 imo=iparticle->GetFirstMother();
608708ba 252 TParticle* mother= (TParticle *) particles->At(imo);
253 imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
110410f9 254 }
e06552ba 255// Put particle on the stack ...
1df2195f 256// printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
608708ba 257
e06552ba 258 gAlice->SetTrack(0,imo,kf,p,origin,polar,
65fb704d 259 tof,kPPrimary,nt);
e06552ba 260// ... and keep it there
261 gAlice->KeepTrack(nt);
262//
263 *(newPos+i)=nt;
264 } // selected
265 } // particle loop parents
266//
267// Now write the final state particles
268//
269
270 for (i = 0; i<np; i++) {
271 TParticle * iparticle = (TParticle *) particles->At(i);
e06552ba 272// Is this a final state particle ?
608708ba 273 if (!Stable(iparticle)) continue;
e06552ba 274//
275 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
276 Bool_t selected = kTRUE;
277 kf = iparticle->GetPdgCode();
2705454f 278 ks = iparticle->GetStatusCode();
279 if (!fSelectAll) {
280 selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
281 if (!fSpectators && selected) selected = (ks != 0 && ks != 10);
282 }
e06552ba 283//
284// Put particle on the stack if selected
285//
286 if (selected) {
287 nc++;
e06552ba 288 p[0]=iparticle->Px();
289 p[1]=iparticle->Py();
290 p[2]=iparticle->Pz();
291 origin[0]=origin0[0]+iparticle->Vx()/10;
292 origin[1]=origin0[1]+iparticle->Vy()/10;
293 origin[2]=origin0[2]+iparticle->Vz()/10;
294 tof=kconv*iparticle->T();
295 imo=-1;
608708ba 296
e06552ba 297 if (hasMother) {
298 imo=iparticle->GetFirstMother();
608708ba 299 TParticle* mother= (TParticle *) particles->At(imo);
300 imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
110410f9 301 }
e06552ba 302// Put particle on the stack
303 gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
65fb704d 304 tof,kPNoProcess,nt);
305// tof,"Secondary",nt);
608708ba 306
1df2195f 307// printf("\n set track final: %d %d %d",imo, kf, nt);
608708ba 308 gAlice->KeepTrack(nt);
110410f9 309 *(newPos+i)=nt;
310 } // selected
e06552ba 311 } // particle loop final state
312
110410f9 313 delete newPos;
e06552ba 314
110410f9 315 printf("\n I've put %i particles on the stack \n",nc);
316 if (nc > 0) {
317 jev+=nc;
318 if (jev >= fNpart || fNpart == -1) {
319 fKineBias=Float_t(fNpart)/Float_t(fTrials);
320 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
321 break;
322 }
323 }
324 } // event loop
b78ff279 325 fHijing->Rluget(50,-1);
110410f9 326}
327
328Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
329{
330// Perform kinematic selection
c007646a 331 Double_t px=particle->Px();
332 Double_t py=particle->Py();
333 Double_t pz=particle->Pz();
334 Double_t e=particle->Energy();
110410f9 335
336//
337// transverse momentum cut
c007646a 338 Double_t pt=TMath::Sqrt(px*px+py*py);
110410f9 339 if (pt > fPtMax || pt < fPtMin)
340 {
341// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
342 return kFALSE;
343 }
344//
345// momentum cut
c007646a 346 Double_t p=TMath::Sqrt(px*px+py*py+pz*pz);
110410f9 347 if (p > fPMax || p < fPMin)
348 {
349// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
350 return kFALSE;
351 }
352
353//
354// theta cut
c007646a 355 Double_t theta = Double_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
110410f9 356 if (theta > fThetaMax || theta < fThetaMin)
357 {
358
359// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
360 return kFALSE;
361 }
362
363//
364// rapidity cut
c007646a 365 Double_t y;
366 if(e<=pz) y = 99;
367 else if (e<=-pz) y = -99;
368 else y = 0.5*TMath::Log((e+pz)/(e-pz));
110410f9 369 if (y > fYMax || y < fYMin)
370 {
371// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
372 return kFALSE;
373 }
374
375//
376// phi cut
c007646a 377 Double_t phi=Double_t(TMath::ATan2(Double_t(py),Double_t(px)));
110410f9 378 if (phi > fPhiMax || phi < fPhiMin)
379 {
380// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
381 return kFALSE;
382 }
383
384 return kTRUE;
385}
386
387void AliGenHijing::KeepFullEvent()
388{
389 fKeep=1;
390}
391
392void AliGenHijing::EvaluateCrossSections()
393{
394// Glauber Calculation of geometrical x-section
395//
396 Float_t xTot=0.; // barn
397 Float_t xTotHard=0.; // barn
398 Float_t xPart=0.; // barn
399 Float_t xPartHard=0.; // barn
400 Float_t sigmaHard=0.1; // mbarn
401 Float_t bMin=0.;
402 Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
403 const Float_t kdib=0.2;
404 Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
405
406
407 printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
408 printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
409 Int_t i;
410 Float_t oldvalue=0.;
411
412 for (i=0; i<kMax; i++)
413 {
414 Float_t xb=bMin+i*kdib;
415 Float_t ov;
416 ov=fHijing->Profile(xb);
417 Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
418 Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
419 xTot+=gb;
420 xTotHard+=gbh;
421 if (xb > fMinImpactParam && xb < fMaxImpactParam)
422 {
423 xPart+=gb;
424 xPartHard+=gbh;
425 }
426
c007646a 427 if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
110410f9 428 oldvalue=xTot;
429 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
430 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
431 }
432 printf("\n Total cross section (barn): %f \n",xTot);
433 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
434 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
435 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
436}
437
438Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
439{
440//
441// Looks recursively if one of the daughters has been selected
442//
443// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
444 Int_t imin=-1;
445 Int_t imax=-1;
446 Int_t i;
447 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
448 Bool_t selected=kFALSE;
449 if (hasDaughters) {
450 imin=iparticle->GetFirstDaughter();
451 imax=iparticle->GetLastDaughter();
110410f9 452 for (i=imin; i<= imax; i++){
453 TParticle * jparticle = (TParticle *) particles->At(i);
eb342d2d 454 Int_t ip=jparticle->GetPdgCode();
608708ba 455 if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {
608708ba 456 selected=kTRUE; break;
457 }
110410f9 458 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
459 }
460 } else {
461 return kFALSE;
462 }
463
464 return selected;
465}
466
467
eb342d2d 468Bool_t AliGenHijing::SelectFlavor(Int_t pid)
469{
470// Select flavor of particle
471// 0: all
472// 4: charm and beauty
473// 5: beauty
474 if (fFlavor == 0) return kTRUE;
475
476 Int_t ifl=TMath::Abs(pid/100);
477 if (ifl > 10) ifl/=10;
608708ba 478 return (fFlavor == ifl);
479}
eb342d2d 480
608708ba 481Bool_t AliGenHijing::Stable(TParticle* particle)
482{
675e9664 483// Return true for a stable particle
484//
608708ba 485 Int_t kf = TMath::Abs(particle->GetPdgCode());
486
487 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
488
489 {
490 return kTRUE;
491 } else {
492 return kFALSE;
493 }
eb342d2d 494}
110410f9 495
3b8bd63f 496void AliGenHijing::MakeHeader()
497{
498// Builds the event header, to be called after each event
499 AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
500// header->SetDate(date);
501// header->SetRunNumber(run);
502// header->SetEventNumber(event);
503 header->SetNProduced(fHijing->GetNATT());
504 header->SetImpactParameter(fHijing->GetHINT1(19));
505 header->SetTotalEnergy(fHijing->GetEATT());
506 header->SetHardScatters(fHijing->GetJATT());
507 header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
508 header->SetCollisions(fHijing->GetN0(),
509 fHijing->GetN01(),
510 fHijing->GetN10(),
511 fHijing->GetN11());
512}
513
110410f9 514AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
515{
516// Assignment operator
517 return *this;
518}
608708ba 519
65fb704d 520#ifndef WIN32
521# define rluget_hijing rluget_hijing_
522# define rluset_hijing rluset_hijing_
523# define rlu_hijing rlu_hijing_
524# define type_of_call
525#else
526# define rluget_hijing RLUGET_HIJING
527# define rluset_hijing RLUSET_HIJING
528# define rlu_hijing RLU_HIJING
529# define type_of_call _stdcall
530#endif
531
532
533extern "C" {
534 void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
535 {printf("Dummy version of rluget_hijing reached\n");}
536
537 void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
538 {printf("Dummy version of rluset_hijing reached\n");}
539
540 Double_t type_of_call rlu_hijing(Int_t & /*idum*/)
7498b3ad 541 {
542 Float_t r;
543 do r=sRandom->Rndm(); while(0 >= r || r >= 1);
544 return r;
545 }
65fb704d 546}