PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx

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