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     fMergingDisabled(false),
  68     fIgnoreESDForAngleCorrection(false)
  69 {
  70   //
  71   // Default Constructor - do not use
  72   //
  73   DGUARD(fDebug,1, "Default CTOR for AliFMDSharingFilter");
  74 }
  75
  76 //____________________________________________________________________
  77 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
  78   : TNamed("fmdSharingFilter", title),
  79     fRingHistos(),
  80     fCorrectAngles(kFALSE),
  81     // fSummed(0),
  82     fHighCuts(0),
  83     fLowCuts(0),
  84     // fOper(0),
  85     fDebug(0),
  86     fZeroSharedHitsBelowThreshold(false),
  87     fLCuts(),
  88     fHCuts(),
  89     fUseSimpleMerging(false),
  90     fThreeStripSharing(true),
  91     fMergingDisabled(false),
  92     fIgnoreESDForAngleCorrection(false)
  93 {
  94   //
  95   // Constructor
  96   //
  97   // Parameters:
  98   //    title Title of object  - not significant
  99   //
 100   DGUARD(fDebug,1, "Named CTOR for AliFMDSharingFilter: %s", title);
 101   fRingHistos.SetName(GetName());
 102   fRingHistos.SetOwner();
 103
 104   fRingHistos.Add(new RingHistos(1, 'I'));
 105   fRingHistos.Add(new RingHistos(2, 'I'));
 106   fRingHistos.Add(new RingHistos(2, 'O'));
 107   fRingHistos.Add(new RingHistos(3, 'I'));
 108   fRingHistos.Add(new RingHistos(3, 'O'));
 109
 110   fHCuts.Set(AliFMDMultCuts::kLandauSigmaWidth, 1);
 111   fLCuts.Set(AliFMDMultCuts::kFixed, .15);
 112
 113   // fExtraDead.Reset(-1);
 114 }
 115
 116 //____________________________________________________________________
 117 AliFMDSharingFilter::~AliFMDSharingFilter()
 118 {
 119   //
 120   // Destructor
 121   //
 122   DGUARD(fDebug,3, "DTOR for AliFMDSharingFilter");
 123   // fRingHistos.Delete();
 124 }
 125
 126 //____________________________________________________________________
 127 AliFMDSharingFilter::RingHistos*
 128 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
 129 {
 130   //
 131   // Get the ring histogram container
 132   //
 133   // Parameters:
 134   //    d Detector
 135   //    r Ring
 136   //
 137   // Return:
 138   //    Ring histogram container
 139   //
 140   Int_t idx = -1;
 141   switch (d) {
 142   case 1: idx = 0; break;
 143   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 144   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 145   }
 146   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 147
 148   return static_cast<RingHistos*>(fRingHistos.At(idx));
 149 }
 150
 151 //____________________________________________________________________
 152 void
 153 AliFMDSharingFilter::SetupForData(const TAxis& axis)
 154 {
 155   // Initialise - called on first event
 156   DGUARD(fDebug,1, "Initialize for AliFMDSharingFilter");
 157   AliForwardCorrectionManager& fcm  = AliForwardCorrectionManager::Instance();
 158   const AliFMDCorrELossFit*    fits = fcm.GetELossFit();
 159
 160   // Get the high cut.  The high cut is defined as the
 161   // most-probably-value peak found from the energy distributions, minus
 162   // 2 times the width of the corresponding Landau.
 163
 164   TAxis eAxis(axis.GetNbins(),
 165               axis.GetXmin(),
 166               axis.GetXmax());
 167   if(fits)
 168     eAxis.Set(fits->GetEtaAxis().GetNbins(),
 169               fits->GetEtaAxis().GetXmin(),
 170               fits->GetEtaAxis().GetXmax());
 171
 172   UShort_t nEta = eAxis.GetNbins();
 173
 174   fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
 175   fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
 176   fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
 177   fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
 178   fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
 179   fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
 180
 181   fLowCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
 182   fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
 183   fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
 184   fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
 185   fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
 186   fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
 187
 188   // Cache our cuts in histograms
 189   fLCuts.FillHistogram(fLowCuts);
 190   fHCuts.FillHistogram(fHighCuts);
 191 }
 192
 193 //____________________________________________________________________
 194 #define ETA2COS(ETA)                                            \
 195   TMath::Cos(2*TMath::ATan(TMath::Exp(-TMath::Abs(ETA))))
 196
 197 Bool_t
 198 AliFMDSharingFilter::Filter(const AliESDFMD& input,
 199                             Bool_t           /*lowFlux*/,
 200                             AliESDFMD&       output,
 201                             Double_t         /*zvtx*/)
 202 {
 203   //
 204   // Filter the input AliESDFMD object
 205   //
 206   // Parameters:
 207   //    input     Input
 208   //    lowFlux   If this is a low-flux event
 209   //    output    Output AliESDFMD object
 210   //
 211   // Return:
 212   //    True on success, false otherwise
 213   //
 214   DGUARD(fDebug,1, "Filter event in AliFMDSharingFilter");
 215   output.Clear();
 216   TIter    next(&fRingHistos);
 217   RingHistos* o      = 0;
 218   while ((o = static_cast<RingHistos*>(next()))) o->Clear();
 219
 220   Int_t nSingle    = 0;
 221   Int_t nDouble    = 0;
 222   Int_t nTriple    = 0;
 223
 224   for(UShort_t d = 1; d <= 3; d++) {
 225     Int_t nRings = (d == 1 ? 1 : 2);
 226     for (UShort_t q = 0; q < nRings; q++) {
 227       Char_t      r      = (q == 0 ? 'I' : 'O');
 228       UShort_t    nsec   = (q == 0 ?  20 :  40);
 229       UShort_t    nstr   = (q == 0 ? 512 : 256);
 230       RingHistos* histos = GetRingHistos(d, r);
 231
 232       for(UShort_t s = 0; s < nsec;  s++) {
 233         // `used' flags if the _current_ strip was used by _previous_
 234         // iteration.
 235         Bool_t   used            = kFALSE;
 236         // `eTotal' contains the current sum of merged signals so far
 237         Double_t eTotal          = -1;
 238         // Int_t    nDistanceBefore = -1;
 239         // Int_t    nDistanceAfter  = -1;
 240         // `twoLow' flags if we saw two consequtive strips with a
 241         // signal between the two cuts.
 242         Bool_t   twoLow          = kFALSE;
 243         Int_t    nStripsAboveCut = 0;
 244
 245         for(UShort_t t = 0; t < nstr; t++) {
 246           // nDistanceBefore++;
 247           // nDistanceAfter++;
 248
 249           output.SetMultiplicity(d,r,s,t,0.);
 250           Float_t mult         = SignalInStrip(input,d,r,s,t);
 251           Float_t multNext     = (t<nstr-1) ? SignalInStrip(input,d,r,s,t+1) :0;
 252           Float_t multNextNext = (t<nstr-2) ? SignalInStrip(input,d,r,s,t+2) :0;
 253           if (multNext     ==  AliESDFMD::kInvalidMult) multNext     = 0;
 254           if (multNextNext ==  AliESDFMD::kInvalidMult) multNextNext = 0;
 255           if(!fThreeStripSharing) multNextNext = 0;
 256
 257           // Get the pseudo-rapidity
 258           Double_t eta = input.Eta(d,r,s,t);
 259           Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
 260           if (s == 0) output.SetEta(d,r,s,t,eta);
 261
 262           // Keep dead-channel information - either from the ESD (but
 263           // see above for older data) or from the settings in the
 264           // ForwardAODConfig.C file.
 265           if (mult == AliESDFMD::kInvalidMult) {
 266             output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
 267             histos->fBefore->Fill(-1);
 268             mult = AliESDFMD::kInvalidMult;
 269           }
 270
 271           Double_t lowCut  = GetLowCut(d, r, eta);
 272           Double_t highCut = GetHighCut(d, r, eta, false);
 273           if (mult != AliESDFMD::kInvalidMult && mult > lowCut) {
 274             // Always fill the ESD sum histogram
 275             histos->fSumESD->Fill(eta, phi, mult);
 276           }
 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             // Flush a possible signal
 282             if (eTotal > 0 && t > 0)
 283               output.SetMultiplicity(d,r,s,t-1,eTotal);
 284             // Reset states so we do not try to merge over a dead strip.
 285             eTotal = -1;
 286             used   = false;
 287             twoLow = false;
 288             if (t > 0)
 289               histos->fNConsecutive->Fill(nStripsAboveCut);
 290             if (mult == AliESDFMD::kInvalidMult)
 291               // Why not fill immidiately here?
 292               nStripsAboveCut = -1;
 293             else
 294               // Why not fill immidiately here?
 295               nStripsAboveCut = 0;
 296             continue;
 297           }
 298
 299           // Fill the diagnostics histogram
 300           histos->fBefore->Fill(mult);
 301
 302           Double_t mergedEnergy = mult;
 303           // it seems to me that this logic could be condensed a bit
 304           if(mult > lowCut) {
 305             if(nStripsAboveCut < 1) {
 306               if(t > 0)
 307                 histos->fNConsecutive->Fill(nStripsAboveCut);
 308               nStripsAboveCut=0;
 309             }
 310             nStripsAboveCut++;
 311           }
 312           else {
 313             if (t > 0)
 314               histos->fNConsecutive->Fill(nStripsAboveCut);
 315             nStripsAboveCut=0;
 316           }
 317
 318           if (!fMergingDisabled) {
 319             mergedEnergy = 0;
 320
 321             // The current sum
 322             Float_t etot = 0;
 323
 324             // Fill in neighbor information
 325             if (t < nstr-1) histos->fNeighborsBefore->Fill(mult,multNext);
 326
 327             Bool_t thisValid = mult     > lowCut;
 328             Bool_t nextValid = multNext > lowCut;
 329             Bool_t thisSmall = mult     < highCut;
 330             Bool_t nextSmall = multNext < highCut;
 331
 332             // If this strips signal is above the high cut, reset distance
 333             // if (!thisSmall) {
 334             //    histos->fDistanceBefore->Fill(nDistanceBefore);
 335             //    nDistanceBefore = -1;
 336             // }
 337
 338             // If the total signal in the past 1 or 2 strips are non-zero
 339             // we need to check
 340             if (eTotal > 0) {
 341               // Here, we have already flagged one strip as a candidate
 342
 343               // If 3-strip merging is enabled, then check the next
 344               // strip to see that it falls within cut, or if we have
 345               // two low signals
 346               if (fThreeStripSharing && nextValid && (nextSmall || twoLow)) {
 347                 eTotal = eTotal + multNext;
 348                 used = kTRUE;
 349                 histos->fTriple->Fill(eTotal);
 350                 nTriple++;
 351                 twoLow = kFALSE;
 352               }
 353               // Otherwise, we got a double hit before, and that
 354               // should be stored.
 355               else {
 356                 used = kFALSE;
 357                 histos->fDouble->Fill(eTotal);
 358                 nDouble++;
 359               }
 360               // Store energy loss and reset sum
 361               etot   = eTotal;
 362               eTotal = -1;
 363             } // if (eTotal>0)
 364             else {
 365               // If we have no current sum
 366
 367               // Check if this is marked as used, and if so, continue
 368               if (used) {used = kFALSE; continue; }
 369
 370               // If the signal is abvoe the cut, set current
 371               if (thisValid) etot = mult;
 372
 373               // If the signal is abiove the cut, and so is the next
 374               // signal and either of them are below the high cut,
 375               if (thisValid  && nextValid  && (thisSmall || nextSmall)) {
 376
 377                 // If this is below the high cut, and the next is too, then
 378                 // we have two low signals
 379                 if (thisSmall && nextSmall) twoLow = kTRUE;
 380
 381                 // If this signal is bigger than the next, and the
 382                 // one after that is below the low-cut, then update
 383                 // the sum
 384                 if (mult>multNext && multNextNext < lowCut) {
 385                   etot = mult + multNext;
 386                   used = kTRUE;
 387                   histos->fDouble->Fill(etot);
 388                   nDouble++;
 389                 }
 390                 // Otherwise, we may need to merge with a third strip
 391                 else {
 392                   etot   = 0;
 393                   eTotal = mult + multNext;
 394                 }
 395               }
 396               // This is a signle hit
 397               else if(etot > 0) {
 398                 histos->fSingle->Fill(etot);
 399                 histos->fSinglePerStrip->Fill(etot,t);
 400                 nSingle++;
 401               }
 402             } // else if (etotal >= 0)
 403
 404             mergedEnergy = etot;
 405             // if (mergedEnergy > GetHighCut(d, r, eta ,false)) {
 406             //   histos->fDistanceAfter->Fill(nDistanceAfter);
 407             //   nDistanceAfter    = -1;
 408             // }
 409             //if(mult>0 && multNext >0)
 410             //  std::cout<<mult<<"  "<<multNext<<"  "<<mergedEnergy<<std::endl;
 411           } // if (!fMergingDisabled)
 412
 413           if (!fCorrectAngles)
 414             mergedEnergy = AngleCorrect(mergedEnergy, eta);
 415           // if (mergedEnergy > 0) histos->Incr();
 416
 417           if (t != 0)
 418             histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1),
 419                                           mergedEnergy);
 420           histos->fBeforeAfter->Fill(mult, mergedEnergy);
 421           if(mergedEnergy > 0)
 422             histos->fAfter->Fill(mergedEnergy);
 423           histos->fSum->Fill(eta,phi,mergedEnergy);
 424
 425           output.SetMultiplicity(d,r,s,t,mergedEnergy);
 426         } // for strip
 427         histos->fNConsecutive->Fill(nStripsAboveCut); // fill the last sector
 428       } // for sector
 429     } // for ring
 430   } // for detector
 431   DMSG(fDebug, 3,"single=%9d, double=%9d, triple=%9d",
 432        nSingle, nDouble, nTriple);
 433   next.Reset();
 434   // while ((o = static_cast<RingHistos*>(next()))) o->Finish();
 435
 436   return kTRUE;
 437 }
 438
 439 //_____________________________________________________________________
 440 Double_t
 441 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
 442                                    UShort_t         d,
 443                                    Char_t           r,
 444                                    UShort_t         s,
 445                                    UShort_t         t) const
 446 {
 447   //
 448   // Get the signal in a strip
 449   //
 450   // Parameters:
 451   //    fmd   ESD object
 452   //    d     Detector
 453   //    r     Ring
 454   //    s     Sector
 455   //    t     Strip
 456   //
 457   // Return:
 458   //    The energy signal
 459   //
 460   Double_t mult = input.Multiplicity(d,r,s,t);
 461   // In case of
 462   //  - bad value (invalid or 0)
 463   //  - we want angle corrected and data is
 464   //  - we don't want angle corrected and data isn't
 465   // just return read value
 466   if (mult == AliESDFMD::kInvalidMult               ||
 467       mult == 0                                     ||
 468       (fCorrectAngles && (fIgnoreESDForAngleCorrection || input.IsAngleCorrected())) ||
 469       (!fCorrectAngles && !fIgnoreESDForAngleCorrection && !input.IsAngleCorrected()))
 470     return mult;
 471
 472   // If we want angle corrected data, correct it,
 473   // otherwise de-correct it
 474   if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
 475   else                mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
 476   return mult;
 477 }
 478
 479 namespace {
 480   Double_t Rng2Cut(UShort_t d, Char_t r, Double_t eta, TH2* h) {
 481     Double_t ret = 1024;
 482     Int_t ybin = 0;
 483     switch(d) {
 484     case 1: ybin = 1; break;
 485     case 2: ybin = (r=='i' || r=='I') ? 2 : 3; break;
 486     case 3: ybin = (r=='i' || r=='I') ? 4 : 5; break;
 487     default: return ret;
 488     }
 489     Int_t xbin = h->GetXaxis()->FindBin(eta);
 490     if (xbin < 1 && xbin > h->GetXaxis()->GetNbins()) return ret;
 491     ret = h->GetBinContent(xbin,ybin);
 492     return ret;
 493   }
 494 }
 495
 496   //_____________________________________________________________________
 497 Double_t
 498 AliFMDSharingFilter::GetLowCut(UShort_t d, Char_t r, Double_t eta) const
 499 {
 500   //
 501   // Get the low cut.  Normally, the low cut is taken to be the lower
 502   // value of the fit range used when generating the energy loss fits.
 503   // However, if fLowCut is set (using SetLowCit) to a value greater
 504   // than 0, then that value is used.
 505   //
 506   return Rng2Cut(d, r, eta, fLowCuts);
 507   // return fLCuts.GetMultCut(d,r,eta,false);
 508 }
 509
 510 //_____________________________________________________________________
 511 Double_t
 512 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r,
 513                                 Double_t eta, Bool_t /*errors*/) const
 514 {
 515   //
 516   // Get the high cut.  The high cut is defined as the
 517   // most-probably-value peak found from the energy distributions, minus
 518   // 2 times the width of the corresponding Landau.
 519   //
 520   return Rng2Cut(d, r, eta, fHighCuts);
 521   // return fHCuts.GetMultCut(d,r,eta,errors);
 522 }
 523
 524 //____________________________________________________________________
 525 Double_t
 526 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
 527 {
 528   //
 529   // Angle correct the signal
 530   //
 531   // Parameters:
 532   //    mult Angle Un-corrected Signal
 533   //    eta  Pseudo-rapidity
 534   //
 535   // Return:
 536   //    Angle corrected signal
 537   //
 538   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 539   if (eta < 0) theta -= TMath::Pi();
 540   return mult * TMath::Cos(theta);
 541 }
 542 //____________________________________________________________________
 543 Double_t
 544 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
 545 {
 546   //
 547   // Angle de-correct the signal
 548   //
 549   // Parameters:
 550   //    mult Angle corrected Signal
 551   //    eta  Pseudo-rapidity
 552   //
 553   // Return:
 554   //    Angle un-corrected signal
 555   //
 556   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 557   if (eta < 0) theta -= TMath::Pi();
 558   return mult / TMath::Cos(theta);
 559 }
 560
 561 //____________________________________________________________________
 562 void
 563 AliFMDSharingFilter::Terminate(const TList* dir, TList* output, Int_t nEvents)
 564 {
 565   //
 566   // Scale the histograms to the total number of events
 567   //
 568   // Parameters:
 569   //    dir     Where the output is
 570   //    nEvents Number of events
 571   //
 572   DGUARD(fDebug,1, "Scale histograms in AliFMDSharingFilter");
 573   if (nEvents <= 0) return;
 574   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 575   if (!d) return;
 576
 577   TList* out = new TList;
 578   out->SetName(d->GetName());
 579   out->SetOwner();
 580
 581   TParameter<int>* nFiles =
 582     static_cast<TParameter<int>*>(d->FindObject("nFiles"));
 583
 584   TH2* lowCuts  = static_cast<TH2*>(d->FindObject("lowCuts"));
 585   TH2* highCuts = static_cast<TH2*>(d->FindObject("highCuts"));
 586   if (lowCuts && nFiles) {
 587     lowCuts->Scale(1. / nFiles->GetVal());
 588     out->Add(lowCuts->Clone());
 589   }
 590   else
 591     AliWarning("low cuts histogram not found in input list");
 592   if (highCuts && nFiles) {
 593     highCuts->Scale(1. / nFiles->GetVal());
 594     out->Add(highCuts->Clone());
 595   }
 596   else
 597     AliWarning("high cuts histogram not found in input list");
 598
 599   TIter    next(&fRingHistos);
 600   RingHistos* o = 0;
 601   THStack* sums = new THStack("sums", "Sum of ring signals");
 602   THStack* sumsESD = new THStack("sumsESD", "Sum of ring ESD signals");
 603   while ((o = static_cast<RingHistos*>(next()))) {
 604     o->Terminate(d, nEvents);
 605     if (!o->fSum) {
 606       Warning("Terminate", "No sum histogram found for ring %s", o->GetName());
 607       continue;
 608     }
 609     TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
 610     sum->Scale(1., "width");
 611     sum->SetTitle(o->GetName());
 612     sum->SetDirectory(0);
 613     sum->SetYTitle("#sum_{c} #Delta/#Delta_{mip}");
 614     sums->Add(sum);
 615
 616     sum = o->fSumESD->ProjectionX(o->GetName(), 1, o->fSumESD->GetNbinsY(),"e");
 617     sum->Scale(1., "width");
 618     sum->SetTitle(o->GetName());
 619     sum->SetDirectory(0);
 620     sum->SetYTitle("#sum_{s} #Delta/#Delta_{mip}");
 621     sumsESD->Add(sum);
 622   }
 623   out->Add(sums);
 624   out->Add(sumsESD);
 625   output->Add(out);
 626 }
 627
 628 //____________________________________________________________________
 629 void
 630 AliFMDSharingFilter::CreateOutputObjects(TList* dir)
 631 {
 632   //
 633   // Define the output histograms.  These are put in a sub list of the
 634   // passed list.   The histograms are merged before the parent task calls
 635   // AliAnalysisTaskSE::Terminate
 636   //
 637   // Parameters:
 638   //    dir Directory to add to
 639   //
 640   DGUARD(fDebug,1, "Define output in AliFMDSharingFilter");
 641   TList* d = new TList;
 642   d->SetName(GetName());
 643   dir->Add(d);
 644
 645 #if 0
 646   fSummed = new TH2I("operations", "Strip operations",
 647                      kMergedInto, kNone-.5, kMergedInto+.5,
 648                      51201, -.5, 51200.5);
 649   fSummed->SetXTitle("Operation");
 650   fSummed->SetYTitle("# of strips");
 651   fSummed->SetZTitle("Events");
 652   fSummed->GetXaxis()->SetBinLabel(kNone,            "None");
 653   fSummed->GetXaxis()->SetBinLabel(kCandidate,       "Candidate");
 654   fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
 655   fSummed->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
 656   fSummed->SetDirectory(0);
 657   d->Add(fSummed);
 658 #endif
 659
 660   fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
 661   fHighCuts->SetXTitle("#eta");
 662   fHighCuts->SetDirectory(0);
 663   d->Add(fHighCuts);
 664
 665   fLowCuts = new TH2D("lowCuts", "Low cuts used", 1,0,1, 1,0,1);
 666   fLowCuts->SetXTitle("#eta");
 667   fLowCuts->SetDirectory(0);
 668   d->Add(fLowCuts);
 669
 670   // d->Add(lowCut);
 671   // d->Add(nXi);
 672   // d->Add(sigma);
 673   d->Add(AliForwardUtil::MakeParameter("angle", fCorrectAngles));
 674   d->Add(AliForwardUtil::MakeParameter("lowSignal",
 675                                        fZeroSharedHitsBelowThreshold));
 676   d->Add(AliForwardUtil::MakeParameter("simple", fUseSimpleMerging));
 677   d->Add(AliForwardUtil::MakeParameter("sumThree", fThreeStripSharing));
 678   d->Add(AliForwardUtil::MakeParameter("disabled", fMergingDisabled));
 679   TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
 680   nFiles->SetMergeMode('+');
 681   d->Add(nFiles);
 682
 683   fLCuts.Output(d,"lCuts");
 684   fHCuts.Output(d,"hCuts");
 685
 686   TIter    next(&fRingHistos);
 687   RingHistos* o = 0;
 688   while ((o = static_cast<RingHistos*>(next()))) {
 689     o->CreateOutputObjects(d);
 690   }
 691 }
 692 #define PF(N,V,...)                                     \
 693   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
 694 #define PFB(N,FLAG)                             \
 695   do {                                                                  \
 696     AliForwardUtil::PrintName(N);                                       \
 697     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 698   } while(false)
 699 #define PFV(N,VALUE)                                    \
 700   do {                                                  \
 701     AliForwardUtil::PrintName(N);                       \
 702     std::cout << (VALUE) << std::endl; } while(false)
 703
 704 //____________________________________________________________________
 705 void
 706 AliFMDSharingFilter::Print(Option_t* /*option*/) const
 707 {
 708   //
 709   // Print information
 710   //
 711   // Parameters:
 712   //    option Not used
 713   //
 714   AliForwardUtil::PrintTask(*this);
 715   gROOT->IncreaseDirLevel();
 716
 717   PFB("Use corrected angles",  fCorrectAngles);
 718   PFB("Zero below threshold",  fZeroSharedHitsBelowThreshold);
 719   PFB("Use simple sharing",    fUseSimpleMerging);
 720   PFB("Allow 3 strip merging", fThreeStripSharing);
 721   PF("Low cuts",        "");
 722   fLCuts.Print();
 723   PF("High cuts",       "");
 724   fHCuts.Print();
 725   gROOT->DecreaseDirLevel();
 726 }
 727
 728 //====================================================================
 729 AliFMDSharingFilter::RingHistos::RingHistos()
 730   : AliForwardUtil::RingHistos(),
 731     fBefore(0),
 732     fAfter(0),
 733     fSingle(0),
 734     fDouble(0),
 735     fTriple(0),
 736     fSinglePerStrip(0),
 737     // fDistanceBefore(0),
 738     // fDistanceAfter(0),
 739     fBeforeAfter(0),
 740     fNeighborsBefore(0),
 741     fNeighborsAfter(0),
 742     fSumESD(0),
 743     fSum(0),
 744     fNConsecutive(0)
 745      // ,
 746     // fHits(0),
 747     // fNHits(0)
 748 {
 749   //
 750   // Default CTOR
 751   //
 752
 753 }
 754
 755 //____________________________________________________________________
 756 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
 757   : AliForwardUtil::RingHistos(d,r),
 758     fBefore(0),
 759     fAfter(0),
 760     fSingle(0),
 761     fDouble(0),
 762     fTriple(0),
 763     fSinglePerStrip(0),
 764     // fDistanceBefore(0),
 765     // fDistanceAfter(0),
 766     fBeforeAfter(0),
 767     fNeighborsBefore(0),
 768     fNeighborsAfter(0),
 769     fSumESD(0),
 770     fSum(0),
 771     fNConsecutive(0)
 772      //,
 773     // fHits(0),
 774     // fNHits(0)
 775 {
 776   //
 777   // Constructor
 778   //
 779   // Parameters:
 780   //    d detector
 781   //    r ring
 782   //
 783   fBefore = new TH1D("esdEloss", Form("Energy loss in %s (reconstruction)",
 784                                       GetName()), 640, -1, 15);
 785   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
 786   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
 787   fBefore->SetFillColor(Color());
 788   fBefore->SetFillStyle(3001);
 789   fBefore->SetLineColor(kBlack);
 790   fBefore->SetLineStyle(2);
 791   fBefore->SetDirectory(0);
 792
 793   fAfter  = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
 794   fAfter->SetTitle(Form("Energy loss in %s (sharing corrected)", GetName()));
 795   fAfter->SetFillColor(Color()+2);
 796   fAfter->SetLineStyle(1);
 797   fAfter->SetDirectory(0);
 798
 799   fSingle = new TH1D("singleEloss", "Energy loss (single strips)",
 800                      600, 0, 15);
 801   fSingle->SetXTitle("#Delta/#Delta_{mip}");
 802   fSingle->SetYTitle("P(#Delta/#Delta_{mip})");
 803   fSingle->SetFillColor(Color());
 804   fSingle->SetFillStyle(3001);
 805   fSingle->SetLineColor(kBlack);
 806   fSingle->SetLineStyle(2);
 807   fSingle->SetDirectory(0);
 808
 809   fDouble = static_cast<TH1D*>(fSingle->Clone("doubleEloss"));
 810   fDouble->SetTitle("Energy loss (two strips)");
 811   fDouble->SetFillColor(Color()+1);
 812   fDouble->SetDirectory(0);
 813
 814   fTriple = static_cast<TH1D*>(fSingle->Clone("tripleEloss"));
 815   fTriple->SetTitle("Energy loss (three strips)");
 816   fTriple->SetFillColor(Color()+2);
 817   fTriple->SetDirectory(0);
 818
 819   //Int_t nBinsForInner = (r == 'I' ? 32 : 16);
 820   Int_t nBinsForInner = (r == 'I' ? 512 : 256);
 821   Int_t nStrips       = (r == 'I' ? 512 : 256);
 822
 823   fSinglePerStrip = new TH2D("singlePerStrip", "SinglePerStrip",
 824                              600,0,15, nBinsForInner,0,nStrips);
 825   fSinglePerStrip->SetXTitle("#Delta/#Delta_{mip}");
 826   fSinglePerStrip->SetYTitle("Strip #");
 827   fSinglePerStrip->SetZTitle("Counts");
 828   fSinglePerStrip->SetDirectory(0);
 829
 830 #if 0
 831   fDistanceBefore = new TH1D("distanceBefore", "Distance before sharing",
 832                              nStrips , 0,nStrips );
 833   fDistanceBefore->SetXTitle("Distance");
 834   fDistanceBefore->SetYTitle("Counts");
 835   fDistanceBefore->SetFillColor(kGreen+2);
 836   fDistanceBefore->SetFillStyle(3001);
 837   fDistanceBefore->SetLineColor(kBlack);
 838   fDistanceBefore->SetLineStyle(2);
 839   fDistanceBefore->SetDirectory(0);
 840
 841   fDistanceAfter = static_cast<TH1D*>(fDistanceBefore->Clone("distanceAfter"));
 842   fDistanceAfter->SetTitle("Distance after sharing");
 843   fDistanceAfter->SetFillColor(kGreen+1);
 844   fDistanceAfter->SetDirectory(0);
 845 #endif
 846
 847   Double_t max = 15;
 848   Double_t min = -1;
 849   Int_t    n   = int((max-min) / (max / 300));
 850   fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation",
 851                           n, min, max, n, min, max);
 852   fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
 853   fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
 854   fBeforeAfter->SetZTitle("Correlation");
 855   fBeforeAfter->SetDirectory(0);
 856
 857   fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
 858   fNeighborsBefore->SetTitle("Correlation of neighbors before");
 859   fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
 860   fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
 861   fNeighborsBefore->SetDirectory(0);
 862
 863   fNeighborsAfter =
 864     static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
 865   fNeighborsAfter->SetTitle("Correlation of neighbors after");
 866   fNeighborsAfter->SetDirectory(0);
 867
 868   fSumESD = new TH2D("summedESD", "Summed ESD signal", 200, -4, 6,
 869                   NSector(), 0, 2*TMath::Pi());
 870   fSumESD->SetDirectory(0);
 871   fSumESD->Sumw2();
 872   fSumESD->SetMarkerColor(Color());
 873   // fSum->SetFillColor(Color());
 874   fSumESD->SetXTitle("#eta");
 875   fSumESD->SetYTitle("#varphi [radians]");
 876   fSumESD->SetZTitle("#sum_{strip} #Delta/#Delta_{mip}(#eta,#varphi) ");
 877
 878   fSum = static_cast<TH2D*>(fSumESD->Clone("summed"));
 879   fSum->SetTitle("Summed cluster signal");
 880   fSum->SetZTitle("#sum_{cluster} #Delta/#Delta_{mip}(#eta,#varphi) ");
 881   fSum->SetDirectory(0);
 882
 883   // Perhaps we need to ensure that this histogram has enough range to
 884   // accommondate all possible ranges - that is, from -1 to the number
 885   // of strips in this ring(-type) - i.e., NStrips().  Perhaps the
 886   // axis should be defined with increasin bin size - e.g.,
 887   //
 888   //   -1.5,-.5,.5,1.5,...,100.5,128.5,192.5,...,NStrips()
 889   //
 890   fNConsecutive = new TH1D("nConsecutive","# consecutive strips above low cut",
 891                            201,-1.5,199.5);
 892   fNConsecutive->SetXTitle("N_{strips}");
 893   fNConsecutive->SetYTitle("N_{entries}");
 894   fNConsecutive->SetFillColor(kYellow+2);
 895   fNConsecutive->SetFillStyle(3001);
 896   fNConsecutive->SetDirectory(0);
 897
 898
 899 #if 0
 900   fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
 901   fHits->SetDirectory(0);
 902   fHits->SetXTitle("# of hits");
 903   fHits->SetFillColor(kGreen+1);
 904 #endif
 905 }
 906 //____________________________________________________________________
 907 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
 908   : AliForwardUtil::RingHistos(o),
 909     fBefore(o.fBefore),
 910     fAfter(o.fAfter),
 911     fSingle(o.fSingle),
 912     fDouble(o.fDouble),
 913     fTriple(o.fTriple),
 914     fSinglePerStrip(o.fSinglePerStrip),
 915     // fDistanceBefore(o.fDistanceBefore),
 916     // fDistanceAfter(o.fDistanceAfter),
 917     fBeforeAfter(o.fBeforeAfter),
 918     fNeighborsBefore(o.fNeighborsBefore),
 919     fNeighborsAfter(o.fNeighborsAfter),
 920     fSumESD(o.fSumESD), //,
 921     fSum(o.fSum),
 922     fNConsecutive(o.fNConsecutive)
 923      //,
 924     // fHits(o.fHits),
 925     // fNHits(o.fNHits)
 926 {
 927   //
 928   // Copy constructor
 929   //
 930   // Parameters:
 931   //    o Object to copy from
 932   //
 933 }
 934
 935 //____________________________________________________________________
 936 AliFMDSharingFilter::RingHistos&
 937 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
 938 {
 939   //
 940   // Assignment operator
 941   //
 942   // Parameters:
 943   //    o Object to assign from
 944   //
 945   // Return:
 946   //    Reference to this
 947   //
 948   if (&o == this) return *this;
 949   AliForwardUtil::RingHistos::operator=(o);
 950   fDet = o.fDet;
 951   fRing = o.fRing;
 952
 953   if (fBefore)         delete  fBefore;
 954   if (fAfter)          delete  fAfter;
 955   if (fSingle)         delete  fSingle;
 956   if (fDouble)         delete  fDouble;
 957   if (fTriple)         delete  fTriple;
 958   if (fSinglePerStrip) delete  fSinglePerStrip;
 959   if (fNConsecutive)   delete  fNConsecutive;
 960   // if (fDistanceBefore) delete fDistanceBefore;
 961   // if (fDistanceAfter)  delete fDistanceAfter;
 962   // if (fHits)                delete fHits;
 963
 964
 965   fBefore          = static_cast<TH1D*>(o.fBefore->Clone());
 966   fAfter           = static_cast<TH1D*>(o.fAfter->Clone());
 967   fSingle          = static_cast<TH1D*>(o.fSingle->Clone());
 968   fDouble          = static_cast<TH1D*>(o.fDouble->Clone());
 969   fTriple          = static_cast<TH1D*>(o.fTriple->Clone());
 970   fSinglePerStrip  = static_cast<TH2D*>(o.fSinglePerStrip->Clone());
 971   // fDistanceBefore  = static_cast<TH1D*>(o.fDistanceBefore->Clone());
 972   // fDistanceAfter   = static_cast<TH1D*>(o.fDistanceAfter->Clone());
 973   fBeforeAfter     = static_cast<TH2D*>(o.fBeforeAfter->Clone());
 974   fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
 975   fNeighborsAfter  = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
 976   // fHits            = static_cast<TH1D*>(o.fHits->Clone());
 977   fSumESD          = static_cast<TH2D*>(o.fSumESD->Clone());
 978   fSum             = static_cast<TH2D*>(o.fSum->Clone());
 979   fNConsecutive    = static_cast<TH1D*>(o.fNConsecutive->Clone());
 980
 981   return *this;
 982 }
 983 //____________________________________________________________________
 984 AliFMDSharingFilter::RingHistos::~RingHistos()
 985 {
 986   //
 987   // Destructor
 988   //
 989 }
 990 #if 0
 991 //____________________________________________________________________
 992 void
 993 AliFMDSharingFilter::RingHistos::Finish()
 994 {
 995   //
 996   // Finish off
 997   //
 998   //
 999   // fHits->Fill(fNHits);
1000 }
1001 #endif
1002 //____________________________________________________________________
1003 void
1004 AliFMDSharingFilter::RingHistos::Terminate(const TList* dir, Int_t nEvents)
1005 {
1006   //
1007   // Scale the histograms to the total number of events
1008   //
1009   // Parameters:
1010   //    nEvents Number of events
1011   //    dir     Where the output is
1012   //
1013   TList* l = GetOutputList(dir);
1014   if (!l) return;
1015
1016 #if 0
1017   TH1D* before = static_cast<TH1D*>(l->FindObject("esdEloss"));
1018   TH1D* after  = static_cast<TH1D*>(l->FindObject("anaEloss"));
1019   if (before) before->Scale(1./nEvents);
1020   if (after)  after->Scale(1./nEvents);
1021 #endif
1022
1023   TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
1024   if (summed) summed->Scale(1./nEvents);
1025   fSum = summed;
1026
1027   TH2D* summedESD = static_cast<TH2D*>(l->FindObject("summedESD"));
1028   if (summedESD) summedESD->Scale(1./nEvents);
1029   fSumESD = summedESD;
1030
1031   TH1D* consecutive = static_cast<TH1D*>(l->FindObject("nConsecutive"));
1032   if (consecutive) consecutive->Scale(1./nEvents);
1033   fNConsecutive= consecutive;
1034 }
1035
1036 //____________________________________________________________________
1037 void
1038 AliFMDSharingFilter::RingHistos::CreateOutputObjects(TList* dir)
1039 {
1040   //
1041   // Make output
1042   //
1043   // Parameters:
1044   //    dir where to store
1045   //
1046   TList* d = DefineOutputList(dir);
1047
1048   d->Add(fBefore);
1049   d->Add(fAfter);
1050   d->Add(fSingle);
1051   d->Add(fDouble);
1052   d->Add(fTriple);
1053   d->Add(fSinglePerStrip);
1054   // d->Add(fDistanceBefore);
1055   // d->Add(fDistanceAfter);
1056   d->Add(fBeforeAfter);
1057   d->Add(fNeighborsBefore);
1058   d->Add(fNeighborsAfter);
1059   // d->Add(fHits);
1060   d->Add(fSumESD);
1061   d->Add(fSum);
1062   d->Add(fNConsecutive);
1063
1064   // Removed to avoid doubly adding the list which destroys
1065   // the merging
1066   //dir->Add(d);
1067 }
1068
1069 //____________________________________________________________________
1070 //
1071 // EOF
1072 //