Fixes for pA indenfication of events
[u/mrichter/AliRoot.git] / PWGLF / 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 <AliAODInputHandler.h>
11 #include <AliAODHandler.h>
12 #include <AliAODEvent.h>
13 #include <AliESDEvent.h>
14 #include <AliAnalysisTaskSE.h>
15 #include <AliPhysicsSelection.h>
16 #include <AliTriggerAnalysis.h>
17 #include <AliMultiplicity.h>
18 #include <TParameter.h>
19 #include <TH2D.h>
20 #include <TH1I.h>
21 #include <TF1.h>
22 #include <TFitResult.h>
23 #include <TMath.h>
24 #include <TError.h>
25 #include <TROOT.h>
26
27 //====================================================================
28 UShort_t
29 AliForwardUtil::ParseCollisionSystem(const char* sys)
30 {
31   // 
32   // Parse a collision system spec given in a string.   Known values are 
33   // 
34   //  - "ppb", "p-pb", "pa", "p-a"  which returns kPPb
35   //  - "pp", "p-p"                 which returns kPP 
36   //  - "PbPb", "Pb-Pb", "A-A",     which returns kPbPb 
37   //  - Everything else gives kUnknown 
38   // 
39   // Parameters:
40   //    sys Collision system spec 
41   // 
42   // Return:
43   //    Collision system id 
44   //
45   TString s(sys);
46   s.ToLower();
47   // we do pA first to avoid pp catch on ppb string (AH)
48   if (s.Contains("p-pb")  || s.Contains("ppb"))   return AliForwardUtil::kPPb;
49   if (s.Contains("p-a")   || s.Contains("pa"))    return AliForwardUtil::kPPb;
50   if (s.Contains("a-p")   || s.Contains("ap"))    return AliForwardUtil::kPPb;
51   if (s.Contains("p-p")   || s.Contains("pp"))    return AliForwardUtil::kPP; 
52   if (s.Contains("pb-pb") || s.Contains("pbpb"))  return AliForwardUtil::kPbPb;
53   if (s.Contains("a-a")   || s.Contains("aa"))    return AliForwardUtil::kPbPb;
54   return AliForwardUtil::kUnknown;
55 }
56 //____________________________________________________________________
57 const char*
58 AliForwardUtil::CollisionSystemString(UShort_t sys)
59 {
60   // 
61   // Get a string representation of the collision system 
62   // 
63   // Parameters:
64   //    sys  Collision system 
65   // - kPP -> "pp"
66   // - kPbPb -> "PbPb" 
67   // - anything else gives "unknown"
68   // 
69   // Return:
70   //    String representation of the collision system 
71   //
72   switch (sys) { 
73   case AliForwardUtil::kPP:   return "pp";
74   case AliForwardUtil::kPbPb: return "PbPb";
75   case AliForwardUtil::kPPb:  return "pPb";
76   }
77   return "unknown";
78 }
79 //____________________________________________________________________
80 UShort_t
81 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t /* sys */, Float_t v)
82 {
83   // 
84   // Parse the center of mass energy given as a float and return known 
85   // values as a unsigned integer
86   // 
87   // Parameters:
88   //    sys  Collision system (needed for AA)
89   //    cms  Center of mass energy * total charge 
90   // 
91   // Return:
92   //    Center of mass energy per nucleon
93   //
94   Float_t energy = v; 
95   // Below no longer needed apparently
96   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
97   if (TMath::Abs(energy - 900.)   < 10)  return 900;
98   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
99   if (TMath::Abs(energy - 2750.)  < 20)  return 2750;
100   if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
101   if (TMath::Abs(energy - 5022.)  < 10)  return 5000;
102   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
103   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
104   if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
105   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
106   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
107   return 0;
108 }
109 //____________________________________________________________________
110 const char* 
111 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
112 {
113   // 
114   // Get a string representation of the center of mass energy per nuclean
115   // 
116   // Parameters:
117   //    cms  Center of mass energy per nucleon
118   // 
119   // Return:
120   //    String representation of the center of mass energy per nuclean
121   //
122   return Form("%04dGeV", cms);
123 }
124 //____________________________________________________________________
125 Short_t
126 AliForwardUtil::ParseMagneticField(Float_t v)
127 {
128   // 
129   // Parse the magnetic field (in kG) as given by a floating point number
130   // 
131   // Parameters:
132   //    field  Magnetic field in kG 
133   // 
134   // Return:
135   //    Short integer value of magnetic field in kG 
136   //
137   if (TMath::Abs(v - 5.) < 1 ) return +5;
138   if (TMath::Abs(v + 5.) < 1 ) return -5;
139   if (TMath::Abs(v) < 1)       return 0;
140   return 999;
141 }
142 //____________________________________________________________________
143 const char* 
144 AliForwardUtil::MagneticFieldString(Short_t f)
145 {
146   // 
147   // Get a string representation of the magnetic field
148   // 
149   // Parameters:
150   //    field Magnetic field in kG
151   // 
152   // Return:
153   //    String representation of the magnetic field
154   //
155   return Form("%01dkG", f);
156 }
157 //_____________________________________________________________________
158 AliAODEvent* AliForwardUtil::GetAODEvent(AliAnalysisTaskSE* task)
159 {
160   // Check if AOD is the output event
161   if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
162
163   AliAODEvent* ret = task->AODEvent();
164   if (ret) return ret; 
165   
166   // Check if AOD is the input event 
167   ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
168   if (!ret) ::Warning("GetAODEvent", "No AOD event found");
169   
170   return ret; 
171 }
172 //_____________________________________________________________________
173 UShort_t AliForwardUtil::CheckForAOD()
174 {
175   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
176   if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
177     ::Info("CheckForAOD", "Found AOD Input handler");
178     return 1;
179   }
180   if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
181     ::Info("CheckForAOD", "Found AOD Output handler");
182     return 2;
183   }
184
185   ::Warning("CheckForAOD", 
186             "Neither and input nor output AOD handler is specified");
187   return 0;
188 }
189 //_____________________________________________________________________
190 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
191 {
192   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
193   if (!cls) { 
194     AliAnalysisTask* t = am->GetTask(clsOrName);
195     if (!t) { 
196       ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
197       return false;
198     }
199     ::Info("CheckForTask", "Found task %s", clsOrName);
200     return true;
201   }
202   TClass* dep = gROOT->GetClass(clsOrName);
203   if (!dep) { 
204     ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
205     return false;
206   }
207   TIter next(am->GetTasks());
208   TObject* o = 0;
209   while ((o = next())) { 
210     if (o->IsA()->InheritsFrom(dep)) {
211       ::Info("CheckForTask", "Found task of class %s: %s", 
212              clsOrName, o->GetName());
213       return true;
214     }
215   }
216   ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
217   return false;
218 }
219
220 //_____________________________________________________________________
221 TObject* AliForwardUtil::MakeParameter(const Char_t* name, UShort_t value)
222 {
223   TParameter<int>* ret = new TParameter<int>(name, value);
224   ret->SetUniqueID(value);
225   return ret;
226 }
227 //_____________________________________________________________________
228 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Int_t value)
229 {
230   TParameter<int>* ret = new TParameter<int>(name, value);
231   ret->SetUniqueID(value);
232   return ret;
233 }
234 //_____________________________________________________________________
235 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Double_t value)
236 {
237   TParameter<double>* ret = new TParameter<double>(name, value);
238   Float_t v = value;
239   ret->SetUniqueID(*reinterpret_cast<UInt_t*>(&v));
240   return ret;
241 }
242 //_____________________________________________________________________
243 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
244 {
245   TParameter<bool>* ret = new TParameter<bool>(name, value);
246   ret->SetUniqueID(value);
247   return ret;
248 }
249
250 //_____________________________________________________________________
251 void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
252 {
253   if (!o) return;
254   value = o->GetUniqueID();
255 }
256 //_____________________________________________________________________
257 void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
258 {
259   if (!o) return;
260   value = o->GetUniqueID();
261 }
262 //_____________________________________________________________________
263 void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
264 {
265   if (!o) return;
266   UInt_t  i = o->GetUniqueID();
267   Float_t v = *reinterpret_cast<Float_t*>(&i);
268   value = v;
269 }
270 //_____________________________________________________________________
271 void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
272 {
273   if (!o) return;
274   value = o->GetUniqueID();
275 }
276   
277 //_____________________________________________________________________
278 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx)
279 {
280   //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta but support displaced vertices)
281   
282   //Get max R of ring
283   Double_t maxR = 0;
284   Double_t minR = 0;
285   Bool_t inner = false;
286   switch (ring) { 
287   case 'i': case 'I': maxR = 17.2; minR =  4.5213; inner = true;  break;
288   case 'o': case 'O': maxR = 28.0; minR = 15.4;    inner = false; break;
289   default: 
290     return -99999;
291   }
292
293   Double_t   rad       =  maxR- minR;
294   Double_t   nStrips   = (ring == 'I' ? 512 : 256);
295   Double_t   segment   = rad / nStrips;
296   Double_t   r         =  minR + segment*strip;
297   Int_t hybrid = sec / 2;
298   
299   Double_t z = 0;
300   switch (det) { 
301   case 1: z = 320.266; break;
302   case 2: z = (inner ? 83.666 : 74.966); break;
303   case 3: z = (inner ? -63.066 : -74.966); break; 
304   default: return -999999;
305   }
306   if ((hybrid % 2) == 0) z -= .5;
307   
308   Double_t   theta = TMath::ATan2(r,z-zvtx);
309   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
310   
311   return eta;
312 }
313
314
315 //====================================================================
316 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
317 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
318 namespace {
319   // 
320   // The shift of the most probable value for the ROOT function TMath::Landau 
321   //
322   const Double_t  mpshift  = -0.22278298;
323   // 
324   // Integration normalisation 
325   //
326   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
327
328   // 
329   // Utility function to use in TF1 defintition 
330   //
331   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
332   {
333     Double_t x        = xp[0];
334     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
335     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
336     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
337     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
338     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
339
340     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
341   }
342
343   // 
344   // Utility function to use in TF1 defintition 
345   //
346   Double_t landauGausN(Double_t* xp, Double_t* pp) 
347   {
348     Double_t  x        = xp[0];
349     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
350     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
351     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
352     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
353     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
354     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
355     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
356
357     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
358                                                   n, a);
359   }
360   // 
361   // Utility function to use in TF1 defintition 
362   //
363   Double_t landauGausI(Double_t* xp, Double_t* pp) 
364   {
365     Double_t x         = xp[0];
366     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
367     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
368     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
369     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
370     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
371     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
372
373     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
374   }
375
376
377 }
378 //____________________________________________________________________
379 Double_t 
380 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
381 {
382   // 
383   // Calculate the shifted Landau
384   // @f[
385   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
386   // @f]
387   //
388   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
389   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
390   // @f$\Delta=0,\xi=1@f$. 
391   // 
392   // Parameters:
393   //    x      Where to evaluate @f$ f'_{L}@f$ 
394   //    delta  Most probable value 
395   //    xi     The 'width' of the distribution 
396   //
397   // Return:
398   //    @f$ f'_{L}(x;\Delta,\xi) @f$
399   //
400   return TMath::Landau(x, delta - xi * mpshift, xi);
401 }
402 //____________________________________________________________________
403 Double_t 
404 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
405                            Double_t sigma, Double_t sigmaN)
406 {
407   // 
408   // Calculate the value of a Landau convolved with a Gaussian 
409   // 
410   // @f[ 
411   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
412   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
413   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
414   // @f]
415   // 
416   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
417   // energy loss, @f$ \xi@f$ the width of the Landau, and 
418   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
419   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
420   // noise in the detector.  
421   //
422   // Note that this function uses the constants fgConvolutionSteps and
423   // fgConvolutionNSigma
424   // 
425   // References: 
426   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
427   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
428   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
429   // 
430   // Parameters:
431   //    x         where to evaluate @f$ f@f$
432   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
433   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
434   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
435   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
436   // 
437   // Return:
438   //    @f$ f@f$ evaluated at @f$ x@f$.  
439   //
440   Double_t deltap = delta - xi * mpshift;
441   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
442   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
443   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
444   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
445   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
446   Double_t sum    = 0;
447   
448   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
449     Double_t x1 = xlow  + (i - .5) * step;
450     Double_t x2 = xhigh - (i - .5) * step;
451     
452     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
453     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
454   }
455   return step * sum * invSq2pi / sigma1;
456 }
457
458 //____________________________________________________________________
459 Double_t 
460 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
461                             Double_t sigma, Double_t sigmaN, Int_t i)
462 {
463   // 
464   // Evaluate 
465   // @f[ 
466   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
467   // @f] 
468   // corresponding to @f$ i@f$ particles i.e., with the substitutions 
469   // @f{eqnarray*}{ 
470   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
471   //    \xi       \rightarrow \xi_i       &=& i \xi
472   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
473   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
474   // @f} 
475   // 
476   // Parameters:
477   //    x        Where to evaluate 
478   //    delta    @f$ \Delta@f$ 
479   //    xi       @f$ \xi@f$ 
480   //    sigma    @f$ \sigma@f$ 
481   //    sigma_n  @f$ \sigma_n@f$
482   //    i        @f$ i @f$
483   // 
484   // Return:
485   //    @f$ f_i @f$ evaluated
486   //  
487   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
488   Double_t xiI    =  i * xi;
489   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
490   if (sigmaI < 1e-10) { 
491     // Fall back to landau 
492     return AliForwardUtil::Landau(x, deltaI, xiI);
493   }
494   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
495 }
496
497 //____________________________________________________________________
498 Double_t 
499 AliForwardUtil::IdLandauGausdPar(Double_t x, 
500                                  UShort_t par,   Double_t dPar, 
501                                  Double_t delta, Double_t xi, 
502                                  Double_t sigma, Double_t sigmaN, 
503                                  Int_t    i)
504 {
505   // 
506   // Numerically evaluate 
507   // @f[ 
508   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
509   // @f] 
510   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
511   // of the parameters is given by 
512   //
513   // - 0: @f$\Delta@f$ 
514   // - 1: @f$\xi@f$ 
515   // - 2: @f$\sigma@f$ 
516   // - 3: @f$\sigma_n@f$ 
517   //
518   // This is the partial derivative with respect to the parameter of
519   // the response function corresponding to @f$ i@f$ particles i.e.,
520   // with the substitutions
521   // @f[ 
522   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
523   //    \xi       \rightarrow \xi_i       = i \xi
524   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
525   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
526   // @f] 
527   // 
528   // Parameters:
529   //    x        Where to evaluate 
530   //    ipar     Parameter number 
531   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
532   //    delta    @f$ \Delta@f$ 
533   //    xi       @f$ \xi@f$ 
534   //    sigma    @f$ \sigma@f$ 
535   //    sigma_n  @f$ \sigma_n@f$
536   //    i        @f$ i@f$
537   // 
538   // Return:
539   //    @f$ f_i@f$ evaluated
540   //  
541   if (dPar == 0) return 0;
542   Double_t dp      = dPar;
543   Double_t d2      = dPar / 2;
544   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
545   Double_t xiI     =  i * xi;
546   Double_t si      =  TMath::Sqrt(Double_t(i));
547   Double_t sigmaI  =  si*sigma;
548   Double_t y1      = 0;
549   Double_t y2      = 0;
550   Double_t y3      = 0;
551   Double_t y4      = 0;
552   switch (par) {
553   case 0: 
554     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
555     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
556     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
557     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
558     break;
559   case 1: 
560     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
561     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
562     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
563     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
564     break;
565   case 2: 
566     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
567     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
568     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
569     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
570     break;
571   case 3: 
572     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
573     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
574     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
575     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
576     break;
577   default:
578     return 0;
579   } 
580   
581   Double_t d0  = y1 - y4;
582   Double_t d1  = 2 * (y2 - y3);
583   
584   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
585    
586   return g;
587 }
588
589 //____________________________________________________________________
590 Double_t 
591 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
592                             Double_t sigma, Double_t sigmaN, Int_t n, 
593                             const Double_t* a)
594 {
595   // 
596   // Evaluate 
597   // @f[ 
598   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
599   // @f] 
600   // 
601   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
602   // Landau with a Gaussian (see LandauGaus).  Note that 
603   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
604   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
605   //  
606   // References: 
607   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
608   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
609   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
610   // 
611   // Parameters:
612   //    x        Where to evaluate @f$ f_N@f$
613   //    delta    @f$ \Delta_1@f$ 
614   //    xi       @f$ \xi_1@f$
615   //    sigma    @f$ \sigma_1@f$ 
616   //    sigma_n  @f$ \sigma_n@f$ 
617   //    n        @f$ N@f$ in the sum above.
618   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
619   //                 @f$ i > 1@f$ 
620   // 
621   // Return:
622   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
623   //
624   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
625   for (Int_t i = 2; i <= n; i++) 
626     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
627   return result;
628 }
629 namespace { 
630   const Int_t kColors[] = { kRed+1, 
631                             kPink+3, 
632                             kMagenta+2, 
633                             kViolet+2, 
634                             kBlue+1, 
635                             kAzure+3, 
636                             kCyan+1, 
637                             kTeal+2, 
638                             kGreen+2, 
639                             kSpring+3, 
640                             kYellow+2, 
641                             kOrange+2 };
642 }
643
644 //____________________________________________________________________
645 TF1*
646 AliForwardUtil::MakeNLandauGaus(Double_t  c, 
647                                 Double_t  delta, Double_t xi, 
648                                 Double_t  sigma, Double_t sigmaN, Int_t n, 
649                                 const Double_t* a, 
650                                 Double_t  xmin, Double_t xmax)
651 {
652   // 
653   // Generate a TF1 object of @f$ f_N@f$ 
654   // 
655   // Parameters:
656   //    c         Constant                             
657   //    delta     @f$ \Delta@f$                        
658   //    xi            @f$ \xi_1@f$                     
659   //    sigma     @f$ \sigma_1@f$                      
660   //    sigma_n   @f$ \sigma_n@f$                      
661   //    n             @f$ N@f$ - how many particles to sum to
662   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
663   //                  @f$ i > 1@f$ 
664   //    xmin      Least value of range  
665   //    xmax      Largest value of range
666   // 
667   // Return:
668   //    Newly allocated TF1 object
669   //
670   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
671   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
672   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
673   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
674   landaun->SetLineWidth(2);
675   landaun->SetNpx(500);
676   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
677
678   // Set the initial parameters from the seed fit 
679   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
680   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
681   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
682   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
683   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
684   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);       
685
686   // Set the range and name of the scale parameters 
687   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
688     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
689     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
690   }
691   return landaun;
692 }
693 //____________________________________________________________________
694 TF1*
695 AliForwardUtil::MakeILandauGaus(Double_t  c, 
696                                 Double_t  delta, Double_t xi, 
697                                 Double_t  sigma, Double_t sigmaN, Int_t i, 
698                                 Double_t  xmin, Double_t xmax)
699 {
700   // 
701   // Generate a TF1 object of @f$ f_I@f$ 
702   // 
703   // Parameters:
704   //    c        Constant
705   //    delta    @f$ \Delta@f$ 
706   //    xi       @f$ \xi_1@f$          
707   //    sigma    @f$ \sigma_1@f$               
708   //    sigma_n  @f$ \sigma_n@f$               
709   //    i            @f$ i@f$ - the number of particles
710   //    xmin     Least value of range
711   //    xmax     Largest value of range
712   // 
713   // Return:
714   //    Newly allocated TF1 object
715   //
716   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
717   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
718   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
719   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
720   landaui->SetLineWidth(1);
721   landaui->SetNpx(500);
722   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
723
724   // Set the initial parameters from the seed fit 
725   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
726   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
727   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
728   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
729   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN); 
730   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);       
731
732   return landaui;
733 }
734
735 //====================================================================
736 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
737                                          Double_t maxRange, 
738                                          UShort_t minusBins) 
739   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
740     fFitResults(0), fFunctions(0)
741 {
742   // 
743   // Constructor 
744   // 
745   // Parameters:
746   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
747   //    maxRange   Maximum range to fit to 
748   //    minusBins  The number of bins below maximum to use 
749   //
750   fFitResults.SetOwner();
751   fFunctions.SetOwner();
752 }
753 //____________________________________________________________________
754 AliForwardUtil::ELossFitter::~ELossFitter()
755 {
756   // 
757   // Destructor
758   // 
759   //
760   fFitResults.Delete();
761   fFunctions.Delete();
762 }
763 //____________________________________________________________________
764 void
765 AliForwardUtil::ELossFitter::Clear()
766 {
767   // 
768   // Clear internal arrays 
769   // 
770   //
771   fFitResults.Clear();
772   fFunctions.Clear();
773 }
774 //____________________________________________________________________
775 TF1*
776 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
777 {
778   // 
779   // Fit a 1-particle signal to the passed energy loss distribution 
780   // 
781   // Note that this function clears the internal arrays first 
782   // 
783   // Parameters:
784   //    dist    Data to fit the function to 
785   //    sigman If larger than zero, the initial guess of the
786   //               detector induced noise. If zero or less, then this 
787   //               parameter is ignored in the fit (fixed at 0)
788   // 
789   // Return:
790   //    The function fitted to the data 
791   //
792
793   // Clear the cache 
794   Clear();
795   
796   // Find the fit range 
797   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
798   
799   // Get the bin with maximum 
800   Int_t    maxBin = dist->GetMaximumBin();
801   Double_t maxE   = dist->GetBinLowEdge(maxBin);
802   
803   // Get the low edge 
804   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
805   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
806   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
807   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
808
809   // Restore the range 
810   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
811   
812   // Define the function to fit 
813   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
814
815   // Set initial guesses, parameter names, and limits  
816   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
817   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
818   landau1->SetNpx(500);
819   landau1->SetParLimits(kDelta, minE, fMaxRange);
820   landau1->SetParLimits(kXi,    0.00, fMaxRange);
821   landau1->SetParLimits(kSigma, 0.01, 0.1);
822   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
823   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
824
825   // Do the fit, getting the result object 
826   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
827   landau1->SetRange(minE, fMaxRange);
828   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
829   fFunctions.AddAtAndExpand(landau1, 0);
830
831   return landau1;
832 }
833 //____________________________________________________________________
834 TF1*
835 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
836                                           Double_t sigman)
837 {
838   // 
839   // Fit a N-particle signal to the passed energy loss distribution 
840   //
841   // If there's no 1-particle fit present, it does that first 
842   // 
843   // Parameters:
844   //    dist   Data to fit the function to 
845   //    n      Number of particle signals to fit 
846   //    sigman If larger than zero, the initial guess of the
847   //               detector induced noise. If zero or less, then this 
848   //               parameter is ignored in the fit (fixed at 0)
849   // 
850   // Return:
851   //    The function fitted to the data 
852   //
853
854   // Get the seed fit result 
855   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
856   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
857   if (!r || !f) { 
858     f = Fit1Particle(dist, sigman);
859     r = static_cast<TFitResult*>(fFitResults.At(0));
860     if (!r || !f) { 
861       ::Warning("FitNLandau", "No first shot at landau fit");
862       return 0;
863     }
864   }
865
866   // Get some parameters from seed fit 
867   Double_t delta1  = r->Parameter(kDelta);
868   Double_t xi1     = r->Parameter(kXi);
869   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
870   Double_t minE    = f->GetXmin();
871
872   // Array of weights 
873   TArrayD a(n-1);
874   for (UShort_t i = 2; i <= n; i++) 
875     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
876   // Make the fit function 
877   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
878                                      r->Parameter(kDelta),
879                                      r->Parameter(kXi),
880                                      r->Parameter(kSigma),
881                                      r->Parameter(kSigmaN),
882                                      n,a.fArray,minE,maxEi);
883   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
884   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
885   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
886   // Check if we're using the noise sigma 
887   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
888   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
889
890   // Set the range and name of the scale parameters 
891   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
892     landaun->SetParLimits(kA+i-2, 0,1);
893   }
894
895   // Do the fit 
896   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
897   
898   landaun->SetRange(minE, fMaxRange);
899   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
900   fFunctions.AddAtAndExpand(landaun, n-1);
901   
902   return landaun;
903 }  
904
905 //====================================================================
906 AliForwardUtil::Histos::~Histos()
907 {
908   // 
909   // Destructor
910   //
911   if (fFMD1i) delete fFMD1i;
912   if (fFMD2i) delete fFMD2i;
913   if (fFMD2o) delete fFMD2o;
914   if (fFMD3i) delete fFMD3i;
915   if (fFMD3o) delete fFMD3o;
916 }
917
918 //____________________________________________________________________
919 TH2D*
920 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, 
921                              const TAxis& etaAxis) const
922 {
923   // 
924   // Make a histogram 
925   // 
926   // Parameters:
927   //    d        Detector
928   //    r        Ring 
929   //    etaAxis  Eta axis to use
930   // 
931   // Return:
932   //    Newly allocated histogram 
933   //
934   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
935   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r), 
936                         Form("FMD%d%c cache", d, r),
937                         etaAxis.GetNbins(), etaAxis.GetXmin(), 
938                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
939   hist->SetXTitle("#eta");
940   hist->SetYTitle("#phi [radians]");
941   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
942   hist->Sumw2();
943   hist->SetDirectory(0);
944
945   return hist;
946 }
947 //____________________________________________________________________
948 void
949 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
950 {
951   // 
952   // Initialize the object 
953   // 
954   // Parameters:
955   //    etaAxis Eta axis to use 
956   //
957   fFMD1i = Make(1, 'I', etaAxis);
958   fFMD2i = Make(2, 'I', etaAxis);
959   fFMD2o = Make(2, 'O', etaAxis);
960   fFMD3i = Make(3, 'I', etaAxis);
961   fFMD3o = Make(3, 'O', etaAxis);
962 }
963 //____________________________________________________________________
964 void
965 AliForwardUtil::Histos::Clear(Option_t* option)
966 {
967   // 
968   // Clear data 
969   // 
970   // Parameters:
971   //    option Not used 
972   //
973   if (fFMD1i) fFMD1i->Reset(option);
974   if (fFMD2i) fFMD2i->Reset(option);
975   if (fFMD2o) fFMD2o->Reset(option);
976   if (fFMD3i) fFMD3i->Reset(option);
977   if (fFMD3o) fFMD3o->Reset(option);
978 }
979
980 //____________________________________________________________________
981 TH2D*
982 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
983 {
984   // 
985   // Get the histogram for a particular detector,ring
986   // 
987   // Parameters:
988   //    d Detector 
989   //    r Ring 
990   // 
991   // Return:
992   //    Histogram for detector,ring or nul 
993   //
994   switch (d) { 
995   case 1: return fFMD1i;
996   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
997   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
998   }
999   return 0;
1000 }
1001 //====================================================================
1002 TList*
1003 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1004 {
1005   // 
1006   // Define the outout list in @a d
1007   // 
1008   // Parameters:
1009   //    d Where to put the output list
1010   // 
1011   // Return:
1012   //    Newly allocated TList object or null
1013   //
1014   if (!d) return 0;
1015   TList* list = new TList;
1016   list->SetOwner();
1017   list->SetName(fName.Data());
1018   d->Add(list);
1019   return list;
1020 }
1021 //____________________________________________________________________
1022 TList*
1023 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1024 {
1025   // 
1026   // Get our output list from the container @a d
1027   // 
1028   // Parameters:
1029   //    d where to get the output list from 
1030   // 
1031   // Return:
1032   //    The found TList or null
1033   //
1034   if (!d) return 0;
1035   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1036   return list;
1037 }
1038
1039 //____________________________________________________________________
1040 TH1*
1041 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1042 {
1043   // 
1044   // Find a specific histogram in the source list @a d
1045   // 
1046   // Parameters:
1047   //    d     (top)-container 
1048   //    name  Name of histogram
1049   // 
1050   // Return:
1051   //    Found histogram or null
1052   //
1053   return static_cast<TH1*>(d->FindObject(name));
1054 }
1055
1056 //====================================================================
1057 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl, 
1058                                        const char* format, ...)
1059   : fMsg("")
1060 {
1061   if (lvl < msgLvl) return; 
1062   va_list ap;
1063   va_start(ap, format);
1064   Format(fMsg, format, ap);
1065   va_end(ap);
1066   Output(+1, fMsg);
1067 }
1068 //____________________________________________________________________
1069 AliForwardUtil::DebugGuard::~DebugGuard()
1070 {
1071   if (fMsg.IsNull()) return;
1072   Output(-1, fMsg);
1073 }
1074 //____________________________________________________________________
1075 void
1076 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl, 
1077                                     const char* format, ...)
1078 {
1079   if (lvl < msgLvl) return; 
1080   TString msg;
1081   va_list ap;
1082   va_start(ap, format);
1083   Format(msg, format, ap);
1084   va_end(ap);
1085   Output(0, msg);
1086 }
1087
1088 //____________________________________________________________________
1089 void
1090 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1091 {
1092   static char buf[512];
1093   Int_t n = gROOT->GetDirLevel() + 2;
1094   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1095   vsnprintf(&(buf[n]), 511-n, format, ap);
1096   buf[511] = '\0';
1097   out = buf;  
1098 }
1099 //____________________________________________________________________
1100 void
1101 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1102 {
1103   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1104   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1105   if      (in > 0) gROOT->IncreaseDirLevel();
1106   else if (in < 0) gROOT->DecreaseDirLevel();
1107 }
1108
1109
1110
1111 //
1112 // EOF
1113 //