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