1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ////////////////////////////////////////////////////////////////////////////////
18 // Glauber MC implementation
20 // origin: PHOBOS experiment
21 // alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
22 // update: You Zhou, Nikhef, yzhou@nikhef.nl
23 ////////////////////////////////////////////////////////////////////////////////
25 #include <Riostream.h>
30 #include <TObjArray.h>
38 #include "AliGlauberNucleon.h"
39 #include "AliGlauberNucleus.h"
40 #include "AliGlauberMC.h"
42 ClassImp(AliGlauberMC)
44 //______________________________________________________________________________
45 AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
101 fMeanr2Cos2PhiColl(0),
102 fMeanr2Sin2PhiColl(0),
103 fMeanr2Cos3PhiColl(0),
104 fMeanr2Sin3PhiColl(0),
105 fMeanr2Cos4PhiColl(0),
106 fMeanr2Sin4PhiColl(0),
107 fMeanr2Cos5PhiColl(0),
108 fMeanr2Sin5PhiColl(0),
109 fMeanr3Cos3PhiColl(0),
110 fMeanr3Sin3PhiColl(0),
111 fMeanr4Cos4PhiColl(0),
112 fMeanr4Sin4PhiColl(0),
113 fMeanr5Cos5PhiColl(0),
114 fMeanr5Sin5PhiColl(0),
126 fdNdEtaParam[0] = 8.0;
127 fdNdEtaParam[1] = 0.13;
128 fdNdEtaGBWParam[0] = 0.79;
129 fdNdEtaGBWParam[1] = 0.288;
130 fdNdEtaGBWParam[2] = 30.25;
131 fdNdEtaNBDParam[0] = 3.0;
132 fdNdEtaNBDParam[1] = 4.0;
133 fdNdEtaNBDParam[2] = 0.13;
134 fdNdEtaTwoNBDParam[0] = 3.0;
135 fdNdEtaTwoNBDParam[1] = 4.0;
136 fdNdEtaTwoNBDParam[2] = 2.0;
137 fdNdEtaTwoNBDParam[3] = 11.0;
138 fdNdEtaTwoNBDParam[4] = 0.4;
139 fdNdEtaTwoNBDParam[5] = 0.13;
141 SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
142 SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
145 //______________________________________________________________________________
146 AliGlauberMC::~AliGlauberMC()
152 //______________________________________________________________________________
153 AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
155 fANucleus(in.fANucleus),
156 fBNucleus(in.fBNucleus),
158 fNucleonsA(in.fNucleonsA),
159 fNucleonsB(in.fNucleonsB),
166 fMeanXParts(in.fMeanXParts),
167 fMeanYParts(in.fMeanYParts),
168 fMeanXColl(in.fMeanXColl),
169 fMeanYColl(in.fMeanYColl),
170 fMeanX2Coll(in.fMeanX2Coll),
171 fMeanY2Coll(in.fMeanY2Coll),
172 fMeanXYColl(in.fMeanXYColl),
173 fMeanXSystem(in.fMeanXSystem),
174 fMeanYSystem(in.fMeanYSystem),
181 fTotalEvents(in.fTotalEvents),
184 fMaxNpartFound(in.fMaxNpartFound),
191 fMeanr2Cos2Phi(in.fMeanr2Cos2Phi),
192 fMeanr2Sin2Phi(in.fMeanr2Sin2Phi),
193 fMeanr2Cos3Phi(in.fMeanr2Cos3Phi),
194 fMeanr2Sin3Phi(in.fMeanr2Sin3Phi),
195 fMeanr2Cos4Phi(in.fMeanr2Cos4Phi),
196 fMeanr2Sin4Phi(in.fMeanr2Sin4Phi),
197 fMeanr2Cos5Phi(in.fMeanr2Cos5Phi),
198 fMeanr2Sin5Phi(in.fMeanr2Sin5Phi),
199 fMeanr3Cos3Phi(in.fMeanr3Cos3Phi),
200 fMeanr3Sin3Phi(in.fMeanr3Sin3Phi),
201 fMeanr4Cos4Phi(in.fMeanr4Cos4Phi),
202 fMeanr4Sin4Phi(in.fMeanr4Sin4Phi),
203 fMeanr5Cos5Phi(in.fMeanr5Cos5Phi),
204 fMeanr5Sin5Phi(in.fMeanr5Sin5Phi),
205 fMeanr2Coll(in.fMeanr2Coll),
206 fMeanr3Coll(in.fMeanr3Coll),
207 fMeanr4Coll(in.fMeanr4Coll),
208 fMeanr5Coll(in.fMeanr5Coll),
209 fMeanr2Cos2PhiColl(in.fMeanr2Cos2PhiColl),
210 fMeanr2Sin2PhiColl(in.fMeanr2Sin2PhiColl),
211 fMeanr2Cos3PhiColl(in.fMeanr2Cos3PhiColl),
212 fMeanr2Sin3PhiColl(in.fMeanr2Sin3PhiColl),
213 fMeanr2Cos4PhiColl(in.fMeanr2Cos4PhiColl),
214 fMeanr2Sin4PhiColl(in.fMeanr2Sin4PhiColl),
215 fMeanr2Cos5PhiColl(in.fMeanr2Cos5PhiColl),
216 fMeanr2Sin5PhiColl(in.fMeanr2Sin5PhiColl),
217 fMeanr3Cos3PhiColl(in.fMeanr3Cos3PhiColl),
218 fMeanr3Sin3PhiColl(in.fMeanr3Sin3PhiColl),
219 fMeanr4Cos4PhiColl(in.fMeanr4Cos4PhiColl),
220 fMeanr4Sin4PhiColl(in.fMeanr4Sin4PhiColl),
221 fMeanr5Cos5PhiColl(in.fMeanr5Cos5PhiColl),
222 fMeanr5Sin5PhiColl(in.fMeanr5Sin5PhiColl),
226 fSx2Coll(in.fSx2Coll),
227 fSy2Coll(in.fSy2Coll),
228 fSxyColl(in.fSxyColl),
234 memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
235 memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
236 memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
237 memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
240 //______________________________________________________________________________
241 AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
244 fANucleus=in.fANucleus;
245 fBNucleus=in.fBNucleus;
247 fNucleonsA=in.fNucleonsA;
248 fNucleonsB=in.fNucleonsB;
259 fMeanXParts=in.fMeanXParts;
260 fMeanYParts=in.fMeanYParts;
261 fMeanXColl=in.fMeanXColl;
262 fMeanYColl=in.fMeanYColl;
263 fMeanX2Coll=in.fMeanX2Coll;
264 fMeanY2Coll=in.fMeanY2Coll;
265 fMeanXYColl=in.fMeanXYColl;
266 fMeanr2Coll=in.fMeanr2Coll;
267 fMeanr3Coll=in.fMeanr3Coll;
268 fMeanr4Coll=in.fMeanr4Coll;
269 fMeanr5Coll=in.fMeanr5Coll;
270 fMeanXSystem=in.fMeanXSystem;
271 fMeanYSystem=in.fMeanYSystem;
278 fTotalEvents=in.fTotalEvents;
281 memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
282 memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
283 memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
284 memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
285 fMaxNpartFound=in.fMaxNpartFound;
288 fMeanr2Cos2Phi=in.fMeanr2Cos2Phi;
289 fMeanr2Sin2Phi=in.fMeanr2Sin2Phi;
290 fMeanr2Cos3Phi=in.fMeanr2Cos3Phi;
291 fMeanr2Sin3Phi=in.fMeanr2Sin3Phi;
292 fMeanr2Cos4Phi=in.fMeanr2Cos4Phi;
293 fMeanr2Sin4Phi=in.fMeanr2Sin4Phi;
294 fMeanr2Cos5Phi=in.fMeanr2Cos5Phi;
295 fMeanr2Sin5Phi=in.fMeanr2Sin5Phi;
296 fMeanr3Cos3Phi=in.fMeanr3Cos3Phi;
297 fMeanr3Sin3Phi=in.fMeanr3Sin3Phi;
298 fMeanr4Cos4Phi=in.fMeanr4Cos4Phi;
299 fMeanr4Sin4Phi=in.fMeanr4Sin4Phi;
300 fMeanr5Cos5Phi=in.fMeanr5Cos5Phi;
301 fMeanr5Sin5Phi=in.fMeanr5Sin5Phi;
302 fMeanr2Cos2PhiColl=in.fMeanr2Cos2PhiColl;
303 fMeanr2Sin2PhiColl=in.fMeanr2Sin2PhiColl;
304 fMeanr2Cos3PhiColl=in.fMeanr2Cos3PhiColl;
305 fMeanr2Sin3PhiColl=in.fMeanr2Sin3PhiColl;
306 fMeanr2Cos4PhiColl=in.fMeanr2Cos4PhiColl;
307 fMeanr2Sin4PhiColl=in.fMeanr2Sin4PhiColl;
308 fMeanr2Cos5PhiColl=in.fMeanr2Cos5PhiColl;
309 fMeanr2Sin5PhiColl=in.fMeanr2Sin5PhiColl;
310 fMeanr3Cos3PhiColl=in.fMeanr3Cos3PhiColl;
311 fMeanr3Sin3PhiColl=in.fMeanr3Sin3PhiColl;
312 fMeanr4Cos4PhiColl=in.fMeanr4Cos4PhiColl;
313 fMeanr4Sin4PhiColl=in.fMeanr4Sin4PhiColl;
314 fMeanr5Cos5PhiColl=in.fMeanr5Cos5PhiColl;
315 fMeanr5Sin5PhiColl=in.fMeanr5Sin5PhiColl;
319 fSx2Coll=in.fSx2Coll;
320 fSy2Coll=in.fSy2Coll;
321 fSxyColl=in.fSxyColl;
327 //______________________________________________________________________________
328 Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
331 fANucleus.ThrowNucleons(-bgen/2.);
332 fNucleonsA = fANucleus.GetNucleons();
333 fAN = fANucleus.GetN();
334 for (Int_t i = 0; i<fAN; i++)
336 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
337 nucleonA->SetInNucleusA();
339 fBNucleus.ThrowNucleons(bgen/2.);
340 fNucleonsB = fBNucleus.GetNucleons();
341 fBN = fBNucleus.GetN();
342 for (Int_t i = 0; i<fBN; i++)
344 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
345 nucleonB->SetInNucleusB();
348 // "ball" diameter = distance at which two balls interact
349 Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
351 // for each of the A nucleons in nucleus B
352 for (Int_t i = 0; i<fBN; i++)
354 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
355 for (Int_t j = 0 ; j < fAN ; j++)
357 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
358 Double_t dx = nucleonB->GetX()-nucleonA->GetX();
359 Double_t dy = nucleonB->GetY()-nucleonA->GetY();
360 Double_t dij = dx*dx+dy*dy;
369 return CalcResults(bgen);
372 //______________________________________________________________________________
373 Bool_t AliGlauberMC::CalcResults(Double_t bgen)
375 // calc results for the given event
376 //return true if we have participants
418 fMeanr2Cos2PhiColl=0.;
419 fMeanr2Sin2PhiColl=0.;
420 fMeanr2Cos3PhiColl=0.;
421 fMeanr2Sin3PhiColl=0.;
422 fMeanr2Cos4PhiColl=0.;
423 fMeanr2Sin4PhiColl=0.;
424 fMeanr2Cos5PhiColl=0.;
425 fMeanr2Sin5PhiColl=0.;
426 fMeanr3Cos3PhiColl=0.;
427 fMeanr3Sin3PhiColl=0.;
428 fMeanr4Cos4PhiColl=0.;
429 fMeanr4Sin4PhiColl=0.;
430 fMeanr5Cos5PhiColl=0.;
431 fMeanr5Sin5PhiColl=0.;
433 for (Int_t i = 0; i<fAN; i++)
435 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
436 Double_t xA = nucleonA->GetX();
437 Double_t yA = nucleonA->GetY();
438 Double_t r2A = xA*xA+yA*yA;
439 Double_t rA = TMath::Sqrt(r2A);
440 Double_t r3A = r2A*rA;
441 Double_t r4A = r3A*rA;
442 Double_t r5A = r4A*rA;
443 Double_t phiA = TMath::ATan2(yA,xA);
444 Double_t sin2PhiA = TMath::Sin(2.*phiA);
445 Double_t cos2PhiA = TMath::Cos(2.*phiA);
446 Double_t sin3PhiA = TMath::Sin(3.*phiA);
447 Double_t cos3PhiA = TMath::Cos(3.*phiA);
448 Double_t sin4PhiA = TMath::Sin(4.*phiA);
449 Double_t cos4PhiA = TMath::Cos(4.*phiA);
450 Double_t sin5PhiA = TMath::Sin(5.*phiA);
451 Double_t cos5PhiA = TMath::Cos(5.*phiA);
452 //printf("x: %.3f, y:%.3f, r:%.3f, phi:%.3f, sp:%.6f,cp:%.6f\n",xA,yA,TMath::Sqrt(r2A),PhiA,Sin2PhiA, Cos2PhiA);
458 if(nucleonA->IsWounded())
470 fMeanr2Cos2Phi += r2A*cos2PhiA;
471 fMeanr2Sin2Phi += r2A*sin2PhiA;
472 fMeanr2Cos3Phi += r2A*cos3PhiA;
473 fMeanr2Sin3Phi += r2A*sin3PhiA;
474 fMeanr2Cos4Phi += r2A*cos4PhiA;
475 fMeanr2Sin4Phi += r2A*sin4PhiA;
476 fMeanr2Cos5Phi += r2A*cos5PhiA;
477 fMeanr2Sin5Phi += r2A*sin5PhiA;
478 fMeanr3Cos3Phi += r3A*cos3PhiA;
479 fMeanr3Sin3Phi += r3A*sin3PhiA;
480 fMeanr4Cos4Phi += r4A*cos4PhiA;
481 fMeanr4Sin4Phi += r4A*sin4PhiA;
482 fMeanr5Cos5Phi += r5A*cos5PhiA;
483 fMeanr5Sin5Phi += r5A*sin5PhiA;
487 for (Int_t i = 0; i<fBN; i++)
489 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
490 Double_t xB=nucleonB->GetX();
491 Double_t yB=nucleonB->GetY();
492 Double_t r2B=xB*xB+yB*yB;
493 Double_t rB = TMath::Sqrt(r2B);
494 Double_t r3B = r2B*rB;
495 Double_t r4B = r3B*rB;
496 Double_t r5B = r4B*rB;
497 Double_t phiB = TMath::ATan2(yB,xB);
498 Double_t sin2PhiB = TMath::Sin(2*phiB);
499 Double_t cos2PhiB = TMath::Cos(2*phiB);
500 Double_t sin3PhiB = TMath::Sin(3*phiB);
501 Double_t cos3PhiB = TMath::Cos(3*phiB);
502 Double_t sin4PhiB = TMath::Sin(4*phiB);
503 Double_t cos4PhiB = TMath::Cos(4*phiB);
504 Double_t sin5PhiB = TMath::Sin(5*phiB);
505 Double_t cos5PhiB = TMath::Cos(5*phiB);
511 if(nucleonB->IsWounded())
513 Int_t ncoll = nucleonB->GetNColl();
516 fMeanXColl += xB*ncoll;
518 fMeanYColl += yB*ncoll;
520 fMeanX2Coll += xB*xB*ncoll*ncoll;
522 fMeanY2Coll += yB*yB*ncoll*ncoll;
524 fMeanXYColl += xB*yB*ncoll*ncoll;
525 fNcoll += nucleonB->GetNColl();
530 fMeanr2Cos2Phi += r2B*cos2PhiB;
531 fMeanr2Sin2Phi += r2B*sin2PhiB;
532 fMeanr2Cos3Phi += r2B*cos3PhiB;
533 fMeanr2Sin3Phi += r2B*sin3PhiB;
534 fMeanr2Cos4Phi += r2B*cos4PhiB;
535 fMeanr2Sin4Phi += r2B*sin4PhiB;
536 fMeanr2Cos5Phi += r2B*cos5PhiB;
537 fMeanr2Sin5Phi += r2B*sin5PhiB;
538 fMeanr3Cos3Phi += r3B*cos3PhiB;
539 fMeanr3Sin3Phi += r3B*sin3PhiB;
540 fMeanr4Cos4Phi += r4B*cos4PhiB;
541 fMeanr4Sin4Phi += r4B*sin4PhiB;
542 fMeanr5Cos5Phi += r5B*cos5PhiB;
543 fMeanr5Sin5Phi += r5B*sin5PhiB;
544 fMeanr2Coll += r2B*ncoll*ncoll;
545 fMeanr3Coll += r3B*ncoll*ncoll;
546 fMeanr4Coll += r4B*ncoll*ncoll;
547 fMeanr5Coll += r5B*ncoll*ncoll;
548 fMeanr2Cos2PhiColl += r2B*cos2PhiB*ncoll*ncoll;
549 fMeanr2Sin2PhiColl += r2B*sin2PhiB*ncoll*ncoll;
550 fMeanr2Cos3PhiColl += r2B*cos3PhiB*ncoll*ncoll;
551 fMeanr2Sin3PhiColl += r2B*sin3PhiB*ncoll*ncoll;
552 fMeanr2Cos4PhiColl += r2B*cos4PhiB*ncoll*ncoll;
553 fMeanr2Sin4PhiColl += r2B*sin4PhiB*ncoll*ncoll;
554 fMeanr2Cos5PhiColl += r2B*cos5PhiB*ncoll*ncoll;
555 fMeanr2Sin5PhiColl += r2B*sin5PhiB*ncoll*ncoll;
556 fMeanr3Cos3PhiColl += r3B*cos3PhiB*ncoll*ncoll;
557 fMeanr3Sin3PhiColl += r3B*sin3PhiB*ncoll*ncoll;
558 fMeanr4Cos4PhiColl += r4B*cos4PhiB*ncoll*ncoll;
559 fMeanr4Sin4PhiColl += r4B*sin4PhiB*ncoll*ncoll;
560 fMeanr5Cos5PhiColl += r5B*cos5PhiB*ncoll*ncoll;
561 fMeanr5Sin5PhiColl += r5B*sin5PhiB*ncoll*ncoll;
567 fMeanXParts /= fNpart;
568 fMeanYParts /= fNpart;
576 fMeanr2Cos2Phi /= fNpart;
577 fMeanr2Sin2Phi /= fNpart;
578 fMeanr2Cos3Phi /= fNpart;
579 fMeanr2Sin3Phi /= fNpart;
580 fMeanr2Cos4Phi /= fNpart;
581 fMeanr2Sin4Phi /= fNpart;
582 fMeanr2Cos5Phi /= fNpart;
583 fMeanr2Sin5Phi /= fNpart;
584 fMeanr3Cos3Phi /= fNpart;
585 fMeanr3Sin3Phi /= fNpart;
586 fMeanr4Cos4Phi /= fNpart;
587 fMeanr4Sin4Phi /= fNpart;
588 fMeanr5Cos5Phi /= fNpart;
589 fMeanr5Sin5Phi /= fNpart;
620 fMeanXColl /= fNcoll;
621 fMeanYColl /= fNcoll;
622 fMeanX2Coll /= fNcoll;
623 fMeanY2Coll /= fNcoll;
624 fMeanXYColl /= fNcoll;
625 fMeanr2Coll /= fNcoll;
626 fMeanr3Coll /= fNcoll;
627 fMeanr4Coll /= fNcoll;
628 fMeanr5Coll /= fNcoll;
629 fMeanr2Cos2PhiColl /= fNcoll;
630 fMeanr2Sin2PhiColl /= fNcoll;
631 fMeanr2Cos3PhiColl /= fNcoll;
632 fMeanr2Sin3PhiColl /= fNcoll;
633 fMeanr2Cos4PhiColl /= fNcoll;
634 fMeanr2Sin4PhiColl /= fNcoll;
635 fMeanr2Cos5PhiColl /= fNcoll;
636 fMeanr2Sin5PhiColl /= fNcoll;
637 fMeanr3Cos3PhiColl /= fNcoll;
638 fMeanr3Sin3PhiColl /= fNcoll;
639 fMeanr4Cos4PhiColl /= fNcoll;
640 fMeanr4Sin4PhiColl /= fNcoll;
641 fMeanr5Cos5PhiColl /= fNcoll;
642 fMeanr5Sin5PhiColl /= fNcoll;
651 fMeanr2Cos2PhiColl =0;
652 fMeanr2Sin2PhiColl =0;
653 fMeanr2Cos3PhiColl =0;
654 fMeanr2Sin3PhiColl =0;
655 fMeanr2Cos4PhiColl =0;
656 fMeanr2Sin4PhiColl =0;
657 fMeanr2Cos5PhiColl =0;
658 fMeanr2Sin5PhiColl =0;
659 fMeanr3Cos3PhiColl =0;
660 fMeanr3Sin3PhiColl =0;
661 fMeanr4Cos4PhiColl =0;
662 fMeanr4Sin4PhiColl =0;
663 fMeanr5Cos5PhiColl =0;
664 fMeanr5Sin5PhiColl =0;
668 fMeanXSystem /= (fAN + fBN);
669 fMeanYSystem /= (fAN + fBN);
698 fSx2=fMeanX2-(fMeanXParts*fMeanXParts);
699 fSy2=fMeanY2-(fMeanYParts*fMeanYParts);
700 fSxy=fMeanXY-fMeanXParts*fMeanYParts;
701 fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl);
702 fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl);
703 fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl;
706 if (fNpart>0) fEvents++;
707 if (fNpart==0) return kFALSE;
708 if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
713 //______________________________________________________________________________
714 void AliGlauberMC::Draw(Option_t* /*option*/)
717 fANucleus.Draw(fXSect, 2);
718 fBNucleus.Draw(fXSect, 4);
725 e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0);
726 e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0);
729 //______________________________________________________________________________
730 Double_t AliGlauberMC::GetTotXSect() const
733 return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
736 //______________________________________________________________________________
737 Double_t AliGlauberMC::GetTotXSectErr() const
739 //total xsection error
740 return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
741 TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
744 //______________________________________________________________________________
745 TObjArray *AliGlauberMC::GetNucleons()
747 //get array of nucleons
748 if(!fNucleonsA || !fNucleonsB) return 0;
749 fNucleonsA->SetOwner(0);
750 fNucleonsB->SetOwner(0);
751 TObjArray *allnucleons=new TObjArray(fAN+fBN);
752 allnucleons->SetOwner();
753 for (Int_t i = 0; i<fAN; i++)
755 allnucleons->Add(fNucleonsA->UncheckedAt(i));
757 for (Int_t i = 0; i<fBN; i++)
759 allnucleons->Add(fNucleonsB->UncheckedAt(i));
764 //______________________________________________________________________________
765 Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
767 //produce a number distributed acc. neg.bin.dist
770 cout << "Error, zero or negative K" << endl;
773 return (TMath::Binomial(x+k-1,x))
774 *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
775 *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
778 //______________________________________________________________________________
779 Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
781 //return random integer from a Negative Binomial Distribution
782 static const Int_t fMaxPlot = 1000;
783 Double_t array[fMaxPlot];
784 array[0] = NegativeBinomialDistribution(0,k,nmean);
785 for (Int_t i=1; i<fMaxPlot; i++)
787 array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
789 Double_t r = gRandom->Uniform(0,1);
790 return TMath::BinarySearch(fMaxPlot,array,r)+2;
793 //______________________________________________________________________________
794 Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
798 Double_t alpha ) const
800 //return random integer from a Double Negative Binomial Distribution
801 static const Int_t fMaxPlot = 1000;
802 Double_t array[fMaxPlot];
803 array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2);
804 for (Int_t i=1; i<fMaxPlot; i++)
806 array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
808 Double_t r = gRandom->Uniform(0,1);
809 return TMath::BinarySearch(fMaxPlot,array,r)+2;
812 //______________________________________________________________________________
813 void AliGlauberMC::SetdNdEtaParam(Double_t nnp, Double_t x)
815 // set parameters used for calculating multiplicity, see GetdNdEta() for commments
820 //______________________________________________________________________________
821 void AliGlauberMC::SetdNdEtaGBWParam(Double_t delta, Double_t lambda, Double_t snn)
823 // set parameters used for calculating multiplicity see GetdNdEtaGBW() for commments
824 fdNdEtaGBWParam[0]=delta;
825 fdNdEtaGBWParam[1]=lambda;
826 fdNdEtaGBWParam[2]=snn;
829 //______________________________________________________________________________
830 void AliGlauberMC::SetdNdEtaNBDParam(Double_t k, Double_t nmean, Double_t beta)
832 // set parameters used for calculating multiplicity see GetdNdEtaNBD() for commments
833 fdNdEtaNBDParam[0]=k;
834 fdNdEtaNBDParam[1]=nmean;
835 fdNdEtaNBDParam[2]=beta;
838 //______________________________________________________________________________
839 void AliGlauberMC::SetdNdEtaTwoNBDParam( Double_t k1,
846 // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
847 fdNdEtaTwoNBDParam[0]=k1;
848 fdNdEtaTwoNBDParam[1]=nmean1;
849 fdNdEtaTwoNBDParam[2]=k2;
850 fdNdEtaTwoNBDParam[3]=nmean2;
851 fdNdEtaTwoNBDParam[4]=alpha;
852 fdNdEtaTwoNBDParam[5]=beta;
855 //______________________________________________________________________________
856 Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x) const
858 //Get particle density per unit of rapidity
859 //using two component model
861 return nnp*((1.-x)*fNpart/2.+x*fNcoll);
864 //______________________________________________________________________________
865 Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta,
869 //Get particle density per unit of rapidity
870 //using the GBW model
871 //Parameters: delta, lambda, snn
872 return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
873 * TMath::Power(fNpart,(1.-delta)/3./delta);
876 //_______________________________________________________________________________
877 Double_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta) const
879 //Get particle density per unit of rapidity
880 //using a aandomized number from a negative binomial distrubution
881 //Parameters: k = related to distribition width
882 // nmean = mean of distribution
883 // beta = set contribution of participants / binary collisions to multiplicity
885 for(Int_t i = 0; i<fNpart; i++)
887 mulnp+=NegativeBinomialRandom(k,nmean);
890 for(Int_t i = 0; i<fNcoll; i++)
892 mulnb+=NegativeBinomialRandom(k,nmean);
894 return (1-beta)*mulnp/2+beta*mulnb;
897 //______________________________________________________________________________
898 Double_t AliGlauberMC::GetdNdEtaTwoNBD ( Int_t k1,
903 Double_t beta ) const
905 //Get particle density per unit of rapidity
906 //using random numbers from two negative binomial distributions
907 //Parameters: k1 = related to distribition width of distribution 1
908 // nmean1 = mean of distribution 1
909 // k2 = related to distribition width of distribution 2
910 // nmean2 = mean of distribution 2
911 // alpha = set contributions of distrubitin 1 / distribution 2
912 // beta = set contribution of participants / binary collisions to multiplicity
914 for(Int_t i = 0; i<fNpart; i++)
916 mulnp+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
919 for(Int_t i = 0; i<fNcoll; i++)
921 mulnb+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
923 Double_t mul=(1-beta)*mulnp/2+beta*mulnb;
927 //______________________________________________________________________________
928 Double_t AliGlauberMC::GetEccentricityPart() const
930 //get participant eccentricity of participants
931 if (fNpart<2) return 0.0;
932 return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
935 //_____________________________________________________________________________
938 Double_t AliGlauberMC::GetEccentricityPartColl() const
940 //get participant eccentricity of binary collisions
941 if (fNcoll<2) return 0.0;
942 return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
945 //______________________________________________________________________________
946 Double_t AliGlauberMC::GetEccentricity() const
948 //get standard eccentricity of participants
949 if (fNpart<2) return 0.0;
950 return ((fSy2-fSx2)/(fSy2+fSx2));
953 //______________________________________________________________________________
954 Double_t AliGlauberMC::GetEccentricityColl() const
956 //get standard eccentricity of binary collisions
957 if (fNcoll<2) return 0.0;
958 return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
961 //______________________________________________________________________________
962 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
965 Int_t nAttempts = 10; // set indices, max attempts and max comparisons
966 Bool_t succes = kFALSE;
967 for(Int_t j=0; j<nAttempts; j++)
969 if(bgen<0||!succes) //get impactparameter
971 bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
973 if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
977 //______________________________________________________________________________
978 Double_t AliGlauberMC::GetEpsilon2Part() const
980 //get participant eccentricity of participants
981 if (fNpart<2) return 0.0;
982 return (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2);
983 //return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
985 //______________________________________________________________________________
986 Double_t AliGlauberMC::GetEpsilon3Part() const
988 //get participant eccentricity of participants
989 if (fNpart<2) return 0.0;
990 return (TMath::Sqrt(fMeanr2Cos3Phi*fMeanr2Cos3Phi+fMeanr2Sin3Phi*fMeanr2Sin3Phi)/fMeanr2);
991 //return (TMath::Sqrt(fMeanr3Cos3Phi*fMeanr3Cos3Phi+fMeanr3Sin3Phi*fMeanr3Sin3Phi)/fMeanr3);
993 //______________________________________________________________________________
994 Double_t AliGlauberMC::GetEpsilon4Part() const
996 //get participant eccentricity of participants
997 if (fNpart<2) return 0.0;
998 return (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2);
999 //return (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4);
1002 //______________________________________________________________________________
1003 Double_t AliGlauberMC::GetEpsilon5Part() const
1005 //get participant eccentricity of participants
1006 if (fNpart<2) return 0.0;
1007 return (TMath::Sqrt(fMeanr2Cos5Phi*fMeanr2Cos5Phi+fMeanr2Sin5Phi*fMeanr2Sin5Phi)/fMeanr2);
1008 //return (TMath::Sqrt(fMeanr5Cos5Phi*fMeanr5Cos5Phi+fMeanr5Sin5Phi*fMeanr5Sin5Phi)/fMeanr5);
1011 //______________________________________________________________________________
1012 Double_t AliGlauberMC::GetEpsilon2Coll() const
1014 //get epsilon2 of binary collisions
1015 if (fNcoll<2) return 0.0;
1016 return (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll);
1019 //______________________________________________________________________________
1020 Double_t AliGlauberMC::GetEpsilon3Coll() const
1022 //get epsilon3 of binary collisions
1023 if (fNcoll<2) return 0.0;
1024 return (TMath::Sqrt(fMeanr2Cos3PhiColl*fMeanr2Cos3PhiColl+fMeanr2Sin3PhiColl*fMeanr2Sin3PhiColl)/fMeanr2Coll);
1025 //return (TMath::Sqrt(fMeanr3Cos3PhiColl*fMeanr3Cos3PhiColl+fMeanr3Sin3PhiColl*fMeanr3Sin3PhiColl)/fMeanr3Coll);
1028 //______________________________________________________________________________
1029 Double_t AliGlauberMC::GetEpsilon4Coll() const
1031 //get epsilon4 of binary collisions
1032 if (fNcoll<2) return 0.0;
1033 return (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll);
1034 //return (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll);
1037 //______________________________________________________________________________
1038 Double_t AliGlauberMC::GetEpsilon5Coll() const
1040 //get epsilon5 of binary collisions
1041 if (fNcoll<2) return 0.0;
1042 return (TMath::Sqrt(fMeanr2Cos5PhiColl*fMeanr2Cos5PhiColl+fMeanr2Sin5PhiColl*fMeanr2Sin5PhiColl)/fMeanr2Coll);
1043 //return (TMath::Sqrt(fMeanr5Cos5PhiColl*fMeanr5Cos5PhiColl+fMeanr5Sin5PhiColl*fMeanr5Sin5PhiColl)/fMeanr5Coll);
1046 //______________________________________________________________________________
1047 void AliGlauberMC::Run(Int_t nevents)
1050 cout << "Generating " << nevents << " events..." << endl;
1051 TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
1052 TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
1055 fnt = new TNtuple(name,title,
1056 "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:VarEColl:VarEPart:VarEPartColl:dNdEta:dNdEtaGBW:dNdEtaNBD:dNdEtaTwoNBD:xsect:tAA:Epsl2:Epsl3:Epsl4:Epsl5:E2Coll:E3Coll:E4Coll:E5Coll");
1057 fnt->SetDirectory(0);
1061 for (Int_t i = 0; i<nevents; i++)
1084 v[11] = fMeanXSystem;
1085 v[12] = fMeanYSystem;
1090 v[17] = GetEccentricity();
1091 v[18] = GetEccentricityColl();
1092 v[19] = GetEccentricityPart();
1093 v[20] = GetEccentricityPartColl();
1096 v[21] = GetdNdEta( fdNdEtaParam[0],fdNdEtaParam[1] );
1097 v[22] = GetdNdEtaGBW( fdNdEtaGBWParam[0],fdNdEtaGBWParam[1],fdNdEtaGBWParam[2] );
1098 v[23] = GetdNdEtaNBD( TMath::Nint(fdNdEtaNBDParam[0]),
1100 fdNdEtaNBDParam[2] );
1101 v[24] = GetdNdEtaTwoNBD( TMath::Nint(fdNdEtaTwoNBDParam[0]),
1102 fdNdEtaTwoNBDParam[1],
1103 TMath::Nint(fdNdEtaTwoNBDParam[2]),
1104 fdNdEtaTwoNBDParam[3],
1105 fdNdEtaTwoNBDParam[4],
1106 fdNdEtaTwoNBDParam[5] );
1118 if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
1120 //_____________epsilon2,3,4,4_______
1121 v[27] = GetEpsilon2Part();
1122 v[28] = GetEpsilon3Part();
1123 v[29] = GetEpsilon4Part();
1124 v[30] = GetEpsilon5Part();
1125 v[31] = GetEpsilon2Coll();
1126 v[32] = GetEpsilon3Coll();
1127 v[33] = GetEpsilon4Coll();
1128 v[34] = GetEpsilon5Coll();
1132 if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
1134 std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events: " << q << " discarded events: " << u <<"."<< endl;
1137 //---------------------------------------------------------------------------------
1138 void AliGlauberMC::RunAndSaveNtuple( Int_t n,
1139 const Option_t *sysA,
1140 const Option_t *sysB,
1148 AliGlauberMC mcg(sysA,sysB,signn);
1149 mcg.SetMinDistance(mind);
1153 TNtuple *nt=mcg.GetNtuple();
1154 TFile out(fname,"recreate",fname,9);
1156 printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
1160 //---------------------------------------------------------------------------------
1161 void AliGlauberMC::RunAndSaveNucleons( Int_t n,
1162 const Option_t *sysA,
1164 const Double_t signn,
1170 AliGlauberMC mcg(sysA,sysB,signn);
1171 mcg.SetMinDistance(mind);
1174 out=new TFile(fname,"recreate",fname,9);
1176 for(Int_t ievent=0; ievent<n; ievent++)
1179 //get an event with at least one collision
1180 while(!mcg.NextEvent()) {}
1182 //access, save and (if wanted) print out nucleons
1183 TObjArray* nucleons=mcg.GetNucleons();
1184 if(!nucleons) continue;
1186 nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
1190 cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
1191 cout<<"B = "<<mcg.GetB()<<" Npart = "<<mcg.GetNpart()<<endl<<endl;
1192 printf("Nucleus\t X\t Y\t Z\tNcoll\n");
1193 Int_t nNucls=nucleons->GetEntries();
1194 for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
1196 AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
1198 if(nucl->IsInNucleusB()) nucleus='B';
1199 Double_t x=nucl->GetX();
1200 Double_t y=nucl->GetY();
1201 Double_t z=nucl->GetZ();
1202 Int_t ncoll=nucl->GetNColl();
1203 printf(" %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
1210 //---------------------------------------------------------------------------------
1211 void AliGlauberMC::Reset()
1214 delete fnt; fnt=NULL;