Improved eloss fitting - see NIM B1, 16
[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 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
19 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
20 namespace {
21   /** 
22    * The shift of the most probable value for the ROOT function TMath::Landau 
23    */
24   const Double_t  mpshift  = -0.22278298;
25   /** 
26    * Integration normalisation 
27    */
28   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
29
30   /** 
31    * Utility function to use in TF1 defintition 
32    */
33   Double_t landauGaus1(Double_t* xp, Double_t* pp) 
34   {
35     Double_t x        = xp[0];
36     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
37     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
38     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
39     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
40     Double_t sigma_n  = pp[AliForwardUtil::ELossFitter::kSigmaN];
41
42     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
43   }
44
45   /** 
46    * Utility function to use in TF1 defintition 
47    */
48   Double_t landauGausN(Double_t* xp, Double_t* pp) 
49   {
50     Double_t  x        = xp[0];
51     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
52     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
53     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
54     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
55     Double_t sigma_n   = pp[AliForwardUtil::ELossFitter::kSigmaN];
56     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
57     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
58
59     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
60                                                   n, a);
61   }
62
63
64 }
65 //____________________________________________________________________
66 Double_t 
67 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
68 {
69   return TMath::Landau(x, delta - xi * mpshift, xi);
70 }
71 //____________________________________________________________________
72 Double_t 
73 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
74                            Double_t sigma, Double_t sigma_n)
75 {
76   Double_t deltap = delta - xi * mpshift;
77   Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
78   Double_t sigma1 = sigma_n == 0 ? sigma : TMath::Sqrt(sigma2);
79   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
80   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
81   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
82   Double_t sum    = 0;
83   
84   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
85     Double_t x1 = xlow  + (i - .5) * step;
86     Double_t x2 = xhigh - (i - .5) * step;
87     
88     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
89     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
90   }
91   return step * sum * invSq2pi / sigma1;
92 }
93
94 //____________________________________________________________________
95 Double_t 
96 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
97                             Double_t sigma, Double_t sigma_n, Int_t n, 
98                             Double_t* a)
99 {
100   Double_t result = LandauGaus(x, delta, xi, sigma, sigma_n);
101   for (Int_t i = 2; i <= n; i++) { 
102     Double_t delta_i =  i * (delta + xi * TMath::Log(i));
103     Double_t xi_i    =  i * xi;
104     Double_t sigma_i =  TMath::Sqrt(Double_t(n))*sigma;
105     Double_t a_i     =  a[i-2];
106     result           += a_i * AliForwardUtil::LandauGaus(x, delta_i, xi_i, 
107                                                          sigma_i, sigma_n);
108   }
109   return result;
110 }
111
112 //====================================================================
113 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
114                                          Double_t maxRange, 
115                                          UShort_t minusBins) 
116   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
117     fFitResults(0), fFunctions(0)
118 {
119   fFitResults.SetOwner();
120   fFunctions.SetOwner();
121 }
122 //____________________________________________________________________
123 AliForwardUtil::ELossFitter::~ELossFitter()
124 {
125   fFitResults.Delete();
126   fFunctions.Delete();
127 }
128 //____________________________________________________________________
129 void
130 AliForwardUtil::ELossFitter::Clear()
131 {
132   fFitResults.Clear();
133   fFunctions.Clear();
134 }
135 //____________________________________________________________________
136 TF1*
137 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
138 {
139   // Clear the cache 
140   Clear();
141   
142   // Find the fit range 
143   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
144   
145   // Get the bin with maximum 
146   Int_t    maxBin = dist->GetMaximumBin();
147   Double_t maxE   = dist->GetBinLowEdge(maxBin);
148   
149   // Get the low edge 
150   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
151   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
152   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
153   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
154
155   // Restore the range 
156   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
157   
158   // Define the function to fit 
159   TF1*          landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
160
161   // Set initial guesses, parameter names, and limits  
162   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
163   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
164   landau1->SetNpx(500);
165   landau1->SetParLimits(kDelta, minE, fMaxRange);
166   landau1->SetParLimits(kXi,    0.00, fMaxRange);
167   landau1->SetParLimits(kSigma, 0.01, 0.1);
168   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
169   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
170
171   // Do the fit, getting the result object 
172   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
173   landau1->SetRange(minE, fMaxRange);
174   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
175   fFunctions.AddAtAndExpand(landau1, 0);
176
177   return landau1;
178 }
179 //____________________________________________________________________
180 TF1*
181 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
182                                           Double_t sigman)
183 {
184   // Get the seed fit result 
185   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
186   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
187   if (!r || !f) { 
188     f = Fit1Particle(dist, sigman);
189     r = static_cast<TFitResult*>(fFitResults.At(0));
190     if (!r || !f) { 
191       ::Warning("FitNLandau", "No first shot at landau fit");
192       return 0;
193     }
194   }
195
196   // Get some parameters from seed fit 
197   Double_t delta1  = r->Parameter(kDelta);
198   Double_t xi1     = r->Parameter(kXi);
199   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
200   Double_t minE    = f->GetXmin();
201
202   // Make the fit function 
203   TF1* landaun     = new TF1(Form("landau%d", n), &landauGausN,minE,maxEi,kN+n);
204   landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
205   landaun->SetLineColor(((n-2) % 10)+2); // start at red
206   landaun->SetLineWidth(1);
207   landaun->SetNpx(500);
208   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
209
210   // Set the initial parameters from the seed fit 
211   landaun->SetParameter(kC,      r->Parameter(kC));      // Constant
212   landaun->SetParameter(kDelta,  r->Parameter(kDelta));  // Delta 
213   landaun->SetParameter(kXi,     r->Parameter(kXi));     // xi
214   landaun->SetParameter(kSigma,  r->Parameter(kSigma));  // sigma
215   landaun->SetParameter(kSigmaN, r->Parameter(kSigmaN)); // sigma_n
216   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
217   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
218   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
219   // Check if we're using the noise sigma 
220   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
221   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
222   // Fix the number parameter 
223   landaun->FixParameter(kN, n);               // N
224
225   // Set the range and name of the scale parameters 
226   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
227     landaun->SetParameter(kA+i-2, n == 2 ? 0.05 : 0.000001);
228     landaun->SetParLimits(kA+i-2, 0,1);
229     landaun->SetParName(kA+i-2, Form("a_{%d}", i));
230   }
231
232   // Do the fit 
233   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
234   
235   landaun->SetRange(minE, fMaxRange);
236   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
237   fFunctions.AddAtAndExpand(landaun, n-1);
238   
239   return landaun;
240 }  
241
242 //====================================================================
243 AliForwardUtil::Histos::~Histos()
244 {
245   if (fFMD1i) delete fFMD1i;
246   if (fFMD2i) delete fFMD2i;
247   if (fFMD2o) delete fFMD2o;
248   if (fFMD3i) delete fFMD3i;
249   if (fFMD3o) delete fFMD3o;
250 }
251
252 //____________________________________________________________________
253 TH2D*
254 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, 
255                              const TAxis& etaAxis) const
256 {
257   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
258   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r), 
259                         Form("FMD%d%c cache", d, r),
260                         etaAxis.GetNbins(), etaAxis.GetXmin(), 
261                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
262   hist->SetXTitle("#eta");
263   hist->SetYTitle("#phi [radians]");
264   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
265   hist->Sumw2();
266   hist->SetDirectory(0);
267
268   return hist;
269 }
270 //____________________________________________________________________
271 void
272 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
273 {
274   fFMD1i = Make(1, 'I', etaAxis);
275   fFMD2i = Make(2, 'I', etaAxis);
276   fFMD2o = Make(2, 'O', etaAxis);
277   fFMD3i = Make(3, 'I', etaAxis);
278   fFMD3o = Make(3, 'O', etaAxis);
279 }
280 //____________________________________________________________________
281 void
282 AliForwardUtil::Histos::Clear(Option_t* option)
283 {
284   fFMD1i->Reset(option);
285   fFMD2i->Reset(option);
286   fFMD2o->Reset(option);
287   fFMD3i->Reset(option);
288   fFMD3o->Reset(option);
289 }
290
291 //____________________________________________________________________
292 TH2D*
293 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
294 {
295   switch (d) { 
296   case 1: return fFMD1i;
297   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
298   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
299   }
300   return 0;
301 }
302 //====================================================================
303 TList*
304 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
305 {
306   if (!d) return 0;
307   TList* list = new TList;
308   list->SetName(fName.Data());
309   d->Add(list);
310   return list;
311 }
312 //____________________________________________________________________
313 TList*
314 AliForwardUtil::RingHistos::GetOutputList(TList* d) const
315 {
316   if (!d) return 0;
317   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
318   return list;
319 }
320
321 //____________________________________________________________________
322 TH1*
323 AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
324 {
325   return static_cast<TH1*>(d->FindObject(name));
326 }
327
328 //
329 // EOF
330 //