e6097e1d8439bdbc9a73b51ada9e48392e477c57
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDHistCollector.cxx
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 //
15 #include "AliFMDHistCollector.h"
16 #include <AliESDFMD.h>
17 #include <TAxis.h>
18 #include <TList.h>
19 #include <TMath.h>
20 #include "AliForwardCorrectionManager.h"
21 #include "AliLog.h"
22 #include <TH2D.h>
23 #include <TH3D.h>
24 #include <TH1I.h>
25 #include <TProfile.h>
26 #include <TProfile2D.h>
27 #include <TArrayI.h>
28 #include <TROOT.h>
29 #include <iostream>
30 #include <iomanip>
31
32 ClassImp(AliFMDHistCollector)
33 #if 0
34 ; // For Emacs
35 #endif 
36
37 //____________________________________________________________________
38 AliFMDHistCollector::AliFMDHistCollector()
39   : fNCutBins(0), 
40     fCorrectionCut(0), 
41     fFirstBins(), 
42     fLastBins(), 
43     fDebug(0),
44     fList(0),
45     fSumRings(0),
46     fCoverage(0),
47     fMergeMethod(kStraightMean),
48     fFiducialMethod(kByCut),
49     fSkipFMDRings(0),
50     fBgAndHitMaps(false)
51 {
52   DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector");
53 }
54
55 //____________________________________________________________________
56 AliFMDHistCollector::AliFMDHistCollector(const char* title)
57   : TNamed("fmdHistCollector", title), 
58     fNCutBins(2), 
59     fCorrectionCut(0.5), 
60     fFirstBins(1), 
61     fLastBins(1), 
62     fDebug(0),
63     fList(0),
64     fSumRings(0),
65     fCoverage(0),
66     fMergeMethod(kStraightMean),
67     fFiducialMethod(kByCut),
68     fSkipFMDRings(0),
69     fBgAndHitMaps(false)
70 {
71   DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title);
72 }
73 //____________________________________________________________________
74 AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
75   : TNamed(o), 
76     fNCutBins(o.fNCutBins), 
77     fCorrectionCut(o.fCorrectionCut),
78     fFirstBins(o.fFirstBins), 
79     fLastBins(o.fLastBins), 
80     fDebug(o.fDebug),
81     fList(o.fList),
82     fSumRings(o.fSumRings),
83     fCoverage(o.fCoverage),
84     fMergeMethod(o.fMergeMethod),
85     fFiducialMethod(o.fFiducialMethod),
86     fSkipFMDRings(o.fSkipFMDRings),
87     fBgAndHitMaps(o.fBgAndHitMaps)
88 {
89   DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
90 }
91
92 //____________________________________________________________________
93 AliFMDHistCollector::~AliFMDHistCollector()
94
95   DGUARD(fDebug, 3, "DTOR of AliFMDHistCollector");
96   // if (fList) delete fList;
97 }
98 //____________________________________________________________________
99 AliFMDHistCollector&
100 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
101 {
102   // 
103   // Assignement operator
104   // 
105   // Parameters:
106   //    o Object to assign from 
107   //
108   // Return:
109   //    Reference to this object
110   //
111   DGUARD(fDebug, 3, "Assignment of AliFMDHistCollector");
112   if (&o == this) return *this; 
113   TNamed::operator=(o);
114
115   fNCutBins       = o.fNCutBins;
116   fCorrectionCut  = o.fCorrectionCut;
117   fFirstBins      = o.fFirstBins;
118   fLastBins       = o.fLastBins;
119   fDebug          = o.fDebug;
120   fList           = o.fList;
121   fSumRings       = o.fSumRings;
122   fCoverage       = o.fCoverage;
123   fMergeMethod    = o.fMergeMethod;
124   fFiducialMethod = o.fFiducialMethod;
125   fSkipFMDRings   = o.fSkipFMDRings;
126   fBgAndHitMaps   = o.fBgAndHitMaps;
127   
128   return *this;
129 }
130
131 //____________________________________________________________________
132 void
133 AliFMDHistCollector::SetupForData(const TAxis& vtxAxis,
134                           const TAxis& etaAxis)
135 {
136   // 
137   // Intialise 
138   // 
139   // Parameters:
140   //    vtxAxis  Vertex axis 
141   //  
142   DGUARD(fDebug, 1, "Initialization of AliFMDHistCollector");
143
144   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
145
146   fSumRings = new TH2D("sumRings", "Sum in individual rings",
147                        etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
148                        5, 1, 6);
149   fSumRings->Sumw2();
150   fSumRings->SetDirectory(0);
151   fSumRings->SetXTitle("#eta");
152   fSumRings->GetYaxis()->SetBinLabel(1,"FMD1i");
153   fSumRings->GetYaxis()->SetBinLabel(2,"FMD2i");
154   fSumRings->GetYaxis()->SetBinLabel(3,"FMD2o");
155   fSumRings->GetYaxis()->SetBinLabel(4,"FMD3i");
156   fSumRings->GetYaxis()->SetBinLabel(5,"FMD3o");
157   fList->Add(fSumRings);
158
159   fCoverage = new TH2D("coverage", "#eta coverage per v_{z}", 
160                        etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
161                        vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
162   fCoverage->SetDirectory(0);
163   fCoverage->SetXTitle("#eta");
164   fCoverage->SetYTitle("v_{z} [cm]");
165   fCoverage->SetZTitle("n_{bins}");
166   fList->Add(fCoverage);
167                        
168   UShort_t nVz = vtxAxis.GetNbins();
169   fFirstBins.Set(5*nVz);
170   fLastBins.Set(5*nVz);
171
172   // --- Add parameters to output ------------------------------------
173   fList->Add(AliForwardUtil::MakeParameter("nCutBins",fNCutBins));
174   fList->Add(AliForwardUtil::MakeParameter("skipRings",fSkipFMDRings));
175   fList->Add(AliForwardUtil::MakeParameter("bgAndHits",fBgAndHitMaps));
176   fList->Add(AliForwardUtil::MakeParameter("fiducial",Int_t(fFiducialMethod)));
177   fList->Add(AliForwardUtil::MakeParameter("fiducialCut",fCorrectionCut));
178   fList->Add(AliForwardUtil::MakeParameter("merge",Int_t(fMergeMethod)));
179
180   // Find the eta bin ranges 
181   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
182     
183     Double_t vMin    = vtxAxis.GetBinLowEdge(iVz);
184     Double_t vMax    = vtxAxis.GetBinUpEdge(iVz);
185     TList*   vtxList=0;
186     if(fBgAndHitMaps) {
187       vtxList = new TList;
188       vtxList->SetName(Form("%c%02d_%c%02d", 
189                             vMin < 0 ? 'm' : 'p', int(TMath::Abs(vMin)),
190                             vMax < 0 ? 'm' : 'p', int(TMath::Abs(vMax))));
191       fList->Add(vtxList);
192     }
193     // Find the first and last eta bin to use for each ring for 
194     // each vertex bin.   This is instead of using the methods 
195     // provided by AliFMDAnaParameters 
196     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
197       UShort_t d = 0;
198       Char_t   r = 0;
199       GetDetRing(iIdx, d, r);
200       
201       // Skipping selected FMD rings 
202       if(d==1 && r=='I' && (fSkipFMDRings & kFMD1I)) continue;
203       if(d==2 && r=='I' && (fSkipFMDRings & kFMD2I)) continue;
204       if(d==2 && r=='O' && (fSkipFMDRings & kFMD2O)) continue;
205       if(d==3 && r=='I' && (fSkipFMDRings & kFMD3I)) continue;
206       if(d==3 && r=='O' && (fSkipFMDRings & kFMD3O)) continue;
207       
208       // Get the background object 
209       // TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
210       TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
211       Int_t nEta  = bg->GetNbinsX();
212       Int_t first = nEta+1;
213       Int_t last  = 0;
214
215       // Loop over the eta bins 
216       for (Int_t ie = 1; ie <= nEta; ie++) { 
217         // Loop over the phi bins to make sure that we 
218         // do not have holes in the coverage 
219         bool ok = true;
220         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) { 
221           if (!CheckCorrection(bg, ie, ip)) {
222             ok = false;
223             continue;
224           }
225         }
226         if (!ok) continue;
227
228         first = TMath::Min(ie, first);
229         last  = TMath::Max(ie, last);
230       }
231       
232       // Store the result for later use 
233       fFirstBins[(iVz-1)*5+iIdx] = first;
234       fLastBins[(iVz-1)*5+iIdx]  = last;
235       TH2D* obg=0;
236       if(fBgAndHitMaps) {
237         obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
238         obg->SetDirectory(0);
239         obg->Reset();
240         vtxList->Add(obg);
241         
242         TH2D* hitmap = static_cast<TH2D*>(bg->Clone(Form("hitMapFMD%d%c", d, r)));
243         if(r == 'O') hitmap->RebinY(2);
244         hitmap->SetDirectory(0);
245         hitmap->GetZaxis()->SetTitle("");
246         hitmap->Reset();
247         vtxList->Add(hitmap);
248       }
249       // Fill diagnostics histograms 
250       for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
251         Double_t old = fCoverage->GetBinContent(ie, iVz);
252         fCoverage->SetBinContent(ie, iVz, old+1);
253         if(fBgAndHitMaps) {
254           for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
255             obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
256             obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
257           }
258         }
259       }
260     } // for j 
261   }
262 }
263
264 //____________________________________________________________________
265 Bool_t
266 AliFMDHistCollector::CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const
267 {
268   // 
269   // Check if we should include the bin in the data range 
270   // 
271   // Parameters:
272   //    bg Secondary map histogram
273   //    ie Eta bin
274   //    ip Phi bin
275   // 
276   // Return:
277   //    True if to be used
278   //
279   Double_t c = bg->GetBinContent(ie,ip);
280   switch (fFiducialMethod) { 
281   case kByCut:
282     return c >= fCorrectionCut;
283   case kDistance: 
284     if (2 * c < bg->GetBinContent(ie+1,ip) ||
285         2 * c < bg->GetBinContent(ie-1,ip)) return false;
286     return true;
287   default: 
288     AliError("No fiducal cut method defined");
289   }
290   return false;
291 }
292     
293 //____________________________________________________________________
294 void
295 AliFMDHistCollector::CreateOutputObjects(TList* dir)
296 {
297   // 
298   // Output diagnostic histograms to directory 
299   // 
300   // Parameters:
301   //    dir List to write in
302   //  
303   DGUARD(fDebug, 1, "Define output of AliFMDHistCollector");
304   fList = new TList;
305   fList->SetOwner();
306   fList->SetName(GetName());
307   dir->Add(fList);
308
309 }
310
311
312 //____________________________________________________________________
313 Int_t
314 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
315 {
316   // 
317   // Get the ring index from detector number and ring identifier 
318   // 
319   // Parameters:
320   //    d Detector
321   //    r Ring identifier 
322   // 
323   // Return:
324   //    ring index or -1 in case of problems 
325   //
326   Int_t idx = -1;
327   switch (d) { 
328   case 1: idx = 0; break;
329   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
330   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
331   }
332   return idx;
333 }
334 //____________________________________________________________________
335 void
336 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
337 {
338   // 
339   // Get the detector and ring from the ring index 
340   // 
341   // Parameters:
342   //    idx Ring index 
343   //    d   On return, the detector or 0 in case of errors 
344   //    r   On return, the ring id or '0' in case of errors 
345   //
346   d = 0; 
347   r = '\0';
348   switch (idx) { 
349   case 0: d = 1; r = 'I'; break;
350   case 1: d = 2; r = 'I'; break;
351   case 2: d = 2; r = 'O'; break;
352   case 3: d = 3; r = 'I'; break;
353   case 4: d = 3; r = 'O'; break;
354   }
355 }
356
357 //____________________________________________________________________
358 void
359 AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin, 
360                                      Int_t& first, Int_t& last) const
361 {
362   // 
363   // Get the first and last eta bin to use for a given ring and vertex 
364   // 
365   // Parameters:
366   //    idx      Ring index as given by GetIdx
367   //    vtxBin   Vertex bin (1 based) 
368   //    first    On return, the first eta bin to use 
369   //    last     On return, the last eta bin to use 
370   //
371   first = 0; 
372   last  = 0;
373   
374   if (idx    <  0) return;
375   if (vtxbin <= 0) return;
376   idx += (vtxbin-1) * 5;
377       
378   if (idx < 0 || idx >= fFirstBins.GetSize()) return;
379   
380   first = fFirstBins.At(idx)+fNCutBins;  
381   last  = fLastBins.At(idx)-fNCutBins;
382 }
383
384 //____________________________________________________________________
385 Int_t
386 AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const 
387 {
388   // 
389   // Get the first eta bin to use for a given ring and vertex 
390   // 
391   // Parameters:
392   //    idx Ring index as given by GetIdx
393   //    v vertex bin (1 based)
394   // 
395   // Return:
396   //    First eta bin to use, or -1 in case of problems 
397   //  
398   Int_t f, l;
399   GetFirstAndLast(idx,v,f,l);
400   return f;
401 }
402
403
404 //____________________________________________________________________
405 Int_t
406 AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const 
407 {
408   // 
409   // Get the last eta bin to use for a given ring and vertex 
410   // 
411   // Parameters:
412   //    idx Ring index as given by GetIdx
413   //    v vertex bin (1 based)
414   // 
415   // Return:
416   //    Last eta bin to use, or -1 in case of problems 
417   //  
418   Int_t f, l;
419   GetFirstAndLast(idx,v,f,l);
420   return l;
421 }
422
423 //____________________________________________________________________
424 Int_t 
425 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r, 
426                                 Int_t bin,  UShort_t vtxbin) const
427 {
428   // 
429   // Get the possibly overlapping histogram of eta bin @a e in 
430   // detector and ring 
431   // 
432   // Parameters:
433   //    d Detector
434   //    r Ring 
435   //    e Eta bin
436   //    v Vertex bin (1 based)
437   //
438   // Return:
439   //    Overlapping histogram index or -1
440   //
441
442   Int_t other = -1;
443   if (d == 1) {
444     if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
445   }
446   else if (d == 2 && r == 'I') {
447     if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
448     else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
449   }
450   else if (d == 2 && r == 'O') {
451     if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
452   }
453   else if (d == 3 && r == 'O') {
454     if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
455   }
456   else if (d == 3 && r == 'I') {
457     if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
458   }
459   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
460   return other;
461 }
462 //____________________________________________________________________
463 Int_t
464 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
465 {
466   // 
467   // Get the possibly overlapping histogram of eta bin @a e in 
468   // detector and ring 
469   // 
470   // Parameters:
471   //    i Ring index
472   //    e Eta bin
473   //    v Vertex bin (1 based)
474   //
475   // Return:
476   //    Overlapping histogram index or -1
477   //
478   UShort_t d = 0; 
479   Char_t   r = '\0';
480   GetDetRing(idx, d, r);
481   if (d == 0 || r == '\0') return 0;
482
483   return GetOverlap(d, r, bin, vtxbin);
484 }
485   
486   
487 //____________________________________________________________________
488 void
489 AliFMDHistCollector::MergeBins(Double_t c,   Double_t e, 
490                                Double_t oc,  Double_t oe,
491                                Double_t& rc, Double_t& re) const
492 {
493   // 
494   // Merge bins accoring to set method
495   // 
496   // Parameters:
497   //    c   Current content
498   //    e   Current error
499   //    oc  Old content
500   //    oe  Old error
501   //    rc  On return, the new content
502   //    re  On return, tne new error
503   //
504   rc = re = 0;
505   switch (fMergeMethod) { 
506   case kStraightMean:
507     // calculate the average of old value (half the original), 
508     // and this value, as well as the summed squared errors 
509     // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2) 
510     // and half the error of this.   
511     //
512     // So, on the first overlapping histogram we get 
513     // 
514     //    c = c_1 / 2
515     //    e = sqrt((e_1 / 2)^2) = e_1/2
516     // 
517     // On the second we get 
518     // 
519     //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2 
520     //    e' = sqrt(e^2 + (e_2/2)^2) 
521     //       = sqrt(e_1^2/4 + e_2^2/4) 
522     //       = sqrt(1/4 * (e_1^2+e_2^2)) 
523     //       = 1/2 * sqrt(e_1^2 + e_2^2)
524     rc = oc + c/2;
525     re = TMath::Sqrt(oe*oe+(e*e)/4);
526     break;
527   case kStraightMeanNoZero:
528     // If there's data in the overlapping histogram, 
529     // calculate the average and add the errors in 
530     // quadrature.  
531     // 
532     // If there's no data in the overlapping histogram, 
533     // then just fill in the data 
534     if (oe > 0) {
535       rc = (oc + c)/2;
536       re = TMath::Sqrt(oe*oe + e*e)/2;
537     }
538     else {
539       rc = c;
540       re = e;
541     }       
542     break;
543   case kWeightedMean: {
544     // Calculate the weighted mean 
545     Double_t w  = 1/(e*e);
546     Double_t sc = w * c;
547     Double_t sw = w;
548     if (oe > 0) {
549       Double_t ow =  1/(oe*oe);
550       sc          += ow * oc;
551       sw          += ow;
552     }
553     rc = sc / sw;
554     re = TMath::Sqrt(1 / sw);
555   }
556     break;
557   case kLeastError:
558     if (e < oe) {
559       rc = c;
560       re = e;
561     }
562     else {
563       rc = oc;
564       re = oe;
565     }
566     break;
567   case kSum:
568     rc = c + oc;
569     re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature 
570     break;
571   default:
572     AliError("No method for defining content of overlapping bins defined");
573     return;
574   }
575 }
576
577 //____________________________________________________________________
578 Bool_t
579 AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
580                              AliForwardUtil::Histos& sums,
581                              UShort_t                vtxbin, 
582                              TH2D&                   out,
583                              TList*                  lout,
584                              Double_t                cent,
585                              TList*      sumsv)
586 {
587   // 
588   // Do the calculations 
589   // 
590   // Parameters:
591   //    hists    Cache of histograms 
592   //    vtxBin   Vertex bin (1 based)
593   //    out      Output histogram
594   // 
595   // Return:
596   //    true on successs 
597   //
598   DGUARD(fDebug, 1, "Collect final histogram of AliFMDHistCollector");
599   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
600   const TAxis* vtxAxis = fcm.GetVertexAxis();
601   Double_t vMin    = vtxAxis->GetBinLowEdge(vtxbin);
602   Double_t vMax    = vtxAxis->GetBinUpEdge(vtxbin);
603   TList* vtxList 
604     = static_cast<TList*>(fList->FindObject(Form("%c%02d_%c%02d", 
605                                                 vMin < 0 ? 'm' : 'p', 
606                                                 int(TMath::Abs(vMin)), 
607                                                 vMax < 0 ? 'm' : 'p', 
608                                                 int(TMath::Abs(vMax)))));
609   
610   for (UShort_t d=1; d<=3; d++) { 
611     UShort_t nr = (d == 1 ? 1 : 2);
612     for (UShort_t q=0; q<nr; q++) { 
613       Char_t      r = (q == 0 ? 'I' : 'O');
614       // Skipping selected FMD rings 
615       if(d==1 && r=='I' && (fSkipFMDRings & kFMD1I)) continue; 
616       if(d==2 && r=='I' && (fSkipFMDRings & kFMD2I)) continue; 
617       if(d==2 && r=='O' && (fSkipFMDRings & kFMD2O)) continue; 
618       if(d==3 && r=='I' && (fSkipFMDRings & kFMD3I)) continue; 
619       if(d==3 && r=='O' && (fSkipFMDRings & kFMD3O)) continue; 
620
621       TH2D*       h = hists.Get(d,r);
622       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
623       Int_t       i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
624       TH2D*       o = sums.Get(d, r);
625       TH2D*       ovrt=0x0;
626       if(sumsv)
627       { 
628         AliForwardUtil::Histos* sumsvhistos=static_cast<AliForwardUtil::Histos*>(sumsv->At(vtxbin-1));
629         if(sumsvhistos)
630         {
631                 ovrt=sumsvhistos->Get(d, r);
632         }       
633       } 
634       TH3D* detavcent=0x0;
635       if(lout)
636       {
637          detavcent=static_cast<TH3D*>(lout->FindObject(Form("FMD%d%cetavcent",d,r)));         
638       } 
639  
640       // Get valid range 
641       Int_t first = 0;
642       Int_t last  = 0;
643       GetFirstAndLast(d, r, vtxbin, first, last);
644       
645       // Zero outside valid range 
646       Int_t nY = t->GetNbinsY();
647       for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) { 
648         // Lower range 
649         for (Int_t iEta = 1; iEta < first; iEta++) { 
650           t->SetBinContent(iEta,iPhi,0);
651           t->SetBinError(iEta,iPhi,0);
652         }
653         for (Int_t iEta = last+1; iEta <= t->GetNbinsX(); iEta++) {
654           t->SetBinContent(iEta,iPhi,0);
655           t->SetBinError(iEta,iPhi,0);
656         }
657       }
658       for (Int_t iEta = first; iEta <= last; iEta++) {
659         // Double_t phiAcc = t->GetBinContent(iEta, nY+1);
660         // if (phiAcc > 1e-12 && phiAcc < 1)
661         //   Info("", "FMD%d%c %3d phi acceptance: %f",d,r,iEta,phiAcc);
662         t->SetBinContent(iEta,0,1);
663       }
664       // Add to our per-ring sum 
665       o->Add(t);
666       if(ovrt)  
667         ovrt->Add(t);   
668       // fillinig the deta v cent histo
669       if(cent>0&&detavcent)
670       {
671                 Int_t nYbins=t->GetYaxis()->GetNbins();
672                 Int_t nXbins=t->GetXaxis()->GetNbins();
673                 Int_t cenbin=detavcent->GetYaxis()->FindBin(cent);
674                 if(cenbin>0&&cenbin<=detavcent->GetYaxis()->GetNbins()) 
675                 {
676                         TH1D* projectionX=(TH1D*)t->ProjectionX("tmp",1,nYbins);
677                         for (int ibineta=1;ibineta<nXbins;ibineta++) 
678                         {
679                                 Double_t v1=projectionX->GetBinContent(ibineta);
680                                 Double_t e1=projectionX->GetBinError(ibineta);
681                                 Double_t v2=detavcent->GetBinContent(ibineta,cenbin,1);
682                                 Double_t e2=detavcent->GetBinError(ibineta,cenbin,1);
683                                 detavcent->SetBinContent(ibineta,cenbin,1,v1+v2);
684                                 detavcent->SetBinError(ibineta,cenbin,1,TMath::Sqrt(e1*e1+e2*e2));
685                                 if (t->GetBinContent(ibineta,0)>0.0)
686                                         detavcent->SetBinContent(ibineta,cenbin,0,detavcent->GetBinContent(ibineta,cenbin,0)+t->GetBinContent(ibineta,0));
687                                 if (t->GetBinContent(ibineta,nYbins+1)>0.0)
688                                         detavcent->SetBinContent(ibineta,cenbin,2,detavcent->GetBinContent(ibineta,cenbin,2)+t->GetBinContent(ibineta,nYbins+1));         
689                         }
690                 }       
691       }                       
692
693
694       // Outer rings have better phi segmentation - rebin to same as inner. 
695       if (q == 1) t->RebinY(2);
696
697       nY = t->GetNbinsY();
698       // Now update profile output 
699       for (Int_t iEta = first; iEta <= last; iEta++) { 
700
701         // Get the possibly overlapping histogram 
702         Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
703
704         // Get factor 
705         Float_t fac      = (fMergeMethod != kSum && overlap >= 0 ? .5 : 1); 
706
707         // Fill eta acceptance for this event into the phi underflow bin
708         Float_t ooc      = out.GetBinContent(iEta,0);
709         out.SetBinContent(iEta, 0, ooc + fac);
710
711         // Fill phi acceptance for this event into the phi overflow bin
712         Float_t oop      = out.GetBinContent(iEta,nY+1);
713         Float_t nop      = t->GetBinContent(iEta,nY+1);
714 #if 0
715         Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f", 
716              iEta, fac, oop, nop, fac*(oop+nop));
717 #endif
718         out.SetBinContent(iEta, nY+1, fac * nop + oop);
719
720         // Should we loop over h or t Y bins - I think it's t
721         for (Int_t iPhi = 1; iPhi <= nY; iPhi++) { 
722           Double_t c  = t->GetBinContent(iEta,iPhi);
723           Double_t e  = t->GetBinError(iEta,iPhi);
724           Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
725           fSumRings->Fill(ee, i, c);
726
727           // If there's no signal, continue 
728           // if (e <= 0) continue;
729           if (c <= 0 || e <= 0)     continue;
730           
731           // If there's no overlapping histogram (ring), then 
732           // fill in data and continue to the next phi bin 
733           if (overlap < 0) { 
734             out.SetBinContent(iEta,iPhi,c);
735             out.SetBinError(iEta,iPhi,e);
736             continue;
737           }
738
739           // Get the current bin content and error 
740           Double_t oc = out.GetBinContent(iEta,iPhi);
741           Double_t oe = out.GetBinError(iEta,iPhi);
742
743           Double_t rc, re;
744           MergeBins(c, e, oc, oe, rc, re);
745           out.SetBinContent(iEta,iPhi, rc);
746           out.SetBinError(iEta,iPhi, re);
747         }
748       }
749       if(fBgAndHitMaps) {
750         TH2D* hRingSumVtx 
751           = static_cast<TH2D*>(vtxList->FindObject(Form("hitMapFMD%d%c", 
752                                                         d, r)));
753         hRingSumVtx->Add(t);
754       }
755       // Remove temporary histogram 
756       delete t;
757     } // for r
758   } // for d 
759  return kTRUE;
760 }
761
762 //____________________________________________________________________
763 void
764 AliFMDHistCollector::Print(Option_t* /* option */) const
765 {
766   // 
767   // Print information 
768   // 
769   // Parameters:
770   //    option Not used
771   //
772   char ind[gROOT->GetDirLevel()+1];
773   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
774   ind[gROOT->GetDirLevel()] = '\0';
775   std::cout << ind << ClassName() << ": " << GetName() << '\n'
776             << ind << " # of cut bins:          " << fNCutBins << '\n'
777             << ind << " Fiducal method:         " 
778             << (fFiducialMethod == kByCut ? "cut" : "distance") << "\n"
779             << ind << " Fiducial cut:           " << fCorrectionCut << "\n"
780             << ind << " Merge method:           ";
781   switch (fMergeMethod) {
782   case kStraightMean:       std::cout << "straight mean\n"; break;
783   case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
784   case kWeightedMean:       std::cout << "weighted mean\n"; break;
785   case kLeastError:         std::cout << "least error\n"; break;
786   case kSum:                std::cout << "straight sum\n"; break;
787   }
788     
789   std::cout << ind << " Bin ranges:\n" << ind << "  rings  ";
790   Int_t nVz = fFirstBins.fN / 5;
791   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
792     UShort_t d = 0;
793     Char_t   r = 0;
794     GetDetRing(iIdx, d, r);
795     std::cout << ind << " | FMD" << d << r << " ";
796   }
797   std::cout << '\n' << ind << "  /vz_bin ";
798   for (Int_t iIdx = 0; iIdx < 5; iIdx++) 
799     std::cout << "-+--------";
800   std::cout << std::endl;
801
802   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
803     std::cout << " " << std::setw(7) << iVz << "   ";
804     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
805       UShort_t d = 0;
806       Char_t   r = 0;
807       GetDetRing(iIdx, d, r);
808     
809       Int_t first, last;
810       GetFirstAndLast(iIdx, iVz, first, last);
811       std::cout << " | " << std::setw(3) << first << "-" 
812                 << std::setw(3) << last;
813     }
814     std::cout << std::endl;
815   }
816 }
817
818 //____________________________________________________________________
819 //
820 // EOF
821 //
822           
823
824