More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.cxx
1 #include "AliForwardUtil.h"
2 #include <AliAnalysisManager.h>
3 #include "AliAODForwardMult.h"
4 #include <AliLog.h>
5 #include <AliInputEventHandler.h>
6 #include <AliESDEvent.h>
7 #include <AliPhysicsSelection.h>
8 #include <AliTriggerAnalysis.h>
9 #include <AliMultiplicity.h>
10 #include <TH2D.h>
11 #include <TH1I.h>
12 #include <TF1.h>
13 #include <TFitResult.h>
14 #include <TMath.h>
15 #include <TError.h>
16
17 //====================================================================
18 UShort_t
19 AliForwardUtil::ParseCollisionSystem(const char* sys)
20 {
21   TString s(sys);
22   s.ToLower();
23   if (s.Contains("p-p")   || s.Contains("pp"))   return AliForwardUtil::kPP; 
24   if (s.Contains("pb-pb") || s.Contains("pbpb")) return AliForwardUtil::kPbPb;
25   if (s.Contains("a-a")   || s.Contains("aa"))   return AliForwardUtil::kPbPb;
26   return AliForwardUtil::kUnknown;
27 }
28 //____________________________________________________________________
29 const char*
30 AliForwardUtil::CollisionSystemString(UShort_t sys)
31 {
32   switch (sys) { 
33   case AliForwardUtil::kPP:   return "pp";
34   case AliForwardUtil::kPbPb: return "PbPb";
35   }
36   return "unknown";
37 }
38 //____________________________________________________________________
39 UShort_t
40 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t v)
41 {
42   Float_t energy = v; 
43   if (sys != AliForwardUtil::kPP) energy = energy / 208 * 82;
44   if (TMath::Abs(energy - 900.)   < 10)  return 900;
45   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
46   if (TMath::Abs(energy - 2750.)  < 10)  return 2750;
47   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
48   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
49   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
50   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
51   return 0;
52 }
53 //____________________________________________________________________
54 const char* 
55 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
56 {
57   return Form("%04dGeV", cms);
58 }
59 //____________________________________________________________________
60 Short_t
61 AliForwardUtil::ParseMagneticField(Float_t v)
62 {
63   if (TMath::Abs(v - 5.) < 1 ) return +5;
64   if (TMath::Abs(v + 5.) < 1 ) return -5;
65   if (TMath::Abs(v) < 1)       return 0;
66   return 999;
67 }
68 //____________________________________________________________________
69 const char* 
70 AliForwardUtil::MagneticFieldString(Short_t f)
71 {
72   return Form("%01dkG", f);
73 }
74
75
76 //====================================================================
77 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
78 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
79 namespace {
80   /** 
81    * The shift of the most probable value for the ROOT function TMath::Landau 
82    */
83   const Double_t  mpshift  = -0.22278298;
84   /** 
85    * Integration normalisation 
86    */
87   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
88
89   /** 
90    * Utility function to use in TF1 defintition 
91    */
92   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
93   {
94     Double_t x        = xp[0];
95     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
96     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
97     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
98     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
99     Double_t sigma_n  = pp[AliForwardUtil::ELossFitter::kSigmaN];
100
101     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
102   }
103
104   /** 
105    * Utility function to use in TF1 defintition 
106    */
107   Double_t landauGausN(Double_t* xp, Double_t* pp) 
108   {
109     Double_t  x        = xp[0];
110     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
111     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
112     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
113     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
114     Double_t sigma_n   = pp[AliForwardUtil::ELossFitter::kSigmaN];
115     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
116     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
117
118     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
119                                                   n, a);
120   }
121   /** 
122    * Utility function to use in TF1 defintition 
123    */
124   Double_t landauGausI(Double_t* xp, Double_t* pp) 
125   {
126     Double_t x         = xp[0];
127     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
128     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
129     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
130     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
131     Double_t sigma_n   = pp[AliForwardUtil::ELossFitter::kSigmaN];
132     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
133
134     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
135   }
136
137
138 }
139 //____________________________________________________________________
140 Double_t 
141 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
142 {
143   return TMath::Landau(x, delta - xi * mpshift, xi);
144 }
145 //____________________________________________________________________
146 Double_t 
147 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
148                            Double_t sigma, Double_t sigma_n)
149 {
150   Double_t deltap = delta - xi * mpshift;
151   Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
152   Double_t sigma1 = sigma_n == 0 ? sigma : TMath::Sqrt(sigma2);
153   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
154   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
155   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
156   Double_t sum    = 0;
157   
158   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
159     Double_t x1 = xlow  + (i - .5) * step;
160     Double_t x2 = xhigh - (i - .5) * step;
161     
162     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
163     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
164   }
165   return step * sum * invSq2pi / sigma1;
166 }
167
168 //____________________________________________________________________
169 Double_t 
170 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
171                             Double_t sigma, Double_t sigma_n, Int_t i)
172 {
173   Double_t delta_i =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
174   Double_t xi_i    =  i * xi;
175   Double_t sigma_i =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
176   if (sigma_i < 1e-10) { 
177     // Fall back to landau 
178     return AliForwardUtil::Landau(x, delta_i, xi_i);
179   }
180   return AliForwardUtil::LandauGaus(x, delta_i, xi_i, sigma_i, sigma_n);
181 }
182
183 //____________________________________________________________________
184 Double_t 
185 AliForwardUtil::IdLandauGausdPar(Double_t x, 
186                                  UShort_t par,   Double_t dPar, 
187                                  Double_t delta, Double_t xi, 
188                                  Double_t sigma, Double_t sigma_n, 
189                                  Int_t    i)
190 {
191   if (dPar == 0) return 0;
192   Double_t dp      = dPar;
193   Double_t d2      = dPar / 2;
194   Double_t delta_i =  i * (delta + xi * TMath::Log(i));
195   Double_t xi_i    =  i * xi;
196   Double_t si      =  TMath::Sqrt(Double_t(i));
197   Double_t sigma_i =  si*sigma;
198   Double_t y1      = 0;
199   Double_t y2      = 0;
200   Double_t y3      = 0;
201   Double_t y4      = 0;
202   switch (par) {
203   case 0: 
204     y1 = ILandauGaus(x, delta_i+i*dp, xi_i, sigma_i, sigma_n, i);
205     y2 = ILandauGaus(x, delta_i+i*d2, xi_i, sigma_i, sigma_n, i);
206     y3 = ILandauGaus(x, delta_i-i*d2, xi_i, sigma_i, sigma_n, i);
207     y4 = ILandauGaus(x, delta_i-i*dp, xi_i, sigma_i, sigma_n, i);
208     break;
209   case 1: 
210     y1 = ILandauGaus(x, delta_i, xi_i+i*dp, sigma_i, sigma_n, i);
211     y2 = ILandauGaus(x, delta_i, xi_i+i*d2, sigma_i, sigma_n, i);
212     y3 = ILandauGaus(x, delta_i, xi_i-i*d2, sigma_i, sigma_n, i);
213     y4 = ILandauGaus(x, delta_i, xi_i-i*dp, sigma_i, sigma_n, i);
214     break;
215   case 2: 
216     y1 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*dp, sigma_n, i);
217     y2 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*d2, sigma_n, i);
218     y3 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*d2, sigma_n, i);
219     y4 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*dp, sigma_n, i);
220     break;
221   case 3: 
222     y1 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+dp, i);
223     y2 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+d2, i);
224     y3 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-d2, i);
225     y4 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-dp, i);
226     break;
227   default:
228     return 0;
229   } 
230   
231   Double_t d0  = y1 - y4;
232   Double_t d1  = 2 * (y2 - y3);
233   
234   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
235    
236   return g;
237 }
238
239 //____________________________________________________________________
240 Double_t 
241 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
242                             Double_t sigma, Double_t sigma_n, Int_t n, 
243                             Double_t* a)
244 {
245   Double_t result = ILandauGaus(x, delta, xi, sigma, sigma_n, 1);
246   for (Int_t i = 2; i <= n; i++) 
247     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
248   return result;
249 }
250 namespace { 
251   const Int_t kColors[] = { kRed+1, 
252                             kPink+3, 
253                             kMagenta+2, 
254                             kViolet+2, 
255                             kBlue+1, 
256                             kAzure+3, 
257                             kCyan+1, 
258                             kTeal+2, 
259                             kGreen+2, 
260                             kSpring+3, 
261                             kYellow+2, 
262                             kOrange+2 };
263 }
264
265 //____________________________________________________________________
266 TF1*
267 AliForwardUtil::MakeNLandauGaus(Double_t  c, 
268                                 Double_t  delta, Double_t xi, 
269                                 Double_t  sigma, Double_t sigma_n, Int_t n, 
270                                 Double_t* a, 
271                                 Double_t  xmin, Double_t xmax)
272 {
273   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
274   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
275   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
276   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
277   landaun->SetLineWidth(2);
278   landaun->SetNpx(500);
279   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
280
281   // Set the initial parameters from the seed fit 
282   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
283   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
284   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
285   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
286   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n); 
287   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);       
288
289   // Set the range and name of the scale parameters 
290   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
291     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
292     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
293   }
294   return landaun;
295 }
296 //____________________________________________________________________
297 TF1*
298 AliForwardUtil::MakeILandauGaus(Double_t  c, 
299                                 Double_t  delta, Double_t xi, 
300                                 Double_t  sigma, Double_t sigma_n, Int_t i, 
301                                 Double_t  xmin, Double_t xmax)
302 {
303   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
304   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
305   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
306   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
307   landaui->SetLineWidth(1);
308   landaui->SetNpx(500);
309   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
310
311   // Set the initial parameters from the seed fit 
312   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);       
313   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);   
314   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);      
315   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);   
316   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n); 
317   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);       
318
319   return landaui;
320 }
321
322 //====================================================================
323 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
324                                          Double_t maxRange, 
325                                          UShort_t minusBins) 
326   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
327     fFitResults(0), fFunctions(0)
328 {
329   fFitResults.SetOwner();
330   fFunctions.SetOwner();
331 }
332 //____________________________________________________________________
333 AliForwardUtil::ELossFitter::~ELossFitter()
334 {
335   fFitResults.Delete();
336   fFunctions.Delete();
337 }
338 //____________________________________________________________________
339 void
340 AliForwardUtil::ELossFitter::Clear()
341 {
342   fFitResults.Clear();
343   fFunctions.Clear();
344 }
345 //____________________________________________________________________
346 TF1*
347 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
348 {
349   // Clear the cache 
350   Clear();
351   
352   // Find the fit range 
353   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
354   
355   // Get the bin with maximum 
356   Int_t    maxBin = dist->GetMaximumBin();
357   Double_t maxE   = dist->GetBinLowEdge(maxBin);
358   
359   // Get the low edge 
360   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
361   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
362   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
363   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
364
365   // Restore the range 
366   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
367   
368   // Define the function to fit 
369   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
370
371   // Set initial guesses, parameter names, and limits  
372   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
373   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
374   landau1->SetNpx(500);
375   landau1->SetParLimits(kDelta, minE, fMaxRange);
376   landau1->SetParLimits(kXi,    0.00, fMaxRange);
377   landau1->SetParLimits(kSigma, 0.01, 0.1);
378   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
379   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
380
381   // Do the fit, getting the result object 
382   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
383   landau1->SetRange(minE, fMaxRange);
384   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
385   fFunctions.AddAtAndExpand(landau1, 0);
386
387   return landau1;
388 }
389 //____________________________________________________________________
390 TF1*
391 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
392                                           Double_t sigman)
393 {
394   // Get the seed fit result 
395   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
396   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
397   if (!r || !f) { 
398     f = Fit1Particle(dist, sigman);
399     r = static_cast<TFitResult*>(fFitResults.At(0));
400     if (!r || !f) { 
401       ::Warning("FitNLandau", "No first shot at landau fit");
402       return 0;
403     }
404   }
405
406   // Get some parameters from seed fit 
407   Double_t delta1  = r->Parameter(kDelta);
408   Double_t xi1     = r->Parameter(kXi);
409   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
410   Double_t minE    = f->GetXmin();
411
412   // Array of weights 
413   TArrayD a(n-1);
414   for (UShort_t i = 2; i <= n; i++) 
415     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
416   // Make the fit function 
417   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
418                                      r->Parameter(kDelta),
419                                      r->Parameter(kXi),
420                                      r->Parameter(kSigma),
421                                      r->Parameter(kSigmaN),
422                                      n,a.fArray,minE,maxEi);
423   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
424   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
425   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
426   // Check if we're using the noise sigma 
427   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
428   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
429
430   // Set the range and name of the scale parameters 
431   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
432     landaun->SetParLimits(kA+i-2, 0,1);
433   }
434
435   // Do the fit 
436   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
437   
438   landaun->SetRange(minE, fMaxRange);
439   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
440   fFunctions.AddAtAndExpand(landaun, n-1);
441   
442   return landaun;
443 }  
444
445 //====================================================================
446 AliForwardUtil::Histos::~Histos()
447 {
448   if (fFMD1i) delete fFMD1i;
449   if (fFMD2i) delete fFMD2i;
450   if (fFMD2o) delete fFMD2o;
451   if (fFMD3i) delete fFMD3i;
452   if (fFMD3o) delete fFMD3o;
453 }
454
455 //____________________________________________________________________
456 TH2D*
457 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, 
458                              const TAxis& etaAxis) const
459 {
460   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
461   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r), 
462                         Form("FMD%d%c cache", d, r),
463                         etaAxis.GetNbins(), etaAxis.GetXmin(), 
464                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
465   hist->SetXTitle("#eta");
466   hist->SetYTitle("#phi [radians]");
467   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
468   hist->Sumw2();
469   hist->SetDirectory(0);
470
471   return hist;
472 }
473 //____________________________________________________________________
474 void
475 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
476 {
477   fFMD1i = Make(1, 'I', etaAxis);
478   fFMD2i = Make(2, 'I', etaAxis);
479   fFMD2o = Make(2, 'O', etaAxis);
480   fFMD3i = Make(3, 'I', etaAxis);
481   fFMD3o = Make(3, 'O', etaAxis);
482 }
483 //____________________________________________________________________
484 void
485 AliForwardUtil::Histos::Clear(Option_t* option)
486 {
487   fFMD1i->Reset(option);
488   fFMD2i->Reset(option);
489   fFMD2o->Reset(option);
490   fFMD3i->Reset(option);
491   fFMD3o->Reset(option);
492 }
493
494 //____________________________________________________________________
495 TH2D*
496 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
497 {
498   switch (d) { 
499   case 1: return fFMD1i;
500   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
501   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
502   }
503   return 0;
504 }
505 //====================================================================
506 TList*
507 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
508 {
509   if (!d) return 0;
510   TList* list = new TList;
511   list->SetName(fName.Data());
512   d->Add(list);
513   return list;
514 }
515 //____________________________________________________________________
516 TList*
517 AliForwardUtil::RingHistos::GetOutputList(TList* d) const
518 {
519   if (!d) return 0;
520   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
521   return list;
522 }
523
524 //____________________________________________________________________
525 TH1*
526 AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
527 {
528   return static_cast<TH1*>(d->FindObject(name));
529 }
530
531 //
532 // EOF
533 //