Minor updates
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDDensityCalculator.cxx
CommitLineData
7e4038b5 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
10ClassImp(AliFMDDensityCalculator)
11#if 0
12; // For Emacs
13#endif
14
15//____________________________________________________________________
16AliFMDDensityCalculator::AliFMDDensityCalculator()
17 : TNamed(),
18 fRingHistos(),
19 fMultCut(0.3)
20{}
21
22//____________________________________________________________________
23AliFMDDensityCalculator::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//____________________________________________________________________
37AliFMDDensityCalculator::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//____________________________________________________________________
49AliFMDDensityCalculator::~AliFMDDensityCalculator()
50{
51 fRingHistos.Delete();
52}
53
54//____________________________________________________________________
55AliFMDDensityCalculator&
56AliFMDDensityCalculator::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//____________________________________________________________________
72AliFMDDensityCalculator::RingHistos*
73AliFMDDensityCalculator::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//____________________________________________________________________
87Bool_t
88AliFMDDensityCalculator::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//_____________________________________________________________________
129Float_t
130AliFMDDensityCalculator::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//_____________________________________________________________________
162Float_t
163AliFMDDensityCalculator::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//_____________________________________________________________________
193Float_t
194AliFMDDensityCalculator::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//____________________________________________________________________
219void
9d99b0dd 220AliFMDDensityCalculator::ScaleHistograms(TList* dir, Int_t nEvents)
7e4038b5 221{
222 if (nEvents <= 0) return;
9d99b0dd 223 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
224 if (!d) return;
7e4038b5 225
226 TIter next(&fRingHistos);
227 RingHistos* o = 0;
9d99b0dd 228 while ((o = static_cast<RingHistos*>(next())))
229 o->ScaleHistograms(d, nEvents);
7e4038b5 230}
231
232//____________________________________________________________________
233void
9d99b0dd 234AliFMDDensityCalculator::DefineOutput(TList* dir)
7e4038b5 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//====================================================================
247AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd 248 : AliForwardUtil::RingHistos(),
7e4038b5 249 fEvsN(0),
250 fEvsM(0),
251 fDensity(0)
252{}
253
254//____________________________________________________________________
255AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 256 : AliForwardUtil::RingHistos(d,r),
7e4038b5 257 fEvsN(0),
258 fEvsM(0),
259 fDensity(0)
260{
9d99b0dd 261 fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()),
7e4038b5 262 Form("Energy loss vs uncorrected inclusive "
9d99b0dd 263 "N_{ch} in %s", fName.Data()),
7e4038b5 264 100, -.5, 24.5, 100, -.5, 24.5);
9d99b0dd 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);
7e4038b5 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
9d99b0dd 277 fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()),
278 Form("Inclusive N_{ch} densitu in %s", fName.Data()),
7e4038b5 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//____________________________________________________________________
287AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd 288 : AliForwardUtil::RingHistos(o),
7e4038b5 289 fEvsN(o.fEvsN),
290 fEvsM(o.fEvsM),
291 fDensity(o.fDensity)
292{}
293
294//____________________________________________________________________
295AliFMDDensityCalculator::RingHistos&
296AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
297{
9d99b0dd 298 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 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//____________________________________________________________________
311AliFMDDensityCalculator::RingHistos::~RingHistos()
312{
313 if (fEvsN) delete fEvsN;
314 if (fEvsM) delete fEvsM;
315 if (fDensity) delete fDensity;
316}
317
318//____________________________________________________________________
319void
320AliFMDDensityCalculator::RingHistos::Output(TList* dir)
321{
9d99b0dd 322 TList* d = DefineOutputList(dir);
7e4038b5 323 d->Add(fEvsN);
324 d->Add(fEvsM);
325 d->Add(fDensity);
9d99b0dd 326}
327
328//____________________________________________________________________
329void
330AliFMDDensityCalculator::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);
7e4038b5 337}
338
339//____________________________________________________________________
340//
341// EOF
342//
343
344
345