More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDSharingFilter.cxx
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 "AliForwardCorrectionManager.h"
11 // #include "AliFMDAnaParameters.h"
12 #include <AliLog.h>
13 #include <TROOT.h>
14 #include <iostream>
15 #include <iomanip>
16
17 ClassImp(AliFMDSharingFilter)
18 #if 0
19 ; // This is for Emacs
20 #endif 
21
22
23 //____________________________________________________________________
24 AliFMDSharingFilter::AliFMDSharingFilter()
25   : TNamed(), 
26     fRingHistos(),
27     fLowCut(0.),
28     fCorrectAngles(kFALSE), 
29     fNXi(1),
30     fDebug(0)
31 {}
32
33 //____________________________________________________________________
34 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
35   : TNamed("fmdSharingFilter", title), 
36     fRingHistos(), 
37     fLowCut(0.),
38     fCorrectAngles(kFALSE), 
39     fNXi(1),
40     fDebug(0)
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'));
47 }
48
49 //____________________________________________________________________
50 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
51   : TNamed(o), 
52     fRingHistos(), 
53     fLowCut(o.fLowCut),
54     fCorrectAngles(o.fCorrectAngles), 
55     fNXi(o.fNXi),
56     fDebug(o.fDebug)
57 {
58   TIter    next(&o.fRingHistos);
59   TObject* obj = 0;
60   while ((obj = next())) fRingHistos.Add(obj);
61 }
62
63 //____________________________________________________________________
64 AliFMDSharingFilter::~AliFMDSharingFilter()
65 {
66   fRingHistos.Delete();
67 }
68
69 //____________________________________________________________________
70 AliFMDSharingFilter&
71 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
72 {
73   TNamed::operator=(o);
74
75   fLowCut        = o.fLowCut;
76   fCorrectAngles = o.fCorrectAngles;
77   fNXi           = o.fNXi;
78   fDebug         = o.fDebug;
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 //____________________________________________________________________
89 AliFMDSharingFilter::RingHistos*
90 AliFMDSharingFilter::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 //____________________________________________________________________
104 Bool_t
105 AliFMDSharingFilter::Filter(const AliESDFMD& input, 
106                             Bool_t           lowFlux,
107                             AliESDFMD&       output)
108 {
109   output.Clear();
110   TIter    next(&fRingHistos);
111   RingHistos* o      = 0;
112   Double_t    lowCut = GetLowCut();
113   while ((o = static_cast<RingHistos*>(next())))
114     o->Clear();
115
116   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
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 
143           Double_t eta = input.Eta(d,r,s,t);
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);
163           if (mergedEnergy > lowCut) histos->Incr();
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 //_____________________________________________________________________
181 Double_t 
182 AliFMDSharingFilter::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 }
197 //_____________________________________________________________________
198 Double_t 
199 AliFMDSharingFilter::GetLowCut() const
200 {
201   if (fLowCut > 0) return fLowCut;
202   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
203   AliFMDCorrELossFit* fits = fcm.GetELossFit();
204   return fits->GetLowCut();
205 }
206                         
207 //_____________________________________________________________________
208 Double_t 
209 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
210 {
211   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
212
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.
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);
232   return highCut;
233 }
234            
235 //_____________________________________________________________________
236 Double_t 
237 AliFMDSharingFilter::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 
251   Double_t    lowCut = GetLowCut();
252   if(mult > 12 || mult < lowCut) {
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   
268   // Get the high cut 
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 
290   if (prevE > lowCut && 
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 
298   if(nextE > lowCut && nextE < highCut ) {
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 //____________________________________________________________________
331 Double_t
332 AliFMDSharingFilter::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 //____________________________________________________________________
339 Double_t
340 AliFMDSharingFilter::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 //____________________________________________________________________
348 void
349 AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
350 {
351   if (nEvents <= 0) return;
352   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
353   if (!d) return;
354
355   TIter    next(&fRingHistos);
356   RingHistos* o = 0;
357   while ((o = static_cast<RingHistos*>(next())))
358     o->ScaleHistograms(d, nEvents);
359 }
360
361 //____________________________________________________________________
362 void
363 AliFMDSharingFilter::DefineOutput(TList* dir)
364 {
365   TList* d = new TList;
366   d->SetName(GetName());
367   dir->Add(d);
368   TIter    next(&fRingHistos);
369   RingHistos* o = 0;
370   while ((o = static_cast<RingHistos*>(next()))) {
371     o->Output(d);
372   }
373 }
374 //____________________________________________________________________
375 void
376 AliFMDSharingFilter::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 }
387   
388 //====================================================================
389 AliFMDSharingFilter::RingHistos::RingHistos()
390   : AliForwardUtil::RingHistos(), 
391     fBefore(0), 
392     fAfter(0), 
393     fHits(0),
394     fNHits(0)
395 {}
396
397 //____________________________________________________________________
398 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
399   : AliForwardUtil::RingHistos(d,r), 
400     fBefore(0), 
401     fAfter(0),
402     fHits(0),
403     fNHits(0)
404 {
405   fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()), 
406                      Form("Energy loss in %s (reconstruction)", fName.Data()), 
407                      1000, 0, 25);
408   fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()), 
409                      Form("Energy loss in %s (sharing corrected)",
410                           fName.Data()), 1000, 0, 25);
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
424   fHits = new TH1D(Form("%s_Hits", fName.Data()), 
425                    Form("Number of hits in %s", fName.Data()), 
426                    200, 0, 200000);
427   fHits->SetDirectory(0);
428   fHits->SetXTitle("# of hits");
429   fHits->SetFillColor(kGreen+1);
430 }
431 //____________________________________________________________________
432 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
433   : AliForwardUtil::RingHistos(o), 
434     fBefore(o.fBefore), 
435     fAfter(o.fAfter),
436     fHits(o.fHits),
437     fNHits(o.fNHits)
438 {}
439
440 //____________________________________________________________________
441 AliFMDSharingFilter::RingHistos&
442 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
443 {
444   AliForwardUtil::RingHistos::operator=(o);
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 //____________________________________________________________________
459 AliFMDSharingFilter::RingHistos::~RingHistos()
460 {
461   if (fBefore) delete fBefore;
462   if (fAfter)  delete fAfter;
463   if (fHits)   delete fHits;
464 }
465 //____________________________________________________________________
466 void
467 AliFMDSharingFilter::RingHistos::Finish()
468 {
469   fHits->Fill(fNHits);
470 }
471
472 //____________________________________________________________________
473 void
474 AliFMDSharingFilter::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
487 //____________________________________________________________________
488 void
489 AliFMDSharingFilter::RingHistos::Output(TList* dir)
490 {
491   TList* d = DefineOutputList(dir);
492   d->Add(fBefore);
493   d->Add(fAfter);
494   d->Add(fHits);
495   dir->Add(d);
496 }
497
498 //____________________________________________________________________
499 //
500 // EOF
501 //