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