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