Various improvements
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.cxx
1 // 
2 // Utilities used in the forward multiplcity analysis 
3 // 
4 //
5 #include "AliForwardUtil.h"
6 #include <AliAnalysisManager.h>
7 #include "AliAODForwardMult.h"
8 #include <AliLog.h>
9 #include <AliInputEventHandler.h>
10 #include <AliESDEvent.h>
11 #include <AliPhysicsSelection.h>
12 #include <AliTriggerAnalysis.h>
13 #include <AliMultiplicity.h>
14 #include <TH2D.h>
15 #include <TH1I.h>
16 #include <TF1.h>
17 #include <TFitResult.h>
18 #include <TMath.h>
19 #include <TError.h>
20
21 //====================================================================
22 UShort_t
23 AliForwardUtil::ParseCollisionSystem(const char* sys)
24 {
25   // 
26   // Parse a collision system spec given in a string.   Known values are 
27   // 
28   //  - "pp", "p-p" which returns kPP 
29   //  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb 
30   //  - Everything else gives kUnknown 
31   // 
32   // Parameters:
33   //    sys Collision system spec 
34   // 
35   // Return:
36   //    Collision system id 
37   //
38   TString s(sys);
39   s.ToLower();
40   if (s.Contains("p-p")   || s.Contains("pp"))   return AliForwardUtil::kPP; 
41   if (s.Contains("pb-pb") || s.Contains("pbpb")) return AliForwardUtil::kPbPb;
42   if (s.Contains("a-a")   || s.Contains("aa"))   return AliForwardUtil::kPbPb;
43   return AliForwardUtil::kUnknown;
44 }
45 //____________________________________________________________________
46 const char*
47 AliForwardUtil::CollisionSystemString(UShort_t sys)
48 {
49   // 
50   // Get a string representation of the collision system 
51   // 
52   // Parameters:
53   //    sys  Collision system 
54   // - kPP -> "pp"
55   // - kPbPb -> "PbPb" 
56   // - anything else gives "unknown"
57   // 
58   // Return:
59   //    String representation of the collision system 
60   //
61   switch (sys) { 
62   case AliForwardUtil::kPP:   return "pp";
63   case AliForwardUtil::kPbPb: return "PbPb";
64   }
65   return "unknown";
66 }
67 //____________________________________________________________________
68 UShort_t
69 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t /* sys */, Float_t v)
70 {
71   // 
72   // Parse the center of mass energy given as a float and return known 
73   // values as a unsigned integer
74   // 
75   // Parameters:
76   //    sys  Collision system (needed for AA)
77   //    cms  Center of mass energy * total charge 
78   // 
79   // Return:
80   //    Center of mass energy per nucleon
81   //
82   Float_t energy = v; 
83   // Below no longer needed apparently
84   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
85   if (TMath::Abs(energy - 900.)   < 10)  return 900;
86   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
87   if (TMath::Abs(energy - 2750.)  < 10)  return 2750;
88   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
89   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
90   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
91   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
92   return 0;
93 }
94 //____________________________________________________________________
95 const char* 
96 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
97 {
98   // 
99   // Get a string representation of the center of mass energy per nuclean
100   // 
101   // Parameters:
102   //    cms  Center of mass energy per nucleon
103   // 
104   // Return:
105   //    String representation of the center of mass energy per nuclean
106   //
107   return Form("%04dGeV", cms);
108 }
109 //____________________________________________________________________
110 Short_t
111 AliForwardUtil::ParseMagneticField(Float_t v)
112 {
113   // 
114   // Parse the magnetic field (in kG) as given by a floating point number
115   // 
116   // Parameters:
117   //    field  Magnetic field in kG 
118   // 
119   // Return:
120   //    Short integer value of magnetic field in kG 
121   //
122   if (TMath::Abs(v - 5.) < 1 ) return +5;
123   if (TMath::Abs(v + 5.) < 1 ) return -5;
124   if (TMath::Abs(v) < 1)       return 0;
125   return 999;
126 }
127 //____________________________________________________________________
128 const char* 
129 AliForwardUtil::MagneticFieldString(Short_t f)
130 {
131   // 
132   // Get a string representation of the magnetic field
133   // 
134   // Parameters:
135   //    field Magnetic field in kG
136   // 
137   // Return:
138   //    String representation of the magnetic field
139   //
140   return Form("%01dkG", f);
141 }
142
143
144 //====================================================================
145 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
146 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
147 namespace {
148   // 
149   // The shift of the most probable value for the ROOT function TMath::Landau 
150   //
151   const Double_t  mpshift  = -0.22278298;
152   // 
153   // Integration normalisation 
154   //
155   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
156
157   // 
158   // Utility function to use in TF1 defintition 
159   //
160   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
161   {
162     Double_t x        = xp[0];
163     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
164     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
165     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
166     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
167     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
168
169     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
170   }
171
172   // 
173   // Utility function to use in TF1 defintition 
174   //
175   Double_t landauGausN(Double_t* xp, Double_t* pp) 
176   {
177     Double_t  x        = xp[0];
178     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
179     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
180     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
181     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
182     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
183     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
184     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
185
186     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
187                                                   n, a);
188   }
189   // 
190   // Utility function to use in TF1 defintition 
191   //
192   Double_t landauGausI(Double_t* xp, Double_t* pp) 
193   {
194     Double_t x         = xp[0];
195     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
196     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
197     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
198     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
199     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
200     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
201
202     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
203   }
204
205
206 }
207 //____________________________________________________________________
208 Double_t 
209 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
210 {
211   // 
212   // Calculate the shifted Landau
213   // @f[
214   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
215   // @f]
216   //
217   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
218   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
219   // @f$\Delta=0,\xi=1@f$. 
220   // 
221   // Parameters:
222   //    x      Where to evaluate @f$ f'_{L}@f$ 
223   //    delta  Most probable value 
224   //    xi     The 'width' of the distribution 
225   //
226   // Return:
227   //    @f$ f'_{L}(x;\Delta,\xi) @f$
228   //
229   return TMath::Landau(x, delta - xi * mpshift, xi);
230 }
231 //____________________________________________________________________
232 Double_t 
233 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
234                            Double_t sigma, Double_t sigmaN)
235 {
236   // 
237   // Calculate the value of a Landau convolved with a Gaussian 
238   // 
239   // @f[ 
240   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
241   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
242   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
243   // @f]
244   // 
245   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
246   // energy loss, @f$ \xi@f$ the width of the Landau, and 
247   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
248   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
249   // noise in the detector.  
250   //
251   // Note that this function uses the constants fgConvolutionSteps and
252   // fgConvolutionNSigma
253   // 
254   // References: 
255   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
256   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
257   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
258   // 
259   // Parameters:
260   //    x         where to evaluate @f$ f@f$
261   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
262   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
263   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
264   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
265   // 
266   // Return:
267   //    @f$ f@f$ evaluated at @f$ x@f$.  
268   //
269   Double_t deltap = delta - xi * mpshift;
270   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
271   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
272   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
273   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
274   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
275   Double_t sum    = 0;
276   
277   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
278     Double_t x1 = xlow  + (i - .5) * step;
279     Double_t x2 = xhigh - (i - .5) * step;
280     
281     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
282     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
283   }
284   return step * sum * invSq2pi / sigma1;
285 }
286
287 //____________________________________________________________________
288 Double_t 
289 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
290                             Double_t sigma, Double_t sigmaN, Int_t i)
291 {
292   // 
293   // Evaluate 
294   // @f[ 
295   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
296   // @f] 
297   // corresponding to @f$ i@f$ particles i.e., with the substitutions 
298   // @f{eqnarray*}{ 
299   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
300   //    \xi       \rightarrow \xi_i       &=& i \xi
301   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
302   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
303   // @f} 
304   // 
305   // Parameters:
306   //    x        Where to evaluate 
307   //    delta    @f$ \Delta@f$ 
308   //    xi       @f$ \xi@f$ 
309   //    sigma    @f$ \sigma@f$ 
310   //    sigma_n  @f$ \sigma_n@f$
311   //    i        @f$ i @f$
312   // 
313   // Return:
314   //    @f$ f_i @f$ evaluated
315   //  
316   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
317   Double_t xiI    =  i * xi;
318   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
319   if (sigmaI < 1e-10) { 
320     // Fall back to landau 
321     return AliForwardUtil::Landau(x, deltaI, xiI);
322   }
323   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
324 }
325
326 //____________________________________________________________________
327 Double_t 
328 AliForwardUtil::IdLandauGausdPar(Double_t x, 
329                                  UShort_t par,   Double_t dPar, 
330                                  Double_t delta, Double_t xi, 
331                                  Double_t sigma, Double_t sigmaN, 
332                                  Int_t    i)
333 {
334   // 
335   // Numerically evaluate 
336   // @f[ 
337   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
338   // @f] 
339   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
340   // of the parameters is given by 
341   //
342   // - 0: @f$\Delta@f$ 
343   // - 1: @f$\xi@f$ 
344   // - 2: @f$\sigma@f$ 
345   // - 3: @f$\sigma_n@f$ 
346   //
347   // This is the partial derivative with respect to the parameter of
348   // the response function corresponding to @f$ i@f$ particles i.e.,
349   // with the substitutions
350   // @f[ 
351   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
352   //    \xi       \rightarrow \xi_i       = i \xi
353   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
354   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
355   // @f] 
356   // 
357   // Parameters:
358   //    x        Where to evaluate 
359   //    ipar     Parameter number 
360   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
361   //    delta    @f$ \Delta@f$ 
362   //    xi       @f$ \xi@f$ 
363   //    sigma    @f$ \sigma@f$ 
364   //    sigma_n  @f$ \sigma_n@f$
365   //    i        @f$ i@f$
366   // 
367   // Return:
368   //    @f$ f_i@f$ evaluated
369   //  
370   if (dPar == 0) return 0;
371   Double_t dp      = dPar;
372   Double_t d2      = dPar / 2;
373   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
374   Double_t xiI     =  i * xi;
375   Double_t si      =  TMath::Sqrt(Double_t(i));
376   Double_t sigmaI  =  si*sigma;
377   Double_t y1      = 0;
378   Double_t y2      = 0;
379   Double_t y3      = 0;
380   Double_t y4      = 0;
381   switch (par) {
382   case 0: 
383     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
384     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
385     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
386     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
387     break;
388   case 1: 
389     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
390     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
391     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
392     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
393     break;
394   case 2: 
395     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
396     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
397     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
398     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
399     break;
400   case 3: 
401     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
402     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
403     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
404     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
405     break;
406   default:
407     return 0;
408   } 
409   
410   Double_t d0  = y1 - y4;
411   Double_t d1  = 2 * (y2 - y3);
412   
413   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
414    
415   return g;
416 }
417
418 //____________________________________________________________________
419 Double_t 
420 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
421                             Double_t sigma, Double_t sigmaN, Int_t n, 
422                             Double_t* a)
423 {
424   // 
425   // Evaluate 
426   // @f[ 
427   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
428   // @f] 
429   // 
430   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
431   // Landau with a Gaussian (see LandauGaus).  Note that 
432   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
433   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
434   //  
435   // References: 
436   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
437   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
438   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
439   // 
440   // Parameters:
441   //    x        Where to evaluate @f$ f_N@f$
442   //    delta    @f$ \Delta_1@f$ 
443   //    xi       @f$ \xi_1@f$
444   //    sigma    @f$ \sigma_1@f$ 
445   //    sigma_n  @f$ \sigma_n@f$ 
446   //    n        @f$ N@f$ in the sum above.
447   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
448   //                 @f$ i > 1@f$ 
449   // 
450   // Return:
451   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
452   //
453   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
454   for (Int_t i = 2; i <= n; i++) 
455     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
456   return result;
457 }
458 namespace { 
459   const Int_t kColors[] = { kRed+1, 
460                             kPink+3, 
461                             kMagenta+2, 
462                             kViolet+2, 
463                             kBlue+1, 
464                             kAzure+3, 
465                             kCyan+1, 
466                             kTeal+2, 
467                             kGreen+2, 
468                             kSpring+3, 
469                             kYellow+2, 
470                             kOrange+2 };
471 }
472
473 //____________________________________________________________________
474 TF1*
475 AliForwardUtil::MakeNLandauGaus(Double_t  c, 
476                                 Double_t  delta, Double_t xi, 
477                                 Double_t  sigma, Double_t sigmaN, Int_t n, 
478                                 Double_t* a, 
479                                 Double_t  xmin, Double_t xmax)
480 {
481   // 
482   // Generate a TF1 object of @f$ f_N@f$ 
483   // 
484   // Parameters:
485   //    c         Constant                             
486   //    delta     @f$ \Delta@f$                        
487   //    xi            @f$ \xi_1@f$                     
488   //    sigma     @f$ \sigma_1@f$                      
489   //    sigma_n   @f$ \sigma_n@f$                      
490   //    n             @f$ N@f$ - how many particles to sum to
491   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
492   //                  @f$ i > 1@f$ 
493   //    xmin      Least value of range  
494   //    xmax      Largest value of range
495   // 
496   // Return:
497   //    Newly allocated TF1 object
498   //
499   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
500   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
501   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
502   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
503   landaun->SetLineWidth(2);
504   landaun->SetNpx(500);
505   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
506
507   // Set the initial parameters from the seed fit 
508   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
509   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
510   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
511   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
512   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
513   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);       
514
515   // Set the range and name of the scale parameters 
516   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
517     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
518     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
519   }
520   return landaun;
521 }
522 //____________________________________________________________________
523 TF1*
524 AliForwardUtil::MakeILandauGaus(Double_t  c, 
525                                 Double_t  delta, Double_t xi, 
526                                 Double_t  sigma, Double_t sigmaN, Int_t i, 
527                                 Double_t  xmin, Double_t xmax)
528 {
529   // 
530   // Generate a TF1 object of @f$ f_I@f$ 
531   // 
532   // Parameters:
533   //    c        Constant
534   //    delta    @f$ \Delta@f$ 
535   //    xi       @f$ \xi_1@f$          
536   //    sigma    @f$ \sigma_1@f$               
537   //    sigma_n  @f$ \sigma_n@f$               
538   //    i            @f$ i@f$ - the number of particles
539   //    xmin     Least value of range
540   //    xmax     Largest value of range
541   // 
542   // Return:
543   //    Newly allocated TF1 object
544   //
545   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
546   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
547   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
548   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
549   landaui->SetLineWidth(1);
550   landaui->SetNpx(500);
551   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
552
553   // Set the initial parameters from the seed fit 
554   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
555   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
556   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
557   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
558   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
559   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);       
560
561   return landaui;
562 }
563
564 //====================================================================
565 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
566                                          Double_t maxRange, 
567                                          UShort_t minusBins) 
568   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
569     fFitResults(0), fFunctions(0)
570 {
571   // 
572   // Constructor 
573   // 
574   // Parameters:
575   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
576   //    maxRange   Maximum range to fit to 
577   //    minusBins  The number of bins below maximum to use 
578   //
579   fFitResults.SetOwner();
580   fFunctions.SetOwner();
581 }
582 //____________________________________________________________________
583 AliForwardUtil::ELossFitter::~ELossFitter()
584 {
585   // 
586   // Destructor
587   // 
588   //
589   fFitResults.Delete();
590   fFunctions.Delete();
591 }
592 //____________________________________________________________________
593 void
594 AliForwardUtil::ELossFitter::Clear()
595 {
596   // 
597   // Clear internal arrays 
598   // 
599   //
600   fFitResults.Clear();
601   fFunctions.Clear();
602 }
603 //____________________________________________________________________
604 TF1*
605 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
606 {
607   // 
608   // Fit a 1-particle signal to the passed energy loss distribution 
609   // 
610   // Note that this function clears the internal arrays first 
611   // 
612   // Parameters:
613   //    dist    Data to fit the function to 
614   //    sigman If larger than zero, the initial guess of the
615   //               detector induced noise. If zero or less, then this 
616   //               parameter is ignored in the fit (fixed at 0)
617   // 
618   // Return:
619   //    The function fitted to the data 
620   //
621
622   // Clear the cache 
623   Clear();
624   
625   // Find the fit range 
626   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
627   
628   // Get the bin with maximum 
629   Int_t    maxBin = dist->GetMaximumBin();
630   Double_t maxE   = dist->GetBinLowEdge(maxBin);
631   
632   // Get the low edge 
633   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
634   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
635   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
636   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
637
638   // Restore the range 
639   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
640   
641   // Define the function to fit 
642   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
643
644   // Set initial guesses, parameter names, and limits  
645   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
646   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
647   landau1->SetNpx(500);
648   landau1->SetParLimits(kDelta, minE, fMaxRange);
649   landau1->SetParLimits(kXi,    0.00, fMaxRange);
650   landau1->SetParLimits(kSigma, 0.01, 0.1);
651   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
652   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
653
654   // Do the fit, getting the result object 
655   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
656   landau1->SetRange(minE, fMaxRange);
657   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
658   fFunctions.AddAtAndExpand(landau1, 0);
659
660   return landau1;
661 }
662 //____________________________________________________________________
663 TF1*
664 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
665                                           Double_t sigman)
666 {
667   // 
668   // Fit a N-particle signal to the passed energy loss distribution 
669   //
670   // If there's no 1-particle fit present, it does that first 
671   // 
672   // Parameters:
673   //    dist   Data to fit the function to 
674   //    n      Number of particle signals to fit 
675   //    sigman If larger than zero, the initial guess of the
676   //               detector induced noise. If zero or less, then this 
677   //               parameter is ignored in the fit (fixed at 0)
678   // 
679   // Return:
680   //    The function fitted to the data 
681   //
682
683   // Get the seed fit result 
684   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
685   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
686   if (!r || !f) { 
687     f = Fit1Particle(dist, sigman);
688     r = static_cast<TFitResult*>(fFitResults.At(0));
689     if (!r || !f) { 
690       ::Warning("FitNLandau", "No first shot at landau fit");
691       return 0;
692     }
693   }
694
695   // Get some parameters from seed fit 
696   Double_t delta1  = r->Parameter(kDelta);
697   Double_t xi1     = r->Parameter(kXi);
698   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
699   Double_t minE    = f->GetXmin();
700
701   // Array of weights 
702   TArrayD a(n-1);
703   for (UShort_t i = 2; i <= n; i++) 
704     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
705   // Make the fit function 
706   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
707                                      r->Parameter(kDelta),
708                                      r->Parameter(kXi),
709                                      r->Parameter(kSigma),
710                                      r->Parameter(kSigmaN),
711                                      n,a.fArray,minE,maxEi);
712   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
713   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
714   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
715   // Check if we're using the noise sigma 
716   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
717   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
718
719   // Set the range and name of the scale parameters 
720   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
721     landaun->SetParLimits(kA+i-2, 0,1);
722   }
723
724   // Do the fit 
725   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
726   
727   landaun->SetRange(minE, fMaxRange);
728   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
729   fFunctions.AddAtAndExpand(landaun, n-1);
730   
731   return landaun;
732 }  
733
734 //====================================================================
735 AliForwardUtil::Histos::~Histos()
736 {
737   // 
738   // Destructor
739   //
740   if (fFMD1i) delete fFMD1i;
741   if (fFMD2i) delete fFMD2i;
742   if (fFMD2o) delete fFMD2o;
743   if (fFMD3i) delete fFMD3i;
744   if (fFMD3o) delete fFMD3o;
745 }
746
747 //____________________________________________________________________
748 TH2D*
749 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, 
750                              const TAxis& etaAxis) const
751 {
752   // 
753   // Make a histogram 
754   // 
755   // Parameters:
756   //    d        Detector
757   //    r        Ring 
758   //    etaAxis  Eta axis to use
759   // 
760   // Return:
761   //    Newly allocated histogram 
762   //
763   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
764   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r), 
765                         Form("FMD%d%c cache", d, r),
766                         etaAxis.GetNbins(), etaAxis.GetXmin(), 
767                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
768   hist->SetXTitle("#eta");
769   hist->SetYTitle("#phi [radians]");
770   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
771   hist->Sumw2();
772   hist->SetDirectory(0);
773
774   return hist;
775 }
776 //____________________________________________________________________
777 void
778 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
779 {
780   // 
781   // Initialize the object 
782   // 
783   // Parameters:
784   //    etaAxis Eta axis to use 
785   //
786   fFMD1i = Make(1, 'I', etaAxis);
787   fFMD2i = Make(2, 'I', etaAxis);
788   fFMD2o = Make(2, 'O', etaAxis);
789   fFMD3i = Make(3, 'I', etaAxis);
790   fFMD3o = Make(3, 'O', etaAxis);
791 }
792 //____________________________________________________________________
793 void
794 AliForwardUtil::Histos::Clear(Option_t* option)
795 {
796   // 
797   // Clear data 
798   // 
799   // Parameters:
800   //    option Not used 
801   //
802   fFMD1i->Reset(option);
803   fFMD2i->Reset(option);
804   fFMD2o->Reset(option);
805   fFMD3i->Reset(option);
806   fFMD3o->Reset(option);
807 }
808
809 //____________________________________________________________________
810 TH2D*
811 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
812 {
813   // 
814   // Get the histogram for a particular detector,ring
815   // 
816   // Parameters:
817   //    d Detector 
818   //    r Ring 
819   // 
820   // Return:
821   //    Histogram for detector,ring or nul 
822   //
823   switch (d) { 
824   case 1: return fFMD1i;
825   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
826   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
827   }
828   return 0;
829 }
830 //====================================================================
831 TList*
832 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
833 {
834   // 
835   // Define the outout list in @a d
836   // 
837   // Parameters:
838   //    d Where to put the output list
839   // 
840   // Return:
841   //    Newly allocated TList object or null
842   //
843   if (!d) return 0;
844   TList* list = new TList;
845   list->SetName(fName.Data());
846   d->Add(list);
847   return list;
848 }
849 //____________________________________________________________________
850 TList*
851 AliForwardUtil::RingHistos::GetOutputList(TList* d) const
852 {
853   // 
854   // Get our output list from the container @a d
855   // 
856   // Parameters:
857   //    d where to get the output list from 
858   // 
859   // Return:
860   //    The found TList or null
861   //
862   if (!d) return 0;
863   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
864   return list;
865 }
866
867 //____________________________________________________________________
868 TH1*
869 AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
870 {
871   // 
872   // Find a specific histogram in the source list @a d
873   // 
874   // Parameters:
875   //    d     (top)-container 
876   //    name  Name of histogram
877   // 
878   // Return:
879   //    Found histogram or null
880   //
881   return static_cast<TH1*>(d->FindObject(name));
882 }
883
884 //
885 // EOF
886 //