PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx

   1 //
   2 // This class collects the event histograms into single histograms,
   3 // one for each ring in each vertex bin.
   4 //
   5 // Input:
   6 //   - AliESDFMD object possibly corrected for sharing
   7 //
   8 // Output:
   9 //   - 5 RingHistos objects - each with a number of vertex dependent
  10 //     2D histograms of the inclusive charge particle density
  11 //
  12 // HistCollector used:
  13 //   - AliFMDCorrSecondaryMap
  14 //
  15 #include "AliFMDHistCollector.h"
  16 #include <AliESDFMD.h>
  17 #include <TAxis.h>
  18 #include <TList.h>
  19 #include <TMath.h>
  20 #include "AliForwardCorrectionManager.h"
  21 #include "AliLog.h"
  22 #include <TH2D.h>
  23 #include <TH1I.h>
  24 #include <TProfile.h>
  25 #include <TProfile2D.h>
  26 #include <TArrayI.h>
  27 #include <TROOT.h>
  28 #include <iostream>
  29 #include <iomanip>
  30
  31 ClassImp(AliFMDHistCollector)
  32 #if 0
  33 ; // For Emacs
  34 #endif
  35
  36 //____________________________________________________________________
  37 AliFMDHistCollector::AliFMDHistCollector()
  38   : fNCutBins(0),
  39     fCorrectionCut(0),
  40     fFirstBins(),
  41     fLastBins(),
  42     fDebug(0),
  43     fList(0),
  44     fSumRings(0),
  45     fCoverage(0),
  46     fMergeMethod(kStraightMean),
  47     fFiducialMethod(kByCut),
  48     fSkipFMDRings(0),
  49     fBgAndHitMaps(false)
  50 {}
  51
  52 //____________________________________________________________________
  53 AliFMDHistCollector::AliFMDHistCollector(const char* title)
  54   : TNamed("fmdHistCollector", title),
  55     fNCutBins(2),
  56     fCorrectionCut(0.5),
  57     fFirstBins(1),
  58     fLastBins(1),
  59     fDebug(0),
  60     fList(0),
  61     fSumRings(0),
  62     fCoverage(0),
  63     fMergeMethod(kStraightMean),
  64     fFiducialMethod(kByCut),
  65     fSkipFMDRings(0),
  66     fBgAndHitMaps(false)
  67 {
  68 }
  69 //____________________________________________________________________
  70 AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
  71   : TNamed(o),
  72     fNCutBins(o.fNCutBins),
  73     fCorrectionCut(o.fCorrectionCut),
  74     fFirstBins(o.fFirstBins),
  75     fLastBins(o.fLastBins),
  76     fDebug(o.fDebug),
  77     fList(o.fList),
  78     fSumRings(o.fSumRings),
  79     fCoverage(o.fCoverage),
  80     fMergeMethod(o.fMergeMethod),
  81     fFiducialMethod(o.fFiducialMethod),
  82     fSkipFMDRings(o.fSkipFMDRings),
  83     fBgAndHitMaps(o.fBgAndHitMaps)
  84 {}
  85
  86 //____________________________________________________________________
  87 AliFMDHistCollector::~AliFMDHistCollector()
  88 {
  89   if (fList) delete fList;
  90 }
  91 //____________________________________________________________________
  92 AliFMDHistCollector&
  93 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
  94 {
  95   //
  96   // Assignement operator
  97   //
  98   // Parameters:
  99   //    o Object to assign from
 100   //
 101   // Return:
 102   //    Reference to this object
 103   //
 104   if (&o == this) return *this;
 105   TNamed::operator=(o);
 106
 107   fNCutBins       = o.fNCutBins;
 108   fCorrectionCut  = o.fCorrectionCut;
 109   fFirstBins      = o.fFirstBins;
 110   fLastBins       = o.fLastBins;
 111   fDebug          = o.fDebug;
 112   fList           = o.fList;
 113   fSumRings       = o.fSumRings;
 114   fCoverage       = o.fCoverage;
 115   fMergeMethod    = o.fMergeMethod;
 116   fFiducialMethod = o.fFiducialMethod;
 117   fSkipFMDRings   = o.fSkipFMDRings;
 118   fBgAndHitMaps   = o.fBgAndHitMaps;
 119
 120   return *this;
 121 }
 122
 123 //____________________________________________________________________
 124 void
 125 AliFMDHistCollector::Init(const TAxis& vtxAxis,
 126                           const TAxis& etaAxis)
 127 {
 128   //
 129   // Intialise
 130   //
 131   // Parameters:
 132   //    vtxAxis  Vertex axis
 133   //
 134
 135   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 136
 137   fSumRings = new TH2D("sumRings", "Sum in individual rings",
 138                        etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
 139                        5, 1, 6);
 140   fSumRings->Sumw2();
 141   fSumRings->SetDirectory(0);
 142   fSumRings->SetXTitle("#eta");
 143   fSumRings->GetYaxis()->SetBinLabel(1,"FMD1i");
 144   fSumRings->GetYaxis()->SetBinLabel(2,"FMD2i");
 145   fSumRings->GetYaxis()->SetBinLabel(3,"FMD2o");
 146   fSumRings->GetYaxis()->SetBinLabel(4,"FMD3i");
 147   fSumRings->GetYaxis()->SetBinLabel(5,"FMD3o");
 148   fList->Add(fSumRings);
 149
 150   fCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
 151                        etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
 152                        vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
 153   fCoverage->SetDirectory(0);
 154   fCoverage->SetXTitle("#eta");
 155   fCoverage->SetYTitle("v_{z} [cm]");
 156   fCoverage->SetZTitle("n_{bins}");
 157   fList->Add(fCoverage);
 158
 159   UShort_t nVz = vtxAxis.GetNbins();
 160   fFirstBins.Set(5*nVz);
 161   fLastBins.Set(5*nVz);
 162
 163   // Find the eta bin ranges
 164   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 165
 166     Double_t vMin    = vtxAxis.GetBinLowEdge(iVz);
 167     Double_t vMax    = vtxAxis.GetBinUpEdge(iVz);
 168     TList*   vtxList=0;
 169     if(fBgAndHitMaps) {
 170       vtxList = new TList;
 171       vtxList->SetName(Form("%c%02d_%c%02d",
 172                             vMin < 0 ? 'm' : 'p', int(TMath::Abs(vMin)),
 173                             vMax < 0 ? 'm' : 'p', int(TMath::Abs(vMax))));
 174       fList->Add(vtxList);
 175     }
 176     // Find the first and last eta bin to use for each ring for
 177     // each vertex bin.   This is instead of using the methods
 178     // provided by AliFMDAnaParameters
 179     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 180       UShort_t d = 0;
 181       Char_t   r = 0;
 182       GetDetRing(iIdx, d, r);
 183
 184       // Skipping selected FMD rings
 185       if(d==1 && r=='I' && (fSkipFMDRings & kFMD1I)) continue;
 186       if(d==2 && r=='I' && (fSkipFMDRings & kFMD2I)) continue;
 187       if(d==2 && r=='O' && (fSkipFMDRings & kFMD2O)) continue;
 188       if(d==3 && r=='I' && (fSkipFMDRings & kFMD3I)) continue;
 189       if(d==3 && r=='O' && (fSkipFMDRings & kFMD3O)) continue;
 190
 191       // Get the background object
 192       // TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
 193       TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
 194       Int_t nEta  = bg->GetNbinsX();
 195       Int_t first = nEta+1;
 196       Int_t last  = 0;
 197
 198       // Loop over the eta bins
 199       for (Int_t ie = 1; ie <= nEta; ie++) {
 200         // Loop over the phi bins to make sure that we
 201         // do not have holes in the coverage
 202         bool ok = true;
 203         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 204           if (!CheckCorrection(bg, ie, ip)) {
 205             ok = false;
 206             continue;
 207           }
 208         }
 209         if (!ok) continue;
 210
 211         first = TMath::Min(ie, first);
 212         last  = TMath::Max(ie, last);
 213       }
 214
 215       // Store the result for later use
 216       fFirstBins[(iVz-1)*5+iIdx] = first;
 217       fLastBins[(iVz-1)*5+iIdx]  = last;
 218       TH2D* obg=0;
 219       if(fBgAndHitMaps) {
 220         obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
 221         obg->SetDirectory(0);
 222         obg->Reset();
 223         vtxList->Add(obg);
 224
 225         TH2D* hitmap = static_cast<TH2D*>(bg->Clone(Form("hitMapFMD%d%c", d, r)));
 226         if(r == 'O') hitmap->RebinY(2);
 227         hitmap->SetDirectory(0);
 228         hitmap->GetZaxis()->SetTitle("");
 229         hitmap->Reset();
 230         vtxList->Add(hitmap);
 231       }
 232       // Fill diagnostics histograms
 233       for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
 234         Double_t old = fCoverage->GetBinContent(ie, iVz);
 235         fCoverage->SetBinContent(ie, iVz, old+1);
 236         if(fBgAndHitMaps) {
 237           for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 238             obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
 239             obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
 240           }
 241         }
 242       }
 243     } // for j
 244   }
 245 }
 246
 247 //____________________________________________________________________
 248 Bool_t
 249 AliFMDHistCollector::CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const
 250 {
 251   //
 252   // Check if we should include the bin in the data range
 253   //
 254   // Parameters:
 255   //    bg Secondary map histogram
 256   //    ie Eta bin
 257   //    ip Phi bin
 258   //
 259   // Return:
 260   //    True if to be used
 261   //
 262   Double_t c = bg->GetBinContent(ie,ip);
 263   switch (fFiducialMethod) {
 264   case kByCut:
 265     return c >= fCorrectionCut;
 266   case kDistance:
 267     if (2 * c < bg->GetBinContent(ie+1,ip) ||
 268         2 * c < bg->GetBinContent(ie-1,ip)) return false;
 269     return true;
 270   default:
 271     AliError("No fiducal cut method defined");
 272   }
 273   return false;
 274 }
 275
 276 //____________________________________________________________________
 277 void
 278 AliFMDHistCollector::DefineOutput(TList* dir)
 279 {
 280   //
 281   // Output diagnostic histograms to directory
 282   //
 283   // Parameters:
 284   //    dir List to write in
 285   //
 286   fList = new TList;
 287   fList->SetOwner();
 288   fList->SetName(GetName());
 289   dir->Add(fList);
 290
 291 }
 292
 293
 294 //____________________________________________________________________
 295 Int_t
 296 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
 297 {
 298   //
 299   // Get the ring index from detector number and ring identifier
 300   //
 301   // Parameters:
 302   //    d Detector
 303   //    r Ring identifier
 304   //
 305   // Return:
 306   //    ring index or -1 in case of problems
 307   //
 308   Int_t idx = -1;
 309   switch (d) {
 310   case 1: idx = 0; break;
 311   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 312   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 313   }
 314   return idx;
 315 }
 316 //____________________________________________________________________
 317 void
 318 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
 319 {
 320   //
 321   // Get the detector and ring from the ring index
 322   //
 323   // Parameters:
 324   //    idx Ring index
 325   //    d   On return, the detector or 0 in case of errors
 326   //    r   On return, the ring id or '0' in case of errors
 327   //
 328   d = 0;
 329   r = '\0';
 330   switch (idx) {
 331   case 0: d = 1; r = 'I'; break;
 332   case 1: d = 2; r = 'I'; break;
 333   case 2: d = 2; r = 'O'; break;
 334   case 3: d = 3; r = 'I'; break;
 335   case 4: d = 3; r = 'O'; break;
 336   }
 337 }
 338
 339 //____________________________________________________________________
 340 void
 341 AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
 342                                      Int_t& first, Int_t& last) const
 343 {
 344   //
 345   // Get the first and last eta bin to use for a given ring and vertex
 346   //
 347   // Parameters:
 348   //    idx      Ring index as given by GetIdx
 349   //    vtxBin   Vertex bin (1 based)
 350   //    first    On return, the first eta bin to use
 351   //    last     On return, the last eta bin to use
 352   //
 353   first = 0;
 354   last  = 0;
 355
 356   if (idx    <  0) return;
 357   if (vtxbin <= 0) return;
 358   idx += (vtxbin-1) * 5;
 359
 360   if (idx < 0 || idx >= fFirstBins.GetSize()) return;
 361
 362   first = fFirstBins.At(idx)+fNCutBins;
 363   last  = fLastBins.At(idx)-fNCutBins;
 364 }
 365
 366 //____________________________________________________________________
 367 Int_t
 368 AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
 369 {
 370   //
 371   // Get the first eta bin to use for a given ring and vertex
 372   //
 373   // Parameters:
 374   //    idx Ring index as given by GetIdx
 375   //    v vertex bin (1 based)
 376   //
 377   // Return:
 378   //    First eta bin to use, or -1 in case of problems
 379   //
 380   Int_t f, l;
 381   GetFirstAndLast(idx,v,f,l);
 382   return f;
 383 }
 384
 385
 386 //____________________________________________________________________
 387 Int_t
 388 AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
 389 {
 390   //
 391   // Get the last eta bin to use for a given ring and vertex
 392   //
 393   // Parameters:
 394   //    idx Ring index as given by GetIdx
 395   //    v vertex bin (1 based)
 396   //
 397   // Return:
 398   //    Last eta bin to use, or -1 in case of problems
 399   //
 400   Int_t f, l;
 401   GetFirstAndLast(idx,v,f,l);
 402   return l;
 403 }
 404
 405 //____________________________________________________________________
 406 Int_t
 407 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
 408                                 Int_t bin,  UShort_t vtxbin) const
 409 {
 410   //
 411   // Get the possibly overlapping histogram of eta bin @a e in
 412   // detector and ring
 413   //
 414   // Parameters:
 415   //    d Detector
 416   //    r Ring
 417   //    e Eta bin
 418   //    v Vertex bin (1 based)
 419   //
 420   // Return:
 421   //    Overlapping histogram index or -1
 422   //
 423
 424   Int_t other = -1;
 425   if (d == 1) {
 426     if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
 427   }
 428   else if (d == 2 && r == 'I') {
 429     if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
 430     else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
 431   }
 432   else if (d == 2 && r == 'O') {
 433     if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
 434   }
 435   else if (d == 3 && r == 'O') {
 436     if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
 437   }
 438   else if (d == 3 && r == 'I') {
 439     if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
 440   }
 441   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
 442   return other;
 443 }
 444 //____________________________________________________________________
 445 Int_t
 446 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
 447 {
 448   //
 449   // Get the possibly overlapping histogram of eta bin @a e in
 450   // detector and ring
 451   //
 452   // Parameters:
 453   //    i Ring index
 454   //    e Eta bin
 455   //    v Vertex bin (1 based)
 456   //
 457   // Return:
 458   //    Overlapping histogram index or -1
 459   //
 460   UShort_t d = 0;
 461   Char_t   r = '\0';
 462   GetDetRing(idx, d, r);
 463   if (d == 0 || r == '\0') return 0;
 464
 465   return GetOverlap(d, r, bin, vtxbin);
 466 }
 467
 468
 469 //____________________________________________________________________
 470 void
 471 AliFMDHistCollector::MergeBins(Double_t c,   Double_t e,
 472                                Double_t oc,  Double_t oe,
 473                                Double_t& rc, Double_t& re) const
 474 {
 475   //
 476   // Merge bins accoring to set method
 477   //
 478   // Parameters:
 479   //    c   Current content
 480   //    e   Current error
 481   //    oc  Old content
 482   //    oe  Old error
 483   //    rc  On return, the new content
 484   //    re  On return, tne new error
 485   //
 486   rc = re = 0;
 487   switch (fMergeMethod) {
 488   case kStraightMean:
 489     // calculate the average of old value (half the original),
 490     // and this value, as well as the summed squared errors
 491     // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
 492     // and half the error of this.
 493     //
 494     // So, on the first overlapping histogram we get
 495     //
 496     //    c = c_1 / 2
 497     //    e = sqrt((e_1 / 2)^2) = e_1/2
 498     //
 499     // On the second we get
 500     //
 501     //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
 502     //    e' = sqrt(e^2 + (e_2/2)^2)
 503     //       = sqrt(e_1^2/4 + e_2^2/4)
 504     //       = sqrt(1/4 * (e_1^2+e_2^2))
 505     //       = 1/2 * sqrt(e_1^2 + e_2^2)
 506     rc = oc + c/2;
 507     re = TMath::Sqrt(oe*oe+(e*e)/4);
 508     break;
 509   case kStraightMeanNoZero:
 510     // If there's data in the overlapping histogram,
 511     // calculate the average and add the errors in
 512     // quadrature.
 513     //
 514     // If there's no data in the overlapping histogram,
 515     // then just fill in the data
 516     if (oe > 0) {
 517       rc = (oc + c)/2;
 518       re = TMath::Sqrt(oe*oe + e*e)/2;
 519     }
 520     else {
 521       rc = c;
 522       re = e;
 523     }
 524     break;
 525   case kWeightedMean: {
 526     // Calculate the weighted mean
 527     Double_t w  = 1/(e*e);
 528     Double_t sc = w * c;
 529     Double_t sw = w;
 530     if (oe > 0) {
 531       Double_t ow =  1/(oe*oe);
 532       sc          += ow * oc;
 533       sw          += ow;
 534     }
 535     rc = sc / sw;
 536     re = TMath::Sqrt(1 / sw);
 537   }
 538     break;
 539   case kLeastError:
 540     if (e < oe) {
 541       rc = c;
 542       re = e;
 543     }
 544     else {
 545       rc = oc;
 546       re = oe;
 547     }
 548     break;
 549   default:
 550     AliError("No method for defining content of overlapping bins defined");
 551     return;
 552   }
 553 }
 554
 555 //____________________________________________________________________
 556 Bool_t
 557 AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
 558                              AliForwardUtil::Histos& sums,
 559                              UShort_t                vtxbin,
 560                              TH2D&                   out)
 561 {
 562   //
 563   // Do the calculations
 564   //
 565   // Parameters:
 566   //    hists    Cache of histograms
 567   //    vtxBin   Vertex bin (1 based)
 568   //    out      Output histogram
 569   //
 570   // Return:
 571   //    true on successs
 572   //
 573   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 574   const TAxis* vtxAxis = fcm.GetVertexAxis();
 575   Double_t vMin    = vtxAxis->GetBinLowEdge(vtxbin);
 576   Double_t vMax    = vtxAxis->GetBinUpEdge(vtxbin);
 577   TList* vtxList
 578     = static_cast<TList*>(fList->FindObject(Form("%c%02d_%c%02d",
 579                                                 vMin < 0 ? 'm' : 'p',
 580                                                 int(TMath::Abs(vMin)),
 581                                                 vMax < 0 ? 'm' : 'p',
 582                                                 int(TMath::Abs(vMax)))));
 583
 584
 585   for (UShort_t d=1; d<=3; d++) {
 586     UShort_t nr = (d == 1 ? 1 : 2);
 587     for (UShort_t q=0; q<nr; q++) {
 588       Char_t      r = (q == 0 ? 'I' : 'O');
 589       TH2D*       h = hists.Get(d,r);
 590       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
 591       Int_t       i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
 592       TH2D*       o = sums.Get(d, r);
 593
 594       // Skipping selected FMD rings
 595       if(d==1 && r=='I' && (fSkipFMDRings & kFMD1I)) { delete t; continue; }
 596       if(d==2 && r=='I' && (fSkipFMDRings & kFMD2I)) { delete t; continue; }
 597       if(d==2 && r=='O' && (fSkipFMDRings & kFMD2O)) { delete t; continue; }
 598       if(d==3 && r=='I' && (fSkipFMDRings & kFMD3I)) { delete t; continue; }
 599       if(d==3 && r=='O' && (fSkipFMDRings & kFMD3O)) { delete t; continue; }
 600
 601       // Get valid range
 602       Int_t first = 0;
 603       Int_t last  = 0;
 604       GetFirstAndLast(d, r, vtxbin, first, last);
 605
 606       // Zero outside valid range
 607       for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) {
 608         // Lower range
 609         for (Int_t iEta = 1; iEta < first; iEta++) {
 610           t->SetBinContent(iEta,iPhi,0);
 611           t->SetBinError(iEta,iPhi,0);
 612         }
 613         for (Int_t iEta = last+1; iEta <= t->GetNbinsX(); iEta++) {
 614           t->SetBinContent(iEta,iPhi,0);
 615           t->SetBinError(iEta,iPhi,0);
 616         }
 617       }
 618       for (Int_t iEta = first; iEta <= last; iEta++)
 619         t->SetBinContent(iEta,0,1);
 620       // Add to our per-ring sum
 621       o->Add(t);
 622
 623       // Outer rings have better phi segmentation - rebin to same as inner.
 624       if (q == 1) t->RebinY(2);
 625
 626       // Now update profile output
 627       for (Int_t iEta = first; iEta <= last; iEta++) {
 628
 629         // Get the possibly overlapping histogram
 630         Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
 631
 632         // Fill eta acceptance for this event into the phi underlow bin
 633         Float_t ooc      = out.GetBinContent(iEta,0);
 634         Float_t noc      = overlap >= 0? 0.5 : 1;
 635         out.SetBinContent(iEta, 0, ooc + noc);
 636
 637         // Should we loop over h or t Y bins - I think it's t
 638         for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) {
 639           Double_t c  = t->GetBinContent(iEta,iPhi);
 640           Double_t e  = t->GetBinError(iEta,iPhi);
 641           Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
 642           fSumRings->Fill(ee, i, c);
 643
 644           // If there's no signal, continue
 645           // if (e <= 0) continue;
 646           if (c <= 0 || e <= 0)     continue;
 647
 648           // If there's no overlapping histogram (ring), then
 649           // fill in data and continue to the next phi bin
 650           if (overlap < 0) {
 651             out.SetBinContent(iEta,iPhi,c);
 652             out.SetBinError(iEta,iPhi,e);
 653             continue;
 654           }
 655
 656           // Get the current bin content and error
 657           Double_t oc = out.GetBinContent(iEta,iPhi);
 658           Double_t oe = out.GetBinError(iEta,iPhi);
 659
 660           Double_t rc, re;
 661           MergeBins(c, e, oc, oe, rc, re);
 662           out.SetBinContent(iEta,iPhi, rc);
 663           out.SetBinError(iEta,iPhi, re);
 664         }
 665       }
 666       // Remove temporary histogram
 667       if(fBgAndHitMaps) {
 668         TH2D* hRingSumVtx
 669           = static_cast<TH2D*>(vtxList->FindObject(Form("hitMapFMD%d%c", d, r)));
 670         hRingSumVtx->Add(t);
 671       }
 672       delete t;
 673     } // for r
 674   } // for d
 675   return kTRUE;
 676 }
 677
 678 //____________________________________________________________________
 679 void
 680 AliFMDHistCollector::Print(Option_t* /* option */) const
 681 {
 682   //
 683   // Print information
 684   //
 685   // Parameters:
 686   //    option Not used
 687   //
 688   char ind[gROOT->GetDirLevel()+1];
 689   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 690   ind[gROOT->GetDirLevel()] = '\0';
 691   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 692             << ind << " # of cut bins:          " << fNCutBins << '\n'
 693             << ind << " Correction cut:         " << fCorrectionCut << '\n'
 694             << ind << " Merge method:           ";
 695   switch (fMergeMethod) {
 696   case kStraightMean:       std::cout << "straight mean\n"; break;
 697   case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
 698   case kWeightedMean:       std::cout << "weighted mean\n"; break;
 699   case kLeastError:         std::cout << "least error\n"; break;
 700   }
 701
 702   std::cout << ind << " Bin ranges:\n" << ind << "  rings  ";
 703   Int_t nVz = fFirstBins.fN / 5;
 704   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 705     UShort_t d = 0;
 706     Char_t   r = 0;
 707     GetDetRing(iIdx, d, r);
 708     std::cout << ind << " | FMD" << d << r << " ";
 709   }
 710   std::cout << '\n' << ind << "  /vz_bin ";
 711   for (Int_t iIdx = 0; iIdx < 5; iIdx++)
 712     std::cout << "-+--------";
 713   std::cout << std::endl;
 714
 715   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 716     std::cout << " " << std::setw(7) << iVz << "   ";
 717     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 718       UShort_t d = 0;
 719       Char_t   r = 0;
 720       GetDetRing(iIdx, d, r);
 721
 722       Int_t first, last;
 723       GetFirstAndLast(iIdx, iVz, first, last);
 724       std::cout << " | " << std::setw(3) << first << "-"
 725                 << std::setw(3) << last;
 726     }
 727     std::cout << std::endl;
 728   }
 729 }
 730
 731 //____________________________________________________________________
 732 //
 733 // EOF
 734 //
 735
 736
 737