AliPhysics  5eaf189 (5eaf189)
AliFMDHistCollector.cxx
Go to the documentation of this file.
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>
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 //____________________________________________________________________
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 //____________________________________________________________________
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),
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 //____________________________________________________________________
81  : TNamed(o),
82  fNCutBins(o.fNCutBins),
84  fDebug(o.fDebug),
85  fList(o.fList),
88  fSkipped(o.fSkipped),
93  fVtxList(o.fVtxList),
94  fByCent(o.fByCent),
96 {
97  DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
98 }
99 
100 //____________________________________________________________________
102 {
103  DGUARD(fDebug, 3, "DTOR of AliFMDHistCollector");
104  // if (fList) delete fList;
105 }
106 //____________________________________________________________________
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;
125  fDebug = o.fDebug;
126  fList = o.fList;
127  fSumRings = o.fSumRings;
128  fCoverage = o.fCoverage;
129  fSkipped = o.fSkipped;
134  fVtxList = o.fVtxList;
135  fByCent = o.fByCent;
136  fDoByCent = o.fDoByCent;
137  return *this;
138 }
139 
140 //____________________________________________________________________
141 void
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 ------------------------------------
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 
210  fCorrectionCut, fList, etaAxis,
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
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
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
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
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
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 //____________________________________________________________________
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 //____________________________________________________________________
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
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
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,
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
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 
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 //____________________________________________________________________
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 //____________________________________________________________________
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),
600 {
601 }
602 //____________________________________________________________________
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*
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
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) {
644  fHitMap->Init(etaAxis);
645  }
646  fFirstBin.Set(5);
647  fLastBin.Set(5);
648 
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
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
738 {
739  Int_t first, last;
740  GetFirstAndLast(idx, first , last);
741  return first;
742 }
743 //____________________________________________________________________
744 Int_t
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
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
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 
VtxBin & operator=(const VtxBin &o)
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
Definition: External.C:244
Int_t GetLast(Int_t idx) const
AliForwardUtil::Histos * fHitMap
TH2D * GetCorrection(UShort_t d, Char_t r, Double_t v) const
char Char_t
Definition: External.C:18
void SetupForData(TH2 *coverage, UShort_t skip, FiducialMethod fiducial, Double_t cut, TList *l, const TAxis &etaAxis, Bool_t doHitMap, Bool_t storeSecMap)
TList * fByCent
Per-vertex list.
TCanvas * c
Definition: TestFitELoss.C:172
Int_t GetFirst(Int_t idx) const
int Int_t
Definition: External.C:63
static Color_t RingColor(UShort_t d, Char_t r)
Int_t GetOverlap(UShort_t d, Char_t r, Int_t bin) const
float Float_t
Definition: External.C:68
const Char_t * GetName() const
Definition: External.C:252
virtual void CreateOutputObjects(TList *dir)
virtual Bool_t Collect(const AliForwardUtil::Histos &hists, AliForwardUtil::Histos &sums, UShort_t vtxBin, TH2D &out, Double_t cent=-1.0, Bool_t eta2phi=false, Bool_t add=true)
Bool_t Collect(const AliForwardUtil::Histos &hists, AliForwardUtil::Histos &sums, TH2D &out, TH2D *sumRings, TH1D *skipped, Double_t cent, MergeMethod m, UShort_t skips, TList *byCent, Bool_t eta2phi, Bool_t add)
void GetFirstAndLast(UShort_t d, UShort_t r, Int_t &first, Int_t &last) const
Definition: External.C:228
Definition: External.C:212
void Print(Option_t *option="") const
static void MergeBins(MergeMethod m, Double_t c, Double_t e, Double_t oc, Double_t oe, Double_t &rc, Double_t &re)
AliFMDHistCollector & operator=(const AliFMDHistCollector &)
#define DGUARD(L, N, F,...)
static void PrintTask(const TObject &o)
TH2D * Get(UShort_t d, Char_t r) const
static TObject * MakeParameter(const char *name, UShort_t value)
VtxBin(Int_t index=0, Double_t minIpZ=999, Double_t maxIpZ=-999, Int_t nCut=0)
#define PFV(N, VALUE)
static void GetDetRing(Int_t idx, UShort_t &d, Char_t &r)
VtxBin * GetVtxBin(Int_t ivtx)
Definition: External.C:220
void Init(const TAxis &etaAxis)
static Int_t GetIdx(UShort_t d, Char_t r)
static Bool_t CheckCorrection(FiducialMethod m, Double_t cut, const TH2D *bg, Int_t ie, Int_t ip)
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
const AliFMDCorrSecondaryMap * GetSecondaryMap() const
virtual void SetupForData(const TAxis &vtxAxis, const TAxis &etaAxis)
Definition: External.C:196
static Bool_t CheckSkip(UShort_t d, Char_t r, UShort_t skips)
static AliForwardCorrectionManager & Instance()
FiducialMethod fFiducialMethod
TDirectoryFile * dir