2 // This class collects the event histograms into single histograms,
3 // one for each ring in each vertex bin.
6 // - AliESDFMD object possibly corrected for sharing
9 // - 5 RingHistos objects - each with a number of vertex dependent
10 // 2D histograms of the inclusive charge particle density
12 // HistCollector used:
13 // - AliFMDCorrSecondaryMap
15 #include "AliFMDHistCollector.h"
16 #include <AliESDFMD.h>
20 #include "AliForwardCorrectionManager.h"
25 #include <TProfile2D.h>
31 ClassImp(AliFMDHistCollector)
37 //____________________________________________________________________
39 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
42 // Assignement operator
45 // o Object to assign from
48 // Reference to this object
52 fNCutBins = o.fNCutBins;
53 fCorrectionCut = o.fCorrectionCut;
54 fFirstBins = o.fFirstBins;
55 fLastBins = o.fLastBins;
61 //____________________________________________________________________
63 AliFMDHistCollector::Init(const TAxis& vtxAxis)
69 // vtxAxis Vertex axis
72 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
74 UShort_t nVz = vtxAxis.GetNbins();
75 fFirstBins.Set(5*nVz);
78 // Find the eta bin ranges
79 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
80 // Find the first and last eta bin to use for each ring for
81 // each vertex bin. This is instead of using the methods
82 // provided by AliFMDAnaParameters
83 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
86 GetDetRing(iIdx, d, r);
88 // Get the background object
89 // TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
90 TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
91 Int_t nEta = bg->GetNbinsX();
95 // Loop over the eta bins
96 for (Int_t ie = 1; ie <= nEta; ie++) {
97 if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
99 // Loop over the phi bins to make sure that we
100 // do not have holes in the coverage
102 for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
103 if (bg->GetBinContent(ie,ip) < 0.001) {
110 first = TMath::Min(ie, first);
111 last = TMath::Max(ie, last);
114 // Store the result for later use
115 fFirstBins[(iVz-1)*5+iIdx] = first;
116 fLastBins[(iVz-1)*5+iIdx] = last;
121 //____________________________________________________________________
123 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
126 // Get the ring index from detector number and ring identifier
133 // ring index or -1 in case of problems
137 case 1: idx = 0; break;
138 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
139 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
143 //____________________________________________________________________
145 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
148 // Get the detector and ring from the ring index
152 // d On return, the detector or 0 in case of errors
153 // r On return, the ring id or '0' in case of errors
158 case 0: d = 1; r = 'I'; break;
159 case 1: d = 2; r = 'I'; break;
160 case 2: d = 2; r = 'O'; break;
161 case 3: d = 3; r = 'I'; break;
162 case 4: d = 3; r = 'O'; break;
166 //____________________________________________________________________
168 AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
169 Int_t& first, Int_t& last) const
172 // Get the first and last eta bin to use for a given ring and vertex
175 // idx Ring index as given by GetIdx
176 // vtxBin Vertex bin (1 based)
177 // first On return, the first eta bin to use
178 // last On return, the last eta bin to use
184 if (vtxbin <= 0) return;
185 idx += (vtxbin-1) * 5;
187 if (idx < 0 || idx >= fFirstBins.GetSize()) return;
189 first = fFirstBins.At(idx)+fNCutBins;
190 last = fLastBins.At(idx)-fNCutBins;
193 //____________________________________________________________________
195 AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
198 // Get the first eta bin to use for a given ring and vertex
201 // idx Ring index as given by GetIdx
202 // v vertex bin (1 based)
205 // First eta bin to use, or -1 in case of problems
208 GetFirstAndLast(idx,v,f,l);
213 //____________________________________________________________________
215 AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
218 // Get the last eta bin to use for a given ring and vertex
221 // idx Ring index as given by GetIdx
222 // v vertex bin (1 based)
225 // Last eta bin to use, or -1 in case of problems
228 GetFirstAndLast(idx,v,f,l);
232 //____________________________________________________________________
234 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
235 Int_t bin, UShort_t vtxbin) const
238 // Get the possibly overlapping histogram of eta bin @a e in
245 // v Vertex bin (1 based)
248 // Overlapping histogram index or -1
253 if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
255 else if (d == 2 && r == 'I') {
256 if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
257 else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
259 else if (d == 2 && r == 'O') {
260 if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
262 else if (d == 3 && r == 'O') {
263 if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
265 else if (d == 3 && r == 'I') {
266 if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
268 // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
271 //____________________________________________________________________
273 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
276 // Get the possibly overlapping histogram of eta bin @a e in
282 // v Vertex bin (1 based)
285 // Overlapping histogram index or -1
289 GetDetRing(idx, d, r);
290 if (d == 0 || r == '\0') return 0;
292 return GetOverlap(d, r, bin, vtxbin);
296 //____________________________________________________________________
298 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
303 // Do the calculations
306 // hists Cache of histograms
307 // vtxBin Vertex bin (1 based)
308 // out Output histogram
313 for (UShort_t d=1; d<=3; d++) {
314 UShort_t nr = (d == 1 ? 1 : 2);
315 for (UShort_t q=0; q<nr; q++) {
316 Char_t r = (q == 0 ? 'I' : 'O');
317 TH2D* h = hists.Get(d,r);
318 TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
320 // Outer rings have better phi segmentation - rebin to same as inner.
321 if (q == 1) t->RebinY(2);
325 GetFirstAndLast(d, r, vtxbin, first, last);
327 // Now update profile output
328 for (Int_t iEta = first; iEta <= last; iEta++) {
330 // Get the possibly overlapping histogram
331 Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
333 // Fill eta acceptance for this event into the phi underlow bin
334 Float_t ooc = out.GetBinContent(iEta,0);
335 Float_t noc = overlap >= 0? 0.5 : 1;
336 out.SetBinContent(iEta, 0, ooc + noc);
338 for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
339 Double_t c = t->GetBinContent(iEta,iPhi);
340 Double_t e = t->GetBinError(iEta,iPhi);
342 // If there's no signal, continue
343 // if (e <= 0) continue;
344 if (c <= 0 || e <= 0) continue;
346 // If there's no overlapping histogram (ring), then
347 // fill in data and continue to the next phi bin
349 out.SetBinContent(iEta,iPhi,c);
350 out.SetBinError(iEta,iPhi,e);
354 // Get the current bin content and error
355 Double_t oc = out.GetBinContent(iEta,iPhi);
356 Double_t oe = out.GetBinError(iEta,iPhi);
358 #define USE_STRAIGHT_MEAN
359 // #define USE_STRAIGHT_MEAN_NONZERO
360 // #define USE_WEIGHTED_MEAN
361 // #define USE_MOST_CERTAIN
362 #if defined(USE_STRAIGHT_MEAN)
363 // calculate the average of old value (half the original),
364 // and this value, as well as the summed squared errors
365 // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
366 // and half the error of this.
368 // So, on the first overlapping histogram we get
371 // e = sqrt((e_1 / 2)^2) = e_1/2
373 // On the second we get
375 // c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
376 // e' = sqrt(e^2 + (e_2/2)^2)
377 // = sqrt(e_1^2/4 + e_2^2/4)
378 // = sqrt(1/4 * (e_1^2+e_2^2))
379 // = 1/2 * sqrt(e_1^2 + e_2^2)
380 out.SetBinContent(iEta,iPhi,oc + c/2);
381 out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
382 #elif defined(USE_STRAIGHT_MEAN_NONZERO)
384 // If there's data in the overlapping histogram,
385 // calculate the average and add the errors in
388 // If there's no data in the overlapping histogram,
389 // then just fill in the data
391 out.SetBinContent(iEta,iPhi,(oc + c)/2);
392 out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
395 out.SetBinContent(iEta,iPhi,c);
396 out.SetBinError(iEta,iPhi,e);
398 #elif defined(USE_WEIGHTED_MEAN)
399 // Calculate the weighted mean
400 Double_t w = 1/(e*e);
404 Double_t ow = 1/(oe*oe);
408 Double_t nc = sc / sw;
409 Double_t ne = TMath::Sqrt(1 / sw);
410 out.SetBinContent(iEta,iPhi,nc,ne);
411 #elif defined(USE_MOST_CERTAIN)
414 out.SetBinContent(iEta,iPhi,c);
415 out.SetBinError(iEta,iPhi,e);
418 out.SetBinContent(iEta,iPhi,oc);
419 out.SetBinError(iEta,iPhi,oe);
422 # error No method for defining content of overlapping bins defined
424 #if defined(ZERO_OTHER)
425 // Get the content of the overlapping histogram,
426 // and zero the content so that we won't use it
428 UShort_t od; Char_t oR;
429 GetDetRing(overlap, od, oR);
430 TH2D* other = hists.Get(od,oR);
431 other->SetBinContent(iEta,iPhi,0);
432 other->SetBinError(iEta,iPhi,0);
436 // Remove temporary histogram
443 //____________________________________________________________________
445 AliFMDHistCollector::Print(Option_t* /* option */) const
453 char ind[gROOT->GetDirLevel()+1];
454 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
455 ind[gROOT->GetDirLevel()] = '\0';
456 std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
457 << ind << " # of cut bins: " << fNCutBins << '\n'
458 << ind << " Correction cut: " << fCorrectionCut << '\n'
459 << ind << " Bin ranges:\n" << ind << " v_z bin";
460 Int_t nVz = fFirstBins.fN / 5;
461 for (UShort_t iVz = 1; iVz <= nVz; iVz++)
462 std::cout << " | " << std::setw(7) << iVz;
463 std::cout << '\n' << ind << " / ring ";
464 for (UShort_t iVz = 1; iVz <= nVz; iVz++)
465 std::cout << "-+--------";
466 std::cout << std::endl;
468 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
471 GetDetRing(iIdx, d, r);
473 std::cout << ind << " FMD" << d << r << " ";
474 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
476 GetFirstAndLast(iIdx, iVz, first, last);
477 std::cout << " | " << std::setw(3) << first << "-"
478 << std::setw(3) << last;
480 std::cout << std::endl;
484 //____________________________________________________________________