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