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