Make sure that histograms are obtained from output list in
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDDensityCalculator.cxx
1 #include "AliFMDDensityCalculator.h"
2 #include <AliESDFMD.h>
3 #include <TAxis.h>
4 #include <TList.h>
5 #include <TMath.h>
6 #include "AliFMDAnaParameters.h"
7 #include "AliLog.h"
8 #include <TH2D.h>
9
10 ClassImp(AliFMDDensityCalculator)
11 #if 0
12 ; // For Emacs
13 #endif 
14
15 //____________________________________________________________________
16 AliFMDDensityCalculator::AliFMDDensityCalculator()
17   : TNamed(), 
18     fRingHistos(),
19     fMultCut(0.3)
20 {}
21
22 //____________________________________________________________________
23 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
24   : TNamed("fmdDensityCalculator", title), 
25     fRingHistos(), 
26     fMultCut(0.3)
27 {
28   fRingHistos.SetName(GetName());
29   fRingHistos.Add(new RingHistos(1, 'I'));
30   fRingHistos.Add(new RingHistos(2, 'I'));
31   fRingHistos.Add(new RingHistos(2, 'O'));
32   fRingHistos.Add(new RingHistos(3, 'I'));
33   fRingHistos.Add(new RingHistos(3, 'O'));
34 }
35
36 //____________________________________________________________________
37 AliFMDDensityCalculator::AliFMDDensityCalculator(const 
38                                                  AliFMDDensityCalculator& o)
39   : TNamed(o), 
40     fRingHistos(), 
41     fMultCut(o.fMultCut)
42 {
43   TIter    next(&o.fRingHistos);
44   TObject* obj = 0;
45   while ((obj = next())) fRingHistos.Add(obj);
46 }
47
48 //____________________________________________________________________
49 AliFMDDensityCalculator::~AliFMDDensityCalculator()
50 {
51   fRingHistos.Delete();
52 }
53
54 //____________________________________________________________________
55 AliFMDDensityCalculator&
56 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
57 {
58   SetName(o.GetName());
59   SetTitle(o.GetTitle());
60
61   fMultCut = o.fMultCut;
62
63   fRingHistos.Delete();
64   TIter    next(&o.fRingHistos);
65   TObject* obj = 0;
66   while ((obj = next())) fRingHistos.Add(obj);
67   
68   return *this;
69 }
70
71 //____________________________________________________________________
72 AliFMDDensityCalculator::RingHistos*
73 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
74 {
75   Int_t idx = -1;
76   switch (d) { 
77   case 1: idx = 0; break;
78   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
79   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
80   }
81   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
82   
83   return static_cast<RingHistos*>(fRingHistos.At(idx));
84 }
85     
86 //____________________________________________________________________
87 Bool_t
88 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
89                                    AliForwardUtil::Histos& hists,
90                                    Int_t                   vtxbin, 
91                                    Bool_t                  lowFlux)
92 {
93   for (UShort_t d=1; d<=3; d++) { 
94     UShort_t nr = (d == 1 ? 1 : 2);
95     for (UShort_t q=0; q<nr; q++) { 
96       Char_t      r = (q == 0 ? 'I' : 'O');
97       UShort_t    ns= (q == 0 ?  20 :  40);
98       UShort_t    nt= (q == 0 ? 512 : 256);
99       TH2D*       h = hists.Get(d,r);
100       RingHistos* rh= GetRingHistos(d,r);
101
102       for (UShort_t s=0; s<ns; s++) { 
103         for (UShort_t t=0; t<nt; t++) {
104           Float_t mult = fmd.Multiplicity(d,r,s,t);
105           
106           if (mult == 0 || mult > 20) continue;
107
108           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
109           Float_t eta = fmd.Eta(d,r,s,t);
110           
111           Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
112           rh->fEvsN->Fill(mult,n);
113
114           Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
115           if (c > 0) n /= c;
116           rh->fEvsM->Fill(mult,n);
117
118           h->Fill(eta,phi,n);
119           rh->fDensity->Fill(eta,phi,n);
120         } // for t
121       } // for s 
122     } // for q
123   } // for d
124   
125   return kTRUE;
126 }
127
128 //_____________________________________________________________________
129 Float_t 
130 AliFMDDensityCalculator::NParticles(Float_t  mult, 
131                                     UShort_t d, 
132                                     Char_t   r, 
133                                     UShort_t /*s*/, 
134                                     UShort_t /*t*/, 
135                                     Int_t    /*v*/, 
136                                     Float_t  eta,
137                                     Bool_t   lowFlux) const
138 {
139   if (mult <= fMultCut) return 0;
140   if (lowFlux) return 1;
141
142   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
143
144   Float_t mpv  = pars->GetMPV(d,r,eta);
145   Float_t w    = pars->GetSigma(d,r,eta);
146   Float_t w2   = pars->Get2MIPWeight(d,r,eta);
147   Float_t w3   = pars->Get3MIPWeight(d,r,eta);
148   Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
149   Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
150   
151   Float_t sum  = (TMath::Landau(mult,mpv,w,kTRUE) +
152                   w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
153                   w3  * TMath::Landau(mult,mpv3,3*w,kTRUE));
154   Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
155                   2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
156                   3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
157   
158   return (sum > 0) ? wsum / sum : 1;
159 }
160
161 //_____________________________________________________________________
162 Float_t 
163 AliFMDDensityCalculator::Correction(UShort_t d, Char_t r, UShort_t /*s*/, 
164                                     UShort_t t, Int_t v, Float_t eta,
165                                     Bool_t lowFlux) const
166 {
167   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
168
169   Float_t correction = AcceptanceCorrection(r,t);
170   if (lowFlux) { 
171     TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
172     if (dblHitCor) {
173       Float_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
174       if (dblC > 0) 
175         correction *= dblC;
176     }
177     else 
178       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
179   }
180   
181   TH1F* deadCor = pars->GetFMDDeadCorrection(v);
182   if (deadCor) { 
183     Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
184     if (deadC > 0) 
185       correction *= deadC; 
186   }
187   
188   return correction;
189 }
190
191
192 //_____________________________________________________________________
193 Float_t 
194 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
195 {
196   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
197   
198   //AliFMDRing fmdring(ring);
199   //fmdring.Init();
200   Float_t   rad       = pars->GetMaxR(r)-pars->GetMinR(r);
201   Float_t   slen      = pars->GetStripLength(r,t);
202   Float_t   sblen     = pars->GetBaseStripLength(r,t);
203   Float_t   nstrips   = (r == 'I' ? 512 : 256);
204   Float_t   segment   = rad / nstrips;
205   Float_t   radius    = pars->GetMinR(r) + segment*t;
206   
207   Float_t   basearea1 = 0.5*sblen*TMath::Power(radius,2);
208   Float_t   basearea2 = 0.5*sblen*TMath::Power((radius-segment),2);
209   Float_t   basearea  = basearea1 - basearea2;
210   
211   Float_t   area1     = 0.5*slen*TMath::Power(radius,2);
212   Float_t   area2     = 0.5*slen*TMath::Power((radius-segment),2);
213   Float_t   area      = area1 - area2;
214   
215   return area/basearea;
216 }
217
218 //____________________________________________________________________
219 void
220 AliFMDDensityCalculator::ScaleHistograms(TList* dir, Int_t nEvents)
221 {
222   if (nEvents <= 0) return;
223   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
224   if (!d) return;
225
226   TIter    next(&fRingHistos);
227   RingHistos* o = 0;
228   while ((o = static_cast<RingHistos*>(next())))
229     o->ScaleHistograms(d, nEvents);
230 }
231
232 //____________________________________________________________________
233 void
234 AliFMDDensityCalculator::DefineOutput(TList* dir)
235 {
236   TList* d = new TList;
237   d->SetName(GetName());
238   dir->Add(d);
239   TIter    next(&fRingHistos);
240   RingHistos* o = 0;
241   while ((o = static_cast<RingHistos*>(next()))) {
242     o->Output(d);
243   }
244 }
245
246 //====================================================================
247 AliFMDDensityCalculator::RingHistos::RingHistos()
248   : AliForwardUtil::RingHistos(),
249     fEvsN(0), 
250     fEvsM(0), 
251     fDensity(0)
252 {}
253
254 //____________________________________________________________________
255 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
256   : AliForwardUtil::RingHistos(d,r),
257     fEvsN(0), 
258     fEvsM(0),
259     fDensity(0)
260 {
261   fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), 
262                    Form("Energy loss vs uncorrected inclusive "
263                         "N_{ch} in %s", fName.Data()), 
264                    100, -.5, 24.5, 100, -.5, 24.5);
265   fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), 
266                    Form("Energy loss vs corrected inclusive N_{ch} in %s",
267                         fName.Data()), 100, -.5, 24.5, 100, -.5, 24.5);
268   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
269   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
270   fEvsN->Sumw2();
271   fEvsN->SetDirectory(0);
272   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
273   fEvsN->SetYTitle("Inclusive N_{ch} (corrected)");
274   fEvsM->Sumw2();
275   fEvsM->SetDirectory(0);
276
277   fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), 
278                       Form("Inclusive N_{ch} densitu in %s", fName.Data()), 
279                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
280                       0, 2*TMath::Pi());
281   fDensity->SetDirectory(0);
282   fDensity->SetXTitle("#eta");
283   fDensity->SetYTitle("#phi [radians]");
284   fDensity->SetZTitle("Inclusive N_{ch} density");
285 }
286 //____________________________________________________________________
287 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
288   : AliForwardUtil::RingHistos(o), 
289     fEvsN(o.fEvsN), 
290     fEvsM(o.fEvsM),
291     fDensity(o.fDensity)
292 {}
293
294 //____________________________________________________________________
295 AliFMDDensityCalculator::RingHistos&
296 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
297 {
298   AliForwardUtil::RingHistos::operator=(o);
299   
300   if (fEvsN)    delete  fEvsN;
301   if (fEvsM)    delete  fEvsM;
302   if (fDensity) delete fDensity;
303   
304   fEvsN    = static_cast<TH2D*>(o.fEvsN->Clone());
305   fEvsM    = static_cast<TH2D*>(o.fEvsM->Clone());
306   fDensity = static_cast<TH2D*>(o.fDensity->Clone());
307   
308   return *this;
309 }
310 //____________________________________________________________________
311 AliFMDDensityCalculator::RingHistos::~RingHistos()
312 {
313   if (fEvsN)    delete fEvsN;
314   if (fEvsM)    delete fEvsM;
315   if (fDensity) delete fDensity;
316 }
317
318 //____________________________________________________________________
319 void
320 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
321 {
322   TList* d = DefineOutputList(dir);
323   d->Add(fEvsN);
324   d->Add(fEvsM);
325   d->Add(fDensity);
326 }
327
328 //____________________________________________________________________
329 void
330 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
331 {
332   TList* l = GetOutputList(dir);
333   if (!l) return; 
334
335   TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
336   if (density) density->Scale(1./nEvents);
337 }
338
339 //____________________________________________________________________
340 //
341 // EOF
342 //
343           
344
345