New implementation of the forward multiplicity analysis.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / older / 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 AliFMDHistCollector::AliFMDHistCollector()
22   : TNamed(), 
23     fList(), 
24     fNEvents(0),
25     fNCutBins(0),
26     fCorrectionCut(0),
27     fFirstBins(),
28     fLastBins(), 
29     fUseEtaFromData(kFALSE),
30     fEtaNorm(0), 
31     fOutput()
32 {}
33
34 //____________________________________________________________________
35 AliFMDHistCollector::AliFMDHistCollector(const char* title)
36   : TNamed("fmdHistCollector", title), 
37     fList(), 
38     fNEvents(0),
39     fNCutBins(1), 
40     fCorrectionCut(0.5),
41     fFirstBins(1),
42     fLastBins(1), 
43     fUseEtaFromData(kFALSE),
44     fEtaNorm(0), 
45     fOutput()
46 {
47   fList.SetName(GetName());
48   fOutput.SetName(GetName());
49 }
50
51 //____________________________________________________________________
52 AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
53   : TNamed(o), 
54     fList(), 
55     fNEvents(o.fNEvents),
56     fNCutBins(o.fNCutBins), 
57     fCorrectionCut(o.fCorrectionCut),
58     fFirstBins(o.fFirstBins),
59     fLastBins(o.fLastBins), 
60     fUseEtaFromData(o.fUseEtaFromData),
61     fEtaNorm(o.fEtaNorm), 
62     fOutput()
63 {
64
65   TObject* obj = 0;
66   TIter nextL(&o.fList);
67   while ((obj = nextL())) fList.Add(obj->Clone());
68   TIter nextO(&o.fOutput);
69   while ((obj = nextO())) fOutput.Add(obj->Clone());
70 }
71
72 //____________________________________________________________________
73 AliFMDHistCollector::~AliFMDHistCollector()
74 {
75   fList.Delete();
76   fOutput.Delete();
77   if (fEtaNorm) delete fEtaNorm;
78 }
79
80 //____________________________________________________________________
81 AliFMDHistCollector&
82 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
83 {
84   SetName(o.GetName());
85   SetTitle(o.GetTitle());
86   if (fEtaNorm) delete fEtaNorm;
87   fList.Delete();
88   fOutput.Delete();
89
90   fNEvents        = o.fNEvents;
91   fNCutBins       = o.fNCutBins;
92   fCorrectionCut  = o.fCorrectionCut;
93   fFirstBins      = o.fFirstBins;
94   fLastBins       = o.fLastBins;
95   fUseEtaFromData = o.fUseEtaFromData;
96
97   TObject* obj = 0;
98   TIter nextL(&o.fList);
99   while ((obj = nextL())) fList.Add(obj->Clone());
100   TIter nextO(&o.fOutput);
101   while ((obj = nextO())) fOutput.Add(obj->Clone());
102   
103   return *this;
104 }
105
106 //____________________________________________________________________
107 void
108 AliFMDHistCollector::Init(const TAxis& vtxAxis, const TAxis& etaAxis)
109 {
110   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
111
112   Int_t nVz = vtxAxis.GetNbins();
113   fFirstBins.Set(5*nVz);
114   fLastBins.Set(5*nVz);
115
116   // Find the eta bin ranges 
117   for (Int_t iVz = 0; iVz < nVz; iVz++) {
118     AliForwardUtil::Histos* h = new AliForwardUtil::Histos;
119     h->Init(etaAxis);
120     fList.AddAt(h, iVz);
121
122     // Find the first and last eta bin to use for each ring for 
123     // each vertex bin.   This is instead of using the methods 
124     // provided by AliFMDAnaParameters 
125     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
126       UShort_t d = 0;
127       Char_t   r = 0;
128       GetDetRing(iIdx, d, r);
129       
130       // Get the background object 
131       TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
132       Int_t nEta  = bg->GetNbinsX();
133       Int_t first = nEta+1;
134       Int_t last  = 0;
135
136       // Loop over the eta bins 
137       for (Int_t ie = 1; ie <= nEta; ie++) { 
138         if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
139
140         // Loop over the phi bins to make sure that we 
141         // do not have holes in the coverage 
142         bool ok = true;
143         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) { 
144           if (bg->GetBinContent(ie,ip) < 0.001) {
145             ok = false;
146             continue;
147           }
148         }
149         if (!ok) continue;
150
151         first = TMath::Min(ie, first);
152         last  = TMath::Max(ie, last);
153       }
154       
155       // Store the result for later use 
156       fFirstBins[iVz*5+iIdx] = first;
157       fLastBins[iVz*5+iIdx]  = last;
158     } // for j 
159   }
160
161   // Next, using the ranges found above, do the eta acceptance
162   // normalisation.
163   fEtaNorm = new TH1F("etaAcceptance", "Acceptance in eta", 
164                       etaAxis.GetNbins(), 
165                       etaAxis.GetXmin(), 
166                       etaAxis.GetXmax());
167   fEtaNorm->SetDirectory(0);
168   fEtaNorm->SetFillColor(kRed+1);
169   fEtaNorm->SetFillStyle(3001);
170   
171   TH1F* tmp = new TH1F("tmp", "tmp", 
172                       etaAxis.GetNbins(), 
173                       etaAxis.GetXmin(), 
174                       etaAxis.GetXmax());
175   tmp->SetDirectory(0);
176   TH2D* bg = new TH2D("bg", "Background map", 
177                       etaAxis.GetNbins(), etaAxis.GetXmin(),etaAxis.GetXmax(),
178                       20, 0, 2*TMath::Pi());
179   bg->SetDirectory(0);
180
181   Int_t colors[] = { kRed,   kPink, kMagenta, kViolet, kBlue, 
182                      kAzure, kCyan, kTeal,    kGreen,  kSpring };
183
184
185   // Loop over the vertex bins 
186   for (Int_t iVz = 0; iVz < nVz; iVz++) { 
187     // Clear cache 
188     tmp->Reset();
189     tmp->SetFillColor(colors[iVz % 10] + iVz/10);
190     tmp->SetFillStyle(3001);
191
192     bg->Reset();
193     TList* vzList = new TList;
194     vzList->SetName(Form("vtxbin%02d", iVz));
195     fOutput.AddAt(vzList, iVz);
196
197     // Loop over rings 
198     for (Int_t d = 1; d <= 3; d++) { 
199       Int_t nq = (d == 1 ? 1 : 2);
200       for (Int_t q = 0; q < nq; q++) { 
201         Char_t r = (q == 0 ? 'I' : 'O'); 
202         
203         // Get the first and last bin to use 
204         Int_t first, last;
205         GetFirstAndLast(d, r, iVz, first, last);
206         
207         TH2F* ibg = static_cast<TH2F*>(pars->GetBackgroundCorrection(d,r,iVz)
208                                        ->Clone("bg"));
209
210         // Zero out-side of selected range 
211         if (q == 1) { ibg->RebinY(2); ibg->Scale(0.5); }
212         for (Int_t iEta = 1; iEta < first; iEta++) 
213           for (Int_t iPhi = 1; iPhi <= 20; iPhi++) 
214             ibg->SetBinContent(iEta,iPhi,0);
215         for (Int_t iEta = last+1; iEta <= etaAxis.GetNbins(); iEta++) 
216           for (Int_t iPhi = 1; iPhi <= 20; iPhi++) 
217             ibg->SetBinContent(iEta,iPhi,0);
218         // Add to output 
219         bg->Add(ibg);
220         delete ibg;
221
222         // Loop over the eta bins 
223         for (Int_t iEta = first; iEta <= last; iEta++) {
224           Int_t   overlap = GetOverlap(d,r,iEta,iVz);
225           Float_t oc      = tmp->GetBinContent(iEta);
226           Float_t nc      = overlap >= 0? 0.5 : 1;
227           tmp->SetBinContent(iEta, oc + nc);
228         }
229       }
230     }
231     vzList->Add(tmp->Clone("etaAcceptance"));
232     vzList->Add(bg->Clone("secondaryMap"));
233     fEtaNorm->Add(tmp);
234   }
235   fEtaNorm->Scale(1. / nVz);
236   fOutput.Add(fEtaNorm);
237   delete tmp;
238   delete bg;
239
240 }
241
242 //____________________________________________________________________
243 Int_t
244 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
245 {
246   Int_t idx = -1;
247   switch (d) { 
248   case 1: idx = 0; break;
249   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
250   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
251   }
252   return idx;
253 }
254 //____________________________________________________________________
255 void
256 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
257 {
258   d = 0; 
259   r = '\0';
260   switch (idx) { 
261   case 0: d = 1; r = 'I'; break;
262   case 1: d = 2; r = 'I'; break;
263   case 2: d = 2; r = 'O'; break;
264   case 3: d = 3; r = 'I'; break;
265   case 4: d = 3; r = 'O'; break;
266   }
267 }
268
269 //____________________________________________________________________
270 void
271 AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin, 
272                                      Int_t& first, Int_t& last) const
273 {
274   first = 0; 
275   last  = 0;
276   
277   if (idx < 0) return;
278   idx += vtxbin * 5;
279       
280   if (idx < 0 || idx >= fFirstBins.GetSize()) return;
281   
282   first = fFirstBins.At(idx)+fNCutBins;  
283   last  = fLastBins.At(idx)-fNCutBins;
284 }
285
286 //____________________________________________________________________
287 Int_t
288 AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const 
289 {
290   Int_t f, l;
291   GetFirstAndLast(idx,v,f,l);
292   return f;
293 }
294
295
296 //____________________________________________________________________
297 Int_t
298 AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const 
299 {
300   Int_t f, l;
301   GetFirstAndLast(idx,v,f,l);
302   return l;
303 }
304
305 //____________________________________________________________________
306 Int_t 
307 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r, 
308                                 Int_t bin, Int_t vtxbin) const
309 {
310   // Get the possibly overlapping histogram 
311   Int_t other = -1;
312   if (d == 1) {
313     if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
314   }
315   else if (d == 2 && r == 'I') {
316     if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
317     else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
318   }
319   else if (d == 2 && r == 'O') {
320     if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
321   }
322   else if (d == 3 && r == 'O') {
323     if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
324   }
325   else if (d == 3 && r == 'I') {
326     if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
327   }
328   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
329   return other;
330 }
331 //____________________________________________________________________
332 Int_t
333 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
334 {
335   UShort_t d = 0; 
336   Char_t   r = '\0';
337   GetDetRing(idx, d, r);
338   if (d == 0 || r == '\0') return 0;
339
340   return GetOverlap(d, r, bin, vtxbin);
341 }
342   
343   
344 //____________________________________________________________________
345 Bool_t
346 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
347                              Int_t                   vtxbin, 
348                              TH2D&                   out)
349 {
350   if (fList.GetEntries() <= vtxbin) {
351     AliWarning(Form("Vertex bin %d out of range [0,%d]", 
352                     vtxbin, fList.GetEntries()));
353     return kFALSE;
354   }
355   
356   Merge(hists, vtxbin, out);
357   Store(hists, vtxbin);
358
359   return kTRUE;
360 }
361
362 //____________________________________________________________________
363 void
364 AliFMDHistCollector::Store(AliForwardUtil::Histos& hists,
365                            Int_t                   vtxbin)
366 {
367   AliForwardUtil::Histos* histos = 
368     static_cast<AliForwardUtil::Histos*>(fList.At(vtxbin));
369   if (!histos) { 
370     AliWarning(Form("No histogram for vertex bin %d", vtxbin));
371     return;
372   }
373   if (!histos->fFMD1i || 
374       !histos->fFMD2i || 
375       !histos->fFMD2o || 
376       !histos->fFMD3i || 
377       !histos->fFMD3o) { 
378     AliWarning(Form("Histograms for vertex bin %d not initialised "
379                     "(%p,%p,%p,%p,%p)", vtxbin, 
380                     histos->fFMD1i, histos->fFMD2i, histos->fFMD2o,
381                     histos->fFMD3i, histos->fFMD3o));
382     return;
383   }
384   if (!hists.fFMD1i || 
385       !hists.fFMD2i || 
386       !hists.fFMD2o || 
387       !hists.fFMD3i || 
388       !hists.fFMD3o) { 
389     AliWarning(Form("Cache histograms not initialised (%p,%p,%p,%p,%p)",
390                     hists.fFMD1i, hists.fFMD2i, hists.fFMD2o,
391                     hists.fFMD3i, hists.fFMD3o));
392     return;
393   }
394
395   histos->fFMD1i->Add(hists.fFMD1i);
396   histos->fFMD2i->Add(hists.fFMD2i);
397   histos->fFMD2o->Add(hists.fFMD2o);
398   histos->fFMD3i->Add(hists.fFMD3i);
399   histos->fFMD3o->Add(hists.fFMD3o);
400 }
401
402
403 //____________________________________________________________________
404 void
405 AliFMDHistCollector::Merge(AliForwardUtil::Histos& hists,
406                            Int_t                   vtxbin, 
407                            TH2D&                   out)
408 {
409   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
410   Bool_t isProfile = (out.InheritsFrom(TProfile2D::Class()));
411
412   for (UShort_t d=1; d<=3; d++) { 
413     UShort_t nr = (d == 1 ? 1 : 2);
414     for (UShort_t q=0; q<nr; q++) { 
415       Char_t      r = (q == 0 ? 'I' : 'O');
416       TH2D*       h = hists.Get(d,r);
417       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
418       
419       // Outer rings have better phi segmentation - rebin to same as inner. 
420       if (q == 1) t->RebinY(2);
421
422       Int_t first = 0;
423       Int_t last  = 0;
424       GetFirstAndLast(d, r, vtxbin, first, last);
425
426       // Now update profile output 
427       for (Int_t iEta = first; iEta <= last; iEta++) { 
428
429         // Get the possibly overlapping histogram 
430         Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
431
432         // Fill eta acceptance for this event into the phi underlow bin
433         Float_t ooc      = out.GetBinContent(iEta,0);
434         Float_t noc      = overlap >= 0? 0.5 : 1;
435         out.SetBinContent(iEta, 0, ooc + noc);
436
437         for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) { 
438           Double_t c = t->GetBinContent(iEta,iPhi);
439           Double_t e = t->GetBinError(iEta,iPhi);
440
441           // If there's no signal, continue 
442           // if (e <= 0) continue;
443           if (c <= 0 || e <= 0)     continue;
444           
445           if (isProfile) {
446             static_cast<TProfile2D&>(out).Fill(h->GetXaxis()->GetBinCenter(iEta), 
447                                                h->GetYaxis()->GetBinCenter(iPhi), 
448                                                c, e);
449             continue;
450           }
451           
452           // If there's no overlapping histogram (ring), then 
453           // fill in data and continue to the next phi bin 
454           if (overlap < 0) { 
455             out.SetBinContent(iEta,iPhi,c);
456             out.SetBinError(iEta,iPhi,e);
457             continue;
458           }
459
460           // Get the current bin content and error 
461           Double_t oc = out.GetBinContent(iEta,iPhi);
462           Double_t oe = out.GetBinError(iEta,iPhi);
463
464 #define USE_STRAIGHT_MEAN
465 // #define USE_STRAIGHT_MEAN_NONZERO
466 // #define USE_WEIGHTED_MEAN
467 // #define USE_MOST_CERTAIN
468 #if defined(USE_STRAIGHT_MEAN)
469           // calculate the average of old value (half the original), 
470           // and this value, as well as the summed squared errors 
471           // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2) 
472           // and half the error of this.   
473           //
474           // So, on the first overlapping histogram we get 
475           // 
476           //    c = c_1 / 2
477           //    e = sqrt((e_1 / 2)^2) = e_1/2
478           // 
479           // On the second we get 
480           // 
481           //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2 
482           //    e' = sqrt(e^2 + (e_2/2)^2) 
483           //       = sqrt(e_1^2/4 + e_2^2/4) 
484           //       = sqrt(1/4 * (e_1^2+e_2^2)) 
485           //       = 1/2 * sqrt(e_1^2 + e_2^2)
486           out.SetBinContent(iEta,iPhi,oc + c/2);
487           out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
488 #elif defined(USE_STRAIGHT_MEAN_NONZERO)
489 # define ZERO_OTHER
490           // If there's data in the overlapping histogram, 
491           // calculate the average and add the errors in 
492           // quadrature.  
493           // 
494           // If there's no data in the overlapping histogram, 
495           // then just fill in the data 
496           if (oe > 0) {
497             out.SetBinContent(iEta,iPhi,(oc + c)/2);
498             out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
499           }
500           else {
501             out.SetBinContent(iEta,iPhi,c);
502             out.SetBinError(iEta,iPhi,e);
503           }         
504 #elif defined(USE_WEIGHTED_MEAN) 
505           // Calculate the weighted mean 
506           Double_t w  = 1/(e*e);
507           Double_t sc = w * c;
508           Double_t sw = w;
509           if (oe > 0) {
510             Double_t ow = 1/(oe*oe);
511             sc          += ow * oc;
512             sw          += ow;
513           }
514           Double_t nc = sc / sw;
515           Double_t ne = TMath::Sqrt(1 / sw);
516           out.SetBinContent(iEta,iPhi,nc,ne);
517 #elif defined(USE_MOST_CERTAIN)
518 # define ZERO_OTHER
519           if (e < oe) {
520             out.SetBinContent(iEta,iPhi,c);
521             out.SetBinError(iEta,iPhi,e);
522           }
523           else {
524             out.SetBinContent(iEta,iPhi,oc);
525             out.SetBinError(iEta,iPhi,oe);
526           }
527 #else 
528 #         error No method for defining content of overlapping bins defined
529 #endif
530 #if defined(ZERO_OTHER)
531           // Get the content of the overlapping histogram, 
532           // and zero the content so that we won't use it 
533           // again 
534           UShort_t od; Char_t oR; 
535           GetDetRing(overlap, od, oR);
536           TH2D* other = hists.Get(od,oR);
537           other->SetBinContent(iEta,iPhi,0);
538           other->SetBinError(iEta,iPhi,0);
539 #endif
540         }
541       }
542       // Remove temporary histogram 
543       delete t;
544     } // for r
545   } // for d 
546
547   // Scale the histogram to the bin size 
548   // Int_t    nEta = out.GetNbinsX();
549   // Int_t    nPhi = out.GetNbinsY();
550   // Double_t rEta = out.GetXaxis()->GetXmax()-out.GetXaxis()->GetXmin();
551   // Double_t rPhi = out.GetYaxis()->GetXmax()-out.GetYaxis()->GetXmin();
552   // out.Scale(1. / (rEta/nEta*rPhi/nPhi));
553   // out.Scale(1., "width"); //
554 }
555
556 //____________________________________________________________________
557 void
558 AliFMDHistCollector::ScaleHistograms(const TH1I& nEvents)
559 {
560   fNEvents = &nEvents;
561
562 #if 0
563   TIter    next(&fList);
564   AliForwardUtil::Histos* histos = 0;
565   Int_t i = 1;
566   while ((histos = static_cast<AliForwardUtil::Histos*>(next()))) {
567     Int_t nev = nEvents.GetBinContent(i++);
568     if (nev <= 0) continue;
569     histos->fFMD1i->Scale(1. / nev);
570     histos->fFMD2i->Scale(1. / nev);
571     histos->fFMD2o->Scale(1. / nev);
572     histos->fFMD3i->Scale(1. / nev);
573     histos->fFMD3o->Scale(1. / nev);
574   }
575 #endif 
576 }
577
578 //____________________________________________________________________
579 void
580 AliFMDHistCollector::Output(TList* dir)
581 {
582   dir->Add(&fOutput);
583
584   // The axis 
585   TAxis* eAxis = static_cast<AliForwardUtil::Histos*>(fList.At(0))
586     ->fFMD1i->GetXaxis();
587   TAxis* vAxis = fNEvents->GetXaxis();
588   Int_t  nVtx  = vAxis->GetNbins();
589
590   // Create profiles of each vertex bin 
591   TList  mults;
592   for (Int_t i = 0; i < nVtx; i++) {
593     TProfile* mult = new TProfile(Form("dndeta_vtx%02d", i), 
594                                   Form("dN_{ch}/d#eta %f<v_{z}<%f", 
595                                        fNEvents->GetXaxis()->GetBinLowEdge(i+1),
596                                        fNEvents->GetXaxis()->GetBinUpEdge(i+1)),
597                                   eAxis->GetNbins(),eAxis->GetXmin(),
598                                   eAxis->GetXmax());
599     mult->Sumw2();
600     mults.AddAt(mult, i);
601     fOutput.Add(mult);
602   }
603   // Loop over vertex bins 
604   for (Int_t iVtx = 0; iVtx < nVtx; iVtx++) {
605     AliForwardUtil::Histos* histos = 
606       static_cast<AliForwardUtil::Histos*>(fList.At(iVtx));
607       // Output list for this vertex 
608     TList* l = 
609       static_cast<TList*>(fOutput.FindObject(Form("vtxbin%02d", iVtx)));
610
611     // Number of events in this vertex bin 
612     Int_t nev = fNEvents->GetBinContent(iVtx+1);
613     if (nev <= 0) { 
614       continue;
615     }
616     
617     // One of the outputs 
618     TProfile* mult = static_cast<TProfile*>(mults.At(iVtx));
619     if (!mult) { 
620       AliWarning(Form("No multiplicity for vertex %d", iVtx));
621       continue;
622     }
623
624     // Array of histograms and iterator 
625     TH2D* hh[] = { histos->fFMD1i, histos->fFMD2i, histos->fFMD2o,
626                    histos->fFMD3i, histos->fFMD3o, 0 };
627     TH2D** pp = hh;
628     Int_t idx = 0;
629     while (*pp) { 
630       // (eta,phi) histogram of the ring in this vertex 
631       (*pp)->SetName(Form("d2ndetadphi_%s", (*pp)->GetName()));
632       // Projection (sum) over phi of each vertex histogram 
633       TH1D* h = (*pp)->ProjectionX(Form("%s_proj", (*pp)->GetName()));
634       // Scale to number of events and bin width 
635       h->Scale(1./nev, "width");
636       // add to output 
637       l->Add(h);
638
639       Int_t firstNonZero = h->GetNbinsX();
640       Int_t lastNonZero  = 0;
641       if (!fUseEtaFromData) {
642         GetFirstAndLast(idx, iVtx, firstNonZero, lastNonZero);
643         // firstNonZero += 2;
644         // lastNonZero -= 2;
645       }
646       else {
647         // Find first and las non-zerio bins 
648         for (Int_t j = 1; j <= h->GetNbinsX(); j++) { 
649           if (h->GetBinContent(j) != 0) { 
650             firstNonZero = TMath::Min(j, firstNonZero);
651             lastNonZero  = TMath::Max(j, lastNonZero);
652           }
653         }
654         firstNonZero += fNCutBins;
655         lastNonZero  -= fNCutBins;
656       }
657
658       // Fill profile histogram 
659       UShort_t d; Char_t r;
660       GetDetRing(idx, d, r);
661       for (Int_t iEta = firstNonZero; iEta <= lastNonZero; iEta++) {
662         Double_t c = h->GetBinContent(iEta);
663         Double_t e = h->GetBinError(iEta);
664         if (e <= 0) { 
665           AliWarning(Form("error=%f in bin %d of %s/vtx=%d", 
666                           e, iEta, h->GetName(), iVtx));
667           continue;
668         }
669         mult->Fill(h->GetBinCenter(iEta), c, e);
670       }
671       idx++;
672       pp++;
673     } // while (*pp)
674
675     // Loop again to output 2D hists
676     pp = hh;
677     while (*pp) { 
678       // (*pp)->Scale(1./nev, "width");
679       l->Add(*pp);
680       pp++;
681     }
682   } // for histos 
683
684   // Result 
685   TProfile* total = new TProfile("dndeta", "#frac{1}{N}#frac{dN_{ch}}{d#eta}", 
686                                  eAxis->GetNbins(),eAxis->GetXmin(),
687                                  eAxis->GetXmax());
688
689   
690   for (Int_t iVtx = 0; iVtx < vAxis->GetNbins(); iVtx++) { 
691     TProfile* mult = static_cast<TProfile*>(mults.At(iVtx));
692     total->Add(mult);
693   }
694   fOutput.Add(total);
695 }
696
697 //____________________________________________________________________
698 //
699 // EOF
700 //
701           
702
703