[u/mrichter/AliRoot.git] / PWGPP / macros / tpcPIDResol.C

/*
  //Analysis of the output of AliTPCtaskPID.
  //3 6D histograms  - THnSparse created in the task:
  //TPC raw dEdx
  //TPC normalized dEdx (dEdx_rec/dNdx_mc)
  //TPC PID probabilities
  //
  //The values are binned in following variables:
  // Some of them are correlated - but THnSpase handle it  
  //                               ~ 14 MBy per object needed
  //
  // 0 - particle specie as defined in the AliPID - negatives+5 <0,9>
  // 1 - momenta - at the entrance of the TPC
  // 2 - tan lambda- fP[3]
  // 3 - betagamma
  // 4 - npoints
  // 5 - measurement - dEdx, dEdx/BB resp.  PID probability
  // 6 - BB


  .x ~/NimStyle.C
  .x ~/UliStyle.C
  .L $ALICE_ROOT/PWGPP/Macros/tpcPIDResol.C+
  Init();
  //
  // 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)
  GetProjection(fTPCsignalNorm,5,2,  -1,11,   5,10,  -1.4,1.4,   120,160,0,4)  ->ProfileX()->Draw();
  //

  //
  // predefined draw
  //                         0          1          2             4         6
  //1. SetRange -             part      p          eta           ncl       par 6
  SetRange(fTPCsignal,      0,10,    0.2,5,     -0.9,0.9,    120,160,  0,6);
  SetRange(fTPCsignalNorm,  0,10,    0.2,5,     -0.9,0.9,    120,160,  0,6);
  SetRange(fTPCr,           0,100,    0.2,5,     -0.9,0.9,    120,160,  0,10);
  
  //
  // Draw
  //
  DrawPIDpt();
  DrawSignalTypeBG();
  DrawSignalTypeMom();
  DrawNormSignalTypeBG();
  DrawNormSignalTypeMom();
  DrawNormSignalTypeBB();

*/
#include "TFile.h"
#include "THnSparse.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TF1.h"

#include "TLegend.h"
#include "TCanvas.h"

#include "AliPID.h"

Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};


THnSparse * fTPCsignal     = 0;
THnSparse * fTPCsignalNorm = 0;
THnSparse * fTPCr = 0;
Float_t     fNorm=50;
TString selName;
Int_t version=0;

void Init(){
  TFile *f = new TFile("OutputPID.root");
  TObjArray *array= (TObjArray*)f->Get("tpcTaskPID");
  delete f;
  fTPCsignal     = (THnSparse*)array->At(0);
  fTPCsignalNorm = (THnSparse*)array->At(1);
  fTPCr          = (THnSparse*)array->At(2);

  fTPCsignal->GetAxis(5)->SetTitle("dEdx_{rec}");
  fTPCsignal->GetAxis(6)->SetTitle("dNdx_{mc}");
  fTPCsignal->GetAxis(5)->SetName("dEdx_{rec}");
  fTPCsignal->GetAxis(6)->SetName("dNdx_{mc}");

  fTPCsignalNorm->GetAxis(5)->SetTitle("dEdx_{rec}/dNdx_{mc}");
  fTPCsignalNorm->GetAxis(6)->SetTitle("dNdx_{mc}");
  fTPCsignalNorm->GetAxis(5)->SetName("dEdx_{rec}/dNdx_{mc}");
  fTPCsignalNorm->GetAxis(6)->SetName("dNdx_{mc}");
  new TCanvas("dEdx study");
}

void ResetRange(THnSparse* his){
  //
  // reset range to full
  for (Int_t idim=0;idim< his->GetNdimensions(); idim++){
    his->GetAxis(idim)->SetRange(0, his->GetAxis(idim)->GetNbins());
  }    
}

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){
  //
  //
  //
  ResetRange(his);
  selName=Form("pic/p%f_%f_eta_%f_%f_ncl_%d_%d_p6_%f_%f",p0,p1,eta0,eta1,ncl0,ncl1,p60,p61);
  his->GetAxis(0)->SetRangeUser(type0,type1);
  his->GetAxis(1)->SetRangeUser(p0,p1);
  his->GetAxis(2)->SetRangeUser(eta0,eta1);
  his->GetAxis(4)->SetRangeUser(ncl0,ncl1);
  his->GetAxis(6)->SetRangeUser(p60,p61);
  version++;
}


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){

  SetRange(his,  type0, type1,   p0, p1,  eta0,eta1 ,   ncl0, ncl1,   p60,p61);
  TH2F * res = (TH2F*) his->Projection(i0,i1);
  res->SetXTitle(his->GetAxis(i1)->GetTitle());
  res->SetYTitle(his->GetAxis(i0)->GetTitle());
  return res;
}

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){

  SetRange(his,  type0, type1,   p0, p1,  eta0,eta1 ,   ncl0, ncl1,   p60,p61);
  TH1F * res = (TH1F*) his->Projection(i0);
  res->SetXTitle(his->GetAxis(i0)->GetTitle());
  return res;
}


