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 {
 474   //
 475   // Do the calculations
 476   //
 477   // Parameters:
 478   //    hists    Cache of histograms
 479   //    vtxBin   Vertex bin (1 based)
 480   //    out      Output histogram
 481   //
 482   // Return:
 483   //    true on successs
 484   //
 485   DGUARD(fDebug, 1, "Collect final histogram of AliFMDHistCollector");
 486   // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 487   // const TAxis* vtxAxis = fcm.GetVertexAxis();
 488   // Double_t vMin    = vtxAxis->GetBinLowEdge(vtxbin);
 489   // Double_t vMax    = vtxAxis->GetBinUpEdge(vtxbin);
 490   VtxBin*  bin     = GetVtxBin(vtxbin);
 491   if (!bin) return false;
 492   Bool_t   ret     = bin->Collect(hists, sums, out, fSumRings, fSkipped, cent,
 493                                   fMergeMethod, fSkipFMDRings,
 494                                   fByCent, eta2phi);
 495
 496   return ret;
 497 }
 498
 499 #define PF(N,V,...)                                     \
 500   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
 501 #define PFB(N,FLAG)                             \
 502   do {                                                                  \
 503     AliForwardUtil::PrintName(N);                                       \
 504     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 505   } while(false)
 506 #define PFV(N,VALUE)                                    \
 507   do {                                                  \
 508     AliForwardUtil::PrintName(N);                       \
 509     std::cout << (VALUE) << std::endl; } while(false)
 510
 511 //____________________________________________________________________
 512 void
 513 AliFMDHistCollector::Print(Option_t* /* option */) const
 514 {
 515   //
 516   // Print information
 517   //
 518   // Parameters:
 519   //    option Not used
 520   //
 521   TString merge("unknown");
 522   switch (fMergeMethod) {
 523   case kStraightMean:       merge = "straight mean"; break;
 524   case kStraightMeanNoZero: merge = "straight mean (no zeros)"; break;
 525   case kWeightedMean:       merge = "weighted mean"; break;
 526   case kLeastError:         merge = "least error"; break;
 527   case kSum:                merge = "straight sum"; break;
 528   case kPreferInner:        merge = "prefer inners"; break;
 529   case kPreferOuter:        merge = "prefer outers"; break;
 530   }
 531
 532   AliForwardUtil::PrintTask(*this);
 533   gROOT->IncreaseDirLevel();
 534   PFV("# of cut bins",  fNCutBins );
 535   PFV("Fiducal method", (fFiducialMethod == kByCut ? "cut" : "distance"));
 536   PFV("Fiducial cut",   fCorrectionCut );
 537   PFV("Merge method",   merge);
 538
 539   if (!fVtxList) {
 540     gROOT->DecreaseDirLevel();
 541     return;
 542   }
 543   char ind[64];
 544   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 545   ind[gROOT->GetDirLevel()] = '\0';
 546
 547   std::cout << ind << "Bin ranges:\n" << ind << "  rings   |   Range  ";
 548   Int_t nVz = fVtxList->GetEntriesFast();
 549   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 550     UShort_t d = 0;
 551     Char_t   r = 0;
 552     GetDetRing(iIdx, d, r);
 553     std::cout << ind << " | FMD" << d << r << " ";
 554   }
 555   std::cout << '\n' << ind << "  /vz_bin |-----------";
 556   for (Int_t iIdx = 0; iIdx < 5; iIdx++)
 557     std::cout << "-+--------";
 558   std::cout << std::endl;
 559
 560   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 561     const VtxBin* bin = GetVtxBin(iVz);
 562     if (!bin) continue;
 563     std::cout << "    " << std::right << std::setw(6) << iVz << " | "
 564               << std::setw(3) << bin->fLow << " - " << std::left
 565               << std::setw(3) << bin->fHigh << " ";
 566     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 567       Int_t first, last;
 568       bin->GetFirstAndLast(iIdx, first, last);
 569       std::cout << " | " << std::setw(3) << first << "-"
 570                 << std::setw(3) << last;
 571     }
 572     std::cout << std::endl;
 573   }
 574     gROOT->DecreaseDirLevel();
 575 }
 576
 577 //____________________________________________________________________
 578 AliFMDHistCollector::VtxBin::VtxBin(Int_t idx, Double_t minIpZ, Double_t maxIpZ,
 579                                     Int_t nCutBins)
 580   : fIndex(idx),
 581     fLow(minIpZ),
 582     fHigh(maxIpZ),
 583     fHitMap(0),
 584     fFirstBin(1),
 585     fLastBin(1),
 586     fNCutBins(nCutBins)
 587 {
 588 }
 589 //____________________________________________________________________
 590 AliFMDHistCollector::VtxBin::VtxBin(const VtxBin& o)
 591   : TObject(o),
 592     fIndex(o.fIndex),
 593     fLow(o.fLow),
 594     fHigh(o.fHigh),
 595     fHitMap(o.fHitMap),
 596     fFirstBin(o.fFirstBin),
 597     fLastBin(o.fLastBin),
 598     fNCutBins(o.fNCutBins)
 599 {
 600 }
 601 //____________________________________________________________________
 602 AliFMDHistCollector::VtxBin&
 603 AliFMDHistCollector::VtxBin::operator=(const VtxBin& o)
 604 {
 605   if (&o == this) return *this;
 606   fIndex    = o.fIndex;
 607   fLow      = o.fLow;
 608   fHigh     = o.fHigh;
 609   fHitMap   = o.fHitMap;
 610   fFirstBin = o.fFirstBin;
 611   fLastBin  = o.fLastBin;
 612   fNCutBins = o.fNCutBins;
 613   return *this;
 614 }
 615 //____________________________________________________________________
 616 const Char_t*
 617 AliFMDHistCollector::VtxBin::GetName() const
 618 {
 619   return Form("%c%03d_%c%03d",
 620               (fLow >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fLow)),
 621               (fHigh >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fHigh)));
 622 }
 623 //____________________________________________________________________
 624 void
 625 AliFMDHistCollector::VtxBin::SetupForData(TH2*           coverage,
 626                                           UShort_t       skips,
 627                                           FiducialMethod fiducial,
 628                                           Double_t       cut,
 629                                           TList*         l,
 630                                           const TAxis&   etaAxis,
 631                                           Bool_t         doHitMaps,
 632                                           Bool_t         storeSecMap)
 633 {
 634   TList* out = 0;
 635   if (doHitMaps || storeSecMap) {
 636     out = new TList;
 637     out->SetName(GetName());
 638     out->SetOwner();
 639     l->Add(out);
 640   }
 641   if (doHitMaps) {
 642     fHitMap = new AliForwardUtil::Histos();
 643     fHitMap->Init(etaAxis);
 644   }
 645   fFirstBin.Set(5);
 646   fLastBin.Set(5);
 647
 648   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 649
 650   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 651     UShort_t d = 0;
 652     Char_t   r = 0;
 653     GetDetRing(iIdx, d, r);
 654
 655     // Skipping selected FMD rings
 656     if (CheckSkip(d, r, skips)) continue;
 657
 658     // Get the background object
 659     TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,UShort_t(fIndex));
 660     Int_t nEta  = bg->GetNbinsX();
 661     Int_t first = nEta+1;
 662     Int_t last  = 0;
 663
 664     // Loop over the eta bins
 665     for (Int_t ie = 1; ie <= nEta; ie++) {
 666       // Loop over the phi bins to make sure that we
 667       // do not have holes in the coverage
 668       bool ok = true;
 669       for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 670         if (!CheckCorrection(fiducial, cut, bg, ie, ip)) {
 671           ok = false;
 672           continue;
 673         }
 674       }
 675       if (!ok) continue;
 676
 677       first = TMath::Min(ie, first);
 678       last  = TMath::Max(ie, last);
 679     }
 680     // Store result of first/last bin for this ring
 681     fFirstBin[iIdx] = first;
 682     fLastBin[iIdx]  = last;
 683
 684     if (fHitMap) {
 685       TH2* h = fHitMap->Get(d, r);
 686       h->SetDirectory(0);
 687       h->SetName(Form("hitMapFMD%d%c", d, r));
 688       // if (r == 'O') h->RebinY(2);
 689       out->Add(h);
 690     }
 691
 692     TH2D* obg=0;
 693     if(storeSecMap) {
 694       obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
 695       obg->SetDirectory(0);
 696       obg->Reset();
 697       out->Add(obg);
 698     }
 699     // Fill diagnostics histograms
 700     for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
 701       Double_t old = coverage->GetBinContent(ie, fIndex);
 702       coverage->SetBinContent(ie, fIndex, old+1);
 703       if(obg) {
 704         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 705           obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
 706           obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
 707         } // for (ip)
 708       } // if (doSecHits)
 709     } // for (ie)
 710   } // for (iIdx)
 711 }
 712
 713 //____________________________________________________________________
 714 void
 715 AliFMDHistCollector::VtxBin::GetFirstAndLast(Int_t  idx,
 716                                              Int_t& first,
 717                                              Int_t& last) const
 718 {
 719   // Get the first and last eta bin to use for a given ring and vertex
 720   //
 721   // Parameters:
 722   //    idx      Ring index as given by GetIdx
 723   //    first    On return, the first eta bin to use
 724   //    last     On return, the last eta bin to use
 725   //
 726   first = 0;
 727   last  = 0;
 728
 729   if (idx < 0 || idx >= fFirstBin.GetSize()) return;
 730
 731   first = fFirstBin.At(idx)+fNCutBins;
 732   last  = fLastBin.At(idx)-fNCutBins;
 733 }
 734 //____________________________________________________________________
 735 Int_t
 736 AliFMDHistCollector::VtxBin::GetFirst(Int_t  idx) const
 737 {
 738   Int_t first, last;
 739   GetFirstAndLast(idx, first , last);
 740   return first;
 741 }
 742 //____________________________________________________________________
 743 Int_t
 744 AliFMDHistCollector::VtxBin::GetLast(Int_t  idx) const
 745 {
 746   Int_t first, last;
 747   GetFirstAndLast(idx, first , last);
 748   return last;
 749 }
 750
 751 //____________________________________________________________________
 752 Bool_t
 753 AliFMDHistCollector::VtxBin::Collect(const AliForwardUtil::Histos& hists,
 754                                      AliForwardUtil::Histos&       sums,
 755                                      TH2D&                         out,
 756                                      TH2D*                         sumRings,
 757                                      TH1D*                         skipped,
 758                                      Double_t                      cent,
 759                                      MergeMethod                   m,
 760                                      UShort_t                      skips,
 761                                      TList*                        byCent,
 762                                      Bool_t                        eta2phi)
 763 {
 764   for (UShort_t d=1; d<=3; d++) {
 765     UShort_t nr = (d == 1 ? 1 : 2);
 766     for (UShort_t q=0; q<nr; q++) {
 767       Char_t      r = (q == 0 ? 'I' : 'O');
 768       Int_t       i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
 769       TH2D*       h = hists.Get(d,r);
 770       if (CheckSkip(d, r, skips) || !h ||
 771           h->TestBit(AliForwardUtil::kSkipRing)) {
 772         // Skipping a ring - either because disable, not there, or
 773         // because of flagged (too many outliers, ...)
 774         skipped->Fill(i);
 775         continue;
 776       }
 777       TH2D*       o = sums.Get(d, r);
 778       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
 779
 780       // Get valid range
 781       Int_t first = 0;
 782       Int_t last  = 0;
 783       GetFirstAndLast(d, r, first, last);
 784
 785       // Zero outside valid range
 786       Int_t nY = t->GetNbinsY();
 787       Int_t nX = t->GetNbinsX();
 788       for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) {
 789         // Lower range
 790         for (Int_t iEta = 1; iEta < first; iEta++) {
 791           t->SetBinContent(iEta,iPhi,0);
 792           t->SetBinError(iEta,iPhi,0);
 793         }
 794         for (Int_t iEta = last+1; iEta <= nX; iEta++) {
 795           t->SetBinContent(iEta,iPhi,0);
 796           t->SetBinError(iEta,iPhi,0);
 797         }
 798       }
 799       // Fill under-flow bins with eta coverage
 800       for (Int_t iEta = first; iEta <= last; iEta++)
 801         t->SetBinContent(iEta,0,1);
 802       if (eta2phi) {
 803         for (Int_t iEta = first; iEta <= last; iEta++)
 804           t->SetBinContent(iEta,nY+1,1);
 805       }
 806
 807       // Add to our per-ring sum
 808       o->Add(t);
 809
 810       // If we store hit maps, update here
 811       if (fHitMap) fHitMap->Get(d, r)->Add(t);
 812
 813       if (byCent) {
 814         TH3* dNdetaCent = static_cast<TH3*>(byCent->At(i-1));
 815         if (cent >= 0 && dNdetaCent) {
 816           Int_t iCent = dNdetaCent->GetYaxis()->FindBin(cent);
 817
 818           if (iCent > 0 && iCent <= dNdetaCent->GetNbinsY()) {
 819             // Make a projection of data
 820             TH1* proj = static_cast<TH1*>(t->ProjectionX("tmp", 1, nY));
 821             proj->SetDirectory(0);
 822             for (Int_t iEta = 1; iEta <= nX; iEta++) {
 823               Double_t v1 = proj->GetBinContent(iEta);
 824               Double_t e1 = proj->GetBinError(iEta);
 825               Double_t v2 = dNdetaCent->GetBinContent(iEta, iCent, 1);
 826               Double_t e2 = dNdetaCent->GetBinError(iEta, iCent, 1);
 827               dNdetaCent->SetBinContent(iEta,iCent,1, v1+v2);
 828               dNdetaCent->SetBinError(iEta,iCent,1, TMath::Sqrt(e1*e1+e2*e2));
 829
 830               // Check under-/overflow bins
 831               Double_t uF = t->GetBinContent(iEta, 0);
 832               Double_t oF = t->GetBinContent(iEta, nY+1);
 833               if (uF > 0) {
 834                 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 0);
 835                 dNdetaCent->SetBinContent(iEta, iCent, 0, old + uF);
 836               }
 837               if (oF > 0) {
 838                 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 2);
 839                 dNdetaCent->SetBinContent(iEta, iCent, 2, old + oF);
 840               }
 841             } // for(iEta)
 842             delete proj;
 843           } // if(iCent)
 844         } // if (cent)
 845       } // if (byCent)
 846
 847       // Outer rings have better phi segmentation - rebin to same as inner.
 848       if (q == 1) t->RebinY(2);
 849
 850       nY = t->GetNbinsY();
 851
 852       // Now update profile output
 853       for (Int_t iEta = first; iEta <= last; iEta++) {
 854
 855         // Get the possibly overlapping histogram
 856         Int_t overlap = GetOverlap(d,r,iEta);
 857
 858         // Get factor
 859         MergeMethod mm  = m; // Possibly override method locally
 860         Float_t     fac = 1;
 861         if (m != kSum && overlap >= 0) {
 862           // Default is to average
 863           fac = 0.5;
 864           if (m == kPreferInner) {
 865             // Current one is an outer overlapping an inner
 866             if ((r == 'o' || r == 'O') &&
 867                 (overlap == 0 || overlap == 1 || overlap == 3))
 868               // Do not use this signal
 869               fac = 0;
 870             // Current one is inner overlapping an outer
 871             else if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
 872               // Prefer this one
 873               fac = 1;
 874             else
 875               // In case of two overlapping inners
 876               mm = kStraightMean;
 877           }
 878           else if (m == kPreferOuter) {
 879             // Current one is an inner overlapping an outer
 880             if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
 881               // Do not use this signal
 882               fac = 0;
 883             else if ((r == 'O' || r == 'o') &&
 884                      (overlap == 0 || overlap == 1 || overlap == 3))
 885               fac = 1;
 886             else
 887               // In case of two overlapping outers
 888               mm = kStraightMean;
 889           }
 890         }
 891
 892         // Fill eta acceptance for this event into the phi underflow bin
 893         Float_t ooc      = out.GetBinContent(iEta,0);
 894         out.SetBinContent(iEta, 0, ooc + fac);
 895
 896         // Fill phi acceptance for this event into the phi overflow bin
 897         Float_t oop      = out.GetBinContent(iEta,nY+1);
 898         Float_t nop      = t->GetBinContent(iEta,nY+1);
 899 #if 0
 900         Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f",
 901              iEta, fac, oop, nop, fac*(oop+nop));
 902 #endif
 903         out.SetBinContent(iEta, nY+1, fac * nop + oop);
 904
 905         // Should we loop over h or t Y bins - I think it's t
 906         for (Int_t iPhi = 1; iPhi <= nY; iPhi++) {
 907           Double_t c  = t->GetBinContent(iEta,iPhi);
 908           Double_t e  = t->GetBinError(iEta,iPhi);
 909           Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
 910           sumRings->Fill(ee, i, c);
 911
 912           // If there's no signal or the signal was ignored because we
 913           // prefer the inners/outers, continue if (e <= 0) continue;
 914           if (fac <= 0 || c <= 0 || e <= 0)     continue;
 915
 916           // If there's no overlapping histogram (ring) or if we
 917           // prefer inner/outer, then fill in data and continue to the
 918           // next phi bin
 919           if (overlap < 0 || fac >= 1) {
 920             out.SetBinContent(iEta,iPhi,c);
 921             out.SetBinError(iEta,iPhi,e);
 922             continue;
 923           }
 924
 925           // Get the current bin content and error
 926           Double_t oc = out.GetBinContent(iEta,iPhi);
 927           Double_t oe = out.GetBinError(iEta,iPhi);
 928
 929           Double_t rc, re;
 930           MergeBins(mm, c, e, oc, oe, rc, re);
 931           out.SetBinContent(iEta,iPhi, rc);
 932           out.SetBinError(iEta,iPhi, re);
 933         }
 934       }
 935       // Remove temporary histogram
 936       delete t;
 937     } // for r
 938   } // for d
 939   return true;
 940 }
 941 //____________________________________________________________________
 942 Int_t
 943 AliFMDHistCollector::VtxBin::GetOverlap(UShort_t d, Char_t r,
 944                                         Int_t bin) const
 945 {
 946   //
 947   // Get the possibly overlapping histogram of eta bin @a e in
 948   // detector and ring
 949   //
 950   // Parameters:
 951   //    d Detector
 952   //    r Ring
 953   //    e Eta bin
 954   //    v Vertex bin (1 based)
 955   //
 956   // Return:
 957   //    Overlapping histogram index or -1
 958   //
 959
 960   Int_t other = -1;
 961   if (d == 1) {
 962     if (bin <= GetLast(2,'I')) other = GetIdx(2,'I');
 963   }
 964   else if (d == 2 && r == 'I') {
 965     if      (bin <= GetLast(2,  'O')) other = GetIdx(2, 'O');
 966     else if (bin >= GetFirst(1, 'I')) other = GetIdx(1, 'I');
 967   }
 968   else if (d == 2 && r == 'O') {
 969     if (bin >= GetFirst(2, 'I'))      other = GetIdx(2,'I');
 970   }
 971   else if (d == 3 && r == 'O') {
 972     if (bin <= GetLast(3, 'I'))       other = GetIdx(3, 'I');
 973   }
 974   else if (d == 3 && r == 'I') {
 975     if (bin >= GetFirst(3, 'O'))      other = GetIdx(3, 'O');
 976   }
 977   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
 978   return other;
 979 }
 980
 981
 982 //____________________________________________________________________
 983 //
 984 // EOF
 985 //
 986
 987
 988