More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.cxx
CommitLineData
7e4038b5 1#include "AliForwardUtil.h"
9d99b0dd 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>
7e4038b5 10#include <TH2D.h>
9d99b0dd 11#include <TH1I.h>
7f759bb7 12#include <TF1.h>
13#include <TFitResult.h>
7e4038b5 14#include <TMath.h>
7f759bb7 15#include <TError.h>
16
0bd4b00f 17//====================================================================
18UShort_t
19AliForwardUtil::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//____________________________________________________________________
29const char*
30AliForwardUtil::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//____________________________________________________________________
39UShort_t
40AliForwardUtil::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//____________________________________________________________________
54const char*
55AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
56{
57 return Form("%04dGeV", cms);
58}
59//____________________________________________________________________
60Short_t
61AliForwardUtil::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//____________________________________________________________________
69const char*
70AliForwardUtil::MagneticFieldString(Short_t f)
71{
72 return Form("%01dkG", f);
73}
74
75
7f759bb7 76//====================================================================
77Int_t AliForwardUtil::fgConvolutionSteps = 100;
78Double_t AliForwardUtil::fgConvolutionNSigma = 5;
79namespace {
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];
c389303e 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];
7f759bb7 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];
c389303e 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]);
7f759bb7 117
118 return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
119 n, a);
120 }
0bd4b00f 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 }
7f759bb7 136
137
138}
139//____________________________________________________________________
140Double_t
141AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
142{
143 return TMath::Landau(x, delta - xi * mpshift, xi);
144}
145//____________________________________________________________________
146Double_t
147AliForwardUtil::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;
c389303e 152 Double_t sigma1 = sigma_n == 0 ? sigma : TMath::Sqrt(sigma2);
7f759bb7 153 Double_t xlow = x - fgConvolutionNSigma * sigma1;
c389303e 154 Double_t xhigh = x + fgConvolutionNSigma * sigma1;
7f759bb7 155 Double_t step = (xhigh - xlow) / fgConvolutionSteps;
156 Double_t sum = 0;
157
158 for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
c389303e 159 Double_t x1 = xlow + (i - .5) * step;
160 Double_t x2 = xhigh - (i - .5) * step;
7f759bb7 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
0bd4b00f 168//____________________________________________________________________
169Double_t
170AliForwardUtil::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//____________________________________________________________________
184Double_t
185AliForwardUtil::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
7f759bb7 239//____________________________________________________________________
240Double_t
241AliForwardUtil::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{
0bd4b00f 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);
7f759bb7 248 return result;
249}
0bd4b00f 250namespace {
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//____________________________________________________________________
266TF1*
267AliForwardUtil::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//____________________________________________________________________
297TF1*
298AliForwardUtil::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}
7f759bb7 321
322//====================================================================
323AliForwardUtil::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//____________________________________________________________________
333AliForwardUtil::ELossFitter::~ELossFitter()
334{
335 fFitResults.Delete();
336 fFunctions.Delete();
337}
338//____________________________________________________________________
339void
340AliForwardUtil::ELossFitter::Clear()
341{
342 fFitResults.Clear();
343 fFunctions.Clear();
344}
345//____________________________________________________________________
346TF1*
347AliForwardUtil::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
7f759bb7 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
0bd4b00f 369 TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
7f759bb7 370
371 // Set initial guesses, parameter names, and limits
c389303e 372 landau1->SetParameters(1,0.5,0.07,0.1,sigman);
7f759bb7 373 landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
c389303e 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);
7f759bb7 380
381 // Do the fit, getting the result object
382 TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
c389303e 383 landau1->SetRange(minE, fMaxRange);
7f759bb7 384 fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
385 fFunctions.AddAtAndExpand(landau1, 0);
386
387 return landau1;
388}
389//____________________________________________________________________
390TF1*
391AliForwardUtil::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
c389303e 407 Double_t delta1 = r->Parameter(kDelta);
408 Double_t xi1 = r->Parameter(kXi);
7f759bb7 409 Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
410 Double_t minE = f->GetXmin();
411
0bd4b00f 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);
7f759bb7 416 // Make the fit function
0bd4b00f 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);
c389303e 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);
7f759bb7 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
c389303e 432 landaun->SetParLimits(kA+i-2, 0,1);
7f759bb7 433 }
434
435 // Do the fit
436 TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
437
c389303e 438 landaun->SetRange(minE, fMaxRange);
7f759bb7 439 fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
440 fFunctions.AddAtAndExpand(landaun, n-1);
441
442 return landaun;
443}
7e4038b5 444
445//====================================================================
446AliForwardUtil::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//____________________________________________________________________
456TH2D*
457AliForwardUtil::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//____________________________________________________________________
474void
475AliForwardUtil::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//____________________________________________________________________
484void
485AliForwardUtil::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//____________________________________________________________________
495TH2D*
496AliForwardUtil::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}
9d99b0dd 505//====================================================================
506TList*
507AliForwardUtil::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//____________________________________________________________________
516TList*
517AliForwardUtil::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//____________________________________________________________________
525TH1*
526AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
527{
528 return static_cast<TH1*>(d->FindObject(name));
529}
530
7e4038b5 531//
532// EOF
533//