void DrawSignalTypeBG(){
  //
  //
  //
  //
  TObjArray array;
  TH2F *hisdEdxBG_part[10];
  TH1  *hisdEdxBG_partProfile[10];
  //  TH1  *hisdEdxBG_partFit[10];
  for (Int_t ipart=0;ipart<10;ipart++){  
    fTPCsignal->GetAxis(0)->SetRange(1+ipart, 1+ipart);
    hisdEdxBG_part[ipart] = (TH2F*)fTPCsignal->Projection(5,3);
    if (!hisdEdxBG_part[ipart]) continue;
    hisdEdxBG_partProfile[ipart]=(TH1F*)(hisdEdxBG_part[ipart]->ProfileX(Form("SignalTypeBG_%d_%dProf",ipart,version)));
    //    hisdEdxBG_partProfile[ipart]->SetName(;
    hisdEdxBG_partProfile[ipart]->SetMaximum(3);
    hisdEdxBG_partProfile[ipart]->SetXTitle("#beta#gamma");
    hisdEdxBG_partProfile[ipart]->SetYTitle("dEdx/MIP");
    hisdEdxBG_partProfile[ipart]->Scale(1/fNorm);
    hisdEdxBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisdEdxBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
    delete hisdEdxBG_part[ipart];
  }
  if (gPad) gPad->SetLogx(kTRUE);
  hisdEdxBG_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<10;ipart++)  {
    hisdEdxBG_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisdEdxBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"_hisdEdxBG.gif");
  gPad->SaveAs(selName+"_hisdEdxBG.eps");
}

void DrawSignalTypeMom(){
  //
  //
  //
  TObjArray array;
  TH2F *hisdEdxMom_part[10];
  TH1  *hisdEdxMom_partProfile[10];
  //  TH1  *hisdEdxMom_partFit[10];
  for (Int_t ipart=0;ipart<10;ipart++){  
    fTPCsignal->GetAxis(0)->SetRange(1+ipart, 1+ipart);
    hisdEdxMom_part[ipart] = (TH2F*)fTPCsignal->Projection(5,1);
    if (!hisdEdxMom_part[ipart]) continue;
    hisdEdxMom_partProfile[ipart]=(TH1F*)(hisdEdxMom_part[ipart]->ProfileX(Form("SignalTypeMom_%d_%dProf",ipart,version)));
    
    hisdEdxMom_partProfile[ipart]->SetMaximum(5);
    hisdEdxMom_partProfile[ipart]->SetMinimum(0);
    hisdEdxMom_partProfile[ipart]->SetXTitle("p (GeV/c)");
    hisdEdxMom_partProfile[ipart]->SetYTitle("dEdx/MIP");
    hisdEdxMom_partProfile[ipart]->Scale(1/fNorm);
    hisdEdxMom_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisdEdxMom_partProfile[ipart]->SetMarkerStyle(22+ipart);
    delete hisdEdxMom_part[ipart];
  } 
  if (gPad) gPad->SetLogx(kTRUE);
  hisdEdxMom_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<10;ipart++)  {
    hisdEdxMom_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisdEdxMom_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"_hisdEdxMom.gif");
  gPad->SaveAs(selName+"_hisdEdxMom.eps");
}


void DrawNormSignalTypeBG(){
  //
  //
  //
  //
  TObjArray array;
  TH2F *hisdEdxNormBG_part[10];
  TH1  *hisdEdxNormBG_partProfile[10];
  //  TH1  *hisdEdxNormBG_partFit[10];
  for (Int_t ipart=0;ipart<10;ipart++){  
    fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
    hisdEdxNormBG_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,3);
    if (!hisdEdxNormBG_part[ipart]) continue;
    hisdEdxNormBG_partProfile[ipart]=(TH1F*)(hisdEdxNormBG_part[ipart]->ProfileX(Form("NormSignalTypeBG_%d_%dProf",ipart,version)));
    //    hisdEdxNormBG_partProfile[ipart]->SetName(;
    hisdEdxNormBG_partProfile[ipart]->SetMaximum(1.5);
    hisdEdxNormBG_partProfile[ipart]->SetMinimum(0.5);
    hisdEdxNormBG_partProfile[ipart]->SetXTitle("#beta#gamma");
    hisdEdxNormBG_partProfile[ipart]->SetYTitle("dEdx/MIP");
    hisdEdxNormBG_partProfile[ipart]->Scale(1/fNorm);
    hisdEdxNormBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisdEdxNormBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
    delete hisdEdxNormBG_part[ipart];
  }
  if (gPad) gPad->SetLogx(kTRUE);
  hisdEdxNormBG_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<10;ipart++)  {
    hisdEdxNormBG_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisdEdxNormBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
  gPad->SaveAs(selName+"_hisdEdxNormBG.eps");    
}


