]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDHistCollector.cxx
CommitLineData
7984e5f7 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//
7e4038b5 15#include "AliFMDHistCollector.h"
16#include <AliESDFMD.h>
17#include <TAxis.h>
18#include <TList.h>
19#include <TMath.h>
0bd4b00f 20#include "AliForwardCorrectionManager.h"
8449e3e0 21#include "AliFMDCorrSecondaryMap.h"
7e4038b5 22#include "AliLog.h"
23#include <TH2D.h>
4f9319f3 24#include <TH3D.h>
7e4038b5 25#include <TH1I.h>
26#include <TProfile.h>
27#include <TProfile2D.h>
8449e3e0 28#include <TObjArray.h>
7e4038b5 29#include <TArrayI.h>
0bd4b00f 30#include <TROOT.h>
31#include <iostream>
32#include <iomanip>
7e4038b5 33
34ClassImp(AliFMDHistCollector)
35#if 0
36; // For Emacs
37#endif
38
5bb5d1f6 39//____________________________________________________________________
40AliFMDHistCollector::AliFMDHistCollector()
41 : fNCutBins(0),
42 fCorrectionCut(0),
5bb5d1f6 43 fDebug(0),
44 fList(0),
45 fSumRings(0),
46 fCoverage(0),
77f97e3f 47 fSkipped(0),
5bb5d1f6 48 fMergeMethod(kStraightMean),
ce731b9f 49 fFiducialMethod(kByCut),
50 fSkipFMDRings(0),
8449e3e0 51 fBgAndHitMaps(false),
52 fVtxList(0),
53 fByCent(0),
54 fDoByCent(false)
6ab100ec 55{
5ca83fee 56 DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector");
6ab100ec 57}
5bb5d1f6 58
59//____________________________________________________________________
60AliFMDHistCollector::AliFMDHistCollector(const char* title)
61 : TNamed("fmdHistCollector", title),
62 fNCutBins(2),
63 fCorrectionCut(0.5),
5bb5d1f6 64 fDebug(0),
65 fList(0),
66 fSumRings(0),
67 fCoverage(0),
77f97e3f 68 fSkipped(0),
5bb5d1f6 69 fMergeMethod(kStraightMean),
ce731b9f 70 fFiducialMethod(kByCut),
71 fSkipFMDRings(0),
8449e3e0 72 fBgAndHitMaps(false),
73 fVtxList(0),
74 fByCent(0),
75 fDoByCent(false)
5bb5d1f6 76{
5ca83fee 77 DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title);
5bb5d1f6 78}
79//____________________________________________________________________
80AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
81 : TNamed(o),
82 fNCutBins(o.fNCutBins),
83 fCorrectionCut(o.fCorrectionCut),
5bb5d1f6 84 fDebug(o.fDebug),
85 fList(o.fList),
86 fSumRings(o.fSumRings),
87 fCoverage(o.fCoverage),
77f97e3f 88 fSkipped(o.fSkipped),
5bb5d1f6 89 fMergeMethod(o.fMergeMethod),
ce731b9f 90 fFiducialMethod(o.fFiducialMethod),
91 fSkipFMDRings(o.fSkipFMDRings),
8449e3e0 92 fBgAndHitMaps(o.fBgAndHitMaps),
93 fVtxList(o.fVtxList),
94 fByCent(o.fByCent),
95 fDoByCent(o.fDoByCent)
6ab100ec 96{
5ca83fee 97 DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
6ab100ec 98}
5bb5d1f6 99
12fffad7 100//____________________________________________________________________
101AliFMDHistCollector::~AliFMDHistCollector()
102{
6ab100ec 103 DGUARD(fDebug, 3, "DTOR of AliFMDHistCollector");
b7ab8a2c 104 // if (fList) delete fList;
12fffad7 105}
7e4038b5 106//____________________________________________________________________
107AliFMDHistCollector&
108AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
109{
7984e5f7 110 //
111 // Assignement operator
112 //
113 // Parameters:
114 // o Object to assign from
115 //
116 // Return:
117 // Reference to this object
118 //
6ab100ec 119 DGUARD(fDebug, 3, "Assignment of AliFMDHistCollector");
d015ecfe 120 if (&o == this) return *this;
ea3e5d95 121 TNamed::operator=(o);
7e4038b5 122
123 fNCutBins = o.fNCutBins;
124 fCorrectionCut = o.fCorrectionCut;
ea3e5d95 125 fDebug = o.fDebug;
12fffad7 126 fList = o.fList;
127 fSumRings = o.fSumRings;
128 fCoverage = o.fCoverage;
77f97e3f 129 fSkipped = o.fSkipped;
e1f47419 130 fMergeMethod = o.fMergeMethod;
131 fFiducialMethod = o.fFiducialMethod;
ce731b9f 132 fSkipFMDRings = o.fSkipFMDRings;
133 fBgAndHitMaps = o.fBgAndHitMaps;
8449e3e0 134 fVtxList = o.fVtxList;
135 fByCent = o.fByCent;
136 fDoByCent = o.fDoByCent;
7e4038b5 137 return *this;
138}
139
140//____________________________________________________________________
141void
5934a3e3 142AliFMDHistCollector::SetupForData(const TAxis& vtxAxis,
8449e3e0 143 const TAxis& etaAxis)
7e4038b5 144{
7984e5f7 145 //
146 // Intialise
147 //
148 // Parameters:
149 // vtxAxis Vertex axis
150 //
6ab100ec 151 DGUARD(fDebug, 1, "Initialization of AliFMDHistCollector");
7984e5f7 152
8449e3e0 153 // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5 154
12fffad7 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}",
e1f47419 169 etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
170 vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
12fffad7 171 fCoverage->SetDirectory(0);
172 fCoverage->SetXTitle("#eta");
173 fCoverage->SetYTitle("v_{z} [cm]");
174 fCoverage->SetZTitle("n_{bins}");
175 fList->Add(fCoverage);
176
77f97e3f
CHC		 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);
7e4038b5 188
6af747ad 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
8449e3e0 197 UShort_t nVz = vtxAxis.GetNbins();
198 fVtxList = new TObjArray(nVz, 1);
199 fVtxList->SetName("histCollectorVtxBins");
200 fVtxList->SetOwner();
201
7e4038b5 202 // Find the eta bin ranges
0bd4b00f 203 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
e1f47419 204 Double_t vMin = vtxAxis.GetBinLowEdge(iVz);
205 Double_t vMax = vtxAxis.GetBinUpEdge(iVz);
8449e3e0 206 VtxBin* bin = new VtxBin(iVz, vMin, vMax, fNCutBins);
207 fVtxList->AddAt(bin, iVz);
7e4038b5 208
8449e3e0 209 bin->SetupForData(fCoverage, fSkipFMDRings, fFiducialMethod,
210 fCorrectionCut, fList, etaAxis,
211 fBgAndHitMaps, fBgAndHitMaps);
212 }
7e4038b5 213
8449e3e0 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);
7e4038b5 240 }
241}
e1f47419 242//____________________________________________________________________
243Bool_t
8449e3e0 244AliFMDHistCollector::CheckCorrection(FiducialMethod m,
245 Double_t cut,
246 const TH2D* bg,
247 Int_t ie,
248 Int_t ip)
e1f47419 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);
8449e3e0 262 switch (m) {
e1f47419 263 case kByCut:
8449e3e0 264 return c >= cut;
e1f47419 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:
8449e3e0 270 AliErrorClass("No fiducal cut method defined");
e1f47419 271 }
272 return false;
273}
274
12fffad7 275//____________________________________________________________________
276void
5934a3e3 277AliFMDHistCollector::CreateOutputObjects(TList* dir)
12fffad7 278{
279 //
280 // Output diagnostic histograms to directory
281 //
282 // Parameters:
283 // dir List to write in
284 //
6ab100ec 285 DGUARD(fDebug, 1, "Define output of AliFMDHistCollector");
12fffad7 286 fList = new TList;
287 fList->SetOwner();
288 fList->SetName(GetName());
289 dir->Add(fList);
5bb5d1f6 290
12fffad7 291}
292
8449e3e0 293//____________________________________________________________________
294Bool_t
295AliFMDHistCollector::CheckSkip(UShort_t d, Char_t r, UShort_t skips)
296{
8449e3e0 297 UShort_t q = (r == 'I' || r == 'i' ? 0 : 1);
bfab35d9 298 UShort_t c = 1 << (d-1);
299 UShort_t t = 1 << (c+q-1);
300
8449e3e0 301 return (t & skips) == t;
302}
5bb5d1f6 303
7e4038b5 304//____________________________________________________________________
305Int_t
8449e3e0 306AliFMDHistCollector::GetIdx(UShort_t d, Char_t r)
7e4038b5 307{
7984e5f7 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 //
7e4038b5 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//____________________________________________________________________
327void
8449e3e0 328AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r)
7e4038b5 329{
7984e5f7 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 //
7e4038b5 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//____________________________________________________________________
8449e3e0 350AliFMDHistCollector::VtxBin*
351AliFMDHistCollector::GetVtxBin(Int_t ivtx)
7e4038b5 352{
7984e5f7 353 // Parameters:
7984e5f7 354 // vtxBin Vertex bin (1 based)
8449e3e0 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;
7e4038b5 359}
7e4038b5 360//____________________________________________________________________
8449e3e0 361const AliFMDHistCollector::VtxBin*
362AliFMDHistCollector::GetVtxBin(Int_t ivtx) const
7e4038b5 363{
7984e5f7 364 // Parameters:
8449e3e0 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;
7e4038b5 370}
7e4038b5 371
e1f47419 372//____________________________________________________________________
373void
8449e3e0 374AliFMDHistCollector::MergeBins(MergeMethod m,
375 Double_t c, Double_t e,
e1f47419 376 Double_t oc, Double_t oe,
8449e3e0 377 Double_t& rc, Double_t& re)
e1f47419 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;
8449e3e0 391 switch (m) {
e1f47419 392 case kStraightMean:
bfab35d9 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
e1f47419 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;
5ca83fee 455 case kSum:
456 rc = c + oc;
457 re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature
458 break;
e1f47419 459 default:
8449e3e0 460 AliErrorClass("No method for defining content of overlapping bins defined");
e1f47419 461 return;
462 }
463}
464
7e4038b5 465//____________________________________________________________________
466Bool_t
5bb5d1f6 467AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
468 AliForwardUtil::Histos& sums,
0bd4b00f 469 UShort_t vtxbin,
4f9319f3 470 TH2D& out,
c8b1a7db 471 Double_t cent,
45533148 472 Bool_t eta2phi,
473 Bool_t add)
7e4038b5 474{
7984e5f7 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 //
6ab100ec 486 DGUARD(fDebug, 1, "Collect final histogram of AliFMDHistCollector");
8449e3e0 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);
717a4c2c 492 if (!bin) return false;
77f97e3f 493 Bool_t ret = bin->Collect(hists, sums, out, fSumRings, fSkipped, cent,
8449e3e0 494 fMergeMethod, fSkipFMDRings,
45533148 495 fByCent, eta2phi, add);
8449e3e0 496
497 return ret;
498}
499
c8b1a7db 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
8449e3e0 512//____________________________________________________________________
513void
514AliFMDHistCollector::Print(Option_t* /* option */) const
515{
516 //
517 // Print information
518 //
519 // Parameters:
520 // option Not used
521 //
c8b1a7db 522 TString merge("unknown");
8449e3e0 523 switch (fMergeMethod) {
c8b1a7db 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;
8449e3e0 531 }
c8b1a7db 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);
8449e3e0 539
c8b1a7db 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';
8449e3e0 547
c8b1a7db 548 std::cout << ind << "Bin ranges:\n" << ind << " rings | Range ";
8449e3e0 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 }
c8b1a7db 575 gROOT->DecreaseDirLevel();
8449e3e0 576}
577
578//____________________________________________________________________
579AliFMDHistCollector::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//____________________________________________________________________
591AliFMDHistCollector::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//____________________________________________________________________
603AliFMDHistCollector::VtxBin&
604AliFMDHistCollector::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//____________________________________________________________________
617const Char_t*
618AliFMDHistCollector::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//____________________________________________________________________
625void
626AliFMDHistCollector::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
73e536c6 649 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
73e536c6 650
8449e3e0 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
a26aa8b2 685 if (doHitMaps) {
8449e3e0 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//____________________________________________________________________
715void
716AliFMDHistCollector::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//____________________________________________________________________
736Int_t
737AliFMDHistCollector::VtxBin::GetFirst(Int_t idx) const
738{
739 Int_t first, last;
740 GetFirstAndLast(idx, first , last);
741 return first;
742}
743//____________________________________________________________________
744Int_t
745AliFMDHistCollector::VtxBin::GetLast(Int_t idx) const
746{
747 Int_t first, last;
748 GetFirstAndLast(idx, first , last);
749 return last;
750}
751
6b7e379a 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
8449e3e0 765//____________________________________________________________________
766Bool_t
767AliFMDHistCollector::VtxBin::Collect(const AliForwardUtil::Histos& hists,
768 AliForwardUtil::Histos& sums,
769 TH2D& out,
770 TH2D* sumRings,
77f97e3f 771 TH1D* skipped,
8449e3e0 772 Double_t cent,
773 MergeMethod m,
774 UShort_t skips,
c8b1a7db 775 TList* byCent,
45533148 776 Bool_t eta2phi,
777 Bool_t add)
8449e3e0 778{
7e4038b5 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');
77f97e3f 783 Int_t i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
7e4038b5 784 TH2D* h = hists.Get(d,r);
77f97e3f
CHC		 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 }
8449e3e0 792 TH2D* o = sums.Get(d, r);
7e4038b5 793 TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
8449e3e0 794
5bb5d1f6 795 // Get valid range
7e4038b5 796 Int_t first = 0;
797 Int_t last = 0;
8449e3e0 798 GetFirstAndLast(d, r, first, last);
5bb5d1f6 799
800 // Zero outside valid range
5ca83fee 801 Int_t nY = t->GetNbinsY();
8449e3e0 802 Int_t nX = t->GetNbinsX();
5ca83fee 803 for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) {
5bb5d1f6 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 }
8449e3e0 809 for (Int_t iEta = last+1; iEta <= nX; iEta++) {
5bb5d1f6 810 t->SetBinContent(iEta,iPhi,0);
811 t->SetBinError(iEta,iPhi,0);
812 }
813 }
bfab35d9 814 // Fill under-flow bins with eta coverage
815 for (Int_t iEta = first; iEta <= last; iEta++)
5bb5d1f6 816 t->SetBinContent(iEta,0,1);
c8b1a7db 817 if (eta2phi) {
818 for (Int_t iEta = first; iEta <= last; iEta++)
819 t->SetBinContent(iEta,nY+1,1);
820 }
bfab35d9 821
45533148 822 if (add) {
823 // Add to our per-ring sum
824 o->Add(t);
bfab35d9 825
3a2b980d 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);
45533148 830
4f9319f3 831
45533148 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);
8449e3e0 836
45533148 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));
8449e3e0 848
45533148 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)
4f9319f3 866
5bb5d1f6 867 // Outer rings have better phi segmentation - rebin to same as inner.
6b7e379a 868 if (q == 1) {
869 // PRINT_OVERFLOW(d, r, "before", t);
870 t->RebinY(2);
871 // PRINT_OVERFLOW(d, r, "after", t);
872 }
7e4038b5 873
5ca83fee 874 nY = t->GetNbinsY();
bfab35d9 875
7e4038b5 876 // Now update profile output
877 for (Int_t iEta = first; iEta <= last; iEta++) {
878
879 // Get the possibly overlapping histogram
8449e3e0 880 Int_t overlap = GetOverlap(d,r,iEta);
7e4038b5 881
5ca83fee 882 // Get factor
bfab35d9 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 }
5ca83fee 915
916 // Fill eta acceptance for this event into the phi underflow bin
7e4038b5 917 Float_t ooc = out.GetBinContent(iEta,0);
5ca83fee 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);
7e4038b5 928
5bb5d1f6 929 // Should we loop over h or t Y bins - I think it's t
5ca83fee 930 for (Int_t iPhi = 1; iPhi <= nY; iPhi++) {
12fffad7 931 Double_t c = t->GetBinContent(iEta,iPhi);
932 Double_t e = t->GetBinError(iEta,iPhi);
933 Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
8449e3e0 934 sumRings->Fill(ee, i, c);
7e4038b5 935
bfab35d9 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;
7e4038b5 939
bfab35d9 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) {
7e4038b5 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
e1f47419 953 Double_t rc, re;
bfab35d9 954 MergeBins(mm, c, e, oc, oe, rc, re);
e1f47419 955 out.SetBinContent(iEta,iPhi, rc);
956 out.SetBinError(iEta,iPhi, re);
7e4038b5 957 }
958 }
36ffcf83 959 // Remove temporary histogram
7e4038b5 960 delete t;
961 } // for r
962 } // for d
8449e3e0 963 return true;
7e4038b5 964}
0bd4b00f 965//____________________________________________________________________
8449e3e0 966Int_t
967AliFMDHistCollector::VtxBin::GetOverlap(UShort_t d, Char_t r,
968 Int_t bin) const
0bd4b00f 969{
7984e5f7 970 //
8449e3e0 971 // Get the possibly overlapping histogram of eta bin @a e in
972 // detector and ring
7984e5f7 973 //
974 // Parameters:
8449e3e0 975 // d Detector
976 // r Ring
977 // e Eta bin
978 // v Vertex bin (1 based)
7984e5f7 979 //
8449e3e0 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');
e1f47419 987 }
8449e3e0 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');
5bb5d1f6 991 }
8449e3e0 992 else if (d == 2 && r == 'O') {
993 if (bin >= GetFirst(2, 'I')) other = GetIdx(2,'I');
0bd4b00f 994 }
8449e3e0 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;
0bd4b00f 1003}
8449e3e0 1004
7e4038b5 1005
1006//____________________________________________________________________
1007//
1008// EOF
1009//
1010
1011
1012