PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx

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