Mega commit of many changes to PWGLFforward
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDDensityCalculator.cxx
1 // This class calculates the inclusive charged particle density
2 // in each for the 5 FMD rings. 
3 //
4 #include "AliFMDDensityCalculator.h"
5 #include <AliESDFMD.h>
6 #include <TAxis.h>
7 #include <TList.h>
8 #include <TMath.h>
9 #include "AliForwardCorrectionManager.h"
10 #include "AliFMDCorrDoubleHit.h"
11 #include "AliFMDCorrELossFit.h"
12 #include "AliLog.h"
13 #include <TH2D.h>
14 #include <TProfile.h>
15 #include <THStack.h>
16 #include <TROOT.h>
17 #include <TVector3.h>
18 #include <iostream>
19 #include <iomanip>
20
21 ClassImp(AliFMDDensityCalculator)
22 #if 0
23 ; // For Emacs
24 #endif 
25
26 //____________________________________________________________________
27 const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator";
28
29 //____________________________________________________________________
30 AliFMDDensityCalculator::AliFMDDensityCalculator()
31   : TNamed(), 
32     fRingHistos(),
33     fSumOfWeights(0),
34     fWeightedSum(0),
35     fCorrections(0),
36     fMaxParticles(5),
37     fUsePoisson(false),
38     fUsePhiAcceptance(kPhiCorrectNch),
39     fAccI(0),
40     fAccO(0),
41     fFMD1iMax(0),
42     fFMD2iMax(0),
43     fFMD2oMax(0),
44     fFMD3iMax(0),
45     fFMD3oMax(0),
46     fMaxWeights(0),
47     fLowCuts(0),
48     fEtaLumping(32), 
49     fPhiLumping(4),    
50     fDebug(0),
51     fCuts(),
52     fRecalculateEta(false),
53     fRecalculatePhi(false),
54     fMinQuality(10)
55 {
56   // 
57   // Constructor 
58   //
59   DGUARD(fDebug, 3, "Default CTOR of FMD density calculator");
60 }
61
62 //____________________________________________________________________
63 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
64   : TNamed(fgkFolderName, title), 
65     fRingHistos(), 
66     fSumOfWeights(0),
67     fWeightedSum(0),
68     fCorrections(0),
69     fMaxParticles(5),
70     fUsePoisson(false),
71     fUsePhiAcceptance(kPhiCorrectNch),
72     fAccI(0),
73     fAccO(0),
74     fFMD1iMax(0),
75     fFMD2iMax(0),
76     fFMD2oMax(0),
77     fFMD3iMax(0),
78     fFMD3oMax(0),
79     fMaxWeights(0),
80     fLowCuts(0),
81     fEtaLumping(32), 
82     fPhiLumping(4),
83     fDebug(0),
84     fCuts(),
85     fRecalculateEta(false),
86     fRecalculatePhi(false),
87     fMinQuality(10)
88 {
89   // 
90   // Constructor 
91   // 
92   // Parameters:
93   //    name Name of object
94   //
95   DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title);
96   fRingHistos.SetName(GetName());
97   fRingHistos.SetOwner();
98   fRingHistos.Add(new RingHistos(1, 'I'));
99   fRingHistos.Add(new RingHistos(2, 'I'));
100   fRingHistos.Add(new RingHistos(2, 'O'));
101   fRingHistos.Add(new RingHistos(3, 'I'));
102   fRingHistos.Add(new RingHistos(3, 'O'));
103   fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
104                            200, 0, 20);
105   fSumOfWeights->SetFillColor(kRed+1);
106   fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
107   fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
108                            200, 0, 20);
109   fWeightedSum->SetFillColor(kBlue+1);
110   fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
111   fCorrections  = new TH1D("corrections", "Distribution of corrections", 
112                            100, 0, 10);
113   fCorrections->SetFillColor(kBlue+1);
114   fCorrections->SetXTitle("correction");
115
116   fAccI = GenerateAcceptanceCorrection('I');
117   fAccO = GenerateAcceptanceCorrection('O');
118
119   fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use", 
120                          1, 0, 1, 1, 0, 1);
121   fMaxWeights->SetXTitle("#eta");
122   fMaxWeights->SetDirectory(0);
123
124   fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
125   fLowCuts->SetXTitle("#eta");
126   fLowCuts->SetDirectory(0);
127
128 }
129
130 //____________________________________________________________________
131 AliFMDDensityCalculator::AliFMDDensityCalculator(const 
132                                                  AliFMDDensityCalculator& o)
133   : TNamed(o), 
134     fRingHistos(), 
135     fSumOfWeights(o.fSumOfWeights),
136     fWeightedSum(o.fWeightedSum),
137     fCorrections(o.fCorrections),
138     fMaxParticles(o.fMaxParticles),
139     fUsePoisson(o.fUsePoisson),
140     fUsePhiAcceptance(o.fUsePhiAcceptance),
141     fAccI(o.fAccI),
142     fAccO(o.fAccO),
143     fFMD1iMax(o.fFMD1iMax),
144     fFMD2iMax(o.fFMD2iMax),
145     fFMD2oMax(o.fFMD2oMax),
146     fFMD3iMax(o.fFMD3iMax),
147     fFMD3oMax(o.fFMD3oMax),
148     fMaxWeights(o.fMaxWeights),
149     fLowCuts(o.fLowCuts),
150     fEtaLumping(o.fEtaLumping), 
151     fPhiLumping(o.fPhiLumping),
152     fDebug(o.fDebug),
153     fCuts(o.fCuts),
154     fRecalculateEta(o.fRecalculateEta),
155     fRecalculatePhi(o.fRecalculatePhi),
156     fMinQuality(o.fMinQuality)
157 {
158   // 
159   // Copy constructor 
160   // 
161   // Parameters:
162   //    o Object to copy from 
163   //
164   DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator");
165   TIter    next(&o.fRingHistos);
166   TObject* obj = 0;
167   while ((obj = next())) fRingHistos.Add(obj);
168 }
169
170 //____________________________________________________________________
171 AliFMDDensityCalculator::~AliFMDDensityCalculator()
172 {
173   // 
174   // Destructor 
175   //
176   DGUARD(fDebug, 3, "DTOR of FMD density calculator");
177   // fRingHistos.Delete();
178 }
179
180 //____________________________________________________________________
181 AliFMDDensityCalculator&
182 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
183 {
184   // 
185   // Assignement operator
186   // 
187   // Parameters:
188   //    o Object to assign from 
189   // 
190   // Return:
191   //    Reference to this object
192   //
193   DGUARD(fDebug, 3, "Assignment of FMD density calculator");
194   if (&o == this) return *this; 
195   TNamed::operator=(o);
196
197   fDebug              = o.fDebug;
198   fMaxParticles       = o.fMaxParticles;
199   fUsePoisson         = o.fUsePoisson;
200   fUsePhiAcceptance   = o.fUsePhiAcceptance;
201   fAccI               = o.fAccI;
202   fAccO               = o.fAccO;
203   fFMD1iMax           = o.fFMD1iMax;
204   fFMD2iMax           = o.fFMD2iMax;
205   fFMD2oMax           = o.fFMD2oMax;
206   fFMD3iMax           = o.fFMD3iMax;
207   fFMD3oMax           = o.fFMD3oMax;
208   fMaxWeights         = o.fMaxWeights;
209   fLowCuts            = o.fLowCuts;
210   fEtaLumping         = o.fEtaLumping;
211   fPhiLumping         = o.fPhiLumping;
212   fCuts               = o.fCuts;
213   fRecalculateEta     = o.fRecalculateEta;
214   fRecalculatePhi     = o.fRecalculatePhi;
215   fMinQuality         = o.fMinQuality;
216
217   fRingHistos.Delete();
218   TIter    next(&o.fRingHistos);
219   TObject* obj = 0;
220   while ((obj = next())) fRingHistos.Add(obj);
221   
222   return *this;
223 }
224
225 //____________________________________________________________________
226 void
227 AliFMDDensityCalculator::SetupForData(const TAxis& axis)
228 {
229   // Intialize this sub-algorithm 
230   //
231   // Parameters:
232   //   etaAxis   Eta axis
233   DGUARD(fDebug, 1, "Initialize FMD density calculator");
234   CacheMaxWeights(axis);
235  
236   TIter    next(&fRingHistos);
237   RingHistos* o = 0;
238   while ((o = static_cast<RingHistos*>(next()))) {
239     o->SetupForData(axis);
240     // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
241     // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
242   }
243 }
244
245 //____________________________________________________________________
246 AliFMDDensityCalculator::RingHistos*
247 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
248 {
249   // 
250   // Get the ring histogram container 
251   // 
252   // Parameters:
253   //    d Detector
254   //    r Ring 
255   // 
256   // Return:
257   //    Ring histogram container 
258   //
259   Int_t idx = -1;
260   switch (d) { 
261   case 1: idx = 0; break;
262   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
263   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
264   }
265   if (idx < 0 || idx >= fRingHistos.GetEntries()) {
266     AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
267     return 0;
268   }
269   
270   return static_cast<RingHistos*>(fRingHistos.At(idx));
271 }
272
273 //____________________________________________________________________
274 Double_t
275 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
276                                     Bool_t errors) const
277 {
278   // 
279   // Get the multiplicity cut.  If the user has set fMultCut (via
280   // SetMultCut) then that value is used.  If not, then the lower
281   // value of the fit range for the enery loss fits is returned.
282   // 
283   // Return:
284   //    Lower cut on multiplicity
285   //
286   return fCuts.GetMultCut(d,r,ieta,errors);
287 }
288     
289 //____________________________________________________________________
290 Double_t
291 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
292                                     Bool_t errors) const
293 {
294   // 
295   // Get the multiplicity cut.  If the user has set fMultCut (via
296   // SetMultCut) then that value is used.  If not, then the lower
297   // value of the fit range for the enery loss fits is returned.
298   // 
299   // Return:
300   //    Lower cut on multiplicity
301   //
302   return fCuts.GetMultCut(d,r,eta,errors);
303 }
304   
305 //____________________________________________________________________
306 Bool_t
307 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
308                                    AliForwardUtil::Histos& hists,
309                                    Bool_t                  lowFlux,
310                                    Double_t                /*cent*/, 
311                                    const TVector3&         ip)
312 {
313   // 
314   // Do the calculations 
315   // 
316   // Parameters:
317   //    fmd      AliESDFMD object (possibly) corrected for sharing
318   //    hists    Histogram cache
319   //    vtxBin   Vertex bin 
320   //    lowFlux  Low flux flag. 
321   // 
322   // Return:
323   //    true on successs 
324   DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
325
326   // --- Loop over detectors -----------------------------------------
327   for (UShort_t d=1; d<=3; d++) { 
328     UShort_t nr = (d == 1 ? 1 : 2);
329     for (UShort_t q=0; q<nr; q++) { 
330       Char_t      r = (q == 0 ? 'I' : 'O');
331       UShort_t    ns= (q == 0 ?  20 :  40);
332       UShort_t    nt= (q == 0 ? 512 : 256);
333       TH2D*       h = hists.Get(d,r);
334       RingHistos* rh= GetRingHistos(d,r);
335       if (!rh) { 
336         AliError(Form("No ring histogram found for FMD%d%c", d, r));
337         fRingHistos.ls();
338         return false;
339       }
340       // rh->fPoisson.SetObject(d,r,vtxbin,cent);
341       rh->fPoisson.Reset(0);
342       rh->fTotal->Reset();
343       rh->fGood->Reset();
344       // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
345       
346       // --- Loop over sectors and strips ----------------------------
347       for (UShort_t s=0; s<ns; s++) { 
348         for (UShort_t t=0; t<nt; t++) {
349           
350           Float_t  mult   = fmd.Multiplicity(d,r,s,t);
351           Float_t  phi    = fmd.Phi(d,r,s,t) * TMath::DegToRad();
352           Float_t  eta    = fmd.Eta(d,r,s,t);
353           Double_t oldPhi = phi;
354
355           // --- Re-calculate eta - needed for satelittes ------------
356           if( fRecalculateEta)  
357             eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
358
359           // --- Check this strip ------------------------------------
360           rh->fTotal->Fill(eta);
361           if (mult == AliESDFMD::kInvalidMult || mult > 20) {
362             // Do not count invalid stuff 
363             rh->fELoss->Fill(-1);
364             // rh->fEvsN->Fill(mult,-1);
365             // rh->fEvsM->Fill(mult,-1);
366             continue;
367           }
368           // --- Automatic calculation of acceptance -----------------
369           rh->fGood->Fill(eta);
370
371           // --- If we asked to re-calculate phi for (x,y) IP --------
372           if (fRecalculatePhi) {
373             oldPhi = phi;
374             phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
375           }
376
377           // --- Apply phi corner correction to eloss ----------------
378           if (fUsePhiAcceptance == kPhiCorrectELoss) 
379             mult *= AcceptanceCorrection(r,t);
380
381           // --- Get the low multiplicity cut ------------------------
382           Double_t cut  = 1024;
383           if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
384
385           // --- Now caluculate Nch for this strip using fits --------
386           Double_t n   = 0;
387           if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
388           rh->fELoss->Fill(mult);
389           // rh->fEvsN->Fill(mult,n);
390           // rh->fEtaVsN->Fill(eta, n);
391           
392           // --- Calculate correction if needed ----------------------
393           Double_t c = Correction(d,r,t,eta,lowFlux);
394           fCorrections->Fill(c);
395           if (c > 0) n /= c;
396           // rh->fEvsM->Fill(mult,n);
397           // rh->fEtaVsM->Fill(eta, n);
398           rh->fCorr->Fill(eta, c);
399
400           // --- Accumulate Poisson statistics -----------------------
401           Bool_t hit = (n > 0.9 && c > 0);
402           if (hit) {
403             rh->fELossUsed->Fill(mult);
404             if (fRecalculatePhi) {
405               rh->fPhiBefore->Fill(oldPhi);
406               rh->fPhiAfter->Fill(phi);
407             }
408           }
409           rh->fPoisson.Fill(t,s,hit,1./c);
410           h->Fill(eta,phi,n);
411
412           // --- If we use ELoss fits, apply now ---------------------
413           if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
414         } // for t
415       } // for s 
416
417       // --- Automatic acceptance - Calculate as an efficiency -------
418       rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
419
420       // --- Make a copy and reset as needed -------------------------
421       TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
422       if (!fUsePoisson) hclone->Reset();
423       if ( fUsePoisson) h->Reset();
424       
425       // --- Store Poisson result ------------------------------------
426       TH2D* poisson = rh->fPoisson.Result();
427       for (Int_t t=0; t < poisson->GetNbinsX(); t++) { 
428         for (Int_t s=0; s < poisson->GetNbinsY(); s++) { 
429           
430           Double_t poissonV = poisson->GetBinContent(t+1,s+1);
431           Double_t  phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
432           Double_t  eta  = fmd.Eta(d,r,s,t);
433           if(fRecalculateEta)  
434             eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
435           if (fUsePoisson)
436             h->Fill(eta,phi,poissonV);
437           else
438             hclone->Fill(eta,phi,poissonV);
439           if (fUsePoisson) rh->fDensity->Fill(eta, phi, poissonV);
440         }
441       }
442       
443       // --- Make diagnostics - eloss vs poisson ---------------------
444       Int_t nY = h->GetNbinsY();
445       for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) { 
446         // Set the overflow bin to contain the phi acceptance 
447         Double_t phiAcc  = rh->fGood->GetBinContent(ieta);
448         Double_t phiAccE = rh->fGood->GetBinError(ieta);
449         h->SetBinContent(ieta, nY+1, phiAcc);
450         h->SetBinError(ieta, nY+1, phiAccE);
451         Double_t eta     = h->GetXaxis()->GetBinCenter(ieta);
452         rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
453 #if 0
454         if (phiAcc > 1e-12) {
455           Info("", "FMD%d%c, eta=%3d phi acceptance: %f (%f)", 
456                d, r, ieta, phiAcc, h->GetBinContent(ieta, nY+1));
457         }
458 #endif
459         for (Int_t iphi=1; iphi<= nY; iphi++) { 
460           
461           Double_t poissonV =  0; //h->GetBinContent(,s+1);
462           Double_t eLossV =  0;
463           if(fUsePoisson) { 
464             poissonV = h->GetBinContent(ieta,iphi);
465             eLossV  = hclone->GetBinContent(ieta,iphi);
466           }
467           else { 
468             poissonV = hclone->GetBinContent(ieta,iphi);
469             eLossV  = h->GetBinContent(ieta,iphi);
470           }
471           
472           rh->fELossVsPoisson->Fill(eLossV, poissonV);
473           rh->fDiffELossPoisson->Fill(poissonV < 1e-12 ? 0 : 
474                                       (eLossV - poissonV) / poissonV);
475                                       
476         }
477       }
478       delete hclone;
479       
480     } // for q
481   } // for d
482
483   return kTRUE;
484 }
485
486 //_____________________________________________________________________
487 Int_t
488 AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
489                                        UShort_t d, Char_t r, Int_t iEta) const
490 {
491   // 
492   // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
493   // 
494   // Parameters:
495   //    cor   Correction
496   //    d     Detector 
497   //    r     Ring 
498   //    iEta  Eta bin 
499   //
500   DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
501   if(!cor) return -1; 
502
503   AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,iEta, -1);
504   if (!fit) { 
505     // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
506     return -1;
507   }
508   return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError, 
509                             AliFMDCorrELossFit::ELossFit::fgLeastWeight, 
510                             fMaxParticles);
511 }
512   
513 //_____________________________________________________________________
514 void
515 AliFMDDensityCalculator::CacheMaxWeights(const TAxis& axis)
516 {
517   // 
518   // Find the max weights and cache them 
519   // 
520   DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
521   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
522   const AliFMDCorrELossFit*     cor = fcm.GetELossFit();
523   cor->CacheBins(fMinQuality);
524
525   TAxis eta(axis.GetNbins(),
526             axis.GetXmin(),
527             axis.GetXmax());
528
529   if(cor) 
530     eta.Set(cor->GetEtaAxis().GetNbins(), 
531             cor->GetEtaAxis().GetXmin(), 
532             cor->GetEtaAxis().GetXmax());
533
534   Int_t nEta = eta.GetNbins();
535   fFMD1iMax.Set(nEta);
536   fFMD2iMax.Set(nEta);
537   fFMD2oMax.Set(nEta);
538   fFMD3iMax.Set(nEta);
539   fFMD3oMax.Set(nEta);
540   
541   fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
542   fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
543   fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
544   fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
545   fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
546   fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
547
548   AliInfo(Form("Get eta axis with %d bins from %f to %f",
549                nEta, eta.GetXmin(), eta.GetXmax()));
550   fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
551   fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
552   fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
553   fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
554   fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
555   fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
556   
557   for (Int_t i = 0; i < nEta; i++) {
558     Double_t w[5];
559     w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
560     w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
561     w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
562     w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
563     w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
564     Double_t l[5];
565     l[0] = GetMultCut(1, 'I', i+1, false);
566     l[1] = GetMultCut(2, 'I', i+1, false);
567     l[2] = GetMultCut(2, 'O', i+1, false);
568     l[3] = GetMultCut(3, 'I', i+1, false);
569     l[4] = GetMultCut(3, 'O', i+1, false);
570     for (Int_t j = 0; j < 5; j++) { 
571       if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
572       if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
573     }
574   }
575 }
576
577 //_____________________________________________________________________
578 Int_t
579 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
580 {
581   // 
582   // Find the (cached) maximum weight for FMD<i>dr</i> in 
583   // @f$\eta@f$ bin @a iEta
584   // 
585   // Parameters:
586   //    d     Detector
587   //    r     Ring
588   //    iEta  Eta bin
589   // 
590   // Return:
591   //    max weight or <= 0 in case of problems 
592   //
593   if (iEta < 0) return -1;
594
595   const TArrayI* max  = 0;
596   switch (d) { 
597   case 1:  max = &fFMD1iMax;                                       break;
598   case 2:  max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
599   case 3:  max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
600   }
601   if (!max) { 
602     AliWarning(Form("No array for FMD%d%c", d, r));
603     return -1;
604   }
605   
606   if (iEta >= max->fN) { 
607     AliWarning(Form("Eta bin %3d out of bounds [0,%d]", 
608                     iEta, max->fN-1));
609     return -1;
610   }
611
612   AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta, 
613                    max->At(iEta)));
614   return max->At(iEta);
615 }
616
617 //_____________________________________________________________________
618 Int_t
619 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
620 {
621   // 
622   // Find the (cached) maximum weight for FMD<i>dr</i> iat
623   // @f$\eta@f$ 
624   // 
625   // Parameters:
626   //    d     Detector
627   //    r     Ring
628   //    eta   Eta bin
629   // 
630   // Return:
631   //    max weight or <= 0 in case of problems 
632   //
633   AliForwardCorrectionManager&  fcm  = AliForwardCorrectionManager::Instance();
634   Int_t                         iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
635
636   return GetMaxWeight(d, r, iEta);
637 }
638
639 //_____________________________________________________________________
640 Float_t 
641 AliFMDDensityCalculator::NParticles(Float_t  mult, 
642                                     UShort_t d, 
643                                     Char_t   r, 
644                                     Float_t  eta,
645                                     Bool_t   lowFlux) const
646 {
647   // 
648   // Get the number of particles corresponding to the signal mult
649   // 
650   // Parameters:
651   //    mult     Signal
652   //    d        Detector
653   //    r        Ring 
654   //    s        Sector 
655   //    t        Strip (not used)
656   //    v        Vertex bin 
657   //    eta      Pseudo-rapidity 
658   //    lowFlux  Low-flux flag 
659   // 
660   // Return:
661   //    The number of particles 
662   //
663   // if (mult <= GetMultCut()) return 0;
664   DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
665   if (lowFlux) return 1;
666   
667   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
668   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta, -1);
669   if (!fit) { 
670     AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f qual=%d", 
671                     d, r, eta, fMinQuality));
672     return 0;
673   }
674   
675   Int_t    m   = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
676   if (m < 1) { 
677     AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
678     return 0;
679   }
680   
681   UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
682   Double_t ret = fit->EvaluateWeighted(mult, n);
683   
684   if (fDebug > 10) {
685     AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
686   }
687     
688   fWeightedSum->Fill(ret);
689   fSumOfWeights->Fill(ret);
690   
691   return ret;
692 }
693
694 //_____________________________________________________________________
695 Float_t 
696 AliFMDDensityCalculator::Correction(UShort_t d, 
697                                     Char_t   r, 
698                                     UShort_t t, 
699                                     Float_t  eta,
700                                     Bool_t   lowFlux) const
701 {
702   // 
703   // Get the inverse correction factor.  This consist of
704   // 
705   // - acceptance correction (corners of sensors) 
706   // - double hit correction (for low-flux events) 
707   // - dead strip correction 
708   // 
709   // Parameters:
710   //    d        Detector
711   //    r        Ring 
712   //    s        Sector 
713   //    t        Strip (not used)
714   //    v        Vertex bin 
715   //    eta      Pseudo-rapidity 
716   //    lowFlux  Low-flux flag 
717   // 
718   // Return:
719   //    
720   //
721   DGUARD(fDebug, 10, "Apply correction in FMD density calculator");
722   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
723
724   Float_t correction = 1; 
725   if (fUsePhiAcceptance == kPhiCorrectNch) 
726     correction = AcceptanceCorrection(r,t);
727   if (lowFlux) { 
728     TH1D* dblHitCor = 0;
729     if (fcm.GetDoubleHit()) 
730       dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
731
732     if (dblHitCor) {
733       Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
734       if (dblC > 0) correction *= dblC;
735     }
736     // else {
737     //   AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
738     // }
739   }
740   return correction;
741 }
742
743 //_____________________________________________________________________
744 TH1D*
745 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
746 {
747   // 
748   // Generate the acceptance corrections 
749   // 
750   // Parameters:
751   //    r Ring to generate for 
752   // 
753   // Return:
754   //    Newly allocated histogram of acceptance corrections
755   //
756   DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
757   const Double_t ic1[] = { 4.9895, 15.3560 };
758   const Double_t ic2[] = { 1.8007, 17.2000 };
759   const Double_t oc1[] = { 4.2231, 26.6638 };
760   const Double_t oc2[] = { 1.8357, 27.9500 };
761   const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
762   const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
763   Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
764   Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
765   Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
766   Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
767   Float_t   basearc    = 2 * TMath::Pi() / nSec;
768   Double_t  rad        = maxR - minR;
769   Float_t   segment    = rad / nStrips;
770   Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
771
772   // Numbers used to find end-point of strip.
773   // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
774   Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
775   Float_t dx           = c2[0] - c1[0];
776   Float_t dy           = c2[1] - c1[1];
777   Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);
778
779   TH1D* ret = new TH1D(Form("acc%c", r), 
780                        Form("Acceptance correction for FMDx%c", r), 
781                        nStrips, -.5, nStrips-.5);
782   ret->SetXTitle("Strip");
783   ret->SetYTitle("#varphi acceptance");
784   ret->SetDirectory(0);
785   ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
786   ret->SetFillStyle(3001);
787
788   for (Int_t t = 0; t < nStrips; t++) { 
789     Float_t   radius     = minR + t * segment;
790     
791     // If the radius of the strip is smaller than the radius corresponding 
792     // to the first corner we have a full strip length 
793     if (radius <= cr) {
794       ret->SetBinContent(t+1, 1);
795       continue;
796     }
797
798     // Next, we should find the end-point of the strip - that is, 
799     // the coordinates where the line from c1 to c2 intersects a circle 
800     // with radius given by the strip. 
801     // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
802     // Calculate the determinant 
803     Float_t det = radius * radius * dr * dr - D*D;
804
805     if (det <=  0) { 
806       // <0 means No intersection
807       // =0 means Exactly tangent
808       ret->SetBinContent(t+1, 1);
809       continue;
810     }
811
812     // Calculate end-point and the corresponding opening angle 
813     Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
814     Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
815     Float_t th  = TMath::ATan2(x, y);
816
817     ret->SetBinContent(t+1, th / basearc);
818   }
819   return ret;
820 }
821
822 //_____________________________________________________________________
823 Float_t 
824 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
825 {
826   // 
827   // Get the acceptance correction for strip @a t in an ring of type @a r
828   // 
829   // Parameters:
830   //    r  Ring type ('I' or 'O')
831   //    t  Strip number 
832   // 
833   // Return:
834   //    Inverse acceptance correction 
835   //
836   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
837   return acc->GetBinContent(t+1);
838 }
839
840 //____________________________________________________________________
841 void
842 AliFMDDensityCalculator::Terminate(const TList* dir, TList* output, 
843                                    Int_t nEvents)
844 {
845   // 
846   // Scale the histograms to the total number of events 
847   // 
848   // Parameters:
849   //    dir     where to put the output
850   //    nEvents Number of events 
851   //
852   DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
853   if (nEvents <= 0) return;
854   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
855   if (!d) return;
856
857   TList* out = new TList;
858   out->SetName(d->GetName());
859   out->SetOwner();
860   
861   TIter    next(&fRingHistos);
862   RingHistos* o = 0;
863   THStack* sums = new THStack("sums", "sums of ring signals");
864   while ((o = static_cast<RingHistos*>(next()))) {
865     o->Terminate(d, nEvents);
866     if (!o->fDensity) { 
867       Warning("Terminate", "No density in %s", o->GetName());
868       continue;
869     }
870     TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, 
871                                          o->fDensity->GetNbinsY(),"e");
872     sum->Scale(1., "width");
873     sum->SetTitle(o->GetName());
874     sum->SetDirectory(0);
875     sum->SetYTitle("#sum N_{ch,incl}");
876     sums->Add(sum);
877   }
878   out->Add(sums);
879   output->Add(out);
880 }
881
882 //____________________________________________________________________
883 void
884 AliFMDDensityCalculator::CreateOutputObjects(TList* dir)
885 {
886   // 
887   // Output diagnostic histograms to directory 
888   // 
889   // Parameters:
890   //    dir List to write in
891   //  
892   DGUARD(fDebug, 1, "Define output FMD density calculator");
893   TList* d = new TList;
894   d->SetOwner();
895   d->SetName(GetName());
896   dir->Add(d);
897   d->Add(fWeightedSum);
898   d->Add(fSumOfWeights);
899   d->Add(fCorrections);
900   d->Add(fAccI);
901   d->Add(fAccO);
902   d->Add(fMaxWeights);
903   d->Add(fLowCuts);
904
905   // TNamed* sigma  = new TNamed("sigma",
906   // (fIncludeSigma ? "included" : "excluded"));
907   TObject* maxP   = AliForwardUtil::MakeParameter("maxParticle", fMaxParticles);
908   TObject* method = AliForwardUtil::MakeParameter("method", fUsePoisson);
909   TObject* phiA   = AliForwardUtil::MakeParameter("phiAcceptance", 
910                                                   fUsePhiAcceptance);
911   TObject* etaL   = AliForwardUtil::MakeParameter("etaLumping", fEtaLumping);
912   TObject* phiL   = AliForwardUtil::MakeParameter("phiLumping", fPhiLumping);
913   TObject* reEt   = AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta);
914   TObject* rePh   = AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi);
915
916   // d->Add(sigma);
917   d->Add(maxP);
918   d->Add(method);
919   d->Add(phiA);
920   d->Add(etaL);
921   d->Add(phiL);
922   d->Add(reEt);
923   d->Add(rePh);
924   // d->Add(nxi);
925   fCuts.Output(d,"lCuts");
926
927   TIter    next(&fRingHistos);
928   RingHistos* o = 0;
929   while ((o = static_cast<RingHistos*>(next()))) {
930     o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
931     o->CreateOutputObjects(d);
932   }
933 }
934 //____________________________________________________________________
935 void
936 AliFMDDensityCalculator::Print(Option_t* option) const
937 {
938   // 
939   // Print information 
940   // 
941   // Parameters:
942   //    option Not used
943   //
944   char ind[gROOT->GetDirLevel()+3];
945   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
946   ind[gROOT->GetDirLevel()] = '\0';
947   std::cout << ind << ClassName() << ": " << GetName() << '\n'
948             << std::boolalpha 
949             << ind << " Max(particles):         " << fMaxParticles << '\n'
950             << ind << " Poisson method:         " << fUsePoisson << '\n'
951             << ind << " Eta lumping:            " << fEtaLumping << '\n'
952             << ind << " Phi lumping:            " << fPhiLumping << '\n'
953             << ind << " Recalculate eta:        " << fRecalculateEta << '\n'
954             << ind << " Recalculate phi:        " << fRecalculatePhi << '\n'
955             << ind << " Use phi acceptance:     "
956             << std::noboolalpha
957             << std::flush;
958   switch (fUsePhiAcceptance) { 
959   case kPhiNoCorrect:    std::cout << "none\n"; break;
960   case kPhiCorrectNch:   std::cout << "correct Nch\n"; break;
961   case kPhiCorrectELoss: std::cout << "correct energy loss\n"; break;
962   }
963   std::cout << ind << " Lower cut:" << std::endl;
964   fCuts.Print();
965   TString opt(option);
966   opt.ToLower();
967   if (opt.Contains("nomax")) return;
968   
969   std::cout << ind << " Max weights:\n";
970
971   for (UShort_t d=1; d<=3; d++) { 
972     UShort_t nr = (d == 1 ? 1 : 2);
973     for (UShort_t q=0; q<nr; q++) { 
974       ind[gROOT->GetDirLevel()]   = ' ';
975       ind[gROOT->GetDirLevel()+1] = '\0';
976       Char_t      r = (q == 0 ? 'I' : 'O');
977       std::cout << ind << " FMD" << d << r << ":";
978       ind[gROOT->GetDirLevel()+1] = ' ';
979       ind[gROOT->GetDirLevel()+2] = '\0';
980       
981       const TArrayI& a = (d == 1 ? fFMD1iMax : 
982                           (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) : 
983                            (r == 'I' ? fFMD3iMax : fFMD3oMax)));
984       Int_t j = 0;
985       for (Int_t i = 0; i < a.fN; i++) { 
986         if (a.fArray[i] < 1) continue; 
987         if (j % 6 == 0)      std::cout << "\n " << ind;
988         j++;
989         std::cout << "  " << std::setw(3) << i << ": " << a.fArray[i];
990       }
991       std::cout << std::endl;
992     }
993   }
994 }
995
996 //====================================================================
997 AliFMDDensityCalculator::RingHistos::RingHistos()
998   : AliForwardUtil::RingHistos(),
999     fList(0),
1000     // fEvsN(0), 
1001     // fEvsM(0), 
1002     // fEtaVsN(0),
1003     // fEtaVsM(0),
1004     fCorr(0),
1005     fDensity(0),
1006     fELossVsPoisson(0),
1007     fDiffELossPoisson(0),
1008     fPoisson(),
1009     fELoss(0),
1010     fELossUsed(0),
1011     fMultCut(0), 
1012     fTotal(0),
1013     fGood(0),
1014     fPhiAcc(0), 
1015     fPhiBefore(0),
1016     fPhiAfter(0)
1017 {
1018   // 
1019   // Default CTOR
1020   //
1021 }
1022
1023 //____________________________________________________________________
1024 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
1025   : AliForwardUtil::RingHistos(d,r),
1026     fList(0),
1027     // fEvsN(0), 
1028     // fEvsM(0),
1029     // fEtaVsN(0),
1030     // fEtaVsM(0),
1031     fCorr(0),
1032     fDensity(0),
1033     fELossVsPoisson(0),
1034     fDiffELossPoisson(0),
1035     fPoisson("ignored"),
1036     fELoss(0),
1037     fELossUsed(0),
1038     fMultCut(0), 
1039     fTotal(0),
1040     fGood(0),
1041     fPhiAcc(0), 
1042     fPhiBefore(0),
1043     fPhiAfter(0)
1044 {
1045   // 
1046   // Constructor
1047   // 
1048   // Parameters:
1049   //    d detector
1050   //    r ring 
1051   //
1052 #if 0
1053   fEvsN = new TH2D("elossVsNnocorr", 
1054                    "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
1055                    250, -.5, 24.5, 251, -1.5, 24.5);
1056   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
1057   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
1058   fEvsN->Sumw2();
1059   fEvsN->SetDirectory(0);
1060
1061   fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
1062   fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
1063   fEvsM->SetDirectory(0);
1064
1065   fEtaVsN = new TProfile("etaVsNnocorr",
1066                          "Average inclusive N_{ch} vs #eta (uncorrected)",
1067                          200, -4, 6);
1068   fEtaVsN->SetXTitle("#eta");
1069   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
1070   fEtaVsN->SetDirectory(0);
1071   fEtaVsN->SetLineColor(Color());
1072   fEtaVsN->SetFillColor(Color());
1073
1074   fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
1075   fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
1076   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
1077   fEtaVsM->SetDirectory(0);
1078 #endif
1079
1080   fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
1081   fCorr->SetXTitle("#eta");
1082   fCorr->SetYTitle("#LT correction#GT");
1083   fCorr->SetDirectory(0);
1084   fCorr->SetLineColor(Color());
1085   fCorr->SetFillColor(Color());
1086
1087   fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
1088                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
1089                       0, 2*TMath::Pi());
1090   fDensity->SetDirectory(0);
1091   fDensity->Sumw2();  fDensity->SetMarkerColor(Color());
1092   fDensity->SetXTitle("#eta");
1093   fDensity->SetYTitle("#phi [radians]");
1094   fDensity->SetZTitle("Inclusive N_{ch} density");
1095
1096   fELossVsPoisson = new TH2D("elossVsPoisson", 
1097                              "N_{ch} from energy loss vs from Poission",
1098                              500, 0, 100, 500, 0, 100);
1099   fELossVsPoisson->SetDirectory(0);
1100   fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
1101   fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
1102   fELossVsPoisson->SetZTitle("Correlation");
1103
1104   fDiffELossPoisson = new TH1D("diffElossPoisson",
1105                                "(N_{ch,#DeltaE}-N_{ch,Poisson})/N_{ch,Poisson}",
1106                                100, -1, 1);
1107   fDiffELossPoisson->SetDirectory(0);
1108   fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
1109   fDiffELossPoisson->SetYTitle("Frequency");
1110   fDiffELossPoisson->SetMarkerColor(Color());
1111   fDiffELossPoisson->SetFillColor(Color());
1112   fDiffELossPoisson->Sumw2();
1113                                
1114   fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips", 
1115                     640, -1, 15);
1116   fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1117   fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1118   fELoss->SetFillColor(Color()-2);
1119   fELoss->SetFillStyle(3003);
1120   fELoss->SetLineColor(kBlack);
1121   fELoss->SetLineStyle(2);
1122   fELoss->SetLineWidth(2);
1123   fELoss->SetDirectory(0);
1124
1125   fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1126   fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1127   fELossUsed->SetFillStyle(3002);
1128   fELossUsed->SetLineStyle(1);
1129   fELossUsed->SetDirectory(0);
1130
1131   fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
1132                         (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1133   fPhiBefore->SetDirectory(0);
1134   fPhiBefore->SetXTitle("#phi");
1135   fPhiBefore->SetYTitle("Events");
1136   fPhiBefore->SetMarkerColor(Color());
1137   fPhiBefore->SetLineColor(Color());
1138   fPhiBefore->SetFillColor(Color());
1139   fPhiBefore->SetMarkerStyle(20);
1140
1141   fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
1142   fPhiAfter->SetTitle("#phi distribution (after re-calc)");
1143   fPhiAfter->SetDirectory(0);
1144   
1145 }
1146 //____________________________________________________________________
1147 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
1148   : AliForwardUtil::RingHistos(o),
1149     fList(o.fList), 
1150     // fEvsN(o.fEvsN), 
1151     // fEvsM(o.fEvsM),
1152     // fEtaVsN(o.fEtaVsN),
1153     // fEtaVsM(o.fEtaVsM),
1154     fCorr(o.fCorr),
1155     fDensity(o.fDensity),
1156     fELossVsPoisson(o.fELossVsPoisson),
1157     fDiffELossPoisson(o.fDiffELossPoisson),
1158     fPoisson(o.fPoisson),
1159     fELoss(o.fELoss),
1160     fELossUsed(o.fELossUsed),
1161     fMultCut(o.fMultCut), 
1162     fTotal(o.fTotal),
1163     fGood(o.fGood),
1164     fPhiAcc(o.fPhiAcc), 
1165     fPhiBefore(o.fPhiBefore),
1166     fPhiAfter(o.fPhiAfter)
1167 {
1168   // 
1169   // Copy constructor 
1170   // 
1171   // Parameters:
1172   //    o Object to copy from 
1173   //
1174 }
1175
1176 //____________________________________________________________________
1177 AliFMDDensityCalculator::RingHistos&
1178 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1179 {
1180   // 
1181   // Assignment operator 
1182   // 
1183   // Parameters:
1184   //    o Object to assign from 
1185   // 
1186   // Return:
1187   //    Reference to this 
1188   //
1189   if (&o == this) return *this; 
1190   AliForwardUtil::RingHistos::operator=(o);
1191   
1192   // if (fEvsN)             delete fEvsN;
1193   // if (fEvsM)             delete fEvsM;
1194   // if (fEtaVsN)           delete fEtaVsN;
1195   // if (fEtaVsM)           delete fEtaVsM;
1196   if (fCorr)             delete fCorr;
1197   if (fDensity)          delete fDensity;
1198   if (fELossVsPoisson)   delete fELossVsPoisson;
1199   if (fDiffELossPoisson) delete fDiffELossPoisson;
1200   if (fTotal)            delete fTotal;
1201   if (fGood)             delete fGood;
1202   if (fPhiAcc)           delete fPhiAcc;
1203   if (fPhiBefore)        delete fPhiBefore;
1204   if (fPhiAfter)         delete fPhiAfter;
1205
1206   // fEvsN             = static_cast<TH2D*>(o.fEvsN->Clone());
1207   // fEvsM             = static_cast<TH2D*>(o.fEvsM->Clone());
1208   // fEtaVsN           = static_cast<TProfile*>(o.fEtaVsN->Clone());
1209   // fEtaVsM           = static_cast<TProfile*>(o.fEtaVsM->Clone());
1210   fCorr             = static_cast<TProfile*>(o.fCorr->Clone());
1211   fDensity          = static_cast<TH2D*>(o.fDensity->Clone());
1212   fELossVsPoisson   = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1213   fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1214   fPoisson          = o.fPoisson;
1215   fELoss            = static_cast<TH1D*>(o.fELoss->Clone());
1216   fELossUsed        = static_cast<TH1D*>(o.fELossUsed->Clone());
1217   fTotal            = static_cast<TH1D*>(o.fTotal->Clone());
1218   fGood             = static_cast<TH1D*>(o.fGood->Clone());
1219   fPhiAcc           = static_cast<TH2D*>(o.fPhiAcc->Clone());
1220   fPhiBefore        = static_cast<TH1D*>(o.fPhiBefore->Clone());
1221   fPhiAfter         = static_cast<TH1D*>(o.fPhiAfter->Clone());
1222   return *this;
1223 }
1224 //____________________________________________________________________
1225 AliFMDDensityCalculator::RingHistos::~RingHistos()
1226 {
1227   // 
1228   // Destructor 
1229   //
1230 }
1231
1232
1233 //____________________________________________________________________
1234 void
1235 AliFMDDensityCalculator::RingHistos::SetupForData(const TAxis& eAxis)
1236 {
1237   // Initialize 
1238   // This is called on first event 
1239   fPoisson.Init(-1,-1);
1240   fTotal = new TH1D("total", "Total # of strips per #eta",
1241                     eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
1242   fTotal->SetDirectory(0);
1243   fTotal->SetXTitle("#eta");
1244   fTotal->SetYTitle("# of strips");
1245   fGood = static_cast<TH1D*>(fTotal->Clone("good"));
1246   fGood->SetTitle("# of good strips per #eta");
1247   fGood->SetDirectory(0);
1248   
1249   fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}", 
1250                      eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
1251                      10, -10, 10);
1252   fPhiAcc->SetXTitle("#eta");
1253   fPhiAcc->SetYTitle("v_{z} [cm]");
1254   fPhiAcc->SetZTitle("#phi acceptance");
1255   fPhiAcc->SetDirectory(0);
1256
1257   if (fList) fList->Add(fPhiAcc);
1258 }
1259
1260 //____________________________________________________________________
1261 void
1262 AliFMDDensityCalculator::RingHistos::CreateOutputObjects(TList* dir)
1263 {
1264   // 
1265   // Make output.  This is called as part of SlaveBegin
1266   // 
1267   // Parameters:
1268   //    dir Where to put it 
1269   //
1270   TList* d = DefineOutputList(dir);
1271   // d->Add(fEvsN);
1272   // d->Add(fEvsM);
1273   // d->Add(fEtaVsN);
1274   // d->Add(fEtaVsM);
1275   d->Add(fCorr);
1276   d->Add(fDensity);
1277   d->Add(fELossVsPoisson);
1278   d->Add(fDiffELossPoisson);
1279   fPoisson.Output(d);
1280   fPoisson.GetOccupancy()->SetFillColor(Color());
1281   fPoisson.GetMean()->SetFillColor(Color());
1282   fPoisson.GetOccupancy()->SetFillColor(Color());
1283   d->Add(fELoss);
1284   d->Add(fELossUsed);
1285   d->Add(fPhiBefore);
1286   d->Add(fPhiAfter);
1287
1288   Bool_t inner = (fRing == 'I' || fRing == 'i');
1289   Int_t nStr = inner ? 512 : 256;
1290   Int_t nSec = inner ?  20 :  40;
1291   TAxis x(nStr, -.5, nStr-.5);
1292   TAxis y(nSec, -.5, nSec-.5);
1293   x.SetTitle("strip");
1294   y.SetTitle("sector");
1295   fPoisson.Define(x, y);
1296
1297   d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
1298   fList = d;
1299 }
1300
1301 //____________________________________________________________________
1302 void
1303 AliFMDDensityCalculator::RingHistos::Terminate(TList* dir, Int_t nEvents)
1304 {
1305   // 
1306   // Scale the histograms to the total number of events 
1307   // 
1308   // Parameters:
1309   //    dir     Where the output is 
1310   //    nEvents Number of events 
1311   //
1312   TList* l = GetOutputList(dir);
1313   if (!l) return; 
1314
1315   TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"inclDensity"));
1316   if (density) density->Scale(1./nEvents);
1317   fDensity = density;
1318 }
1319
1320 //____________________________________________________________________
1321 //
1322 // EOF
1323 //
1324           
1325
1326