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