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