PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx

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