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