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