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