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