26176494dd34d877b201aac3e327dee1b5d74158
[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     fEtaCorr(0)
26 {}
27
28 //____________________________________________________________________
29 AliFMDSharingFilter::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 //____________________________________________________________________
48 AliFMDSharingFilter::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 //____________________________________________________________________
61 AliFMDSharingFilter::~AliFMDSharingFilter()
62 {
63   fRingHistos.Delete();
64 }
65
66 //____________________________________________________________________
67 AliFMDSharingFilter&
68 AliFMDSharingFilter::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 //____________________________________________________________________
85 AliFMDSharingFilter::RingHistos*
86 AliFMDSharingFilter::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 //____________________________________________________________________
100 Bool_t
101 AliFMDSharingFilter::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 //_____________________________________________________________________
180 Double_t 
181 AliFMDSharingFilter::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 //_____________________________________________________________________
198 Double_t 
199 AliFMDSharingFilter::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 //_____________________________________________________________________
215 Double_t 
216 AliFMDSharingFilter::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 //____________________________________________________________________
309 Double_t
310 AliFMDSharingFilter::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 //____________________________________________________________________
317 Double_t
318 AliFMDSharingFilter::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 //____________________________________________________________________
326 void
327 AliFMDSharingFilter::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 //____________________________________________________________________
340 void
341 AliFMDSharingFilter::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 //====================================================================
355 AliFMDSharingFilter::RingHistos::RingHistos()
356   : fDet(0),
357     fRing('\0'),
358     fBefore(0), 
359     fAfter(0), 
360     fHits(0),
361     fNHits(0)
362 {}
363
364 //____________________________________________________________________
365 AliFMDSharingFilter::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 //____________________________________________________________________
400 AliFMDSharingFilter::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 //____________________________________________________________________
411 AliFMDSharingFilter::RingHistos&
412 AliFMDSharingFilter::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 //____________________________________________________________________
428 AliFMDSharingFilter::RingHistos::~RingHistos()
429 {
430   if (fBefore) delete fBefore;
431   if (fAfter)  delete fAfter;
432   if (fHits)   delete fHits;
433 }
434 //____________________________________________________________________
435 void
436 AliFMDSharingFilter::RingHistos::Finish()
437 {
438   fHits->Fill(fNHits);
439 }
440
441 //____________________________________________________________________
442 void
443 AliFMDSharingFilter::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 //