Adding the option to zero shared entries below threshold
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDSharingFilter.cxx
1 //
2 // Class to do the sharing correction.  That is, a filter that merges 
3 // adjacent strip signals presumably originating from a single particle 
4 // that impinges on the detector in such a way that it deposite energy 
5 // into two or more strips. 
6 //
7 // Input: 
8 //    - AliESDFMD object  - from reconstruction
9 //
10 // Output: 
11 //    - AliESDFMD object  - copy of input, but with signals merged 
12 //
13 // Corrections used: 
14 //    - AliFMDCorrELossFit
15 //
16 // Histograms: 
17 //    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of 
18 //      signals before and after the filter.  
19 //    - For each ring (see above), an array of distributions of number of 
20 //      hit strips for each vertex bin (if enabled - see Init method)
21 // 
22 //
23 //
24 #include "AliFMDSharingFilter.h"
25 #include <AliESDFMD.h>
26 #include <TAxis.h>
27 #include <TList.h>
28 #include <TH1.h>
29 #include <TMath.h>
30 #include "AliForwardCorrectionManager.h"
31 #include "AliFMDCorrELossFit.h"
32 #include <AliLog.h>
33 #include <TROOT.h>
34 #include <THStack.h>
35 #include <iostream>
36 #include <iomanip>
37
38 ClassImp(AliFMDSharingFilter)
39 #if 0
40 ; // This is for Emacs
41 #endif 
42
43 #define DBG(L,M) \
44   do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false) 
45 #define DBGL(L,M) \
46   do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false) 
47                               
48
49
50 //____________________________________________________________________
51 AliFMDSharingFilter::AliFMDSharingFilter()
52   : TNamed(), 
53     fRingHistos(),
54     fLowCut(0.),
55     fCorrectAngles(kFALSE), 
56     fNXi(1),
57     fIncludeSigma(true),
58     fSummed(0),
59     fHighCuts(0),
60     fOper(0),
61     fDebug(0),
62     fZeroSharedHitsBelowThreshold(0)
63 {
64   // 
65   // Default Constructor - do not use 
66   //
67 }
68
69 //____________________________________________________________________
70 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
71   : TNamed("fmdSharingFilter", title), 
72     fRingHistos(), 
73     fLowCut(0.),
74     fCorrectAngles(kFALSE), 
75     fNXi(1),
76     fIncludeSigma(true),
77     fSummed(0),
78     fHighCuts(0),
79     fOper(0),
80     fDebug(0),
81     fZeroSharedHitsBelowThreshold(0)
82 {
83   // 
84   // Constructor 
85   // 
86   // Parameters:
87   //    title Title of object  - not significant 
88   //
89   fRingHistos.Add(new RingHistos(1, 'I'));
90   fRingHistos.Add(new RingHistos(2, 'I'));
91   fRingHistos.Add(new RingHistos(2, 'O'));
92   fRingHistos.Add(new RingHistos(3, 'I'));
93   fRingHistos.Add(new RingHistos(3, 'O'));
94 }
95
96 //____________________________________________________________________
97 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
98   : TNamed(o), 
99     fRingHistos(), 
100     fLowCut(o.fLowCut),
101     fCorrectAngles(o.fCorrectAngles), 
102     fNXi(o.fNXi),
103     fIncludeSigma(o.fIncludeSigma),
104     fSummed(o.fSummed),
105     fHighCuts(o.fHighCuts),
106     fOper(o.fOper),
107     fDebug(o.fDebug),
108     fZeroSharedHitsBelowThreshold(o.fZeroSharedHitsBelowThreshold)
109 {
110   // 
111   // Copy constructor 
112   // 
113   // Parameters:
114   //    o Object to copy from 
115   //
116   TIter    next(&o.fRingHistos);
117   TObject* obj = 0;
118   while ((obj = next())) fRingHistos.Add(obj);
119 }
120
121 //____________________________________________________________________
122 AliFMDSharingFilter::~AliFMDSharingFilter()
123 {
124   // 
125   // Destructor
126   //
127   fRingHistos.Delete();
128 }
129
130 //____________________________________________________________________
131 AliFMDSharingFilter&
132 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
133 {
134   // 
135   // Assignment operator 
136   // 
137   // Parameters:
138   //    o Object to assign from 
139   // 
140   // Return:
141   //    Reference to this 
142   //
143   TNamed::operator=(o);
144
145   fLowCut                       = o.fLowCut;
146   fCorrectAngles                = o.fCorrectAngles;
147   fNXi                          = o.fNXi;
148   fDebug                        = o.fDebug;
149   fOper                         = o.fOper;
150   fSummed                       = o.fSummed;
151   fHighCuts                     = o.fHighCuts;
152   fIncludeSigma                 = o.fIncludeSigma;
153   fZeroSharedHitsBelowThreshold = o.fZeroSharedHitsBelowThreshold;
154
155   fRingHistos.Delete();
156   TIter    next(&o.fRingHistos);
157   TObject* obj = 0;
158   while ((obj = next())) fRingHistos.Add(obj);
159   
160   return *this;
161 }
162
163 //____________________________________________________________________
164 AliFMDSharingFilter::RingHistos*
165 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
166 {
167   // 
168   // Get the ring histogram container 
169   // 
170   // Parameters:
171   //    d Detector
172   //    r Ring 
173   // 
174   // Return:
175   //    Ring histogram container 
176   //
177   Int_t idx = -1;
178   switch (d) { 
179   case 1: idx = 0; break;
180   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
181   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
182   }
183   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
184   
185   return static_cast<RingHistos*>(fRingHistos.At(idx));
186 }
187
188 //____________________________________________________________________
189 void
190 AliFMDSharingFilter::Init()
191 {
192   // Initialise 
193   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::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   AliFMDCorrELossFit* fits = fcm.GetELossFit();
199   const TAxis& eAxis = fits->GetEtaAxis();
200
201   UShort_t nEta = eAxis.GetNbins();
202   fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
203   fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
204   fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
205   fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
206   fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
207   fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
208
209   UShort_t ybin = 0;
210   for (UShort_t d = 1; d <= 3; d++) {
211     UShort_t nr = (d == 1 ? 1 : 2);
212     for (UShort_t q = 0; q < nr; q++) { 
213       Char_t r = (q == 0 ? 'I' : 'O');
214       ybin++;
215       for (UShort_t e = 1; e <= nEta; e++) { 
216         Double_t eta = eAxis.GetBinCenter(e);
217         Double_t cut = GetHighCut(d, r, eta, false);
218         if (cut <= 0) continue;
219         fHighCuts->SetBinContent(e, ybin, cut);
220       }
221     }
222   }
223 }
224
225 //____________________________________________________________________
226 Bool_t
227 AliFMDSharingFilter::Filter(const AliESDFMD& input, 
228                             Bool_t           lowFlux,
229                             AliESDFMD&       output)
230 {
231   // 
232   // Filter the input AliESDFMD object
233   // 
234   // Parameters:
235   //    input     Input 
236   //    lowFlux   If this is a low-flux event 
237   //    output    Output AliESDFMD object 
238   // 
239   // Return:
240   //    True on success, false otherwise 
241   //
242   output.Clear();
243   TIter    next(&fRingHistos);
244   RingHistos* o      = 0;
245   while ((o = static_cast<RingHistos*>(next())))
246     o->Clear();
247
248   if (fOper) fOper->Reset(0);
249   Int_t nNone      = 0;
250   Int_t nCandidate = 0;
251   Int_t nMerged    = 0;
252   Int_t nSummed    = 0;
253
254   Status status[512];
255   
256   for(UShort_t d = 1; d <= 3; d++) {
257     Int_t nRings = (d == 1 ? 1 : 2);
258     for (UShort_t q = 0; q < nRings; q++) {
259       Char_t      r      = (q == 0 ? 'I' : 'O');
260       UShort_t    nsec   = (q == 0 ?  20 :  40);
261       UShort_t    nstr   = (q == 0 ? 512 : 256);
262       RingHistos* histos = GetRingHistos(d, r);
263       
264       for(UShort_t s = 0; s < nsec;  s++) {
265         for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate;
266 #ifdef USE_OLDER_MERGING
267         Bool_t usedThis   = kFALSE;
268         Bool_t usedPrev   = kFALSE;
269 #endif  
270         for(UShort_t t = 0; t < nstr; t++) {
271           output.SetMultiplicity(d,r,s,t,0.);
272           Float_t mult = SignalInStrip(input,d,r,s,t);
273           // Get the pseudo-rapidity 
274           Double_t eta = input.Eta(d,r,s,t);
275           Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
276           if (s == 0) output.SetEta(d,r,s,t,eta);
277           
278           // Keep dead-channel information. 
279           if(mult == AliESDFMD::kInvalidMult)
280             output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
281
282           // If no signal or dead strip, go on. 
283           if (mult == AliESDFMD::kInvalidMult || mult == 0) {
284             if (mult == 0) histos->fSum->Fill(eta,phi,mult);
285             status[t] = kNone;
286             continue;
287           }
288
289           // Fill the diagnostics histogram 
290           histos->fBefore->Fill(mult);
291
292           // Get next and previous signal - if any 
293           Double_t prevE = 0;
294           Double_t nextE = 0;
295           Status   prevStatus = (t == 0      ? kNone : status[t-1]);
296           Status   thisStatus = status[t];
297           Status   nextStatus = (t == nstr-1 ? kNone : status[t+1]);
298           if (t != 0) {
299             prevE = SignalInStrip(input,d,r,s,t-1);
300             if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
301           }
302           if (t != nstr - 1) {
303             nextE = SignalInStrip(input,d,r,s,t+1);
304             if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
305           }
306           if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult);
307           
308 #ifdef USE_OLDER_MERGING
309           Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
310                                                       lowFlux,d,r,s,t, 
311                                                       usedPrev,usedThis);
312           status[t] = (usedPrev ? kMergedWithOther : kNone);
313           if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone);
314 #else 
315           Double_t mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, 
316                                                       eta, lowFlux, 
317                                                       d, r, s, t, 
318                                                       prevStatus, 
319                                                       thisStatus, 
320                                                       nextStatus);
321           if (t != 0)      status[t-1] = prevStatus;
322           if (t != nstr-1) status[t+1] = nextStatus;
323           status[t] = thisStatus;
324
325 #endif
326           // If we're processing on non-angle corrected data, we
327           // should do the angle correction here
328           if (!fCorrectAngles)
329             mergedEnergy = AngleCorrect(mergedEnergy, eta);
330           if (mergedEnergy > 0) histos->Incr();
331
332           if (t != 0) 
333             histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1), 
334                                           mergedEnergy);
335           histos->fBeforeAfter->Fill(mult, mergedEnergy);
336           histos->fAfter->Fill(mergedEnergy);
337           histos->fSum->Fill(eta,phi,mergedEnergy);
338
339           output.SetMultiplicity(d,r,s,t,mergedEnergy);
340         } // for strip
341         for (UShort_t t = 0; t < nstr; t++) {
342           if (fOper) fOper->operator()(d, r, s, t) = status[t];
343           switch (status[t]) { 
344           case kNone:            nNone++;      break;
345           case kCandidate:       nCandidate++; break;
346           case kMergedWithOther: nMerged++;    break;
347           case kMergedInto:      nSummed++;    break;
348           }
349         }
350       } // for sector
351     } // for ring 
352   } // for detector
353   fSummed->Fill(kNone,            nNone);
354   fSummed->Fill(kCandidate,       nCandidate);
355   fSummed->Fill(kMergedWithOther, nMerged);
356   fSummed->Fill(kMergedInto,      nSummed);
357
358   DBGL(1, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d", 
359                nNone, nCandidate, nMerged, nSummed));
360   next.Reset();
361   while ((o = static_cast<RingHistos*>(next())))
362     o->Finish();
363
364   return kTRUE;
365 }
366
367 //_____________________________________________________________________
368 Double_t 
369 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input, 
370                                    UShort_t         d,
371                                    Char_t           r,
372                                    UShort_t         s,
373                                    UShort_t         t) const
374 {
375   // 
376   // Get the signal in a strip 
377   // 
378   // Parameters:
379   //    fmd   ESD object
380   //    d     Detector
381   //    r     Ring 
382   //    s     Sector 
383   //    t     Strip
384   // 
385   // Return:
386   //    The energy signal 
387   //
388   Double_t mult = input.Multiplicity(d,r,s,t);
389   // In case of 
390   //  - bad value (invalid or 0) 
391   //  - we want angle corrected and data is 
392   //  - we don't want angle corrected and data isn't 
393   // just return read value  
394   if (mult == AliESDFMD::kInvalidMult               || 
395       mult == 0                                     ||
396       (fCorrectAngles && input.IsAngleCorrected()) || 
397       (!fCorrectAngles && !input.IsAngleCorrected()))
398     return mult;
399
400   // If we want angle corrected data, correct it, 
401   // otherwise de-correct it 
402   if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
403   else                mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
404   return mult;
405 }
406 //_____________________________________________________________________
407 Double_t 
408 AliFMDSharingFilter::GetLowCut(UShort_t, Char_t, Double_t) const
409 {
410   //
411   // Get the low cut.  Normally, the low cut is taken to be the lower
412   // value of the fit range used when generating the energy loss fits.
413   // However, if fLowCut is set (using SetLowCit) to a value greater
414   // than 0, then that value is used.
415   //
416   if (fLowCut > 0) return fLowCut;
417   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
418   AliFMDCorrELossFit* fits = fcm.GetELossFit();
419   return fits->GetLowCut();
420 }
421                         
422 //_____________________________________________________________________
423 Double_t 
424 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, 
425                                 Double_t eta, Bool_t errors) const
426 {
427   //
428   // Get the high cut.  The high cut is defined as the 
429   // most-probably-value peak found from the energy distributions, minus 
430   // 2 times the width of the corresponding Landau.
431   //
432   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
433
434  
435   // Get the high cut.  The high cut is defined as the 
436   // most-probably-value peak found from the energy distributions, minus 
437   // 2 times the width of the corresponding Landau.
438   AliFMDCorrELossFit* fits = fcm.GetELossFit();
439   
440   return fits->GetLowerBound(d, r, eta, fNXi, errors, fIncludeSigma);
441 }
442
443 //_____________________________________________________________________
444 namespace { 
445   const char* status2String(AliFMDSharingFilter::Status s)
446   {
447     switch(s) { 
448     case AliFMDSharingFilter::kNone:            return "none     "; 
449     case AliFMDSharingFilter::kCandidate:       return "candidate"; 
450     case AliFMDSharingFilter::kMergedWithOther: return "merged   "; 
451     case AliFMDSharingFilter::kMergedInto:      return "summed   "; 
452     }
453     return "unknown  ";
454   }
455 }
456
457 //_____________________________________________________________________
458 Double_t 
459 AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE,
460                                          Double_t prevE,
461                                          Double_t nextE,
462                                          Double_t eta,
463                                          Bool_t   lowFlux,
464                                          UShort_t d,
465                                          Char_t   r,
466                                          UShort_t s,
467                                          UShort_t t,
468                                          Status&  prevStatus, 
469                                          Status&  thisStatus, 
470                                          Status&  nextStatus) const
471 {
472   Double_t lowCut = GetLowCut(d, r, eta);
473
474   DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)", 
475              d, r, s, t,
476              thisE, status2String(thisStatus), 
477              prevE, status2String(prevStatus), 
478              nextE, status2String(nextStatus)));
479
480   // If below cut, then modify zero signal and make sure the next
481   // strip is considered a candidate.
482   if (thisE < lowCut || thisE > 20) { 
483     thisStatus = kNone;
484     DBGL(3,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE));
485     if (prevStatus == kCandidate) prevStatus = kNone;
486     return 0;  
487   }
488   // It this strip was merged with the previous strip, then 
489   // make the next strip a candidate and zero the value in this strip. 
490   if (thisStatus == kMergedWithOther) { 
491     DBGL(3,Form(" Merged with other, 0'ed"));
492     return 0;
493   }
494
495   // Get the upper cut on the sharing 
496   Double_t highCut = GetHighCut(d, r, eta ,false);
497   if (highCut <= 0) {
498     thisStatus = kNone;
499     return 0;
500   }
501
502   // Only for low-flux  events 
503   if (lowFlux) { 
504     // If this is smaller than the next, and the next is larger 
505     // then the cut, mark both as candidates for sharing. 
506     // They should be be merged on the next iteration 
507     if (thisE < nextE && nextE > highCut) { 
508       thisStatus = kCandidate;
509       nextStatus = kCandidate;
510       return 0;
511     }
512   }
513   
514   // Variable to sum signal in 
515   Double_t totalE = thisE;
516
517   // If the previous signal was between the two cuts, and is still
518   // marked as a candidate , then merge it in here.
519   if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) {
520     totalE     += prevE;
521     prevStatus =  kMergedWithOther;
522     DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f", 
523                 lowCut, prevE, highCut, totalE));
524   }
525
526   // If the next signal is between the two cuts, then merge it here
527   if (nextE > lowCut && nextE < highCut) { 
528     totalE     += nextE;
529     nextStatus =  kMergedWithOther;
530     DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE));
531   }
532   
533   // If the signal is bigger than the high-cut, then return the value 
534   if (totalE > highCut) { 
535     if (totalE > thisE) 
536       thisStatus = kMergedInto;
537     else 
538       thisStatus = kNone;
539     DBGL(3, Form(" %9f>%f9 (was %9f) -> %9f %s", 
540                  totalE, highCut, thisE, totalE,status2String(thisStatus)));
541     return totalE;
542   }
543   else {
544     // This is below cut, so flag next as a candidate 
545     DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut));
546     nextStatus = kCandidate;
547   }
548   
549   // If the total signal is smaller than low cut then zero this and kill this 
550   if (totalE < lowCut)  {
551     DBGL(3, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE));
552     thisStatus = kNone;
553     return 0;
554   }
555
556   // If total signal not above high cut or lower than low cut, 
557   // mark this as a candidate for merging into the next, and zero signal 
558   DBGL(3, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate", 
559                    lowCut, totalE, highCut, thisE));
560   thisStatus = kCandidate;
561   if(fZeroSharedHitsBelowThreshold)
562     return 0;
563   else return totalE; 
564 }
565            
566 //_____________________________________________________________________
567 Double_t 
568 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
569                                          Double_t eta,
570                                          Double_t prevE,
571                                          Double_t nextE,
572                                          Bool_t   lowFlux,
573                                          UShort_t d,
574                                          Char_t   r,
575                                          UShort_t /*s*/,
576                                          UShort_t /*t*/,
577                                          Bool_t&  usedPrev, 
578                                          Bool_t&  usedThis) const
579 {
580   // 
581   // The actual algorithm 
582   // 
583   // Parameters:
584   //    mult      The unfiltered signal in the strip
585   //    eta       Psuedo rapidity 
586   //    prevE     Previous strip signal (or 0)
587   //    nextE     Next strip signal (or 0) 
588   //    lowFlux   Whether this is a low flux event 
589   //    d         Detector
590   //    r         Ring 
591   //    s         Sector 
592   //    t         Strip
593   //    usedPrev  Whether the previous strip was used in sharing or not
594   //    usedThis  Wether this strip was used in sharing or not. 
595   // 
596   // Return:
597   //    The filtered signal in the strip
598   //
599
600   // IF the multiplicity is very large, or below the cut, reject it, and 
601   // flag it as candidate for sharing 
602   Double_t    lowCut = GetLowCut(d,r,eta);
603   if(mult > 12 || mult < lowCut) {
604     usedThis      = kFALSE;
605     usedPrev      = kFALSE;
606     return 0;
607   }
608
609   // If this was shared with the previous signal, return 0 and mark it as 
610   // not a candiate for sharing 
611   if(usedThis) {
612     usedThis      = kFALSE;
613     usedPrev      = kTRUE;
614     return 0.;
615   }
616
617   //analyse and perform sharing on one strip
618   
619   // Get the high cut 
620   Double_t highCut = GetHighCut(d, r, eta, false);
621   if (highCut <= 0) {
622     usedThis = kFALSE;
623     usedPrev = kTRUE;
624     return 0;
625   }
626
627   // If this signal is smaller than the next, and the next signal is smaller 
628   // than then the high cut, and it's a low-flux event, then mark this and 
629   // the next signal as candidates for sharing 
630   // 
631   // This is the test if the signal is the smaller of two possibly
632   // shared signals.   
633   if (mult  < nextE   && 
634       nextE > highCut && 
635       lowFlux) {
636     usedThis      = kFALSE;
637     usedPrev      = kFALSE;
638     return 0;
639   }
640   
641   Double_t totalE  = mult;
642   
643   // If the previous signal was larger than the low cut, and 
644   // the previous was smaller than high cut, and the previous signal 
645   // isn't marked as used, then add it's energy to this signal 
646   if (prevE > lowCut && 
647       prevE < highCut && 
648       !usedPrev) 
649     totalE += prevE;
650
651   // If the next signal is larger than the low cut, and 
652   // the next signal is smaller than the high cut, then add the next signal 
653   // to this, and marked the next signal as used 
654   if(nextE > lowCut && nextE < highCut ) {
655     totalE      += nextE;
656     usedThis =  kTRUE;
657   }
658
659   // If we're processing on non-angle corrected data, we should do the 
660   // angle correction here 
661   // if (!fCorrectAngles) 
662   //   totalE = AngleCorrect(totalE, eta);
663
664   // Fill post processing histogram
665   // if(totalE > fLowCut)
666   //   GetRingHistos(d,r)->fAfter->Fill(totalE);
667
668   // Return value 
669   Double_t mergedEnergy = 0;
670   
671   if(totalE > 0) {
672     // If we have a signal, then this is used (sharedPrev=true) and
673     // the signal is set to the result
674     mergedEnergy = totalE;
675     usedPrev     = kTRUE;
676   }
677   else {
678     // If we do not have a signal (too low), then this is not shared, 
679     // and the next strip is not shared either 
680     usedThis  = kFALSE;
681     usedPrev  = kFALSE;
682   }
683   
684   return mergedEnergy; 
685 }
686 //____________________________________________________________________
687 Double_t
688 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
689 {
690   // 
691   // Angle correct the signal 
692   // 
693   // Parameters:
694   //    mult Angle Un-corrected Signal 
695   //    eta  Pseudo-rapidity 
696   // 
697   // Return:
698   //    Angle corrected signal 
699   //
700   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
701   if (eta < 0) theta -= TMath::Pi();
702   return mult * TMath::Cos(theta);
703 }
704 //____________________________________________________________________
705 Double_t
706 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
707 {
708   // 
709   // Angle de-correct the signal 
710   // 
711   // Parameters:
712   //    mult Angle corrected Signal 
713   //    eta  Pseudo-rapidity 
714   // 
715   // Return:
716   //    Angle un-corrected signal 
717   //
718   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
719   if (eta < 0) theta -= TMath::Pi();
720   return mult / TMath::Cos(theta);
721 }
722
723 //____________________________________________________________________
724 void
725 AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
726 {
727   // 
728   // Scale the histograms to the total number of events 
729   // 
730   // Parameters:
731   //    dir     Where the output is 
732   //    nEvents Number of events 
733   //
734   if (nEvents <= 0) return;
735   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
736   if (!d) return;
737
738   TIter    next(&fRingHistos);
739   RingHistos* o = 0;
740   THStack* sums = new THStack("sums", "Sum of ring signals");
741   while ((o = static_cast<RingHistos*>(next()))) {
742     o->ScaleHistograms(d, nEvents);
743     TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
744     sum->Scale(1., "width");
745     sum->SetTitle(o->GetName());
746     sum->SetDirectory(0);
747     sum->SetYTitle("#sum #Delta/#Delta_{mip}");
748     sums->Add(sum);
749   }
750   d->Add(sums);
751 }
752
753 //____________________________________________________________________
754 void
755 AliFMDSharingFilter::DefineOutput(TList* dir)
756 {
757   // 
758   // Define the output histograms.  These are put in a sub list of the
759   // passed list.   The histograms are merged before the parent task calls 
760   // AliAnalysisTaskSE::Terminate 
761   // 
762   // Parameters:
763   //    dir Directory to add to 
764   //
765   TList* d = new TList;
766   d->SetName(GetName());
767   dir->Add(d);
768
769   fSummed = new TH2I("operations", "Strip operations", 
770                      kMergedInto, kNone-.5, kMergedInto+.5,
771                      51201, -.5, 51200.5);
772   fSummed->SetXTitle("Operation");
773   fSummed->SetYTitle("# of strips");
774   fSummed->SetZTitle("Events");
775   fSummed->GetXaxis()->SetBinLabel(kNone,            "None");
776   fSummed->GetXaxis()->SetBinLabel(kCandidate,       "Candidate");
777   fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
778   fSummed->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
779   fSummed->SetDirectory(0);
780   d->Add(fSummed);
781
782   fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
783   fHighCuts->SetXTitle("#eta");
784   fHighCuts->SetDirectory(0);
785   d->Add(fHighCuts);
786
787
788   TIter    next(&fRingHistos);
789   RingHistos* o = 0;
790   while ((o = static_cast<RingHistos*>(next()))) {
791     o->Output(d);
792   }
793 }
794 //____________________________________________________________________
795 void
796 AliFMDSharingFilter::Print(Option_t* /*option*/) const
797 {
798   // 
799   // Print information
800   // 
801   // Parameters:
802   //    option Not used 
803   //
804   char ind[gROOT->GetDirLevel()+1];
805   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
806   ind[gROOT->GetDirLevel()] = '\0';
807   std::cout << ind << ClassName() << ": " << GetName() << '\n'
808             << std::boolalpha 
809             << ind << " Debug:                  " << fDebug << "\n"
810             << ind << " Low cut:                " << fLowCut << '\n'
811             << ind << " N xi factor:            " << fNXi    << '\n'
812             << ind << " Include sigma in cut:   " << fIncludeSigma << '\n'
813             << ind << " Use corrected angles:   " << fCorrectAngles 
814             << std::noboolalpha << std::endl;
815 }
816   
817 //====================================================================
818 AliFMDSharingFilter::RingHistos::RingHistos()
819   : AliForwardUtil::RingHistos(), 
820     fBefore(0), 
821     fAfter(0), 
822     fBeforeAfter(0),
823     fNeighborsBefore(0),
824     fNeighborsAfter(0),
825     fSum(0),
826     fHits(0),
827     fNHits(0)
828 {
829   // 
830   // Default CTOR
831   //
832
833 }
834
835 //____________________________________________________________________
836 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
837   : AliForwardUtil::RingHistos(d,r), 
838     fBefore(0), 
839     fAfter(0),
840     fBeforeAfter(0),
841     fNeighborsBefore(0),
842     fNeighborsAfter(0),
843     fSum(0),
844     fHits(0),
845     fNHits(0)
846 {
847   // 
848   // Constructor
849   // 
850   // Parameters:
851   //    d detector
852   //    r ring 
853   //
854   fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)", 
855                      1000, 0, 25);
856   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
857   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
858   fBefore->SetFillColor(Color());
859   fBefore->SetFillStyle(3001);
860   fBefore->SetDirectory(0);
861
862   fAfter  = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
863   fAfter->SetTitle("Energy loss in %s (sharing corrected)");
864   fAfter->SetFillColor(Color()+2);
865   fAfter->SetDirectory(0);
866
867   Double_t max = 15;
868   Double_t min = -1;
869   Int_t    n   = int((max-min) / (max / 300));
870   fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation", 
871                           n, min, max, n, min, max);
872   fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
873   fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
874   fBeforeAfter->SetZTitle("Correlation");
875   fBeforeAfter->SetDirectory(0);
876
877   fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
878   fNeighborsBefore->SetTitle("Correlation of neighbors befire");
879   fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
880   fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
881   fNeighborsBefore->SetDirectory(0);
882   
883   fNeighborsAfter = 
884     static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
885   fNeighborsAfter->SetTitle("Correlation of neighbors after");
886   fNeighborsAfter->SetDirectory(0);
887
888   fSum = new TH2D("summed", "Summed signal", 200, -4, 6, 
889                   (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi());
890   fSum->SetDirectory(0);
891   fSum->Sumw2();
892   fSum->SetMarkerColor(Color());
893   // fSum->SetFillColor(Color());
894   fSum->SetXTitle("#eta");
895   fSum->SetYTitle("#varphi [radians]");
896   fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) ");
897   
898   fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
899   fHits->SetDirectory(0);
900   fHits->SetXTitle("# of hits");
901   fHits->SetFillColor(kGreen+1);
902 }
903 //____________________________________________________________________
904 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
905   : AliForwardUtil::RingHistos(o), 
906     fBefore(o.fBefore), 
907     fAfter(o.fAfter),
908     fBeforeAfter(o.fBeforeAfter),
909     fNeighborsBefore(o.fNeighborsBefore),
910     fNeighborsAfter(o.fNeighborsAfter),
911     fSum(o.fSum),
912     fHits(o.fHits),
913     fNHits(o.fNHits)
914 {
915   // 
916   // Copy constructor 
917   // 
918   // Parameters:
919   //    o Object to copy from 
920   //
921 }
922
923 //____________________________________________________________________
924 AliFMDSharingFilter::RingHistos&
925 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
926 {
927   // 
928   // Assignment operator 
929   // 
930   // Parameters:
931   //    o Object to assign from 
932   // 
933   // Return:
934   //    Reference to this 
935   //
936   AliForwardUtil::RingHistos::operator=(o);
937   fDet = o.fDet;
938   fRing = o.fRing;
939   
940   if (fBefore) delete  fBefore;
941   if (fAfter)  delete  fAfter;
942   if (fHits)   delete fHits;
943   
944   fBefore          = static_cast<TH1D*>(o.fBefore->Clone());
945   fAfter           = static_cast<TH1D*>(o.fAfter->Clone());
946   fBeforeAfter     = static_cast<TH2D*>(o.fBeforeAfter->Clone());
947   fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
948   fNeighborsAfter  = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
949   fHits            = static_cast<TH1D*>(o.fHits->Clone());
950   fSum             = static_cast<TH2D*>(o.fSum->Clone());
951
952   return *this;
953 }
954 //____________________________________________________________________
955 AliFMDSharingFilter::RingHistos::~RingHistos()
956 {
957   // 
958   // Destructor 
959   //
960 }
961 //____________________________________________________________________
962 void
963 AliFMDSharingFilter::RingHistos::Finish()
964 {
965   // 
966   // Finish off 
967   // 
968   //
969   fHits->Fill(fNHits);
970 }
971
972 //____________________________________________________________________
973 void
974 AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir, 
975                                                  Int_t nEvents)
976 {
977   // 
978   // Scale the histograms to the total number of events 
979   // 
980   // Parameters:
981   //    nEvents Number of events 
982   //    dir     Where the output is 
983   //
984   TList* l = GetOutputList(dir);
985   if (!l) return; 
986
987   TH1D* before = static_cast<TH1D*>(l->FindObject("esdELoss"));
988   TH1D* after  = static_cast<TH1D*>(l->FindObject("anaELoss"));
989   TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
990   if (before) before->Scale(1./nEvents);
991   if (after)  after->Scale(1./nEvents);
992   if (summed) summed->Scale(1./nEvents);
993   
994 }
995
996 //____________________________________________________________________
997 void
998 AliFMDSharingFilter::RingHistos::Output(TList* dir)
999 {
1000   // 
1001   // Make output 
1002   // 
1003   // Parameters:
1004   //    dir where to store 
1005   //
1006   TList* d = DefineOutputList(dir);
1007
1008   d->Add(fBefore);
1009   d->Add(fAfter);
1010   d->Add(fBeforeAfter);
1011   d->Add(fNeighborsBefore);
1012   d->Add(fNeighborsAfter);
1013   d->Add(fHits);
1014   d->Add(fSum);
1015
1016   dir->Add(d);
1017 }
1018
1019 //____________________________________________________________________
1020 //
1021 // EOF
1022 //