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