Fixes
[u/mrichter/AliRoot.git] / PWG2 / 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 <TROOT.h>
16 #include <iostream>
17 #include <iomanip>
18
19 ClassImp(AliFMDDensityCalculator)
20 #if 0
21 ; // For Emacs
22 #endif 
23
24 //____________________________________________________________________
25 AliFMDDensityCalculator::AliFMDDensityCalculator()
26   : TNamed(), 
27     fRingHistos(),
28     fMultCut(0),
29     fSumOfWeights(0),
30     fWeightedSum(0),
31     fCorrections(0),
32     fMaxParticles(5),
33     fUsePoisson(false),
34     fAccI(0),
35     fAccO(0),
36     fFMD1iMax(0),
37     fFMD2iMax(0),
38     fFMD2oMax(0),
39     fFMD3iMax(0),
40     fFMD3oMax(0),
41     fEtaLumping(1), 
42     fPhiLumping(1),    
43     fDebug(0)
44 {
45   // 
46   // Constructor 
47   //
48    
49 }
50
51 //____________________________________________________________________
52 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
53   : TNamed("fmdDensityCalculator", title), 
54     fRingHistos(), 
55     fMultCut(0),
56     fSumOfWeights(0),
57     fWeightedSum(0),
58     fCorrections(0),
59     fMaxParticles(5),
60     fUsePoisson(false),
61     fAccI(0),
62     fAccO(0),
63     fFMD1iMax(0),
64     fFMD2iMax(0),
65     fFMD2oMax(0),
66     fFMD3iMax(0),
67     fFMD3oMax(0),
68     fEtaLumping(5), 
69     fPhiLumping(1),
70     fDebug(0)
71 {
72   // 
73   // Constructor 
74   // 
75   // Parameters:
76   //    name Name of object
77   //
78   fRingHistos.SetName(GetName());
79   fRingHistos.SetOwner();
80   fRingHistos.Add(new RingHistos(1, 'I'));
81   fRingHistos.Add(new RingHistos(2, 'I'));
82   fRingHistos.Add(new RingHistos(2, 'O'));
83   fRingHistos.Add(new RingHistos(3, 'I'));
84   fRingHistos.Add(new RingHistos(3, 'O'));
85   fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
86                            200, 0, 20);
87   fSumOfWeights->SetFillColor(kRed+1);
88   fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
89   fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
90                            200, 0, 20);
91   fWeightedSum->SetFillColor(kBlue+1);
92   fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
93   fCorrections  = new TH1D("corrections", "Distribution of corrections", 
94                            100, 0, 10);
95   fCorrections->SetFillColor(kBlue+1);
96   fCorrections->SetXTitle("correction");
97
98   fAccI = GenerateAcceptanceCorrection('I');
99   fAccO = GenerateAcceptanceCorrection('O');
100 }
101
102 //____________________________________________________________________
103 AliFMDDensityCalculator::AliFMDDensityCalculator(const 
104                                                  AliFMDDensityCalculator& o)
105   : TNamed(o), 
106     fRingHistos(), 
107     fMultCut(o.fMultCut),
108     fSumOfWeights(o.fSumOfWeights),
109     fWeightedSum(o.fWeightedSum),
110     fCorrections(o.fCorrections),
111     fMaxParticles(o.fMaxParticles),
112     fUsePoisson(o.fUsePoisson),
113     fAccI(o.fAccI),
114     fAccO(o.fAccO),
115     fFMD1iMax(o.fFMD1iMax),
116     fFMD2iMax(o.fFMD2iMax),
117     fFMD2oMax(o.fFMD2oMax),
118     fFMD3iMax(o.fFMD3iMax),
119     fFMD3oMax(o.fFMD3oMax),
120     fEtaLumping(o.fEtaLumping), 
121     fPhiLumping(o.fPhiLumping),
122     fDebug(o.fDebug)
123 {
124   // 
125   // Copy constructor 
126   // 
127   // Parameters:
128   //    o Object to copy from 
129   //
130   TIter    next(&o.fRingHistos);
131   TObject* obj = 0;
132   while ((obj = next())) fRingHistos.Add(obj);
133 }
134
135 //____________________________________________________________________
136 AliFMDDensityCalculator::~AliFMDDensityCalculator()
137 {
138   // 
139   // Destructor 
140   //
141   fRingHistos.Delete();
142 }
143
144 //____________________________________________________________________
145 AliFMDDensityCalculator&
146 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
147 {
148   // 
149   // Assignement operator
150   // 
151   // Parameters:
152   //    o Object to assign from 
153   // 
154   // Return:
155   //    Reference to this object
156   //
157   TNamed::operator=(o);
158
159   fMultCut      = o.fMultCut;
160   fDebug        = o.fDebug;
161   fMaxParticles = o.fMaxParticles;
162   fUsePoisson   = o.fUsePoisson;
163   fAccI         = o.fAccI;
164   fAccO         = o.fAccO;
165   fFMD1iMax     = o.fFMD1iMax;
166   fFMD2iMax     = o.fFMD2iMax;
167   fFMD2oMax     = o.fFMD2oMax;
168   fFMD3iMax     = o.fFMD3iMax;
169   fFMD3oMax     = o.fFMD3oMax;
170   fEtaLumping   = o.fEtaLumping;
171   fPhiLumping   = o.fPhiLumping;
172   fRingHistos.Delete();
173   TIter    next(&o.fRingHistos);
174   TObject* obj = 0;
175   while ((obj = next())) fRingHistos.Add(obj);
176   
177   return *this;
178 }
179
180 //____________________________________________________________________
181 void
182 AliFMDDensityCalculator::Init(const TAxis& axis)
183 {
184   // Intialize this sub-algorithm 
185   //
186   // Parameters:
187   //   etaAxis   Not used
188   CacheMaxWeights();
189
190   TIter    next(&fRingHistos);
191   RingHistos* o = 0;
192   while ((o = static_cast<RingHistos*>(next())))
193     o->Init(axis);
194 }
195
196 //____________________________________________________________________
197 AliFMDDensityCalculator::RingHistos*
198 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
199 {
200   // 
201   // Get the ring histogram container 
202   // 
203   // Parameters:
204   //    d Detector
205   //    r Ring 
206   // 
207   // Return:
208   //    Ring histogram container 
209   //
210   Int_t idx = -1;
211   switch (d) { 
212   case 1: idx = 0; break;
213   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
214   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
215   }
216   if (idx < 0 || idx >= fRingHistos.GetEntries()) {
217     AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
218     return 0;
219   }
220   
221   return static_cast<RingHistos*>(fRingHistos.At(idx));
222 }
223     
224 //____________________________________________________________________
225 Double_t
226 AliFMDDensityCalculator::GetMultCut() const
227 {
228   // 
229   // Get the multiplicity cut.  If the user has set fMultCut (via
230   // SetMultCut) then that value is used.  If not, then the lower
231   // value of the fit range for the enery loss fits is returned.
232   // 
233   // Return:
234   //    Lower cut on multiplicity
235   //
236   if (fMultCut > 0) return fMultCut;
237
238   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
239   AliFMDCorrELossFit* fits = fcm.GetELossFit();
240   return fits->GetLowCut();
241 }
242   
243 //____________________________________________________________________
244 Bool_t
245 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
246                                    AliForwardUtil::Histos& hists,
247                                    UShort_t                vtxbin, 
248                                    Bool_t                  lowFlux)
249 {
250   // 
251   // Do the calculations 
252   // 
253   // Parameters:
254   //    fmd      AliESDFMD object (possibly) corrected for sharing
255   //    hists    Histogram cache
256   //    vtxBin   Vertex bin 
257   //    lowFlux  Low flux flag. 
258   // 
259   // Return:
260   //    true on successs 
261   //
262   for (UShort_t d=1; d<=3; d++) { 
263     UShort_t nr = (d == 1 ? 1 : 2);
264     for (UShort_t q=0; q<nr; q++) { 
265       Char_t      r = (q == 0 ? 'I' : 'O');
266       UShort_t    ns= (q == 0 ?  20 :  40);
267       UShort_t    nt= (q == 0 ? 512 : 256);
268       TH2D*       h = hists.Get(d,r);
269       RingHistos* rh= GetRingHistos(d,r);
270       if (!rh) { 
271         AliError(Form("No ring histogram found for FMD%d%c", d, r));
272         fRingHistos.ls();
273         return false;
274       }
275       rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
276       
277       for (UShort_t s=0; s<ns; s++) { 
278         for (UShort_t t=0; t<nt; t++) {
279           Float_t mult = fmd.Multiplicity(d,r,s,t);
280           Float_t phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
281           Float_t eta  = fmd.Eta(d,r,s,t);
282           rh->fTotalStrips->Fill(eta, phi);
283           if (mult < GetMultCut() || mult > 20) rh->fEmptyStrips->Fill(eta,phi);
284           if (mult == AliESDFMD::kInvalidMult || mult > 20) continue;
285           
286           //if (mult == 0) { 
287           //  rh->fEmptyStrips->Fill(eta,phi);
288           //  continue;
289           // }
290           
291           Float_t n   = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
292           
293           rh->fEvsN->Fill(mult,n);
294           rh->fEtaVsN->Fill(eta, n);
295           
296           Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
297           fCorrections->Fill(c);
298           if (c > 0) n /= c;
299           rh->fEvsM->Fill(mult,n);
300           rh->fEtaVsM->Fill(eta, n);
301           rh->fCorr->Fill(eta, c);
302           
303           if (n > 0.9 && c > 0) rh->fBasicHits->Fill(eta,phi, 1./c);
304           else                  rh->fEmptyStrips->Fill(eta,phi);
305             
306           h->Fill(eta,phi,n);
307           rh->fDensity->Fill(eta,phi,n);
308         } // for t
309       } // for s 
310
311
312       // --- Loop over poisson histograms ----------------------------
313       for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) { 
314         for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
315           Double_t eLossV   = h->GetBinContent(ieta, iphi);
316           Float_t  eta      = h->GetXaxis()->GetBinCenter(ieta);
317           Float_t  phi      = h->GetYaxis()->GetBinCenter(iphi);
318           Int_t    jeta     = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
319           Int_t    jphi     = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
320           Double_t empty    = rh->fEmptyStrips->GetBinContent(jeta, jphi);
321           Double_t total    = rh->fTotalStrips->GetBinContent(jeta, jphi);
322           Double_t hits     = rh->fBasicHits->GetBinContent(ieta,iphi);
323           // Mean in region of interest 
324           Double_t poissonM = (total <= 0 || empty <= 0 ? 0 : 
325                                -TMath::Log(empty / total));
326           Double_t poissonV = 0;
327           if(poissonM > 0)
328             // Correct for counting statistics and weight by counts 
329             poissonV = (hits * poissonM) / (1 - TMath::Exp(-1*poissonM));
330           Double_t poissonE = 0 ;
331           if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
332           
333           rh->fELossVsPoisson->Fill(eLossV, poissonV);
334           rh->fEmptyVsTotal->Fill(total, empty);
335           if (fUsePoisson) {
336             h->SetBinContent(ieta,iphi,poissonV);
337             h->SetBinError(ieta,iphi,poissonE);
338           }
339         }
340       }
341         
342     } // for q
343   } // for d
344   
345   return kTRUE;
346 }
347
348 //_____________________________________________________________________
349 Int_t
350 AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
351                                        UShort_t d, Char_t r, Int_t iEta) const
352 {
353   // 
354   // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
355   // 
356   // Parameters:
357   //    cor   Correction
358   //    d     Detector 
359   //    r     Ring 
360   //    iEta  Eta bin 
361   //
362   AliFMDCorrELossFit::ELossFit* fit = cor->GetFit(d,r,iEta);
363   if (!fit) { 
364     // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
365     return -1;
366   }
367   return fit->FindMaxWeight();
368 }
369   
370 //_____________________________________________________________________
371 void
372 AliFMDDensityCalculator::CacheMaxWeights()
373 {
374   // 
375   // Find the max weights and cache them 
376   // 
377   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
378   AliFMDCorrELossFit*           cor = fcm.GetELossFit();
379   const TAxis&                  eta = cor->GetEtaAxis();
380
381   Int_t nEta = eta.GetNbins();
382   fFMD1iMax.Set(nEta);
383   fFMD2iMax.Set(nEta);
384   fFMD2oMax.Set(nEta);
385   fFMD3iMax.Set(nEta);
386   fFMD3oMax.Set(nEta);
387   
388   for (Int_t i = 0; i < nEta; i++) {
389     fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
390     fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
391     fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
392     fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
393     fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
394   }
395 }
396
397 //_____________________________________________________________________
398 Int_t
399 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
400 {
401   // 
402   // Find the (cached) maximum weight for FMD<i>dr</i> in 
403   // @f$\eta@f$ bin @a iEta
404   // 
405   // Parameters:
406   //    d     Detector
407   //    r     Ring
408   //    iEta  Eta bin
409   // 
410   // Return:
411   //    max weight or <= 0 in case of problems 
412   //
413   if (iEta < 0) return -1;
414
415   const TArrayI* max  = 0;
416   switch (d) { 
417   case 1:  max = &fFMD1iMax;                                       break;
418   case 2:  max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
419   case 3:  max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
420   }
421   if (!max) { 
422     AliWarning(Form("No array for FMD%d%c", d, r));
423     return -1;
424   }
425   
426   if (iEta >= max->fN) { 
427     AliWarning(Form("Eta bin %3d out of bounds [0,%d]", 
428                     iEta, max->fN-1));
429     return -1;
430   }
431
432   AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta, 
433                    max->At(iEta)));
434   return max->At(iEta);
435 }
436
437 //_____________________________________________________________________
438 Int_t
439 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
440 {
441   // 
442   // Find the (cached) maximum weight for FMD<i>dr</i> iat
443   // @f$\eta@f$ 
444   // 
445   // Parameters:
446   //    d     Detector
447   //    r     Ring
448   //    eta   Eta bin
449   // 
450   // Return:
451   //    max weight or <= 0 in case of problems 
452   //
453   AliForwardCorrectionManager&  fcm  = AliForwardCorrectionManager::Instance();
454   Int_t                         iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
455
456   return GetMaxWeight(d, r, iEta);
457 }
458
459 //_____________________________________________________________________
460 Float_t 
461 AliFMDDensityCalculator::NParticles(Float_t  mult, 
462                                     UShort_t d, 
463                                     Char_t   r, 
464                                     UShort_t /*s*/, 
465                                     UShort_t /*t*/, 
466                                     UShort_t /*v*/, 
467                                     Float_t  eta,
468                                     Bool_t   lowFlux) const
469 {
470   // 
471   // Get the number of particles corresponding to the signal mult
472   // 
473   // Parameters:
474   //    mult     Signal
475   //    d        Detector
476   //    r        Ring 
477   //    s        Sector 
478   //    t        Strip (not used)
479   //    v        Vertex bin 
480   //    eta      Pseudo-rapidity 
481   //    lowFlux  Low-flux flag 
482   // 
483   // Return:
484   //    The number of particles 
485   //
486   if (mult <= GetMultCut()) return 0;
487   if (lowFlux) return 1;
488   
489   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
490   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
491   if (!fit) { 
492     AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
493     return 0;
494   }
495   
496   Int_t    m   = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
497   if (m < 1) { 
498     AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
499     return 0;
500   }
501   UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
502   Double_t ret = fit->EvaluateWeighted(mult, n);
503
504   if (fDebug > 10) {
505     AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
506   }
507   
508   fWeightedSum->Fill(ret);
509   fSumOfWeights->Fill(ret);
510
511   return ret;
512 }
513
514 //_____________________________________________________________________
515 Float_t 
516 AliFMDDensityCalculator::Correction(UShort_t d, 
517                                     Char_t   r, 
518                                     UShort_t /*s*/, 
519                                     UShort_t t, 
520                                     UShort_t /*v*/, 
521                                     Float_t  eta,
522                                     Bool_t   lowFlux) const
523 {
524   // 
525   // Get the inverse correction factor.  This consist of
526   // 
527   // - acceptance correction (corners of sensors) 
528   // - double hit correction (for low-flux events) 
529   // - dead strip correction 
530   // 
531   // Parameters:
532   //    d        Detector
533   //    r        Ring 
534   //    s        Sector 
535   //    t        Strip (not used)
536   //    v        Vertex bin 
537   //    eta      Pseudo-rapidity 
538   //    lowFlux  Low-flux flag 
539   // 
540   // Return:
541   //    
542   //
543   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
544
545   Float_t correction = AcceptanceCorrection(r,t);
546   if (lowFlux) { 
547     TH1D* dblHitCor = 0;
548     if (fcm.GetDoubleHit()) 
549       dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
550
551     if (dblHitCor) {
552       Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
553       if (dblC > 0) correction *= dblC;
554     }
555     else {
556       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
557     }
558   }
559   return correction;
560 }
561
562 //_____________________________________________________________________
563 TH1D*
564 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
565 {
566   // 
567   // Generate the acceptance corrections 
568   // 
569   // Parameters:
570   //    r Ring to generate for 
571   // 
572   // Return:
573   //    Newly allocated histogram of acceptance corrections
574   //
575   const Double_t ic1[] = { 4.9895, 15.3560 };
576   const Double_t ic2[] = { 1.8007, 17.2000 };
577   const Double_t oc1[] = { 4.2231, 26.6638 };
578   const Double_t oc2[] = { 1.8357, 27.9500 };
579   const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
580   const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
581   Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
582   Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
583   Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
584   Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
585   Float_t   basearc    = 2 * TMath::Pi() / nSec;
586   Double_t  rad        = maxR - minR;
587   Float_t   segment    = rad / nStrips;
588   Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
589
590   // Numbers used to find end-point of strip.
591   // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
592   Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
593   Float_t dx           = c2[0] - c1[0];
594   Float_t dy           = c2[1] - c1[1];
595   Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);
596
597   TH1D* ret = new TH1D(Form("acc%c", r), 
598                        Form("Acceptance correction for FMDx%c", r), 
599                        nStrips, -.5, nStrips-.5);
600   ret->SetXTitle("Strip");
601   ret->SetYTitle("#varphi acceptance");
602   ret->SetDirectory(0);
603   ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
604   ret->SetFillStyle(3001);
605
606   for (Int_t t = 0; t < nStrips; t++) { 
607     Float_t   radius     = minR + t * segment;
608     
609     // If the radius of the strip is smaller than the radius corresponding 
610     // to the first corner we have a full strip length 
611     if (radius <= cr) {
612       ret->SetBinContent(t+1, 1);
613       continue;
614     }
615
616     // Next, we should find the end-point of the strip - that is, 
617     // the coordinates where the line from c1 to c2 intersects a circle 
618     // with radius given by the strip. 
619     // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
620     // Calculate the determinant 
621     Float_t det = radius * radius * dr * dr - D*D;
622
623     if (det <=  0) { 
624       // <0 means No intersection
625       // =0 means Exactly tangent
626       ret->SetBinContent(t+1, 1);
627       continue;
628     }
629
630     // Calculate end-point and the corresponding opening angle 
631     Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
632     Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
633     Float_t th  = TMath::ATan2(x, y);
634
635     ret->SetBinContent(t+1, th / basearc);
636   }
637   return ret;
638 }
639
640 //_____________________________________________________________________
641 Float_t 
642 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
643 {
644   // 
645   // Get the acceptance correction for strip @a t in an ring of type @a r
646   // 
647   // Parameters:
648   //    r  Ring type ('I' or 'O')
649   //    t  Strip number 
650   // 
651   // Return:
652   //    Inverse acceptance correction 
653   //
654   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
655   return acc->GetBinContent(t+1);
656 }
657
658 //____________________________________________________________________
659 void
660 AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents)
661 {
662   // 
663   // Scale the histograms to the total number of events 
664   // 
665   // Parameters:
666   //    dir     where to put the output
667   //    nEvents Number of events 
668   //
669   if (nEvents <= 0) return;
670   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
671   if (!d) return;
672
673   TIter    next(&fRingHistos);
674   RingHistos* o = 0;
675   while ((o = static_cast<RingHistos*>(next())))
676     o->ScaleHistograms(d, nEvents);
677 }
678
679 //____________________________________________________________________
680 void
681 AliFMDDensityCalculator::DefineOutput(TList* dir)
682 {
683   // 
684   // Output diagnostic histograms to directory 
685   // 
686   // Parameters:
687   //    dir List to write in
688   //  
689   TList* d = new TList;
690   d->SetName(GetName());
691   dir->Add(d);
692   d->Add(fWeightedSum);
693   d->Add(fSumOfWeights);
694   d->Add(fCorrections);
695   d->Add(fAccI);
696   d->Add(fAccO);
697
698   TIter    next(&fRingHistos);
699   RingHistos* o = 0;
700   while ((o = static_cast<RingHistos*>(next()))) {
701     o->Output(d);
702   }
703 }
704 //____________________________________________________________________
705 void
706 AliFMDDensityCalculator::Print(Option_t* option) const
707 {
708   // 
709   // Print information 
710   // 
711   // Parameters:
712   //    option Not used
713   //
714   char ind[gROOT->GetDirLevel()+3];
715   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
716   ind[gROOT->GetDirLevel()] = '\0';
717   std::cout << ind << "AliFMDDensityCalculator: " << GetName() << '\n'
718             << ind << " Multiplicity cut:       " << fMultCut << '\n'
719             << ind << " Max(particles):         " << fMaxParticles 
720             << std::endl;
721   TString opt(option);
722   opt.ToLower();
723   if (opt.Contains("nomax")) return;
724   
725   std::cout << ind << " Max weights:\n";
726
727   for (UShort_t d=1; d<=3; d++) { 
728     UShort_t nr = (d == 1 ? 1 : 2);
729     for (UShort_t q=0; q<nr; q++) { 
730       ind[gROOT->GetDirLevel()]   = ' ';
731       ind[gROOT->GetDirLevel()+1] = '\0';
732       Char_t      r = (q == 0 ? 'I' : 'O');
733       std::cout << ind << " FMD" << d << r << ":";
734       ind[gROOT->GetDirLevel()+1] = ' ';
735       ind[gROOT->GetDirLevel()+2] = '\0';
736       
737       const TArrayI& a = (d == 1 ? fFMD1iMax : 
738                           (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) : 
739                            (r == 'I' ? fFMD3iMax : fFMD3oMax)));
740       Int_t j = 0;
741       for (Int_t i = 0; i < a.fN; i++) { 
742         if (a.fArray[i] < 1) continue; 
743         if (j % 6 == 0)      std::cout << "\n " << ind;
744         j++;
745         std::cout << "  " << std::setw(3) << i << ": " << a.fArray[i];
746       }
747       std::cout << std::endl;
748     }
749   }
750 }
751
752 //====================================================================
753 AliFMDDensityCalculator::RingHistos::RingHistos()
754   : AliForwardUtil::RingHistos(),
755     fEvsN(0), 
756     fEvsM(0), 
757     fEtaVsN(0),
758     fEtaVsM(0),
759     fCorr(0),
760     fDensity(0),
761     fELossVsPoisson(0),
762     fTotalStrips(0),
763     fEmptyStrips(0),
764     fBasicHits(0),
765     fEmptyVsTotal(0)
766 {
767   // 
768   // Default CTOR
769   //
770 }
771
772 //____________________________________________________________________
773 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
774   : AliForwardUtil::RingHistos(d,r),
775     fEvsN(0), 
776     fEvsM(0),
777     fEtaVsN(0),
778     fEtaVsM(0),
779     fCorr(0),
780     fDensity(0),
781     fELossVsPoisson(0),
782     fTotalStrips(0),
783     fEmptyStrips(0),
784     fBasicHits(0),
785     fEmptyVsTotal(0)
786 {
787   // 
788   // Constructor
789   // 
790   // Parameters:
791   //    d detector
792   //    r ring 
793   //
794   fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), 
795                    Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
796                         "N_{ch} in %s", fName.Data()), 
797                    2500, -.5, 24.5, 2500, -.5, 24.5);
798   fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), 
799                    Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
800                         "N_{ch} in %s", fName.Data()), 
801                    2500, -.5, 24.5, 2500, -.5, 24.5);
802   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
803   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
804   fEvsN->Sumw2();
805   fEvsN->SetDirectory(0);
806   fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
807   fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
808   fEvsM->Sumw2();
809   fEvsM->SetDirectory(0);
810
811   fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
812                          Form("Average inclusive N_{ch} vs #eta (uncorrected) "
813                               "in %s", fName.Data()), 200, -4, 6);
814   fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
815                          Form("Average inclusive N_{ch} vs #eta (corrected) "
816                               "in %s", fName.Data()), 200, -4, 6);
817   fEtaVsN->SetXTitle("#eta");
818   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
819   fEtaVsN->SetDirectory(0);
820   fEtaVsN->SetLineColor(Color());
821   fEtaVsN->SetFillColor(Color());
822   fEtaVsM->SetXTitle("#eta");
823   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
824   fEtaVsM->SetDirectory(0);
825   fEtaVsM->SetLineColor(Color());
826   fEtaVsM->SetFillColor(Color());
827
828
829   fCorr = new TProfile(Form("%s_corr", fName.Data()),
830                          Form("Average correction in %s", fName.Data()), 
831                        200, -4, 6);
832   fCorr->SetXTitle("#eta");
833   fCorr->SetYTitle("#LT correction#GT");
834   fCorr->SetDirectory(0);
835   fCorr->SetLineColor(Color());
836   fCorr->SetFillColor(Color());
837
838   fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), 
839                       Form("Inclusive N_{ch} density in %s", fName.Data()), 
840                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
841                       0, 2*TMath::Pi());
842   fDensity->SetDirectory(0);
843   fDensity->SetXTitle("#eta");
844   fDensity->SetYTitle("#phi [radians]");
845   fDensity->SetZTitle("Inclusive N_{ch} density");
846
847   fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()),
848                              Form("N_{ch} from energy loss vs from Poission %s",
849                                   fName.Data()), 100, 0, 20, 100, 0, 20);
850   fELossVsPoisson->SetDirectory(0);
851   fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
852   fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
853   fELossVsPoisson->SetZTitle("Correlation");
854
855   fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data()), 
856                            Form("# of empty strips vs. total @ # strips in %s", 
857                                 fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5);
858   fEmptyVsTotal->SetDirectory(0);
859   fEmptyVsTotal->SetXTitle("Total # strips");
860   fEmptyVsTotal->SetYTitle("Empty # strips");
861   fEmptyVsTotal->SetZTitle("Correlation");  
862 }
863 //____________________________________________________________________
864 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
865   : AliForwardUtil::RingHistos(o), 
866     fEvsN(o.fEvsN), 
867     fEvsM(o.fEvsM),
868     fEtaVsN(o.fEtaVsN),
869     fEtaVsM(o.fEtaVsM),
870     fCorr(o.fCorr),
871     fDensity(o.fDensity),
872     fELossVsPoisson(o.fELossVsPoisson),
873     fTotalStrips(o.fTotalStrips),
874     fEmptyStrips(o.fEmptyStrips),
875     fBasicHits(o.fBasicHits),
876     fEmptyVsTotal(o.fEmptyVsTotal)
877 {
878   // 
879   // Copy constructor 
880   // 
881   // Parameters:
882   //    o Object to copy from 
883   //
884 }
885
886 //____________________________________________________________________
887 AliFMDDensityCalculator::RingHistos&
888 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
889 {
890   // 
891   // Assignment operator 
892   // 
893   // Parameters:
894   //    o Object to assign from 
895   // 
896   // Return:
897   //    Reference to this 
898   //
899   AliForwardUtil::RingHistos::operator=(o);
900   
901   if (fEvsN)           delete fEvsN;
902   if (fEvsM)           delete fEvsM;
903   if (fEtaVsN)         delete fEtaVsN;
904   if (fEtaVsM)         delete fEtaVsM;
905   if (fCorr)           delete fCorr;
906   if (fDensity)        delete fDensity;
907   if (fELossVsPoisson) delete fELossVsPoisson;
908   if (fTotalStrips)    delete fTotalStrips;
909   if (fEmptyStrips)    delete fEmptyStrips;
910   
911   fEvsN           = static_cast<TH2D*>(o.fEvsN->Clone());
912   fEvsM           = static_cast<TH2D*>(o.fEvsM->Clone());
913   fEtaVsN         = static_cast<TProfile*>(o.fEtaVsN->Clone());
914   fEtaVsM         = static_cast<TProfile*>(o.fEtaVsM->Clone());
915   fCorr           = static_cast<TProfile*>(o.fCorr->Clone());
916   fDensity        = static_cast<TH2D*>(o.fDensity->Clone());
917   fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson);
918   fTotalStrips    = static_cast<TH2D*>(o.fTotalStrips);
919   fEmptyStrips    = static_cast<TH2D*>(o.fEmptyStrips);
920   
921   return *this;
922 }
923 //____________________________________________________________________
924 AliFMDDensityCalculator::RingHistos::~RingHistos()
925 {
926   // 
927   // Destructor 
928   //
929   if (fEvsN)           delete fEvsN;
930   if (fEvsM)           delete fEvsM;
931   if (fEtaVsN)         delete fEtaVsN;
932   if (fEtaVsM)         delete fEtaVsM;
933   if (fCorr)           delete fCorr;
934   if (fDensity)        delete fDensity;
935   if (fELossVsPoisson) delete fELossVsPoisson;
936   if (fTotalStrips)    delete fTotalStrips;
937   if (fEmptyStrips)    delete fEmptyStrips;
938 }
939
940 //____________________________________________________________________
941 void
942 AliFMDDensityCalculator::RingHistos::ResetPoissonHistos(const TH2D* h,
943                                                         Int_t etaLumping,
944                                                         Int_t phiLumping)
945 {
946   if (fTotalStrips) { 
947     fTotalStrips->Reset();
948     fEmptyStrips->Reset();
949     fBasicHits->Reset();
950     return;
951   }
952   //Inserted by HHD
953   
954   Int_t    nEta    = h->GetNbinsX() / etaLumping;
955   Int_t    nEtaF   = h->GetNbinsX();
956   Double_t etaMin  = h->GetXaxis()->GetXmin();
957   Double_t etaMax  = h->GetXaxis()->GetXmax();
958   Int_t    nPhi    = h->GetNbinsY() / phiLumping;
959   Int_t    nPhiF   = h->GetNbinsY();
960   Double_t phiMin  = h->GetYaxis()->GetXmin();
961   Double_t phiMax  = h->GetYaxis()->GetXmax();
962
963   fTotalStrips = new TH2D(Form("total%s", fName.Data()), 
964                           Form("Total number of strips in %s", fName.Data()), 
965                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
966   fEmptyStrips = new TH2D(Form("empty%s", fName.Data()), 
967                           Form("Empty number of strips in %s", fName.Data()), 
968                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
969   fBasicHits   = new TH2D(Form("basichits%s", fName.Data()), 
970                           Form("Basic number of hits in %s", fName.Data()), 
971                           nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
972   
973   fTotalStrips->SetDirectory(0);
974   fEmptyStrips->SetDirectory(0);
975   fBasicHits->SetDirectory(0);
976   fTotalStrips->SetXTitle("#eta");
977   fEmptyStrips->SetXTitle("#eta");
978   fBasicHits->SetXTitle("#eta");
979   fTotalStrips->SetYTitle("#varphi [radians]");
980   fEmptyStrips->SetYTitle("#varphi [radians]");
981   fBasicHits->SetYTitle("#varphi [radians]");
982   fTotalStrips->Sumw2();
983   fEmptyStrips->Sumw2();
984   fBasicHits->Sumw2();
985 }
986
987 //____________________________________________________________________
988 void
989 AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
990 {
991 }
992
993 //____________________________________________________________________
994 void
995 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
996 {
997   // 
998   // Make output 
999   // 
1000   // Parameters:
1001   //    dir Where to put it 
1002   //
1003   TList* d = DefineOutputList(dir);
1004   d->Add(fEvsN);
1005   d->Add(fEvsM);
1006   d->Add(fEtaVsN);
1007   d->Add(fEtaVsM);
1008   d->Add(fCorr);
1009   d->Add(fDensity);
1010   d->Add(fELossVsPoisson);
1011   d->Add(fEmptyVsTotal);
1012   d->Add(fTotalStrips);
1013   d->Add(fEmptyStrips);
1014   d->Add(fBasicHits);
1015   
1016   
1017 }
1018
1019 //____________________________________________________________________
1020 void
1021 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1022 {
1023   // 
1024   // Scale the histograms to the total number of events 
1025   // 
1026   // Parameters:
1027   //    dir     Where the output is 
1028   //    nEvents Number of events 
1029   //
1030   TList* l = GetOutputList(dir);
1031   if (!l) return; 
1032
1033   TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
1034   if (density) density->Scale(1./nEvents);
1035 }
1036
1037 //____________________________________________________________________
1038 //
1039 // EOF
1040 //
1041           
1042
1043