PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx

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