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