More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDDensityCalculator.cxx
1 #include "AliFMDDensityCalculator.h"
2 #include <AliESDFMD.h>
3 #include <TAxis.h>
4 #include <TList.h>
5 #include <TMath.h>
6 #include "AliForwardCorrectionManager.h"
7 // #include "AliFMDAnaParameters.h"
8 #include "AliLog.h"
9 #include <TH2D.h>
10 #include <TProfile.h>
11 #include <TROOT.h>
12 #include <iostream>
13 #include <iomanip>
14
15 ClassImp(AliFMDDensityCalculator)
16 #if 0
17 ; // For Emacs
18 #endif 
19
20 //____________________________________________________________________
21 AliFMDDensityCalculator::AliFMDDensityCalculator()
22   : TNamed(), 
23     fRingHistos(),
24     fMultCut(0),
25     fSumOfWeights(0),
26     fWeightedSum(0),
27     fCorrections(0),
28     fMaxParticles(5),
29     fAccI(0),
30     fAccO(0),
31     fDebug(0)
32 {}
33
34 //____________________________________________________________________
35 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
36   : TNamed("fmdDensityCalculator", title), 
37     fRingHistos(), 
38     fMultCut(0),
39     fSumOfWeights(0),
40     fWeightedSum(0),
41     fCorrections(0),
42     fMaxParticles(5),
43     fAccI(0),
44     fAccO(0),
45     fDebug(0)
46 {
47   fRingHistos.SetName(GetName());
48   fRingHistos.Add(new RingHistos(1, 'I'));
49   fRingHistos.Add(new RingHistos(2, 'I'));
50   fRingHistos.Add(new RingHistos(2, 'O'));
51   fRingHistos.Add(new RingHistos(3, 'I'));
52   fRingHistos.Add(new RingHistos(3, 'O'));
53   fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
54                            200, 0, 20);
55   fSumOfWeights->SetFillColor(kRed+1);
56   fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
57   fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
58                            200, 0, 20);
59   fWeightedSum->SetFillColor(kBlue+1);
60   fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
61   fCorrections  = new TH1D("corrections", "Distribution of corrections", 
62                            100, 0, 10);
63   fCorrections->SetFillColor(kBlue+1);
64   fCorrections->SetXTitle("correction");
65
66   fAccI = GenerateAcceptanceCorrection('I');
67   fAccO = GenerateAcceptanceCorrection('O');
68 }
69
70 //____________________________________________________________________
71 AliFMDDensityCalculator::AliFMDDensityCalculator(const 
72                                                  AliFMDDensityCalculator& o)
73   : TNamed(o), 
74     fRingHistos(), 
75     fMultCut(o.fMultCut),
76     fSumOfWeights(o.fSumOfWeights),
77     fWeightedSum(o.fWeightedSum),
78     fCorrections(o.fCorrections),
79     fMaxParticles(o.fMaxParticles),
80     fAccI(o.fAccI),
81     fAccO(o.fAccO),
82     fDebug(o.fDebug)
83 {
84   TIter    next(&o.fRingHistos);
85   TObject* obj = 0;
86   while ((obj = next())) fRingHistos.Add(obj);
87 }
88
89 //____________________________________________________________________
90 AliFMDDensityCalculator::~AliFMDDensityCalculator()
91 {
92   fRingHistos.Delete();
93 }
94
95 //____________________________________________________________________
96 AliFMDDensityCalculator&
97 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
98 {
99   TNamed::operator=(o);
100
101   fMultCut      = o.fMultCut;
102   fDebug        = o.fDebug;
103   fMaxParticles = o.fMaxParticles;
104   fAccI         = o.fAccI;
105   fAccO         = o.fAccO;
106
107   fRingHistos.Delete();
108   TIter    next(&o.fRingHistos);
109   TObject* obj = 0;
110   while ((obj = next())) fRingHistos.Add(obj);
111   
112   return *this;
113 }
114
115 //____________________________________________________________________
116 AliFMDDensityCalculator::RingHistos*
117 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
118 {
119   Int_t idx = -1;
120   switch (d) { 
121   case 1: idx = 0; break;
122   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
123   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
124   }
125   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
126   
127   return static_cast<RingHistos*>(fRingHistos.At(idx));
128 }
129     
130 //____________________________________________________________________
131 Double_t
132 AliFMDDensityCalculator::GetMultCut() const
133 {
134   if (fMultCut > 0) return fMultCut;
135
136   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
137   AliFMDCorrELossFit* fits = fcm.GetELossFit();
138   return fits->GetLowCut();
139 }
140   
141 //____________________________________________________________________
142 Bool_t
143 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
144                                    AliForwardUtil::Histos& hists,
145                                    UShort_t                vtxbin, 
146                                    Bool_t                  lowFlux)
147 {
148   for (UShort_t d=1; d<=3; d++) { 
149     UShort_t nr = (d == 1 ? 1 : 2);
150     for (UShort_t q=0; q<nr; q++) { 
151       Char_t      r = (q == 0 ? 'I' : 'O');
152       UShort_t    ns= (q == 0 ?  20 :  40);
153       UShort_t    nt= (q == 0 ? 512 : 256);
154       TH2D*       h = hists.Get(d,r);
155       RingHistos* rh= GetRingHistos(d,r);
156
157       for (UShort_t s=0; s<ns; s++) { 
158         for (UShort_t t=0; t<nt; t++) {
159           Float_t mult = fmd.Multiplicity(d,r,s,t);
160           
161           if (mult == 0 || mult > 20) continue;
162
163           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
164           Float_t eta = fmd.Eta(d,r,s,t);
165           
166           Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
167           rh->fEvsN->Fill(mult,n);
168           rh->fEtaVsN->Fill(eta, n);
169
170           Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
171           fCorrections->Fill(c);
172           if (c > 0) n /= c;
173           rh->fEvsM->Fill(mult,n);
174           rh->fEtaVsM->Fill(eta, n);
175           rh->fCorr->Fill(eta, c);
176
177           h->Fill(eta,phi,n);
178           rh->fDensity->Fill(eta,phi,n);
179         } // for t
180       } // for s 
181     } // for q
182   } // for d
183   
184   return kTRUE;
185 }
186
187 //_____________________________________________________________________
188 Float_t 
189 AliFMDDensityCalculator::NParticles(Float_t  mult, 
190                                     UShort_t d, 
191                                     Char_t   r, 
192                                     UShort_t /*s*/, 
193                                     UShort_t /*t*/, 
194                                     UShort_t /*v*/, 
195                                     Float_t  eta,
196                                     Bool_t   lowFlux) const
197 {
198   if (mult <= GetMultCut()) return 0;
199   if (lowFlux) return 1;
200   
201   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
202   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
203   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
204   if (!fit) { 
205     AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
206     return 0;
207   }
208   
209   Int_t    m   = fit->FindMaxWeight();
210   if (m < 1) { 
211     AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
212     return 0;
213   }
214   UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
215   Double_t ret = fit->EvaluateWeighted(mult, n);
216
217   if (fDebug > 10) {
218     AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
219   }
220
221   fWeightedSum->Fill(ret);
222   fSumOfWeights->Fill(ret);
223
224   return ret;
225 #if 0
226   Float_t mpv  = pars->GetMPV(d,r,eta);
227   Float_t w    = pars->GetSigma(d,r,eta);
228   Float_t w2   = pars->Get2MIPWeight(d,r,eta);
229   Float_t w3   = pars->Get3MIPWeight(d,r,eta);
230   Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
231   Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
232   
233   Float_t sum  = (TMath::Landau(mult,mpv,w,kTRUE) +
234                   w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
235                   w3  * TMath::Landau(mult,mpv3,3*w,kTRUE));
236   Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
237                   2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
238                   3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
239   
240   fWeightedSum->Fill(wsum);
241   fSumOfWeights->Fill(sum);
242   return (sum > 0) ? wsum / sum : 1;
243 #endif
244 }
245
246 //_____________________________________________________________________
247 Float_t 
248 AliFMDDensityCalculator::Correction(UShort_t d, 
249                                     Char_t   r, 
250                                     UShort_t /*s*/, 
251                                     UShort_t t, 
252                                     UShort_t /*v*/, 
253                                     Float_t  eta,
254                                     Bool_t   lowFlux) const
255 {
256   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
257   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
258
259   Float_t correction = AcceptanceCorrection(r,t);
260   if (lowFlux) { 
261     // TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
262     TH1D* dblHitCor = 0;
263     if (fcm.GetDoubleHit()) 
264       dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
265
266     if (dblHitCor) {
267       Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
268       if (dblC > 0) correction *= dblC;
269     }
270     else {
271       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
272     }
273   }
274   
275 #if 0
276   TH1F* deadCor = pars->GetFMDDeadCorrection(v);
277   if (deadCor) { 
278     Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
279     if (deadC > 0) 
280       correction *= deadC; 
281   }
282 #endif  
283
284   return correction;
285 }
286
287 //_____________________________________________________________________
288 TH1D*
289 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
290 {
291   const Double_t ic1[] = { 4.9895, 15.3560 };
292   const Double_t ic2[] = { 1.8007, 17.2000 };
293   const Double_t oc1[] = { 4.2231, 26.6638 };
294   const Double_t oc2[] = { 1.8357, 27.9500 };
295   const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
296   const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
297   Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
298   Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
299   Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
300   Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
301   Float_t   basearc    = 2 * TMath::Pi() / nSec;
302   Double_t  rad        = maxR - minR;
303   Float_t   segment    = rad / nStrips;
304   Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
305
306   // Numbers used to find end-point of strip.
307   // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
308   Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
309   Float_t dx           = c2[0] - c1[0];
310   Float_t dy           = c2[1] - c1[1];
311   Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);
312
313   TH1D* ret = new TH1D(Form("acc%c", r), 
314                        Form("Acceptance correction for FMDx%c", r), 
315                        nStrips, -.5, nStrips-.5);
316   ret->SetXTitle("Strip");
317   ret->SetYTitle("#varphi acceptance");
318   ret->SetDirectory(0);
319   ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
320   ret->SetFillStyle(3001);
321
322   for (Int_t t = 0; t < nStrips; t++) { 
323     Float_t   radius     = minR + t * segment;
324     
325     // If the radius of the strip is smaller than the radius corresponding 
326     // to the first corner we have a full strip length 
327     if (radius <= cr) {
328       ret->SetBinContent(t+1, 1);
329       continue;
330     }
331
332     // Next, we should find the end-point of the strip - that is, 
333     // the coordinates where the line from c1 to c2 intersects a circle 
334     // with radius given by the strip. 
335     // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
336     // Calculate the determinant 
337     Float_t det = radius * radius * dr * dr - D*D;
338
339     if (det <=  0) { 
340       // <0 means No intersection
341       // =0 means Exactly tangent
342       ret->SetBinContent(t+1, 1);
343       continue;
344     }
345
346     // Calculate end-point and the corresponding opening angle 
347     Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
348     Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
349     Float_t th  = TMath::ATan2(x, y);
350
351     ret->SetBinContent(t+1, th / basearc);
352   }
353   return ret;
354 }
355
356 //_____________________________________________________________________
357 Float_t 
358 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
359 {
360   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
361   return acc->GetBinContent(t+1);
362
363 #if 0
364   const Double_t ic1[] = { 4.9895, 15.3560 };
365   const Double_t ic2[] = { 1.8007, 17.2000 };
366   const Double_t oc1[] = { 4.2231, 26.6638 };
367   const Double_t oc2[] = { 1.8357, 27.9500 };
368   const Double_t* c1   = (r == 'I' ? ic1      : oc1);
369   const Double_t* c2   = (r == 'I' ? ic2      : oc2);
370   Double_t  minR       = (r == 'I' ?   4.5213 :  15.4);
371   Double_t  maxR       = (r == 'I' ?  17.2    :  28.0);
372   Int_t     nStrips    = (r == 'I' ? 512      : 256);
373   Int_t     nSec       = (r == 'I' ?  20      :  40);
374   Float_t   basearc    = 2 * TMath::Pi() / nSec;
375   Double_t  rad        = maxR - minR;
376   Float_t   segment    = rad / nStrips;
377   Float_t   radius     = minR + t * segment;
378
379   // Old method - calculate full strip area and take ratio to extended
380   // strip area
381   Float_t   baselen    = basearc * radius;
382   Float_t   slope      = (c1[1] - c2[1]) / (c1[0] - c2[0]);
383   Float_t   constant   = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
384   Float_t   d          = (TMath::Power(TMath::Abs(radius*slope),2) + 
385                           TMath::Power(radius,2) - TMath::Power(constant,2));
386
387   // If below corners return 1
388   if (d >= 0) return 1;
389  
390   Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
391                          (1+TMath::Power(slope, 2)));
392   Float_t   y         = slope*x + constant;
393
394   // If x is larger than corner x or y less than corner y, we have full
395   // length strip
396   if(x >= c1[0] || y <= c1[1]) return 1;
397
398   //One sector since theta is by definition half-hybrid
399   Float_t   theta     = TMath::ATan2(x,y);
400   Float_t   arclen    = radius * theta;
401   
402   // Calculate the area of a strip with no cut
403   Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
404   Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
405   Float_t   basearea  = basearea1 - basearea2;
406
407   // Calculate the area of a strip with cut
408   Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
409   Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
410   Float_t   area      = area1 - area2;
411   
412   // Acceptance is ratio 
413   return area/basearea;
414 #endif
415 }
416
417 //____________________________________________________________________
418 void
419 AliFMDDensityCalculator::ScaleHistograms(TList* dir, Int_t nEvents)
420 {
421   if (nEvents <= 0) return;
422   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
423   if (!d) return;
424
425   TIter    next(&fRingHistos);
426   RingHistos* o = 0;
427   while ((o = static_cast<RingHistos*>(next())))
428     o->ScaleHistograms(d, nEvents);
429 }
430
431 //____________________________________________________________________
432 void
433 AliFMDDensityCalculator::DefineOutput(TList* dir)
434 {
435   TList* d = new TList;
436   d->SetName(GetName());
437   dir->Add(d);
438   d->Add(fWeightedSum);
439   d->Add(fSumOfWeights);
440   d->Add(fCorrections);
441   d->Add(fAccI);
442   d->Add(fAccO);
443
444   TIter    next(&fRingHistos);
445   RingHistos* o = 0;
446   while ((o = static_cast<RingHistos*>(next()))) {
447     o->Output(d);
448   }
449 }
450 //____________________________________________________________________
451 void
452 AliFMDDensityCalculator::Print(Option_t*) const
453 {
454   char ind[gROOT->GetDirLevel()+1];
455   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
456   ind[gROOT->GetDirLevel()] = '\0';
457   std::cout << ind << "AliFMDDensityCalculator: " << GetName() << '\n'
458             << ind << " Multiplicity cut:       " << fMultCut << '\n'
459             << ind << " Max(particles):         " << fMaxParticles 
460             << std::endl;
461 }
462
463 //====================================================================
464 AliFMDDensityCalculator::RingHistos::RingHistos()
465   : AliForwardUtil::RingHistos(),
466     fEvsN(0), 
467     fEvsM(0), 
468     fEtaVsN(0),
469     fEtaVsM(0),
470     fCorr(0),
471     fDensity(0)
472 {}
473
474 //____________________________________________________________________
475 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
476   : AliForwardUtil::RingHistos(d,r),
477     fEvsN(0), 
478     fEvsM(0),
479     fEtaVsN(0),
480     fEtaVsM(0),
481     fCorr(0),
482     fDensity(0)
483 {
484   fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), 
485                    Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
486                         "N_{ch} in %s", fName.Data()), 
487                    2500, -.5, 24.5, 2500, -.5, 24.5);
488   fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), 
489                    Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
490                         "N_{ch} in %s", fName.Data()), 
491                    2500, -.5, 24.5, 2500, -.5, 24.5);
492   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
493   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
494   fEvsN->Sumw2();
495   fEvsN->SetDirectory(0);
496   fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
497   fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
498   fEvsM->Sumw2();
499   fEvsM->SetDirectory(0);
500
501   fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
502                          Form("Average inclusive N_{ch} vs #eta (uncorrected) "
503                               "in %s", fName.Data()), 200, -4, 6);
504   fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
505                          Form("Average inclusive N_{ch} vs #eta (corrected) "
506                               "in %s", fName.Data()), 200, -4, 6);
507   fEtaVsN->SetXTitle("#eta");
508   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
509   fEtaVsN->SetDirectory(0);
510   fEtaVsN->SetLineColor(Color());
511   fEtaVsN->SetFillColor(Color());
512   fEtaVsM->SetXTitle("#eta");
513   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
514   fEtaVsM->SetDirectory(0);
515   fEtaVsM->SetLineColor(Color());
516   fEtaVsM->SetFillColor(Color());
517
518
519   fCorr = new TProfile(Form("%s_corr", fName.Data()),
520                          Form("Average correction in %s", fName.Data()), 
521                        200, -4, 6);
522   fCorr->SetXTitle("#eta");
523   fCorr->SetYTitle("#LT correction#GT");
524   fCorr->SetDirectory(0);
525   fCorr->SetLineColor(Color());
526   fCorr->SetFillColor(Color());
527
528   fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), 
529                       Form("Inclusive N_{ch} density in %s", fName.Data()), 
530                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
531                       0, 2*TMath::Pi());
532   fDensity->SetDirectory(0);
533   fDensity->SetXTitle("#eta");
534   fDensity->SetYTitle("#phi [radians]");
535   fDensity->SetZTitle("Inclusive N_{ch} density");
536 }
537 //____________________________________________________________________
538 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
539   : AliForwardUtil::RingHistos(o), 
540     fEvsN(o.fEvsN), 
541     fEvsM(o.fEvsM),
542     fEtaVsN(o.fEtaVsN),
543     fEtaVsM(o.fEtaVsM),
544     fCorr(o.fCorr),
545     fDensity(o.fDensity)
546 {}
547
548 //____________________________________________________________________
549 AliFMDDensityCalculator::RingHistos&
550 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
551 {
552   AliForwardUtil::RingHistos::operator=(o);
553   
554   if (fEvsN)    delete  fEvsN;
555   if (fEvsM)    delete  fEvsM;
556   if (fEtaVsN)  delete  fEtaVsN;
557   if (fEtaVsM)  delete  fEtaVsM;
558   if (fCorr)    delete  fCorr;
559   if (fDensity) delete fDensity;
560   
561   fEvsN    = static_cast<TH2D*>(o.fEvsN->Clone());
562   fEvsM    = static_cast<TH2D*>(o.fEvsM->Clone());
563   fEtaVsN  = static_cast<TProfile*>(o.fEtaVsN->Clone());
564   fEtaVsM  = static_cast<TProfile*>(o.fEtaVsM->Clone());
565   fCorr    = static_cast<TProfile*>(o.fCorr->Clone());
566   fDensity = static_cast<TH2D*>(o.fDensity->Clone());
567   
568   return *this;
569 }
570 //____________________________________________________________________
571 AliFMDDensityCalculator::RingHistos::~RingHistos()
572 {
573   if (fEvsN)    delete fEvsN;
574   if (fEvsM)    delete fEvsM;
575   if (fEtaVsN)  delete fEtaVsN;
576   if (fEtaVsM)  delete fEtaVsM;
577   if (fCorr)    delete fCorr;
578   if (fDensity) delete fDensity;
579 }
580
581 //____________________________________________________________________
582 void
583 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
584 {
585   TList* d = DefineOutputList(dir);
586   d->Add(fEvsN);
587   d->Add(fEvsM);
588   d->Add(fEtaVsN);
589   d->Add(fEtaVsM);
590   d->Add(fCorr);
591   d->Add(fDensity);
592 }
593
594 //____________________________________________________________________
595 void
596 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
597 {
598   TList* l = GetOutputList(dir);
599   if (!l) return; 
600
601   TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
602   if (density) density->Scale(1./nEvents);
603 }
604
605 //____________________________________________________________________
606 //
607 // EOF
608 //
609           
610
611