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