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