]>
Commit | Line | Data |
---|---|---|
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 | 18 | Revision 1.13 2000/11/09 17:40:27 morsch |
19 | Possibility to select/unselect spectator protons and neutrons. | |
20 | Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano) | |
21 | ||
2705454f | 22 | Revision 1.12 2000/10/20 13:38:38 morsch |
23 | Debug printouts commented. | |
24 | ||
1df2195f | 25 | Revision 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 | 30 | Revision 1.10 2000/10/17 15:10:20 morsch |
31 | Write first all the parent particles to the stack and then the final state particles. | |
32 | ||
e06552ba | 33 | Revision 1.9 2000/10/17 13:38:59 morsch |
34 | Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..) (FCA) | |
35 | ||
c007646a | 36 | Revision 1.8 2000/10/17 12:46:31 morsch |
37 | Protect EvaluateCrossSections() against division by zero. | |
38 | ||
73bc6a56 | 39 | Revision 1.7 2000/10/02 21:28:06 fca |
40 | Removal of useless dependecies via forward declarations | |
41 | ||
94de3818 | 42 | Revision 1.6 2000/09/11 13:23:37 morsch |
43 | Write last seed to file (fortran lun 50) and reed back from same lun using calls to | |
44 | luget_hijing and luset_hijing. | |
45 | ||
b78ff279 | 46 | Revision 1.5 2000/09/07 16:55:40 morsch |
47 | fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko) | |
48 | ||
ba244bac | 49 | Revision 1.4 2000/07/11 18:24:56 fca |
50 | Coding convention corrections + few minor bug fixes | |
51 | ||
aee8290b | 52 | Revision 1.3 2000/06/30 12:08:36 morsch |
53 | In member data: char* replaced by TString, Init takes care of resizing the strings to | |
54 | 8 characters required by Hijing. | |
55 | ||
3b8bd63f | 56 | Revision 1.2 2000/06/15 14:15:05 morsch |
57 | Add possibility for heavy flavor selection: charm and beauty. | |
58 | ||
eb342d2d | 59 | Revision 1.1 2000/06/09 20:47:27 morsch |
60 | AliGenerator 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 | ||
74 | AliGenHijing::AliGenHijing() | |
75 | :AliGenerator() | |
76 | { | |
77 | // Constructor | |
78 | } | |
79 | ||
80 | AliGenHijing::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 | ||
100 | AliGenHijing::AliGenHijing(const AliGenHijing & Hijing) | |
101 | { | |
102 | // copy constructor | |
103 | } | |
104 | ||
105 | ||
106 | AliGenHijing::~AliGenHijing() | |
107 | { | |
108 | // Destructor | |
109 | } | |
110 | ||
111 | void 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 | ||
137 | void 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 | ||
298 | Bool_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 | ||
357 | void AliGenHijing::KeepFullEvent() | |
358 | { | |
359 | fKeep=1; | |
360 | } | |
361 | ||
362 | void 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 | ||
408 | Bool_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 | 438 | Bool_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 | 451 | Bool_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 | 464 | void 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 | 482 | AliGenHijing& 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 | ||
501 | extern "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 |