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