GetNumberOfCollisionsPerEvent added (suggested by G. Martinez)
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
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 /* $Id$ */
17 //
18 // Utility class to make simple Glauber type calculations 
19 //           for SYMMTRIC collision geometries (AA):
20 // Impact parameter, production points, reaction plane dependence
21 //
22 // The SimulateTrigger method can be used for simple MB and hard-process
23 // (binary scaling) trigger studies.
24 // 
25 // Some basic quantities can be visualized directly.
26 //
27 // The default set-up for PbPb or AUAu collisions can be read from a file 
28 // calling Init(1) or Init(2) if you want to read Almonds too.
29 //
30 // ***** If you change settings dont forget to call init afterwards, *****
31 // ***** in order to update the formulas with the new parameters.    *****
32 // 
33 // Author: andreas.morsch@cern.ch
34 //=================== Added by A. Dainese 11/02/04 ===========================
35 // Calculate path length for a parton with production point (x0,y0)
36 // and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
37 // in a collision with impact parameter b and functions that make use
38 // of it.
39 //=================== Added by A. Dainese 05/03/04 ===========================
40 // Calculation of line integrals I0 and I1
41 //  integral0 =    \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
42 //  integral1 =    \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
43 // mostly for use in the Quenching class
44 //=================== Added by C. Loizdes 27/03/04 ===========================
45 // Handling of AuAu collisions
46 // More get/set functions
47 // Comments, units and clearing of code
48 //
49
50 // from AliRoot
51 #include "AliFastGlauber.h"
52 // from root
53 #include <TStyle.h>
54 #include <TH1F.h>
55 #include <TH2F.h>
56 #include <TF1.h>
57 #include <TF2.h>
58 #include <TLegend.h>
59 #include <TCanvas.h>
60 #include <TRandom.h>
61 #include <TFile.h>
62 #include <TROOT.h>
63 #include <Riostream.h>
64
65 ClassImp(AliFastGlauber)
66
67 Float_t AliFastGlauber::fgBMax           = 0.;
68 TF1*    AliFastGlauber::fgWSb            = NULL;     
69 TF2*    AliFastGlauber::fgWSbz           = NULL;    
70 TF1*    AliFastGlauber::fgWSz            = NULL;     
71 TF1*    AliFastGlauber::fgWSta           = NULL;    
72 TF2*    AliFastGlauber::fgWStarfi        = NULL; 
73 TF2*    AliFastGlauber::fgWAlmond        = NULL; 
74 TF1*    AliFastGlauber::fgWStaa          = NULL;   
75 TF1*    AliFastGlauber::fgWSgeo          = NULL;   
76 TF1*    AliFastGlauber::fgWSbinary       = NULL;   
77 TF1*    AliFastGlauber::fgWSN            = NULL;   
78 TF1*    AliFastGlauber::fgWPathLength0   = NULL;   
79 TF1*    AliFastGlauber::fgWPathLength    = NULL;
80 TF1*    AliFastGlauber::fgWEnergyDensity = NULL;   
81 TF1*    AliFastGlauber::fgWIntRadius     = NULL;   
82 TF2*    AliFastGlauber::fgWKParticipants = NULL; 
83 TF1*    AliFastGlauber::fgWParticipants  = NULL; 
84 TF2*    AliFastGlauber::fgWAlmondCurrent = NULL;    
85 TF2     AliFastGlauber::fgWAlmondFixedB[40]; 
86 const Int_t AliFastGlauber::fgkMCInts = 100000;
87 Int_t AliFastGlauber::fgCounter = 0;       
88
89 AliFastGlauber::AliFastGlauber() : fName()
90 {
91   //  Default Constructor
92   //
93   fgCounter++;
94   if(fgCounter>1)
95     Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter);
96
97   //  Defaults for Pb
98   SetMaxImpact();
99   SetLengthDefinition();
100   SetPbPbLHC();
101 }
102
103 AliFastGlauber::~AliFastGlauber()
104 {
105   fgCounter--;
106   //if(fgCounter==0) Reset();
107 }
108
109 void AliFastGlauber::SetAuAuRhic()
110 {
111   //Set all parameters for RHIC
112   SetWoodSaxonParametersAu();
113   SetHardCrossSection();
114   SetNNCrossSection(42);
115   SetNucleus(197);
116   SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
117 }
118
119 void AliFastGlauber::SetPbPbLHC()
120 {
121   //Set all parameters for LHC
122   SetWoodSaxonParametersPb();
123   SetHardCrossSection();
124   SetNNCrossSection();
125   SetNucleus();
126   SetFileName();
127 }
128
129 void AliFastGlauber::Init(Int_t mode)
130 {
131   // Initialisation
132   //    mode = 0; all functions are calculated 
133   //    mode = 1; overlap function is read from file (for Pb-Pb only)
134   //    mode = 2; interaction almond functions are read from file 
135   //              USE THIS FOR PATH LENGTH CALC.!
136   //
137
138   // 
139   Reset(); 
140   //
141
142   //
143   //  Wood-Saxon
144   //
145   fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
146   fgWSb->SetParameter(0, fWSr0);
147   fgWSb->SetParameter(1, fWSd);
148   fgWSb->SetParameter(2, fWSw);
149   fgWSb->SetParameter(3, fWSn);
150
151   fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
152   fgWSbz->SetParameter(0, fWSr0);
153   fgWSbz->SetParameter(1, fWSd);
154   fgWSbz->SetParameter(2, fWSw);
155   fgWSbz->SetParameter(3, fWSn);
156
157   fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
158   fgWSz->SetParameter(0, fWSr0);
159   fgWSz->SetParameter(1, fWSd);
160   fgWSz->SetParameter(2, fWSw);
161   fgWSz->SetParameter(3, fWSn);
162
163   //
164   //  Thickness
165   //
166   fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
167     
168   //
169   //  Overlap Kernel
170   //
171   fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
172   fgWStarfi->SetParameter(0, 0.);     
173   fgWStarfi->SetNpx(200);     
174   fgWStarfi->SetNpy(20);    
175
176   //
177   //  Participants Kernel
178   //
179   fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
180   fgWKParticipants->SetParameter(0, 0.);     
181   fgWKParticipants->SetParameter(1, fSigmaNN);     
182   fgWKParticipants->SetParameter(2, fA);     
183   fgWKParticipants->SetNpx(200);     
184   fgWKParticipants->SetNpy(20);      
185
186   //
187   //  Overlap and Participants
188   //
189   if (!mode) {
190     fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
191     fgWStaa->SetNpx(100);
192     fgWStaa->SetParameter(0,fA);
193     fgWStaa->SetNpx(100);
194     fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
195     fgWParticipants->SetParameter(0, fSigmaNN);     
196     fgWParticipants->SetParameter(1, fA);     
197     fgWParticipants->SetNpx(100);
198   } else {
199     Info("Init","Reading overlap function from file %s",fName.Data()); 
200     TFile* f = new TFile(fName.Data());
201     if(!f->IsOpen()){
202       Fatal("Init", "Could not open file %s",fName.Data());
203     }
204     fgWStaa = (TF1*) f->Get("WStaa");
205     fgWParticipants = (TF1*) f->Get("WParticipants");
206     delete f;
207   }
208
209   //
210   //  Energy Density
211   //
212   fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
213   fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
214     
215   //
216   //  Geometrical Cross-Section
217   //
218   fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
219   fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
220   fgWSgeo->SetNpx(100);
221
222   //
223   //  Hard cross section (binary collisions)
224   //
225   fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
226   fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
227   fgWSbinary->SetNpx(100);
228
229   //
230   // Hard collisions per event
231   //
232   fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
233   fgWSN->SetNpx(100);
234
235   //
236   //  Almond shaped interaction region
237   //
238   fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
239   fgWAlmond->SetParameter(0, 0.);     
240   fgWAlmond->SetNpx(200);     
241   fgWAlmond->SetNpy(200);  
242   
243   if(mode==2) {
244     Info("Init","Reading interaction almonds from file: %s",fName.Data());
245     Char_t almondName[100];
246     TFile* ff = new TFile(fName.Data());
247     for(Int_t k=0; k<40; k++) {
248       sprintf(almondName,"WAlmondFixedB%d",k);
249       fgWAlmondCurrent = (TF2*)ff->Get(almondName);
250       new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
251     }
252     delete ff;
253   }
254
255   fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
256   fgWIntRadius->SetParameter(0, 0.);    
257
258   //
259   //  Path Length as a function of Phi
260   //    
261   fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
262   fgWPathLength0->SetParameter(0, 0.);
263   fgWPathLength0->SetParameter(1, 0.); //Pathlength definition     
264
265   fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
266   fgWPathLength->SetParameter(0, 0.);    //Impact Parameter
267   fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
268   fgWPathLength->SetParameter(2, 0);     //Pathlength definition
269 }
270
271 void AliFastGlauber::Reset()
272 {
273   //
274   // Reset dynamic allocated formulas
275   // in case init is called twice
276
277   if(fgWSb)            delete fgWSb;     
278   if(fgWSbz)           delete fgWSbz;
279   if(fgWSz)            delete fgWSz;
280   if(fgWSta)           delete fgWSta;
281   if(fgWStarfi)        delete fgWStarfi;
282   if(fgWAlmond)        delete fgWAlmond;
283   if(fgWStaa)          delete fgWStaa;
284   if(fgWSgeo)          delete fgWSgeo;
285   if(fgWSbinary)       delete fgWSbinary;
286   if(fgWSN)            delete fgWSN;
287   if(fgWPathLength0)   delete fgWPathLength0;
288   if(fgWPathLength)    delete fgWPathLength;
289   if(fgWEnergyDensity) delete fgWEnergyDensity;
290   if(fgWIntRadius)     delete fgWIntRadius; 
291   if(fgWKParticipants) delete fgWKParticipants; 
292   if(fgWParticipants)  delete fgWParticipants;
293 }
294
295 void AliFastGlauber::DrawWSb() const
296 {
297   //
298   //  Draw Wood-Saxon Nuclear Density Function
299   //
300   TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
301   c1->cd();
302   Double_t max=fgWSb->GetMaximum(0,fgBMax)*1.01;
303   TH2F *h2f=new TH2F("h2fwsb","Wood Saxon: #rho(r) = n (1-#omega(r/r_{0})^2)/(1+exp((r-r_{0})/d)) [fm^{-3}]",2,0,fgBMax,2,0,max);
304   h2f->SetStats(0);
305   h2f->GetXaxis()->SetTitle("r [fm]");
306   h2f->GetYaxis()->SetNoExponent(kTRUE);
307   h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
308   h2f->Draw(); 
309   fgWSb->Draw("same");
310   TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
311   l1a->SetFillStyle(0);
312   l1a->SetBorderSize(0);
313   Char_t label[100];
314   sprintf(label,"r_{0} = %.2f fm",fWSr0);
315   l1a->AddEntry(fgWSb,label,"");
316   sprintf(label,"d = %.2f fm",fWSd);
317   l1a->AddEntry(fgWSb,label,"");
318   sprintf(label,"n = %.2e fm^{-3}",fWSn);
319   l1a->AddEntry(fgWSb,label,"");
320   sprintf(label,"#omega = %.2f",fWSw);
321   l1a->AddEntry(fgWSb,label,"");
322   l1a->Draw();
323   c1->Update();
324 }
325
326 void AliFastGlauber::DrawOverlap() const
327 {
328   //
329   //  Draw Overlap Function
330   //
331   TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
332   c2->cd();
333   Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
334   TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0,max);
335   h2f->SetStats(0);
336   h2f->GetXaxis()->SetTitle("b [fm]");
337   h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
338   h2f->Draw(); 
339   fgWStaa->Draw("same");
340 }
341
342 void AliFastGlauber::DrawParticipants() const
343 {
344   //
345   //  Draw Number of Participants Npart
346   //
347   TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
348   c3->cd();
349   Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
350   TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
351   h2f->SetStats(0);
352   h2f->GetXaxis()->SetTitle("b [fm]");
353   h2f->GetYaxis()->SetTitle("N_{part}");
354   h2f->Draw(); 
355   fgWParticipants->Draw("same");
356   TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
357   l1a->SetFillStyle(0);
358   l1a->SetBorderSize(0);
359   Char_t label[100];
360   sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
361   l1a->AddEntry(fgWParticipants,label,"");
362   l1a->Draw();
363   c3->Update();
364 }
365
366 void AliFastGlauber::DrawThickness() const
367 {
368   //
369   //  Draw Thickness Function
370   //
371   TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
372   c4->cd();
373   Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
374   TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
375   h2f->SetStats(0);
376   h2f->GetXaxis()->SetTitle("b [fm]");
377   h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
378   h2f->Draw(); 
379   fgWSta->Draw("same");
380 }
381
382 void AliFastGlauber::DrawGeo() const
383 {
384   //
385   //  Draw Geometrical Cross-Section
386   //
387   TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
388   c5->cd();
389   Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
390   TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
391   h2f->SetStats(0);
392   h2f->GetXaxis()->SetTitle("b [fm]");
393   h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
394   h2f->Draw(); 
395   fgWSgeo->Draw("same");
396   TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
397   l1a->SetFillStyle(0);
398   l1a->SetBorderSize(0);
399   Char_t label[100];
400   sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
401   l1a->AddEntry(fgWSgeo,label,"");
402   l1a->Draw();
403   c5->Update();
404 }
405
406 void AliFastGlauber::DrawBinary() const
407 {
408   //
409   //  Draw Binary Cross-Section
410   //
411   TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
412   c6->cd();
413   Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
414   TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
415   h2f->SetStats(0);
416   h2f->GetXaxis()->SetTitle("b [fm]");
417   h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
418   h2f->Draw(); 
419   fgWSbinary->Draw("same");
420   TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
421   l1a->SetFillStyle(0);
422   l1a->SetBorderSize(0);
423   Char_t label[100];
424   sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
425   l1a->AddEntry(fgWSb,label,"");
426   l1a->Draw();
427   c6->Update();
428 }
429
430 void AliFastGlauber::DrawN() const
431 {
432   //
433   //  Draw Binaries per event (Ncoll)
434   //
435   TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
436   c7->cd();
437   Double_t max=fgWSN->GetMaximum(0,fgBMax)*1.01;
438   TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
439   h2f->SetStats(0);
440   h2f->GetXaxis()->SetTitle("b [fm]");
441   h2f->GetYaxis()->SetTitle("N_{coll}");
442   h2f->Draw(); 
443   fgWSN->Draw("same");
444   TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
445   l1a->SetFillStyle(0);
446   l1a->SetBorderSize(0);
447   Char_t label[100];
448   sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
449   l1a->AddEntry(fgWSN,label,"");
450   sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
451   l1a->AddEntry(fgWSN,label,"");
452   l1a->Draw();
453   c7->Update();
454 }
455
456 void AliFastGlauber::DrawKernel(Double_t b) const
457 {
458   //
459   //  Draw Kernel
460   //
461   TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
462   c8->cd();
463   fgWStarfi->SetParameter(0, b);
464   TH2F *h2f=new TH2F("h2fkernel","Kernel of Overlap function: d^{2}T_{AB}/dr/d#phi [fm^{-3}]",2,0,fgBMax,2,0,TMath::Pi());
465   h2f->SetStats(0);
466   h2f->GetXaxis()->SetTitle("r [fm]");
467   h2f->GetYaxis()->SetTitle("#phi [rad]");
468   h2f->Draw(); 
469   fgWStarfi->Draw("same");
470   TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
471   l1a->SetFillStyle(0);
472   l1a->SetBorderSize(0);
473   Char_t label[100];
474   sprintf(label,"b = %.1f fm",b);
475   l1a->AddEntry(fgWStarfi,label,"");
476   l1a->Draw();
477   c8->Update();
478 }
479
480 void AliFastGlauber::DrawAlmond(Double_t b) const
481 {
482   //
483   //  Draw Interaction Almond
484   //
485   TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
486   c9->cd();
487   fgWAlmond->SetParameter(0, b);
488   TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,0,fgBMax,2,0,fgBMax);
489   h2f->SetStats(0);
490   h2f->GetXaxis()->SetTitle("x [fm]");
491   h2f->GetYaxis()->SetTitle("y [fm]");
492   h2f->Draw(); 
493   fgWAlmond->Draw("same");
494   TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
495   l1a->SetFillStyle(0);
496   l1a->SetBorderSize(0);
497   Char_t label[100];
498   sprintf(label,"b = %.1f fm",b);
499   l1a->AddEntry(fgWAlmond,label,"");
500   l1a->Draw();
501   c9->Update();
502 }
503
504 void AliFastGlauber::DrawEnergyDensity() const
505 {
506   //
507   //  Draw energy density
508   //
509   TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
510   c10->cd();
511   fgWEnergyDensity->SetMinimum(0.);
512   Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
513   TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
514   h2f->SetStats(0);
515   h2f->GetXaxis()->SetTitle("b [fm]");
516   h2f->GetYaxis()->SetTitle("fm^{-4}");
517   h2f->Draw(); 
518   fgWEnergyDensity->Draw("same");
519   c10->Update();
520 }
521
522 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
523 {
524   //
525   //  Draw Path Length
526   //
527   TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
528   c11->cd();
529   fgWPathLength0->SetParameter(0, b);
530   fgWPathLength0->SetParameter(1, Double_t(iopt));
531   fgWPathLength0->SetMinimum(0.); 
532   fgWPathLength0->SetMaximum(10.); 
533   TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
534   h2f->SetStats(0);
535   h2f->GetXaxis()->SetTitle("#phi [rad]");
536   h2f->GetYaxis()->SetTitle("l [fm]");
537   h2f->Draw(); 
538   fgWPathLength0->Draw("same");
539 }
540
541 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
542 {
543   //
544   //  Draw Path Length
545   //
546   TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
547   c12->cd();
548   fgWAlmond->SetParameter(0, b);
549   fgWPathLength->SetParameter(0, b);
550   fgWPathLength->SetParameter(1, Double_t (ni));
551   fgWPathLength->SetParameter(2, Double_t (iopt));
552   fgWPathLength->SetMinimum(0.); 
553   fgWPathLength->SetMaximum(10.); 
554   TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
555   h2f->SetStats(0);
556   h2f->GetXaxis()->SetTitle("#phi [rad]");
557   h2f->GetYaxis()->SetTitle("l [fm]");
558   h2f->Draw(); 
559   fgWPathLength->Draw("same");
560 }
561
562 void AliFastGlauber::DrawIntRadius(Double_t b) const
563 {
564   //
565   //  Draw Interaction Radius
566   //
567   TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
568   c13->cd();
569   fgWIntRadius->SetParameter(0, b);
570   fgWIntRadius->SetMinimum(0);
571   Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
572   TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
573   h2f->SetStats(0);
574   h2f->GetXaxis()->SetTitle("r [fm]");
575   h2f->GetYaxis()->SetTitle("[fm^{-3}]");
576   h2f->Draw(); 
577   fgWIntRadius->Draw("same");
578 }
579
580 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
581 {
582   //
583   //  Woods-Saxon Parameterisation
584   //  as a function of radius (xx)
585   //
586   const Double_t kxx  = x[0];   //fm
587   const Double_t kr0  = par[0]; //fm
588   const Double_t kd   = par[1]; //fm   
589   const Double_t kw   = par[2]; //no units
590   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
591   Double_t y   = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
592   return y; //fm^-3
593 }
594
595 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
596 {
597   //
598   //  Wood Saxon Parameterisation
599   //  as a function of z and  b
600   //
601   const Double_t kbb  = x[0];   //fm
602   const Double_t kzz  = x[1];   //fm
603   const Double_t kr0  = par[0]; //fm
604   const Double_t kd   = par[1]; //fm
605   const Double_t kw   = par[2]; //no units
606   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
607   const Double_t kxx  = TMath::Sqrt(kbb*kbb+kzz*kzz);
608   Double_t y  = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
609   return y; //fm^-3
610 }
611
612 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
613 {
614   //
615   //  Wood Saxon Parameterisation
616   //  as a function of z for fixed b
617   //
618   const Double_t kzz  = x[0];   //fm
619   const Double_t kr0  = par[0]; //fm
620   const Double_t kd   = par[1]; //fm
621   const Double_t kw   = par[2]; //no units
622   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
623   const Double_t kbb  = par[4]; //fm
624   const Double_t kxx  = TMath::Sqrt(kbb*kbb+kzz*kzz);
625   Double_t y  = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
626   return y; //fm^-3
627 }
628
629 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
630 {
631   //
632   //  Thickness function T_A
633   //  as a function of b
634   //
635   const Double_t kb  = x[0];
636   fgWSz->SetParameter(4, kb);
637   Double_t y  = 2. * fgWSz->Integral(0., fgBMax);
638   return y; //fm^-2
639 }
640
641 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
642 {
643   //
644   //  Kernel for overlap function: T_A(s)*T_A(s-b)
645   //  as a function of r and phi
646   const Double_t kr1  = x[0];
647   const Double_t kphi = x[1];
648   const Double_t kb   = par[0];
649   const Double_t kr2  = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi)); 
650   Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
651   return y; //fm^-3
652 }
653
654 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
655 {
656   //
657   //  Overlap function 
658   //  T_{AB}=Int d2s T_A(s)*T_B(s-b)
659   //  as a function of b
660   // (normalized to fA*fB)
661   //
662   const Double_t kb  = x[0];
663   const Double_t ka = par[0];
664   fgWStarfi->SetParameter(0, kb);
665
666   // root integration seems to fail
667   /* 
668      Double_t al[2];
669      Double_t bl[2];
670      al[0] = 1e-6;
671      al[1] = fgBMax;
672      bl[0] = 0.;
673      bl[1] = TMath::Pi();
674      Double_t err;
675      
676      Double_t y =  2. *  208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
677      printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
678   */
679
680   //
681   //  MC Integration
682   //
683   Double_t y = 0;
684   for (Int_t i = 0; i < fgkMCInts; i++)
685     {
686       const Double_t kphi = TMath::Pi() * gRandom->Rndm();
687       const Double_t kb1  = fgBMax      * gRandom->Rndm();      
688       y += fgWStarfi->Eval(kb1, kphi);
689     }
690   y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
691   y *= ka * ka * 0.1; //mbarn^-1
692   return y;
693 }
694
695 Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
696 {
697   //
698   //  Kernel for number of participants
699   //  as a function of r and phi
700   //
701   const Double_t kr1   = x[0];
702   const Double_t kphi  = x[1];
703   const Double_t kb    = par[0]; //fm
704   const Double_t ksig  = par[1]; //mbarn
705   const Double_t ka    = par[2]; //mass number
706   const Double_t kr2   = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi)); 
707   const Double_t kxsi  = fgWSta->Eval(kr2) * ksig * 0.1; //no units
708   /*
709     Double_t y=(1-TMath::Power((1-xsi),aa))
710    */
711   Double_t a   = ka;
712   Double_t sum = ka * kxsi;
713   Double_t y   = sum;
714   for (Int_t i = 1; i <= ka; i++)
715     {
716       a--;
717       sum *= (-kxsi) * a / Float_t(i+1);
718       y  += sum;
719     }
720   y *= kr1 * fgWSta->Eval(kr1);
721   return y; //fm^-1
722 }
723
724 Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
725 {
726   //
727   //  Number of Participants as 
728   //  a function of b
729   //
730   const Double_t kb = x[0];
731   const Double_t ksig  = par[0]; //mbarn
732   const Double_t ka   = par[1];  //mass number
733   fgWKParticipants->SetParameter(0, kb);
734   fgWKParticipants->SetParameter(1, ksig);
735   fgWKParticipants->SetParameter(2, ka);
736
737   //
738   //  MC Integration
739   //
740   Double_t y = 0;
741   for (Int_t i = 0; i < fgkMCInts; i++)
742     {
743       const Double_t kphi = TMath::Pi() * gRandom->Rndm();
744       const Double_t kb1  = fgBMax      * gRandom->Rndm();      
745       y += fgWKParticipants->Eval(kb1, kphi);
746     }
747   y *= 2. *  ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
748   return y; //no units
749 }
750
751 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
752 {
753   //
754   //  Geometrical Cross-Section
755   //  as a function of b
756   //
757   const Double_t kb     = x[0];              //fm
758   const Double_t ksigNN = par[0];            //mbarn
759   const Double_t ktaa   = fgWStaa->Eval(kb); //mbarn^-1
760   Double_t y     = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa)); 
761   return y; //fm
762 }
763
764 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
765 {
766   //
767   //  Number of binary hard collisions
768   //  as a function of b
769   //
770   const Double_t kb     = x[0];              //fm
771   const Double_t ksig   = par[0];            //mbarn
772   const Double_t ktaa   = fgWStaa->Eval(kb); //mbarn^-1
773   Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa; 
774   return y; //fm
775 }
776
777 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
778 {
779   //
780   //  Number of hard processes per event
781   //  as a function of b
782   const Double_t kb = x[0];
783   Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
784   return y; //no units
785 }
786
787 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
788 {
789   //
790   //  Initial energy density 
791   //  as a function of the impact parameter
792   //
793   const Double_t kb     = x[0];
794   const Double_t krA    = par[0];
795   //
796   //  Attention: area of transverse reaction zone in hard-sphere approximation !     
797   const Double_t krA2=krA*krA;
798   const Double_t kb2=kb*kb;  
799   const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2 
800                       - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
801   const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
802   Double_t y=ktaa/ksaa*10;
803   return y; //fm^-4
804 }
805
806 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
807 {
808   //
809   //  Almond shaped interaction region
810   //  as a function of cartesian x,y.
811   //
812   const Double_t kb    = par[0];
813   const Double_t kxx   = x[0] + kb/2.;
814   const Double_t kyy   = x[1];
815   const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
816   const Double_t kphi  = TMath::ATan2(kyy,kxx);
817   const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi)); 
818   //
819   //  Interaction probability calculated as product of thicknesses
820   //
821   Double_t y    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
822   return y; //fm^-4
823 }
824
825 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
826 {
827   //
828   //  Average interaction density over radius 
829   //  at which interaction takes place
830   //  as a function of radius
831   //
832   const Double_t kr    = x[0];
833   const Double_t kb    = par[0];
834   fgWAlmond->SetParameter(0, kb);
835   //  Average over phi in small steps   
836   const Double_t kdphi = 2. * TMath::Pi() / 100.;
837   Double_t phi  = 0.;
838   Double_t y    = 0.;
839   for (Int_t i = 0; i < 100; i++) {
840     const Double_t kxx = kr * TMath::Cos(phi);
841     const Double_t kyy = kr * TMath::Sin(phi);
842     y   += fgWAlmond->Eval(kxx,kyy);
843     phi += kdphi;
844   } // phi loop
845   // Result multiplied by Jacobian (2 pi r)     
846   y *= 2. * TMath::Pi() * kr / 100.;
847   return y; //fm^-3
848 }
849
850 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
851 {
852   //
853   //  Path Length as a function of phi 
854   //  for interaction point fixed at (0,0)
855   //  as a function of phi-direction
856   //
857   //  Phi direction in Almond
858   const Double_t kphi0   = x[0];
859   const Double_t kb      = par[0];
860   //  Path Length definition
861   const Int_t    kiopt   = Int_t(par[1]);
862
863   //  Step along radial direction phi   
864   const Int_t    kNp  = 100; // Steps in r 
865   const Double_t kDr  = fgBMax/kNp;
866   Double_t r  = 0.;
867   Double_t rw = 0.;
868   Double_t w  = 0.;
869   for (Int_t i = 0; i < kNp; i++) {
870     //
871     //  Transform into target frame
872     //
873     const Double_t kxx   = r * TMath::Cos(kphi0) + kb / 2.;
874     const Double_t kyy   = r * TMath::Sin(kphi0);
875     const Double_t kphi  = TMath::ATan2(kyy, kxx);
876     const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
877     // Radius in projectile frame
878     const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi)); 
879     const Double_t ky    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
880
881     rw += ky * r;
882     w  += ky;
883     r  += kDr;
884   } // radial steps
885
886   Double_t y=0.;
887   if (!kiopt)  // My length definition (is exact for hard disk)
888     if(w) y= 2. * rw / w; 
889   else {
890     const Double_t knorm=fgWSta->Eval(1e-4);
891     if(knorm) y =  TMath::Sqrt(2. * rw * kDr / knorm / knorm);
892   }
893   return y; //fm
894 }
895
896 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
897 {
898   //
899   //  Path Length as a function of phi 
900   //  Interaction point from random distribution
901   //  as a function of the phi-direction
902   const Double_t kphi0   = x[0];
903   const Double_t kb      = par[0];
904   fgWAlmond->SetParameter(0, kb); 
905   const Int_t    kNpi    = Int_t (par[1]); //Number of interactions
906   const Int_t    kiopt   = Int_t(par[2]);  //Path Length definition 
907
908   //
909   //  r-steps
910   // 
911   const Int_t    kNp   = 100;
912   const Double_t kDr  = fgBMax/Double_t(kNp);
913   Double_t l = 0.;  //  Path length 
914   for (Int_t in = 0; in < kNpi; in ++) {
915     Double_t rw = 0.;
916     Double_t w  = 0.;
917     // Interaction point
918     Double_t x0, y0;
919     fgWAlmond->GetRandom2(x0, y0);
920     // Initial radius
921     const Double_t kr0  = TMath::Sqrt(x0*x0 + y0*y0);
922     const Int_t    knps = Int_t ((fgBMax - kr0)/kDr) - 1;
923         
924     // Radial steps
925     Double_t r  = 0.;
926     for (Int_t i = 0; (i < knps ); i++) {
927       // Transform into target frame
928       const Double_t kxx   = x0 + r * TMath::Cos(kphi0) + kb / 2.;
929       const Double_t kyy   = y0 + r * TMath::Sin(kphi0);
930       const Double_t kphi  = TMath::ATan2(kyy, kxx);
931       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
932       // Radius in projectile frame
933       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi)); 
934       const Double_t ky    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
935             
936       rw += ky * r;
937       w  += ky;
938       r  += kDr;
939     } // steps
940     // Average over interactions
941     if (!kiopt) {
942       if(w) l += (2. * rw / w);
943     } else {
944       const Double_t knorm=fgWSta->Eval(1e-4);
945       if(knorm) l+= 2. * rw * kDr / knorm / knorm;
946     }
947   } // interactions
948   Double_t ret=0;
949   if (!kiopt) 
950     ret= l / kNpi;
951   else 
952     ret=TMath::Sqrt( l / kNpi);
953   return ret; //fm
954 }
955
956 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
957 {
958   //
959   // Return the geometrical cross-section integrated from b1 to b2 
960   //
961   return fgWSgeo->Integral(b1, b2)*10.; //mbarn
962 }
963
964 Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
965 {
966   //
967   // Return the hard cross-section integrated from b1 to b2 
968   //
969   return fgWSbinary->Integral(b1, b2)*10.; //mbarn
970 }
971
972 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
973 {
974   //
975   // Return fraction of hard cross-section integrated from b1 to b2 
976   //
977   return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
978 }
979
980 Double_t AliFastGlauber::Binaries(Double_t b) const
981 {
982   //
983   // Return number of binary hard collisions normalized to 1 at b=0
984   //
985   if(b==0) b=1e-4;
986   return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
987 }
988
989 Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
990 {
991   //
992   // Return number of binary hard collisions at b
993   //
994   if(b==0) b=1e-4;
995   return fgWSN->Eval(b);
996 }
997
998 Double_t AliFastGlauber::Participants(Double_t  b) const
999 {
1000   //
1001   // Return the number of participants normalized to 1 at b=0
1002   //
1003   if(b==0) b=1e-4;
1004   return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
1005 }
1006
1007 Double_t AliFastGlauber::GetNumberOfParticipants(Double_t  b) const
1008 {
1009   //
1010   // Return the number of participants for impact parameter b
1011   //
1012   if(b==0) b=1e-4;
1013   return (fgWParticipants->Eval(b));
1014 }
1015
1016 Double_t AliFastGlauber::GetNumberOfCollisions(Double_t  b) const
1017 {
1018   //
1019   // Return the number of collisions for impact parameter b
1020   //
1021   if(b==0) b=1e-4;
1022   return (fgWStaa->Eval(b)*fSigmaNN);
1023 }
1024
1025 Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t  b) const
1026 {
1027   //
1028   // Return the number of collisions per event (at least one collision)
1029   // for impact parameter b
1030   //
1031     Double_t n = GetNumberOfCollisions(b);
1032     if (n > 0.) {
1033         return (n / (1. - TMath::Exp(- n)));
1034     } else {
1035         return (0.);
1036     }
1037 }
1038
1039 void AliFastGlauber::SimulateTrigger(Int_t n)
1040 {
1041   //
1042   //  Simulates Trigger
1043   //
1044   TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
1045   TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);   
1046   TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
1047   TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);   
1048
1049   mbtH->SetXTitle("b [fm]");
1050   hdtH->SetXTitle("b [fm]");    
1051   mbmH->SetXTitle("Multiplicity");
1052   hdmH->SetXTitle("Multiplicity");    
1053
1054   TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);    
1055   c0->Divide(2,1);
1056   TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);    
1057   c1->Divide(1,2);
1058
1059   //
1060   //
1061   Init(1);
1062   for (Int_t iev = 0; iev < n; iev++)
1063     {
1064       Float_t b, p, mult;
1065       GetRandom(b, p, mult);
1066       mbtH->Fill(b,1.);
1067       hdtH->Fill(b, p);
1068       mbmH->Fill(mult, 1.);
1069       hdmH->Fill(mult, p);
1070
1071       c0->cd(1);
1072       mbtH->Draw();
1073       c0->cd(2);
1074       hdtH->Draw();     
1075       c0->Update();
1076
1077       c1->cd(1);
1078       mbmH->Draw();
1079       c1->cd(2);
1080       hdmH->Draw();     
1081       c1->Update();
1082     }
1083 }
1084
1085 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1086 {
1087   //
1088   // Gives back a random impact parameter, hard trigger probability and multiplicity
1089   //
1090   b = fgWSgeo->GetRandom();
1091   const Float_t kmu = fgWSN->Eval(b);
1092   p = 1.-TMath::Exp(-kmu);
1093   mult = 6000./fgWSN->Eval(1.) * kmu;
1094 }
1095
1096 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1097 {
1098   //
1099   // Gives back a random impact parameter bin, and hard trigger decission
1100   //
1101   const Float_t kb  = fgWSgeo->GetRandom();
1102   const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1103   const Float_t kp  = 1.-TMath::Exp(-kmu);
1104   if (kb < 5.) {
1105     bin = 1;
1106   } else if (kb <  8.6) {
1107     bin = 2;
1108   } else if (kb < 11.2) {
1109     bin = 3;
1110   } else if (kb < 13.2) {
1111     bin = 4;
1112   } else if (kb < 15.0) {
1113     bin = 5;
1114   } else {
1115     bin = 6;
1116   }
1117   hard = kFALSE;
1118   const Float_t kr = gRandom->Rndm();
1119   if (kr < kp) hard = kTRUE;
1120 }
1121
1122 Double_t  AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
1123 {
1124   //
1125   // Gives back a random impact parameter in the range bmin .. bmax
1126   //
1127   Float_t b = -1.;
1128   while(b < bmin || b > bmax)
1129     b = fgWSgeo->GetRandom();
1130   return b;
1131 }
1132
1133 void AliFastGlauber::StoreFunctions() const
1134 {
1135   //
1136   // Store in file functions
1137   //
1138   TFile* ff = new TFile(fName.Data(),"recreate");
1139   fgWStaa->Write("WStaa");
1140   fgWParticipants->Write("WParticipants");
1141   ff->Close();
1142   return;
1143 }
1144
1145 //=================== Added by A. Dainese 11/02/04 ===========================
1146
1147 void AliFastGlauber::StoreAlmonds() const
1148 {
1149   //
1150   // Store in file 
1151   // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1152   //
1153   Char_t almondName[100];
1154   TFile* ff = new TFile(fName.Data(),"update");
1155   for(Int_t k=0; k<40; k++) {
1156     sprintf(almondName,"WAlmondFixedB%d",k);
1157     Double_t b = 0.25+k*0.5;
1158     Info("StoreAlmonds"," b = %f\n",b); 
1159     fgWAlmond->SetParameter(0,b);
1160     fgWAlmond->Write(almondName);
1161   }
1162   ff->Close();
1163   return;
1164 }
1165
1166 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1167 {
1168   //
1169   // Set limits of centrality class as fractions 
1170   // of the geomtrical cross section  
1171   //
1172   if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
1173     Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
1174     return;
1175   }
1176
1177   Double_t bLow=0.,bUp=0.;
1178   Double_t xsecFr=0.;
1179   const Double_t knorm=fgWSgeo->Integral(0.,100.);
1180   while(xsecFr<xsecFrLow) {
1181     xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
1182     bLow += 0.1;
1183   }
1184   bUp = bLow;
1185   while(xsecFr<xsecFrUp) {
1186     xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
1187     bUp += 0.1;
1188   }
1189
1190   Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
1191          xsecFrLow,xsecFrUp,bLow,bUp);
1192   fgWSbinary->SetRange(bLow,bUp);
1193   fBmin=bLow;
1194   fBmax=bUp;
1195   return;
1196 }
1197
1198 void AliFastGlauber::GetRandomBHard(Double_t& b)
1199 {
1200   //
1201   // Get random impact parameter according to distribution of 
1202   // hard (binary) cross-section, in the range defined by the centrality class
1203   //
1204   b = fgWSbinary->GetRandom();
1205   Int_t bin = 2*(Int_t)b;
1206   if( (b-(Int_t)b) > 0.5) bin++;
1207   fgWAlmondCurrent = &fgWAlmondFixedB[bin]; 
1208   return;
1209 }
1210
1211 void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1212 {
1213   //
1214   // Get random position of parton production point according to 
1215   // product of thickness functions
1216   //
1217   fgWAlmondCurrent->GetRandom2(x,y);
1218   return;
1219 }
1220
1221 void AliFastGlauber::GetRandomPhi(Double_t& phi) 
1222 {
1223   //
1224   // Get random parton azimuthal propagation direction
1225   //
1226   phi = 2.*TMath::Pi()*gRandom->Rndm();
1227   return;
1228 }
1229
1230 Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
1231 {
1232   // 
1233   // Calculate path length for a parton with production point (x0,y0)
1234   // and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
1235   // in a collision with impact parameter b
1236   //
1237
1238   // number of steps in l
1239   const Int_t    kNp  = 100;
1240   const Double_t kDl  = fgBMax/Double_t(kNp);
1241
1242   if(fEllDef==1) {
1243     //
1244     // Definition 1:
1245     // 
1246     // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
1247     //           \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1248     //
1249
1250     // Initial radius
1251     const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1252     const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1253     Double_t l  = 0.;
1254     Double_t integral1 = 0.;
1255     Double_t integral2 = 0.;
1256     // Radial steps
1257     for (Int_t i = 0; i < knps; i++) {
1258       
1259       // Transform into target frame
1260       const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1261       const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1262       const Double_t kphi  = TMath::ATan2(kyy, kxx);
1263       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1264       // Radius in projectile frame
1265       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi)); 
1266       const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1267       
1268       integral1 += kprodTATB * l * kDl;
1269       integral2 += kprodTATB * kDl;
1270       l  += kDl;
1271     } // steps
1272     
1273     Double_t ell=0.;
1274     if(integral2)
1275       ell = (2. * integral1 / integral2);
1276     return ell;
1277   } else if(fEllDef==2) {
1278     //
1279     // Definition 2:
1280     // 
1281     // ell = \int_0^\infty dl*
1282     //       \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
1283     //
1284
1285     // Initial radius
1286     const Double_t kr0  = TMath::Sqrt(x0*x0 + y0*y0);
1287     const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1288     const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
1289     // Radial steps
1290     Double_t l  = 0.;
1291     Double_t integral = 0.;
1292     for (Int_t i = 0; i < knps; i++) {
1293       // Transform into target frame
1294       const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1295       const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1296       const Double_t kphi  = TMath::ATan2(kyy, kxx);
1297       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1298       // Radius in projectile frame
1299       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi)); 
1300       const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1301       if(kprodTATB>kprodTATBHalfMax) integral += kDl;
1302       l  += kDl;
1303     } // steps
1304     Double_t ell = integral;
1305     return ell;
1306   } else {
1307     Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
1308     return -1.;
1309   }
1310 }
1311
1312 void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
1313 {
1314   //
1315   // Return length from random b, x0, y0, phi0 
1316   // Return also phi0
1317   //
1318   Double_t x0,y0,phi0;
1319   if(b<0.) GetRandomBHard(b);
1320   GetRandomXY(x0,y0);
1321   GetRandomPhi(phi0);
1322   phi = phi0;
1323   ell = CalculateLength(b,x0,y0,phi0);
1324   return;
1325 }
1326
1327 void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1328 {
1329   //
1330   // Return length from random b, x0, y0, phi0 
1331   //
1332   Double_t phi;
1333   GetLengthAndPhi(ell,phi,b);
1334   return;
1335 }
1336
1337 void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1338 {
1339   //
1340   // Return 2 lengths back to back from random b, x0, y0, phi0 
1341   // Return also phi0 
1342  // 
1343   Double_t x0,y0,phi0;
1344   if(b<0.) GetRandomBHard(b);
1345   GetRandomXY(x0,y0);
1346   GetRandomPhi(phi0);
1347   const Double_t kphi0plusPi = phi0+TMath::Pi();
1348   phi = phi0;
1349   ell1 = CalculateLength(b,x0,y0,phi0);
1350   ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
1351   return;
1352 }
1353
1354 void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1355                                           Double_t b)
1356 {
1357   //
1358   // Return 2 lengths back to back from random b, x0, y0, phi0 
1359   // 
1360   Double_t phi;
1361   GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1362   return;
1363 }
1364
1365 void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
1366 {
1367   //
1368   // Returns lenghts for n partons with azimuthal angles phi[n] 
1369   // from random b, x0, y0
1370   //
1371   Double_t x0, y0;
1372   if(b < 0.) GetRandomBHard(b);
1373   GetRandomXY(x0,y0);
1374   for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
1375   return;
1376 }
1377
1378 void AliFastGlauber::PlotBDistr(Int_t n)
1379 {  
1380   // 
1381   // Plot distribution of n impact parameters
1382   //
1383   Double_t b;
1384   TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax); 
1385   hB->SetXTitle("b [fm]");
1386   hB->SetYTitle("dN/db [a.u.]");
1387   hB->SetFillColor(3);
1388   for(Int_t i=0; i<n; i++) {
1389     GetRandomBHard(b);
1390     hB->Fill(b);
1391   }
1392   TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1393   cB->cd();
1394   hB->Draw();
1395   return;
1396 }
1397
1398 void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
1399 {
1400   //
1401   // Plot length distribution
1402   //
1403   Double_t ell;
1404   TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15); 
1405   hEll->SetXTitle("Transverse path length, L [fm]");
1406   hEll->SetYTitle("Probability");
1407   hEll->SetFillColor(2);
1408   for(Int_t i=0; i<n; i++) {
1409     GetLength(ell);
1410     hEll->Fill(ell);
1411   }
1412   hEll->Scale(1/(Double_t)n);
1413   TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1414   cL->cd();
1415   hEll->Draw();
1416
1417   if(save) {
1418     TFile *f = new TFile(fname,"recreate");
1419     hEll->Write();
1420     f->Close();
1421   }
1422   return;
1423 }
1424
1425 void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
1426 {
1427   //
1428   // Plot lengths back-to-back distributions
1429   //
1430   Double_t ell1,ell2;
1431   TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15); 
1432   hElls->SetXTitle("Transverse path length, L [fm]");
1433   hElls->SetYTitle("Transverse path length, L [fm]");
1434   for(Int_t i=0; i<n; i++) {
1435     GetLengthsBackToBack(ell1,ell2);
1436     hElls->Fill(ell1,ell2);
1437   }
1438   hElls->Scale(1/(Double_t)n);
1439   TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1440   gStyle->SetPalette(1,0);
1441   cLs->cd();
1442   hElls->Draw("col,Z");
1443   if(save) {
1444     TFile *f = new TFile(fname,"recreate");
1445     hElls->Write();
1446     f->Close();
1447   }
1448   return;
1449 }
1450
1451 void AliFastGlauber::PlotAlmonds() const
1452 {
1453   //
1454   // Plot almonds for some impact parameters
1455   //
1456   TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1457   gStyle->SetPalette(1,0);
1458   c->Divide(2,2);
1459   c->cd(1);
1460   fgWAlmondFixedB[0].Draw("cont1");
1461   c->cd(2);
1462   fgWAlmondFixedB[10].Draw("cont1");
1463   c->cd(3);
1464   fgWAlmondFixedB[20].Draw("cont1");
1465   c->cd(4);
1466   fgWAlmondFixedB[30].Draw("cont1");
1467   return;
1468 }
1469
1470 //=================== Added by A. Dainese 05/03/04 ===========================
1471
1472 void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
1473                                    Double_t b,Double_t x0,Double_t y0,
1474                                    Double_t phi0,Double_t ellCut) const
1475 {
1476   // 
1477   // Calculate integrals: 
1478   //  integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1479   //  integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
1480   //
1481   // for a parton with production point (x0,y0)
1482   // and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
1483   // in a collision with impact parameter b
1484   // 
1485
1486   // number of steps in l
1487   const Int_t    kNp  = 100;
1488   const Double_t kDl  = fgBMax/Double_t(kNp);
1489     
1490   // Initial radius
1491   const Double_t kr0  = TMath::Sqrt(x0 * x0 + y0 * y0);
1492   const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1493     
1494   // Radial steps
1495   Double_t l  = 0.;
1496   integral0 = 0.;
1497   integral1 = 0.;
1498   Int_t i = 0;
1499   while((i < knps) && (l < ellCut)) {
1500     // Transform into target frame
1501     const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1502     const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1503     const Double_t kphi  = TMath::ATan2(kyy, kxx);
1504     const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1505     // Radius in projectile frame
1506     const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi)); 
1507     const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1508     integral0 += kprodTATB * kDl;
1509     integral1 += kprodTATB * l * kDl;
1510     l  += kDl;
1511     i++;
1512   } // steps
1513   return;  
1514 }
1515
1516 void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1517                                    Double_t& phi,
1518                                    Double_t ellCut,Double_t b)
1519 {
1520   //
1521   // Return I0 and I1 from random b, x0, y0, phi0 
1522   // Return also phi
1523   //
1524   Double_t x0,y0,phi0;
1525   if(b<0.) GetRandomBHard(b);
1526   GetRandomXY(x0,y0);
1527   GetRandomPhi(phi0);
1528   phi = phi0;
1529   CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1530   return;
1531 }
1532
1533 void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1534                              Double_t ellCut,Double_t b)
1535 {
1536   //
1537   // Return I0 and I1 from random b, x0, y0, phi0 
1538   //
1539   Double_t phi;
1540   GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1541   return;
1542 }
1543
1544 void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1545                                              Double_t& integral02,Double_t& integral12,
1546                                              Double_t& phi,
1547                                              Double_t ellCut,Double_t b)
1548 {
1549   //
1550   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 
1551   // Return also phi0
1552   //
1553   Double_t x0,y0,phi0;
1554   if(b<0.) GetRandomBHard(b);
1555   GetRandomXY(x0,y0);
1556   GetRandomPhi(phi0);
1557   phi = phi0;
1558   const Double_t kphi0plusPi = phi0+TMath::Pi();
1559   CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1560   CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1561   return;
1562 }
1563
1564 void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11,
1565                                                   Double_t& integral02,Double_t& integral12,
1566                                                   Double_t& phi,Double_t &x,Double_t &y,
1567                                                   Double_t ellCut,Double_t b)
1568 {
1569   //
1570   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 
1571   // Return also phi0
1572   //
1573   Double_t x0,y0,phi0;
1574   if(b<0.) GetRandomBHard(b);
1575   GetRandomXY(x0,y0);
1576   GetRandomPhi(phi0);
1577   phi = phi0; x=x0; y=y0;
1578   const Double_t kphi0plusPi = phi0+TMath::Pi();
1579   CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1580   CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1581   return;
1582 }
1583
1584 void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1585                                        Double_t& integral02,Double_t& integral12,
1586                                        Double_t ellCut,Double_t b)
1587 {
1588   //
1589   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 
1590   //
1591   Double_t phi;
1592   GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1593                           phi,ellCut,b);
1594   return;
1595 }
1596
1597 void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1598                                       Double_t* integral0,Double_t* integral1,
1599                                       Double_t ellCut,Double_t b)
1600 {
1601   //
1602   // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n] 
1603   // from random b, x0, y0
1604   //
1605   Double_t x0,y0;
1606   if(b<0.) GetRandomBHard(b);
1607   GetRandomXY(x0,y0);
1608   for(Int_t i=0; i<n; i++) 
1609     CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1610   return;
1611 }
1612
1613 void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi,
1614                                       Double_t* integral0,Double_t* integral1,
1615                                       Double_t &x,Double_t& y,
1616                                       Double_t ellCut,Double_t b)
1617 {
1618   //
1619   // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n] 
1620   // from random b, x0, y0 and return x0,y0
1621   //
1622   Double_t x0,y0;
1623   if(b<0.) GetRandomBHard(b);
1624   GetRandomXY(x0,y0);
1625   for(Int_t i=0; i<n; i++) 
1626     CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1627   x=x0;
1628   y=y0;
1629   return;
1630 }
1631
1632 void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
1633                                    Bool_t save,Char_t *fname)
1634 {
1635   //
1636   // Plot I0-I1 distribution
1637   //
1638   Double_t i0,i1;
1639   TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01); 
1640   hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1641   hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1642
1643   TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
1644                        1000,0,0.001); 
1645   hI0->SetXTitle("I_{0} [fm^{-3}]");
1646   hI0->SetYTitle("Probability");
1647   hI0->SetFillColor(3);
1648   TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
1649                        1000,0,0.01); 
1650   hI1->SetXTitle("I_{1} [fm^{-2}]");
1651   hI1->SetYTitle("Probability");
1652   hI1->SetFillColor(4);
1653   TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1654                       100,0,0.02); 
1655   h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1656   h2->SetYTitle("Probability");
1657   h2->SetFillColor(2);
1658   TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1659                       100,0,15); 
1660   h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1661   h3->SetYTitle("Probability");
1662   h3->SetFillColor(5);
1663   TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1664                       100,0,0.00015); 
1665   h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1666   h4->SetYTitle("Probability");
1667   h4->SetFillColor(7);
1668
1669   for(Int_t i=0; i<n; i++) {
1670     GetI0I1(i0,i1,ellCut);
1671     hI0I1s->Fill(i0,i1);
1672     hI0->Fill(i0);
1673     hI1->Fill(i1);
1674     h2->Fill(2.*i1*i1/i0);
1675     h3->Fill(2.*i1/i0);
1676     h4->Fill(i0*i0/2./i1);
1677   }
1678   hI0->Scale(1/(Double_t)n);
1679   hI1->Scale(1/(Double_t)n);
1680   h2->Scale(1/(Double_t)n);
1681   h3->Scale(1/(Double_t)n);
1682   h4->Scale(1/(Double_t)n);
1683   hI0I1s->Scale(1/(Double_t)n);
1684
1685   TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1686   cI0I1->Divide(3,2);
1687   cI0I1->cd(1);
1688   hI0->Draw();
1689   cI0I1->cd(2);
1690   hI1->Draw();
1691   cI0I1->cd(3);
1692   h2->Draw();
1693   cI0I1->cd(4);
1694   h3->Draw();
1695   cI0I1->cd(5);
1696   h4->Draw();
1697   cI0I1->cd(6);
1698   gStyle->SetPalette(1,0);
1699   hI0I1s->Draw("col,Z");
1700
1701   if(save) {
1702     TFile *f = new TFile(fname,"recreate");
1703     hI0I1s->Write();
1704     hI0->Write();
1705     hI1->Write();
1706     h2->Write();
1707     h3->Write();
1708     h4->Write();
1709     f->Close();
1710   }
1711   return;
1712 }
1713
1714 void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
1715                                       Bool_t save,Char_t *fname)
1716 {
1717   //
1718   // Plot I0-I1 back-to-back distributions
1719   //
1720   Double_t i01,i11,i02,i12;
1721   TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100); 
1722   hI0s->SetXTitle("I_{0} [fm^{-3}]");
1723   hI0s->SetYTitle("I_{0} [fm^{-3}]");
1724   TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100); 
1725   hI1s->SetXTitle("I_{1} [fm^{-2}]");
1726   hI1s->SetYTitle("I_{1} [fm^{-2}]");
1727
1728   for(Int_t i=0; i<n; i++) {
1729     GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1730     hI0s->Fill(i01,i02);
1731     hI1s->Fill(i11,i12);
1732   }
1733   hI0s->Scale(1/(Double_t)n);
1734   hI1s->Scale(1/(Double_t)n);
1735
1736   TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
1737   gStyle->SetPalette(1,0);
1738   cI0I1s->Divide(2,1);
1739   cI0I1s->cd(1);
1740   hI0s->Draw("col,Z");
1741   cI0I1s->cd(2);
1742   hI1s->Draw("col,Z");
1743
1744   if(save) {
1745     TFile *f = new TFile(fname,"recreate");
1746     hI0s->Write();
1747     hI1s->Write();
1748     f->Close();
1749   }
1750   return;
1751 }
1752