PWG1/macros/tpcPIDResol.C

   1 /*
   2   //Analysis of the output of AliTPCtaskPID.
   3   //3 6D histograms  - THnSparse created in the task:
   4   //TPC raw dEdx
   5   //TPC normalized dEdx (dEdx_rec/dNdx_mc)
   6   //TPC PID probabilities
   7   //
   8   //The values are binned in following variables:
   9   // Some of them are correlated - but THnSpase handle it
  10   //                               ~ 14 MBy per object needed
  11   //
  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]
  15   // 3 - betagamma
  16   // 4 - npoints
  17   // 5 - measurement - dEdx, dEdx/BB resp.  PID probability
  18   // 6 - BB
  19
  20
  21
  22   .x ~/NimStyle.C
  23   .x ~/UliStyle.C
  24   .L $ALICE_ROOT/PWG1/Macros/tpcPIDResol.C+
  25   Init();
  26   //
  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();
  29   //
  30
  31   //
  32   // predefined draw
  33   //                         0          1          2             4         6
  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);
  38
  39   //
  40   // Draw
  41   //
  42   DrawPIDpt();
  43   DrawSignalTypeBG();
  44   DrawSignalTypeMom();
  45   DrawNormSignalTypeBG();
  46   DrawNormSignalTypeMom();
  47   DrawNormSignalTypeBB();
  48
  49 */
  50 #include "TFile.h"
  51 #include "THnSparse.h"
  52 #include "TH1F.h"
  53 #include "TH2F.h"
  54 #include "TProfile.h"
  55 #include "TF1.h"
  56
  57 #include "TLegend.h"
  58 #include "TCanvas.h"
  59
  60 #include "AliPID.h"
  61
  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};
  64
  65
  66
  67 THnSparse * fTPCsignal     = 0;
  68 THnSparse * fTPCsignalNorm = 0;
  69 THnSparse * fTPCr = 0;
  70 Float_t     fNorm=50;
  71 TString selName;
  72 Int_t version=0;
  73
  74 void Init(){
  75   TFile *f = new TFile("OutputPID.root");
  76   TObjArray *array= (TObjArray*)f->Get("tpcTaskPID");
  77   delete f;
  78   fTPCsignal     = (THnSparse*)array->At(0);
  79   fTPCsignalNorm = (THnSparse*)array->At(1);
  80   fTPCr          = (THnSparse*)array->At(2);
  81
  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}");
  86
  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");
  92 }
  93
  94 void ResetRange(THnSparse* his){
  95   //
  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());
  99   }
 100 }
 101
 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){
 103   //
 104   //
 105   //
 106   ResetRange(his);
 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);
 113   version++;
 114 }
 115
 116
 117
 118
 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){
 120
 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());
 125   return res;
 126 }
 127
 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){
 129
 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());
 133   return res;
 134 }
 135
 136
 137
 138 void DrawSignalTypeBG(){
 139   //
 140   //
 141   //
 142   //
 143   TObjArray array;
 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];
 160   }
 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");
 167   }
 168   legend->Draw();
 169   gPad->SaveAs(selName+"_hisdEdxBG.gif");
 170   gPad->SaveAs(selName+"_hisdEdxBG.eps");
 171 }
 172
 173 void DrawSignalTypeMom(){
 174   //
 175   //
 176   //
 177   TObjArray array;
 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)));
 186
 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];
 195   }
 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");
 202   }
 203   legend->Draw();
 204   gPad->SaveAs(selName+"_hisdEdxMom.gif");
 205   gPad->SaveAs(selName+"_hisdEdxMom.eps");
 206 }
 207
 208
 209 void DrawNormSignalTypeBG(){
 210   //
 211   //
 212   //
 213   //
 214   TObjArray array;
 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];
 232   }
 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");
 239   }
 240   legend->Draw();
 241   gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
 242   gPad->SaveAs(selName+"_hisdEdxNormBG.eps");
 243 }
 244
 245
 246 void DrawNormSignalTypeMom(){
 247   //
 248   //
 249   //
 250   //
 251   //
 252   TObjArray array;
 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];
 269   }
 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");
 276   }
 277   legend->Draw();
 278   gPad->SaveAs(selName+"_hisdEdxNormMom.gif");
 279   gPad->SaveAs(selName+"_hisdEdxNormMom.eps");
 280 }
 281
 282
 283
 284 void DrawNormSignalTypeBB(){
 285   //
 286   //
 287   //
 288   //
 289   //
 290   TObjArray array;
 291   TH2F *hisdEdxNormBG_part[10];
 292   TH1  *hisdEdxNormBG_partProfile[10];
 293   //  TH1  *hisdEdxNormBG_partFit[10];
 294   Double_t itypes[10];
 295   Double_t ival[10];
 296   Double_t ierr[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);
 309     itypes[ipart]=ipart;
 310     itypes[ipart]=0;
 311     itypes[ipart]=0;
 312     if (ipart%5!=0) {
 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);
 317       }
 318     }
 319     printf("%d\t\t%s%f\t%f\n",ipart,AliPID::ParticleName(ipart%5), ival[ipart], ierr[ipart]);
 320     delete hisdEdxNormBG_part[ipart];
 321   }
 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");
 328   }
 329   legend->Draw();
 330   gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
 331   gPad->SaveAs(selName+"_hisdEdxNormBG.eps");
 332 }
 333
 334
 335 void DrawPIDpt(){
 336   //
 337   //
 338   //
 339   //
 340   TObjArray array;
 341   TH2F *hisPIDMCRCSame_part[10];
 342   TH1  *hisPIDMCRCSame_partProfile[10];
 343
 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];
 359   }
 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");
 366   }
 367   legend->Draw();
 368   gPad->SaveAs(selName+"pidPIDMCRCSame.gif");
 369   gPad->SaveAs(selName+"pidPIDMCRCSame.eps");
 370 }
 371
 372
 373 /*
 374
 375
 376 TObjArray arr;
 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);
 382 hsigma->Scale(100);
 383 hsigma->SetMaximum(10);
 384 hsigma->SetMinimum(0.0);
 385 hsigma->SetMarkerStyle(22);
 386 hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
 387 hsigma->Draw();
 388 gPad->SaveAs("pic/resoldEdx_pt.eps");
 389 gPad->SaveAs("pic/resoldEdx_pt.gif");
 390
 391
 392 */
 393 /*
 394   TObjArray arr;
 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);
 400   hsigma->Scale(100);
 401   hsigma->SetMaximum(10);
 402   hsigma->SetMinimum(0.0);
 403   hsigma->SetMarkerStyle(22);
 404   hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
 405   hsigma->Draw();
 406   gPad->SaveAs("pic/resoldEdx_theta.eps");
 407   gPad->SaveAs("pic/resoldEdx_theta.gif");
 408
 409 */