]>
Commit | Line | Data |
---|---|---|
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$ | |
18 | Revision 1.16 2000/12/02 11:41:39 morsch | |
19 | Use SetRandom() to initialize random number generator in constructor. | |
20 | ||
21 | Revision 1.15 2000/11/30 20:29:02 morsch | |
22 | Initialise static variable sRandom in constructor: sRandom = fRandom; | |
23 | ||
24 | Revision 1.14 2000/11/30 07:12:50 alibrary | |
25 | Introducing new Rndm and QA classes | |
26 | ||
27 | Revision 1.13 2000/11/09 17:40:27 morsch | |
28 | Possibility to select/unselect spectator protons and neutrons. | |
29 | Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano) | |
30 | ||
31 | Revision 1.12 2000/10/20 13:38:38 morsch | |
32 | Debug printouts commented. | |
33 | ||
34 | Revision 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 | ||
39 | Revision 1.10 2000/10/17 15:10:20 morsch | |
40 | Write first all the parent particles to the stack and then the final state particles. | |
41 | ||
42 | Revision 1.9 2000/10/17 13:38:59 morsch | |
43 | Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..) (FCA) | |
44 | ||
45 | Revision 1.8 2000/10/17 12:46:31 morsch | |
46 | Protect EvaluateCrossSections() against division by zero. | |
47 | ||
48 | Revision 1.7 2000/10/02 21:28:06 fca | |
49 | Removal of useless dependecies via forward declarations | |
50 | ||
51 | Revision 1.6 2000/09/11 13:23:37 morsch | |
52 | Write last seed to file (fortran lun 50) and reed back from same lun using calls to | |
53 | luget_hijing and luset_hijing. | |
54 | ||
55 | Revision 1.5 2000/09/07 16:55:40 morsch | |
56 | fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko) | |
57 | ||
58 | Revision 1.4 2000/07/11 18:24:56 fca | |
59 | Coding convention corrections + few minor bug fixes | |
60 | ||
61 | Revision 1.3 2000/06/30 12:08:36 morsch | |
62 | In member data: char* replaced by TString, Init takes care of resizing the strings to | |
63 | 8 characters required by Hijing. | |
64 | ||
65 | Revision 1.2 2000/06/15 14:15:05 morsch | |
66 | Add possibility for heavy flavor selection: charm and beauty. | |
67 | ||
68 | Revision 1.1 2000/06/09 20:47:27 morsch | |
69 | AliGenerator interface class to HIJING using THijing (test version) | |
70 | ||
71 | */ | |
72 | ||
73 | #include "AliGenHijing.h" | |
74 | #include "AliGenHijingEventHeader.h" | |
75 | #include "AliRun.h" | |
76 | ||
77 | #include <TArrayI.h> | |
78 | #include <TParticle.h> | |
79 | #include <THijing.h> | |
80 | ||
81 | ||
82 | ClassImp(AliGenHijing) | |
83 | ||
84 | AliGenHijing::AliGenHijing() | |
85 | :AliGenerator() | |
86 | { | |
87 | // Constructor | |
88 | } | |
89 | ||
90 | AliGenHijing::AliGenHijing(Int_t npart) | |
91 | :AliGenerator(npart) | |
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; | |
106 | fFlavor=0; | |
107 | fSpectators=1; | |
108 | // | |
109 | // Set random number generator | |
110 | sRandom = fRandom; | |
111 | } | |
112 | ||
113 | AliGenHijing::AliGenHijing(const AliGenHijing & Hijing) | |
114 | { | |
115 | // copy constructor | |
116 | } | |
117 | ||
118 | ||
119 | AliGenHijing::~AliGenHijing() | |
120 | { | |
121 | // Destructor | |
122 | } | |
123 | ||
124 | void AliGenHijing::Init() | |
125 | { | |
126 | // Initialisation | |
127 | fFrame.Resize(8); | |
128 | fTarget.Resize(8); | |
129 | fProjectile.Resize(8); | |
130 | ||
131 | SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, | |
132 | fAProjectile, fZProjectile, fATarget, fZTarget, | |
133 | fMinImpactParam, fMaxImpactParam)); | |
134 | ||
135 | fHijing=(THijing*) fgMCEvGen; | |
136 | ||
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); | |
142 | fHijing->Rluset(50,0); | |
143 | fHijing->Initialize(); | |
144 | ||
145 | ||
146 | // | |
147 | if (fEvaluate) EvaluateCrossSections(); | |
148 | // | |
149 | // | |
150 | // Initialize random generator | |
151 | } | |
152 | ||
153 | void 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; | |
170 | kf=0; | |
171 | ||
172 | if(!particles) particles=new TClonesArray("TParticle",10000); | |
173 | ||
174 | fTrials=0; | |
175 | for (j=0;j<3;j++) origin0[j]=fOrigin[j]; | |
176 | if(fVertexSmear==kPerEvent) { | |
177 | Rndm(random,6); | |
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 | } | |
183 | } else if (fVertexSmear==kPerTrack) { | |
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]; | |
201 | ||
202 | for (i = 0; i<np; i++) *(newPos+i)=i; | |
203 | // | |
204 | // First write parent particles | |
205 | // | |
206 | ||
207 | for (i = 0; i<np; i++) { | |
208 | TParticle * iparticle = (TParticle *) particles->At(i); | |
209 | // Is this a parent particle ? | |
210 | if (Stable(iparticle)) continue; | |
211 | // | |
212 | Bool_t hasMother = (iparticle->GetFirstMother() >=0); | |
213 | Bool_t selected = kTRUE; | |
214 | Bool_t hasSelectedDaughters = kFALSE; | |
215 | ||
216 | ||
217 | kf = iparticle->GetPdgCode(); | |
218 | ks = iparticle->GetStatusCode(); | |
219 | if (kf == 92) continue; | |
220 | ||
221 | if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf); | |
222 | hasSelectedDaughters = DaughtersSelection(iparticle, particles); | |
223 | // | |
224 | // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle | |
225 | // | |
226 | if (selected || hasSelectedDaughters) { | |
227 | nc++; | |
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(); | |
238 | TParticle* mother= (TParticle *) particles->At(imo); | |
239 | imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1; | |
240 | } | |
241 | // Put particle on the stack ... | |
242 | // printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters); | |
243 | ||
244 | gAlice->SetTrack(0,imo,kf,p,origin,polar, | |
245 | tof,kPPrimary,nt); | |
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); | |
258 | // Is this a final state particle ? | |
259 | if (!Stable(iparticle)) continue; | |
260 | // | |
261 | Bool_t hasMother = (iparticle->GetFirstMother() >=0); | |
262 | Bool_t selected = kTRUE; | |
263 | kf = iparticle->GetPdgCode(); | |
264 | ks = iparticle->GetStatusCode(); | |
265 | if (!fSelectAll) { | |
266 | selected = KinematicSelection(iparticle)&&SelectFlavor(kf); | |
267 | if (!fSpectators && selected) selected = (ks != 0 && ks != 10); | |
268 | } | |
269 | // | |
270 | // Put particle on the stack if selected | |
271 | // | |
272 | if (selected) { | |
273 | nc++; | |
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; | |
282 | ||
283 | if (hasMother) { | |
284 | imo=iparticle->GetFirstMother(); | |
285 | TParticle* mother= (TParticle *) particles->At(imo); | |
286 | imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1; | |
287 | } | |
288 | // Put particle on the stack | |
289 | gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar, | |
290 | tof,kPNoProcess,nt); | |
291 | // tof,"Secondary",nt); | |
292 | ||
293 | // printf("\n set track final: %d %d %d",imo, kf, nt); | |
294 | gAlice->KeepTrack(nt); | |
295 | *(newPos+i)=nt; | |
296 | } // selected | |
297 | } // particle loop final state | |
298 | ||
299 | delete newPos; | |
300 | ||
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 | |
311 | fHijing->Rluget(50,-1); | |
312 | } | |
313 | ||
314 | Bool_t AliGenHijing::KinematicSelection(TParticle *particle) | |
315 | { | |
316 | // Perform kinematic selection | |
317 | Double_t px=particle->Px(); | |
318 | Double_t py=particle->Py(); | |
319 | Double_t pz=particle->Pz(); | |
320 | Double_t e=particle->Energy(); | |
321 | ||
322 | // | |
323 | // transverse momentum cut | |
324 | Double_t pt=TMath::Sqrt(px*px+py*py); | |
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 | |
332 | Double_t p=TMath::Sqrt(px*px+py*py+pz*pz); | |
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 | |
341 | Double_t theta = Double_t(TMath::ATan2(Double_t(pt),Double_t(pz))); | |
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 | |
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)); | |
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 | |
363 | Double_t phi=Double_t(TMath::ATan2(Double_t(py),Double_t(px))); | |
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 | ||
373 | void AliGenHijing::KeepFullEvent() | |
374 | { | |
375 | fKeep=1; | |
376 | } | |
377 | ||
378 | void 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 | ||
413 | if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break; | |
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 | ||
424 | Bool_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(); | |
438 | for (i=imin; i<= imax; i++){ | |
439 | TParticle * jparticle = (TParticle *) particles->At(i); | |
440 | Int_t ip=jparticle->GetPdgCode(); | |
441 | if (KinematicSelection(jparticle)&&SelectFlavor(ip)) { | |
442 | selected=kTRUE; break; | |
443 | } | |
444 | if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; } | |
445 | } | |
446 | } else { | |
447 | return kFALSE; | |
448 | } | |
449 | ||
450 | return selected; | |
451 | } | |
452 | ||
453 | ||
454 | Bool_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; | |
464 | return (fFlavor == ifl); | |
465 | } | |
466 | ||
467 | Bool_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 | } | |
478 | } | |
479 | ||
480 | void 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 | ||
498 | AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs) | |
499 | { | |
500 | // Assignment operator | |
501 | return *this; | |
502 | } | |
503 | ||
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 | ||
517 | extern "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 | } |