New implementation of the forward multiplicity analysis.
[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
327AliFMDSharingFilter::ScaleHistograms(Int_t nEvents)
328{
329 if (nEvents <= 0) return;
330
331 TIter next(&fRingHistos);
332 RingHistos* o = 0;
333 while ((o = static_cast<RingHistos*>(next()))) {
334 o->fBefore->Scale(1. / nEvents);
335 o->fAfter->Scale(1. / nEvents);
336 }
337}
338
339//____________________________________________________________________
340void
341AliFMDSharingFilter::Output(TList* dir)
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()
356 : fDet(0),
357 fRing('\0'),
358 fBefore(0),
359 fAfter(0),
360 fHits(0),
361 fNHits(0)
362{}
363
364//____________________________________________________________________
365AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
366 : fDet(d),
367 fRing(r),
368 fBefore(0),
369 fAfter(0),
370 fHits(0),
371 fNHits(0)
372{
373 fBefore = new TH1D(Form("FMD%d%c_ESD_Eloss", d, r),
374 Form("Energy loss in FMD%d%c (reconstruction)", d, r),
375 1000, 0, 25);
376 fAfter = new TH1D(Form("FMD%d%c_Ana_Eloss", d, r),
377 Form("Energy loss in FMD%d%c (sharing corrected)", d, r),
378 1000, 0, 25);
379 fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
380 fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
381 fBefore->SetFillColor(kRed+1);
382 fBefore->SetFillStyle(3001);
383 // fBefore->Sumw2();
384 fBefore->SetDirectory(0);
385 fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
386 fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
387 fAfter->SetFillColor(kBlue+1);
388 fAfter->SetFillStyle(3001);
389 // fAfter->Sumw2();
390 fAfter->SetDirectory(0);
391
392 fHits = new TH1D(Form("FMD%d%c_Hits", d, r),
393 Form("Number of hits in FMD%d%c", d, r),
394 200, 0, 200000);
395 fHits->SetDirectory(0);
396 fHits->SetXTitle("# of hits");
397 fHits->SetFillColor(kGreen+1);
398}
399//____________________________________________________________________
400AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
401 : TObject(o),
402 fDet(o.fDet),
403 fRing(o.fRing),
404 fBefore(o.fBefore),
405 fAfter(o.fAfter),
406 fHits(o.fHits),
407 fNHits(o.fNHits)
408{}
409
410//____________________________________________________________________
411AliFMDSharingFilter::RingHistos&
412AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
413{
414 fDet = o.fDet;
415 fRing = o.fRing;
416
417 if (fBefore) delete fBefore;
418 if (fAfter) delete fAfter;
419 if (fHits) delete fHits;
420
421 fBefore = static_cast<TH1D*>(o.fBefore->Clone());
422 fAfter = static_cast<TH1D*>(o.fAfter->Clone());
423 fHits = static_cast<TH1D*>(o.fHits->Clone());
424
425 return *this;
426}
427//____________________________________________________________________
428AliFMDSharingFilter::RingHistos::~RingHistos()
429{
430 if (fBefore) delete fBefore;
431 if (fAfter) delete fAfter;
432 if (fHits) delete fHits;
433}
434//____________________________________________________________________
435void
436AliFMDSharingFilter::RingHistos::Finish()
437{
438 fHits->Fill(fNHits);
439}
440
441//____________________________________________________________________
442void
443AliFMDSharingFilter::RingHistos::Output(TList* dir)
444{
445 TList* d = new TList;
446 d->SetName(Form("FMD%d%c", fDet, fRing));
447 d->Add(fBefore);
448 d->Add(fAfter);
449 d->Add(fHits);
450 dir->Add(d);
451}
452
453//____________________________________________________________________
454//
455// EOF
456//