void DrawNormSignalTypeMom(){
  //
  //
  //
  //
  //
  TObjArray array;
  TH2F *hisdEdxNormMom_part[10];
  TH1  *hisdEdxNormMom_partProfile[10];
  //  TH1  *hisdEdxNormMom_partFit[10];
  for (Int_t ipart=0;ipart<10;ipart++){  
    fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
    hisdEdxNormMom_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,1);
    if (!hisdEdxNormMom_part[ipart]) continue;
    hisdEdxNormMom_partProfile[ipart]=(TH1F*)(hisdEdxNormMom_part[ipart]->ProfileX(Form("NormSignalTypeMom_%d_%dProf",ipart,version)));
    hisdEdxNormMom_partProfile[ipart]->Scale(1/fNorm);
    hisdEdxNormMom_partProfile[ipart]->SetMaximum(1.5);
    hisdEdxNormMom_partProfile[ipart]->SetMinimum(0.5);
    hisdEdxNormMom_partProfile[ipart]->SetXTitle("p (GeV/c)");
    hisdEdxNormMom_partProfile[ipart]->SetYTitle("dEdx_{rec}/dNdx_{mc}");
    hisdEdxNormMom_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisdEdxNormMom_partProfile[ipart]->SetMarkerStyle(22+ipart);
    delete hisdEdxNormMom_part[ipart];
  }
  if (gPad) gPad->SetLogx(kTRUE);
  hisdEdxNormMom_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<10;ipart++)  {
    hisdEdxNormMom_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisdEdxNormMom_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"_hisdEdxNormMom.gif");
  gPad->SaveAs(selName+"_hisdEdxNormMom.eps");  
}


void DrawNormSignalTypeBB(){
  //
  //
  //
  //
  //
  TObjArray array;
  TH2F *hisdEdxNormBG_part[10];
  TH1  *hisdEdxNormBG_partProfile[10];
  //  TH1  *hisdEdxNormBG_partFit[10];
  Double_t itypes[10];
  Double_t ival[10];
  Double_t ierr[10];
  for (Int_t ipart=0;ipart<10;ipart++){  
    fTPCsignalNorm->GetAxis(0)->SetRange(1+ipart, 1+ipart);
    hisdEdxNormBG_part[ipart] = (TH2F*)fTPCsignalNorm->Projection(5,6);
    if (!hisdEdxNormBG_part[ipart]) continue;
    hisdEdxNormBG_partProfile[ipart]=(TH1F*)(hisdEdxNormBG_part[ipart]->ProfileX(Form("NormSignalTypeBB_%d_%dProf",ipart,version)));
    hisdEdxNormBG_partProfile[ipart]->Scale(1/fNorm);
    hisdEdxNormBG_partProfile[ipart]->SetMaximum(1.5);
    hisdEdxNormBG_partProfile[ipart]->SetMinimum(0.9);
    hisdEdxNormBG_partProfile[ipart]->SetXTitle("dNdx_{mc}");
    hisdEdxNormBG_partProfile[ipart]->SetYTitle("dEdx_{rec}/dNdx_{mc}");
    hisdEdxNormBG_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisdEdxNormBG_partProfile[ipart]->SetMarkerStyle(22+ipart);
    itypes[ipart]=ipart;
    itypes[ipart]=0;
    itypes[ipart]=0;
    if (ipart%5!=0) {
      hisdEdxNormBG_partProfile[ipart]->Fit("pol1","q","",1,2);
      if (hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")){
	ival[ipart]= hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")->GetParameter(1);
	ierr[ipart]= hisdEdxNormBG_partProfile[ipart]->GetFunction("pol1")->GetParError(1);
      }
    }
    printf("%d\t\t%s%f\t%f\n",ipart,AliPID::ParticleName(ipart%5), ival[ipart], ierr[ipart]);
    delete hisdEdxNormBG_part[ipart];
  }
  if (gPad) gPad->SetLogx(kFALSE);
  hisdEdxNormBG_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<10;ipart++)  {
    hisdEdxNormBG_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisdEdxNormBG_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"_hisdEdxNormBG.gif");
  gPad->SaveAs(selName+"_hisdEdxNormBG.eps");  
}


