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