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