]>
Commit | Line | Data |
---|---|---|
ec852657 | 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 | // | |
18 | // Glauber MC implementation | |
19 | // | |
20 | // origin: PHOBOS experiment | |
21 | // alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch | |
106d1ca1 | 22 | // update: You Zhou, Nikhef, yzhou@nikhef.nl |
ec852657 | 23 | //////////////////////////////////////////////////////////////////////////////// |
24 | ||
25 | #include <Riostream.h> | |
26 | #include <TMath.h> | |
27 | #include <TEllipse.h> | |
28 | #include <TRandom.h> | |
29 | #include <TNamed.h> | |
30 | #include <TObjArray.h> | |
31 | #include <TNtuple.h> | |
32 | #include <TFile.h> | |
bb442091 | 33 | #include <TTree.h> |
34 | #include <TF1.h> | |
35 | #include <TH1F.h> | |
36 | #include <TArray.h> | |
ec852657 | 37 | |
38 | #include "AliGlauberNucleon.h" | |
39 | #include "AliGlauberNucleus.h" | |
40 | #include "AliGlauberMC.h" | |
41 | ||
42 | ClassImp(AliGlauberMC) | |
43 | ||
44 | //______________________________________________________________________________ | |
45 | AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) : | |
566fc3d0 | 46 | TNamed(), |
bb442091 | 47 | fANucleus(NA), |
48 | fBNucleus(NB), | |
566fc3d0 | 49 | fXSect(xsect), |
bb442091 | 50 | fNucleonsA(0), |
51 | fNucleonsB(0), | |
566fc3d0 | 52 | fAN(0), |
53 | fBN(0), | |
bb442091 | 54 | fnt(0), |
55 | fMeanX2(0), | |
56 | fMeanY2(0), | |
57 | fMeanXY(0), | |
58 | fMeanXParts(0), | |
59 | fMeanYParts(0), | |
60 | fMeanXColl(0), | |
61 | fMeanYColl(0), | |
62 | fMeanX2Coll(0), | |
63 | fMeanY2Coll(0), | |
64 | fMeanXYColl(0), | |
65 | fMeanXSystem(0), | |
66 | fMeanYSystem(0), | |
43215add | 67 | fMeanXA(0), |
68 | fMeanYA(0), | |
69 | fMeanXB(0), | |
70 | fMeanYB(0), | |
71 | fBMC(0), | |
bb442091 | 72 | fEvents(0), |
73 | fTotalEvents(0), | |
566fc3d0 | 74 | fBMin(0.), |
75 | fBMax(20.), | |
bb442091 | 76 | fMaxNpartFound(0), |
77 | fNpart(0), | |
78 | fNcoll(0), | |
106d1ca1 | 79 | fMeanr2(0), |
80 | fMeanr3(0), | |
81 | fMeanr4(0), | |
82 | fMeanr5(0), | |
83 | fMeanr2Cos2Phi(0), | |
84 | fMeanr2Sin2Phi(0), | |
85 | fMeanr2Cos3Phi(0), | |
86 | fMeanr2Sin3Phi(0), | |
87 | fMeanr2Cos4Phi(0), | |
88 | fMeanr2Sin4Phi(0), | |
89 | fMeanr2Cos5Phi(0), | |
90 | fMeanr2Sin5Phi(0), | |
91 | fMeanr3Cos3Phi(0), | |
92 | fMeanr3Sin3Phi(0), | |
93 | fMeanr4Cos4Phi(0), | |
94 | fMeanr4Sin4Phi(0), | |
95 | fMeanr5Cos5Phi(0), | |
96 | fMeanr5Sin5Phi(0), | |
566fc3d0 | 97 | fSx2(0.), |
98 | fSy2(0.), | |
99 | fSxy(0.), | |
100 | fSx2Coll(0.), | |
101 | fSy2Coll(0.), | |
102 | fSxyColl(0.), | |
bb442091 | 103 | fX(0.13), |
14244e14 | 104 | fNpp(8.), |
fee95d0f | 105 | fDoPartProd(kFALSE) |
ec852657 | 106 | { |
43215add | 107 | //ctor |
bb442091 | 108 | fdNdEtaParam[0] = 8.0; |
109 | fdNdEtaParam[1] = 0.13; | |
110 | fdNdEtaGBWParam[0] = 0.79; | |
111 | fdNdEtaGBWParam[1] = 0.288; | |
112 | fdNdEtaGBWParam[2] = 30.25; | |
113 | fdNdEtaNBDParam[0] = 3.0; | |
114 | fdNdEtaNBDParam[1] = 4.0; | |
115 | fdNdEtaNBDParam[2] = 0.13; | |
566fc3d0 | 116 | fdNdEtaTwoNBDParam[0] = 3.0; |
117 | fdNdEtaTwoNBDParam[1] = 4.0; | |
118 | fdNdEtaTwoNBDParam[2] = 2.0; | |
119 | fdNdEtaTwoNBDParam[3] = 11.0; | |
120 | fdNdEtaTwoNBDParam[4] = 0.4; | |
bb442091 | 121 | fdNdEtaTwoNBDParam[5] = 0.13; |
122 | ||
566fc3d0 | 123 | SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName())); |
124 | SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName())); | |
125 | } | |
126 | ||
127 | //______________________________________________________________________________ | |
128 | AliGlauberMC::~AliGlauberMC() | |
129 | { | |
130 | //dtor | |
131 | delete fnt; | |
132 | } | |
133 | ||
134 | //______________________________________________________________________________ | |
135 | AliGlauberMC::AliGlauberMC(const AliGlauberMC& in): | |
136 | TNamed(in), | |
137 | fANucleus(in.fANucleus), | |
138 | fBNucleus(in.fBNucleus), | |
139 | fXSect(in.fXSect), | |
140 | fNucleonsA(in.fNucleonsA), | |
141 | fNucleonsB(in.fNucleonsB), | |
142 | fAN(in.fAN), | |
143 | fBN(in.fBN), | |
144 | fnt(in.fnt), | |
145 | fMeanX2(in.fMeanX2), | |
146 | fMeanY2(in.fMeanY2), | |
147 | fMeanXY(in.fMeanXY), | |
148 | fMeanXParts(in.fMeanXParts), | |
149 | fMeanYParts(in.fMeanYParts), | |
150 | fMeanXColl(in.fMeanXColl), | |
151 | fMeanYColl(in.fMeanYColl), | |
152 | fMeanX2Coll(in.fMeanX2Coll), | |
153 | fMeanY2Coll(in.fMeanY2Coll), | |
154 | fMeanXYColl(in.fMeanXYColl), | |
155 | fMeanXSystem(in.fMeanXSystem), | |
156 | fMeanYSystem(in.fMeanYSystem), | |
43215add | 157 | fMeanXA(in.fMeanXA), |
158 | fMeanYA(in.fMeanYA), | |
159 | fMeanXB(in.fMeanXB), | |
160 | fMeanYB(in.fMeanYB), | |
161 | fBMC(in.fBMC), | |
566fc3d0 | 162 | fEvents(in.fEvents), |
163 | fTotalEvents(in.fTotalEvents), | |
164 | fBMin(in.fBMin), | |
165 | fBMax(in.fBMax), | |
166 | fMaxNpartFound(in.fMaxNpartFound), | |
167 | fNpart(in.fNpart), | |
168 | fNcoll(in.fNcoll), | |
106d1ca1 | 169 | fMeanr2(in.fMeanr2), |
170 | fMeanr3(in.fMeanr3), | |
171 | fMeanr4(in.fMeanr4), | |
172 | fMeanr5(in.fMeanr5), | |
173 | fMeanr2Cos2Phi(in.fMeanr2Cos2Phi), | |
174 | fMeanr2Sin2Phi(in.fMeanr2Sin2Phi), | |
175 | fMeanr2Cos3Phi(in.fMeanr2Cos3Phi), | |
176 | fMeanr2Sin3Phi(in.fMeanr2Sin3Phi), | |
177 | fMeanr2Cos4Phi(in.fMeanr2Cos4Phi), | |
178 | fMeanr2Sin4Phi(in.fMeanr2Sin4Phi), | |
179 | fMeanr2Cos5Phi(in.fMeanr2Cos5Phi), | |
180 | fMeanr2Sin5Phi(in.fMeanr2Sin5Phi), | |
181 | fMeanr3Cos3Phi(in.fMeanr3Cos3Phi), | |
182 | fMeanr3Sin3Phi(in.fMeanr3Sin3Phi), | |
183 | fMeanr4Cos4Phi(in.fMeanr4Cos4Phi), | |
184 | fMeanr4Sin4Phi(in.fMeanr4Sin4Phi), | |
185 | fMeanr5Cos5Phi(in.fMeanr5Cos5Phi), | |
186 | fMeanr5Sin5Phi(in.fMeanr5Sin5Phi), | |
566fc3d0 | 187 | fSx2(in.fSx2), |
188 | fSy2(in.fSy2), | |
189 | fSxy(in.fSxy), | |
190 | fSx2Coll(in.fSx2Coll), | |
191 | fSy2Coll(in.fSy2Coll), | |
192 | fSxyColl(in.fSxyColl), | |
193 | fX(in.fX), | |
97d98e48 | 194 | fNpp(in.fNpp), |
fee95d0f | 195 | fDoPartProd(kFALSE) |
566fc3d0 | 196 | { |
197 | //copy ctor | |
198 | memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t)); | |
199 | memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t)); | |
200 | memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t)); | |
201 | memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t)); | |
202 | } | |
203 | ||
204 | //______________________________________________________________________________ | |
205 | AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in) | |
206 | { | |
207 | //assignment | |
208 | fANucleus=in.fANucleus; | |
209 | fBNucleus=in.fBNucleus; | |
210 | fXSect=in.fXSect; | |
211 | fNucleonsA=in.fNucleonsA; | |
212 | fNucleonsB=in.fNucleonsB; | |
213 | fAN=in.fAN; | |
214 | fBN=in.fBN; | |
215 | fnt=in.fnt; | |
216 | fMeanX2=in.fMeanX2; | |
217 | fMeanY2=in.fMeanY2; | |
106d1ca1 | 218 | fMeanXY=in.fMeanXY; |
219 | fMeanr2=in.fMeanr2; | |
220 | fMeanr3=in.fMeanr3; | |
221 | fMeanr4=in.fMeanr4; | |
222 | fMeanr5=in.fMeanr5; | |
566fc3d0 | 223 | fMeanXParts=in.fMeanXParts; |
224 | fMeanYParts=in.fMeanYParts; | |
225 | fMeanXColl=in.fMeanXColl; | |
226 | fMeanYColl=in.fMeanYColl; | |
227 | fMeanX2Coll=in.fMeanX2Coll; | |
228 | fMeanY2Coll=in.fMeanY2Coll; | |
229 | fMeanXYColl=in.fMeanXYColl; | |
230 | fMeanXSystem=in.fMeanXSystem; | |
231 | fMeanYSystem=in.fMeanYSystem; | |
43215add | 232 | fMeanXA=in.fMeanXA; |
233 | fMeanYA=in.fMeanYA; | |
234 | fMeanXB=in.fMeanXB; | |
235 | fMeanYB=in.fMeanYB; | |
236 | fBMC=in.fBMC; | |
566fc3d0 | 237 | fEvents=in.fEvents; |
238 | fTotalEvents=in.fTotalEvents; | |
239 | fBMin=in.fBMin; | |
240 | fBMax=in.fBMax; | |
241 | memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t)); | |
242 | memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t)); | |
243 | memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t)); | |
244 | memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t)); | |
245 | fMaxNpartFound=in.fMaxNpartFound; | |
246 | fNpart=in.fNpart; | |
106d1ca1 | 247 | fNcoll=in.fNcoll; |
248 | fMeanr2Cos2Phi=in.fMeanr2Cos2Phi; | |
249 | fMeanr2Sin2Phi=in.fMeanr2Sin2Phi; | |
250 | fMeanr2Cos3Phi=in.fMeanr2Cos3Phi; | |
251 | fMeanr2Sin3Phi=in.fMeanr2Sin3Phi; | |
252 | fMeanr2Cos4Phi=in.fMeanr2Cos4Phi; | |
253 | fMeanr2Sin4Phi=in.fMeanr2Sin4Phi; | |
254 | fMeanr2Cos5Phi=in.fMeanr2Cos5Phi; | |
255 | fMeanr2Sin5Phi=in.fMeanr2Sin5Phi; | |
256 | fMeanr3Cos3Phi=in.fMeanr3Cos3Phi; | |
257 | fMeanr3Sin3Phi=in.fMeanr3Sin3Phi; | |
258 | fMeanr4Cos4Phi=in.fMeanr4Cos4Phi; | |
259 | fMeanr4Sin4Phi=in.fMeanr4Sin4Phi; | |
260 | fMeanr5Cos5Phi=in.fMeanr5Cos5Phi; | |
261 | fMeanr5Sin5Phi=in.fMeanr5Sin5Phi; | |
566fc3d0 | 262 | fSx2=in.fSx2; |
263 | fSy2=in.fSy2; | |
264 | fSxy=in.fSxy; | |
265 | fSx2Coll=in.fSx2Coll; | |
266 | fSy2Coll=in.fSy2Coll; | |
267 | fSxyColl=in.fSxyColl; | |
268 | fX=in.fX; | |
269 | fNpp=in.fNpp; | |
270 | return *this; | |
ec852657 | 271 | } |
272 | ||
273 | //______________________________________________________________________________ | |
274 | Bool_t AliGlauberMC::CalcEvent(Double_t bgen) | |
275 | { | |
bb442091 | 276 | // prepare event |
277 | fANucleus.ThrowNucleons(-bgen/2.); | |
278 | fNucleonsA = fANucleus.GetNucleons(); | |
279 | fAN = fANucleus.GetN(); | |
280 | for (Int_t i = 0; i<fAN; i++) | |
106d1ca1 | 281 | { |
282 | AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i)); | |
283 | nucleonA->SetInNucleusA(); | |
284 | } | |
bb442091 | 285 | fBNucleus.ThrowNucleons(bgen/2.); |
286 | fNucleonsB = fBNucleus.GetNucleons(); | |
287 | fBN = fBNucleus.GetN(); | |
288 | for (Int_t i = 0; i<fBN; i++) | |
106d1ca1 | 289 | { |
290 | AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i)); | |
291 | nucleonB->SetInNucleusB(); | |
292 | } | |
293 | ||
bb442091 | 294 | // "ball" diameter = distance at which two balls interact |
295 | Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2 | |
106d1ca1 | 296 | |
bb442091 | 297 | // for each of the A nucleons in nucleus B |
298 | for (Int_t i = 0; i<fBN; i++) | |
bb442091 | 299 | { |
106d1ca1 | 300 | AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i)); |
301 | for (Int_t j = 0 ; j < fAN ; j++) | |
302 | { | |
303 | AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j)); | |
304 | Double_t dx = nucleonB->GetX()-nucleonA->GetX(); | |
305 | Double_t dy = nucleonB->GetY()-nucleonA->GetY(); | |
306 | Double_t dij = dx*dx+dy*dy; | |
307 | if (dij < d2) | |
308 | { | |
309 | nucleonB->Collide(); | |
310 | nucleonA->Collide(); | |
311 | } | |
312 | } | |
bb442091 | 313 | } |
106d1ca1 | 314 | |
ec852657 | 315 | return CalcResults(bgen); |
316 | } | |
317 | ||
318 | //______________________________________________________________________________ | |
319 | Bool_t AliGlauberMC::CalcResults(Double_t bgen) | |
320 | { | |
bb442091 | 321 | // calc results for the given event |
322 | //return true if we have participants | |
106d1ca1 | 323 | |
324 | fNpart=0.; | |
325 | fNcoll=0.; | |
326 | fMeanX2=0.; | |
327 | fMeanY2=0.; | |
328 | fMeanXY=0.; | |
329 | fMeanXParts=0.; | |
330 | fMeanYParts=0.; | |
331 | fMeanXColl=0.; | |
332 | fMeanYColl=0.; | |
333 | fMeanX2Coll=0.; | |
334 | fMeanY2Coll=0.; | |
335 | fMeanXYColl=0.; | |
336 | fMeanXSystem=0.; | |
337 | fMeanYSystem=0.; | |
338 | fMeanXA=0.; | |
339 | fMeanYA=0.; | |
340 | fMeanXB=0.; | |
341 | fMeanYB=0.; | |
342 | fMeanr2=0.; | |
343 | fMeanr3=0.; | |
344 | fMeanr4=0.; | |
345 | fMeanr5=0.; | |
346 | fMeanr2Cos2Phi=0.; | |
347 | fMeanr2Sin2Phi=0.; | |
348 | fMeanr2Cos3Phi=0.; | |
349 | fMeanr2Sin3Phi=0.; | |
350 | fMeanr2Cos4Phi=0.; | |
351 | fMeanr2Sin4Phi=0.; | |
352 | fMeanr2Cos5Phi=0.; | |
353 | fMeanr2Sin5Phi=0.; | |
354 | fMeanr3Cos3Phi=0.; | |
355 | fMeanr3Sin3Phi=0.; | |
356 | fMeanr4Cos4Phi=0.; | |
357 | fMeanr4Sin4Phi=0.; | |
358 | fMeanr5Cos5Phi=0.; | |
359 | fMeanr5Sin5Phi=0.; | |
bb442091 | 360 | |
106d1ca1 | 361 | |
bb442091 | 362 | for (Int_t i = 0; i<fAN; i++) |
bb442091 | 363 | { |
106d1ca1 | 364 | AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i)); |
365 | Double_t xA = nucleonA->GetX(); | |
366 | Double_t yA = nucleonA->GetY(); | |
367 | Double_t r2A = xA*xA+yA*yA; | |
368 | Double_t rA = TMath::Sqrt(r2A); | |
369 | Double_t r3A = r2A*rA; | |
370 | Double_t r4A = r3A*rA; | |
371 | Double_t r5A = r4A*rA; | |
372 | Double_t phiA = TMath::ATan2(yA,xA); | |
373 | Double_t sin2PhiA = TMath::Sin(2.*phiA); | |
374 | Double_t cos2PhiA = TMath::Cos(2.*phiA); | |
375 | Double_t sin3PhiA = TMath::Sin(3.*phiA); | |
376 | Double_t cos3PhiA = TMath::Cos(3.*phiA); | |
377 | Double_t sin4PhiA = TMath::Sin(4.*phiA); | |
378 | Double_t cos4PhiA = TMath::Cos(4.*phiA); | |
379 | Double_t sin5PhiA = TMath::Sin(5.*phiA); | |
380 | Double_t cos5PhiA = TMath::Cos(5.*phiA); | |
381 | //printf("x: %.3f, y:%.3f, r:%.3f, phi:%.3f, sp:%.6f,cp:%.6f\n",xA,yA,TMath::Sqrt(r2A),PhiA,Sin2PhiA, Cos2PhiA); | |
382 | fMeanXSystem += xA; | |
383 | fMeanYSystem += yA; | |
384 | fMeanXA += xA; | |
385 | fMeanYA += yA; | |
386 | ||
387 | if(nucleonA->IsWounded()) | |
388 | { | |
389 | fNpart++; | |
390 | fMeanXParts += xA; | |
391 | fMeanYParts += yA; | |
392 | fMeanX2 += xA * xA; | |
393 | fMeanY2 += yA * yA; | |
394 | fMeanXY += xA * yA; | |
395 | fMeanr2 += r2A; | |
396 | fMeanr3 += r3A; | |
397 | fMeanr4 += r4A; | |
398 | fMeanr5 += r5A; | |
399 | fMeanr2Cos2Phi += r2A*cos2PhiA; | |
400 | fMeanr2Sin2Phi += r2A*sin2PhiA; | |
401 | fMeanr2Cos3Phi += r2A*cos3PhiA; | |
402 | fMeanr2Sin3Phi += r2A*sin3PhiA; | |
403 | fMeanr2Cos4Phi += r2A*cos4PhiA; | |
404 | fMeanr2Sin4Phi += r2A*sin4PhiA; | |
405 | fMeanr2Cos5Phi += r2A*cos5PhiA; | |
406 | fMeanr2Sin5Phi += r2A*sin5PhiA; | |
407 | fMeanr3Cos3Phi += r3A*cos3PhiA; | |
408 | fMeanr3Sin3Phi += r3A*sin3PhiA; | |
409 | fMeanr4Cos4Phi += r4A*cos4PhiA; | |
410 | fMeanr4Sin4Phi += r4A*sin4PhiA; | |
411 | fMeanr5Cos5Phi += r5A*cos5PhiA; | |
412 | fMeanr5Sin5Phi += r5A*sin5PhiA; | |
413 | } | |
bb442091 | 414 | } |
106d1ca1 | 415 | |
bb442091 | 416 | for (Int_t i = 0; i<fBN; i++) |
bb442091 | 417 | { |
106d1ca1 | 418 | AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i)); |
419 | Double_t xB=nucleonB->GetX(); | |
420 | Double_t yB=nucleonB->GetY(); | |
421 | Double_t r2B=xB*xB+yB*yB; | |
422 | Double_t rB = TMath::Sqrt(r2B); | |
423 | Double_t r3B = r2B*rB; | |
424 | Double_t r4B = r3B*rB; | |
425 | Double_t r5B = r4B*rB; | |
426 | Double_t phiB = TMath::ATan2(yB,xB); | |
427 | Double_t sin2PhiB = TMath::Sin(2*phiB); | |
428 | Double_t cos2PhiB = TMath::Cos(2*phiB); | |
429 | Double_t sin3PhiB = TMath::Sin(3*phiB); | |
430 | Double_t cos3PhiB = TMath::Cos(3*phiB); | |
431 | Double_t sin4PhiB = TMath::Sin(4*phiB); | |
432 | Double_t cos4PhiB = TMath::Cos(4*phiB); | |
433 | Double_t sin5PhiB = TMath::Sin(5*phiB); | |
434 | Double_t cos5PhiB = TMath::Cos(5*phiB); | |
435 | fMeanXSystem += xB; | |
436 | fMeanYSystem += yB; | |
437 | fMeanXB += xB; | |
438 | fMeanYB += yB; | |
439 | ||
440 | if(nucleonB->IsWounded()) | |
441 | { | |
442 | Int_t ncoll = nucleonB->GetNColl(); | |
443 | fNpart++; | |
444 | fMeanXParts += xB; | |
445 | fMeanXColl += xB*ncoll; | |
446 | fMeanYParts += yB; | |
447 | fMeanYColl += yB*ncoll; | |
448 | fMeanX2 += xB * xB; | |
449 | fMeanX2Coll += xB*xB*ncoll*ncoll; | |
450 | fMeanY2 += yB * yB; | |
451 | fMeanY2Coll += yB*yB*ncoll*ncoll; | |
452 | fMeanXY += xB * yB; | |
453 | fMeanXYColl += xB*yB*ncoll*ncoll; | |
454 | fNcoll += nucleonB->GetNColl(); | |
455 | fMeanr2 += r2B; | |
456 | fMeanr2Cos2Phi += r2B*cos2PhiB; | |
457 | fMeanr2Sin2Phi += r2B*sin2PhiB; | |
458 | fMeanr2Cos3Phi += r2B*cos3PhiB; | |
459 | fMeanr2Sin3Phi += r2B*sin3PhiB; | |
460 | fMeanr2Cos4Phi += r2B*cos4PhiB; | |
461 | fMeanr2Sin4Phi += r2B*sin4PhiB; | |
462 | fMeanr2Cos5Phi += r2B*cos5PhiB; | |
463 | fMeanr2Sin5Phi += r2B*sin5PhiB; | |
464 | fMeanr3Cos3Phi += r3B*cos3PhiB; | |
465 | fMeanr3Sin3Phi += r3B*sin3PhiB; | |
466 | fMeanr4Cos4Phi += r4B*cos4PhiB; | |
467 | fMeanr4Sin4Phi += r4B*sin4PhiB; | |
468 | fMeanr5Cos5Phi += r5B*cos5PhiB; | |
469 | fMeanr5Sin5Phi += r5B*sin5PhiB; | |
470 | } | |
bb442091 | 471 | } |
106d1ca1 | 472 | |
bb442091 | 473 | if (fNpart>0) |
106d1ca1 | 474 | { |
475 | fMeanXParts /= fNpart; | |
476 | fMeanYParts /= fNpart; | |
477 | fMeanX2 /= fNpart; | |
478 | fMeanY2 /= fNpart; | |
479 | fMeanXY /= fNpart; | |
480 | fMeanr2 /= fNpart; | |
481 | fMeanr2Cos2Phi /= fNpart; | |
482 | fMeanr2Sin2Phi /= fNpart; | |
483 | fMeanr2Cos3Phi /= fNpart; | |
484 | fMeanr2Sin3Phi /= fNpart; | |
485 | fMeanr2Cos4Phi /= fNpart; | |
486 | fMeanr2Sin4Phi /= fNpart; | |
487 | fMeanr2Cos5Phi /= fNpart; | |
488 | fMeanr2Sin5Phi /= fNpart; | |
489 | fMeanr3Cos3Phi /= fNpart; | |
490 | fMeanr3Sin3Phi /= fNpart; | |
491 | fMeanr4Cos4Phi /= fNpart; | |
492 | fMeanr4Sin4Phi /= fNpart; | |
493 | fMeanr5Cos5Phi /= fNpart; | |
494 | fMeanr5Sin5Phi /= fNpart; | |
495 | } | |
bb442091 | 496 | else |
106d1ca1 | 497 | { |
498 | fMeanXParts = 0; | |
499 | fMeanYParts = 0; | |
500 | fMeanX2 = 0; | |
501 | fMeanY2 = 0; | |
502 | fMeanXY = 0; | |
503 | fMeanr2 = 0; | |
504 | fMeanr2Cos2Phi = 0; | |
505 | fMeanr2Sin2Phi = 0; | |
506 | fMeanr2Cos3Phi = 0; | |
507 | fMeanr2Sin3Phi = 0; | |
508 | fMeanr2Cos4Phi = 0; | |
509 | fMeanr2Sin4Phi = 0; | |
510 | fMeanr2Cos5Phi = 0; | |
511 | fMeanr2Sin5Phi = 0; | |
512 | fMeanr3Cos3Phi = 0; | |
513 | fMeanr3Sin3Phi = 0; | |
514 | fMeanr4Cos4Phi = 0; | |
515 | fMeanr4Sin4Phi = 0; | |
516 | fMeanr5Cos5Phi = 0; | |
517 | fMeanr5Sin5Phi = 0; | |
518 | } | |
519 | ||
bb442091 | 520 | if (fNcoll>0) |
106d1ca1 | 521 | { |
522 | fMeanXColl /= fNcoll; | |
523 | fMeanYColl /= fNcoll; | |
524 | fMeanX2Coll /= fNcoll; | |
525 | fMeanY2Coll /= fNcoll; | |
526 | fMeanXYColl /= fNcoll; | |
527 | } | |
bb442091 | 528 | else |
106d1ca1 | 529 | { |
530 | fMeanXColl = 0; | |
531 | fMeanYColl = 0; | |
532 | fMeanX2Coll = 0; | |
533 | fMeanY2Coll = 0; | |
534 | fMeanXYColl = 0; | |
535 | } | |
536 | ||
bb442091 | 537 | if(fAN+fBN>0) |
106d1ca1 | 538 | { |
539 | fMeanXSystem /= (fAN + fBN); | |
540 | fMeanYSystem /= (fAN + fBN); | |
541 | } | |
bb442091 | 542 | else |
106d1ca1 | 543 | { |
544 | fMeanXSystem = 0; | |
545 | fMeanYSystem = 0; | |
546 | } | |
bb442091 | 547 | if(fAN>0) |
106d1ca1 | 548 | { |
549 | fMeanXA /= fAN; | |
550 | fMeanYA /= fAN; | |
551 | } | |
bb442091 | 552 | else |
106d1ca1 | 553 | { |
554 | fMeanXA = 0; | |
555 | fMeanYA = 0; | |
556 | } | |
557 | ||
bb442091 | 558 | if(fBN>0) |
106d1ca1 | 559 | { |
560 | fMeanXB /= fBN; | |
561 | fMeanYB /= fBN; | |
562 | } | |
bb442091 | 563 | else |
106d1ca1 | 564 | { |
565 | fMeanXB = 0; | |
566 | fMeanYB = 0; | |
567 | } | |
568 | ||
569 | //fEcc = Math::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2; | |
bb442091 | 570 | fSx2=fMeanX2-(fMeanXParts*fMeanXParts); |
571 | fSy2=fMeanY2-(fMeanYParts*fMeanYParts); | |
572 | fSxy=fMeanXY-fMeanXParts*fMeanYParts; | |
106d1ca1 | 573 | |
bb442091 | 574 | fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl); |
575 | fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl); | |
576 | fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl; | |
106d1ca1 | 577 | //fPsi3 = (TMath::Atan2(fMeanr2Sin3Phi,fMeanr2Cos3Phi)+3.1415926)/3.; |
578 | //printf("e:%.5f\n",TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2); | |
43215add | 579 | fBMC = bgen; |
bb442091 | 580 | |
581 | fTotalEvents++; | |
582 | if (fNpart>0) fEvents++; | |
583 | ||
584 | if (fNpart==0) return kFALSE; | |
585 | if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart; | |
586 | ||
106d1ca1 | 587 | |
bb442091 | 588 | return kTRUE; |
ec852657 | 589 | } |
590 | ||
591 | //______________________________________________________________________________ | |
592 | void AliGlauberMC::Draw(Option_t* /*option*/) | |
593 | { | |
43215add | 594 | //draw method |
bb442091 | 595 | fANucleus.Draw(fXSect, 2); |
596 | fBNucleus.Draw(fXSect, 4); | |
597 | ||
598 | TEllipse e; | |
599 | e.SetFillColor(0); | |
600 | e.SetLineColor(1); | |
601 | e.SetLineStyle(2); | |
602 | e.SetLineWidth(1); | |
603 | e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0); | |
604 | e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0); | |
ec852657 | 605 | } |
606 | ||
607 | //______________________________________________________________________________ | |
608 | Double_t AliGlauberMC::GetTotXSect() const | |
609 | { | |
43215add | 610 | //total xsection |
bb442091 | 611 | return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100; |
ec852657 | 612 | } |
613 | ||
614 | //______________________________________________________________________________ | |
615 | Double_t AliGlauberMC::GetTotXSectErr() const | |
616 | { | |
43215add | 617 | //total xsection error |
bb442091 | 618 | return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) * |
106d1ca1 | 619 | TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents)); |
ec852657 | 620 | } |
621 | ||
622 | //______________________________________________________________________________ | |
bb442091 | 623 | TObjArray *AliGlauberMC::GetNucleons() |
ec852657 | 624 | { |
43215add | 625 | //get array of nucleons |
bb442091 | 626 | if(!fNucleonsA || !fNucleonsB) return 0; |
627 | fNucleonsA->SetOwner(0); | |
628 | fNucleonsB->SetOwner(0); | |
629 | TObjArray *allnucleons=new TObjArray(fAN+fBN); | |
630 | allnucleons->SetOwner(); | |
631 | for (Int_t i = 0; i<fAN; i++) | |
106d1ca1 | 632 | { |
633 | allnucleons->Add(fNucleonsA->UncheckedAt(i)); | |
634 | } | |
bb442091 | 635 | for (Int_t i = 0; i<fBN; i++) |
106d1ca1 | 636 | { |
637 | allnucleons->Add(fNucleonsB->UncheckedAt(i)); | |
638 | } | |
bb442091 | 639 | return allnucleons; |
640 | } | |
641 | //______________________________________________________________________________ | |
642 | Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean) | |
643 | { | |
43215add | 644 | //produce a number distributed acc. neg.bin.dist |
bb442091 | 645 | if(k<=0) |
106d1ca1 | 646 | { |
647 | cout << "Error, zero or negative K" << endl; | |
648 | return 0; | |
649 | } | |
bb442091 | 650 | return (TMath::Binomial(x+k-1,x)) |
106d1ca1 | 651 | *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x)) |
652 | *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k)))); | |
bb442091 | 653 | } |
654 | //______________________________________________________________________________ | |
43215add | 655 | Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const |
bb442091 | 656 | { |
43215add | 657 | //return random integer from a Negative Binomial Distribution |
bb442091 | 658 | static const Int_t fMaxPlot = 1000; |
659 | Double_t array[fMaxPlot]; | |
660 | array[0] = NegativeBinomialDistribution(0,k,nmean); | |
661 | for (Int_t i=1; i<fMaxPlot; i++) | |
106d1ca1 | 662 | { |
663 | array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1]; | |
664 | } | |
bb442091 | 665 | Double_t r = gRandom->Uniform(0,1); |
666 | return TMath::BinarySearch(fMaxPlot,array,r)+2; | |
667 | } | |
668 | //______________________________________________________________________________ | |
669 | Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k, | |
670 | Double_t nmean, | |
671 | Int_t k2, | |
672 | Double_t nmean2, | |
43215add | 673 | Double_t alpha ) const |
bb442091 | 674 | { |
675 | //return random integer from a Double Negative Binomial Distribution | |
676 | static const Int_t fMaxPlot = 1000; | |
677 | Double_t array[fMaxPlot]; | |
678 | array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2); | |
679 | for (Int_t i=1; i<fMaxPlot; i++) | |
106d1ca1 | 680 | { |
681 | array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1]; | |
682 | } | |
bb442091 | 683 | Double_t r = gRandom->Uniform(0,1); |
684 | return TMath::BinarySearch(fMaxPlot,array,r)+2; | |
685 | } | |
686 | ||
bb442091 | 687 | //______________________________________________________________________________ |
688 | void AliGlauberMC::SetdNdEtaParam(Double_t nnp, Double_t x) | |
689 | { | |
690 | // set parameters used for calculating multiplicity, see GetdNdEta() for commments | |
691 | fdNdEtaParam[0]=nnp; | |
692 | fdNdEtaParam[1]=x; | |
ec852657 | 693 | } |
694 | ||
1b4934d5 | 695 | //______________________________________________________________________________ |
bb442091 | 696 | void AliGlauberMC::SetdNdEtaGBWParam(Double_t delta, Double_t lambda, Double_t snn) |
697 | { | |
698 | // set parameters used for calculating multiplicity see GetdNdEtaGBW() for commments | |
699 | fdNdEtaGBWParam[0]=delta; | |
700 | fdNdEtaGBWParam[1]=lambda; | |
701 | fdNdEtaGBWParam[2]=snn; | |
702 | } | |
14244e14 | 703 | |
bb442091 | 704 | //______________________________________________________________________________ |
705 | void AliGlauberMC::SetdNdEtaNBDParam(Double_t k, Double_t nmean, Double_t beta) | |
706 | { | |
707 | // set parameters used for calculating multiplicity see GetdNdEtaNBD() for commments | |
708 | fdNdEtaNBDParam[0]=k; | |
709 | fdNdEtaNBDParam[1]=nmean; | |
710 | fdNdEtaNBDParam[2]=beta; | |
711 | } | |
14244e14 | 712 | |
bb442091 | 713 | //______________________________________________________________________________ |
566fc3d0 | 714 | void AliGlauberMC::SetdNdEtaTwoNBDParam( Double_t k1, |
715 | Double_t nmean1, | |
716 | Double_t k2, | |
717 | Double_t nmean2, | |
718 | Double_t alpha, | |
719 | Double_t beta) | |
bb442091 | 720 | { |
721 | // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments | |
566fc3d0 | 722 | fdNdEtaTwoNBDParam[0]=k1; |
723 | fdNdEtaTwoNBDParam[1]=nmean1; | |
724 | fdNdEtaTwoNBDParam[2]=k2; | |
725 | fdNdEtaTwoNBDParam[3]=nmean2; | |
726 | fdNdEtaTwoNBDParam[4]=alpha; | |
bb442091 | 727 | fdNdEtaTwoNBDParam[5]=beta; |
728 | } | |
14244e14 | 729 | |
bb442091 | 730 | //______________________________________________________________________________ |
43215add | 731 | Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x) const |
1b4934d5 | 732 | { |
733 | //Get particle density per unit of rapidity | |
bb442091 | 734 | //using two component model |
735 | //Parameters: npp, x | |
736 | return nnp*((1.-x)*fNpart/2.+x*fNcoll); | |
0448b9f3 | 737 | } |
738 | ||
739 | //______________________________________________________________________________ | |
bb442091 | 740 | Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta, |
741 | Double_t lambda, | |
43215add | 742 | Double_t snn) const |
0448b9f3 | 743 | { |
744 | //Get particle density per unit of rapidity | |
745 | //using the GBW model | |
bb442091 | 746 | //Parameters: delta, lambda, snn |
747 | return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda)) | |
106d1ca1 | 748 | * TMath::Power(fNpart,(1.-delta)/3./delta); |
1b4934d5 | 749 | } |
750 | ||
14244e14 | 751 | //_______________________________________________________________________________ |
43215add | 752 | Double_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta) const |
bb442091 | 753 | { |
754 | //Get particle density per unit of rapidity | |
755 | //using a aandomized number from a negative binomial distrubution | |
756 | //Parameters: k = related to distribition width | |
757 | // nmean = mean of distribution | |
758 | // beta = set contribution of participants / binary collisions to multiplicity | |
759 | Double_t mulnp=0.0; | |
760 | for(Int_t i = 0; i<fNpart; i++) | |
106d1ca1 | 761 | { |
762 | mulnp+=NegativeBinomialRandom(k,nmean); | |
763 | } | |
bb442091 | 764 | Double_t mulnb=0.0; |
765 | for(Int_t i = 0; i<fNcoll; i++) | |
106d1ca1 | 766 | { |
767 | mulnb+=NegativeBinomialRandom(k,nmean); | |
768 | } | |
bb442091 | 769 | return (1-beta)*mulnp/2+beta*mulnb; |
770 | } | |
bb442091 | 771 | |
14244e14 | 772 | //______________________________________________________________________________ |
566fc3d0 | 773 | Double_t AliGlauberMC::GetdNdEtaTwoNBD ( Int_t k1, |
774 | Double_t nmean1, | |
775 | Int_t k2, | |
776 | Double_t nmean2, | |
777 | Double_t alpha, | |
43215add | 778 | Double_t beta ) const |
bb442091 | 779 | { |
780 | //Get particle density per unit of rapidity | |
781 | //using random numbers from two negative binomial distributions | |
782 | //Parameters: k1 = related to distribition width of distribution 1 | |
783 | // nmean1 = mean of distribution 1 | |
784 | // k2 = related to distribition width of distribution 2 | |
785 | // nmean2 = mean of distribution 2 | |
786 | // alpha = set contributions of distrubitin 1 / distribution 2 | |
787 | // beta = set contribution of participants / binary collisions to multiplicity | |
788 | Double_t mulnp=0.0; | |
789 | for(Int_t i = 0; i<fNpart; i++) | |
790 | { | |
791 | mulnp+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha); | |
792 | } | |
793 | Double_t mulnb=0.0; | |
794 | for(Int_t i = 0; i<fNcoll; i++) | |
795 | { | |
796 | mulnb+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha); | |
797 | } | |
798 | Double_t mul=(1-beta)*mulnp/2+beta*mulnb; | |
799 | return mul; | |
bb442091 | 800 | } |
14244e14 | 801 | |
1b4934d5 | 802 | //______________________________________________________________________________ |
803 | Double_t AliGlauberMC::GetEccentricityPart() const | |
804 | { | |
bb442091 | 805 | //get participant eccentricity of participants |
a0278ee7 | 806 | if (fNpart<2) return 0.0; |
1b4934d5 | 807 | return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2)); |
808 | } | |
809 | ||
106d1ca1 | 810 | //_____________________________________________________________________________ |
811 | ||
812 | ||
bb442091 | 813 | Double_t AliGlauberMC::GetEccentricityPartColl() const |
814 | { | |
815 | //get participant eccentricity of binary collisions | |
a0278ee7 | 816 | if (fNcoll<2) return 0.0; |
bb442091 | 817 | return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll)); |
818 | } | |
819 | ||
1b4934d5 | 820 | //______________________________________________________________________________ |
821 | Double_t AliGlauberMC::GetEccentricity() const | |
822 | { | |
bb442091 | 823 | //get standard eccentricity of participants |
74d0740e | 824 | if (fNpart<2) return 0.0; |
1b4934d5 | 825 | return ((fSy2-fSx2)/(fSy2+fSx2)); |
826 | } | |
bb442091 | 827 | |
828 | //______________________________________________________________________________ | |
829 | Double_t AliGlauberMC::GetEccentricityColl() const | |
830 | { | |
831 | //get standard eccentricity of binary collisions | |
b7166644 | 832 | if (fNcoll<2) return 0.0; |
bb442091 | 833 | return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll)); |
834 | } | |
14244e14 | 835 | |
ec852657 | 836 | //______________________________________________________________________________ |
837 | Bool_t AliGlauberMC::NextEvent(Double_t bgen) | |
838 | { | |
1b4934d5 | 839 | //Make a new event |
bb442091 | 840 | Int_t nAttempts = 10; // set indices, max attempts and max comparisons |
841 | Bool_t succes = kFALSE; | |
842 | for(Int_t j=0; j<nAttempts; j++) | |
843 | { | |
844 | if(bgen<0||!succes) //get impactparameter | |
845 | { | |
846 | bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin); | |
847 | } | |
97d98e48 | 848 | if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts |
bb442091 | 849 | } |
850 | return succes; | |
ec852657 | 851 | } |
852 | ||
853 | //______________________________________________________________________________ | |
106d1ca1 | 854 | |
855 | Double_t AliGlauberMC::GetEpsilon2Part() const | |
856 | { | |
857 | //get participant eccentricity of participants | |
858 | if (fNpart<2) return 0.0; | |
859 | return (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2); | |
860 | } | |
861 | ||
862 | //______________________________________________________________________________ | |
863 | ||
864 | Double_t AliGlauberMC::GetEpsilon3Part() const | |
865 | { | |
866 | //get participant eccentricity of participants | |
867 | if (fNpart<2) return 0.0; | |
868 | return (TMath::Sqrt(fMeanr2Cos3Phi*fMeanr2Cos3Phi+fMeanr2Sin3Phi*fMeanr2Sin3Phi)/fMeanr2); | |
869 | //return (TMath::Sqrt(fMeanr3Cos3Phi*fMeanr3Cos3Phi+fMeanr3Sin3Phi*fMeanr3Sin3Phi)/fMeanr3); | |
870 | } | |
871 | ||
872 | //______________________________________________________________________________ | |
873 | ||
874 | Double_t AliGlauberMC::GetEpsilon4Part() const | |
875 | { | |
876 | //get participant eccentricity of participants | |
877 | if (fNpart<2) return 0.0; | |
878 | return (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2); | |
879 | //return (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4); | |
880 | } | |
881 | ||
882 | //______________________________________________________________________________ | |
883 | ||
884 | Double_t AliGlauberMC::GetEpsilon5Part() const | |
885 | { | |
886 | //get participant eccentricity of participants | |
887 | if (fNpart<2) return 0.0; | |
888 | return (TMath::Sqrt(fMeanr2Cos5Phi*fMeanr2Cos5Phi+fMeanr2Sin5Phi*fMeanr2Sin5Phi)/fMeanr2); | |
889 | //return (TMath::Sqrt(fMeanr5Cos5Phi*fMeanr5Cos5Phi+fMeanr5Sin5Phi*fMeanr5Sin5Phi)/fMeanr5); | |
890 | } | |
891 | //______________________________________________________________________________ | |
892 | ||
ec852657 | 893 | void AliGlauberMC::Run(Int_t nevents) |
894 | { | |
43215add | 895 | //example run |
bb442091 | 896 | cout << "Generating " << nevents << " events..." << endl; |
897 | TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName())); | |
898 | TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect)); | |
899 | if (fnt == 0) | |
900 | { | |
901 | fnt = new TNtuple(name,title, | |
106d1ca1 | 902 | "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:VarEColl:VarEPart:VarEPartColl:dNdEta:dNdEtaGBW:dNdEtaNBD:dNdEtaTwoNBD:xsect:tAA:Epsl2:Epsl3:Epsl4:Epsl5"); |
bb442091 | 903 | fnt->SetDirectory(0); |
904 | } | |
905 | Int_t q = 0; | |
906 | Int_t u = 0; | |
907 | for (Int_t i = 0; i<nevents; i++) | |
bb442091 | 908 | { |
106d1ca1 | 909 | |
910 | if(!NextEvent()) | |
911 | { | |
912 | u++; | |
913 | continue; | |
914 | } | |
915 | ||
916 | q++; | |
917 | //Float_t v[27]; | |
918 | Float_t v[31]; | |
919 | v[0] = GetNpart(); | |
920 | v[1] = GetNcoll(); | |
921 | v[2] = fBMC; | |
922 | v[3] = fMeanXParts; | |
923 | v[4] = fMeanYParts; | |
924 | v[5] = fMeanX2; | |
925 | v[6] = fMeanY2; | |
926 | v[7] = fMeanXY; | |
927 | v[8] = fSx2; | |
928 | v[9] = fSy2; | |
929 | v[10] = fSxy; | |
930 | v[11] = fMeanXSystem; | |
931 | v[12] = fMeanYSystem; | |
932 | v[13] = fMeanXA; | |
933 | v[14] = fMeanYA; | |
934 | v[15] = fMeanXB; | |
935 | v[16] = fMeanYB; | |
936 | v[17] = GetEccentricity(); | |
937 | v[18] = GetEccentricityColl(); | |
938 | v[19] = GetEccentricityPart(); | |
939 | v[20] = GetEccentricityPartColl(); | |
940 | if (fDoPartProd) | |
941 | { | |
942 | v[21] = GetdNdEta( fdNdEtaParam[0],fdNdEtaParam[1] ); | |
943 | v[22] = GetdNdEtaGBW( fdNdEtaGBWParam[0],fdNdEtaGBWParam[1],fdNdEtaGBWParam[2] ); | |
944 | v[23] = GetdNdEtaNBD( TMath::Nint(fdNdEtaNBDParam[0]), | |
945 | fdNdEtaNBDParam[1], | |
946 | fdNdEtaNBDParam[2] ); | |
947 | v[24] = GetdNdEtaTwoNBD( TMath::Nint(fdNdEtaTwoNBDParam[0]), | |
948 | fdNdEtaTwoNBDParam[1], | |
949 | TMath::Nint(fdNdEtaTwoNBDParam[2]), | |
950 | fdNdEtaTwoNBDParam[3], | |
951 | fdNdEtaTwoNBDParam[4], | |
952 | fdNdEtaTwoNBDParam[5] ); | |
953 | } | |
954 | else | |
955 | { | |
956 | v[21] = 0; | |
957 | v[22] = 0; | |
958 | v[23] = 0; | |
959 | v[24] = 0; | |
960 | } | |
961 | v[25]=fXSect; | |
962 | ||
963 | Float_t mytAA=-999; | |
964 | if (GetNcoll()>0) mytAA=GetNcoll()/fXSect; | |
965 | v[26]=mytAA; | |
966 | //_____________epsilon2,3_______ | |
967 | v[27] = GetEpsilon2Part(); | |
968 | v[28] = GetEpsilon3Part(); | |
969 | v[29] = GetEpsilon4Part(); | |
970 | v[30] = GetEpsilon5Part(); | |
971 | //always at the end | |
a0278ee7 | 972 | |
106d1ca1 | 973 | fnt->Fill(v); |
14244e14 | 974 | |
106d1ca1 | 975 | if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush; |
976 | } | |
bb442091 | 977 | std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events: " << q << " discarded events: " << u <<"."<< endl; |
ec852657 | 978 | } |
979 | ||
980 | //--------------------------------------------------------------------------------- | |
43215add | 981 | void AliGlauberMC::RunAndSaveNtuple( Int_t n, |
982 | const Option_t *sysA, | |
983 | const Option_t *sysB, | |
ec852657 | 984 | Double_t signn, |
985 | Double_t mind, | |
106d1ca1 | 986 | Double_t r, |
987 | Double_t a, | |
ec852657 | 988 | const char *fname) |
989 | { | |
43215add | 990 | //example run |
566fc3d0 | 991 | AliGlauberMC mcg(sysA,sysB,signn); |
992 | mcg.SetMinDistance(mind); | |
14244e14 | 993 | mcg.Setr(r); |
994 | mcg.Seta(a); | |
566fc3d0 | 995 | mcg.Run(n); |
996 | TNtuple *nt=mcg.GetNtuple(); | |
bb442091 | 997 | TFile out(fname,"recreate",fname,9); |
998 | if(nt) nt->Write(); | |
14244e14 | 999 | printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect()); |
bb442091 | 1000 | out.Close(); |
ec852657 | 1001 | } |
1002 | ||
1003 | //--------------------------------------------------------------------------------- | |
43215add | 1004 | void AliGlauberMC::RunAndSaveNucleons( Int_t n, |
1005 | const Option_t *sysA, | |
bb442091 | 1006 | Option_t *sysB, |
43215add | 1007 | const Double_t signn, |
ec852657 | 1008 | Double_t mind, |
1009 | Bool_t verbose, | |
1010 | const char *fname) | |
1011 | { | |
43215add | 1012 | //example run |
566fc3d0 | 1013 | AliGlauberMC mcg(sysA,sysB,signn); |
1014 | mcg.SetMinDistance(mind); | |
bb442091 | 1015 | TFile *out=0; |
1016 | if(fname) | |
1017 | out=new TFile(fname,"recreate",fname,9); | |
1018 | ||
1019 | for(Int_t ievent=0; ievent<n; ievent++) | |
bb442091 | 1020 | { |
106d1ca1 | 1021 | |
1022 | //get an event with at least one collision | |
1023 | while(!mcg.NextEvent()) {} | |
1024 | ||
1025 | //access, save and (if wanted) print out nucleons | |
1026 | TObjArray* nucleons=mcg.GetNucleons(); | |
1027 | if(!nucleons) continue; | |
1028 | if(out) | |
1029 | nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey); | |
1030 | ||
1031 | if(verbose) | |
1032 | { | |
1033 | cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl; | |
1034 | cout<<"B = "<<mcg.GetB()<<" Npart = "<<mcg.GetNpart()<<endl<<endl; | |
1035 | printf("Nucleus\t X\t Y\t Z\tNcoll\n"); | |
1036 | Int_t nNucls=nucleons->GetEntries(); | |
1037 | for(Int_t iNucl=0; iNucl<nNucls; iNucl++) | |
1038 | { | |
1039 | AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl); | |
1040 | Char_t nucleus='A'; | |
1041 | if(nucl->IsInNucleusB()) nucleus='B'; | |
1042 | Double_t x=nucl->GetX(); | |
1043 | Double_t y=nucl->GetY(); | |
1044 | Double_t z=nucl->GetZ(); | |
1045 | Int_t ncoll=nucl->GetNColl(); | |
1046 | printf(" %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll); | |
1047 | } | |
1048 | } | |
bb442091 | 1049 | } |
bb442091 | 1050 | if(out) delete out; |
ec852657 | 1051 | } |
1052 | ||
566fc3d0 | 1053 | //--------------------------------------------------------------------------------- |
1054 | void AliGlauberMC::Reset() | |
1055 | { | |
1056 | //delete the ntuple | |
1057 | delete fnt; fnt=NULL; | |
1058 | } |