More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDHistCollector.cxx
1 // Calculate the total energy loss 
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 "AliForwardCorrectionManager.h"
9 #include "AliLog.h"
10 #include <TH2D.h>
11 #include <TH1I.h>
12 #include <TProfile.h>
13 #include <TProfile2D.h>
14 #include <TArrayI.h>
15 #include <TROOT.h>
16 #include <iostream>
17 #include <iomanip>
18
19 ClassImp(AliFMDHistCollector)
20 #if 0
21 ; // For Emacs
22 #endif 
23
24
25 //____________________________________________________________________
26 AliFMDHistCollector&
27 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
28 {
29   TNamed::operator=(o);
30
31   fNCutBins       = o.fNCutBins;
32   fCorrectionCut  = o.fCorrectionCut;
33   fFirstBins      = o.fFirstBins;
34   fLastBins       = o.fLastBins;
35   fDebug          = o.fDebug;
36
37   return *this;
38 }
39
40 //____________________________________________________________________
41 void
42 AliFMDHistCollector::Init(const TAxis& vtxAxis)
43 {
44   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
45   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
46
47   UShort_t nVz = vtxAxis.GetNbins();
48   fFirstBins.Set(5*nVz);
49   fLastBins.Set(5*nVz);
50
51   // Find the eta bin ranges 
52   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
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 
62       // TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
63       TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
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 
88       fFirstBins[(iVz-1)*5+iIdx] = first;
89       fLastBins[(iVz-1)*5+iIdx]  = last;
90     } // for j 
91   }
92 }
93
94 //____________________________________________________________________
95 Int_t
96 AliFMDHistCollector::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 //____________________________________________________________________
107 void
108 AliFMDHistCollector::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 //____________________________________________________________________
122 void
123 AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin, 
124                                      Int_t& first, Int_t& last) const
125 {
126   first = 0; 
127   last  = 0;
128   
129   if (idx    <  0) return;
130   if (vtxbin <= 0) return;
131   idx += (vtxbin-1) * 5;
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 //____________________________________________________________________
140 Int_t
141 AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const 
142 {
143   Int_t f, l;
144   GetFirstAndLast(idx,v,f,l);
145   return f;
146 }
147
148
149 //____________________________________________________________________
150 Int_t
151 AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const 
152 {
153   Int_t f, l;
154   GetFirstAndLast(idx,v,f,l);
155   return l;
156 }
157
158 //____________________________________________________________________
159 Int_t 
160 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r, 
161                                 Int_t bin,  UShort_t vtxbin) const
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 //____________________________________________________________________
185 Int_t
186 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
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 //____________________________________________________________________
198 Bool_t
199 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
200                              UShort_t                vtxbin, 
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
333 //____________________________________________________________________
334 void
335 AliFMDHistCollector::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 }
367
368 //____________________________________________________________________
369 //
370 // EOF
371 //
372           
373
374