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