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 //