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