An effective FD corretion
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / SpectraV0CutVariations.C
CommitLineData
2c8e2a0d 1#if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TMath.h>
3 #include <TROOT.h>
4 #include <Riostream.h>
5 #include <TCanvas.h>
6
7 #include <TString.h>
8
9 #include <TFile.h>
10 #include <TList.h>
11 #include <TH1F.h>
12 #include <TH1D.h>
13 #include <TF2.h>
14 #include <TFitResult.h>
15 #include <TFitResultPtr.h>
16 #include <TH2F.h>
17 #include <TH3F.h>
18#endif
19
2c8e2a0d 20TFile *fAss;
21TFile *fGen;
22TFile *fRaw;
23
24void Generated();
ebebf970 25void Corrections(Float_t cMin, Float_t cMax, TString centr);
2c8e2a0d 26void RawYields(Float_t cMin, Float_t cMax, TString centr);
27void Spectra(TString centr);
28
29void SpectraV0CutVariations(Float_t cMin, Float_t cMax, TString centr) {
30
31 fRaw=TFile::Open((centr+"/AliV0CutVariations.root").Data());
32 fAss=TFile::Open((centr+"/AliV0CutVariationsMC.root").Data());
33 fGen=TFile::Open("Generated.root");
34 if (!fGen) {
35 Generated();
36 fGen=TFile::Open("Generated.root");
37 }
38
ebebf970 39 Corrections(cMin,cMax,centr);
2c8e2a0d 40 RawYields(cMin,cMax,centr);
41 Spectra(centr);
42
43 fAss->Close();
44 fGen->Close();
45 fRaw->Close();
46}
47
48TH1 *GetEfficiency(Float_t, Float_t, const Char_t *, const Char_t *);
ebebf970 49TH1 *GetFeedDown(Float_t, Float_t, TString);
2c8e2a0d 50
ebebf970 51void Corrections(Float_t cmin, Float_t cmax, TString centr) {
2c8e2a0d 52
53 TString name;
54
55 TH1 *effK0s=
56 GetEfficiency(cmin, cmax, "fK0sAs","f3dHistPrimRawPtVsYVsMultK0Short");
57 name="eff_K0s_";
58 name+=centr;
59 effK0s->SetName(name.Data());
60
61 TH1 *effLambda=
62 GetEfficiency(cmin, cmax, "fLambdaAs","f3dHistPrimRawPtVsYVsMultLambda");
63 name="eff_Lambda_";
64 name+=centr;
65 effLambda->SetName(name.Data());
ebebf970 66 TH1 *fdLambda=GetFeedDown(cmin, cmax, centr);
67 name="fd_Lambda_";
68 name+=centr;
69 fdLambda->SetName(name.Data());
70
2c8e2a0d 71
72 TH1 *effLambdaBar=
73 GetEfficiency(cmin,cmax,"fLambdaBarAs","f3dHistPrimRawPtVsYVsMultAntiLambda");
74 name="eff_LambdaBar_";
75 name+=centr;
76 effLambdaBar->SetName(name.Data());
77
ea54d2e3 78 TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d 79 effK0s->Write();
ebebf970 80 effLambda->Write(); fdLambda->Write();
2c8e2a0d 81 effLambdaBar->Write();
82 f->Close();
83}
84
85TH1 *GetEfficiency(Float_t cMin, Float_t cMax,
86 const Char_t *chis, const Char_t *znam) {
9f1ae35f 87 Double_t xbins[]={
88 0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
89 1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
90 2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
91 4.5,5.0,5.5,6.5,8.0,10.0,12.0
92 };
93 Int_t nb=sizeof(xbins)/sizeof(Double_t);
94 Int_t nb1=nb-1;
95
2c8e2a0d 96 // Numerator
97 fAss->cd();
98 TH2F *f2d=(TH2F*)gDirectory->Get(chis); f2d->Sumw2();
99 TH1D *hAs=f2d->ProjectionX("hAs",0,-1,"e");
100
101 // Denominator
102 fGen->cd();
103 TH3F *f3d = (TH3F*)gDirectory->Get(znam);
104 f3d->Sumw2();
105
106 TH1D *fpt = f3d->ProjectionX("fpt",
107 f3d->GetYaxis()->FindBin(-0.5+1e-2),
108 f3d->GetYaxis()->FindBin(+0.5-1e-2),
d58451d7 109 f3d->GetZaxis()->FindBin(cMin+1e-2),
110 f3d->GetZaxis()->FindBin(cMax-1e-2)
2c8e2a0d 111 );
112 TH1 *hMc=fpt->Rebin(nb1,"hMc",xbins);
113
114 //Efficiency
115 TH1 *eff = (TH1*)hAs->Clone();
116 eff->Divide(hAs,hMc,1,1,"b");
117
118 return eff;
119}
120
ebebf970 121Bool_t GetBinContentError(const TH1 *h, Double_t x, Double_t &c, Double_t &e) {
122 Int_t i1=h->GetXaxis()->FindFixBin(x);
123 if (i1 <=0) return kFALSE;
124 Int_t nb=h->GetNbinsX();
125 if (i1 >nb) return kFALSE;
126
127 Double_t x1=h->GetBinCenter(i1);
128
129 if (TMath::Abs(x1-x) < 1e-13) {
130 c=h->GetBinContent(i1);
131 e=h->GetBinError(i1);
132 return kTRUE;
133 }
134
135 Int_t i2 = (x1 < x) ? i1+1 : i1-1;
136 if (i2 <=0) i2=i1+1;
137 if (i2 >nb) i2=i1-1;
138 Double_t x2=h->GetBinCenter(i2);
139
140 Double_t c1=h->GetBinContent(i1);
141 Double_t c2=h->GetBinContent(i2);
142 c = c1 + (x - x1)*(c2 - c1)/(x2 - x1);
143 e = 0.5*(h->GetBinError(i1) + h->GetBinError(i2));
144
145 return kTRUE;
146}
147
148TH1 *GetFeedDown(Float_t cMin, Float_t cMax, TString cent) {
149 // Get the FD matrix
150 fAss->cd();
151 TH3F *f3d = (TH3F*)gDirectory->Get("fLambdaFromXi");
152 f3d->Sumw2();
153 TH2D *fdMat = (TH2D*)f3d->Project3D("zxe");
154 const Int_t nbx=fdMat->GetNbinsX();
155 const Int_t nby=fdMat->GetNbinsY();
156 TAxis *xiPtAxis=fdMat->GetYaxis();
157
158 // Re-normalise the FD matrix with MC Xi spectrum
159 fGen->cd();
160 f3d = (TH3F*)gDirectory->Get("f3dHistGenPtVsYVsMultXiMinus");
161 f3d->Sumw2();
162
163 TH1D *hMcXi = f3d->ProjectionX("fpt",
164 f3d->GetYaxis()->FindBin(-0.5+1e-2),
165 f3d->GetYaxis()->FindBin(+0.5-1e-2),
d58451d7 166 f3d->GetZaxis()->FindBin(cMin+1e-2),
167 f3d->GetZaxis()->FindBin(cMax-1e-2)
ebebf970 168 );
169
170 TH2D *h2McXi=new TH2D(*fdMat);
171 h2McXi->Reset();
172 h2McXi->Sumw2();
173
174 for (Int_t i=1; i<=nbx; i++) {
175 for (Int_t j=1; j<=nby; j++) {
176 Double_t pt=xiPtAxis->GetBinCenter(j);
177 Double_t c=0.,e=0.;
178 if (!GetBinContentError(hMcXi,pt,c,e)) continue;
179 h2McXi->SetBinContent(i,j,c);
180 h2McXi->SetBinError(i,j,e);
181 }
182 }
183 fdMat->Divide(h2McXi);
184 delete h2McXi;
185
186
187 // Multiply the re-normalised matrix with the real Xi spectrum
188 TString fileName("systcorrectedpt_cent");
189 //*** "Dictionary" ***
190 if (cent=="0005") cent="0010";
191 if (cent=="6080") cent="6090";
192 if (cent=="8090") cent="6090";
193 TFile *f=TFile::Open((fileName+cent+".root").Data());
194 TH1F *hReXi=(TH1F*)gDirectory->Get("correctedpt_0");
195
196 TH2D *h2ReXi=new TH2D(*fdMat);
197 h2ReXi->Reset();
198 h2ReXi->Sumw2();
199
200 for (Int_t i=1; i<=nbx; i++) {
201 for (Int_t j=1; j<=nby; j++) {
202 Double_t pt=xiPtAxis->GetBinCenter(j);
203 Double_t c=0.,e=0.;
204 if (!GetBinContentError(hReXi,pt,c,e)) continue;
205 h2ReXi->SetBinContent(i,j,c);
206 h2ReXi->SetBinError(i,j,e);
207 }
208 }
209 fdMat->Multiply(h2ReXi);
210 delete h2ReXi;
211 f->Close();
212
213 TH1 *fd=fdMat->ProjectionX("_px",0,-1,"e");
214 fd->Scale(0.1,"width");
215 return fd;
216}
217
218
2c8e2a0d 219void RawYields(Float_t cMin, Float_t cMax, TString centr) {
220 TString name;
221 TH2F *f2d=0;
222
223 //+++ Number of events for normalisation
ea54d2e3 224 TFile *file=TFile::Open("LHC10h_pass2/Merged.root");
2c8e2a0d 225 TList *v0list=(TList *)gFile->Get("PWGLFExtractV0_PP/clistV0");
226 TH1F *fMult=(TH1F*)v0list->FindObject("fHistMultiplicity");
227 Int_t i1=fMult->GetXaxis()->FindBin(cMin+1e-2);
d58451d7 228 Int_t i2=fMult->GetXaxis()->FindBin(cMax-1e-2);
2c8e2a0d 229 Float_t nEvents=fMult->Integral(i1,i2);
230 file->Close();
231
232 fRaw->cd();
233
234 name="raw_K0s_";
235 name+=centr;
236 f2d=(TH2F*)gDirectory->Get("fK0sSi"); f2d->Sumw2();
237 TH1D *rawK0s=f2d->ProjectionX(name,0,-1,"e");
238 rawK0s->Scale(1/nEvents,"width");
239
240 name="raw_Lambda_";
241 name+=centr;
242 f2d=(TH2F*)gDirectory->Get("fLambdaSi"); f2d->Sumw2();
243 TH1D *rawLambda=f2d->ProjectionX(name,0,-1,"e");
244 rawLambda->Scale(1/nEvents,"width");
245
246 name="raw_LambdaBar_";
247 name+=centr;
248 f2d=(TH2F*)gDirectory->Get("fLambdaBarSi"); f2d->Sumw2();
249 TH1D *rawLambdaBar=f2d->ProjectionX(name,0,-1,"e");
250 rawLambdaBar->Scale(1/nEvents,"width");
251
ea54d2e3 252 TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d 253 rawK0s->Write();
254 rawLambda->Write();
255 rawLambdaBar->Write();
256 f->Close();
257}
258
8fb8d5db 259Double_t fd(Double_t x) {
260 //Effective FD correction
261 return 0.1619 + 0.05295*x - 0.01749*x*x + 0.001425*x*x*x - 3.446e-05*x*x*x*x;
262}
263
2c8e2a0d 264void Spectra(TString centr) {
265 TString name;
266 TH1 *eff=0;
267 TH1D *raw=0;
268 TH1D *spe=0;
269
ea54d2e3 270 TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d 271 name="eff_K0s_";
272 name+=centr;
273 eff = (TH1*)gDirectory->Get(name.Data());
274 name="raw_K0s_";
275 name+=centr;
276 raw = (TH1D*)gDirectory->Get(name.Data());
277 spe = new TH1D(*raw);
278 spe->Divide(eff);
279 name="K0s_";
280 name+=centr;
281 spe->SetName(name.Data());
282 spe->Write();
283
284 name="eff_Lambda_";
285 name+=centr;
286 eff = (TH1*)gDirectory->Get(name.Data());
287 name="raw_Lambda_";
288 name+=centr;
289 raw = (TH1D*)gDirectory->Get(name.Data());
290 spe = new TH1D(*raw);
8fb8d5db 291 for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
292 Double_t pt=spe->GetBinCenter(i);
293 Double_t c=spe->GetBinContent(i);
294 c -= (c*fd(pt));
295 spe->SetBinContent(i,c);
296 }
2c8e2a0d 297 spe->Divide(eff);
298 name="Lambda_";
299 name+=centr;
300 spe->SetName(name.Data());
301 spe->Write();
302
303 name="eff_LambdaBar_";
304 name+=centr;
305 eff = (TH1*)gDirectory->Get(name.Data());
306 name="raw_LambdaBar_";
307 name+=centr;
308 raw = (TH1D*)gDirectory->Get(name.Data());
309 spe = new TH1D(*raw);
8fb8d5db 310 for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
311 Double_t pt=spe->GetBinCenter(i);
312 Double_t c=spe->GetBinContent(i);
313 c -= (c*fd(pt));
314 spe->SetBinContent(i,c);
315 }
2c8e2a0d 316 spe->Divide(eff);
317 name="LambdaBar_";
318 name+=centr;
319 spe->SetName(name.Data());
320 spe->Write();
321 f->Close();
322}
323
324void Generated() {
325 TList *v0listMC=0;
326 TH3F *h3=0;
327 //
328 TFile::Open("LHC11a10b_plus/Merged.root");
329 v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
330
331 TH3F *k0s =
332 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
333 TH3F *k0s_nonInj =
334 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
335
336 TH3F *lam =
337 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
338 TH3F *lam_nonInj =
339 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
340
341 TH3F *alam =
342 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
343 TH3F *alam_nonInj =
344 (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
345
ebebf970 346 TH3F *xiMinus =
347 (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
348 TH3F *xiPlus =
349 (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
350
2c8e2a0d 351
352 //
353 TFile::Open("LHC11a10b_bis/Merged.root");
354 v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
355
356 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
357 k0s->Add(h3);
358 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
359 k0s_nonInj->Add(h3);
360
361 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
362 lam->Add(h3);
363 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
364 lam_nonInj->Add(h3);
365
366 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
367 alam->Add(h3);
368 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
369 alam_nonInj->Add(h3);
370
ebebf970 371 h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
372 xiMinus->Add(h3);
373 h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
374 xiPlus->Add(h3);
375
376
2c8e2a0d 377 //
378 TFile::Open("LHC11a10a_bis/Merged.root");
7e9df831 379 v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
380
2c8e2a0d 381 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
382 k0s_nonInj->Add(h3);
383 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
384 lam_nonInj->Add(h3);
385 h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
386 alam_nonInj->Add(h3);
ebebf970 387 h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
388 xiMinus->Add(h3);
389 h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
390 xiPlus->Add(h3);
2c8e2a0d 391
392 TFile *f=TFile::Open("Generated.root","new");
393 k0s->Write(); k0s_nonInj->Write();
394 lam->Write(); lam_nonInj->Write();
395 alam->Write(); alam_nonInj->Write();
ebebf970 396 xiMinus->Write();
397 xiPlus->Write();
2c8e2a0d 398 f->Close();
399}
ea54d2e3 400