PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx

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