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