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 "AliFMDCorrSecondaryMap.h"
  22 #include "AliLog.h"
  23 #include <TH2D.h>
  24 #include <TH3D.h>
  25 #include <TH1I.h>
  26 #include <TProfile.h>
  27 #include <TProfile2D.h>
  28 #include <TObjArray.h>
  29 #include <TArrayI.h>
  30 #include <TROOT.h>
  31 #include <iostream>
  32 #include <iomanip>
  33
  34 ClassImp(AliFMDHistCollector)
  35 #if 0
  36 ; // For Emacs
  37 #endif
  38
  39 //____________________________________________________________________
  40 AliFMDHistCollector::AliFMDHistCollector()
  41   : fNCutBins(0),
  42     fCorrectionCut(0),
  43     fDebug(0),
  44     fList(0),
  45     fSumRings(0),
  46     fCoverage(0),
  47     fMergeMethod(kStraightMean),
  48     fFiducialMethod(kByCut),
  49     fSkipFMDRings(0),
  50     fBgAndHitMaps(false),
  51     fVtxList(0),
  52     fByCent(0),
  53     fDoByCent(false)
  54 {
  55   DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector");
  56 }
  57
  58 //____________________________________________________________________
  59 AliFMDHistCollector::AliFMDHistCollector(const char* title)
  60   : TNamed("fmdHistCollector", title),
  61     fNCutBins(2),
  62     fCorrectionCut(0.5),
  63     fDebug(0),
  64     fList(0),
  65     fSumRings(0),
  66     fCoverage(0),
  67     fMergeMethod(kStraightMean),
  68     fFiducialMethod(kByCut),
  69     fSkipFMDRings(0),
  70     fBgAndHitMaps(false),
  71     fVtxList(0),
  72     fByCent(0),
  73     fDoByCent(false)
  74 {
  75   DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title);
  76 }
  77 //____________________________________________________________________
  78 AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
  79   : TNamed(o),
  80     fNCutBins(o.fNCutBins),
  81     fCorrectionCut(o.fCorrectionCut),
  82     fDebug(o.fDebug),
  83     fList(o.fList),
  84     fSumRings(o.fSumRings),
  85     fCoverage(o.fCoverage),
  86     fMergeMethod(o.fMergeMethod),
  87     fFiducialMethod(o.fFiducialMethod),
  88     fSkipFMDRings(o.fSkipFMDRings),
  89     fBgAndHitMaps(o.fBgAndHitMaps),
  90     fVtxList(o.fVtxList),
  91     fByCent(o.fByCent),
  92     fDoByCent(o.fDoByCent)
  93 {
  94   DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
  95 }
  96
  97 //____________________________________________________________________
  98 AliFMDHistCollector::~AliFMDHistCollector()
  99 {
 100   DGUARD(fDebug, 3, "DTOR of AliFMDHistCollector");
 101   // if (fList) delete fList;
 102 }
 103 //____________________________________________________________________
 104 AliFMDHistCollector&
 105 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
 106 {
 107   //
 108   // Assignement operator
 109   //
 110   // Parameters:
 111   //    o Object to assign from
 112   //
 113   // Return:
 114   //    Reference to this object
 115   //
 116   DGUARD(fDebug, 3, "Assignment of AliFMDHistCollector");
 117   if (&o == this) return *this;
 118   TNamed::operator=(o);
 119
 120   fNCutBins       = o.fNCutBins;
 121   fCorrectionCut  = o.fCorrectionCut;
 122   fDebug          = o.fDebug;
 123   fList           = o.fList;
 124   fSumRings       = o.fSumRings;
 125   fCoverage       = o.fCoverage;
 126   fMergeMethod    = o.fMergeMethod;
 127   fFiducialMethod = o.fFiducialMethod;
 128   fSkipFMDRings   = o.fSkipFMDRings;
 129   fBgAndHitMaps   = o.fBgAndHitMaps;
 130   fVtxList        = o.fVtxList;
 131   fByCent         = o.fByCent;
 132   fDoByCent       = o.fDoByCent;
 133   return *this;
 134 }
 135
 136 //____________________________________________________________________
 137 void
 138 AliFMDHistCollector::SetupForData(const TAxis& vtxAxis,
 139                                   const TAxis& etaAxis)
 140 {
 141   //
 142   // Intialise
 143   //
 144   // Parameters:
 145   //    vtxAxis  Vertex axis
 146   //
 147   DGUARD(fDebug, 1, "Initialization of AliFMDHistCollector");
 148
 149   // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 150
 151   fSumRings = new TH2D("sumRings", "Sum in individual rings",
 152                        etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
 153                        5, 1, 6);
 154   fSumRings->Sumw2();
 155   fSumRings->SetDirectory(0);
 156   fSumRings->SetXTitle("#eta");
 157   fSumRings->GetYaxis()->SetBinLabel(1,"FMD1i");
 158   fSumRings->GetYaxis()->SetBinLabel(2,"FMD2i");
 159   fSumRings->GetYaxis()->SetBinLabel(3,"FMD2o");
 160   fSumRings->GetYaxis()->SetBinLabel(4,"FMD3i");
 161   fSumRings->GetYaxis()->SetBinLabel(5,"FMD3o");
 162   fList->Add(fSumRings);
 163
 164   fCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
 165                        etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
 166                        vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
 167   fCoverage->SetDirectory(0);
 168   fCoverage->SetXTitle("#eta");
 169   fCoverage->SetYTitle("v_{z} [cm]");
 170   fCoverage->SetZTitle("n_{bins}");
 171   fList->Add(fCoverage);
 172
 173
 174   // --- Add parameters to output ------------------------------------
 175   fList->Add(AliForwardUtil::MakeParameter("nCutBins",fNCutBins));
 176   fList->Add(AliForwardUtil::MakeParameter("skipRings",fSkipFMDRings));
 177   fList->Add(AliForwardUtil::MakeParameter("bgAndHits",fBgAndHitMaps));
 178   fList->Add(AliForwardUtil::MakeParameter("fiducial",Int_t(fFiducialMethod)));
 179   fList->Add(AliForwardUtil::MakeParameter("fiducialCut",fCorrectionCut));
 180   fList->Add(AliForwardUtil::MakeParameter("merge",Int_t(fMergeMethod)));
 181
 182   UShort_t nVz = vtxAxis.GetNbins();
 183   fVtxList = new TObjArray(nVz, 1);
 184   fVtxList->SetName("histCollectorVtxBins");
 185   fVtxList->SetOwner();
 186
 187   // Find the eta bin ranges
 188   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 189     Double_t vMin    = vtxAxis.GetBinLowEdge(iVz);
 190     Double_t vMax    = vtxAxis.GetBinUpEdge(iVz);
 191     VtxBin*  bin     = new VtxBin(iVz, vMin, vMax, fNCutBins);
 192     fVtxList->AddAt(bin, iVz);
 193
 194     bin->SetupForData(fCoverage, fSkipFMDRings, fFiducialMethod,
 195                       fCorrectionCut, fList, etaAxis,
 196                       fBgAndHitMaps, fBgAndHitMaps);
 197   }
 198
 199   if (!fDoByCent) return;
 200
 201   fByCent = new TList;
 202   fByCent->SetName("byCentrality");
 203   fByCent->SetOwner();
 204   fList->Add(fByCent);
 205
 206   Int_t    nCent   = 101;
 207   Double_t minCent = -.5;
 208   Double_t maxCent = 100.5;
 209   for (Int_t i = 0; i < 5; i++) {
 210     UShort_t d;
 211     Char_t   r;
 212     GetDetRing(i, d, r);
 213
 214     TH3* h = new TH3D(Form("FMD%d%c", d, r),
 215                       Form("dN/d#eta per centrality for FMD%d%c", d, r),
 216                       etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
 217                       nCent, minCent, maxCent, 1, 0, 1);
 218     h->SetXTitle("#eta");
 219     h->SetYTitle("Centrality [%]");
 220     h->SetZTitle("dN/d#eta");
 221     h->SetDirectory(0);
 222     h->SetMarkerColor(AliForwardUtil::RingColor(d, r));
 223     h->SetMarkerStyle(20);
 224     fByCent->Add(h);
 225   }
 226 }
 227 //____________________________________________________________________
 228 Bool_t
 229 AliFMDHistCollector::CheckCorrection(FiducialMethod m,
 230                                      Double_t       cut,
 231                                      const TH2D*    bg,
 232                                      Int_t          ie,
 233                                      Int_t          ip)
 234 {
 235   //
 236   // Check if we should include the bin in the data range
 237   //
 238   // Parameters:
 239   //    bg Secondary map histogram
 240   //    ie Eta bin
 241   //    ip Phi bin
 242   //
 243   // Return:
 244   //    True if to be used
 245   //
 246   Double_t c = bg->GetBinContent(ie,ip);
 247   switch (m) {
 248   case kByCut:
 249     return c >= cut;
 250   case kDistance:
 251     if (2 * c < bg->GetBinContent(ie+1,ip) ||
 252         2 * c < bg->GetBinContent(ie-1,ip)) return false;
 253     return true;
 254   default:
 255     AliErrorClass("No fiducal cut method defined");
 256   }
 257   return false;
 258 }
 259
 260 //____________________________________________________________________
 261 void
 262 AliFMDHistCollector::CreateOutputObjects(TList* dir)
 263 {
 264   //
 265   // Output diagnostic histograms to directory
 266   //
 267   // Parameters:
 268   //    dir List to write in
 269   //
 270   DGUARD(fDebug, 1, "Define output of AliFMDHistCollector");
 271   fList = new TList;
 272   fList->SetOwner();
 273   fList->SetName(GetName());
 274   dir->Add(fList);
 275
 276 }
 277
 278 //____________________________________________________________________
 279 Bool_t
 280 AliFMDHistCollector::CheckSkip(UShort_t d, Char_t r, UShort_t skips)
 281 {
 282   UShort_t q  = (r == 'I' || r == 'i' ? 0 : 1);
 283   UShort_t c = 1 << (d-1);
 284   UShort_t t = 1 << (c+q-1);
 285
 286   return (t & skips) == t;
 287 }
 288
 289 //____________________________________________________________________
 290 Int_t
 291 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r)
 292 {
 293   //
 294   // Get the ring index from detector number and ring identifier
 295   //
 296   // Parameters:
 297   //    d Detector
 298   //    r Ring identifier
 299   //
 300   // Return:
 301   //    ring index or -1 in case of problems
 302   //
 303   Int_t idx = -1;
 304   switch (d) {
 305   case 1: idx = 0; break;
 306   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 307   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 308   }
 309   return idx;
 310 }
 311 //____________________________________________________________________
 312 void
 313 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r)
 314 {
 315   //
 316   // Get the detector and ring from the ring index
 317   //
 318   // Parameters:
 319   //    idx Ring index
 320   //    d   On return, the detector or 0 in case of errors
 321   //    r   On return, the ring id or '0' in case of errors
 322   //
 323   d = 0;
 324   r = '\0';
 325   switch (idx) {
 326   case 0: d = 1; r = 'I'; break;
 327   case 1: d = 2; r = 'I'; break;
 328   case 2: d = 2; r = 'O'; break;
 329   case 3: d = 3; r = 'I'; break;
 330   case 4: d = 3; r = 'O'; break;
 331   }
 332 }
 333
 334 //____________________________________________________________________
 335 AliFMDHistCollector::VtxBin*
 336 AliFMDHistCollector::GetVtxBin(Int_t ivtx)
 337 {
 338   // Parameters:
 339   //    vtxBin   Vertex bin (1 based)
 340   if (!fVtxList) return 0;
 341   if (ivtx < 1 || ivtx > fVtxList->GetEntriesFast()) return 0;
 342   VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivtx));
 343   return bin;
 344 }
 345 //____________________________________________________________________
 346 const AliFMDHistCollector::VtxBin*
 347 AliFMDHistCollector::GetVtxBin(Int_t ivtx) const
 348 {
 349   // Parameters:
 350   //    vtxBin   Vertex bin (1 based)
 351   if (!fVtxList) return 0;
 352   if (ivtx < 1 || ivtx > fVtxList->GetEntriesFast()) return 0;
 353   VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivtx));
 354   return bin;
 355 }
 356
 357 //____________________________________________________________________
 358 void
 359 AliFMDHistCollector::MergeBins(MergeMethod   m,
 360                                Double_t c,   Double_t e,
 361                                Double_t oc,  Double_t oe,
 362                                Double_t& rc, Double_t& re)
 363 {
 364   //
 365   // Merge bins accoring to set method
 366   //
 367   // Parameters:
 368   //    c   Current content
 369   //    e   Current error
 370   //    oc  Old content
 371   //    oe  Old error
 372   //    rc  On return, the new content
 373   //    re  On return, tne new error
 374   //
 375   rc = re = 0;
 376   switch (m) {
 377   case kStraightMean:
 378   case kPreferInner:  // We only get these two when there's an overlap
 379   case kPreferOuter:  // between like rings, and we should do a straight mean
 380     // calculate the average of old value (half the original),
 381     // and this value, as well as the summed squared errors
 382     // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
 383     // and half the error of this.
 384     //
 385     // So, on the first overlapping histogram we get
 386     //
 387     //    c = c_1 / 2
 388     //    e = sqrt((e_1 / 2)^2) = e_1/2
 389     //
 390     // On the second we get
 391     //
 392     //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
 393     //    e' = sqrt(e^2 + (e_2/2)^2)
 394     //       = sqrt(e_1^2/4 + e_2^2/4)
 395     //       = sqrt(1/4 * (e_1^2+e_2^2))
 396     //       = 1/2 * sqrt(e_1^2 + e_2^2)
 397     rc = oc + c/2;
 398     re = TMath::Sqrt(oe*oe+(e*e)/4);
 399     break;
 400   case kStraightMeanNoZero:
 401     // If there's data in the overlapping histogram,
 402     // calculate the average and add the errors in
 403     // quadrature.
 404     //
 405     // If there's no data in the overlapping histogram,
 406     // then just fill in the data
 407     if (oe > 0) {
 408       rc = (oc + c)/2;
 409       re = TMath::Sqrt(oe*oe + e*e)/2;
 410     }
 411     else {
 412       rc = c;
 413       re = e;
 414     }
 415     break;
 416   case kWeightedMean: {
 417     // Calculate the weighted mean
 418     Double_t w  = 1/(e*e);
 419     Double_t sc = w * c;
 420     Double_t sw = w;
 421     if (oe > 0) {
 422       Double_t ow =  1/(oe*oe);
 423       sc          += ow * oc;
 424       sw          += ow;
 425     }
 426     rc = sc / sw;
 427     re = TMath::Sqrt(1 / sw);
 428   }
 429     break;
 430   case kLeastError:
 431     if (e < oe) {
 432       rc = c;
 433       re = e;
 434     }
 435     else {
 436       rc = oc;
 437       re = oe;
 438     }
 439     break;
 440   case kSum:
 441     rc = c + oc;
 442     re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature
 443     break;
 444   default:
 445     AliErrorClass("No method for defining content of overlapping bins defined");
 446     return;
 447   }
 448 }
 449
 450 //____________________________________________________________________
 451 Bool_t
 452 AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
 453                              AliForwardUtil::Histos& sums,
 454                              UShort_t                vtxbin,
 455                              TH2D&                   out,
 456                              Double_t                cent)
 457 {
 458   //
 459   // Do the calculations
 460   //
 461   // Parameters:
 462   //    hists    Cache of histograms
 463   //    vtxBin   Vertex bin (1 based)
 464   //    out      Output histogram
 465   //
 466   // Return:
 467   //    true on successs
 468   //
 469   DGUARD(fDebug, 1, "Collect final histogram of AliFMDHistCollector");
 470   // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 471   // const TAxis* vtxAxis = fcm.GetVertexAxis();
 472   // Double_t vMin    = vtxAxis->GetBinLowEdge(vtxbin);
 473   // Double_t vMax    = vtxAxis->GetBinUpEdge(vtxbin);
 474   VtxBin*  bin     = GetVtxBin(vtxbin);
 475   if (!bin) return false;
 476   Bool_t   ret     = bin->Collect(hists, sums, out, fSumRings, cent,
 477                                   fMergeMethod, fSkipFMDRings,
 478                                   fByCent);
 479
 480   return ret;
 481 }
 482
 483 //____________________________________________________________________
 484 void
 485 AliFMDHistCollector::Print(Option_t* /* option */) const
 486 {
 487   //
 488   // Print information
 489   //
 490   // Parameters:
 491   //    option Not used
 492   //
 493   char ind[gROOT->GetDirLevel()+1];
 494   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 495   ind[gROOT->GetDirLevel()] = '\0';
 496   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 497             << ind << " # of cut bins:          " << fNCutBins << '\n'
 498             << ind << " Fiducal method:         "
 499             << (fFiducialMethod == kByCut ? "cut" : "distance") << "\n"
 500             << ind << " Fiducial cut:           " << fCorrectionCut << "\n"
 501             << ind << " Merge method:           ";
 502   switch (fMergeMethod) {
 503   case kStraightMean:       std::cout << "straight mean\n"; break;
 504   case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
 505   case kWeightedMean:       std::cout << "weighted mean\n"; break;
 506   case kLeastError:         std::cout << "least error\n"; break;
 507   case kSum:                std::cout << "straight sum\n"; break;
 508   case kPreferInner:        std::cout << "prefer inners\n"; break;
 509   case kPreferOuter:        std::cout << "prefer outers\n"; break;
 510   }
 511
 512   if (!fVtxList) return;
 513
 514   std::cout << ind << " Bin ranges:\n" << ind << "  rings   |   Range  ";
 515   Int_t nVz = fVtxList->GetEntriesFast();
 516   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 517     UShort_t d = 0;
 518     Char_t   r = 0;
 519     GetDetRing(iIdx, d, r);
 520     std::cout << ind << " | FMD" << d << r << " ";
 521   }
 522   std::cout << '\n' << ind << "  /vz_bin |-----------";
 523   for (Int_t iIdx = 0; iIdx < 5; iIdx++)
 524     std::cout << "-+--------";
 525   std::cout << std::endl;
 526
 527   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 528     const VtxBin* bin = GetVtxBin(iVz);
 529     if (!bin) continue;
 530     std::cout << "    " << std::right << std::setw(6) << iVz << " | "
 531               << std::setw(3) << bin->fLow << " - " << std::left
 532               << std::setw(3) << bin->fHigh << " ";
 533     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 534       Int_t first, last;
 535       bin->GetFirstAndLast(iIdx, first, last);
 536       std::cout << " | " << std::setw(3) << first << "-"
 537                 << std::setw(3) << last;
 538     }
 539     std::cout << std::endl;
 540   }
 541 }
 542
 543 //____________________________________________________________________
 544 AliFMDHistCollector::VtxBin::VtxBin(Int_t idx, Double_t minIpZ, Double_t maxIpZ,
 545                                     Int_t nCutBins)
 546   : fIndex(idx),
 547     fLow(minIpZ),
 548     fHigh(maxIpZ),
 549     fHitMap(0),
 550     fFirstBin(1),
 551     fLastBin(1),
 552     fNCutBins(nCutBins)
 553 {
 554 }
 555 //____________________________________________________________________
 556 AliFMDHistCollector::VtxBin::VtxBin(const VtxBin& o)
 557   : TObject(o),
 558     fIndex(o.fIndex),
 559     fLow(o.fLow),
 560     fHigh(o.fHigh),
 561     fHitMap(o.fHitMap),
 562     fFirstBin(o.fFirstBin),
 563     fLastBin(o.fLastBin),
 564     fNCutBins(o.fNCutBins)
 565 {
 566 }
 567 //____________________________________________________________________
 568 AliFMDHistCollector::VtxBin&
 569 AliFMDHistCollector::VtxBin::operator=(const VtxBin& o)
 570 {
 571   if (&o == this) return *this;
 572   fIndex    = o.fIndex;
 573   fLow      = o.fLow;
 574   fHigh     = o.fHigh;
 575   fHitMap   = o.fHitMap;
 576   fFirstBin = o.fFirstBin;
 577   fLastBin  = o.fLastBin;
 578   fNCutBins = o.fNCutBins;
 579   return *this;
 580 }
 581 //____________________________________________________________________
 582 const Char_t*
 583 AliFMDHistCollector::VtxBin::GetName() const
 584 {
 585   return Form("%c%03d_%c%03d",
 586               (fLow >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fLow)),
 587               (fHigh >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fHigh)));
 588 }
 589 //____________________________________________________________________
 590 void
 591 AliFMDHistCollector::VtxBin::SetupForData(TH2*           coverage,
 592                                           UShort_t       skips,
 593                                           FiducialMethod fiducial,
 594                                           Double_t       cut,
 595                                           TList*         l,
 596                                           const TAxis&   etaAxis,
 597                                           Bool_t         doHitMaps,
 598                                           Bool_t         storeSecMap)
 599 {
 600   TList* out = 0;
 601   if (doHitMaps || storeSecMap) {
 602     out = new TList;
 603     out->SetName(GetName());
 604     out->SetOwner();
 605     l->Add(out);
 606   }
 607   if (doHitMaps) {
 608     fHitMap = new AliForwardUtil::Histos();
 609     fHitMap->Init(etaAxis);
 610   }
 611   fFirstBin.Set(5);
 612   fLastBin.Set(5);
 613
 614   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 615
 616   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 617     UShort_t d = 0;
 618     Char_t   r = 0;
 619     GetDetRing(iIdx, d, r);
 620
 621     // Skipping selected FMD rings
 622     if (CheckSkip(d, r, skips)) continue;
 623
 624     // Get the background object
 625     TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,UShort_t(fIndex));
 626     Int_t nEta  = bg->GetNbinsX();
 627     Int_t first = nEta+1;
 628     Int_t last  = 0;
 629
 630     // Loop over the eta bins
 631     for (Int_t ie = 1; ie <= nEta; ie++) {
 632       // Loop over the phi bins to make sure that we
 633       // do not have holes in the coverage
 634       bool ok = true;
 635       for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 636         if (!CheckCorrection(fiducial, cut, bg, ie, ip)) {
 637           ok = false;
 638           continue;
 639         }
 640       }
 641       if (!ok) continue;
 642
 643       first = TMath::Min(ie, first);
 644       last  = TMath::Max(ie, last);
 645     }
 646     // Store result of first/last bin for this ring
 647     fFirstBin[iIdx] = first;
 648     fLastBin[iIdx]  = last;
 649
 650     if (fHitMap) {
 651       TH2* h = fHitMap->Get(d, r);
 652       h->SetDirectory(0);
 653       h->SetName(Form("hitMapFMD%d%c", d, r));
 654       // if (r == 'O') h->RebinY(2);
 655       out->Add(h);
 656     }
 657
 658     TH2D* obg=0;
 659     if(storeSecMap) {
 660       obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
 661       obg->SetDirectory(0);
 662       obg->Reset();
 663       out->Add(obg);
 664     }
 665     // Fill diagnostics histograms
 666     for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
 667       Double_t old = coverage->GetBinContent(ie, fIndex);
 668       coverage->SetBinContent(ie, fIndex, old+1);
 669       if(obg) {
 670         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 671           obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
 672           obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
 673         } // for (ip)
 674       } // if (doSecHits)
 675     } // for (ie)
 676   } // for (iIdx)
 677 }
 678
 679 //____________________________________________________________________
 680 void
 681 AliFMDHistCollector::VtxBin::GetFirstAndLast(Int_t  idx,
 682                                              Int_t& first,
 683                                              Int_t& last) const
 684 {
 685   // Get the first and last eta bin to use for a given ring and vertex
 686   //
 687   // Parameters:
 688   //    idx      Ring index as given by GetIdx
 689   //    first    On return, the first eta bin to use
 690   //    last     On return, the last eta bin to use
 691   //
 692   first = 0;
 693   last  = 0;
 694
 695   if (idx < 0 || idx >= fFirstBin.GetSize()) return;
 696
 697   first = fFirstBin.At(idx)+fNCutBins;
 698   last  = fLastBin.At(idx)-fNCutBins;
 699 }
 700 //____________________________________________________________________
 701 Int_t
 702 AliFMDHistCollector::VtxBin::GetFirst(Int_t  idx) const
 703 {
 704   Int_t first, last;
 705   GetFirstAndLast(idx, first , last);
 706   return first;
 707 }
 708 //____________________________________________________________________
 709 Int_t
 710 AliFMDHistCollector::VtxBin::GetLast(Int_t  idx) const
 711 {
 712   Int_t first, last;
 713   GetFirstAndLast(idx, first , last);
 714   return last;
 715 }
 716
 717 //____________________________________________________________________
 718 Bool_t
 719 AliFMDHistCollector::VtxBin::Collect(const AliForwardUtil::Histos& hists,
 720                                      AliForwardUtil::Histos&       sums,
 721                                      TH2D&                         out,
 722                                      TH2D*                         sumRings,
 723                                      Double_t                      cent,
 724                                      MergeMethod                   m,
 725                                      UShort_t                      skips,
 726                                      TList*                        byCent)
 727 {
 728   for (UShort_t d=1; d<=3; d++) {
 729     UShort_t nr = (d == 1 ? 1 : 2);
 730     for (UShort_t q=0; q<nr; q++) {
 731       Char_t      r = (q == 0 ? 'I' : 'O');
 732       if (CheckSkip(d, r, skips)) continue;
 733
 734       TH2D*       h = hists.Get(d,r);
 735       TH2D*       o = sums.Get(d, r);
 736       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
 737       Int_t       i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
 738
 739       // Get valid range
 740       Int_t first = 0;
 741       Int_t last  = 0;
 742       GetFirstAndLast(d, r, first, last);
 743
 744       // Zero outside valid range
 745       Int_t nY = t->GetNbinsY();
 746       Int_t nX = t->GetNbinsX();
 747       for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) {
 748         // Lower range
 749         for (Int_t iEta = 1; iEta < first; iEta++) {
 750           t->SetBinContent(iEta,iPhi,0);
 751           t->SetBinError(iEta,iPhi,0);
 752         }
 753         for (Int_t iEta = last+1; iEta <= nX; iEta++) {
 754           t->SetBinContent(iEta,iPhi,0);
 755           t->SetBinError(iEta,iPhi,0);
 756         }
 757       }
 758       // Fill under-flow bins with eta coverage
 759       for (Int_t iEta = first; iEta <= last; iEta++)
 760         t->SetBinContent(iEta,0,1);
 761
 762       // Add to our per-ring sum
 763       o->Add(t);
 764
 765       // If we store hit maps, update here
 766       if (fHitMap) fHitMap->Get(d, r)->Add(t);
 767
 768       if (byCent) {
 769         TH3* dNdetaCent = static_cast<TH3*>(byCent->At(i-1));
 770         if (cent >= 0 && dNdetaCent) {
 771           Int_t iCent = dNdetaCent->GetYaxis()->FindBin(cent);
 772
 773           if (iCent > 0 && iCent <= dNdetaCent->GetNbinsY()) {
 774             // Make a projection of data
 775             TH1* proj = static_cast<TH1*>(t->ProjectionX("tmp", 1, nY));
 776             proj->SetDirectory(0);
 777             for (Int_t iEta = 1; iEta <= nX; iEta++) {
 778               Double_t v1 = proj->GetBinContent(iEta);
 779               Double_t e1 = proj->GetBinError(iEta);
 780               Double_t v2 = dNdetaCent->GetBinContent(iEta, iCent, 1);
 781               Double_t e2 = dNdetaCent->GetBinError(iEta, iCent, 1);
 782               dNdetaCent->SetBinContent(iEta,iCent,1, v1+v2);
 783               dNdetaCent->SetBinError(iEta,iCent,1, TMath::Sqrt(e1*e1+e2*e2));
 784
 785               // Check under-/overflow bins
 786               Double_t uF = t->GetBinContent(iEta, 0);
 787               Double_t oF = t->GetBinContent(iEta, nY+1);
 788               if (uF > 0) {
 789                 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 0);
 790                 dNdetaCent->SetBinContent(iEta, iCent, 0, old + uF);
 791               }
 792               if (oF > 0) {
 793                 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 2);
 794                 dNdetaCent->SetBinContent(iEta, iCent, 2, old + oF);
 795               }
 796             } // for(iEta)
 797             delete proj;
 798           } // if(iCent)
 799         } // if (cent)
 800       } // if (byCent)
 801
 802       // Outer rings have better phi segmentation - rebin to same as inner.
 803       if (q == 1) t->RebinY(2);
 804
 805       nY = t->GetNbinsY();
 806
 807       // Now update profile output
 808       for (Int_t iEta = first; iEta <= last; iEta++) {
 809
 810         // Get the possibly overlapping histogram
 811         Int_t overlap = GetOverlap(d,r,iEta);
 812
 813         // Get factor
 814         MergeMethod mm  = m; // Possibly override method locally
 815         Float_t     fac = 1;
 816         if (m != kSum && overlap >= 0) {
 817           // Default is to average
 818           fac = 0.5;
 819           if (m == kPreferInner) {
 820             // Current one is an outer overlapping an inner
 821             if ((r == 'o' || r == 'O') &&
 822                 (overlap == 0 || overlap == 1 || overlap == 3))
 823               // Do not use this signal
 824               fac = 0;
 825             // Current one is inner overlapping an outer
 826             else if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
 827               // Prefer this one
 828               fac = 1;
 829             else
 830               // In case of two overlapping inners
 831               mm = kStraightMean;
 832           }
 833           else if (m == kPreferOuter) {
 834             // Current one is an inner overlapping an outer
 835             if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
 836               // Do not use this signal
 837               fac = 0;
 838             else if ((r == 'O' || r == 'o') &&
 839                      (overlap == 0 || overlap == 1 || overlap == 3))
 840               fac = 1;
 841             else
 842               // In case of two overlapping outers
 843               mm = kStraightMean;
 844           }
 845         }
 846
 847         // Fill eta acceptance for this event into the phi underflow bin
 848         Float_t ooc      = out.GetBinContent(iEta,0);
 849         out.SetBinContent(iEta, 0, ooc + fac);
 850
 851         // Fill phi acceptance for this event into the phi overflow bin
 852         Float_t oop      = out.GetBinContent(iEta,nY+1);
 853         Float_t nop      = t->GetBinContent(iEta,nY+1);
 854 #if 0
 855         Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f",
 856              iEta, fac, oop, nop, fac*(oop+nop));
 857 #endif
 858         out.SetBinContent(iEta, nY+1, fac * nop + oop);
 859
 860         // Should we loop over h or t Y bins - I think it's t
 861         for (Int_t iPhi = 1; iPhi <= nY; iPhi++) {
 862           Double_t c  = t->GetBinContent(iEta,iPhi);
 863           Double_t e  = t->GetBinError(iEta,iPhi);
 864           Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
 865           sumRings->Fill(ee, i, c);
 866
 867           // If there's no signal or the signal was ignored because we
 868           // prefer the inners/outers, continue if (e <= 0) continue;
 869           if (fac <= 0 || c <= 0 || e <= 0)     continue;
 870
 871           // If there's no overlapping histogram (ring) or if we
 872           // prefer inner/outer, then fill in data and continue to the
 873           // next phi bin
 874           if (overlap < 0 || fac >= 1) {
 875             out.SetBinContent(iEta,iPhi,c);
 876             out.SetBinError(iEta,iPhi,e);
 877             continue;
 878           }
 879
 880           // Get the current bin content and error
 881           Double_t oc = out.GetBinContent(iEta,iPhi);
 882           Double_t oe = out.GetBinError(iEta,iPhi);
 883
 884           Double_t rc, re;
 885           MergeBins(mm, c, e, oc, oe, rc, re);
 886           out.SetBinContent(iEta,iPhi, rc);
 887           out.SetBinError(iEta,iPhi, re);
 888         }
 889       }
 890       // Remove temporary histogram
 891       delete t;
 892     } // for r
 893   } // for d
 894   return true;
 895 }
 896 //____________________________________________________________________
 897 Int_t
 898 AliFMDHistCollector::VtxBin::GetOverlap(UShort_t d, Char_t r,
 899                                         Int_t bin) const
 900 {
 901   //
 902   // Get the possibly overlapping histogram of eta bin @a e in
 903   // detector and ring
 904   //
 905   // Parameters:
 906   //    d Detector
 907   //    r Ring
 908   //    e Eta bin
 909   //    v Vertex bin (1 based)
 910   //
 911   // Return:
 912   //    Overlapping histogram index or -1
 913   //
 914
 915   Int_t other = -1;
 916   if (d == 1) {
 917     if (bin <= GetLast(2,'I')) other = GetIdx(2,'I');
 918   }
 919   else if (d == 2 && r == 'I') {
 920     if      (bin <= GetLast(2,  'O')) other = GetIdx(2, 'O');
 921     else if (bin >= GetFirst(1, 'I')) other = GetIdx(1, 'I');
 922   }
 923   else if (d == 2 && r == 'O') {
 924     if (bin >= GetFirst(2, 'I'))      other = GetIdx(2,'I');
 925   }
 926   else if (d == 3 && r == 'O') {
 927     if (bin <= GetLast(3, 'I'))       other = GetIdx(3, 'I');
 928   }
 929   else if (d == 3 && r == 'I') {
 930     if (bin >= GetFirst(3, 'O'))      other = GetIdx(3, 'O');
 931   }
 932   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
 933   return other;
 934 }
 935
 936
 937 //____________________________________________________________________
 938 //
 939 // EOF
 940 //
 941
 942
 943