void DrawPIDpt(){
  //
  //
  //
  //
  TObjArray array;
  TH2F *hisPIDMCRCSame_part[10];
  TH1  *hisPIDMCRCSame_partProfile[10];

  //  TH1  *hisPIDMCRCSame_partFit[10];
  for (Int_t ipart=0;ipart<5;ipart++){  
    fTPCr->GetAxis(0)->SetRange(1+ipart,1+ipart);
    fTPCr->GetAxis(6)->SetRange(1+ipart,1+ipart);
    hisPIDMCRCSame_part[ipart] = (TH2F*)fTPCr->Projection(5,1);
    if (!hisPIDMCRCSame_part[ipart]) continue;
    hisPIDMCRCSame_partProfile[ipart]=(hisPIDMCRCSame_part[ipart]->ProfileX(Form("PIDpt_%d_%dProf",ipart,version)));
    hisPIDMCRCSame_partProfile[ipart]->SetDirectory(0);
    hisPIDMCRCSame_partProfile[ipart]->SetMaximum(1.0);
    hisPIDMCRCSame_partProfile[ipart]->SetMinimum(0.0);
    hisPIDMCRCSame_partProfile[ipart]->SetXTitle("p (GeV/c)");
    hisPIDMCRCSame_partProfile[ipart]->SetYTitle("TPCr");
    hisPIDMCRCSame_partProfile[ipart]->SetMarkerColor(kmicolors[ipart%5]);
    hisPIDMCRCSame_partProfile[ipart]->SetMarkerStyle(22+ipart);
    delete hisPIDMCRCSame_part[ipart];
  }
  if (gPad) gPad->SetLogx(kTRUE);
  hisPIDMCRCSame_partProfile[0]->Draw("");
  TLegend * legend = new TLegend(.7,.7, .99, .99, "TPC dEdx");
  for (Int_t ipart=0;ipart<5;ipart++)  {
    hisPIDMCRCSame_partProfile[ipart]->Draw("same");
    legend->AddEntry(hisPIDMCRCSame_partProfile[ipart],AliPID::ParticleName(ipart%5),"p");
  }
  legend->Draw();
  gPad->SaveAs(selName+"pidPIDMCRCSame.gif");
  gPad->SaveAs(selName+"pidPIDMCRCSame.eps");    
}


/*


TObjArray arr;   
TH2F * his = (TH2F*) GetProjection(fTPCsignalNorm,5,1 , 7,7,   0.2,5,  -0.6,0.6,   130,160,0.5,1.3);
his->FitSlicesY(0,0,-1,0,"QNR",&arr)
  TH1F * hsigma = (TH1F*)arr.At(2);
TH1F * hmean = (TH1F*)arr.At(1);
hsigma->Divide(hmean);
hsigma->Scale(100);
hsigma->SetMaximum(10);
hsigma->SetMinimum(0.0);
hsigma->SetMarkerStyle(22);
hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
hsigma->Draw();
gPad->SaveAs("pic/resoldEdx_pt.eps");
gPad->SaveAs("pic/resoldEdx_pt.gif");


*/
/*
  TObjArray arr;   
  TH2F * his = (TH2F*) GetProjection(fTPCsignalNorm,5,2 , 7,7,   0.6,1,  -0.9,0.9,   0,160,0,1.3);
  his->FitSlicesY(0,0,-1,0,"QNR",&arr)
  TH1F * hsigma = (TH1F*)arr.At(2);
  TH1F * hmean = (TH1F*)arr.At(1);
  hsigma->Divide(hmean);
  hsigma->Scale(100);
  hsigma->SetMaximum(10);
  hsigma->SetMinimum(0.0);
  hsigma->SetMarkerStyle(22);
  hsigma->SetYTitle("#sigma_{dEdx}/dEdx (%)");
  hsigma->Draw();
  gPad->SaveAs("pic/resoldEdx_theta.eps");
  gPad->SaveAs("pic/resoldEdx_theta.gif");

*/
Commit	Line	Data
e7babc5b	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
2bfe5463	24	.L $ALICE_ROOT/PWGPP/Macros/tpcPIDResol.C+
e7babc5b	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	*/