2 //Analysis of the output of AliTPCtaskPID.
3 //3 6D histograms - THnSparse created in the task:
5 //TPC normalized dEdx (dEdx_rec/dNdx_mc)
6 //TPC PID probabilities
8 //The values are binned in following variables:
9 // Some of them are correlated - but THnSpase handle it
10 // ~ 14 MBy per object needed
12 // 0 - particle specie as defined in the AliPID - negatives+5 <0,9>
13 // 1 - momenta - at the entrance of the TPC
14 // 2 - tan lambda- fP[3]
17 // 5 - measurement - dEdx, dEdx/BB resp. PID probability
24 .L $ALICE_ROOT/PWG1/Macros/tpcPIDResol.C+
27 // custom draw GetProjection(THnSparse*his, Int_t i0, Int_t i1, Float_t type0, Float_t type1, Float_t p0, Float_t p1, Float_t eta0, Float_t eta1, Int_t ncl0, Int_t ncl1)
28 GetProjection(fTPCsignalNorm,5,2, -1,11, 5,10, -1.4,1.4, 120,160,0,4) ->ProfileX()->Draw();
34 //1. SetRange - part p eta ncl par 6
35 SetRange(fTPCsignal, 0,10, 0.2,5, -0.9,0.9, 120,160, 0,6);
36 SetRange(fTPCsignalNorm, 0,10, 0.2,5, -0.9,0.9, 120,160, 0,6);
37 SetRange(fTPCr, 0,100, 0.2,5, -0.9,0.9, 120,160, 0,10);
45 DrawNormSignalTypeBG();
46 DrawNormSignalTypeMom();
47 DrawNormSignalTypeBB();
51 #include "THnSparse.h"
62 Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
63 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
67 THnSparse * fTPCsignal = 0;
68 THnSparse * fTPCsignalNorm = 0;
69 THnSparse * fTPCr = 0;
75 TFile *f = new TFile("OutputPID.root");
76 TObjArray *array= (TObjArray*)f->Get("tpcTaskPID");
78 fTPCsignal = (THnSparse*)array->At(0);
79 fTPCsignalNorm = (THnSparse*)array->At(1);
80 fTPCr = (THnSparse*)array->At(2);
82 fTPCsignal->GetAxis(5)->SetTitle("dEdx_{rec}");
83 fTPCsignal->GetAxis(6)->SetTitle("dNdx_{mc}");
84 fTPCsignal->GetAxis(5)->SetName("dEdx_{rec}");
85 fTPCsignal->GetAxis(6)->SetName("dNdx_{mc}");
87 fTPCsignalNorm->GetAxis(5)->SetTitle("dEdx_{rec}/dNdx_{mc}");
88 fTPCsignalNorm->GetAxis(6)->SetTitle("dNdx_{mc}");
89 fTPCsignalNorm->GetAxis(5)->SetName("dEdx_{rec}/dNdx_{mc}");
90 fTPCsignalNorm->GetAxis(6)->SetName("dNdx_{mc}");
91 new TCanvas("dEdx study");
94 void ResetRange(THnSparse* his){
96 // reset range to full
97 for (Int_t idim=0;idim< his->GetNdimensions(); idim++){
98 his->GetAxis(idim)->SetRange(0, his->GetAxis(idim)->GetNbins());
102 void SetRange(THnSparse*his, Float_t type0, Float_t type1, Float_t p0, Float_t p1, Float_t eta0, Float_t eta1, Int_t ncl0, Int_t ncl1, Float_t p60, Float_t p61){
107 selName=Form("pic/p%f_%f_eta_%f_%f_ncl_%d_%d_p6_%f_%f",p0,p1,eta0,eta1,ncl0,ncl1,p60,p61);
108 his->GetAxis(0)->SetRangeUser(type0,type1);
109 his->GetAxis(1)->SetRangeUser(p0,p1);
110 his->GetAxis(2)->SetRangeUser(eta0,eta1);
111 his->GetAxis(4)->SetRangeUser(ncl0,ncl1);
112 his->GetAxis(6)->SetRangeUser(p60,p61);
119 TH2F* GetProjection(THnSparse*his, Int_t i0, Int_t i1, Float_t type0, Float_t type1, Float_t p0, Float_t p1, Float_t eta0, Float_t eta1, Int_t ncl0, Int_t ncl1, Float_t p60, Float_t p61){
121 SetRange(his, type0, type1, p0, p1, eta0,eta1 , ncl0, ncl1, p60,p61);
122 TH2F * res = (TH2F*) his->Projection(i0,i1);
123 res->SetXTitle(his->GetAxis(i1)->GetTitle());
124 res->SetYTitle(his->GetAxis(i0)->GetTitle());
128 TH1F* GetProjection(THnSparse*his, Int_t i0, Float_t type0, Float_t type1, Float_t p0, Float_t p1, Float_t eta0, Float_t eta1, Int_t ncl0, Int_t ncl1, Float_t p60, Float_t p61){
130 SetRange(his, type0, type1, p0, p1, eta0,eta1 , ncl0, ncl1, p60,p61);
131 TH1F * res = (TH1F*) his->Projection(i0);
132 res->SetXTitle(his->GetAxis(i0)->GetTitle());
138 void DrawSignalTypeBG(){
144 TH2F *hisdEdxBG_part[10];
145 TH1 *hisdEdxBG_partProfile[10];
146 // TH1 *hisdEdxBG_partFit[10];
147 for (Int_t ipart=0;ipart<10;ipart++){
148 fTPCsignal->GetAxis(0)->SetRange(1+ipart, 1+ipart);
149 hisdEdxBG_part[ipart] = (TH2F*)fTPCsignal->Projection(5,3);
150 if (!hisdEdxBG_part[ipart]) continue;
151 hisdEdxBG_partProfile[ipart]=(TH1F*)(hisdEdxBG_part[ipart]->ProfileX(Form("SignalTypeBG_%d_%dProf",ipart,version)));
152 // hisdEdxBG_partProfile[ipart]->SetName(;
153 hisdEdxBG_partProfile[ipart]->SetMaximum(3);
154 hisdEdxBG_partProfile[ipart]->SetXTitle("#beta#gamma");
155 hisdEdxBG_partProfile[ipart]->SetYTitle("dEdx/MIP");
156 hisdEdxBG_partProfile[ipart]->Scale(1/fNorm);
157 hisdEdxBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
158 hisdEdxBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
159 delete hisdEdxBG_part[ipart];
161 if (gPad) gPad->SetLogx(kTRUE);
162 hisdEdxBG_partProfile[0]->Draw("");
163 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
164 for (Int_t ipart=0;ipart<10;ipart++) {
165 hisdEdxBG_partProfile[ipart]->Draw("same");
166 legend->AddEntry(hisdEdxBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
169 gPad->SaveAs(selName+"_hisdEdxBG.gif");
170 gPad->SaveAs(selName+"_hisdEdxBG.eps");
173 void DrawSignalTypeMom(){
178 TH2F *hisdEdxMom_part[10];
179 TH1 *hisdEdxMom_partProfile[10];
180 // TH1 *hisdEdxMom_partFit[10];
181 for (Int_t ipart=0;ipart<10;ipart++){
182 fTPCsignal->GetAxis(0)->SetRange(1+ipart, 1+ipart);
183 hisdEdxMom_part[ipart] = (TH2F*)fTPCsignal->Projection(5,1);
184 if (!hisdEdxMom_part[ipart]) continue;
185 hisdEdxMom_partProfile[ipart]=(TH1F*)(hisdEdxMom_part[ipart]->ProfileX(Form("SignalTypeMom_%d_%dProf",ipart,version)));
187 hisdEdxMom_partProfile[ipart]->SetMaximum(5);
188 hisdEdxMom_partProfile[ipart]->SetMinimum(0);
189 hisdEdxMom_partProfile[ipart]->SetXTitle("p (GeV/c)");
190 hisdEdxMom_partProfile[ipart]->SetYTitle("dEdx/MIP");
191 hisdEdxMom_partProfile[ipart]->Scale(1/fNorm);
192 hisdEdxMom_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
193 hisdEdxMom_partProfile[ipart]->SetMarkerStyle(22+ipart);
194 delete hisdEdxMom_part[ipart];
196 if (gPad) gPad->SetLogx(kTRUE);
197 hisdEdxMom_partProfile[0]->Draw("");
198 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
199 for (Int_t ipart=0;ipart<10;ipart++) {
200 hisdEdxMom_partProfile[ipart]->Draw("same");
201 legend->AddEntry(hisdEdxMom_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
204 gPad->SaveAs(selName+"_hisdEdxMom.gif");
205 gPad->SaveAs(selName+"_hisdEdxMom.eps");
209 void DrawNormSignalTypeBG(){
215 TH2F *hisdEdxNormBG_part[10];
216 TH1 *hisdEdxNormBG_partProfile[10];
217 // TH1 *hisdEdxNormBG_partFit[10];
218 for (Int_t ipart=0;ipart<10;ipart++){
219 fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
220 hisdEdxNormBG_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,3);
221 if (!hisdEdxNormBG_part[ipart]) continue;
222 hisdEdxNormBG_partProfile[ipart]=(TH1F*)(hisdEdxNormBG_part[ipart]->ProfileX(Form("NormSignalTypeBG_%d_%dProf",ipart,version)));
223 // hisdEdxNormBG_partProfile[ipart]->SetName(;
224 hisdEdxNormBG_partProfile[ipart]->SetMaximum(1.5);
225 hisdEdxNormBG_partProfile[ipart]->SetMinimum(0.5);
226 hisdEdxNormBG_partProfile[ipart]->SetXTitle("#beta#gamma");
227 hisdEdxNormBG_partProfile[ipart]->SetYTitle("dEdx/MIP");
228 hisdEdxNormBG_partProfile[ipart]->Scale(1/fNorm);
229 hisdEdxNormBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
230 hisdEdxNormBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
231 delete hisdEdxNormBG_part[ipart];
233 if (gPad) gPad->SetLogx(kTRUE);
234 hisdEdxNormBG_partProfile[0]->Draw("");
235 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
236 for (Int_t ipart=0;ipart<10;ipart++) {
237 hisdEdxNormBG_partProfile[ipart]->Draw("same");
238 legend->AddEntry(hisdEdxNormBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
241 gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
242 gPad->SaveAs(selName+"_hisdEdxNormBG.eps");
246 void DrawNormSignalTypeMom(){
253 TH2F *hisdEdxNormMom_part[10];
254 TH1 *hisdEdxNormMom_partProfile[10];
255 // TH1 *hisdEdxNormMom_partFit[10];
256 for (Int_t ipart=0;ipart<10;ipart++){
257 fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
258 hisdEdxNormMom_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,1);
259 if (!hisdEdxNormMom_part[ipart]) continue;
260 hisdEdxNormMom_partProfile[ipart]=(TH1F*)(hisdEdxNormMom_part[ipart]->ProfileX(Form("NormSignalTypeMom_%d_%dProf",ipart,version)));
261 hisdEdxNormMom_partProfile[ipart]->Scale(1/fNorm);
262 hisdEdxNormMom_partProfile[ipart]->SetMaximum(1.5);
263 hisdEdxNormMom_partProfile[ipart]->SetMinimum(0.5);
264 hisdEdxNormMom_partProfile[ipart]->SetXTitle("p (GeV/c)");
265 hisdEdxNormMom_partProfile[ipart]->SetYTitle("dEdx_{rec}/dNdx_{mc}");
266 hisdEdxNormMom_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
267 hisdEdxNormMom_partProfile[ipart]->SetMarkerStyle(22+ipart);
268 delete hisdEdxNormMom_part[ipart];
270 if (gPad) gPad->SetLogx(kTRUE);
271 hisdEdxNormMom_partProfile[0]->Draw("");
272 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
273 for (Int_t ipart=0;ipart<10;ipart++) {
274 hisdEdxNormMom_partProfile[ipart]->Draw("same");
275 legend->AddEntry(hisdEdxNormMom_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
278 gPad->SaveAs(selName+"_hisdEdxNormMom.gif");
279 gPad->SaveAs(selName+"_hisdEdxNormMom.eps");
284 void DrawNormSignalTypeBB(){
291 TH2F *hisdEdxNormBG_part[10];
292 TH1 *hisdEdxNormBG_partProfile[10];
293 // TH1 *hisdEdxNormBG_partFit[10];
297 for (Int_t ipart=0;ipart<10;ipart++){
298 fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
299 hisdEdxNormBG_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,6);
300 if (!hisdEdxNormBG_part[ipart]) continue;
301 hisdEdxNormBG_partProfile[ipart]=(TH1F*)(hisdEdxNormBG_part[ipart]->ProfileX(Form("NormSignalTypeBB_%d_%dProf",ipart,version)));
302 hisdEdxNormBG_partProfile[ipart]->Scale(1/fNorm);
303 hisdEdxNormBG_partProfile[ipart]->SetMaximum(1.5);
304 hisdEdxNormBG_partProfile[ipart]->SetMinimum(0.9);
305 hisdEdxNormBG_partProfile[ipart]->SetXTitle("dNdx_{mc}");
306 hisdEdxNormBG_partProfile[ipart]->SetYTitle("dEdx_{rec}/dNdx_{mc}");
307 hisdEdxNormBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
308 hisdEdxNormBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
313 hisdEdxNormBG_partProfile[ipart]->Fit("pol1","q","",1,2);
314 if (hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")){
315 ival[ipart]= hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")->GetParameter(1);
316 ierr[ipart]= hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")->GetParError(1);
319 printf("%d\t\t%s%f\t%f\n",ipart,AliPID::ParticleName(ipart%5), ival[ipart], ierr[ipart]);
320 delete hisdEdxNormBG_part[ipart];
322 if (gPad) gPad->SetLogx(kFALSE);
323 hisdEdxNormBG_partProfile[0]->Draw("");
324 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
325 for (Int_t ipart=0;ipart<10;ipart++) {
326 hisdEdxNormBG_partProfile[ipart]->Draw("same");
327 legend->AddEntry(hisdEdxNormBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
330 gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
331 gPad->SaveAs(selName+"_hisdEdxNormBG.eps");
341 TH2F *hisPIDMCRCSame_part[10];
342 TH1 *hisPIDMCRCSame_partProfile[10];
344 // TH1 *hisPIDMCRCSame_partFit[10];
345 for (Int_t ipart=0;ipart<5;ipart++){
346 fTPCr->GetAxis(0)->SetRange(1+ipart,1+ipart);
347 fTPCr->GetAxis(6)->SetRange(1+ipart,1+ipart);
348 hisPIDMCRCSame_part[ipart] = (TH2F*)fTPCr->Projection(5,1);
349 if (!hisPIDMCRCSame_part[ipart]) continue;
350 hisPIDMCRCSame_partProfile[ipart]=(hisPIDMCRCSame_part[ipart]->ProfileX(Form("PIDpt_%d_%dProf",ipart,version)));
351 hisPIDMCRCSame_partProfile[ipart]->SetDirectory(0);
352 hisPIDMCRCSame_partProfile[ipart]->SetMaximum(1.0);
353 hisPIDMCRCSame_partProfile[ipart]->SetMinimum(0.0);
354 hisPIDMCRCSame_partProfile[ipart]->SetXTitle("p (GeV/c)");
355 hisPIDMCRCSame_partProfile[ipart]->SetYTitle("TPCr");
356 hisPIDMCRCSame_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
357 hisPIDMCRCSame_partProfile[ipart]->SetMarkerStyle(22+ipart);
358 delete hisPIDMCRCSame_part[ipart];
360 if (gPad) gPad->SetLogx(kTRUE);
361 hisPIDMCRCSame_partProfile[0]->Draw("");
362 TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
363 for (Int_t ipart=0;ipart<5;ipart++) {
364 hisPIDMCRCSame_partProfile[ipart]->Draw("same");
365 legend->AddEntry(hisPIDMCRCSame_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
368 gPad->SaveAs(selName+"pidPIDMCRCSame.gif");
369 gPad->SaveAs(selName+"pidPIDMCRCSame.eps");
377 TH2F * his = (TH2F*) GetProjection(fTPCsignalNorm,5,1 , 7,7, 0.2,5, -0.6,0.6, 130,160,0.5,1.3);
378 his->FitSlicesY(0,0,-1,0,"QNR",&arr)
379 TH1F * hsigma = (TH1F*)arr.At(2);
380 TH1F * hmean = (TH1F*)arr.At(1);
381 hsigma->Divide(hmean);
383 hsigma->SetMaximum(10);
384 hsigma->SetMinimum(0.0);
385 hsigma->SetMarkerStyle(22);
386 hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
388 gPad->SaveAs("pic/resoldEdx_pt.eps");
389 gPad->SaveAs("pic/resoldEdx_pt.gif");
395 TH2F * his = (TH2F*) GetProjection(fTPCsignalNorm,5,2 , 7,7, 0.6,1, -0.9,0.9, 0,160,0,1.3);
396 his->FitSlicesY(0,0,-1,0,"QNR",&arr)
397 TH1F * hsigma = (TH1F*)arr.At(2);
398 TH1F * hmean = (TH1F*)arr.At(1);
399 hsigma->Divide(hmean);
401 hsigma->SetMaximum(10);
402 hsigma->SetMinimum(0.0);
403 hsigma->SetMarkerStyle(22);
404 hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
406 gPad->SaveAs("pic/resoldEdx_theta.eps");
407 gPad->SaveAs("pic/resoldEdx_theta.gif");