]>
Commit | Line | Data |
---|---|---|
f8715167 | 1 | // |
2 | // Class to do the energy correction of FMD ESD data | |
3 | // | |
4 | #include "AliFMDEnergyFitter.h" | |
0bd4b00f | 5 | #include "AliForwardCorrectionManager.h" |
f8715167 | 6 | #include <AliESDFMD.h> |
7 | #include <TAxis.h> | |
8 | #include <TList.h> | |
9 | #include <TH1.h> | |
10 | #include <TF1.h> | |
11 | #include <TMath.h> | |
f8715167 | 12 | #include <AliLog.h> |
13 | #include <TClonesArray.h> | |
14 | #include <TFitResult.h> | |
15 | #include <THStack.h> | |
0bd4b00f | 16 | #include <TROOT.h> |
17 | #include <iostream> | |
18 | #include <iomanip> | |
f8715167 | 19 | |
20 | ClassImp(AliFMDEnergyFitter) | |
21 | #if 0 | |
22 | ; // This is for Emacs | |
23 | #endif | |
0bd4b00f | 24 | namespace { |
25 | const char* fgkEDistFormat = "%s_etabin%03d"; | |
26 | } | |
f8715167 | 27 | |
f8715167 | 28 | |
29 | //____________________________________________________________________ | |
30 | AliFMDEnergyFitter::AliFMDEnergyFitter() | |
31 | : TNamed(), | |
32 | fRingHistos(), | |
9d05ffeb | 33 | fLowCut(0.4), |
34 | fNParticles(5), | |
35 | fMinEntries(1000), | |
c389303e | 36 | fFitRangeBinWidth(4), |
f8715167 | 37 | fDoFits(false), |
0bd4b00f | 38 | fDoMakeObject(false), |
f8715167 | 39 | fEtaAxis(), |
5e4d905e | 40 | fCentralityAxis(), |
4b9857f3 | 41 | fMaxE(10), |
42 | fNEbins(300), | |
43 | fUseIncreasingBins(true), | |
c389303e | 44 | fMaxRelParError(.25), |
45 | fMaxChi2PerNDF(20), | |
0bd4b00f | 46 | fMinWeight(1e-7), |
f8715167 | 47 | fDebug(0) |
7984e5f7 | 48 | { |
49 | // | |
50 | // Default Constructor - do not use | |
51 | // | |
241cca4d | 52 | DGUARD(fDebug, 1, "Default CTOR of AliFMDEnergyFitter"); |
7984e5f7 | 53 | } |
f8715167 | 54 | |
55 | //____________________________________________________________________ | |
56 | AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title) | |
57 | : TNamed("fmdEnergyFitter", title), | |
58 | fRingHistos(), | |
9d05ffeb | 59 | fLowCut(0.4), |
60 | fNParticles(5), | |
61 | fMinEntries(1000), | |
c389303e | 62 | fFitRangeBinWidth(4), |
f8715167 | 63 | fDoFits(false), |
0bd4b00f | 64 | fDoMakeObject(false), |
4b9857f3 | 65 | fEtaAxis(0,0,0), |
5e4d905e | 66 | fCentralityAxis(0,0,0), |
4b9857f3 | 67 | fMaxE(10), |
68 | fNEbins(300), | |
69 | fUseIncreasingBins(true), | |
c389303e | 70 | fMaxRelParError(.25), |
0bd4b00f | 71 | fMaxChi2PerNDF(20), |
72 | fMinWeight(1e-7), | |
f8715167 | 73 | fDebug(3) |
74 | { | |
7984e5f7 | 75 | // |
76 | // Constructor | |
77 | // | |
78 | // Parameters: | |
79 | // title Title of object - not significant | |
80 | // | |
241cca4d | 81 | DGUARD(fDebug, 1, "Named CTOR of AliFMDEnergyFitter: %s", title); |
f8715167 | 82 | fEtaAxis.SetName("etaAxis"); |
83 | fEtaAxis.SetTitle("#eta"); | |
5e4d905e | 84 | fCentralityAxis.SetName("centralityAxis"); |
85 | fCentralityAxis.SetTitle("Centrality [%]"); | |
f8715167 | 86 | fRingHistos.Add(new RingHistos(1, 'I')); |
87 | fRingHistos.Add(new RingHistos(2, 'I')); | |
88 | fRingHistos.Add(new RingHistos(2, 'O')); | |
89 | fRingHistos.Add(new RingHistos(3, 'I')); | |
90 | fRingHistos.Add(new RingHistos(3, 'O')); | |
91 | } | |
92 | ||
93 | //____________________________________________________________________ | |
94 | AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o) | |
95 | : TNamed(o), | |
96 | fRingHistos(), | |
97 | fLowCut(o.fLowCut), | |
c389303e | 98 | fNParticles(o.fNParticles), |
f8715167 | 99 | fMinEntries(o.fMinEntries), |
c389303e | 100 | fFitRangeBinWidth(o.fFitRangeBinWidth), |
f8715167 | 101 | fDoFits(o.fDoFits), |
0bd4b00f | 102 | fDoMakeObject(o.fDoMakeObject), |
f8715167 | 103 | fEtaAxis(o.fEtaAxis), |
5e4d905e | 104 | fCentralityAxis(o.fCentralityAxis), |
4b9857f3 | 105 | fMaxE(o.fMaxE), |
106 | fNEbins(o.fNEbins), | |
107 | fUseIncreasingBins(o.fUseIncreasingBins), | |
c389303e | 108 | fMaxRelParError(o.fMaxRelParError), |
109 | fMaxChi2PerNDF(o.fMaxChi2PerNDF), | |
0bd4b00f | 110 | fMinWeight(o.fMinWeight), |
f8715167 | 111 | fDebug(o.fDebug) |
112 | { | |
7984e5f7 | 113 | // |
114 | // Copy constructor | |
115 | // | |
116 | // Parameters: | |
117 | // o Object to copy from | |
118 | // | |
241cca4d | 119 | DGUARD(fDebug, 1, "Copy CTOR of AliFMDEnergyFitter"); |
f8715167 | 120 | TIter next(&o.fRingHistos); |
121 | TObject* obj = 0; | |
122 | while ((obj = next())) fRingHistos.Add(obj); | |
123 | } | |
124 | ||
125 | //____________________________________________________________________ | |
126 | AliFMDEnergyFitter::~AliFMDEnergyFitter() | |
127 | { | |
7984e5f7 | 128 | // |
129 | // Destructor | |
130 | // | |
f8715167 | 131 | fRingHistos.Delete(); |
132 | } | |
133 | ||
134 | //____________________________________________________________________ | |
135 | AliFMDEnergyFitter& | |
136 | AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o) | |
137 | { | |
7984e5f7 | 138 | // |
139 | // Assignment operator | |
140 | // | |
141 | // Parameters: | |
142 | // o Object to assign from | |
143 | // | |
144 | // Return: | |
145 | // Reference to this | |
146 | // | |
6ab100ec | 147 | DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter"); |
d015ecfe | 148 | if (&o == this) return *this; |
f8715167 | 149 | TNamed::operator=(o); |
150 | ||
151 | fLowCut = o.fLowCut; | |
c389303e | 152 | fNParticles = o.fNParticles; |
f8715167 | 153 | fMinEntries = o.fMinEntries; |
c389303e | 154 | fFitRangeBinWidth= o.fFitRangeBinWidth; |
f8715167 | 155 | fDoFits = o.fDoFits; |
0bd4b00f | 156 | fDoMakeObject = o.fDoMakeObject; |
5e4d905e | 157 | fEtaAxis.Set(o.fEtaAxis.GetNbins(), |
158 | o.fEtaAxis.GetXmin(), | |
159 | o.fEtaAxis.GetXmax()); | |
160 | if (o.fCentralityAxis.GetXbins()) { | |
161 | const TArrayD* bins = o.fCentralityAxis.GetXbins(); | |
162 | fCentralityAxis.Set(bins->GetSize()-1,bins->GetArray()); | |
163 | } | |
164 | else | |
165 | fCentralityAxis.Set(o.fCentralityAxis.GetNbins(), | |
166 | o.fCentralityAxis.GetXmin(), | |
167 | o.fCentralityAxis.GetXmax()); | |
f8715167 | 168 | fDebug = o.fDebug; |
0bd4b00f | 169 | fMaxRelParError= o.fMaxRelParError; |
170 | fMaxChi2PerNDF = o.fMaxChi2PerNDF; | |
171 | fMinWeight = o.fMinWeight; | |
f8715167 | 172 | |
173 | fRingHistos.Delete(); | |
174 | TIter next(&o.fRingHistos); | |
175 | TObject* obj = 0; | |
176 | while ((obj = next())) fRingHistos.Add(obj); | |
177 | ||
178 | return *this; | |
179 | } | |
180 | ||
181 | //____________________________________________________________________ | |
182 | AliFMDEnergyFitter::RingHistos* | |
183 | AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const | |
184 | { | |
7984e5f7 | 185 | // |
186 | // Get the ring histogram container | |
187 | // | |
188 | // Parameters: | |
189 | // d Detector | |
190 | // r Ring | |
191 | // | |
192 | // Return: | |
193 | // Ring histogram container | |
194 | // | |
f8715167 | 195 | Int_t idx = -1; |
196 | switch (d) { | |
197 | case 1: idx = 0; break; | |
198 | case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break; | |
199 | case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break; | |
200 | } | |
201 | if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0; | |
202 | ||
203 | return static_cast<RingHistos*>(fRingHistos.At(idx)); | |
204 | } | |
205 | ||
206 | //____________________________________________________________________ | |
207 | void | |
208 | AliFMDEnergyFitter::Init(const TAxis& eAxis) | |
209 | { | |
7984e5f7 | 210 | // |
211 | // Initialise the task | |
212 | // | |
213 | // Parameters: | |
214 | // etaAxis The eta axis to use. Note, that if the eta axis | |
215 | // has already been set (using SetEtaAxis), then this parameter will be | |
216 | // ignored | |
217 | // | |
6ab100ec | 218 | DGUARD(fDebug, 1, "Initialize of AliFMDEnergyFitter"); |
4b9857f3 | 219 | if (fEtaAxis.GetNbins() == 0 || |
7c1a1f1d | 220 | TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7) |
4b9857f3 | 221 | SetEtaAxis(eAxis); |
5e4d905e | 222 | if (fCentralityAxis.GetNbins() == 0) { |
223 | UShort_t n = 12; | |
224 | Double_t bins[] = { 0., 5., 10., 15., 20., 30., | |
225 | 40., 50., 60., 70., 80., 100. }; | |
226 | SetCentralityAxis(n, bins); | |
227 | } | |
f8715167 | 228 | TIter next(&fRingHistos); |
229 | RingHistos* o = 0; | |
230 | while ((o = static_cast<RingHistos*>(next()))) | |
5e4d905e | 231 | o->Init(fEtaAxis, fCentralityAxis, fMaxE, fNEbins, fUseIncreasingBins); |
f8715167 | 232 | } |
4b9857f3 | 233 | //____________________________________________________________________ |
234 | void | |
235 | AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis) | |
236 | { | |
7984e5f7 | 237 | // |
238 | // Set the eta axis to use. This will force the code to use this | |
239 | // eta axis definition - irrespective of whatever axis is passed to | |
240 | // the Init member function. Therefore, this member function can be | |
241 | // used to force another eta axis than one found in the correction | |
242 | // objects. | |
243 | // | |
244 | // Parameters: | |
245 | // etaAxis Eta axis to use | |
246 | // | |
4b9857f3 | 247 | SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax()); |
248 | } | |
249 | //____________________________________________________________________ | |
250 | void | |
251 | AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax) | |
252 | { | |
7984e5f7 | 253 | // |
254 | // Set the eta axis to use. This will force the code to use this | |
255 | // eta axis definition - irrespective of whatever axis is passed to | |
256 | // the Init member function. Therefore, this member function can be | |
257 | // used to force another eta axis than one found in the correction | |
258 | // objects. | |
259 | // | |
260 | // Parameters: | |
261 | // nBins Number of bins | |
262 | // etaMin Minimum of the eta axis | |
263 | // etaMax Maximum of the eta axis | |
264 | // | |
4b9857f3 | 265 | fEtaAxis.Set(nBins,etaMin,etaMax); |
266 | } | |
5e4d905e | 267 | //____________________________________________________________________ |
268 | void | |
269 | AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins) | |
270 | { | |
271 | fCentralityAxis.Set(n-1, bins); | |
272 | } | |
273 | ||
4b9857f3 | 274 | |
f8715167 | 275 | //____________________________________________________________________ |
276 | Bool_t | |
5e4d905e | 277 | AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, |
278 | Double_t cent, | |
f8715167 | 279 | Bool_t empty) |
280 | { | |
7984e5f7 | 281 | // |
282 | // Fitter the input AliESDFMD object | |
283 | // | |
284 | // Parameters: | |
285 | // input Input | |
5e4d905e | 286 | // cent Centrality |
7984e5f7 | 287 | // empty Whether the event is 'empty' |
288 | // | |
289 | // Return: | |
290 | // True on success, false otherwise | |
291 | // | |
6ab100ec | 292 | DGUARD(fDebug, 1, "Accumulate statistics in AliFMDEnergyFitter"); |
5e4d905e | 293 | Int_t icent = fCentralityAxis.FindBin(cent); |
294 | if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0; | |
295 | ||
241cca4d | 296 | UShort_t nFills = 0; |
f8715167 | 297 | for(UShort_t d = 1; d <= 3; d++) { |
298 | Int_t nRings = (d == 1 ? 1 : 2); | |
299 | for (UShort_t q = 0; q < nRings; q++) { | |
300 | Char_t r = (q == 0 ? 'I' : 'O'); | |
301 | UShort_t nsec = (q == 0 ? 20 : 40); | |
302 | UShort_t nstr = (q == 0 ? 512 : 256); | |
303 | ||
304 | RingHistos* histos = GetRingHistos(d, r); | |
305 | ||
306 | for(UShort_t s = 0; s < nsec; s++) { | |
307 | for(UShort_t t = 0; t < nstr; t++) { | |
308 | Float_t mult = input.Multiplicity(d,r,s,t); | |
309 | ||
310 | // Keep dead-channel information. | |
311 | if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0) | |
312 | continue; | |
313 | ||
314 | // Get the pseudo-rapidity | |
315 | Double_t eta1 = input.Eta(d,r,s,t); | |
316 | Int_t ieta = fEtaAxis.FindBin(eta1); | |
317 | if (ieta < 1 || ieta > fEtaAxis.GetNbins()) continue; | |
318 | ||
5e4d905e | 319 | histos->Fill(empty, ieta-1, icent, mult); |
241cca4d | 320 | nFills++; |
f8715167 | 321 | } // for strip |
322 | } // for sector | |
323 | } // for ring | |
324 | } // for detector | |
325 | ||
241cca4d | 326 | DMSG(fDebug, 1, "Found a total of %d signals for c=%f, and %sempty event", |
327 | nFills, cent, (empty ? "" : "non-")); | |
f8715167 | 328 | return kTRUE; |
329 | } | |
330 | ||
331 | //____________________________________________________________________ | |
332 | void | |
7c1a1f1d | 333 | AliFMDEnergyFitter::Fit(const TList* dir) |
f8715167 | 334 | { |
7984e5f7 | 335 | // |
336 | // Scale the histograms to the total number of events | |
337 | // | |
338 | // Parameters: | |
339 | // dir Where the histograms are | |
340 | // | |
6ab100ec | 341 | DGUARD(fDebug, 1, "Fit distributions in AliFMDEnergyFitter"); |
f8715167 | 342 | if (!fDoFits) return; |
343 | ||
344 | TList* d = static_cast<TList*>(dir->FindObject(GetName())); | |
345 | if (!d) return; | |
346 | ||
347 | // +1 for chi^2 | |
348 | // +3 for 1 landau | |
349 | // +1 for N | |
c389303e | 350 | // +fNParticles-1 for weights |
351 | Int_t nStack = kN+fNParticles; | |
f8715167 | 352 | THStack* stack[nStack]; |
c389303e | 353 | stack[0] = new THStack("chi2", "#chi^{2}/#nu"); |
354 | stack[kC +1] = new THStack("c", "Constant"); | |
355 | stack[kDelta +1] = new THStack("delta", "#Delta_{p}"); | |
356 | stack[kXi +1] = new THStack("xi", "#xi"); | |
357 | stack[kSigma +1] = new THStack("sigma", "#sigma"); | |
358 | stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}"); | |
359 | stack[kN +1] = new THStack("n", "# Particles"); | |
360 | for (Int_t i = 2; i <= fNParticles; i++) | |
361 | stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i)); | |
f8715167 | 362 | for (Int_t i = 0; i < nStack; i++) |
363 | d->Add(stack[i]); | |
364 | ||
365 | TIter next(&fRingHistos); | |
366 | RingHistos* o = 0; | |
367 | while ((o = static_cast<RingHistos*>(next()))) { | |
c389303e | 368 | TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles, |
369 | fMinEntries, fFitRangeBinWidth, | |
370 | fMaxRelParError, fMaxChi2PerNDF); | |
f8715167 | 371 | if (!l) continue; |
372 | for (Int_t i = 0; i < l->GetEntriesFast(); i++) { | |
373 | stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); | |
374 | } | |
375 | } | |
0bd4b00f | 376 | |
377 | if (!fDoMakeObject) return; | |
378 | ||
379 | MakeCorrectionsObject(d); | |
380 | } | |
381 | ||
382 | //____________________________________________________________________ | |
383 | void | |
384 | AliFMDEnergyFitter::MakeCorrectionsObject(TList* d) | |
385 | { | |
7984e5f7 | 386 | // |
387 | // Generate the corrections object | |
388 | // | |
389 | // Parameters: | |
390 | // dir List to analyse | |
391 | // | |
6ab100ec | 392 | DGUARD(fDebug, 1, "Make the correction objec in AliFMDEnergyFitter"); |
0bd4b00f | 393 | AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance(); |
394 | ||
395 | AliFMDCorrELossFit* obj = new AliFMDCorrELossFit; | |
396 | obj->SetEtaAxis(fEtaAxis); | |
397 | obj->SetLowCut(fLowCut); | |
398 | ||
399 | TIter next(&fRingHistos); | |
400 | RingHistos* o = 0; | |
401 | while ((o = static_cast<RingHistos*>(next()))) { | |
402 | o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError, | |
403 | fMaxChi2PerNDF, fMinWeight); | |
404 | } | |
405 | ||
406 | TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits)); | |
407 | TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits)); | |
408 | AliInfo(Form("Object %s created in output - should be extracted and copied " | |
409 | "to %s", oName.Data(), fileName.Data())); | |
f7cfc454 | 410 | d->Add(new TNamed("filename", fileName.Data())); |
0bd4b00f | 411 | d->Add(obj, oName.Data()); |
f8715167 | 412 | } |
413 | ||
414 | //____________________________________________________________________ | |
415 | void | |
416 | AliFMDEnergyFitter::DefineOutput(TList* dir) | |
417 | { | |
7984e5f7 | 418 | // |
419 | // Define the output histograms. These are put in a sub list of the | |
420 | // passed list. The histograms are merged before the parent task calls | |
421 | // AliAnalysisTaskSE::Terminate | |
422 | // | |
423 | // Parameters: | |
424 | // dir Directory to add to | |
425 | // | |
6ab100ec | 426 | DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter"); |
f8715167 | 427 | TList* d = new TList; |
428 | d->SetName(GetName()); | |
429 | dir->Add(d); | |
430 | d->Add(&fEtaAxis); | |
431 | TIter next(&fRingHistos); | |
432 | RingHistos* o = 0; | |
433 | while ((o = static_cast<RingHistos*>(next()))) { | |
434 | o->Output(d); | |
435 | } | |
436 | } | |
437 | //____________________________________________________________________ | |
438 | void | |
439 | AliFMDEnergyFitter::SetDebug(Int_t dbg) | |
440 | { | |
7984e5f7 | 441 | // |
442 | // Set the debug level. The higher the value the more output | |
443 | // | |
444 | // Parameters: | |
445 | // dbg Debug level | |
446 | // | |
f8715167 | 447 | fDebug = dbg; |
448 | TIter next(&fRingHistos); | |
449 | RingHistos* o = 0; | |
450 | while ((o = static_cast<RingHistos*>(next()))) | |
451 | o->fDebug = dbg; | |
452 | } | |
0bd4b00f | 453 | |
454 | //____________________________________________________________________ | |
455 | void | |
456 | AliFMDEnergyFitter::Print(Option_t*) const | |
457 | { | |
7984e5f7 | 458 | // |
459 | // Print information | |
460 | // | |
461 | // Parameters: | |
462 | // option Not used | |
463 | // | |
0bd4b00f | 464 | char ind[gROOT->GetDirLevel()+1]; |
465 | for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; | |
466 | ind[gROOT->GetDirLevel()] = '\0'; | |
1f480471 | 467 | std::cout << ind << ClassName() << ": " << GetName() << '\n' |
0bd4b00f | 468 | << ind << " Low cut: " << fLowCut << " E/E_mip\n" |
469 | << ind << " Max(particles): " << fNParticles << '\n' | |
470 | << ind << " Min(entries): " << fMinEntries << '\n' | |
471 | << ind << " Fit range: " | |
472 | << fFitRangeBinWidth << " bins\n" | |
473 | << ind << " Make fits: " | |
474 | << (fDoFits ? "yes\n" : "no\n") | |
475 | << ind << " Make object: " | |
476 | << (fDoMakeObject ? "yes\n" : "no\n") | |
477 | << ind << " Max E/E_mip: " << fMaxE << '\n' | |
478 | << ind << " N bins: " << fNEbins << '\n' | |
479 | << ind << " Increasing bins: " | |
480 | << (fUseIncreasingBins ?"yes\n" : "no\n") | |
481 | << ind << " max(delta p/p): " << fMaxRelParError << '\n' | |
482 | << ind << " max(chi^2/nu): " << fMaxChi2PerNDF << '\n' | |
483 | << ind << " min(a_i): " << fMinWeight << std::endl; | |
484 | } | |
f8715167 | 485 | |
486 | //==================================================================== | |
487 | AliFMDEnergyFitter::RingHistos::RingHistos() | |
488 | : AliForwardUtil::RingHistos(), | |
489 | fEDist(0), | |
490 | fEmpty(0), | |
491 | fEtaEDists(), | |
492 | fList(0), | |
0bd4b00f | 493 | fFits("AliFMDCorrELossFit::ELossFit"), |
f8715167 | 494 | fDebug(0) |
7984e5f7 | 495 | { |
496 | // | |
497 | // Default CTOR | |
498 | // | |
6ab100ec | 499 | DGUARD(fDebug, 0, "Default CTOR AliFMDEnergyFitter::RingHistos"); |
7984e5f7 | 500 | } |
f8715167 | 501 | |
502 | //____________________________________________________________________ | |
503 | AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r) | |
504 | : AliForwardUtil::RingHistos(d,r), | |
505 | fEDist(0), | |
506 | fEmpty(0), | |
507 | fEtaEDists(), | |
508 | fList(0), | |
0bd4b00f | 509 | fFits("AliFMDCorrELossFit::ELossFit"), |
f8715167 | 510 | fDebug(0) |
511 | { | |
7984e5f7 | 512 | // |
513 | // Constructor | |
514 | // | |
515 | // Parameters: | |
516 | // d detector | |
517 | // r ring | |
518 | // | |
f8715167 | 519 | fEtaEDists.SetName("EDists"); |
6ab100ec | 520 | DGUARD(fDebug, 0, "Named CTOR AliFMDEnergyFitter::RingHistos: FMD%d%c", |
521 | d, r); | |
f8715167 | 522 | } |
523 | //____________________________________________________________________ | |
524 | AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o) | |
525 | : AliForwardUtil::RingHistos(o), | |
526 | fEDist(o.fEDist), | |
527 | fEmpty(o.fEmpty), | |
528 | fEtaEDists(), | |
529 | fList(0), | |
0bd4b00f | 530 | fFits("AliFMDCorrELossFit::ELossFit"), |
f8715167 | 531 | fDebug(0) |
532 | { | |
7984e5f7 | 533 | // |
534 | // Copy constructor | |
535 | // | |
536 | // Parameters: | |
537 | // o Object to copy from | |
538 | // | |
6ab100ec | 539 | DGUARD(fDebug, 0, "Copy CTOR AliFMDEnergyFitter::RingHistos"); |
0bd4b00f | 540 | fFits.Clear(); |
f8715167 | 541 | TIter next(&o.fEtaEDists); |
542 | TObject* obj = 0; | |
543 | while ((obj = next())) fEtaEDists.Add(obj->Clone()); | |
544 | if (o.fList) { | |
545 | fList = new TList; | |
546 | fList->SetName(fName); | |
547 | TIter next2(o.fList); | |
548 | while ((obj = next2())) fList->Add(obj->Clone()); | |
549 | } | |
550 | } | |
551 | ||
552 | //____________________________________________________________________ | |
553 | AliFMDEnergyFitter::RingHistos& | |
554 | AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o) | |
555 | { | |
7984e5f7 | 556 | // |
557 | // Assignment operator | |
558 | // | |
559 | // Parameters: | |
560 | // o Object to assign from | |
561 | // | |
562 | // Return: | |
563 | // Reference to this | |
564 | // | |
6ab100ec | 565 | DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter::RingHistos"); |
d015ecfe | 566 | if (&o == this) return *this; |
f8715167 | 567 | AliForwardUtil::RingHistos::operator=(o); |
568 | ||
569 | if (fEDist) delete fEDist; | |
570 | if (fEmpty) delete fEmpty; | |
571 | fEtaEDists.Delete(); | |
572 | if (fList) fList->Delete(); | |
0bd4b00f | 573 | fFits.Clear(); |
f8715167 | 574 | |
575 | fEDist = static_cast<TH1D*>(o.fEDist->Clone()); | |
576 | fEmpty = static_cast<TH1D*>(o.fEmpty->Clone()); | |
577 | ||
578 | TIter next(&o.fEtaEDists); | |
579 | TObject* obj = 0; | |
580 | while ((obj = next())) fEtaEDists.Add(obj->Clone()); | |
581 | ||
582 | if (o.fList) { | |
583 | fList = new TList; | |
584 | fList->SetName(fName); | |
585 | TIter next2(o.fList); | |
586 | while ((obj = next2())) fList->Add(obj->Clone()); | |
587 | } | |
588 | ||
589 | return *this; | |
590 | } | |
591 | //____________________________________________________________________ | |
592 | AliFMDEnergyFitter::RingHistos::~RingHistos() | |
593 | { | |
7984e5f7 | 594 | // |
595 | // Destructor | |
596 | // | |
6ab100ec | 597 | DGUARD(fDebug, 3, "DTOR of AliFMDEnergyFitter::RingHistos"); |
f8715167 | 598 | if (fEDist) delete fEDist; |
599 | fEtaEDists.Delete(); | |
600 | } | |
601 | ||
602 | //____________________________________________________________________ | |
603 | void | |
5e4d905e | 604 | AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, |
605 | Int_t /* icent */, Double_t mult) | |
f8715167 | 606 | { |
7984e5f7 | 607 | // |
608 | // Fill histogram | |
609 | // | |
610 | // Parameters: | |
611 | // empty True if event is empty | |
5e4d905e | 612 | // ieta Eta bin (0-based) |
613 | // icent Centrality bin (1-based) | |
7984e5f7 | 614 | // mult Signal |
615 | // | |
241cca4d | 616 | DGUARD(fDebug, 10, "Filling in AliFMDEnergyFitter::RingHistos"); |
f8715167 | 617 | if (empty) { |
618 | fEmpty->Fill(mult); | |
619 | return; | |
620 | } | |
621 | fEDist->Fill(mult); | |
622 | ||
5e4d905e | 623 | // Here, we should first look up in a centrality array, and then in |
624 | // that array look up the eta bin - or perhaps we should do it the | |
625 | // other way around? | |
f8715167 | 626 | if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return; |
627 | ||
628 | TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta)); | |
629 | if (!h) return; | |
630 | ||
631 | h->Fill(mult); | |
632 | } | |
633 | ||
634 | //__________________________________________________________________ | |
635 | TArrayD | |
66b34429 | 636 | AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, |
637 | Double_t max) const | |
f8715167 | 638 | { |
7984e5f7 | 639 | // |
640 | // Make an axis with increasing bins | |
641 | // | |
642 | // Parameters: | |
643 | // n Number of bins | |
644 | // min Minimum | |
645 | // max Maximum | |
646 | // | |
647 | // Return: | |
648 | // An axis with quadratically increasing bin size | |
649 | // | |
650 | ||
f8715167 | 651 | TArrayD bins(n+1); |
652 | Double_t dx = (max-min) / n; | |
653 | bins[0] = min; | |
654 | Int_t i = 1; | |
655 | for (i = 1; i < n+1; i++) { | |
656 | Double_t dI = float(i)/n; | |
66b34429 | 657 | Double_t next = bins[i-1] + dx + dI * dI * dx; |
f8715167 | 658 | bins[i] = next; |
659 | if (next > max) break; | |
660 | } | |
661 | bins.Set(i+1); | |
662 | return bins; | |
663 | } | |
664 | ||
665 | //____________________________________________________________________ | |
666 | void | |
c389303e | 667 | AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, |
668 | Double_t emin, | |
669 | Double_t emax, | |
670 | Double_t deMax, | |
671 | Int_t nDeBins, | |
672 | Bool_t incr) | |
f8715167 | 673 | { |
7984e5f7 | 674 | // |
675 | // Make E/E_mip histogram | |
676 | // | |
677 | // Parameters: | |
678 | // ieta Eta bin | |
679 | // eMin Least signal | |
680 | // eMax Largest signal | |
681 | // deMax Maximum energy loss | |
682 | // nDeBins Number energy loss bins | |
683 | // incr Whether to make bins of increasing size | |
684 | // | |
f8715167 | 685 | TH1D* h = 0; |
0bd4b00f | 686 | TString name = Form(fgkEDistFormat, fName.Data(), ieta); |
f8715167 | 687 | TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f " |
688 | "(#eta bin %d)", fName.Data(), emin, emax, ieta); | |
689 | if (!incr) | |
690 | h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax); | |
691 | else { | |
692 | TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax); | |
693 | h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray); | |
694 | } | |
695 | ||
696 | h->SetDirectory(0); | |
697 | h->SetXTitle("#DeltaE/#DeltaE_{mip}"); | |
698 | h->SetFillColor(Color()); | |
699 | h->SetMarkerColor(Color()); | |
700 | h->SetLineColor(Color()); | |
701 | h->SetFillStyle(3001); | |
702 | h->Sumw2(); | |
703 | ||
704 | fEtaEDists.AddAt(h, ieta-1); | |
705 | } | |
706 | //____________________________________________________________________ | |
c389303e | 707 | TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, |
708 | const char* title, | |
709 | const TAxis& eta) const | |
f8715167 | 710 | { |
7984e5f7 | 711 | // |
712 | // Make a parameter histogram | |
713 | // | |
714 | // Parameters: | |
715 | // name Name of histogram. | |
716 | // title Title of histogram. | |
717 | // eta Eta axis | |
718 | // | |
719 | // Return: | |
720 | // | |
721 | // | |
f8715167 | 722 | TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), |
723 | Form("%s for %s", title, fName.Data()), | |
724 | eta.GetNbins(), eta.GetXmin(), eta.GetXmax()); | |
725 | h->SetXTitle("#eta"); | |
726 | h->SetYTitle(title); | |
727 | h->SetDirectory(0); | |
728 | h->SetFillColor(Color()); | |
729 | h->SetMarkerColor(Color()); | |
730 | h->SetLineColor(Color()); | |
731 | h->SetFillStyle(3001); | |
732 | h->Sumw2(); | |
733 | ||
734 | return h; | |
735 | } | |
736 | //____________________________________________________________________ | |
737 | TH1D* | |
738 | AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, | |
739 | const char* title, | |
740 | const TAxis& eta, | |
741 | Int_t low, | |
742 | Int_t high, | |
743 | Double_t val, | |
744 | Double_t err) const | |
745 | { | |
7984e5f7 | 746 | // |
747 | // Make a histogram that contains the results of the fit over the full ring | |
748 | // | |
749 | // Parameters: | |
750 | // name Name | |
751 | // title Title | |
752 | // eta Eta axis | |
753 | // low Least bin | |
754 | // high Largest bin | |
755 | // val Value of parameter | |
756 | // err Error on parameter | |
757 | // | |
758 | // Return: | |
759 | // The newly allocated histogram | |
760 | // | |
f8715167 | 761 | Double_t xlow = eta.GetBinLowEdge(low); |
762 | Double_t xhigh = eta.GetBinUpEdge(high); | |
763 | TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), | |
764 | Form("%s for %s", title, fName.Data()), | |
765 | 1, xlow, xhigh); | |
766 | h->SetBinContent(1, val); | |
767 | h->SetBinError(1, err); | |
768 | h->SetXTitle("#eta"); | |
769 | h->SetYTitle(title); | |
770 | h->SetDirectory(0); | |
771 | h->SetFillColor(0); | |
772 | h->SetMarkerColor(Color()); | |
773 | h->SetLineColor(Color()); | |
774 | h->SetFillStyle(0); | |
775 | h->SetLineStyle(2); | |
776 | h->SetLineWidth(2); | |
777 | ||
778 | return h; | |
779 | } | |
780 | ||
781 | //____________________________________________________________________ | |
782 | void | |
783 | AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, | |
5e4d905e | 784 | const TAxis& /* cAxis */, |
f8715167 | 785 | Double_t maxDE, Int_t nDEbins, |
786 | Bool_t useIncrBin) | |
787 | { | |
7984e5f7 | 788 | // |
789 | // Initialise object | |
790 | // | |
791 | // Parameters: | |
792 | // eAxis Eta axis | |
793 | // maxDE Max energy loss to consider | |
794 | // nDEbins Number of bins | |
795 | // useIncrBin Whether to use an increasing bin size | |
796 | // | |
6ab100ec | 797 | DGUARD(fDebug, 2, "Initialize in AliFMDEnergyFitter::RingHistos"); |
f8715167 | 798 | fEDist = new TH1D(Form("%s_edist", fName.Data()), |
799 | Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), | |
800 | 200, 0, 6); | |
801 | fEmpty = new TH1D(Form("%s_empty", fName.Data()), | |
802 | Form("#DeltaE/#DeltaE_{mip} for %s (empty events)", | |
803 | fName.Data()), 200, 0, 6); | |
804 | fList->Add(fEDist); | |
805 | fList->Add(fEmpty); | |
5e4d905e | 806 | // Here, we should make an array of centrality bins ... |
807 | // In that case, the structure will be | |
808 | // | |
809 | // -+- Ring1 -+- Centrality1 -+- Eta1 | |
810 | // | | +- Eta2 | |
811 | // ... ... ... | |
812 | // | | +- EtaM | |
813 | // | +- Centrality2 -+- Eta1 | |
814 | // ... ... ... | |
815 | // | | +- EtaM | |
816 | // ... ... | |
817 | // | +- CentralityN -+- Eta1 | |
818 | // ... ... ... | |
819 | // | +- EtaM | |
820 | // -+- Ring2 -+- Centrality1 -+- Eta1 | |
821 | // | | +- Eta2 | |
822 | // ... ... ... | |
823 | // | | +- EtaM | |
824 | // | +- Centrality2 -+- Eta1 | |
825 | // ... ... ... | |
826 | // | | +- EtaM | |
827 | // ... ... | |
828 | // | +- CentralityN -+- Eta1 | |
829 | // ... ... ... | |
830 | // | +- EtaM | |
831 | // ... ... ... | |
832 | // -+- RingO -+- Centrality1 -+- Eta1 | |
833 | // | +- Eta2 | |
834 | // ... ... | |
835 | // | +- EtaM | |
836 | // +- Centrality2 -+- Eta1 | |
837 | // ... ... | |
838 | // | +- EtaM | |
839 | // ... | |
840 | // +- CentralityN -+- Eta1 | |
841 | // ... ... | |
842 | // +- EtaM | |
843 | // | |
f8715167 | 844 | // fEtaEDists.Expand(eAxis.GetNbins()); |
845 | for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { | |
846 | Double_t min = eAxis.GetBinLowEdge(i); | |
847 | Double_t max = eAxis.GetBinUpEdge(i); | |
5e4d905e | 848 | // Or may we should do it here? In that case we'd have the |
849 | // following structure: | |
850 | // Ring1 -+- Eta1 -+- Centrality1 | |
851 | // | +- Centrality2 | |
852 | // ... ... | |
853 | // | +- CentralityN | |
854 | // +- Eta2 -+- Centrality1 | |
855 | // | +- Centrality2 | |
856 | // ... ... | |
857 | // | +- CentralityN | |
858 | // ... | |
859 | // +- EtaM -+- Centrality1 | |
860 | // +- Centrality2 | |
861 | // ... | |
862 | // +- CentralityN | |
f8715167 | 863 | Make(i, min, max, maxDE, nDEbins, useIncrBin); |
864 | } | |
865 | fList->Add(&fEtaEDists); | |
e2213ed5 | 866 | // fEtaEDists.ls(); |
f8715167 | 867 | } |
868 | //____________________________________________________________________ | |
869 | TObjArray* | |
0bd4b00f | 870 | AliFMDEnergyFitter::RingHistos::Fit(TList* dir, |
871 | const TAxis& eta, | |
872 | Double_t lowCut, | |
873 | UShort_t nParticles, | |
874 | UShort_t minEntries, | |
875 | UShort_t minusBins, | |
876 | Double_t relErrorCut, | |
877 | Double_t chi2nuCut) const | |
f8715167 | 878 | { |
7984e5f7 | 879 | // |
880 | // Fit each histogram to up to @a nParticles particle responses. | |
881 | // | |
882 | // Parameters: | |
883 | // dir Output list | |
884 | // eta Eta axis | |
885 | // lowCut Lower cut | |
886 | // nParticles Max number of convolved landaus to fit | |
887 | // minEntries Minimum number of entries | |
888 | // minusBins Number of bins from peak to subtract to | |
889 | // get the fit range | |
890 | // relErrorCut Cut applied to relative error of parameter. | |
891 | // Note, for multi-particle weights, the cut | |
892 | // is loosend by a factor of 2 | |
893 | // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - | |
894 | // the reduced @f$\chi^2@f$ | |
895 | // | |
6ab100ec | 896 | DGUARD(fDebug, 2, "Fit in AliFMDEnergyFitter::RingHistos"); |
7984e5f7 | 897 | |
c389303e | 898 | // Get our ring list from the output container |
f8715167 | 899 | TList* l = GetOutputList(dir); |
900 | if (!l) return 0; | |
901 | ||
c389303e | 902 | // Get the energy distributions from the output container |
f8715167 | 903 | TList* dists = static_cast<TList*>(l->FindObject("EDists")); |
904 | if (!dists) { | |
905 | AliWarning(Form("Didn't find %s_EtaEDists in %s", | |
906 | fName.Data(), l->GetName())); | |
907 | l->ls(); | |
908 | return 0; | |
909 | } | |
910 | ||
c389303e | 911 | // Container of the fit results histograms |
912 | TObjArray* pars = new TObjArray(3+nParticles+1); | |
f8715167 | 913 | pars->SetName("FitResults"); |
914 | l->Add(pars); | |
915 | ||
c389303e | 916 | // Result objects stored in separate list on the output |
917 | TH1* hChi2 = 0; | |
918 | TH1* hC = 0; | |
919 | TH1* hDelta = 0; | |
920 | TH1* hXi = 0; | |
921 | TH1* hSigma = 0; | |
922 | TH1* hSigmaN = 0; | |
923 | TH1* hN = 0; | |
924 | TH1* hA[nParticles-1]; | |
925 | pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta)); | |
926 | pars->Add(hC = MakePar("c", "Constant", eta)); | |
927 | pars->Add(hDelta = MakePar("delta", "#Delta_{p}", eta)); | |
928 | pars->Add(hXi = MakePar("xi", "#xi", eta)); | |
929 | pars->Add(hSigma = MakePar("sigma", "#sigma", eta)); | |
930 | pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}", eta)); | |
931 | pars->Add(hN = MakePar("n", "N", eta)); | |
932 | for (UShort_t i = 1; i < nParticles; i++) | |
f8715167 | 933 | pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta)); |
934 | ||
935 | ||
936 | Int_t nDists = dists->GetEntries(); | |
937 | Int_t low = nDists; | |
938 | Int_t high = 0; | |
9b2f2e39 | 939 | Int_t nEmpty = 0; |
940 | Int_t nLow = 0; | |
941 | Int_t nFitted= 0; | |
f8715167 | 942 | for (Int_t i = 0; i < nDists; i++) { |
943 | TH1D* dist = static_cast<TH1D*>(dists->At(i)); | |
c389303e | 944 | // Ignore empty histograms altoghether |
9b2f2e39 | 945 | if (!dist || dist->GetEntries() <= 0) { |
946 | nEmpty++; | |
947 | continue; | |
948 | } | |
f8715167 | 949 | |
c389303e | 950 | // Scale to the bin-width |
951 | dist->Scale(1., "width"); | |
952 | ||
953 | // Normalize peak to 1 | |
954 | Double_t max = dist->GetMaximum(); | |
0bd4b00f | 955 | if (max <= 0) continue; |
c389303e | 956 | dist->Scale(1/max); |
957 | ||
958 | // Check that we have enough entries | |
9b2f2e39 | 959 | Int_t nEntries = Int_t(dist->GetEntries()); |
960 | if (nEntries <= minEntries) { | |
0bd4b00f | 961 | AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)", |
9b2f2e39 | 962 | i, dist->GetName(), nEntries, minEntries)); |
963 | nLow++; | |
0bd4b00f | 964 | continue; |
965 | } | |
f8715167 | 966 | |
c389303e | 967 | // Now fit |
968 | TF1* res = FitHist(dist,lowCut,nParticles,minusBins, | |
969 | relErrorCut,chi2nuCut); | |
f8715167 | 970 | if (!res) continue; |
9b2f2e39 | 971 | nFitted++; |
c389303e | 972 | // dist->GetListOfFunctions()->Add(res); |
0bd4b00f | 973 | |
c389303e | 974 | // Store eta limits |
f8715167 | 975 | low = TMath::Min(low,i+1); |
976 | high = TMath::Max(high,i+1); | |
977 | ||
c389303e | 978 | // Get the reduced chi square |
f8715167 | 979 | Double_t chi2 = res->GetChisquare(); |
980 | Int_t ndf = res->GetNDF(); | |
c389303e | 981 | |
982 | // Store results of best fit in output histograms | |
983 | res->SetLineWidth(4); | |
984 | hChi2 ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0); | |
985 | hC ->SetBinContent(i+1, res->GetParameter(kC)); | |
986 | hDelta ->SetBinContent(i+1, res->GetParameter(kDelta)); | |
987 | hXi ->SetBinContent(i+1, res->GetParameter(kXi)); | |
988 | hSigma ->SetBinContent(i+1, res->GetParameter(kSigma)); | |
989 | hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN)); | |
990 | hN ->SetBinContent(i+1, res->GetParameter(kN)); | |
991 | ||
992 | hC ->SetBinError(i+1, res->GetParError(kC)); | |
993 | hDelta ->SetBinError(i+1, res->GetParError(kDelta)); | |
994 | hXi ->SetBinError(i+1, res->GetParError(kXi)); | |
995 | hSigma ->SetBinError(i+1, res->GetParError(kSigma)); | |
996 | hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN)); | |
997 | hN ->SetBinError(i+1, res->GetParError(kN)); | |
998 | ||
999 | for (UShort_t j = 0; j < nParticles-1; j++) { | |
1000 | hA[j]->SetBinContent(i+1, res->GetParameter(kA+j)); | |
1001 | hA[j]->SetBinError(i+1, res->GetParError(kA+j)); | |
f8715167 | 1002 | } |
1003 | } | |
9b2f2e39 | 1004 | printf("%s: Out of %d histograms, %d where empty, %d had too little data," |
1005 | "leaving %d to be fitted, of which %d succeeded\n", | |
1006 | GetName(), nDists, nEmpty, nLow, nDists-nEmpty-nLow, nFitted); | |
1007 | ||
1008 | TH1* status = new TH1I("status", "Status of Fits", 5, 0, 5); | |
1009 | status->GetXaxis()->SetBinLabel(1, "Total"); | |
1010 | status->GetXaxis()->SetBinLabel(2, "Empty"); | |
1011 | status->GetXaxis()->SetBinLabel(3, Form("<%d", minEntries)); | |
1012 | status->GetXaxis()->SetBinLabel(4, "Candidates"); | |
1013 | status->GetXaxis()->SetBinLabel(5, "Fitted"); | |
1014 | status->SetXTitle("Status"); | |
1015 | status->SetYTitle("# of #Delta distributions"); | |
1016 | status->SetBinContent(1, nDists); | |
1017 | status->SetBinContent(2, nEmpty); | |
1018 | status->SetBinContent(3, nLow); | |
1019 | status->SetBinContent(4, nDists-nLow-nEmpty); | |
1020 | status->SetBinContent(5, nFitted); | |
1021 | status->SetFillColor(Color()); | |
1022 | status->SetFillStyle(3001); | |
1023 | status->SetLineColor(Color()); | |
1024 | status->SetDirectory(0); | |
1025 | status->SetStats(0); | |
1026 | pars->Add(status); | |
f8715167 | 1027 | |
c389303e | 1028 | // Fit the full-ring histogram |
f8715167 | 1029 | TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data())); |
4b9857f3 | 1030 | if (total && total->GetEntries() >= minEntries) { |
0bd4b00f | 1031 | |
1032 | // Scale to the bin-width | |
1033 | total->Scale(1., "width"); | |
1034 | ||
1035 | // Normalize peak to 1 | |
1036 | Double_t max = total->GetMaximum(); | |
1037 | if (max > 0) total->Scale(1/max); | |
1038 | ||
c389303e | 1039 | TF1* res = FitHist(total,lowCut,nParticles,minusBins, |
1040 | relErrorCut,chi2nuCut); | |
f8715167 | 1041 | if (res) { |
c389303e | 1042 | // Make histograms for the result of this fit |
f8715167 | 1043 | Double_t chi2 = res->GetChisquare(); |
1044 | Int_t ndf = res->GetNDF(); | |
c389303e | 1045 | pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high, |
f8715167 | 1046 | ndf > 0 ? chi2/ndf : 0, 0)); |
c389303e | 1047 | pars->Add(MakeTotal("t_c", "Constant", eta, low, high, |
1048 | res->GetParameter(kC), | |
1049 | res->GetParError(kC))); | |
1050 | pars->Add(MakeTotal("t_delta", "#Delta_{p}", eta, low, high, | |
1051 | res->GetParameter(kDelta), | |
1052 | res->GetParError(kDelta))); | |
1053 | pars->Add(MakeTotal("t_xi", "#xi", eta, low, high, | |
1054 | res->GetParameter(kXi), | |
1055 | res->GetParError(kXi))); | |
1056 | pars->Add(MakeTotal("t_sigma", "#sigma", eta, low, high, | |
1057 | res->GetParameter(kSigma), | |
1058 | res->GetParError(kSigma))); | |
1059 | pars->Add(MakeTotal("t_sigman", "#sigma_{n}", eta, low, high, | |
1060 | res->GetParameter(kSigmaN), | |
1061 | res->GetParError(kSigmaN))); | |
1062 | pars->Add(MakeTotal("t_n", "N", eta, low, high, | |
1063 | res->GetParameter(kN),0)); | |
1064 | for (UShort_t j = 0; j < nParticles-1; j++) | |
1065 | pars->Add(MakeTotal(Form("t_a%d",j+2), | |
1066 | Form("a_{%d}",j+2), eta, low, high, | |
1067 | res->GetParameter(kA+j), | |
1068 | res->GetParError(kA+j))); | |
f8715167 | 1069 | } |
1070 | } | |
1071 | ||
c389303e | 1072 | // Clean up list of histogram. Histograms with no entries or |
1073 | // no functions are deleted. We have to do this using the TObjLink | |
1074 | // objects stored in the list since ROOT cannot guaranty the validity | |
1075 | // of iterators when removing from a list - tsck. Should just implement | |
1076 | // TIter::Remove(). | |
f8715167 | 1077 | TObjLink* lnk = dists->FirstLink(); |
1078 | while (lnk) { | |
1079 | TH1* h = static_cast<TH1*>(lnk->GetObject()); | |
0bd4b00f | 1080 | bool remove = false; |
1081 | if (h->GetEntries() <= 0) { | |
1082 | // AliWarning(Form("No entries in %s - removing", h->GetName())); | |
1083 | remove = true; | |
1084 | } | |
1085 | else if (h->GetListOfFunctions()->GetEntries() <= 0) { | |
1086 | // AliWarning(Form("No fuctions associated with %s - removing", | |
1087 | // h->GetName())); | |
1088 | remove = true; | |
1089 | } | |
1090 | if (remove) { | |
f8715167 | 1091 | TObjLink* keep = lnk->Next(); |
1092 | dists->Remove(lnk); | |
1093 | lnk = keep; | |
1094 | continue; | |
1095 | } | |
1096 | lnk = lnk->Next(); | |
1097 | } | |
1098 | ||
1099 | return pars; | |
1100 | } | |
1101 | ||
1102 | //____________________________________________________________________ | |
1103 | TF1* | |
4b9857f3 | 1104 | AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, |
1105 | Double_t lowCut, | |
c389303e | 1106 | UShort_t nParticles, |
1107 | UShort_t minusBins, | |
1108 | Double_t relErrorCut, | |
1109 | Double_t chi2nuCut) const | |
f8715167 | 1110 | { |
7984e5f7 | 1111 | // |
1112 | // Fit a signal histogram. First, the bin @f$ b_{min}@f$ with | |
1113 | // maximum bin content in the range @f$ [E_{min},\infty]@f$ is | |
1114 | // found. Then the fit range is set to the bin range | |
1115 | // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1 | |
1116 | // particle signal is fitted to that. The parameters of that fit | |
1117 | // is then used as seeds for a fit of the @f$ N@f$ particle response | |
1118 | // to the data in the range | |
1119 | // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$ | |
1120 | // | |
1121 | // Parameters: | |
1122 | // dist Histogram to fit | |
1123 | // lowCut Lower cut @f$ E_{min}@f$ on signal | |
1124 | // nParticles Max number @f$ N@f$ of convolved landaus to fit | |
1125 | // minusBins Number of bins @f$ \Delta b@f$ from peak to | |
1126 | // subtract to get the fit range | |
1127 | // relErrorCut Cut applied to relative error of parameter. | |
1128 | // Note, for multi-particle weights, the cut | |
1129 | // is loosend by a factor of 2 | |
1130 | // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - | |
1131 | // the reduced @f$\chi^2@f$ | |
1132 | // | |
1133 | // Return: | |
1134 | // The best fit function | |
1135 | // | |
6ab100ec | 1136 | DGUARD(fDebug, 3, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s", |
1137 | dist->GetName()); | |
f8715167 | 1138 | Double_t maxRange = 10; |
1139 | ||
c389303e | 1140 | // Create a fitter object |
4b9857f3 | 1141 | AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); |
1142 | f.Clear(); | |
f8715167 | 1143 | |
c389303e | 1144 | |
4b9857f3 | 1145 | // If we are only asked to fit a single particle, return this fit, |
1146 | // no matter what. | |
c389303e | 1147 | if (nParticles == 1) { |
4b9857f3 | 1148 | TF1* r = f.Fit1Particle(dist, 0); |
1149 | if (!r) return 0; | |
e2213ed5 | 1150 | TF1* ret = new TF1(*r); |
1151 | dist->GetListOfFunctions()->Add(ret); | |
1152 | return ret; | |
f8715167 | 1153 | } |
f8715167 | 1154 | |
4b9857f3 | 1155 | // Fit from 2 upto n particles |
c389303e | 1156 | for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0); |
f8715167 | 1157 | |
f8715167 | 1158 | |
4b9857f3 | 1159 | // Now, we need to select the best fit |
fb3430ac | 1160 | Int_t nFits = f.GetFitResults().GetEntriesFast(); |
4b9857f3 | 1161 | TF1* good[nFits]; |
1162 | for (Int_t i = nFits-1; i >= 0; i--) { | |
c389303e | 1163 | good[i] = 0; |
fb3430ac | 1164 | TF1* ff = static_cast<TF1*>(f.GetFunctions().At(i)); |
c389303e | 1165 | // ff->SetLineColor(Color()); |
1166 | ff->SetRange(0, maxRange); | |
1167 | dist->GetListOfFunctions()->Add(new TF1(*ff)); | |
fb3430ac | 1168 | if (CheckResult(static_cast<TFitResult*>(f.GetFitResults().At(i)), |
c389303e | 1169 | relErrorCut, chi2nuCut)) { |
1170 | good[i] = ff; | |
1171 | ff->SetLineWidth(2); | |
1172 | // f.fFitResults.At(i)->Print("V"); | |
1173 | } | |
4b9857f3 | 1174 | } |
1175 | // If all else fails, use the 1 particle fit | |
fb3430ac | 1176 | TF1* ret = static_cast<TF1*>(f.GetFunctions().At(0)); |
c389303e | 1177 | |
1178 | // Find the fit with the most valid particles | |
4b9857f3 | 1179 | for (Int_t i = nFits-1; i > 0; i--) { |
1180 | if (!good[i]) continue; | |
c389303e | 1181 | if (fDebug > 1) { |
1182 | AliInfo(Form("Choosing fit with n=%d", i+1)); | |
fb3430ac | 1183 | f.GetFitResults().At(i)->Print(); |
c389303e | 1184 | } |
4b9857f3 | 1185 | ret = good[i]; |
f8715167 | 1186 | break; |
1187 | } | |
c389303e | 1188 | // Give a warning if we're using fall-back |
fb3430ac | 1189 | if (ret == f.GetFunctions().At(0)) { |
c389303e | 1190 | AliWarning("Choosing fall-back 1 particle fit"); |
0bd4b00f | 1191 | } |
c389303e | 1192 | // Copy our result and return (the functions are owned by the fitter) |
1193 | TF1* fret = new TF1(*ret); | |
1194 | return fret; | |
f8715167 | 1195 | } |
1196 | ||
1197 | //____________________________________________________________________ | |
1198 | Bool_t | |
c389303e | 1199 | AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, |
1200 | Double_t relErrorCut, | |
1201 | Double_t chi2nuCut) const | |
f8715167 | 1202 | { |
7984e5f7 | 1203 | // |
1204 | // Check the result of the fit. Returns true if | |
1205 | // - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$ | |
1206 | // - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters. Note, | |
1207 | // for multi-particle fits, this requirement is relaxed by a | |
1208 | // factor of 2 | |
1209 | // - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$ | |
1210 | // particle response | |
1211 | // | |
1212 | // Parameters: | |
1213 | // r Result to check | |
1214 | // relErrorCut Cut @f$ \delta_e@f$ applied to relative error | |
1215 | // of parameter. | |
1216 | // chi2nuCut Cut @f$ \max{\chi^2/\nu}@f$ | |
1217 | // | |
1218 | // Return: | |
1219 | // true if fit is good. | |
1220 | // | |
c389303e | 1221 | if (fDebug > 10) r->Print(); |
1222 | TString n = r->GetName(); | |
1223 | n.ReplaceAll("TFitResult-", ""); | |
4b9857f3 | 1224 | Double_t chi2 = r->Chi2(); |
1225 | Int_t ndf = r->Ndf(); | |
f8715167 | 1226 | // Double_t prob = r.Prob(); |
c389303e | 1227 | Bool_t ret = kTRUE; |
1228 | ||
1229 | // Check that the reduced chi square isn't larger than 5 | |
1230 | if (ndf <= 0 || chi2 / ndf > chi2nuCut) { | |
0bd4b00f | 1231 | if (fDebug > 2) { |
c389303e | 1232 | AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f", |
1233 | n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf), | |
0bd4b00f | 1234 | chi2nuCut)); } |
c389303e | 1235 | ret = kFALSE; |
f8715167 | 1236 | } |
1237 | ||
c389303e | 1238 | // Check each parameter |
4b9857f3 | 1239 | UShort_t nPar = r->NPar(); |
f8715167 | 1240 | for (UShort_t i = 0; i < nPar; i++) { |
c389303e | 1241 | if (i == kN) continue; // Skip the number parameter |
f8715167 | 1242 | |
c389303e | 1243 | // Get value and error and check value |
1244 | Double_t v = r->Parameter(i); | |
1245 | Double_t e = r->ParError(i); | |
f8715167 | 1246 | if (v == 0) continue; |
c389303e | 1247 | |
1248 | // Calculate the relative error and check it | |
1249 | Double_t re = e / v; | |
1250 | Double_t cut = relErrorCut * (i < kN ? 1 : 2); | |
1251 | if (re > cut) { | |
0bd4b00f | 1252 | if (fDebug > 2) { |
c389303e | 1253 | AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%", |
1254 | n.Data(), r->ParName(i).c_str(), | |
1255 | r->ParName(i).c_str(), e, v, | |
1256 | 100*(v == 0 ? 0 : e/v), | |
0bd4b00f | 1257 | 100*(cut))); } |
c389303e | 1258 | ret = kFALSE; |
f8715167 | 1259 | } |
1260 | } | |
c389303e | 1261 | |
1262 | // Check if we have scale parameters | |
1263 | if (nPar > kN) { | |
1264 | ||
1265 | // Check that the last particle has a significant contribution | |
4b9857f3 | 1266 | Double_t lastScale = r->Parameter(nPar-1); |
1267 | if (lastScale <= 1e-7) { | |
0bd4b00f | 1268 | if (fDebug) { |
c389303e | 1269 | AliWarning(Form("%s: %s=%9.6f<1e-7", |
0bd4b00f | 1270 | n.Data(), r->ParName(nPar-1).c_str(), lastScale)); } |
c389303e | 1271 | ret = kFALSE; |
f8715167 | 1272 | } |
1273 | } | |
c389303e | 1274 | return ret; |
f8715167 | 1275 | } |
1276 | ||
1277 | ||
0bd4b00f | 1278 | //__________________________________________________________________ |
1279 | void | |
fb3430ac | 1280 | AliFMDEnergyFitter::RingHistos::FindBestFits(const TList* d, |
0bd4b00f | 1281 | AliFMDCorrELossFit& obj, |
1282 | const TAxis& eta, | |
1283 | Double_t relErrorCut, | |
1284 | Double_t chi2nuCut, | |
1285 | Double_t minWeightCut) | |
1286 | { | |
7984e5f7 | 1287 | // |
1288 | // Find the best fits | |
1289 | // | |
1290 | // Parameters: | |
1291 | // d Parent list | |
1292 | // obj Object to add fits to | |
1293 | // eta Eta axis | |
1294 | // relErrorCut Cut applied to relative error of parameter. | |
1295 | // Note, for multi-particle weights, the cut | |
1296 | // is loosend by a factor of 2 | |
1297 | // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - | |
1298 | // the reduced @f$\chi^2@f$ | |
1299 | // minWeightCut Least valid @f$ a_i@f$ | |
1300 | // | |
1301 | ||
0bd4b00f | 1302 | // Get our ring list from the output container |
1303 | TList* l = GetOutputList(d); | |
1304 | if (!l) return; | |
1305 | ||
1306 | // Get the energy distributions from the output container | |
1307 | TList* dists = static_cast<TList*>(l->FindObject("EDists")); | |
1308 | if (!dists) { | |
1309 | AliWarning(Form("Didn't find %s_EtaEDists in %s", | |
1310 | fName.Data(), l->GetName())); | |
1311 | l->ls(); | |
1312 | return; | |
1313 | } | |
1314 | Int_t nBin = eta.GetNbins(); | |
1315 | ||
1316 | for (Int_t b = 1; b <= nBin; b++) { | |
1317 | TString n(Form(fgkEDistFormat, fName.Data(), b)); | |
1318 | TH1D* dist = static_cast<TH1D*>(dists->FindObject(n)); | |
1319 | // Ignore empty histograms altoghether | |
1320 | if (!dist || dist->GetEntries() <= 0) continue; | |
1321 | ||
1322 | AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist, | |
1323 | relErrorCut, | |
1324 | chi2nuCut, | |
1325 | minWeightCut); | |
1326 | best->fDet = fDet; | |
1327 | best->fRing = fRing; | |
1328 | best->fBin = b; // | |
1329 | // Double_t eta = fAxis->GetBinCenter(b); | |
1330 | obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best)); | |
1331 | } | |
1332 | } | |
1333 | ||
1334 | //__________________________________________________________________ | |
1335 | AliFMDCorrELossFit::ELossFit* | |
fb3430ac | 1336 | AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist, |
0bd4b00f | 1337 | Double_t relErrorCut, |
1338 | Double_t chi2nuCut, | |
1339 | Double_t minWeightCut) | |
1340 | { | |
7984e5f7 | 1341 | // |
1342 | // Find the best fit | |
1343 | // | |
1344 | // Parameters: | |
1345 | // dist Histogram | |
1346 | // relErrorCut Cut applied to relative error of parameter. | |
1347 | // Note, for multi-particle weights, the cut | |
1348 | // is loosend by a factor of 2 | |
1349 | // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - | |
1350 | // the reduced @f$\chi^2@f$ | |
1351 | // minWeightCut Least valid @f$ a_i@f$ | |
1352 | // | |
1353 | // Return: | |
1354 | // Best fit | |
1355 | // | |
0bd4b00f | 1356 | TList* funcs = dist->GetListOfFunctions(); |
1357 | TIter next(funcs); | |
1358 | TF1* func = 0; | |
1359 | fFits.Clear(); | |
1360 | Int_t i = 0; | |
1361 | // Info("FindBestFit", "%s", dist->GetName()); | |
1362 | while ((func = static_cast<TF1*>(next()))) { | |
1363 | AliFMDCorrELossFit::ELossFit* fit = | |
1364 | new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func); | |
1365 | fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut); | |
1366 | } | |
1367 | fFits.Sort(false); | |
1368 | // fFits.Print(); | |
1369 | return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0)); | |
1370 | } | |
1371 | ||
1372 | ||
f8715167 | 1373 | //____________________________________________________________________ |
1374 | void | |
1375 | AliFMDEnergyFitter::RingHistos::Output(TList* dir) | |
1376 | { | |
7984e5f7 | 1377 | // |
1378 | // Define outputs | |
1379 | // | |
1380 | // Parameters: | |
1381 | // dir | |
1382 | // | |
6ab100ec | 1383 | DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos"); |
f8715167 | 1384 | fList = DefineOutputList(dir); |
1385 | } | |
1386 | ||
1387 | //____________________________________________________________________ | |
1388 | // | |
1389 | // EOF | |
1390 | // |