More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDHistCollector.cxx
CommitLineData
0bd4b00f 1// Calculate the total energy loss
7e4038b5 2#include "AliFMDHistCollector.h"
3#include <AliESDFMD.h>
4#include <TAxis.h>
5#include <TList.h>
6#include <TMath.h>
0bd4b00f 7// #include "AliFMDAnaParameters.h"
8#include "AliForwardCorrectionManager.h"
7e4038b5 9#include "AliLog.h"
10#include <TH2D.h>
11#include <TH1I.h>
12#include <TProfile.h>
13#include <TProfile2D.h>
14#include <TArrayI.h>
0bd4b00f 15#include <TROOT.h>
16#include <iostream>
17#include <iomanip>
7e4038b5 18
19ClassImp(AliFMDHistCollector)
20#if 0
21; // For Emacs
22#endif
23
24
25//____________________________________________________________________
26AliFMDHistCollector&
27AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
28{
ea3e5d95 29 TNamed::operator=(o);
7e4038b5 30
31 fNCutBins = o.fNCutBins;
32 fCorrectionCut = o.fCorrectionCut;
33 fFirstBins = o.fFirstBins;
34 fLastBins = o.fLastBins;
ea3e5d95 35 fDebug = o.fDebug;
36
7e4038b5 37 return *this;
38}
39
40//____________________________________________________________________
41void
42AliFMDHistCollector::Init(const TAxis& vtxAxis)
43{
0bd4b00f 44 // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
45 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5 46
0bd4b00f 47 UShort_t nVz = vtxAxis.GetNbins();
7e4038b5 48 fFirstBins.Set(5*nVz);
49 fLastBins.Set(5*nVz);
50
51 // Find the eta bin ranges
0bd4b00f 52 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
7e4038b5 53 // Find the first and last eta bin to use for each ring for
54 // each vertex bin. This is instead of using the methods
55 // provided by AliFMDAnaParameters
56 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
57 UShort_t d = 0;
58 Char_t r = 0;
59 GetDetRing(iIdx, d, r);
60
61 // Get the background object
0bd4b00f 62 // TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
63 TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
7e4038b5 64 Int_t nEta = bg->GetNbinsX();
65 Int_t first = nEta+1;
66 Int_t last = 0;
67
68 // Loop over the eta bins
69 for (Int_t ie = 1; ie <= nEta; ie++) {
70 if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
71
72 // Loop over the phi bins to make sure that we
73 // do not have holes in the coverage
74 bool ok = true;
75 for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
76 if (bg->GetBinContent(ie,ip) < 0.001) {
77 ok = false;
78 continue;
79 }
80 }
81 if (!ok) continue;
82
83 first = TMath::Min(ie, first);
84 last = TMath::Max(ie, last);
85 }
86
87 // Store the result for later use
0bd4b00f 88 fFirstBins[(iVz-1)*5+iIdx] = first;
89 fLastBins[(iVz-1)*5+iIdx] = last;
7e4038b5 90 } // for j
91 }
92}
93
94//____________________________________________________________________
95Int_t
96AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
97{
98 Int_t idx = -1;
99 switch (d) {
100 case 1: idx = 0; break;
101 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
102 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
103 }
104 return idx;
105}
106//____________________________________________________________________
107void
108AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
109{
110 d = 0;
111 r = '\0';
112 switch (idx) {
113 case 0: d = 1; r = 'I'; break;
114 case 1: d = 2; r = 'I'; break;
115 case 2: d = 2; r = 'O'; break;
116 case 3: d = 3; r = 'I'; break;
117 case 4: d = 3; r = 'O'; break;
118 }
119}
120
121//____________________________________________________________________
122void
0bd4b00f 123AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
7e4038b5 124 Int_t& first, Int_t& last) const
125{
126 first = 0;
127 last = 0;
128
0bd4b00f 129 if (idx < 0) return;
130 if (vtxbin <= 0) return;
131 idx += (vtxbin-1) * 5;
7e4038b5 132
133 if (idx < 0 || idx >= fFirstBins.GetSize()) return;
134
135 first = fFirstBins.At(idx)+fNCutBins;
136 last = fLastBins.At(idx)-fNCutBins;
137}
138
139//____________________________________________________________________
140Int_t
0bd4b00f 141AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
7e4038b5 142{
143 Int_t f, l;
144 GetFirstAndLast(idx,v,f,l);
145 return f;
146}
147
148
149//____________________________________________________________________
150Int_t
0bd4b00f 151AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
7e4038b5 152{
153 Int_t f, l;
154 GetFirstAndLast(idx,v,f,l);
155 return l;
156}
157
158//____________________________________________________________________
159Int_t
160AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
0bd4b00f 161 Int_t bin, UShort_t vtxbin) const
7e4038b5 162{
163 // Get the possibly overlapping histogram
164 Int_t other = -1;
165 if (d == 1) {
166 if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
167 }
168 else if (d == 2 && r == 'I') {
169 if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
170 else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
171 }
172 else if (d == 2 && r == 'O') {
173 if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
174 }
175 else if (d == 3 && r == 'O') {
176 if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
177 }
178 else if (d == 3 && r == 'I') {
179 if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
180 }
181 // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
182 return other;
183}
184//____________________________________________________________________
185Int_t
0bd4b00f 186AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
7e4038b5 187{
188 UShort_t d = 0;
189 Char_t r = '\0';
190 GetDetRing(idx, d, r);
191 if (d == 0 || r == '\0') return 0;
192
193 return GetOverlap(d, r, bin, vtxbin);
194}
195
196
197//____________________________________________________________________
198Bool_t
199AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
0bd4b00f 200 UShort_t vtxbin,
7e4038b5 201 TH2D& out)
202{
203 for (UShort_t d=1; d<=3; d++) {
204 UShort_t nr = (d == 1 ? 1 : 2);
205 for (UShort_t q=0; q<nr; q++) {
206 Char_t r = (q == 0 ? 'I' : 'O');
207 TH2D* h = hists.Get(d,r);
208 TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
209
210 // Outer rings have better phi segmentation - rebin to same as inner.
211 if (q == 1) t->RebinY(2);
212
213 Int_t first = 0;
214 Int_t last = 0;
215 GetFirstAndLast(d, r, vtxbin, first, last);
216
217 // Now update profile output
218 for (Int_t iEta = first; iEta <= last; iEta++) {
219
220 // Get the possibly overlapping histogram
221 Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
222
223 // Fill eta acceptance for this event into the phi underlow bin
224 Float_t ooc = out.GetBinContent(iEta,0);
225 Float_t noc = overlap >= 0? 0.5 : 1;
226 out.SetBinContent(iEta, 0, ooc + noc);
227
228 for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
229 Double_t c = t->GetBinContent(iEta,iPhi);
230 Double_t e = t->GetBinError(iEta,iPhi);
231
232 // If there's no signal, continue
233 // if (e <= 0) continue;
234 if (c <= 0 || e <= 0) continue;
235
236 // If there's no overlapping histogram (ring), then
237 // fill in data and continue to the next phi bin
238 if (overlap < 0) {
239 out.SetBinContent(iEta,iPhi,c);
240 out.SetBinError(iEta,iPhi,e);
241 continue;
242 }
243
244 // Get the current bin content and error
245 Double_t oc = out.GetBinContent(iEta,iPhi);
246 Double_t oe = out.GetBinError(iEta,iPhi);
247
248#define USE_STRAIGHT_MEAN
249// #define USE_STRAIGHT_MEAN_NONZERO
250// #define USE_WEIGHTED_MEAN
251// #define USE_MOST_CERTAIN
252#if defined(USE_STRAIGHT_MEAN)
253 // calculate the average of old value (half the original),
254 // and this value, as well as the summed squared errors
255 // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
256 // and half the error of this.
257 //
258 // So, on the first overlapping histogram we get
259 //
260 // c = c_1 / 2
261 // e = sqrt((e_1 / 2)^2) = e_1/2
262 //
263 // On the second we get
264 //
265 // c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
266 // e' = sqrt(e^2 + (e_2/2)^2)
267 // = sqrt(e_1^2/4 + e_2^2/4)
268 // = sqrt(1/4 * (e_1^2+e_2^2))
269 // = 1/2 * sqrt(e_1^2 + e_2^2)
270 out.SetBinContent(iEta,iPhi,oc + c/2);
271 out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
272#elif defined(USE_STRAIGHT_MEAN_NONZERO)
273# define ZERO_OTHER
274 // If there's data in the overlapping histogram,
275 // calculate the average and add the errors in
276 // quadrature.
277 //
278 // If there's no data in the overlapping histogram,
279 // then just fill in the data
280 if (oe > 0) {
281 out.SetBinContent(iEta,iPhi,(oc + c)/2);
282 out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
283 }
284 else {
285 out.SetBinContent(iEta,iPhi,c);
286 out.SetBinError(iEta,iPhi,e);
287 }
288#elif defined(USE_WEIGHTED_MEAN)
289 // Calculate the weighted mean
290 Double_t w = 1/(e*e);
291 Double_t sc = w * c;
292 Double_t sw = w;
293 if (oe > 0) {
294 Double_t ow = 1/(oe*oe);
295 sc += ow * oc;
296 sw += ow;
297 }
298 Double_t nc = sc / sw;
299 Double_t ne = TMath::Sqrt(1 / sw);
300 out.SetBinContent(iEta,iPhi,nc,ne);
301#elif defined(USE_MOST_CERTAIN)
302# define ZERO_OTHER
303 if (e < oe) {
304 out.SetBinContent(iEta,iPhi,c);
305 out.SetBinError(iEta,iPhi,e);
306 }
307 else {
308 out.SetBinContent(iEta,iPhi,oc);
309 out.SetBinError(iEta,iPhi,oe);
310 }
311#else
312# error No method for defining content of overlapping bins defined
313#endif
314#if defined(ZERO_OTHER)
315 // Get the content of the overlapping histogram,
316 // and zero the content so that we won't use it
317 // again
318 UShort_t od; Char_t oR;
319 GetDetRing(overlap, od, oR);
320 TH2D* other = hists.Get(od,oR);
321 other->SetBinContent(iEta,iPhi,0);
322 other->SetBinError(iEta,iPhi,0);
323#endif
324 }
325 }
326 // Remove temporary histogram
327 delete t;
328 } // for r
329 } // for d
330 return kTRUE;
331}
332
0bd4b00f 333//____________________________________________________________________
334void
335AliFMDHistCollector::Print(Option_t* /* option */) const
336{
337 char ind[gROOT->GetDirLevel()+1];
338 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
339 ind[gROOT->GetDirLevel()] = '\0';
340 std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
341 << ind << " # of cut bins: " << fNCutBins << '\n'
342 << ind << " Correction cut: " << fCorrectionCut << '\n'
343 << ind << " Bin ranges:\n" << ind << " v_z bin";
344 Int_t nVz = fFirstBins.fN / 5;
345 for (UShort_t iVz = 1; iVz <= nVz; iVz++)
346 std::cout << " | " << std::setw(7) << iVz;
347 std::cout << '\n' << ind << " / ring ";
348 for (UShort_t iVz = 1; iVz <= nVz; iVz++)
349 std::cout << "-+--------";
350 std::cout << std::endl;
351
352 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
353 UShort_t d = 0;
354 Char_t r = 0;
355 GetDetRing(iIdx, d, r);
356
357 std::cout << ind << " FMD" << d << r << " ";
358 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
359 Int_t first, last;
360 GetFirstAndLast(iIdx, iVz, first, last);
361 std::cout << " | " << std::setw(3) << first << "-"
362 << std::setw(3) << last;
363 }
364 std::cout << std::endl;
365 }
366}
7e4038b5 367
368//____________________________________________________________________
369//
370// EOF
371//
372
373
374