]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/Glauber/AliGlauberMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
:x
[u/mrichter/AliRoot.git] / PWG / Glauber / AliGlauberMC.cxx
CommitLineData
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
7288f660 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>
22e47719 34#include <TF1.h>
ec852657 35
36#include "AliGlauberNucleon.h"
37#include "AliGlauberNucleus.h"
38#include "AliGlauberMC.h"
39
3a7af7bd 40using std::flush;
ec852657 41ClassImp(AliGlauberMC)
42
43//______________________________________________________________________________
44AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
566fc3d0 45 TNamed(),
bb442091 46 fANucleus(NA),
47 fBNucleus(NB),
566fc3d0 48 fXSect(xsect),
bb442091 49 fNucleonsA(0),
50 fNucleonsB(0),
566fc3d0 51 fAN(0),
998de9b9 52 fQAN(0),
566fc3d0 53 fBN(0),
998de9b9 54 fQBN(0),
bb442091 55 fnt(0),
56 fMeanX2(0),
57 fMeanY2(0),
58 fMeanXY(0),
2c0ff1e6 59 fMeanX2Parts(0),
60 fMeanY2Parts(0),
61 fMeanXYParts(0),
bb442091 62 fMeanXParts(0),
63 fMeanYParts(0),
2c0ff1e6 64 fMeanOXParts(0),
65 fMeanOYParts(0),
bb442091 66 fMeanXColl(0),
67 fMeanYColl(0),
2c0ff1e6 68 fMeanOXColl(0),
69 fMeanOYColl(0),
bb442091 70 fMeanX2Coll(0),
2c0ff1e6 71 fMeanY2Coll(0),
bb442091 72 fMeanXYColl(0),
4d264dfa 73 fMeanXCom(0),
74 fMeanYCom(0),
75 fMeanOXCom(0),
76 fMeanOYCom(0),
77 fMeanX2Com(0),
78 fMeanY2Com(0),
79 fMeanXYCom(0),
bb442091 80 fMeanXSystem(0),
81 fMeanYSystem(0),
43215add 82 fMeanXA(0),
83 fMeanYA(0),
84 fMeanXB(0),
85 fMeanYB(0),
2c0ff1e6 86 fMeanOXA(0),
87 fMeanOYA(0),
88 fMeanOXB(0),
89 fMeanOYB(0),
43215add 90 fBMC(0),
bb442091 91 fEvents(0),
92 fTotalEvents(0),
566fc3d0 93 fBMin(0.),
94 fBMax(20.),
7288f660 95 fMultType(kNBDSV),
bb442091 96 fMaxNpartFound(0),
2c0ff1e6 97 fONpart(0),
98 fONcoll(0),
4d264dfa 99 fONcom(0),
bb442091 100 fNpart(0),
101 fNcoll(0),
95916d44 102 fNcollw(0),
4d264dfa 103 fNcom(0),
106d1ca1 104 fMeanr2(0),
105 fMeanr3(0),
106 fMeanr4(0),
107 fMeanr5(0),
108 fMeanr2Cos2Phi(0),
109 fMeanr2Sin2Phi(0),
110 fMeanr2Cos3Phi(0),
111 fMeanr2Sin3Phi(0),
112 fMeanr2Cos4Phi(0),
113 fMeanr2Sin4Phi(0),
114 fMeanr2Cos5Phi(0),
115 fMeanr2Sin5Phi(0),
116 fMeanr3Cos3Phi(0),
117 fMeanr3Sin3Phi(0),
118 fMeanr4Cos4Phi(0),
119 fMeanr4Sin4Phi(0),
120 fMeanr5Cos5Phi(0),
121 fMeanr5Sin5Phi(0),
9a5780fa 122 fMeanr2Coll(0),
123 fMeanr3Coll(0),
124 fMeanr4Coll(0),
125 fMeanr5Coll(0),
126 fMeanr2Cos2PhiColl(0),
127 fMeanr2Sin2PhiColl(0),
128 fMeanr2Cos3PhiColl(0),
129 fMeanr2Sin3PhiColl(0),
130 fMeanr2Cos4PhiColl(0),
131 fMeanr2Sin4PhiColl(0),
132 fMeanr2Cos5PhiColl(0),
133 fMeanr2Sin5PhiColl(0),
134 fMeanr3Cos3PhiColl(0),
135 fMeanr3Sin3PhiColl(0),
136 fMeanr4Cos4PhiColl(0),
137 fMeanr4Sin4PhiColl(0),
138 fMeanr5Cos5PhiColl(0),
139 fMeanr5Sin5PhiColl(0),
4d264dfa 140 fMeanr2Com(0),
141 fMeanr3Com(0),
142 fMeanr4Com(0),
143 fMeanr5Com(0),
144 fMeanr2Cos2PhiCom(0),
145 fMeanr2Sin2PhiCom(0),
146 fMeanr2Cos3PhiCom(0),
147 fMeanr2Sin3PhiCom(0),
148 fMeanr2Cos4PhiCom(0),
149 fMeanr2Sin4PhiCom(0),
150 fMeanr2Cos5PhiCom(0),
151 fMeanr2Sin5PhiCom(0),
152 fMeanr3Cos3PhiCom(0),
153 fMeanr3Sin3PhiCom(0),
154 fMeanr4Cos4PhiCom(0),
155 fMeanr4Sin4PhiCom(0),
156 fMeanr5Cos5PhiCom(0),
157 fMeanr5Sin5PhiCom(0),
2c0ff1e6 158 fSx2Parts(0.),
159 fSy2Parts(0.),
160 fSxyParts(0.),
566fc3d0 161 fSx2Coll(0.),
162 fSy2Coll(0.),
163 fSxyColl(0.),
4d264dfa 164 fSx2Com(0.),
165 fSy2Com(0.),
166 fSxyCom(0.),
bb442091 167 fX(0.13),
14244e14 168 fNpp(8.),
22e47719 169 fDoPartProd(kFALSE),
170 fBNN(0.),
171 fDoFluc(kFALSE),
172 fOmega(0),
173 fSig0(0),
cab4b5c8 174 fLambda(0),
175 fSigFluc(0)
ec852657 176{
43215add 177 //ctor
7288f660 178 for (UInt_t i=0; i<(sizeof(fdNdEtaParam)/sizeof(fdNdEtaParam[0])); i++)
179 {
180 fdNdEtaParam[i]=0.0;
181 }
bb442091 182
566fc3d0 183 SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
184 SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
185}
186
187//______________________________________________________________________________
188AliGlauberMC::~AliGlauberMC()
189{
190 //dtor
191 delete fnt;
192}
193
194//______________________________________________________________________________
195AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
196 TNamed(in),
197 fANucleus(in.fANucleus),
198 fBNucleus(in.fBNucleus),
199 fXSect(in.fXSect),
200 fNucleonsA(in.fNucleonsA),
201 fNucleonsB(in.fNucleonsB),
202 fAN(in.fAN),
998de9b9 203 fQAN(in.fQAN),
566fc3d0 204 fBN(in.fBN),
998de9b9 205 fQBN(in.fQBN),
566fc3d0 206 fnt(in.fnt),
207 fMeanX2(in.fMeanX2),
208 fMeanY2(in.fMeanY2),
209 fMeanXY(in.fMeanXY),
2c0ff1e6 210 fMeanX2Parts(in.fMeanX2Parts),
211 fMeanY2Parts(in.fMeanY2Parts),
212 fMeanXYParts(in.fMeanXYParts),
566fc3d0 213 fMeanXParts(in.fMeanXParts),
214 fMeanYParts(in.fMeanYParts),
2c0ff1e6 215 fMeanOXParts(in.fMeanOXParts),
216 fMeanOYParts(in.fMeanOYParts),
566fc3d0 217 fMeanXColl(in.fMeanXColl),
218 fMeanYColl(in.fMeanYColl),
2c0ff1e6 219 fMeanOXColl(in.fMeanOXColl),
220 fMeanOYColl(in.fMeanOYColl),
566fc3d0 221 fMeanX2Coll(in.fMeanX2Coll),
222 fMeanY2Coll(in.fMeanY2Coll),
223 fMeanXYColl(in.fMeanXYColl),
4d264dfa 224 fMeanXCom(in.fMeanXCom),
225 fMeanYCom(in.fMeanYCom),
226 fMeanOXCom(in.fMeanOXCom),
227 fMeanOYCom(in.fMeanOYCom),
228 fMeanX2Com(in.fMeanX2Com),
229 fMeanY2Com(in.fMeanY2Com),
230 fMeanXYCom(in.fMeanXYCom),
566fc3d0 231 fMeanXSystem(in.fMeanXSystem),
232 fMeanYSystem(in.fMeanYSystem),
43215add 233 fMeanXA(in.fMeanXA),
234 fMeanYA(in.fMeanYA),
235 fMeanXB(in.fMeanXB),
236 fMeanYB(in.fMeanYB),
2c0ff1e6 237 fMeanOXA(in.fMeanOXA),
238 fMeanOYA(in.fMeanOYA),
239 fMeanOXB(in.fMeanOXB),
240 fMeanOYB(in.fMeanOYB),
43215add 241 fBMC(in.fBMC),
566fc3d0 242 fEvents(in.fEvents),
243 fTotalEvents(in.fTotalEvents),
244 fBMin(in.fBMin),
245 fBMax(in.fBMax),
7288f660 246 fMultType(in.fMultType),
566fc3d0 247 fMaxNpartFound(in.fMaxNpartFound),
2c0ff1e6 248 fONpart(in.fONpart),
249 fONcoll(in.fONcoll),
4d264dfa 250 fONcom(in.fONcom),
566fc3d0 251 fNpart(in.fNpart),
252 fNcoll(in.fNcoll),
4d264dfa 253 fNcom(in.fNcom),
106d1ca1 254 fMeanr2(in.fMeanr2),
255 fMeanr3(in.fMeanr3),
256 fMeanr4(in.fMeanr4),
257 fMeanr5(in.fMeanr5),
258 fMeanr2Cos2Phi(in.fMeanr2Cos2Phi),
259 fMeanr2Sin2Phi(in.fMeanr2Sin2Phi),
260 fMeanr2Cos3Phi(in.fMeanr2Cos3Phi),
261 fMeanr2Sin3Phi(in.fMeanr2Sin3Phi),
262 fMeanr2Cos4Phi(in.fMeanr2Cos4Phi),
263 fMeanr2Sin4Phi(in.fMeanr2Sin4Phi),
264 fMeanr2Cos5Phi(in.fMeanr2Cos5Phi),
265 fMeanr2Sin5Phi(in.fMeanr2Sin5Phi),
266 fMeanr3Cos3Phi(in.fMeanr3Cos3Phi),
267 fMeanr3Sin3Phi(in.fMeanr3Sin3Phi),
268 fMeanr4Cos4Phi(in.fMeanr4Cos4Phi),
269 fMeanr4Sin4Phi(in.fMeanr4Sin4Phi),
270 fMeanr5Cos5Phi(in.fMeanr5Cos5Phi),
271 fMeanr5Sin5Phi(in.fMeanr5Sin5Phi),
9a5780fa 272 fMeanr2Coll(in.fMeanr2Coll),
273 fMeanr3Coll(in.fMeanr3Coll),
274 fMeanr4Coll(in.fMeanr4Coll),
275 fMeanr5Coll(in.fMeanr5Coll),
276 fMeanr2Cos2PhiColl(in.fMeanr2Cos2PhiColl),
277 fMeanr2Sin2PhiColl(in.fMeanr2Sin2PhiColl),
278 fMeanr2Cos3PhiColl(in.fMeanr2Cos3PhiColl),
279 fMeanr2Sin3PhiColl(in.fMeanr2Sin3PhiColl),
280 fMeanr2Cos4PhiColl(in.fMeanr2Cos4PhiColl),
281 fMeanr2Sin4PhiColl(in.fMeanr2Sin4PhiColl),
282 fMeanr2Cos5PhiColl(in.fMeanr2Cos5PhiColl),
283 fMeanr2Sin5PhiColl(in.fMeanr2Sin5PhiColl),
284 fMeanr3Cos3PhiColl(in.fMeanr3Cos3PhiColl),
285 fMeanr3Sin3PhiColl(in.fMeanr3Sin3PhiColl),
286 fMeanr4Cos4PhiColl(in.fMeanr4Cos4PhiColl),
287 fMeanr4Sin4PhiColl(in.fMeanr4Sin4PhiColl),
288 fMeanr5Cos5PhiColl(in.fMeanr5Cos5PhiColl),
289 fMeanr5Sin5PhiColl(in.fMeanr5Sin5PhiColl),
4d264dfa 290 fMeanr2Com(in.fMeanr2Com),
291 fMeanr3Com(in.fMeanr3Com),
292 fMeanr4Com(in.fMeanr4Com),
293 fMeanr5Com(in.fMeanr5Com),
294 fMeanr2Cos2PhiCom(in.fMeanr2Cos2PhiCom),
295 fMeanr2Sin2PhiCom(in.fMeanr2Sin2PhiCom),
296 fMeanr2Cos3PhiCom(in.fMeanr2Cos3PhiCom),
297 fMeanr2Sin3PhiCom(in.fMeanr2Sin3PhiCom),
298 fMeanr2Cos4PhiCom(in.fMeanr2Cos4PhiCom),
299 fMeanr2Sin4PhiCom(in.fMeanr2Sin4PhiCom),
300 fMeanr2Cos5PhiCom(in.fMeanr2Cos5PhiCom),
301 fMeanr2Sin5PhiCom(in.fMeanr2Sin5PhiCom),
302 fMeanr3Cos3PhiCom(in.fMeanr3Cos3PhiCom),
303 fMeanr3Sin3PhiCom(in.fMeanr3Sin3PhiCom),
304 fMeanr4Cos4PhiCom(in.fMeanr4Cos4PhiCom),
305 fMeanr4Sin4PhiCom(in.fMeanr4Sin4PhiCom),
306 fMeanr5Cos5PhiCom(in.fMeanr5Cos5PhiCom),
307 fMeanr5Sin5PhiCom(in.fMeanr5Sin5PhiCom),
2c0ff1e6 308 fSx2Parts(in.fSx2Parts),
309 fSy2Parts(in.fSy2Parts),
310 fSxyParts(in.fSxyParts),
566fc3d0 311 fSx2Coll(in.fSx2Coll),
312 fSy2Coll(in.fSy2Coll),
313 fSxyColl(in.fSxyColl),
4d264dfa 314 fSx2Com(in.fSx2Com),
315 fSy2Com(in.fSy2Com),
316 fSxyCom(in.fSxyCom),
566fc3d0 317 fX(in.fX),
97d98e48 318 fNpp(in.fNpp),
22e47719 319 fDoPartProd(in.fDoPartProd),
320 fBNN(in.fBNN),
321 fDoFluc(in.fDoFluc),
322 fOmega(in.fOmega),
323 fSig0(in.fSig0),
cab4b5c8 324 fLambda(in.fLambda),
325 fSigFluc(in.fSigFluc)
566fc3d0 326{
327 //copy ctor
7288f660 328 memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
566fc3d0 329}
330
331//______________________________________________________________________________
332AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
333{
334 //assignment
2047680a 335 if (&in==this) return *this;
7288f660 336 fANucleus=in.fANucleus;
337 fBNucleus=in.fBNucleus;
338 fXSect=in.fXSect;
566fc3d0 339 fNucleonsA=in.fNucleonsA;
340 fNucleonsB=in.fNucleonsB;
7288f660 341 fAN=in.fAN;
998de9b9 342 fQAN=in.fQAN;
7288f660 343 fBN=in.fBN;
998de9b9 344 fQBN=in.fQBN;
7288f660 345 fnt=in.fnt;
346 fMeanX2=in.fMeanX2;
347 fMeanY2=in.fMeanY2;
348 fMeanXY=in.fMeanXY;
106d1ca1 349 fMeanr2=in.fMeanr2;
350 fMeanr3=in.fMeanr3;
351 fMeanr4=in.fMeanr4;
352 fMeanr5=in.fMeanr5;
566fc3d0 353 fMeanXParts=in.fMeanXParts;
354 fMeanYParts=in.fMeanYParts;
2c0ff1e6 355 fMeanOXParts=in.fMeanOXParts;
356 fMeanOYParts=in.fMeanOYParts;
566fc3d0 357 fMeanXColl=in.fMeanXColl;
358 fMeanYColl=in.fMeanYColl;
2c0ff1e6 359 fMeanOXColl=in.fMeanOXColl;
360 fMeanOYColl=in.fMeanOYColl;
566fc3d0 361 fMeanX2Coll=in.fMeanX2Coll;
362 fMeanY2Coll=in.fMeanY2Coll;
363 fMeanXYColl=in.fMeanXYColl;
9a5780fa 364 fMeanr2Coll=in.fMeanr2Coll;
365 fMeanr3Coll=in.fMeanr3Coll;
366 fMeanr4Coll=in.fMeanr4Coll;
367 fMeanr5Coll=in.fMeanr5Coll;
4d264dfa 368 fMeanXCom=in.fMeanXColl;
369 fMeanYCom=in.fMeanYColl;
370 fMeanOXCom=in.fMeanOXCom;
371 fMeanOYCom=in.fMeanOYCom;
372 fMeanX2Com=in.fMeanX2Com;
373 fMeanY2Com=in.fMeanY2Com;
374 fMeanXYCom=in.fMeanXYCom;
375 fMeanr2Com=in.fMeanr2Com;
376 fMeanr3Com=in.fMeanr3Com;
377 fMeanr4Com=in.fMeanr4Com;
378 fMeanr5Com=in.fMeanr5Com;
566fc3d0 379 fMeanXSystem=in.fMeanXSystem;
2c0ff1e6 380 fMeanYSystem=in.fMeanYSystem;
7288f660 381 fMeanXA=in.fMeanXA;
382 fMeanYA=in.fMeanYA;
383 fMeanXB=in.fMeanXB;
384 fMeanYB=in.fMeanYB;
2c0ff1e6 385 fMeanOXA=in.fMeanOXA;
386 fMeanOYA=in.fMeanOYA;
387 fMeanOXB=in.fMeanOXB;
388 fMeanOYB=in.fMeanOYB;
7288f660 389 fBMC=in.fBMC;
390 fEvents=in.fEvents;
566fc3d0 391 fTotalEvents=in.fTotalEvents;
7288f660 392 fBMin=in.fBMin;
393 fBMax=in.fBMax;
394 fMultType=in.fMultType,
395 memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
566fc3d0 396 fMaxNpartFound=in.fMaxNpartFound;
7288f660 397 fNpart=in.fNpart;
398 fNcoll=in.fNcoll;
4d264dfa 399 fNcom=in.fNcom;
2c0ff1e6 400 fONpart=in.fONpart;
401 fONcoll=in.fONcoll;
4d264dfa 402 fONcom=in.fONcom;
106d1ca1 403 fMeanr2Cos2Phi=in.fMeanr2Cos2Phi;
404 fMeanr2Sin2Phi=in.fMeanr2Sin2Phi;
405 fMeanr2Cos3Phi=in.fMeanr2Cos3Phi;
406 fMeanr2Sin3Phi=in.fMeanr2Sin3Phi;
407 fMeanr2Cos4Phi=in.fMeanr2Cos4Phi;
408 fMeanr2Sin4Phi=in.fMeanr2Sin4Phi;
409 fMeanr2Cos5Phi=in.fMeanr2Cos5Phi;
410 fMeanr2Sin5Phi=in.fMeanr2Sin5Phi;
411 fMeanr3Cos3Phi=in.fMeanr3Cos3Phi;
412 fMeanr3Sin3Phi=in.fMeanr3Sin3Phi;
413 fMeanr4Cos4Phi=in.fMeanr4Cos4Phi;
414 fMeanr4Sin4Phi=in.fMeanr4Sin4Phi;
415 fMeanr5Cos5Phi=in.fMeanr5Cos5Phi;
416 fMeanr5Sin5Phi=in.fMeanr5Sin5Phi;
9a5780fa 417 fMeanr2Cos2PhiColl=in.fMeanr2Cos2PhiColl;
418 fMeanr2Sin2PhiColl=in.fMeanr2Sin2PhiColl;
419 fMeanr2Cos3PhiColl=in.fMeanr2Cos3PhiColl;
420 fMeanr2Sin3PhiColl=in.fMeanr2Sin3PhiColl;
421 fMeanr2Cos4PhiColl=in.fMeanr2Cos4PhiColl;
422 fMeanr2Sin4PhiColl=in.fMeanr2Sin4PhiColl;
423 fMeanr2Cos5PhiColl=in.fMeanr2Cos5PhiColl;
424 fMeanr2Sin5PhiColl=in.fMeanr2Sin5PhiColl;
425 fMeanr3Cos3PhiColl=in.fMeanr3Cos3PhiColl;
426 fMeanr3Sin3PhiColl=in.fMeanr3Sin3PhiColl;
427 fMeanr4Cos4PhiColl=in.fMeanr4Cos4PhiColl;
428 fMeanr4Sin4PhiColl=in.fMeanr4Sin4PhiColl;
429 fMeanr5Cos5PhiColl=in.fMeanr5Cos5PhiColl;
430 fMeanr5Sin5PhiColl=in.fMeanr5Sin5PhiColl;
4d264dfa 431 fMeanr2Cos2PhiCom=in.fMeanr2Cos2PhiCom;
432 fMeanr2Sin2PhiCom=in.fMeanr2Sin2PhiCom;
433 fMeanr2Cos3PhiCom=in.fMeanr2Cos3PhiCom;
434 fMeanr2Sin3PhiCom=in.fMeanr2Sin3PhiCom;
435 fMeanr2Cos4PhiCom=in.fMeanr2Cos4PhiCom;
436 fMeanr2Sin4PhiCom=in.fMeanr2Sin4PhiCom;
437 fMeanr2Cos5PhiCom=in.fMeanr2Cos5PhiCom;
438 fMeanr2Sin5PhiCom=in.fMeanr2Sin5PhiCom;
439 fMeanr3Cos3PhiCom=in.fMeanr3Cos3PhiCom;
440 fMeanr3Sin3PhiCom=in.fMeanr3Sin3PhiCom;
441 fMeanr4Cos4PhiCom=in.fMeanr4Cos4PhiCom;
442 fMeanr4Sin4PhiCom=in.fMeanr4Sin4PhiCom;
443 fMeanr5Cos5PhiCom=in.fMeanr5Cos5PhiCom;
444 fMeanr5Sin5PhiCom=in.fMeanr5Sin5PhiCom;
2c0ff1e6 445 fSx2Parts=in.fSx2Parts;
446 fSy2Parts=in.fSy2Parts;
447 fSxyParts=in.fSxyParts;
7288f660 448 fSx2Coll=in.fSx2Coll;
449 fSy2Coll=in.fSy2Coll;
450 fSxyColl=in.fSxyColl;
4d264dfa 451 fSx2Com=in.fSx2Com;
452 fSy2Com=in.fSy2Com;
453 fSxyCom=in.fSxyCom;
7288f660 454 fX=in.fX;
455 fNpp=in.fNpp;
566fc3d0 456 return *this;
ec852657 457}
458
459//______________________________________________________________________________
460Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
461{
bb442091 462 // prepare event
95916d44
CL		 463
464 if (fDoFluc) {
465 if (!fSigFluc) {
466 fSigFluc = new TF1("fSigFluc","[0]*x/[3]/(x/[3]+[1])*exp(-((x/[1]/[3]-1)/[2])^2)",0,250);
467 fSigFluc->SetParameters(1,fSig0,fOmega,fLambda);
468 cout << "Setting fluc: " << fSig0 << " " << fOmega << " " << fLambda << endl;
469 }
470 }
471
bb442091 472 fANucleus.ThrowNucleons(-bgen/2.);
473 fNucleonsA = fANucleus.GetNucleons();
474 fAN = fANucleus.GetN();
998de9b9 475 fQAN = fAN * 3;
476 //fAN = 3 * fANucleus.GetN(); // for Pb, Number of quark = 3*208;
bb442091 477 for (Int_t i = 0; i<fAN; i++)
7288f660 478 {
479 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
480 nucleonA->SetInNucleusA();
95916d44
CL		 481 nucleonA->SetSigNN(fXSect);
482 if (fDoFluc)
483 nucleonA->SetSigNN(fSigFluc->GetRandom());
7288f660 484 }
bb442091 485 fBNucleus.ThrowNucleons(bgen/2.);
486 fNucleonsB = fBNucleus.GetNucleons();
998de9b9 487 //fBN = 3 * fBNucleus.GetN(); // Number of quark = number of nucleus*3;
bb442091 488 fBN = fBNucleus.GetN();
998de9b9 489 fQBN = fBN * 3;
bb442091 490 for (Int_t i = 0; i<fBN; i++)
7288f660 491 {
492 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
493 nucleonB->SetInNucleusB();
95916d44
CL		 494 nucleonB->SetSigNN(fXSect);
495 if (fDoFluc)
496 nucleonB->SetSigNN(fSigFluc->GetRandom());
7288f660 497 }
498
22e47719 499 if (fDoFluc) {
500 if (!fSigFluc) {
501 fSigFluc = new TF1("fSigFluc","[0]*x/[3]/(x/[3]+[1])*exp(-((x/[1]/[3]-1)/[2])^2)",0,250);
502 fSigFluc->SetParameters(1,fSig0,fOmega,fLambda);
503 cout << "Setting fluc: " << fSig0 << " " << fOmega << " " << fLambda << endl;
504 }
505 fXSect = fSigFluc->GetRandom();
506 }
bb442091 507 // "ball" diameter = distance at which two balls interact
508 Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
7288f660 509
95916d44
CL		 510 Double_t bNN = 0;
511 Double_t Nco = 0;
512 Double_t Ncohc = 0; // hard core
513
bb442091 514 // for each of the A nucleons in nucleus B
515 for (Int_t i = 0; i<fBN; i++)
7288f660 516 {
517 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
518 for (Int_t j = 0 ; j < fAN ; j++)
bb442091 519 {
7288f660 520 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
521 Double_t dx = nucleonB->GetX()-nucleonA->GetX();
522 Double_t dy = nucleonB->GetY()-nucleonA->GetY();
523 Double_t dij = dx*dx+dy*dy;
95916d44
CL		 524 if (fDoFluc) {
525 //fXSect = nucleonA->GetSigNN();
526 //fXSect = (nucleonA->GetSigNN()+nucleonB->GetSigNN())/2.;
527 fXSect = TMath::Max(nucleonA->GetSigNN(),nucleonB->GetSigNN());
528 d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
529 }
7288f660 530 if (dij < d2)
531 {
22e47719 532 bNN += dij;
533 ++Nco;
7288f660 534 nucleonB->Collide();
535 nucleonA->Collide();
95916d44
CL		 536 if (dij<d2/4)
537 ++Ncohc;
7288f660 538 }
bb442091 539 }
7288f660 540 }
95916d44
CL		 541
542 if (Nco>0) {
543 fNcollw = Ncohc;
544 fBNN = bNN/Nco;
545 } else {
546 fNcollw = 0;
547 fBNN = 0.;
548 }
549
22e47719 550 if (Nco>0)
551 fBNN = bNN/Nco;
552 else
553 fBNN = 0.;
ec852657 554 return CalcResults(bgen);
555}
556
557//______________________________________________________________________________
558Bool_t AliGlauberMC::CalcResults(Double_t bgen)
559{
bb442091 560 // calc results for the given event
561 //return true if we have participants
7288f660 562
f98ea3f1 563 fNpart=0;
564 fNcoll=0;
4d264dfa 565 fNcom=0;
2c0ff1e6 566 fONpart=0;
567 fONcoll=0;
4d264dfa 568 fONcom=0;
106d1ca1 569 fMeanX2=0.;
570 fMeanY2=0.;
571 fMeanXY=0.;
572 fMeanXParts=0.;
573 fMeanYParts=0.;
2c0ff1e6 574 fMeanOXParts=0.;
575 fMeanOYParts=0.;
106d1ca1 576 fMeanXColl=0.;
577 fMeanYColl=0.;
2c0ff1e6 578 fMeanOXColl=0.;
579 fMeanOYColl=0.;
106d1ca1 580 fMeanX2Coll=0.;
581 fMeanY2Coll=0.;
582 fMeanXYColl=0.;
4d264dfa 583 fMeanXCom=0.;
584 fMeanYCom=0.;
585 fMeanOXCom=0.;
586 fMeanOYCom=0.;
587 fMeanX2Com=0.;
588 fMeanY2Com=0.;
589 fMeanXYCom=0.;
106d1ca1 590 fMeanXSystem=0.;
591 fMeanYSystem=0.;
592 fMeanXA=0.;
593 fMeanYA=0.;
594 fMeanXB=0.;
595 fMeanYB=0.;
2c0ff1e6 596 fMeanOXA=0.;
597 fMeanOYA=0.;
598 fMeanOXB=0.;
599 fMeanOYB=0.;
106d1ca1 600 fMeanr2=0.;
601 fMeanr3=0.;
602 fMeanr4=0.;
603 fMeanr5=0.;
604 fMeanr2Cos2Phi=0.;
605 fMeanr2Sin2Phi=0.;
606 fMeanr2Cos3Phi=0.;
607 fMeanr2Sin3Phi=0.;
608 fMeanr2Cos4Phi=0.;
609 fMeanr2Sin4Phi=0.;
610 fMeanr2Cos5Phi=0.;
611 fMeanr2Sin5Phi=0.;
612 fMeanr3Cos3Phi=0.;
613 fMeanr3Sin3Phi=0.;
614 fMeanr4Cos4Phi=0.;
615 fMeanr4Sin4Phi=0.;
616 fMeanr5Cos5Phi=0.;
617 fMeanr5Sin5Phi=0.;
9a5780fa 618 fMeanr2Coll=0.;
619 fMeanr3Coll=0.;
620 fMeanr4Coll=0.;
621 fMeanr5Coll=0.;
622 fMeanr2Cos2PhiColl=0.;
623 fMeanr2Sin2PhiColl=0.;
624 fMeanr2Cos3PhiColl=0.;
625 fMeanr2Sin3PhiColl=0.;
626 fMeanr2Cos4PhiColl=0.;
627 fMeanr2Sin4PhiColl=0.;
628 fMeanr2Cos5PhiColl=0.;
629 fMeanr2Sin5PhiColl=0.;
630 fMeanr3Cos3PhiColl=0.;
631 fMeanr3Sin3PhiColl=0.;
632 fMeanr4Cos4PhiColl=0.;
633 fMeanr4Sin4PhiColl=0.;
634 fMeanr5Cos5PhiColl=0.;
635 fMeanr5Sin5PhiColl=0.;
4d264dfa 636 fMeanr2Com=0.;
637 fMeanr3Com=0.;
638 fMeanr4Com=0.;
639 fMeanr5Com=0.;
640 fMeanr2Cos2PhiCom=0.;
641 fMeanr2Sin2PhiCom=0.;
642 fMeanr2Cos3PhiCom=0.;
643 fMeanr2Sin3PhiCom=0.;
644 fMeanr2Cos4PhiCom=0.;
645 fMeanr2Sin4PhiCom=0.;
646 fMeanr2Cos5PhiCom=0.;
647 fMeanr2Sin5PhiCom=0.;
648 fMeanr3Cos3PhiCom=0.;
649 fMeanr3Sin3PhiCom=0.;
650 fMeanr4Cos4PhiCom=0.;
651 fMeanr4Sin4PhiCom=0.;
652 fMeanr5Cos5PhiCom=0.;
653 fMeanr5Sin5PhiCom=0.;
7288f660 654
bb442091 655 for (Int_t i = 0; i<fAN; i++)
7288f660 656 {
657 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
2c0ff1e6 658 Double_t oXA = nucleonA->GetX();
659 Double_t oYA = nucleonA->GetY();
660 //fMeanOXSystem += oXA;
661 //fMeanOYSystem += oYA;
662 fMeanOXA += oXA;
663 fMeanOYA += oYA;
664
665 if(nucleonA->IsWounded())
666 {
667 fONpart++;
668 fMeanOXParts += oXA;
669 fMeanOYParts += oYA;
998de9b9 670
671 fONcom += (1-0.150);
672 fMeanOXCom += oXA*(1-0.150);
673 fMeanOYCom += oXA*(1-0.150);
2c0ff1e6 674 }
675 }
676
677 for (Int_t i = 0; i<fBN; i++)
678 {
679 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
680 Double_t oXB=nucleonB->GetX();
681 Double_t oYB=nucleonB->GetY();
682
683 if(nucleonB->IsWounded())
684 {
685 Int_t oNcoll = nucleonB->GetNColl();
686 fONpart++;
687 fMeanOXParts += oXB;
688 fMeanOXColl += oXB*oNcoll;
998de9b9 689 fMeanOXCom += oXB*((1-0.150)+0.150*oNcoll);
2c0ff1e6 690 fMeanOYParts += oYB;
691 fMeanOYColl += oYB*oNcoll;
998de9b9 692 fMeanOYCom += oYB*((1-0.150)+0.150*oNcoll);
2c0ff1e6 693 fONcoll += oNcoll;
998de9b9 694 fONcom += ((1-0.150)+0.150*oNcoll);
2c0ff1e6 695 }
696 }
697
698 if (fONpart>0)
699 {
700 fMeanOXParts /= fONpart;
701 fMeanOYParts /= fONpart;
702 }
703 else
704 {
705 fMeanOXParts = 0;
706 fMeanOYParts = 0;
707 }
708
709 if (fONcoll>0)
710 {
711 fMeanOXColl /= fONcoll;
712 fMeanOYColl /= fONcoll;
713 }
714 else
715 {
716 fMeanOXColl = 0;
717 fMeanOYColl = 0;
718 }
4d264dfa 719
720 if (fONcom>0)
721 {
722 fMeanOXCom /= fONcom;
723 fMeanOYCom /= fONcom;
724 }
725 else
726 {
727 fMeanOXCom = 0;
728 fMeanOYCom = 0;
729 }
2c0ff1e6 730
731 //////////////////////////////////////////////////////////////////
732 for (Int_t i = 0; i<fAN; i++)
733 {
734 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
735 Double_t xAA = nucleonA->GetX(); // X
736 Double_t yAA = nucleonA->GetY(); // Y
ddf04a47 737 Double_t xAPart = xAA - fMeanOXParts; // X'
738 Double_t yAPart = yAA - fMeanOYParts; // Y'
739 Double_t r2APart = xAPart *xAPart+yAPart*yAPart; // r'^2
740 Double_t rAPart = TMath::Sqrt(r2APart); // r'
741 Double_t r3APart = r2APart*rAPart;
742 Double_t r4APart = r3APart*rAPart;
743 Double_t r5APart = r4APart*rAPart;
744 Double_t phiAPart = TMath::ATan2(yAPart,xAPart);
745 Double_t sin2PhiAPart = TMath::Sin(2.*phiAPart);
746 Double_t cos2PhiAPart = TMath::Cos(2.*phiAPart);
747 Double_t sin3PhiAPart = TMath::Sin(3.*phiAPart);
748 Double_t cos3PhiAPart = TMath::Cos(3.*phiAPart);
749 Double_t sin4PhiAPart = TMath::Sin(4.*phiAPart);
750 Double_t cos4PhiAPart = TMath::Cos(4.*phiAPart);
751 Double_t sin5PhiAPart = TMath::Sin(5.*phiAPart);
752 Double_t cos5PhiAPart = TMath::Cos(5.*phiAPart);
753
998de9b9 754 // for combine
755 Double_t xACom = xAA - fMeanOXCom; // X'-Combine
756 Double_t yACom = yAA - fMeanOYCom; // Y'-C
757 Double_t r2ACom = xACom *xACom+yACom*yACom; // r'^2-C
758 Double_t rACom = TMath::Sqrt(r2ACom); // r'-C
759 Double_t r3ACom = r2ACom*rACom;
760 Double_t r4ACom = r3ACom*rACom;
761 Double_t r5ACom = r4ACom*rACom;
762 Double_t phiACom = TMath::ATan2(yACom,xACom);
763 Double_t sin2PhiACom = TMath::Sin(2.*phiACom);
764 Double_t cos2PhiACom = TMath::Cos(2.*phiACom);
765 Double_t sin3PhiACom = TMath::Sin(3.*phiACom);
766 Double_t cos3PhiACom = TMath::Cos(3.*phiACom);
767 Double_t sin4PhiACom = TMath::Sin(4.*phiACom);
768 Double_t cos4PhiACom = TMath::Cos(4.*phiACom);
769 Double_t sin5PhiACom = TMath::Sin(5.*phiACom);
770 Double_t cos5PhiACom = TMath::Cos(5.*phiACom);
771
2c0ff1e6 772 fMeanXSystem += xAA;
773 fMeanYSystem += yAA;
774 fMeanXA += xAA;
775 fMeanYA += yAA;
776 fMeanX2 += xAA * xAA;
777 fMeanY2 += yAA * yAA;
778 fMeanXY += xAA * yAA;
779
7288f660 780 if(nucleonA->IsWounded())
2c0ff1e6 781 {
998de9b9 782 //Wounded
7288f660 783 fNpart++;
ddf04a47 784 fMeanXParts += xAPart;
785 fMeanYParts += yAPart;
786 fMeanX2Parts += xAPart * xAPart;
787 fMeanY2Parts += yAPart * yAPart;
788 fMeanXYParts += xAPart * yAPart;
789 fMeanr2 += r2APart;
790 fMeanr3 += r3APart;
791 fMeanr4 += r4APart;
792 fMeanr5 += r5APart;
793 fMeanr2Cos2Phi += r2APart*cos2PhiAPart;
794 fMeanr2Sin2Phi += r2APart*sin2PhiAPart;
795 fMeanr2Cos3Phi += r2APart*cos3PhiAPart;
796 fMeanr2Sin3Phi += r2APart*sin3PhiAPart;
797 fMeanr2Cos4Phi += r2APart*cos4PhiAPart;
798 fMeanr2Sin4Phi += r2APart*sin4PhiAPart;
799 fMeanr2Cos5Phi += r2APart*cos5PhiAPart;
800 fMeanr2Sin5Phi += r2APart*sin5PhiAPart;
801 fMeanr3Cos3Phi += r3APart*cos3PhiAPart;
802 fMeanr3Sin3Phi += r3APart*sin3PhiAPart;
803 fMeanr4Cos4Phi += r4APart*cos4PhiAPart;
804 fMeanr4Sin4Phi += r4APart*sin4PhiAPart;
805 fMeanr5Cos5Phi += r5APart*cos5PhiAPart;
806 fMeanr5Sin5Phi += r5APart*sin5PhiAPart;
998de9b9 807
808 // Combined
809 fMeanXCom += xACom*(1-0.150);
810 fMeanYCom += yACom*(1-0.150);
811 fMeanX2Com += xACom*xACom*(1-0.150);
812 fMeanY2Com += yACom*yACom*(1-0.150);
813 fMeanXYCom += xACom*yACom*(1-0.150);
814 fNcom += (1-0.150);
815 fMeanr2Com += r2ACom*(1-0.150);
816 fMeanr3Com += r3ACom*(1-0.150);
817 fMeanr4Com += r4ACom*(1-0.150);
818 fMeanr5Com += r5ACom*(1-0.150);
819 fMeanr2Cos2PhiCom += r2ACom*cos2PhiACom*(1-0.150);
820 fMeanr2Sin2PhiCom += r2ACom*sin2PhiACom*(1-0.150);
821 fMeanr2Cos3PhiCom += r2ACom*cos3PhiACom*(1-0.150);
822 fMeanr2Sin3PhiCom += r2ACom*sin3PhiACom*(1-0.150);
823 fMeanr2Cos4PhiCom += r2ACom*cos4PhiACom*(1-0.150);
824 fMeanr2Sin4PhiCom += r2ACom*sin4PhiACom*(1-0.150);
825 fMeanr2Cos5PhiCom += r2ACom*cos5PhiACom*(1-0.150);
826 fMeanr2Sin5PhiCom += r2ACom*sin5PhiACom*(1-0.150);
827 fMeanr3Cos3PhiCom += r3ACom*cos3PhiACom*(1-0.150);
828 fMeanr3Sin3PhiCom += r3ACom*sin3PhiACom*(1-0.150);
829 fMeanr4Cos4PhiCom += r4ACom*cos4PhiACom*(1-0.150);
830 fMeanr4Sin4PhiCom += r4ACom*sin4PhiACom*(1-0.150);
831 fMeanr5Cos5PhiCom += r5ACom*cos5PhiACom*(1-0.150);
832 fMeanr5Sin5PhiCom += r5ACom*sin5PhiACom*(1-0.150);
833
bb442091 834 }
7288f660 835 }
ddf04a47 836
bb442091 837 for (Int_t i = 0; i<fBN; i++)
bb442091 838 {
2c0ff1e6 839 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
840 Double_t xBB = nucleonB->GetX();
841 Double_t yBB = nucleonB->GetY();
ddf04a47 842 // for Wounded
843 Double_t xBPart = xBB - fMeanOXParts; // X'
844 Double_t yBPart = yBB - fMeanOYParts; // Y'
845 Double_t r2BPart = xBPart*xBPart+yBPart*yBPart; // r'^2
846 Double_t rBPart = TMath::Sqrt(r2BPart); // r'
847 Double_t r3BPart = r2BPart*rBPart;
848 Double_t r4BPart = r3BPart*rBPart;
849 Double_t r5BPart = r4BPart*rBPart;
850 Double_t phiBPart = TMath::ATan2(yBPart,xBPart);
851 Double_t sin2PhiBPart = TMath::Sin(2.*phiBPart);
852 Double_t cos2PhiBPart = TMath::Cos(2.*phiBPart);
853 Double_t sin3PhiBPart = TMath::Sin(3.*phiBPart);
854 Double_t cos3PhiBPart = TMath::Cos(3.*phiBPart);
855 Double_t sin4PhiBPart = TMath::Sin(4.*phiBPart);
856 Double_t cos4PhiBPart = TMath::Cos(4.*phiBPart);
857 Double_t sin5PhiBPart = TMath::Sin(5.*phiBPart);
858 Double_t cos5PhiBPart = TMath::Cos(5.*phiBPart);
859 // for Binary
860 Double_t xBColl = xBB - fMeanOXColl; // X'-Binary
861 Double_t yBColl = yBB - fMeanOYColl; // Y'-B
862 Double_t r2BColl = xBColl*xBColl+yBColl*yBColl; // r'^2-B
863 Double_t rBColl = TMath::Sqrt(r2BColl); // r'-B
864 Double_t r3BColl = r2BColl*rBColl;
865 Double_t r4BColl = r3BColl*rBColl;
866 Double_t r5BColl = r4BColl*rBColl;
867 Double_t phiBColl = TMath::ATan2(yBColl,xBColl);
868 Double_t sin2PhiBColl = TMath::Sin(2.*phiBColl);
869 Double_t cos2PhiBColl = TMath::Cos(2.*phiBColl);
870 Double_t sin3PhiBColl = TMath::Sin(3.*phiBColl);
871 Double_t cos3PhiBColl = TMath::Cos(3.*phiBColl);
872 Double_t sin4PhiBColl = TMath::Sin(4.*phiBColl);
873 Double_t cos4PhiBColl = TMath::Cos(4.*phiBColl);
874 Double_t sin5PhiBColl = TMath::Sin(5.*phiBColl);
875 Double_t cos5PhiBColl = TMath::Cos(5.*phiBColl);
876 // for combine
877 Double_t xBCom = xBB - fMeanOXCom; // X'-Combine
878 Double_t yBCom = yBB - fMeanOYCom; // Y'-C
879 Double_t r2BCom = xBCom *xBCom+yBCom*yBCom; // r'^2-C
880 Double_t rBCom = TMath::Sqrt(r2BCom); // r'-C
881 Double_t r3BCom = r2BCom*rBCom;
882 Double_t r4BCom = r3BCom*rBCom;
883 Double_t r5BCom = r4BCom*rBCom;
884 Double_t phiBCom = TMath::ATan2(yBCom,xBCom);
885 Double_t sin2PhiBCom = TMath::Sin(2.*phiBCom);
886 Double_t cos2PhiBCom = TMath::Cos(2.*phiBCom);
887 Double_t sin3PhiBCom = TMath::Sin(3.*phiBCom);
888 Double_t cos3PhiBCom = TMath::Cos(3.*phiBCom);
889 Double_t sin4PhiBCom = TMath::Sin(4.*phiBCom);
890 Double_t cos4PhiBCom = TMath::Cos(4.*phiBCom);
891 Double_t sin5PhiBCom = TMath::Sin(5.*phiBCom);
892 Double_t cos5PhiBCom = TMath::Cos(5.*phiBCom);
893
2c0ff1e6 894 fMeanXSystem += xBB;
895 fMeanYSystem += yBB;
896 fMeanXB += xBB;
897 fMeanYB += yBB;
898 fMeanX2 += xBB*xBB;
899 fMeanY2 += yBB*yBB;
900 fMeanXY += xBB*yBB;
901
902 if(nucleonB->IsWounded())
903 {
904 Int_t ncoll = nucleonB->GetNColl();
905 fNpart++;
ddf04a47 906 fMeanXParts += xBPart;
907 fMeanXColl += xBColl*ncoll;
998de9b9 908 fMeanXCom += xBCom*((1-0.150)+0.150*ncoll);
ddf04a47 909 fMeanYParts += yBPart;
910 fMeanYColl += yBColl*ncoll;
998de9b9 911 fMeanYCom += yBCom*((1-0.150)+0.150*ncoll);
ddf04a47 912 fMeanX2Parts += xBPart * xBPart;
913 fMeanX2Coll += xBColl*xBColl*ncoll;
998de9b9 914 fMeanX2Com += xBCom*xBCom*((1-0.150)+0.150*ncoll);
ddf04a47 915 fMeanY2Parts += yBPart * yBPart;
916 fMeanY2Coll += yBColl*yBColl*ncoll;
998de9b9 917 fMeanY2Com += yBCom*yBCom*((1-0.150)+0.150*ncoll);
ddf04a47 918 fMeanXYParts += xBPart * yBPart;
919 fMeanXYColl += xBColl*yBColl*ncoll;
998de9b9 920 fMeanXYCom += xBCom*yBCom*((1-0.150)+0.150*ncoll);
2c0ff1e6 921 fNcoll += ncoll;
998de9b9 922 fNcom += ((1-0.150)+0.150*ncoll);
ddf04a47 923 fMeanr2 += r2BPart;
924 fMeanr3 += r3BPart;
925 fMeanr4 += r4BPart;
926 fMeanr5 += r5BPart;
927 fMeanr2Cos2Phi += r2BPart*cos2PhiBPart;
928 fMeanr2Sin2Phi += r2BPart*sin2PhiBPart;
929 fMeanr2Cos3Phi += r2BPart*cos3PhiBPart;
930 fMeanr2Sin3Phi += r2BPart*sin3PhiBPart;
931 fMeanr2Cos4Phi += r2BPart*cos4PhiBPart;
932 fMeanr2Sin4Phi += r2BPart*sin4PhiBPart;
933 fMeanr2Cos5Phi += r2BPart*cos5PhiBPart;
934 fMeanr2Sin5Phi += r2BPart*sin5PhiBPart;
935 fMeanr3Cos3Phi += r3BPart*cos3PhiBPart;
936 fMeanr3Sin3Phi += r3BPart*sin3PhiBPart;
937 fMeanr4Cos4Phi += r4BPart*cos4PhiBPart;
938 fMeanr4Sin4Phi += r4BPart*sin4PhiBPart;
939 fMeanr5Cos5Phi += r5BPart*cos5PhiBPart;
940 fMeanr5Sin5Phi += r5BPart*sin5PhiBPart;
941 fMeanr2Coll += r2BColl*ncoll;
942 fMeanr3Coll += r3BColl*ncoll;
943 fMeanr4Coll += r4BColl*ncoll;
944 fMeanr5Coll += r5BColl*ncoll;
945 fMeanr2Cos2PhiColl += r2BColl*cos2PhiBColl*ncoll;
946 fMeanr2Sin2PhiColl += r2BColl*sin2PhiBColl*ncoll;
947 fMeanr2Cos3PhiColl += r2BColl*cos3PhiBColl*ncoll;
948 fMeanr2Sin3PhiColl += r2BColl*sin3PhiBColl*ncoll;
949 fMeanr2Cos4PhiColl += r2BColl*cos4PhiBColl*ncoll;
950 fMeanr2Sin4PhiColl += r2BColl*sin4PhiBColl*ncoll;
951 fMeanr2Cos5PhiColl += r2BColl*cos5PhiBColl*ncoll;
952 fMeanr2Sin5PhiColl += r2BColl*sin5PhiBColl*ncoll;
953 fMeanr3Cos3PhiColl += r3BColl*cos3PhiBColl*ncoll;
954 fMeanr3Sin3PhiColl += r3BColl*sin3PhiBColl*ncoll;
955 fMeanr4Cos4PhiColl += r4BColl*cos4PhiBColl*ncoll;
956 fMeanr4Sin4PhiColl += r4BColl*sin4PhiBColl*ncoll;
957 fMeanr5Cos5PhiColl += r5BColl*cos5PhiBColl*ncoll;
958 fMeanr5Sin5PhiColl += r5BColl*sin5PhiBColl*ncoll;
998de9b9 959 fMeanr2Com += r2BCom*((1-0.150)+0.150*ncoll);
960 fMeanr3Com += r3BCom*((1-0.150)+0.150*ncoll);
961 fMeanr4Com += r4BCom*((1-0.150)+0.150*ncoll);
962 fMeanr5Com += r5BCom*((1-0.150)+0.150*ncoll);
963 fMeanr2Cos2PhiCom += r2BCom*cos2PhiBCom*((1-0.150)+0.150*ncoll);
964 fMeanr2Sin2PhiCom += r2BCom*sin2PhiBCom*((1-0.150)+0.150*ncoll);
965 fMeanr2Cos3PhiCom += r2BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
966 fMeanr2Sin3PhiCom += r2BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
967 fMeanr2Cos4PhiCom += r2BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
968 fMeanr2Sin4PhiCom += r2BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
969 fMeanr2Cos5PhiCom += r2BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
970 fMeanr2Sin5PhiCom += r2BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
971 fMeanr3Cos3PhiCom += r3BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
972 fMeanr3Sin3PhiCom += r3BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
973 fMeanr4Cos4PhiCom += r4BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
974 fMeanr4Sin4PhiCom += r4BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
975 fMeanr5Cos5PhiCom += r5BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
976 fMeanr5Sin5PhiCom += r5BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
2c0ff1e6 977 }
bb442091 978 }
2c0ff1e6 979
bb442091 980 if (fNpart>0)
2c0ff1e6 981 {
982 fMeanXParts /= fNpart;
983 fMeanYParts /= fNpart;
984 fMeanX2Parts /= fNpart;
985 fMeanY2Parts /= fNpart;
986 fMeanXYParts /= fNpart;
987 fMeanr2 /= fNpart;
988 fMeanr3 /= fNpart;
989 fMeanr4 /= fNpart;
990 fMeanr5 /= fNpart;
991 fMeanr2Cos2Phi /= fNpart;
992 fMeanr2Sin2Phi /= fNpart;
993 fMeanr2Cos3Phi /= fNpart;
994 fMeanr2Sin3Phi /= fNpart;
995 fMeanr2Cos4Phi /= fNpart;
996 fMeanr2Sin4Phi /= fNpart;
997 fMeanr2Cos5Phi /= fNpart;
998 fMeanr2Sin5Phi /= fNpart;
999 fMeanr3Cos3Phi /= fNpart;
1000 fMeanr3Sin3Phi /= fNpart;
1001 fMeanr4Cos4Phi /= fNpart;
1002 fMeanr4Sin4Phi /= fNpart;
1003 fMeanr5Cos5Phi /= fNpart;
1004 fMeanr5Sin5Phi /= fNpart;
1005 }
1006
bb442091 1007 else
2c0ff1e6 1008 {
1009 fMeanXParts = 0;
1010 fMeanYParts = 0;
1011 fMeanX2Parts = 0;
1012 fMeanY2Parts = 0;
1013 fMeanXYParts = 0;
1014 fMeanr2 = 0;
1015 fMeanr3 = 0;
1016 fMeanr4 = 0;
1017 fMeanr5 = 0;
1018 fMeanr2Cos2Phi = 0;
1019 fMeanr2Sin2Phi = 0;
1020 fMeanr2Cos3Phi = 0;
1021 fMeanr2Sin3Phi = 0;
1022 fMeanr2Cos4Phi = 0;
1023 fMeanr2Sin4Phi = 0;
1024 fMeanr2Cos5Phi = 0;
1025 fMeanr2Sin5Phi = 0;
1026 fMeanr3Cos3Phi = 0;
1027 fMeanr3Sin3Phi = 0;
1028 fMeanr4Cos4Phi = 0;
1029 fMeanr4Sin4Phi = 0;
1030 fMeanr5Cos5Phi = 0;
1031 fMeanr5Sin5Phi = 0;
1032 }
1033
bb442091 1034 if (fNcoll>0)
2c0ff1e6 1035 {
1036 fMeanXColl /= fNcoll;
1037 fMeanYColl /= fNcoll;
1038 fMeanX2Coll /= fNcoll;
1039 fMeanY2Coll /= fNcoll;
1040 fMeanXYColl /= fNcoll;
1041 fMeanr2Coll /= fNcoll;
1042 fMeanr3Coll /= fNcoll;
1043 fMeanr4Coll /= fNcoll;
1044 fMeanr5Coll /= fNcoll;
1045 fMeanr2Cos2PhiColl /= fNcoll;
1046 fMeanr2Sin2PhiColl /= fNcoll;
1047 fMeanr2Cos3PhiColl /= fNcoll;
1048 fMeanr2Sin3PhiColl /= fNcoll;
1049 fMeanr2Cos4PhiColl /= fNcoll;
1050 fMeanr2Sin4PhiColl /= fNcoll;
1051 fMeanr2Cos5PhiColl /= fNcoll;
1052 fMeanr2Sin5PhiColl /= fNcoll;
1053 fMeanr3Cos3PhiColl /= fNcoll;
1054 fMeanr3Sin3PhiColl /= fNcoll;
1055 fMeanr4Cos4PhiColl /= fNcoll;
1056 fMeanr4Sin4PhiColl /= fNcoll;
1057 fMeanr5Cos5PhiColl /= fNcoll;
1058 fMeanr5Sin5PhiColl /= fNcoll;
7288f660 1059 }
2c0ff1e6 1060
bb442091 1061 else
7288f660 1062 {
1063 fMeanXColl = 0;
1064 fMeanYColl = 0;
1065 fMeanX2Coll = 0;
1066 fMeanY2Coll = 0;
1067 fMeanXYColl = 0;
1068 fMeanr2Cos2PhiColl =0;
1069 fMeanr2Sin2PhiColl =0;
1070 fMeanr2Cos3PhiColl =0;
1071 fMeanr2Sin3PhiColl =0;
1072 fMeanr2Cos4PhiColl =0;
1073 fMeanr2Sin4PhiColl =0;
1074 fMeanr2Cos5PhiColl =0;
1075 fMeanr2Sin5PhiColl =0;
1076 fMeanr3Cos3PhiColl =0;
1077 fMeanr3Sin3PhiColl =0;
1078 fMeanr4Cos4PhiColl =0;
1079 fMeanr4Sin4PhiColl =0;
1080 fMeanr5Cos5PhiColl =0;
1081 fMeanr5Sin5PhiColl =0;
1082 }
2c0ff1e6 1083
4d264dfa 1084 if (fNcom>0)
1085 {
1086 fMeanXCom /= fNcom;
1087 fMeanYCom /= fNcom;
1088 fMeanX2Com /= fNcom;
1089 fMeanY2Com /= fNcom;
1090 fMeanXYCom /= fNcom;
1091 fMeanr2Com /= fNcom;
1092 fMeanr3Com /= fNcom;
1093 fMeanr4Com /= fNcom;
1094 fMeanr5Com /= fNcom;
1095 fMeanr2Cos2PhiCom /= fNcom;
1096 fMeanr2Sin2PhiCom /= fNcom;
1097 fMeanr2Cos3PhiCom /= fNcom;
1098 fMeanr2Sin3PhiCom /= fNcom;
1099 fMeanr2Cos4PhiCom /= fNcom;
1100 fMeanr2Sin4PhiCom /= fNcom;
1101 fMeanr2Cos5PhiCom /= fNcom;
1102 fMeanr2Sin5PhiCom /= fNcom;
1103 fMeanr3Cos3PhiCom /= fNcom;
1104 fMeanr3Sin3PhiCom /= fNcom;
1105 fMeanr4Cos4PhiCom /= fNcom;
1106 fMeanr4Sin4PhiCom /= fNcom;
1107 fMeanr5Cos5PhiCom /= fNcom;
1108 fMeanr5Sin5PhiCom /= fNcom;
1109 }
1110
1111 else
1112 {
1113 fMeanXCom = 0;
1114 fMeanYCom = 0;
1115 fMeanX2Com = 0;
1116 fMeanY2Com = 0;
1117 fMeanXYCom = 0;
1118 fMeanr2Cos2PhiCom =0;
1119 fMeanr2Sin2PhiCom =0;
1120 fMeanr2Cos3PhiCom =0;
1121 fMeanr2Sin3PhiCom =0;
1122 fMeanr2Cos4PhiCom =0;
1123 fMeanr2Sin4PhiCom =0;
1124 fMeanr2Cos5PhiCom =0;
1125 fMeanr2Sin5PhiCom =0;
1126 fMeanr3Cos3PhiCom =0;
1127 fMeanr3Sin3PhiCom =0;
1128 fMeanr4Cos4PhiCom =0;
1129 fMeanr4Sin4PhiCom =0;
1130 fMeanr5Cos5PhiCom =0;
1131 fMeanr5Sin5PhiCom =0;
1132 }
1133
bb442091 1134 if(fAN+fBN>0)
7288f660 1135 {
1136 fMeanXSystem /= (fAN + fBN);
1137 fMeanYSystem /= (fAN + fBN);
2c0ff1e6 1138 fMeanX2 /= (fAN + fBN);
1139 fMeanY2 /= (fAN + fBN);
1140 fMeanXY /= (fAN + fBN);
7288f660 1141 }
bb442091 1142 else
7288f660 1143 {
1144 fMeanXSystem = 0;
1145 fMeanYSystem = 0;
2c0ff1e6 1146 fMeanX2 = 0;
1147 fMeanY2 = 0;
1148 fMeanXY = 0;
7288f660 1149 }
bb442091 1150 if(fAN>0)
7288f660 1151 {
1152 fMeanXA /= fAN;
1153 fMeanYA /= fAN;
1154 }
bb442091 1155 else
7288f660 1156 {
1157 fMeanXA = 0;
1158 fMeanYA = 0;
1159 }
1160
bb442091 1161 if(fBN>0)
7288f660 1162 {
1163 fMeanXB /= fBN;
1164 fMeanYB /= fBN;
1165 }
bb442091 1166 else
7288f660 1167 {
1168 fMeanXB = 0;
1169 fMeanYB = 0;
1170 }
2c0ff1e6 1171
1172
7288f660 1173
2c0ff1e6 1174
1175 //////////////////////////////////////////////////////////////////
1176 fSx2Parts=fMeanX2Parts-(fMeanXParts*fMeanXParts);
1177 fSy2Parts=fMeanY2Parts-(fMeanYParts*fMeanYParts);
1178 fSxyParts=fMeanXYParts-fMeanXParts*fMeanYParts;
bb442091 1179 fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl);
1180 fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl);
1181 fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl;
4d264dfa 1182 fSx2Com=fMeanX2Com-(fMeanXCom*fMeanXCom);
1183 fSy2Com=fMeanY2Com-(fMeanYCom*fMeanYCom);
1184 fSxyCom=fMeanXYCom-fMeanXCom*fMeanYCom;
43215add 1185 fBMC = bgen;
bb442091 1186 fTotalEvents++;
1187 if (fNpart>0) fEvents++;
bb442091 1188 if (fNpart==0) return kFALSE;
1189 if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
1190
1191 return kTRUE;
ec852657 1192}
1193
1194//______________________________________________________________________________
1195void AliGlauberMC::Draw(Option_t* /*option*/)
1196{
43215add 1197 //draw method
bb442091 1198 fANucleus.Draw(fXSect, 2);
1199 fBNucleus.Draw(fXSect, 4);
1200
1201 TEllipse e;
1202 e.SetFillColor(0);
1203 e.SetLineColor(1);
1204 e.SetLineStyle(2);
1205 e.SetLineWidth(1);
1206 e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0);
1207 e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0);
ec852657 1208}
1209
1210//______________________________________________________________________________
1211Double_t AliGlauberMC::GetTotXSect() const
1212{
43215add 1213 //total xsection
bb442091 1214 return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
ec852657 1215}
1216
1217//______________________________________________________________________________
1218Double_t AliGlauberMC::GetTotXSectErr() const
1219{
43215add 1220 //total xsection error
bb442091 1221 return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
7288f660 1222 TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
ec852657 1223}
1224
1225//______________________________________________________________________________
bb442091 1226TObjArray *AliGlauberMC::GetNucleons()
ec852657 1227{
43215add 1228 //get array of nucleons
bb442091 1229 if(!fNucleonsA || !fNucleonsB) return 0;
1230 fNucleonsA->SetOwner(0);
1231 fNucleonsB->SetOwner(0);
1232 TObjArray *allnucleons=new TObjArray(fAN+fBN);
1233 allnucleons->SetOwner();
1234 for (Int_t i = 0; i<fAN; i++)
7288f660 1235 {
1236 allnucleons->Add(fNucleonsA->UncheckedAt(i));
1237 }
bb442091 1238 for (Int_t i = 0; i<fBN; i++)
7288f660 1239 {
1240 allnucleons->Add(fNucleonsB->UncheckedAt(i));
1241 }
bb442091 1242 return allnucleons;
1243}
9a5780fa 1244
bb442091 1245//______________________________________________________________________________
1246Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
1247{
43215add 1248 //produce a number distributed acc. neg.bin.dist
bb442091 1249 if(k<=0)
7288f660 1250 {
1251 cout << "Error, zero or negative K" << endl;
1252 return 0;
1253 }
bb442091 1254 return (TMath::Binomial(x+k-1,x))
7288f660 1255 *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
1256 *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
bb442091 1257}
9a5780fa 1258
bb442091 1259//______________________________________________________________________________
43215add 1260Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
bb442091 1261{
43215add 1262 //return random integer from a Negative Binomial Distribution
bb442091 1263 static const Int_t fMaxPlot = 1000;
1264 Double_t array[fMaxPlot];
1265 array[0] = NegativeBinomialDistribution(0,k,nmean);
1266 for (Int_t i=1; i<fMaxPlot; i++)
7288f660 1267 {
1268 array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
1269 }
bb442091 1270 Double_t r = gRandom->Uniform(0,1);
1271 return TMath::BinarySearch(fMaxPlot,array,r)+2;
9a5780fa 1272
bb442091 1273}
7288f660 1274
1275//______________________________________________________________________________
1276Int_t AliGlauberMC::NegativeBinomialRandomSV(Double_t k, Double_t nbar) const
1277{
1278 // negative binomial distribution generator, S. Voloshin, 09-May-2007
1279 Double_t sum=0.;
1280 Int_t i=0;
1281 Double_t ran=gRandom->Rndm();
1282 Double_t trm=1./pow(1.+nbar/k,k);
1283 if (trm==0.)
1284 {
1285 cout<<"NBD overflow"<<" nbar="<<nbar<<" k="<<k<<endl;
1286 return -1;
1287 }
1288 for(i=0; i<2000 && sum<ran ; i++)
1289 {
1290 sum += trm;
1291 trm *= (k+i)/(i+1.)*(nbar/(nbar+k));
1292 }
1293 return i-1;
1294}
1295
bb442091 1296//______________________________________________________________________________
1297Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
7288f660 1298 Double_t nmean,
1299 Int_t k2,
1300 Double_t nmean2,
1301 Double_t alpha ) const
bb442091 1302{
1303 //return random integer from a Double Negative Binomial Distribution
1304 static const Int_t fMaxPlot = 1000;
1305 Double_t array[fMaxPlot];
1306 array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2);
1307 for (Int_t i=1; i<fMaxPlot; i++)
7288f660 1308 {
1309 array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
1310 }
bb442091 1311 Double_t r = gRandom->Uniform(0,1);
1312 return TMath::BinarySearch(fMaxPlot,array,r)+2;
1313}
1314
bb442091 1315//______________________________________________________________________________
7288f660 1316Double_t AliGlauberMC::GetdNdEta() const
bb442091 1317{
7288f660 1318 switch (fMultType)
1319 {
1320 case kSimple:
1321 return GetdNdEtaSimple(fdNdEtaParam);
1322 case kNBD:
1323 return GetdNdEtaNBD(fdNdEtaParam);
1324 case kNBDSV:
1325 return GetdNdEtaNBDSV(fdNdEtaParam);
1326 case kTwoNBD:
1327 return GetdNdEtaTwoNBD(fdNdEtaParam);
1328 case kGBW:
1329 return GetdNdEtaGBW(fdNdEtaParam);
1330 default:
1331 return 0.0;
1332 }
bb442091 1333}
14244e14 1334
bb442091 1335//______________________________________________________________________________
7288f660 1336Double_t AliGlauberMC::GetdNdEtaSimple(const Double_t* p) const
1b4934d5 1337{
1338 //Get particle density per unit of rapidity
bb442091 1339 //using two component model
1340 //Parameters: npp, x
7288f660 1341 Double_t nnp = p[0]; //=8.0
1342 Double_t x = p[1]; //=0.13
1343
bb442091 1344 return nnp*((1.-x)*fNpart/2.+x*fNcoll);
0448b9f3 1345}
1346
1347//______________________________________________________________________________
7288f660 1348Double_t AliGlauberMC::GetdNdEtaGBW( const Double_t* p ) const
0448b9f3 1349{
1350 //Get particle density per unit of rapidity
1351 //using the GBW model
bb442091 1352 //Parameters: delta, lambda, snn
7288f660 1353 Double_t delta = p[0]; //=0.79
1354 Double_t lambda = p[1]; //=0.288
1355 Double_t snn = p[2]; //=30.25
1356
bb442091 1357 return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
7288f660 1358 * TMath::Power(fNpart,(1.-delta)/3./delta);
1b4934d5 1359}
1360
14244e14 1361//_______________________________________________________________________________
7288f660 1362Double_t AliGlauberMC::GetdNdEtaNBDSV ( const Double_t* p ) const
1363{
1364 //Get particle density per unit of rapidity (from Sergei)
1365 Double_t npp = p[0]; //=2.49
1366 Double_t ratioSgm2Mu = p[1]; //=1.7
1367 Double_t xhard = p[2]; //=0.13
1368
1369 double knb=npp/(ratioSgm2Mu-1.); // sgm^2/mu = 1+ mu/k
1370 double scale = (1.-xhard)*fNpart/2.+xhard*fNcoll; // effectively get number of pp collisions
1371 float nch1=-99.;
1372 if (knb*scale <1000.)
1373 {
1374 nch1=(float) NegativeBinomialRandomSV( npp*scale, knb*scale );
1375 }
1376 else
1377 {
1378 nch1=(float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. ) +
1379 (float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. );
1380 }
1381 return nch1;
1382}
1383
1384//_______________________________________________________________________________
1385Double_t AliGlauberMC::GetdNdEtaNBD ( const Double_t* p ) const
bb442091 1386{
1387 //Get particle density per unit of rapidity
1388 //using a aandomized number from a negative binomial distrubution
7288f660 1389 //Parameters: k = related to distribition width=3
1390 // nmean = mean of distribution=4
1391 // beta = set contribution of participants / binary collisions to multiplicity=0.13
1392 Int_t k = TMath::Nint(p[0]);
1393 Double_t nmean = p[1];
1394 Double_t beta = p[2];
1395
bb442091 1396 Double_t mulnp=0.0;
1397 for(Int_t i = 0; i<fNpart; i++)
7288f660 1398 {
1399 mulnp+=NegativeBinomialRandom(k,nmean);
1400 }
bb442091 1401 Double_t mulnb=0.0;
1402 for(Int_t i = 0; i<fNcoll; i++)
7288f660 1403 {
1404 mulnb+=NegativeBinomialRandom(k,nmean);
1405 }
bb442091 1406 return (1-beta)*mulnp/2+beta*mulnb;
1407}
bb442091 1408
14244e14 1409//______________________________________________________________________________
7288f660 1410Double_t AliGlauberMC::GetdNdEtaTwoNBD ( const Double_t* p ) const
bb442091 1411{
1412 //Get particle density per unit of rapidity
1413 //using random numbers from two negative binomial distributions
7288f660 1414 //Parameters: k1 = related to distribition width of distribution 1=3
1415 // nmean1 = mean of distribution 1=4
1416 // k2 = related to distribition width of distribution 2=2
1417 // nmean2 = mean of distribution 2=11
1418 // alpha = set contributions of distrubitin 1 / distribution 2=0.4
1419 // beta = set contribution of participants / binary collisions to multiplicity =0.13
1420 Int_t k1 = TMath::Nint(p[0]);
1421 Double_t nmean1 = p[1];
1422 Int_t k2 = TMath::Nint(p[2]);
1423 Double_t nmean2 = p[3];
1424 Double_t alpha = p[4];
1425 Double_t beta = p[6];
1426
bb442091 1427 Double_t mulnp=0.0;
1428 for(Int_t i = 0; i<fNpart; i++)
1429 {
1430 mulnp+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
1431 }
1432 Double_t mulnb=0.0;
1433 for(Int_t i = 0; i<fNcoll; i++)
1434 {
1435 mulnb+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
1436 }
1437 Double_t mul=(1-beta)*mulnp/2+beta*mulnb;
1438 return mul;
bb442091 1439}
14244e14 1440
1b4934d5 1441//______________________________________________________________________________
1442Double_t AliGlauberMC::GetEccentricityPart() const
1443{
bb442091 1444 //get participant eccentricity of participants
a0278ee7 1445 if (fNpart<2) return 0.0;
2c0ff1e6 1446 return (TMath::Sqrt((fSy2Parts-fSx2Parts)*(fSy2Parts-fSx2Parts)+4*fSxyParts*fSxyParts)/(fSy2Parts+fSx2Parts));
1b4934d5 1447}
1448
106d1ca1 1449//_____________________________________________________________________________
bb442091 1450Double_t AliGlauberMC::GetEccentricityPartColl() const
1451{
1452 //get participant eccentricity of binary collisions
998de9b9 1453 if (fNcoll<0) return 0.0;
8b85b325 1454 if (fSy2Coll==0.0) return 0.0;
bb442091 1455 return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
1456}
1457
4d264dfa 1458//_____________________________________________________________________________
1459Double_t AliGlauberMC::GetEccentricityPartCom() const
1460{
1461 //get participant eccentricity of binary collisions
998de9b9 1462 if (fNcom<0) return 0.0;
4d264dfa 1463 return (TMath::Sqrt((fSy2Com-fSx2Com)*(fSy2Com-fSx2Com)+4*fSxyCom*fSxyCom)/(fSy2Com+fSx2Com));
1464}
1b4934d5 1465//______________________________________________________________________________
1466Double_t AliGlauberMC::GetEccentricity() const
1467{
bb442091 1468 //get standard eccentricity of participants
74d0740e 1469 if (fNpart<2) return 0.0;
2c0ff1e6 1470 return ((fSy2Parts-fSx2Parts)/(fSy2Parts+fSx2Parts));
1471}
1472
1473//______________________________________________________________________________
4d264dfa 1474Double_t AliGlauberMC::GetEccentricityColl() const
2c0ff1e6 1475{
4d264dfa 1476 //get standard eccentricity of binary collisions
998de9b9 1477 if (fNcoll<0) return 0.0;
8b85b325 1478 if (fSy2Coll==0.0) return 0.0;
4d264dfa 1479 return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
1b4934d5 1480}
bb442091 1481
1482//______________________________________________________________________________
4d264dfa 1483Double_t AliGlauberMC::GetEccentricityCom() const
bb442091 1484{
ddf04a47 1485 //get standard eccentricity of combined weight
998de9b9 1486 if (fNcom<0) return 0.0;
1487 return ((fSy2Com-fSx2Com)/(fSy2Com+fSx2Com));
1488 //return (((fMeanY2Com-(fMeanYCom*fMeanYCom)) - (fMeanX2Com-(fMeanXCom*fMeanXCom))) / ((fMeanY2Com-(fMeanYCom*fMeanYCom)) + (fMeanX2Com-(fMeanXCom*fMeanXCom))) );
4d264dfa 1489}
1490
1491//______________________________________________________________________________
1492Double_t AliGlauberMC::GetStoa() const
1493{
1494 //get standard Transverse Overlap Area
1495 if (fNpart<2) return 0.0;
1496 return ( TMath::Pi()*(TMath::Sqrt(fSx2Parts))*(TMath::Sqrt(fSy2Parts)));
bb442091 1497}
14244e14 1498
ec852657 1499//______________________________________________________________________________
1500Bool_t AliGlauberMC::NextEvent(Double_t bgen)
1501{
1b4934d5 1502 //Make a new event
bb442091 1503 Int_t nAttempts = 10; // set indices, max attempts and max comparisons
1504 Bool_t succes = kFALSE;
1505 for(Int_t j=0; j<nAttempts; j++)
1506 {
1507 if(bgen<0||!succes) //get impactparameter
1508 {
1509 bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
1510 }
97d98e48 1511 if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
bb442091 1512 }
1513 return succes;
ec852657 1514}
ec852657 1515//______________________________________________________________________________
106d1ca1 1516Double_t AliGlauberMC::GetEpsilon2Part() const
1517{
1518 //get participant eccentricity of participants
1519 if (fNpart<2) return 0.0;
1520 return (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2);
1521}
106d1ca1 1522//______________________________________________________________________________
106d1ca1 1523Double_t AliGlauberMC::GetEpsilon3Part() const
1524{
1525 //get participant eccentricity of participants
1526 if (fNpart<2) return 0.0;
998de9b9 1527 return (TMath::Sqrt(fMeanr2Cos3Phi*fMeanr2Cos3Phi+fMeanr2Sin3Phi*fMeanr2Sin3Phi)/fMeanr2);
1528 //return (TMath::Sqrt(fMeanr3Cos3Phi*fMeanr3Cos3Phi+fMeanr3Sin3Phi*fMeanr3Sin3Phi)/fMeanr3);
106d1ca1 1529}
106d1ca1 1530//______________________________________________________________________________
106d1ca1 1531Double_t AliGlauberMC::GetEpsilon4Part() const
1532{
1533 //get participant eccentricity of participants
1534 if (fNpart<2) return 0.0;
998de9b9 1535 return (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2);
1536 //return (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4);
106d1ca1 1537}
1538
1539//______________________________________________________________________________
106d1ca1 1540Double_t AliGlauberMC::GetEpsilon5Part() const
1541{
1542 //get participant eccentricity of participants
1543 if (fNpart<2) return 0.0;
998de9b9 1544 return (TMath::Sqrt(fMeanr2Cos5Phi*fMeanr2Cos5Phi+fMeanr2Sin5Phi*fMeanr2Sin5Phi)/fMeanr2);
1545 //return (TMath::Sqrt(fMeanr5Cos5Phi*fMeanr5Cos5Phi+fMeanr5Sin5Phi*fMeanr5Sin5Phi)/fMeanr5);
106d1ca1 1546}
9a5780fa 1547
106d1ca1 1548//______________________________________________________________________________
9a5780fa 1549Double_t AliGlauberMC::GetEpsilon2Coll() const
1550{
1551 //get epsilon2 of binary collisions
998de9b9 1552 if (fNcoll<0) return 0.0;
8b85b325 1553 if (fMeanr2Coll==0.0) return 0.0;
9a5780fa 1554 return (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll);
1555}
1556
1557//______________________________________________________________________________
1558Double_t AliGlauberMC::GetEpsilon3Coll() const
1559{
7288f660 1560 //get epsilon3 of binary collisions
998de9b9 1561 if (fNcoll<0) return 0.0;
8b85b325 1562 if (fMeanr2Coll==0.0) return 0.0;
998de9b9 1563 return (TMath::Sqrt(fMeanr2Cos3PhiColl*fMeanr2Cos3PhiColl+fMeanr2Sin3PhiColl*fMeanr2Sin3PhiColl)/fMeanr2Coll);
1564 //return (TMath::Sqrt(fMeanr3Cos3PhiColl*fMeanr3Cos3PhiColl+fMeanr3Sin3PhiColl*fMeanr3Sin3PhiColl)/fMeanr3Coll);
9a5780fa 1565}
1566
1567//______________________________________________________________________________
1568Double_t AliGlauberMC::GetEpsilon4Coll() const
1569{
7288f660 1570 //get epsilon4 of binary collisions
998de9b9 1571 if (fNcoll<0) return 0.0;
8b85b325 1572 if (fMeanr2Coll==0.0) return 0.0;
998de9b9 1573 return (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll);
1574 //return (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll);
9a5780fa 1575}
106d1ca1 1576
9a5780fa 1577//______________________________________________________________________________
1578Double_t AliGlauberMC::GetEpsilon5Coll() const
1579{
1580 //get epsilon5 of binary collisions
998de9b9 1581 if (fNcoll<0) return 0.0;
8b85b325 1582 if (fMeanr2Coll==0.0) return 0.0;
998de9b9 1583 return (TMath::Sqrt(fMeanr2Cos5PhiColl*fMeanr2Cos5PhiColl+fMeanr2Sin5PhiColl*fMeanr2Sin5PhiColl)/fMeanr2Coll);
1584 //return (TMath::Sqrt(fMeanr5Cos5PhiColl*fMeanr5Cos5PhiColl+fMeanr5Sin5PhiColl*fMeanr5Sin5PhiColl)/fMeanr5Coll);
4d264dfa 1585}
1586
1587//______________________________________________________________________________
1588Double_t AliGlauberMC::GetEpsilon2Com() const
1589{
1590 //get epsilon2 of binary collisions
998de9b9 1591 if (fNcom<0) return 0.0;
4d264dfa 1592 return (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com);
1593}
1594
1595//______________________________________________________________________________
1596Double_t AliGlauberMC::GetEpsilon3Com() const
1597{
1598 //get epsilon3 of binary collisions
998de9b9 1599 if (fNcom<0) return 0.0;
1600 return (TMath::Sqrt(fMeanr2Cos3PhiCom*fMeanr2Cos3PhiCom+fMeanr2Sin3PhiCom*fMeanr2Sin3PhiCom)/fMeanr2Com);
1601 //return (TMath::Sqrt(fMeanr3Cos3PhiCom*fMeanr3Cos3PhiCom+fMeanr3Sin3PhiCom*fMeanr3Sin3PhiCom)/fMeanr3Com);
4d264dfa 1602}
1603
1604//______________________________________________________________________________
1605Double_t AliGlauberMC::GetEpsilon4Com() const
1606{
1607 //get epsilon4 of binary collisions
998de9b9 1608 if (fNcom<0) return 0.0;
1609 return (TMath::Sqrt(fMeanr2Cos4PhiCom*fMeanr2Cos4PhiCom+fMeanr2Sin4PhiCom*fMeanr2Sin4PhiCom)/fMeanr2Com);
1610 //return (TMath::Sqrt(fMeanr4Cos4PhiCom*fMeanr4Cos4PhiCom+fMeanr4Sin4PhiCom*fMeanr4Sin4PhiCom)/fMeanr4Com);
4d264dfa 1611}
1612
1613//______________________________________________________________________________
1614Double_t AliGlauberMC::GetEpsilon5Com() const
1615{
1616 //get epsilon5 of binary collisions
998de9b9 1617 if (fNcom<0) return 0.0;
1618 return (TMath::Sqrt(fMeanr2Cos5PhiCom*fMeanr2Cos5PhiCom+fMeanr2Sin5PhiCom*fMeanr2Sin5PhiCom)/fMeanr2Com);
1619 //return (TMath::Sqrt(fMeanr5Cos5PhiCom*fMeanr5Cos5PhiCom+fMeanr5Sin5PhiCom*fMeanr5Sin5PhiCom)/fMeanr5Com);
9a5780fa 1620}
1621
8f9ee6e7 1622//______________________________________________________________________________
1623Double_t AliGlauberMC::GetPsi2() const
1624{
d10259b9 1625 return ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2);
8f9ee6e7 1626}
1627
1628//______________________________________________________________________________
1629Double_t AliGlauberMC::GetPsi3() const
1630{
998de9b9 1631 return ((TMath::ATan2(fMeanr2Sin3Phi,fMeanr2Cos3Phi)+TMath::Pi())/3);
1632 //return ((TMath::ATan2(fMeanr3Sin3Phi,fMeanr3Cos3Phi)+TMath::Pi())/3);
8f9ee6e7 1633}
1634
1635//______________________________________________________________________________
1636Double_t AliGlauberMC::GetPsi4() const
1637{
998de9b9 1638 return ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4);
1639 //return ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4);
8f9ee6e7 1640}
1641
1642//______________________________________________________________________________
1643Double_t AliGlauberMC::GetPsi5() const
1644{
998de9b9 1645 return ((TMath::ATan2(fMeanr2Sin5Phi,fMeanr2Cos5Phi)+TMath::Pi())/5);
1646 //return ((TMath::ATan2(fMeanr5Sin5Phi,fMeanr5Cos5Phi)+TMath::Pi())/5);
8f9ee6e7 1647}
1648
4d264dfa 1649/*_____________________________________________________________________________
1650Double_t AliGlauberMC::GetE43Part() const
1651{
1652 //get participant eccentricity of participants
1653 if (fNpart<2) return 0.0;
1654 //return ((TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2) * TMath::Cos (((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1655 return ((TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1656}
1657
1658
1659//______________________________________________________________________________
1660Double_t AliGlauberMC::GetE43Coll() const
1661{
1662 //get epsilon5 of binary collisions
1663 if (fNcoll<2) return 0.0;
1664 //return ((TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1665 return ((TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1666}
1667
1668//______________________________________________________________________________
1669Double_t AliGlauberMC::GetE43Com() const
1670{
1671 //get epsilon5 of binary collisions
1672 if (fNcom<2) return 0.0;
1673 //return ((TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos4PhiCom*fMeanr2Cos4PhiCom+fMeanr2Sin4PhiCom*fMeanr2Sin4PhiCom)/fMeanr2Com) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1674 return ((TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr4Cos4PhiCom*fMeanr4Cos4PhiCom+fMeanr4Sin4PhiCom*fMeanr4Sin4PhiCom)/fMeanr4Com) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1675}
1676//___________________________________________________________________________
1677
1678Double_t AliGlauberMC::GetPsi4m2() const
1679{
1680 //return (TMath::Cos(4*((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4) - ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2)));
1681 return (TMath::Cos(4*(((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)-((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2))));
1682}
1683*/
9a5780fa 1684//______________________________________________________________________________
ec852657 1685void AliGlauberMC::Run(Int_t nevents)
1686{
43215add 1687 //example run
bb442091 1688 cout << "Generating " << nevents << " events..." << endl;
1689 TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
1690 TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
1691 if (fnt == 0)
1692 {
1693 fnt = new TNtuple(name,title,
95916d44 1694 "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:Stoa:VarEColl:VarECom:VarEPart:VarEPartColl:VarEPartCom:dNdEta:dNdEtaGBW:dNdEtaTwoNBD:xsect:tAA:Epsl2:Epsl3:Epsl4:Epsl5:E2Coll:E3Coll:E4Coll:E5Coll:E2Com:E3Com:E4Com:E5Com:Psi2:Psi3:Psi4:Psi5:BNN:signn:Ncollw");
bb442091 1695 fnt->SetDirectory(0);
1696 }
1697 Int_t q = 0;
1698 Int_t u = 0;
1699 for (Int_t i = 0; i<nevents; i++)
7288f660 1700 {
1701
1702 if(!NextEvent())
1703 {
1704 u++;
1705 continue;
1706 }
1707
1708 q++;
95916d44 1709 Float_t v[48];
7288f660 1710 v[0] = GetNpart();
1711 v[1] = GetNcoll();
1712 v[2] = fBMC;
1713 v[3] = fMeanXParts;
1714 v[4] = fMeanYParts;
2c0ff1e6 1715 v[5] = fMeanX2Parts;
1716 v[6] = fMeanY2Parts;
1717 v[7] = fMeanXYParts;
1718 v[8] = fSx2Parts;
1719 v[9] = fSy2Parts;
1720 v[10] = fSxyParts;
7288f660 1721 v[11] = fMeanXSystem;
1722 v[12] = fMeanYSystem;
1723 v[13] = fMeanXA;
1724 v[14] = fMeanYA;
1725 v[15] = fMeanXB;
1726 v[16] = fMeanYB;
1727 v[17] = GetEccentricity();
2c0ff1e6 1728 v[18] = GetStoa();
1729 v[19] = GetEccentricityColl();
4d264dfa 1730 v[20] = GetEccentricityCom();
1731 v[21] = GetEccentricityPart();
1732 v[22] = GetEccentricityPartColl();
1733 v[23] = GetEccentricityPartCom();
7288f660 1734 if (fDoPartProd)
1735 {
4d264dfa 1736 v[24] = GetdNdEta();
1737 v[25] = GetdNdEta();
1738 v[26] = v[24]+v[25];
7288f660 1739 }
1740 else
bb442091 1741 {
2c0ff1e6 1742 v[24] = 0;
4d264dfa 1743 v[25] = 0;
1744 v[26] = 0;
106d1ca1 1745 }
4d264dfa 1746 v[27]=fXSect;
7288f660 1747
1748 Float_t mytAA=-999;
1749 if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
4d264dfa 1750 v[28]=mytAA;
7288f660 1751 //_____________epsilon2,3,4,4_______
4d264dfa 1752 v[29] = GetEpsilon2Part();
1753 v[30] = GetEpsilon3Part();
1754 v[31] = GetEpsilon4Part();
1755 v[32] = GetEpsilon5Part();
1756 v[33] = GetEpsilon2Coll();
1757 v[34] = GetEpsilon3Coll();
1758 v[35] = GetEpsilon4Coll();
1759 v[36] = GetEpsilon5Coll();
1760 v[37] = GetEpsilon2Com();
1761 v[38] = GetEpsilon3Com();
1762 v[39] = GetEpsilon4Com();
1763 v[40] = GetEpsilon5Com();
1764 v[41] = GetPsi2();
1765 v[42] = GetPsi3();
1766 v[43] = GetPsi4();
1767 v[44] = GetPsi5();
22e47719 1768 v[45] = fBNN;
1769 v[46] = fXSect;
95916d44
CL		 1770 v[47] = fNcollw;
1771
7288f660 1772 //always at the end
1773 fnt->Fill(v);
1774
1775 if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
1776 }
bb442091 1777 std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events: " << q << " discarded events: " << u <<"."<< endl;
ec852657 1778}
1779
1780//---------------------------------------------------------------------------------
43215add 1781void AliGlauberMC::RunAndSaveNtuple( Int_t n,
1782 const Option_t *sysA,
1783 const Option_t *sysB,
ec852657 1784 Double_t signn,
1785 Double_t mind,
7288f660 1786 Double_t r,
1787 Double_t a,
ec852657 1788 const char *fname)
1789{
43215add 1790 //example run
566fc3d0 1791 AliGlauberMC mcg(sysA,sysB,signn);
1792 mcg.SetMinDistance(mind);
14244e14 1793 mcg.Setr(r);
1794 mcg.Seta(a);
566fc3d0 1795 mcg.Run(n);
1796 TNtuple *nt=mcg.GetNtuple();
bb442091 1797 TFile out(fname,"recreate",fname,9);
1798 if(nt) nt->Write();
14244e14 1799 printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
bb442091 1800 out.Close();
ec852657 1801}
1802
1803//---------------------------------------------------------------------------------
43215add 1804void AliGlauberMC::RunAndSaveNucleons( Int_t n,
1805 const Option_t *sysA,
bb442091 1806 Option_t *sysB,
43215add 1807 const Double_t signn,
ec852657 1808 Double_t mind,
1809 Bool_t verbose,
1810 const char *fname)
1811{
43215add 1812 //example run
566fc3d0 1813 AliGlauberMC mcg(sysA,sysB,signn);
1814 mcg.SetMinDistance(mind);
bb442091 1815 TFile *out=0;
1816 if(fname)
1817 out=new TFile(fname,"recreate",fname,9);
1818
1819 for(Int_t ievent=0; ievent<n; ievent++)
7288f660 1820 {
1821
1822 //get an event with at least one collision
1823 while(!mcg.NextEvent()) {}
1824
1825 //access, save and (if wanted) print out nucleons
1826 TObjArray* nucleons=mcg.GetNucleons();
1827 if(!nucleons) continue;
1828 if(out)
1829 nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
1830
1831 if(verbose)
bb442091 1832 {
7288f660 1833 cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
1834 cout<<"B = "<<mcg.GetB()<<" Npart = "<<mcg.GetNpart()<<endl<<endl;
1835 printf("Nucleus\t X\t Y\t Z\tNcoll\n");
1836 Int_t nNucls=nucleons->GetEntries();
1837 for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
1838 {
1839 AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
1840 Char_t nucleus='A';
1841 if(nucl->IsInNucleusB()) nucleus='B';
1842 Double_t x=nucl->GetX();
1843 Double_t y=nucl->GetY();
1844 Double_t z=nucl->GetZ();
1845 Int_t ncoll=nucl->GetNColl();
1846 printf(" %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
1847 }
bb442091 1848 }
7288f660 1849 }
bb442091 1850 if(out) delete out;
ec852657 1851}
1852
566fc3d0 1853//---------------------------------------------------------------------------------
1854void AliGlauberMC::Reset()
1855{
1856 //delete the ntuple
7288f660 1857 delete fnt;
1858 fnt=NULL;
566fc3d0 1859}