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