Added debug flags
[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 "AliFMDAnaParameters.h"
11 #include <AliLog.h>
12
13 ClassImp(AliFMDSharingFilter)
14 #if 0
15 ; // This is for Emacs
16 #endif 
17
18
19 //____________________________________________________________________
20 AliFMDSharingFilter::AliFMDSharingFilter()
21   : TNamed(), 
22     fRingHistos(),
23     fLowCut(0.3),
24     fCorrectAngles(kFALSE), 
25     fDebug(0)
26 {}
27
28 //____________________________________________________________________
29 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
30   : TNamed("fmdSharingFilter", title), 
31     fRingHistos(), 
32     fLowCut(0.3),
33     fCorrectAngles(kFALSE), 
34     fDebug(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 }
42
43 //____________________________________________________________________
44 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
45   : TNamed(o), 
46     fRingHistos(), 
47     fLowCut(o.fLowCut),
48     fCorrectAngles(o.fCorrectAngles), 
49     fDebug(o.fDebug)
50 {
51   TIter    next(&o.fRingHistos);
52   TObject* obj = 0;
53   while ((obj = next())) fRingHistos.Add(obj);
54 }
55
56 //____________________________________________________________________
57 AliFMDSharingFilter::~AliFMDSharingFilter()
58 {
59   fRingHistos.Delete();
60 }
61
62 //____________________________________________________________________
63 AliFMDSharingFilter&
64 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
65 {
66   TNamed::operator=(o);
67
68   fLowCut        = o.fLowCut;
69   fCorrectAngles = o.fCorrectAngles;
70   fDebug         = o.fDebug;
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 //____________________________________________________________________
81 AliFMDSharingFilter::RingHistos*
82 AliFMDSharingFilter::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 //____________________________________________________________________
96 Bool_t
97 AliFMDSharingFilter::Filter(const AliESDFMD& input, 
98                             Bool_t           lowFlux,
99                             AliESDFMD&       output)
100 {
101   output.Clear();
102   TIter    next(&fRingHistos);
103   RingHistos* o = 0;
104   while ((o = static_cast<RingHistos*>(next())))
105     o->Clear();
106
107   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
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 
134           Double_t eta = input.Eta(d,r,s,t);
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 //_____________________________________________________________________
172 Double_t 
173 AliFMDSharingFilter::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 //_____________________________________________________________________
190 Double_t 
191 AliFMDSharingFilter::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 //_____________________________________________________________________
207 Double_t 
208 AliFMDSharingFilter::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 //____________________________________________________________________
301 Double_t
302 AliFMDSharingFilter::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 //____________________________________________________________________
309 Double_t
310 AliFMDSharingFilter::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 //____________________________________________________________________
318 void
319 AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
320 {
321   if (nEvents <= 0) return;
322   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
323   if (!d) return;
324
325   TIter    next(&fRingHistos);
326   RingHistos* o = 0;
327   while ((o = static_cast<RingHistos*>(next())))
328     o->ScaleHistograms(d, nEvents);
329 }
330
331 //____________________________________________________________________
332 void
333 AliFMDSharingFilter::DefineOutput(TList* dir)
334 {
335   TList* d = new TList;
336   d->SetName(GetName());
337   dir->Add(d);
338   TIter    next(&fRingHistos);
339   RingHistos* o = 0;
340   while ((o = static_cast<RingHistos*>(next()))) {
341     o->Output(d);
342   }
343 }
344   
345 //====================================================================
346 AliFMDSharingFilter::RingHistos::RingHistos()
347   : AliForwardUtil::RingHistos(), 
348     fBefore(0), 
349     fAfter(0), 
350     fHits(0),
351     fNHits(0)
352 {}
353
354 //____________________________________________________________________
355 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
356   : AliForwardUtil::RingHistos(d,r), 
357     fBefore(0), 
358     fAfter(0),
359     fHits(0),
360     fNHits(0)
361 {
362   fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()), 
363                      Form("Energy loss in %s (reconstruction)", fName.Data()), 
364                      1000, 0, 25);
365   fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()), 
366                      Form("Energy loss in %s (sharing corrected)",
367                           fName.Data()), 1000, 0, 25);
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
381   fHits = new TH1D(Form("%s_Hits", fName.Data()), 
382                    Form("Number of hits in %s", fName.Data()), 
383                    200, 0, 200000);
384   fHits->SetDirectory(0);
385   fHits->SetXTitle("# of hits");
386   fHits->SetFillColor(kGreen+1);
387 }
388 //____________________________________________________________________
389 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
390   : AliForwardUtil::RingHistos(o), 
391     fBefore(o.fBefore), 
392     fAfter(o.fAfter),
393     fHits(o.fHits),
394     fNHits(o.fNHits)
395 {}
396
397 //____________________________________________________________________
398 AliFMDSharingFilter::RingHistos&
399 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
400 {
401   AliForwardUtil::RingHistos::operator=(o);
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 //____________________________________________________________________
416 AliFMDSharingFilter::RingHistos::~RingHistos()
417 {
418   if (fBefore) delete fBefore;
419   if (fAfter)  delete fAfter;
420   if (fHits)   delete fHits;
421 }
422 //____________________________________________________________________
423 void
424 AliFMDSharingFilter::RingHistos::Finish()
425 {
426   fHits->Fill(fNHits);
427 }
428
429 //____________________________________________________________________
430 void
431 AliFMDSharingFilter::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
444 //____________________________________________________________________
445 void
446 AliFMDSharingFilter::RingHistos::Output(TList* dir)
447 {
448   TList* d = DefineOutputList(dir);
449   d->Add(fBefore);
450   d->Add(fAfter);
451   d->Add(fHits);
452   dir->Add(d);
453 }
454
455 //____________________________________________________________________
456 //
457 // EOF
458 //