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::SetupForData(const TAxis& axis)
 286 {
 287   // Initialise - called on first event
 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             mult = AliESDFMD::kInvalidMult;
 434           }
 435
 436           // If no signal or dead strip, go on.
 437           if (mult == AliESDFMD::kInvalidMult || mult == 0) {
 438             if (mult == 0) histos->fSum->Fill(eta,phi,mult);
 439             status[t] = kNone;
 440             continue;
 441           }
 442
 443           // Fill the diagnostics histogram
 444           histos->fBefore->Fill(mult);
 445
 446           Double_t mergedEnergy = 0;
 447
 448           if(fUseSimpleMerging) {
 449             Float_t etot = 0;
 450             if (t < nstr-1) histos->fNeighborsBefore->Fill(mult,multNext);
 451             if(mult > GetHighCut(d, r, eta ,false)) {
 452               histos->fDistanceBefore->Fill(nDistanceBefore);
 453               nDistanceBefore = -1;
 454             }
 455
 456             if(eTotal > 0) {
 457               if(fThreeStripSharing && multNext > GetLowCut(d, r, eta) &&
 458                  (multNext < GetHighCut(d, r, eta ,false) || twoLow)) {
 459                 eTotal = eTotal + multNext;
 460                 used = kTRUE;
 461                 histos->fTriple->Fill(eTotal);
 462                 nTriple++;
 463                 twoLow = kFALSE;
 464               }
 465               else {
 466                 used = kFALSE;
 467                 histos->fDouble->Fill(eTotal);
 468                 nDouble++;
 469               }
 470               etot   = eTotal;
 471               eTotal = -1;
 472             }
 473             else {
 474               if(used) {used = kFALSE; continue; }
 475               if(mult > GetLowCut(d, r, eta)) etot = mult;
 476
 477               if(mult > GetLowCut(d, r, eta) &&
 478                  multNext > GetLowCut(d, r, eta) &&
 479                  (mult < GetHighCut(d, r, eta ,false) ||
 480                   multNext < GetHighCut(d, r, eta ,false))) {
 481
 482                 if(mult < GetHighCut(d, r, eta ,false) &&
 483                    multNext < GetHighCut(d, r, eta ,false) )
 484                   twoLow = kTRUE;
 485
 486                 if(mult>multNext && multNextNext < GetLowCut(d, r, eta))
 487                   {
 488                     etot = mult + multNext;
 489                     used=kTRUE;
 490                     histos->fDouble->Fill(etot);
 491                     nDouble++;
 492                   }
 493                 else {
 494                   etot   = 0;
 495                   eTotal = mult + multNext;
 496                 }
 497               }
 498               else {
 499                 if(etot > 0) {
 500                   histos->fSingle->Fill(etot);
 501                   histos->fSinglePerStrip->Fill(etot,t);
 502                   nSingle++;
 503                 }
 504               }
 505             }
 506
 507             mergedEnergy = etot;
 508             if(mergedEnergy > GetHighCut(d, r, eta ,false) ) {
 509               histos->fDistanceAfter->Fill(nDistanceAfter);
 510               nDistanceAfter    = -1;
 511             }
 512             //if(mult>0 && multNext >0)
 513             //  std::cout<<mult<<"  "<<multNext<<"  "<<mergedEnergy<<std::endl;
 514           } // End of simple merge
 515           else { // Not simple
 516             // Get next and previous signal - if any
 517             Double_t prevE = 0;
 518             Double_t nextE = 0;
 519             Status   prevStatus = (t == 0      ? kNone : status[t-1]);
 520             Status   thisStatus = status[t];
 521             Status   nextStatus = (t == nstr-1 ? kNone : status[t+1]);
 522             if (t != 0) {
 523               prevE = SignalInStrip(input,d,r,s,t-1);
 524               if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
 525             }
 526             if (t != nstr - 1) {
 527               nextE = SignalInStrip(input,d,r,s,t+1);
 528               if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
 529             }
 530             if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult);
 531
 532             mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE,
 533                                                eta, lowFlux,
 534                                                d, r, s, t,
 535                                                prevStatus,
 536                                                thisStatus,
 537                                                nextStatus);
 538             if (t != 0)      status[t-1] = prevStatus;
 539             if (t != nstr-1) status[t+1] = nextStatus;
 540             status[t] = thisStatus;
 541             // If we're processing on non-angle corrected data, we
 542             // should do the angle correction here
 543           } // End of non-simple
 544           if (!fCorrectAngles)
 545             mergedEnergy = AngleCorrect(mergedEnergy, eta);
 546           if (mergedEnergy > 0) histos->Incr();
 547
 548           if (t != 0)
 549             histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1),
 550                                           mergedEnergy);
 551           histos->fBeforeAfter->Fill(mult, mergedEnergy);
 552           if(mergedEnergy > 0)
 553             histos->fAfter->Fill(mergedEnergy);
 554           histos->fSum->Fill(eta,phi,mergedEnergy);
 555
 556           output.SetMultiplicity(d,r,s,t,mergedEnergy);
 557         } // for strip
 558         for (UShort_t t = 0; t < nstr; t++) {
 559           if (fOper) fOper->operator()(d, r, s, t) = status[t];
 560           switch (status[t]) {
 561           case kNone:            nNone++;      break;
 562           case kCandidate:       nCandidate++; break;
 563           case kMergedWithOther: nMerged++;    break;
 564           case kMergedInto:      nSummed++;    break;
 565           }
 566         }
 567       } // for sector
 568     } // for ring
 569   } // for detector
 570   fSummed->Fill(kNone,            nNone);
 571   fSummed->Fill(kCandidate,       nCandidate);
 572   fSummed->Fill(kMergedWithOther, nMerged);
 573   fSummed->Fill(kMergedInto,      nSummed);
 574
 575   DMSG(fDebug, 3,"single=%9d, double=%9d, triple=%9d",
 576        nSingle, nDouble, nTriple);
 577   next.Reset();
 578   while ((o = static_cast<RingHistos*>(next())))
 579     o->Finish();
 580
 581   return kTRUE;
 582 }
 583
 584 //_____________________________________________________________________
 585 Double_t
 586 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
 587                                    UShort_t         d,
 588                                    Char_t           r,
 589                                    UShort_t         s,
 590                                    UShort_t         t) const
 591 {
 592   //
 593   // Get the signal in a strip
 594   //
 595   // Parameters:
 596   //    fmd   ESD object
 597   //    d     Detector
 598   //    r     Ring
 599   //    s     Sector
 600   //    t     Strip
 601   //
 602   // Return:
 603   //    The energy signal
 604   //
 605   Double_t mult = input.Multiplicity(d,r,s,t);
 606   // In case of
 607   //  - bad value (invalid or 0)
 608   //  - we want angle corrected and data is
 609   //  - we don't want angle corrected and data isn't
 610   // just return read value
 611   if (mult == AliESDFMD::kInvalidMult               ||
 612       mult == 0                                     ||
 613       (fCorrectAngles && input.IsAngleCorrected()) ||
 614       (!fCorrectAngles && !input.IsAngleCorrected()))
 615     return mult;
 616
 617   // If we want angle corrected data, correct it,
 618   // otherwise de-correct it
 619   if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
 620   else                mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
 621   return mult;
 622 }
 623 //_____________________________________________________________________
 624 Double_t
 625 AliFMDSharingFilter::GetLowCut(UShort_t d, Char_t r, Double_t eta) const
 626 {
 627   //
 628   // Get the low cut.  Normally, the low cut is taken to be the lower
 629   // value of the fit range used when generating the energy loss fits.
 630   // However, if fLowCut is set (using SetLowCit) to a value greater
 631   // than 0, then that value is used.
 632   //
 633   return fLCuts.GetMultCut(d,r,eta,false);
 634 #if 0
 635   if (!fCutAtFractionOfMPV && fLowCut > 0) return fLowCut;
 636
 637   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 638   AliFMDCorrELossFit* fits = fcm.GetELossFit();
 639
 640   if (fCutAtFractionOfMPV) {
 641     AliFMDCorrELossFit::ELossFit* func = fits->GetFit(d,r,eta);
 642     return fFractionOfMPV*func->GetDelta() ;
 643   }
 644   return fits->GetLowCut();
 645 #endif
 646 }
 647
 648 //_____________________________________________________________________
 649 Double_t
 650 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r,
 651                                 Double_t eta, Bool_t errors) const
 652 {
 653   //
 654   // Get the high cut.  The high cut is defined as the
 655   // most-probably-value peak found from the energy distributions, minus
 656   // 2 times the width of the corresponding Landau.
 657   //
 658   return fHCuts.GetMultCut(d,r,eta,errors);
 659 }
 660
 661 //_____________________________________________________________________
 662 namespace {
 663   const char* status2String(AliFMDSharingFilter::Status s)
 664   {
 665     switch(s) {
 666     case AliFMDSharingFilter::kNone:            return "none     ";
 667     case AliFMDSharingFilter::kCandidate:       return "candidate";
 668     case AliFMDSharingFilter::kMergedWithOther: return "merged   ";
 669     case AliFMDSharingFilter::kMergedInto:      return "summed   ";
 670     }
 671     return "unknown  ";
 672   }
 673 }
 674
 675 //_____________________________________________________________________
 676 Double_t
 677 AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE,
 678                                          Double_t prevE,
 679                                          Double_t nextE,
 680                                          Double_t eta,
 681                                          Bool_t   lowFlux,
 682                                          UShort_t d,
 683                                          Char_t   r,
 684                                          UShort_t s,
 685                                          UShort_t t,
 686                                          Status&  prevStatus,
 687                                          Status&  thisStatus,
 688                                          Status&  nextStatus) const
 689 {
 690   Double_t lowCut = GetLowCut(d, r, eta);
 691
 692   DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)",
 693              d, r, s, t,
 694              thisE, status2String(thisStatus),
 695              prevE, status2String(prevStatus),
 696              nextE, status2String(nextStatus)));
 697
 698   // If below cut, then modify zero signal and make sure the next
 699   // strip is considered a candidate.
 700   if (thisE < lowCut || thisE > 20) {
 701     thisStatus = kNone;
 702     DBGL(5,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE));
 703     if (prevStatus == kCandidate) prevStatus = kNone;
 704     return 0;
 705   }
 706   // It this strip was merged with the previous strip, then
 707   // make the next strip a candidate and zero the value in this strip.
 708   if (thisStatus == kMergedWithOther) {
 709     DBGL(5,Form(" Merged with other, 0'ed"));
 710     return 0;
 711   }
 712
 713   // Get the upper cut on the sharing
 714   Double_t highCut = GetHighCut(d, r, eta ,false);
 715   if (highCut <= 0) {
 716     thisStatus = kNone;
 717     return 0;
 718   }
 719
 720   // Only for low-flux  events
 721   if (lowFlux) {
 722     // If this is smaller than the next, and the next is larger
 723     // then the cut, mark both as candidates for sharing.
 724     // They should be be merged on the next iteration
 725     if (thisE < nextE && nextE > highCut) {
 726       thisStatus = kCandidate;
 727       nextStatus = kCandidate;
 728       return 0;
 729     }
 730   }
 731
 732   // Variable to sum signal in
 733   Double_t totalE = thisE;
 734
 735   // If the previous signal was between the two cuts, and is still
 736   // marked as a candidate , then merge it in here.
 737   if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) {
 738     totalE     += prevE;
 739     prevStatus =  kMergedWithOther;
 740     DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f",
 741                 lowCut, prevE, highCut, totalE));
 742   }
 743
 744   // If the next signal is between the two cuts, then merge it here
 745   if (nextE > lowCut && nextE < highCut) {
 746     totalE     += nextE;
 747     nextStatus =  kMergedWithOther;
 748     DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE));
 749   }
 750
 751   // If the signal is bigger than the high-cut, then return the value
 752   if (totalE > highCut) {
 753     if (totalE > thisE)
 754       thisStatus = kMergedInto;
 755     else
 756       thisStatus = kNone;
 757     DBGL(5, Form(" %9f>%f9 (was %9f) -> %9f %s",
 758                  totalE, highCut, thisE, totalE,status2String(thisStatus)));
 759     return totalE;
 760   }
 761   else {
 762     // This is below cut, so flag next as a candidate
 763     DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut));
 764     nextStatus = kCandidate;
 765   }
 766
 767   // If the total signal is smaller than low cut then zero this and kill this
 768   if (totalE < lowCut)  {
 769     DBGL(5, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE));
 770     thisStatus = kNone;
 771     return 0;
 772   }
 773
 774   // If total signal not above high cut or lower than low cut,
 775   // mark this as a candidate for merging into the next, and zero signal
 776   DBGL(5, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate",
 777                    lowCut, totalE, highCut, thisE));
 778   thisStatus = kCandidate;
 779   return (fZeroSharedHitsBelowThreshold ? 0 : totalE);
 780 }
 781
 782 //_____________________________________________________________________
 783 Double_t
 784 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 785                                          Double_t eta,
 786                                          Double_t prevE,
 787                                          Double_t nextE,
 788                                          Bool_t   lowFlux,
 789                                          UShort_t d,
 790                                          Char_t   r,
 791                                          UShort_t /*s*/,
 792                                          UShort_t /*t*/,
 793                                          Bool_t&  usedPrev,
 794                                          Bool_t&  usedThis) const
 795 {
 796   //
 797   // The actual algorithm
 798   //
 799   // Parameters:
 800   //    mult      The unfiltered signal in the strip
 801   //    eta       Psuedo rapidity
 802   //    prevE     Previous strip signal (or 0)
 803   //    nextE     Next strip signal (or 0)
 804   //    lowFlux   Whether this is a low flux event
 805   //    d         Detector
 806   //    r         Ring
 807   //    s         Sector
 808   //    t         Strip
 809   //    usedPrev  Whether the previous strip was used in sharing or not
 810   //    usedThis  Wether this strip was used in sharing or not.
 811   //
 812   // Return:
 813   //    The filtered signal in the strip
 814   //
 815
 816   // IF the multiplicity is very large, or below the cut, reject it, and
 817   // flag it as candidate for sharing
 818   Double_t    lowCut = GetLowCut(d,r,eta);
 819   if(mult > 12 || mult < lowCut) {
 820     usedThis      = kFALSE;
 821     usedPrev      = kFALSE;
 822     return 0;
 823   }
 824
 825   // If this was shared with the previous signal, return 0 and mark it as
 826   // not a candiate for sharing
 827   if(usedThis) {
 828     usedThis      = kFALSE;
 829     usedPrev      = kTRUE;
 830     return 0.;
 831   }
 832
 833   //analyse and perform sharing on one strip
 834
 835   // Get the high cut
 836   Double_t highCut = GetHighCut(d, r, eta, false);
 837   if (highCut <= 0) {
 838     usedThis = kFALSE;
 839     usedPrev = kTRUE;
 840     return 0;
 841   }
 842
 843   // If this signal is smaller than the next, and the next signal is smaller
 844   // than then the high cut, and it's a low-flux event, then mark this and
 845   // the next signal as candidates for sharing
 846   //
 847   // This is the test if the signal is the smaller of two possibly
 848   // shared signals.
 849   if (mult  < nextE   &&
 850       nextE > highCut &&
 851       lowFlux) {
 852     usedThis      = kFALSE;
 853     usedPrev      = kFALSE;
 854     return 0;
 855   }
 856
 857   Double_t totalE  = mult;
 858
 859   // If the previous signal was larger than the low cut, and
 860   // the previous was smaller than high cut, and the previous signal
 861   // isn't marked as used, then add it's energy to this signal
 862   if (prevE > lowCut &&
 863       prevE < highCut &&
 864       !usedPrev)
 865     totalE += prevE;
 866
 867   // If the next signal is larger than the low cut, and
 868   // the next signal is smaller than the high cut, then add the next signal
 869   // to this, and marked the next signal as used
 870   if(nextE > lowCut && nextE < highCut ) {
 871     totalE      += nextE;
 872     usedThis =  kTRUE;
 873   }
 874
 875   // If we're processing on non-angle corrected data, we should do the
 876   // angle correction here
 877   // if (!fCorrectAngles)
 878   //   totalE = AngleCorrect(totalE, eta);
 879
 880   // Fill post processing histogram
 881   // if(totalE > fLowCut)
 882   //   GetRingHistos(d,r)->fAfter->Fill(totalE);
 883
 884   // Return value
 885   Double_t mergedEnergy = 0;
 886
 887   if(totalE > 0) {
 888     // If we have a signal, then this is used (sharedPrev=true) and
 889     // the signal is set to the result
 890     mergedEnergy = totalE;
 891     usedPrev     = kTRUE;
 892   }
 893   else {
 894     // If we do not have a signal (too low), then this is not shared,
 895     // and the next strip is not shared either
 896     usedThis  = kFALSE;
 897     usedPrev  = kFALSE;
 898   }
 899
 900   return mergedEnergy;
 901 }
 902 //____________________________________________________________________
 903 Double_t
 904 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
 905 {
 906   //
 907   // Angle correct the signal
 908   //
 909   // Parameters:
 910   //    mult Angle Un-corrected Signal
 911   //    eta  Pseudo-rapidity
 912   //
 913   // Return:
 914   //    Angle corrected signal
 915   //
 916   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 917   if (eta < 0) theta -= TMath::Pi();
 918   return mult * TMath::Cos(theta);
 919 }
 920 //____________________________________________________________________
 921 Double_t
 922 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
 923 {
 924   //
 925   // Angle de-correct the signal
 926   //
 927   // Parameters:
 928   //    mult Angle corrected Signal
 929   //    eta  Pseudo-rapidity
 930   //
 931   // Return:
 932   //    Angle un-corrected signal
 933   //
 934   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 935   if (eta < 0) theta -= TMath::Pi();
 936   return mult / TMath::Cos(theta);
 937 }
 938
 939 //____________________________________________________________________
 940 void
 941 AliFMDSharingFilter::Terminate(const TList* dir, TList* output, Int_t nEvents)
 942 {
 943   //
 944   // Scale the histograms to the total number of events
 945   //
 946   // Parameters:
 947   //    dir     Where the output is
 948   //    nEvents Number of events
 949   //
 950   DGUARD(fDebug,1, "Scale histograms in AliFMDSharingFilter");
 951   if (nEvents <= 0) return;
 952   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 953   if (!d) return;
 954
 955   TList* out = new TList;
 956   out->SetName(d->GetName());
 957   out->SetOwner();
 958
 959   TIter    next(&fRingHistos);
 960   RingHistos* o = 0;
 961   THStack* sums = new THStack("sums", "Sum of ring signals");
 962   while ((o = static_cast<RingHistos*>(next()))) {
 963     o->Terminate(d, nEvents);
 964     if (!o->fSum) {
 965       Warning("Terminate", "No sum histogram found for ring %s", o->GetName());
 966       continue;
 967     }
 968     TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
 969     sum->Scale(1., "width");
 970     sum->SetTitle(o->GetName());
 971     sum->SetDirectory(0);
 972     sum->SetYTitle("#sum #Delta/#Delta_{mip}");
 973     sums->Add(sum);
 974   }
 975   out->Add(sums);
 976   output->Add(out);
 977 }
 978
 979 //____________________________________________________________________
 980 void
 981 AliFMDSharingFilter::CreateOutputObjects(TList* dir)
 982 {
 983   //
 984   // Define the output histograms.  These are put in a sub list of the
 985   // passed list.   The histograms are merged before the parent task calls
 986   // AliAnalysisTaskSE::Terminate
 987   //
 988   // Parameters:
 989   //    dir Directory to add to
 990   //
 991   DGUARD(fDebug,1, "Define output in AliFMDSharingFilter");
 992   TList* d = new TList;
 993   d->SetName(GetName());
 994   dir->Add(d);
 995
 996   fSummed = new TH2I("operations", "Strip operations",
 997                      kMergedInto, kNone-.5, kMergedInto+.5,
 998                      51201, -.5, 51200.5);
 999   fSummed->SetXTitle("Operation");
1000   fSummed->SetYTitle("# of strips");
1001   fSummed->SetZTitle("Events");
1002   fSummed->GetXaxis()->SetBinLabel(kNone,            "None");
1003   fSummed->GetXaxis()->SetBinLabel(kCandidate,       "Candidate");
1004   fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
1005   fSummed->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
1006   fSummed->SetDirectory(0);
1007   d->Add(fSummed);
1008
1009   fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
1010   fHighCuts->SetXTitle("#eta");
1011   fHighCuts->SetDirectory(0);
1012   d->Add(fHighCuts);
1013
1014   fLowCuts = new TH2D("lowCuts", "Low cuts used", 1,0,1, 1,0,1);
1015   fLowCuts->SetXTitle("#eta");
1016   fLowCuts->SetDirectory(0);
1017   d->Add(fLowCuts);
1018
1019   // d->Add(lowCut);
1020   // d->Add(nXi);
1021   // d->Add(sigma);
1022   d->Add(AliForwardUtil::MakeParameter("angle", fCorrectAngles));
1023   d->Add(AliForwardUtil::MakeParameter("lowSignal",
1024                                        fZeroSharedHitsBelowThreshold));
1025   d->Add(AliForwardUtil::MakeParameter("simple", fUseSimpleMerging));
1026
1027   TObjArray* extraDead = new TObjArray;
1028   extraDead->SetOwner();
1029   extraDead->SetName("extraDead");
1030   for (Int_t i = 0; i < fExtraDead.GetSize(); i++) {
1031     if (fExtraDead.At(i) < 0) break;
1032     UShort_t dd, s, t;
1033     Char_t  r;
1034     Int_t   id = fExtraDead.At(i);
1035     AliFMDStripIndex::Unpack(id, dd, r, s, t);
1036     extraDead->Add(AliForwardUtil::MakeParameter(Form("FMD%d%c[%02d,%03d]",
1037                                                       dd, r, s, t), id));
1038   }
1039   d->Add(extraDead);
1040   fLCuts.Output(d,"lCuts");
1041   fHCuts.Output(d,"hCuts");
1042
1043   TIter    next(&fRingHistos);
1044   RingHistos* o = 0;
1045   while ((o = static_cast<RingHistos*>(next()))) {
1046     o->CreateOutputObjects(d);
1047   }
1048 }
1049 //____________________________________________________________________
1050 void
1051 AliFMDSharingFilter::Print(Option_t* /*option*/) const
1052 {
1053   //
1054   // Print information
1055   //
1056   // Parameters:
1057   //    option Not used
1058   //
1059   char ind[gROOT->GetDirLevel()+1];
1060   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1061   ind[gROOT->GetDirLevel()] = '\0';
1062   std::cout << ind << ClassName() << ": " << GetName() << '\n'
1063             << std::boolalpha
1064             << ind << " Debug:                  " << fDebug << "\n"
1065             << ind << " Use corrected angles:   " << fCorrectAngles << '\n'
1066             << ind << " Zero below threshold:   "
1067             << fZeroSharedHitsBelowThreshold << '\n'
1068             << ind << " Use simple sharing:     " << fUseSimpleMerging << '\n'
1069             << std::noboolalpha << std::endl;
1070   std::cout << ind << " Low cuts: " << std::endl;
1071   fLCuts.Print();
1072   std::cout << ind << " High cuts: " << std::endl;
1073   fHCuts.Print();
1074 }
1075
1076 //====================================================================
1077 AliFMDSharingFilter::RingHistos::RingHistos()
1078   : AliForwardUtil::RingHistos(),
1079     fBefore(0),
1080     fAfter(0),
1081     fSingle(0),
1082     fDouble(0),
1083     fTriple(0),
1084     fSinglePerStrip(0),
1085     fDistanceBefore(0),
1086     fDistanceAfter(0),
1087     fBeforeAfter(0),
1088     fNeighborsBefore(0),
1089     fNeighborsAfter(0),
1090     fSum(0),
1091     fHits(0),
1092     fNHits(0)
1093 {
1094   //
1095   // Default CTOR
1096   //
1097
1098 }
1099
1100 //____________________________________________________________________
1101 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
1102   : AliForwardUtil::RingHistos(d,r),
1103     fBefore(0),
1104     fAfter(0),
1105     fSingle(0),
1106     fDouble(0),
1107     fTriple(0),
1108     fSinglePerStrip(0),
1109     fDistanceBefore(0),
1110     fDistanceAfter(0),
1111     fBeforeAfter(0),
1112     fNeighborsBefore(0),
1113     fNeighborsAfter(0),
1114     fSum(0),
1115     fHits(0),
1116     fNHits(0)
1117 {
1118   //
1119   // Constructor
1120   //
1121   // Parameters:
1122   //    d detector
1123   //    r ring
1124   //
1125   fBefore = new TH1D("esdEloss", Form("Energy loss in %s (reconstruction)",
1126                                       GetName()), 600, 0, 15);
1127   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
1128   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
1129   fBefore->SetFillColor(Color());
1130   fBefore->SetFillStyle(3001);
1131   fBefore->SetLineColor(kBlack);
1132   fBefore->SetLineStyle(2);
1133   fBefore->SetDirectory(0);
1134
1135   fAfter  = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
1136   fAfter->SetTitle(Form("Energy loss in %s (sharing corrected)", GetName()));
1137   fAfter->SetFillColor(Color()+2);
1138   fAfter->SetLineStyle(1);
1139   fAfter->SetDirectory(0);
1140
1141   fSingle = new TH1D("singleEloss", "Energy loss (single strips)",
1142                      600, 0, 15);
1143   fSingle->SetXTitle("#Delta/#Delta_{mip}");
1144   fSingle->SetYTitle("P(#Delta/#Delta_{mip})");
1145   fSingle->SetFillColor(Color());
1146   fSingle->SetFillStyle(3001);
1147   fSingle->SetLineColor(kBlack);
1148   fSingle->SetLineStyle(2);
1149   fSingle->SetDirectory(0);
1150
1151   fDouble = static_cast<TH1D*>(fSingle->Clone("doubleEloss"));
1152   fDouble->SetTitle("Energy loss (two strips)");
1153   fDouble->SetFillColor(Color()+1);
1154   fDouble->SetDirectory(0);
1155
1156   fTriple = static_cast<TH1D*>(fSingle->Clone("tripleEloss"));
1157   fTriple->SetTitle("Energy loss (three strips)");
1158   fTriple->SetFillColor(Color()+2);
1159   fTriple->SetDirectory(0);
1160
1161   //Int_t nBinsForInner = (r == 'I' ? 32 : 16);
1162   Int_t nBinsForInner = (r == 'I' ? 512 : 256);
1163   Int_t nStrips       = (r == 'I' ? 512 : 256);
1164
1165   fSinglePerStrip = new TH2D("singlePerStrip", "SinglePerStrip",
1166                              600,0,15, nBinsForInner,0,nStrips);
1167   fSinglePerStrip->SetXTitle("#Delta/#Delta_{mip}");
1168   fSinglePerStrip->SetYTitle("Strip #");
1169   fSinglePerStrip->SetZTitle("Counts");
1170   fSinglePerStrip->SetDirectory(0);
1171
1172   fDistanceBefore = new TH1D("distanceBefore", "Distance before sharing",
1173                              nStrips , 0,nStrips );
1174   fDistanceBefore->SetXTitle("Distance");
1175   fDistanceBefore->SetYTitle("Counts");
1176   fDistanceBefore->SetFillColor(kGreen+2);
1177   fDistanceBefore->SetFillStyle(3001);
1178   fDistanceBefore->SetLineColor(kBlack);
1179   fDistanceBefore->SetLineStyle(2);
1180   fDistanceBefore->SetDirectory(0);
1181
1182   fDistanceAfter = static_cast<TH1D*>(fDistanceBefore->Clone("distanceAfter"));
1183   fDistanceAfter->SetTitle("Distance after sharing");
1184   fDistanceAfter->SetFillColor(kGreen+1);
1185   fDistanceAfter->SetDirectory(0);
1186
1187
1188   Double_t max = 15;
1189   Double_t min = -1;
1190   Int_t    n   = int((max-min) / (max / 300));
1191   fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation",
1192                           n, min, max, n, min, max);
1193   fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
1194   fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
1195   fBeforeAfter->SetZTitle("Correlation");
1196   fBeforeAfter->SetDirectory(0);
1197
1198   fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
1199   fNeighborsBefore->SetTitle("Correlation of neighbors before");
1200   fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
1201   fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
1202   fNeighborsBefore->SetDirectory(0);
1203
1204   fNeighborsAfter =
1205     static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
1206   fNeighborsAfter->SetTitle("Correlation of neighbors after");
1207   fNeighborsAfter->SetDirectory(0);
1208
1209   fSum = new TH2D("summed", "Summed signal", 200, -4, 6,
1210                   (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1211   fSum->SetDirectory(0);
1212   fSum->Sumw2();
1213   fSum->SetMarkerColor(Color());
1214   // fSum->SetFillColor(Color());
1215   fSum->SetXTitle("#eta");
1216   fSum->SetYTitle("#varphi [radians]");
1217   fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) ");
1218
1219   fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
1220   fHits->SetDirectory(0);
1221   fHits->SetXTitle("# of hits");
1222   fHits->SetFillColor(kGreen+1);
1223 }
1224 //____________________________________________________________________
1225 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
1226   : AliForwardUtil::RingHistos(o),
1227     fBefore(o.fBefore),
1228     fAfter(o.fAfter),
1229     fSingle(o.fSingle),
1230     fDouble(o.fDouble),
1231     fTriple(o.fTriple),
1232     fSinglePerStrip(o.fSinglePerStrip),
1233     fDistanceBefore(o.fDistanceBefore),
1234     fDistanceAfter(o.fDistanceAfter),
1235     fBeforeAfter(o.fBeforeAfter),
1236     fNeighborsBefore(o.fNeighborsBefore),
1237     fNeighborsAfter(o.fNeighborsAfter),
1238     fSum(o.fSum),
1239     fHits(o.fHits),
1240     fNHits(o.fNHits)
1241 {
1242   //
1243   // Copy constructor
1244   //
1245   // Parameters:
1246   //    o Object to copy from
1247   //
1248 }
1249
1250 //____________________________________________________________________
1251 AliFMDSharingFilter::RingHistos&
1252 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
1253 {
1254   //
1255   // Assignment operator
1256   //
1257   // Parameters:
1258   //    o Object to assign from
1259   //
1260   // Return:
1261   //    Reference to this
1262   //
1263   if (&o == this) return *this;
1264   AliForwardUtil::RingHistos::operator=(o);
1265   fDet = o.fDet;
1266   fRing = o.fRing;
1267
1268   if (fBefore)         delete  fBefore;
1269   if (fAfter)          delete  fAfter;
1270   if (fSingle)         delete  fSingle;
1271   if (fDouble)         delete  fDouble;
1272   if (fTriple)         delete  fTriple;
1273   if (fSinglePerStrip) delete fSinglePerStrip;
1274   if (fDistanceBefore) delete fDistanceBefore;
1275   if (fDistanceAfter)  delete fDistanceAfter;
1276   if (fHits)           delete fHits;
1277
1278
1279   fBefore          = static_cast<TH1D*>(o.fBefore->Clone());
1280   fAfter           = static_cast<TH1D*>(o.fAfter->Clone());
1281   fSingle          = static_cast<TH1D*>(o.fSingle->Clone());
1282   fDouble          = static_cast<TH1D*>(o.fDouble->Clone());
1283   fTriple          = static_cast<TH1D*>(o.fTriple->Clone());
1284   fSinglePerStrip  = static_cast<TH2D*>(o.fSinglePerStrip->Clone());
1285   fDistanceBefore  = static_cast<TH1D*>(o.fDistanceBefore->Clone());
1286   fDistanceAfter   = static_cast<TH1D*>(o.fDistanceAfter->Clone());
1287   fBeforeAfter     = static_cast<TH2D*>(o.fBeforeAfter->Clone());
1288   fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
1289   fNeighborsAfter  = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
1290   fHits            = static_cast<TH1D*>(o.fHits->Clone());
1291   fSum             = static_cast<TH2D*>(o.fSum->Clone());
1292
1293   return *this;
1294 }
1295 //____________________________________________________________________
1296 AliFMDSharingFilter::RingHistos::~RingHistos()
1297 {
1298   //
1299   // Destructor
1300   //
1301 }
1302 //____________________________________________________________________
1303 void
1304 AliFMDSharingFilter::RingHistos::Finish()
1305 {
1306   //
1307   // Finish off
1308   //
1309   //
1310   fHits->Fill(fNHits);
1311 }
1312
1313 //____________________________________________________________________
1314 void
1315 AliFMDSharingFilter::RingHistos::Terminate(const TList* dir, Int_t nEvents)
1316 {
1317   //
1318   // Scale the histograms to the total number of events
1319   //
1320   // Parameters:
1321   //    nEvents Number of events
1322   //    dir     Where the output is
1323   //
1324   TList* l = GetOutputList(dir);
1325   if (!l) return;
1326
1327 #if 0
1328   TH1D* before = static_cast<TH1D*>(l->FindObject("esdEloss"));
1329   TH1D* after  = static_cast<TH1D*>(l->FindObject("anaEloss"));
1330   if (before) before->Scale(1./nEvents);
1331   if (after)  after->Scale(1./nEvents);
1332 #endif
1333
1334   TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
1335   if (summed) summed->Scale(1./nEvents);
1336   fSum = summed;
1337 }
1338
1339 //____________________________________________________________________
1340 void
1341 AliFMDSharingFilter::RingHistos::CreateOutputObjects(TList* dir)
1342 {
1343   //
1344   // Make output
1345   //
1346   // Parameters:
1347   //    dir where to store
1348   //
1349   TList* d = DefineOutputList(dir);
1350
1351   d->Add(fBefore);
1352   d->Add(fAfter);
1353   d->Add(fSingle);
1354   d->Add(fDouble);
1355   d->Add(fTriple);
1356   d->Add(fSinglePerStrip);
1357   d->Add(fDistanceBefore);
1358   d->Add(fDistanceAfter);
1359   d->Add(fBeforeAfter);
1360   d->Add(fNeighborsBefore);
1361   d->Add(fNeighborsAfter);
1362   d->Add(fHits);
1363   d->Add(fSum);
1364
1365   // Removed to avoid doubly adding the list which destroys
1366   // the merging
1367   //dir->Add(d);
1368 }
1369
1370 //____________________________________________________________________
1371 //
1372 // EOF
1373 //