New implementation of the forward multiplicity analysis.
[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
220AliFMDDensityCalculator::ScaleHistograms(Int_t nEvents)
221{
222 if (nEvents <= 0) return;
223
224 TIter next(&fRingHistos);
225 RingHistos* o = 0;
226 while ((o = static_cast<RingHistos*>(next()))) {
227 o->fDensity->Scale(1. / nEvents);
228 }
229}
230
231//____________________________________________________________________
232void
233AliFMDDensityCalculator::Output(TList* dir)
234{
235 TList* d = new TList;
236 d->SetName(GetName());
237 dir->Add(d);
238 TIter next(&fRingHistos);
239 RingHistos* o = 0;
240 while ((o = static_cast<RingHistos*>(next()))) {
241 o->Output(d);
242 }
243}
244
245//====================================================================
246AliFMDDensityCalculator::RingHistos::RingHistos()
247 : fDet(0),
248 fRing('\0'),
249 fEvsN(0),
250 fEvsM(0),
251 fDensity(0)
252{}
253
254//____________________________________________________________________
255AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
256 : fDet(d),
257 fRing(r),
258 fEvsN(0),
259 fEvsM(0),
260 fDensity(0)
261{
262 fEvsN = new TH2D(Form("FMD%d%c_Eloss_N_nocorr", d, r),
263 Form("Energy loss vs uncorrected inclusive "
264 "N_{ch} in FMD%d%c", d, r),
265 100, -.5, 24.5, 100, -.5, 24.5);
266 fEvsM = new TH2D(Form("FMD%d%c_Eloss_N_corr", d, r),
267 Form("Energy loss vs corrected inclusive N_{ch} in FMD%d%c",
268 d, r), 100, -.5, 24.5, 100, -.5, 24.5);
269 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
270 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
271 fEvsN->Sumw2();
272 fEvsN->SetDirectory(0);
273 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
274 fEvsN->SetYTitle("Inclusive N_{ch} (corrected)");
275 fEvsM->Sumw2();
276 fEvsM->SetDirectory(0);
277
278 fDensity = new TH2D(Form("FMD%d%c_Incl_Density", d, r),
279 Form("in FMD%d%c", d, r),
280 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
281 0, 2*TMath::Pi());
282 fDensity->SetDirectory(0);
283 fDensity->SetXTitle("#eta");
284 fDensity->SetYTitle("#phi [radians]");
285 fDensity->SetZTitle("Inclusive N_{ch} density");
286}
287//____________________________________________________________________
288AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
289 : TObject(o),
290 fDet(o.fDet),
291 fRing(o.fRing),
292 fEvsN(o.fEvsN),
293 fEvsM(o.fEvsM),
294 fDensity(o.fDensity)
295{}
296
297//____________________________________________________________________
298AliFMDDensityCalculator::RingHistos&
299AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
300{
301 fDet = o.fDet;
302 fRing = o.fRing;
303
304 if (fEvsN) delete fEvsN;
305 if (fEvsM) delete fEvsM;
306 if (fDensity) delete fDensity;
307
308 fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
309 fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
310 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
311
312 return *this;
313}
314//____________________________________________________________________
315AliFMDDensityCalculator::RingHistos::~RingHistos()
316{
317 if (fEvsN) delete fEvsN;
318 if (fEvsM) delete fEvsM;
319 if (fDensity) delete fDensity;
320}
321
322//____________________________________________________________________
323void
324AliFMDDensityCalculator::RingHistos::Output(TList* dir)
325{
326 TList* d = new TList;
327 d->SetName(Form("FMD%d%c", fDet, fRing));
328 d->Add(fEvsN);
329 d->Add(fEvsM);
330 d->Add(fDensity);
331 dir->Add(d);
332}
333
334//____________________________________________________________________
335//
336// EOF
337//
338
339
340