New sub-structure for fitting energy loss spectra with
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDSharingFilter.cxx
CommitLineData
7e4038b5 1//
2// Class to do the sharing correction of FMD ESD data
3//
4#include "AliFMDSharingFilter.h"
5#include <AliESDFMD.h>
6#include <TAxis.h>
7#include <TList.h>
8#include <TH1.h>
9#include <TMath.h>
10#include "AliFMDAnaParameters.h"
11#include <AliLog.h>
12
13ClassImp(AliFMDSharingFilter)
14#if 0
15; // This is for Emacs
16#endif
17
18
19//____________________________________________________________________
20AliFMDSharingFilter::AliFMDSharingFilter()
21 : TNamed(),
22 fRingHistos(),
23 fLowCut(0.3),
24 fCorrectAngles(kFALSE),
25 fEtaCorr(0)
26{}
27
28//____________________________________________________________________
29AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
30 : TNamed("fmdSharingFilter", title),
31 fRingHistos(),
32 fLowCut(0.3),
33 fCorrectAngles(kFALSE),
34 fEtaCorr(0)
35{
36 fRingHistos.Add(new RingHistos(1, 'I'));
37 fRingHistos.Add(new RingHistos(2, 'I'));
38 fRingHistos.Add(new RingHistos(2, 'O'));
39 fRingHistos.Add(new RingHistos(3, 'I'));
40 fRingHistos.Add(new RingHistos(3, 'O'));
41 fEtaCorr = new TH2F("etaCorr", "Correction of #eta",
42 200, -4, 6, 200, -4, 6);
43 fEtaCorr->SetXTitle("#eta from ESD");
44 fEtaCorr->SetYTitle("#eta from AnaParameters");
45}
46
47//____________________________________________________________________
48AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
49 : TNamed(o),
50 fRingHistos(),
51 fLowCut(o.fLowCut),
52 fCorrectAngles(o.fCorrectAngles),
53 fEtaCorr(o.fEtaCorr)
54{
55 TIter next(&o.fRingHistos);
56 TObject* obj = 0;
57 while ((obj = next())) fRingHistos.Add(obj);
58}
59
60//____________________________________________________________________
61AliFMDSharingFilter::~AliFMDSharingFilter()
62{
63 fRingHistos.Delete();
64}
65
66//____________________________________________________________________
67AliFMDSharingFilter&
68AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
69{
70 SetName(o.GetName());
71 SetTitle(o.GetTitle());
72
73 fLowCut = o.fLowCut;
74 fCorrectAngles = o.fCorrectAngles;
75
76 fRingHistos.Delete();
77 TIter next(&o.fRingHistos);
78 TObject* obj = 0;
79 while ((obj = next())) fRingHistos.Add(obj);
80
81 return *this;
82}
83
84//____________________________________________________________________
85AliFMDSharingFilter::RingHistos*
86AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
87{
88 Int_t idx = -1;
89 switch (d) {
90 case 1: idx = 0; break;
91 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
92 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
93 }
94 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
95
96 return static_cast<RingHistos*>(fRingHistos.At(idx));
97}
98
99//____________________________________________________________________
100Bool_t
101AliFMDSharingFilter::Filter(const AliESDFMD& input,
102 Bool_t lowFlux,
103 AliESDFMD& output,
104 Double_t vz)
105{
106 output.Clear();
107 TIter next(&fRingHistos);
108 RingHistos* o = 0;
109 while ((o = static_cast<RingHistos*>(next())))
110 o->Clear();
111
112 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
113
114 for(UShort_t d = 1; d <= 3; d++) {
115 Int_t nRings = (d == 1 ? 1 : 2);
116 for (UShort_t q = 0; q < nRings; q++) {
117 Char_t r = (q == 0 ? 'I' : 'O');
118 UShort_t nsec = (q == 0 ? 20 : 40);
119 UShort_t nstr = (q == 0 ? 512 : 256);
120
121 RingHistos* histos = GetRingHistos(d, r);
122
123 for(UShort_t s = 0; s < nsec; s++) {
124 Bool_t usedThis = kFALSE;
125 Bool_t usedPrev = kFALSE;
126
127 for(UShort_t t = 0; t < nstr; t++) {
128 output.SetMultiplicity(d,r,s,t,0.);
129 Float_t mult = SignalInStrip(input,d,r,s,t);
130
131 // Keep dead-channel information.
132 if(mult == AliESDFMD::kInvalidMult)
133 output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
134
135 // If no signal or dead strip, go on.
136 if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
137
138 // Get the pseudo-rapidity
139 Double_t eta1 = input.Eta(d,r,s,t);
140 Double_t eta2 = pars->GetEtaFromStrip(d,r,s,t,vz);
141 fEtaCorr->Fill(eta1, eta2);
142 Double_t eta = eta2;
143
144 // Fill the diagnostics histogram
145 histos->fBefore->Fill(mult);
146
147 // Get next and previous signal - if any
148 Double_t prevE = 0;
149 Double_t nextE = 0;
150 if (t != 0) {
151 prevE = SignalInStrip(input,d,r,s,t-1);
152 if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
153 }
154 if (t != nstr - 1) {
155 nextE = SignalInStrip(input,d,r,s,t+1);
156 if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
157 }
158
159 Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
160 lowFlux,d,r,s,t,
161 usedPrev,usedThis);
162 if (mergedEnergy > fLowCut) histos->Incr();
163 histos->fAfter->Fill(mergedEnergy);
164
165 output.SetMultiplicity(d,r,s,t,mergedEnergy);
166 output.SetEta(d,r,s,t,eta);
167 } // for strip
168 } // for sector
169 } // for ring
170 } // for detector
171
172 next.Reset();
173 while ((o = static_cast<RingHistos*>(next())))
174 o->Finish();
175
176 return kTRUE;
177}
178
179//_____________________________________________________________________
180Double_t
181AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
182 UShort_t d,
183 Char_t r,
184 UShort_t s,
185 UShort_t t) const
186{
187 Double_t mult = input.Multiplicity(d,r,s,t);
188 if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
189
190 if (fCorrectAngles && !input.IsAngleCorrected())
191 mult = AngleCorrect(mult, input.Eta(d,r,s,t));
192 else if (!fCorrectAngles && input.IsAngleCorrected())
193 mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
194 return mult;
195}
196
197//_____________________________________________________________________
198Double_t
199AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
200{
201 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
202
203 // Get the high cut. The high cut is defined as the
204 // most-probably-value peak found from the energy distributions, minus
205 // 2 times the width of the corresponding Landau.
206 Double_t mpv = pars->GetMPV(d,r,eta);
207 Double_t w = pars->GetSigma(d,r,eta);
208 if (mpv > 100)
209 AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
210 Double_t highCut = (mpv - 2 * w);
211 return highCut;
212}
213
214//_____________________________________________________________________
215Double_t
216AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
217 Double_t eta,
218 Double_t prevE,
219 Double_t nextE,
220 Bool_t lowFlux,
221 UShort_t d,
222 Char_t r,
223 UShort_t /*s*/,
224 UShort_t /*t*/,
225 Bool_t& usedPrev,
226 Bool_t& usedThis)
227{
228 // IF the multiplicity is very large, or below the cut, reject it, and
229 // flag it as candidate for sharing
230 if(mult > 12 || mult < fLowCut) {
231 usedThis = kFALSE;
232 usedPrev = kFALSE;
233 return 0;
234 }
235
236 // If this was shared with the previous signal, return 0 and mark it as
237 // not a candiate for sharing
238 if(usedThis) {
239 usedThis = kFALSE;
240 usedPrev = kTRUE;
241 return 0.;
242 }
243
244 //analyse and perform sharing on one strip
245
246 // Calculate the total energy loss
247 Double_t highCut = GetHighCut(d, r, eta);
248
249 // If this signal is smaller than the next, and the next signal is smaller
250 // than then the high cut, and it's a low-flux event, then mark this and
251 // the next signal as candidates for sharing
252 //
253 // This is the test if the signal is the smaller of two possibly
254 // shared signals.
255 if (mult < nextE &&
256 nextE > highCut &&
257 lowFlux) {
258 usedThis = kFALSE;
259 usedPrev = kFALSE;
260 return 0;
261 }
262
263 Double_t totalE = mult;
264
265 // If the previous signal was larger than the low cut, and
266 // the previous was smaller than high cut, and the previous signal
267 // isn't marked as used, then add it's energy to this signal
268 if (prevE > fLowCut &&
269 prevE < highCut &&
270 !usedPrev)
271 totalE += prevE;
272
273 // If the next signal is larger than the low cut, and
274 // the next signal is smaller than the low cut, then add the next signal
275 // to this, and marked the next signal as used
276 if(nextE > fLowCut && nextE < highCut ) {
277 totalE += nextE;
278 usedThis = kTRUE;
279 }
280
281 // If we're processing on non-angle corrected data, we should do the
282 // angle correction here
283 if (!fCorrectAngles)
284 totalE = AngleCorrect(totalE, eta);
285
286 // Fill post processing histogram
287 // if(totalE > fLowCut)
288 // GetRingHistos(d,r)->fAfter->Fill(totalE);
289
290 // Return value
291 Double_t mergedEnergy = 0;
292
293 if(totalE > 0) {
294 // If we have a signal, then this is used (sharedPrev=true) and
295 // the signal is set to the result
296 mergedEnergy = totalE;
297 usedPrev = kTRUE;
298 }
299 else {
300 // If we do not have a signal (too low), then this is not shared,
301 // and the next strip is not shared either
302 usedThis = kFALSE;
303 usedPrev = kFALSE;
304 }
305
306 return mergedEnergy;
307}
308//____________________________________________________________________
309Double_t
310AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
311{
312 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
313 if (eta < 0) theta -= TMath::Pi();
314 return mult * TMath::Cos(theta);
315}
316//____________________________________________________________________
317Double_t
318AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
319{
320 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
321 if (eta < 0) theta -= TMath::Pi();
322 return mult / TMath::Cos(theta);
323}
324
325//____________________________________________________________________
326void
9d99b0dd 327AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
7e4038b5 328{
329 if (nEvents <= 0) return;
9d99b0dd 330 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
331 if (!d) return;
7e4038b5 332
333 TIter next(&fRingHistos);
334 RingHistos* o = 0;
9d99b0dd 335 while ((o = static_cast<RingHistos*>(next())))
336 o->ScaleHistograms(d, nEvents);
7e4038b5 337}
338
339//____________________________________________________________________
340void
9d99b0dd 341AliFMDSharingFilter::DefineOutput(TList* dir)
7e4038b5 342{
343 TList* d = new TList;
344 d->SetName(GetName());
345 dir->Add(d);
346 d->Add(fEtaCorr);
347 TIter next(&fRingHistos);
348 RingHistos* o = 0;
349 while ((o = static_cast<RingHistos*>(next()))) {
350 o->Output(d);
351 }
352}
353
354//====================================================================
355AliFMDSharingFilter::RingHistos::RingHistos()
9d99b0dd 356 : AliForwardUtil::RingHistos(),
7e4038b5 357 fBefore(0),
358 fAfter(0),
359 fHits(0),
360 fNHits(0)
361{}
362
363//____________________________________________________________________
364AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 365 : AliForwardUtil::RingHistos(d,r),
7e4038b5 366 fBefore(0),
367 fAfter(0),
368 fHits(0),
369 fNHits(0)
370{
9d99b0dd 371 fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
372 Form("Energy loss in %s (reconstruction)", fName.Data()),
7e4038b5 373 1000, 0, 25);
9d99b0dd 374 fAfter = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
375 Form("Energy loss in %s (sharing corrected)",
376 fName.Data()), 1000, 0, 25);
7e4038b5 377 fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
378 fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
379 fBefore->SetFillColor(kRed+1);
380 fBefore->SetFillStyle(3001);
381 // fBefore->Sumw2();
382 fBefore->SetDirectory(0);
383 fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
384 fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
385 fAfter->SetFillColor(kBlue+1);
386 fAfter->SetFillStyle(3001);
387 // fAfter->Sumw2();
388 fAfter->SetDirectory(0);
389
9d99b0dd 390 fHits = new TH1D(Form("%s_Hits", fName.Data()),
391 Form("Number of hits in %s", fName.Data()),
7e4038b5 392 200, 0, 200000);
393 fHits->SetDirectory(0);
394 fHits->SetXTitle("# of hits");
395 fHits->SetFillColor(kGreen+1);
396}
397//____________________________________________________________________
398AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd 399 : AliForwardUtil::RingHistos(o),
7e4038b5 400 fBefore(o.fBefore),
401 fAfter(o.fAfter),
402 fHits(o.fHits),
403 fNHits(o.fNHits)
404{}
405
406//____________________________________________________________________
407AliFMDSharingFilter::RingHistos&
408AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
409{
9d99b0dd 410 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 411 fDet = o.fDet;
412 fRing = o.fRing;
413
414 if (fBefore) delete fBefore;
415 if (fAfter) delete fAfter;
416 if (fHits) delete fHits;
417
418 fBefore = static_cast<TH1D*>(o.fBefore->Clone());
419 fAfter = static_cast<TH1D*>(o.fAfter->Clone());
420 fHits = static_cast<TH1D*>(o.fHits->Clone());
421
422 return *this;
423}
424//____________________________________________________________________
425AliFMDSharingFilter::RingHistos::~RingHistos()
426{
427 if (fBefore) delete fBefore;
428 if (fAfter) delete fAfter;
429 if (fHits) delete fHits;
430}
431//____________________________________________________________________
432void
433AliFMDSharingFilter::RingHistos::Finish()
434{
435 fHits->Fill(fNHits);
436}
437
9d99b0dd 438//____________________________________________________________________
439void
440AliFMDSharingFilter::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
441{
442 TList* l = GetOutputList(dir);
443 if (!l) return;
444
445 TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
446 fName.Data())));
447 TH1D* after = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
448 fName.Data())));
449 if (before) before->Scale(1./nEvents);
450 if (after) after->Scale(1./nEvents);
451}
452
7e4038b5 453//____________________________________________________________________
454void
455AliFMDSharingFilter::RingHistos::Output(TList* dir)
456{
9d99b0dd 457 TList* d = DefineOutputList(dir);
7e4038b5 458 d->Add(fBefore);
459 d->Add(fAfter);
460 d->Add(fHits);
461 dir->Add(d);
462}
463
464//____________________________________________________________________
465//
466// EOF
467//