[u/mrichter/AliRoot.git] / TPC / macros / data2011 / makeBBFit.C

/*
   Macro to make the MC based dEdx fits, and study the influence of thre track selection to the PID 
   
   .x $HOME/rootlogon.C
   .x $HOME/NimStyle.C
    .L $ALICE_ROOT/TPC/macros/data2011/makeBBFit.C+
    Init();
    MakeESDCutsPID();
    makeBBfit(10000000);

*/

#include "TCut.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TDatabasePDG.h"
#include "TStatToolkit.h"
#include "TGraphErrors.h"
#include "TStopwatch.h"
#include "TLegend.h"
#include "TLatex.h"
//
#include "TTreeStream.h"
#include "AliTPCParam.h"
#include "AliTPCcalibBase.h"

//
// Global variables
//

TCut cutGeom, cutNcr, cutNcl;
TCut cutFiducial;
Int_t pdgs[4]={11,211,321,2212};
Double_t massPDG[4]={0};
const char *partName[4]={"Electron","Pion","Kaon","Proton"};
const Int_t colors[5]={kGreen,kBlack,kRed,kBlue,5};
const Int_t markers[5]={24,20,21,25,25};
const Float_t markerSize[5]={1.2,1.2,1.0,1.2,1};

void Init();            // init (ALICE) BB parameterization

void FitTransferFunctionAllTheta();
void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream);


void Init(){
   //
  // 1.) Register default ALICE parametrization
  //
  AliTPCParam param;
  param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
  for (Int_t ihis=0; ihis<4; ihis++)     massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
}


void MakeESDCutsPID(){
  //
  // Cuts to be used in the sequenece
  // 1.)  cutGeom - stronger cut on geometry  - perfect agreement data-MC        
  //              - Length factor 1
  // 2.)  cutNcr  - relativally stong cut on the number of crossed rows to gaurantee tracking performance 
  //                relativally good agrement in the MC except of the week decay propability
  //              - Length factor 0.8
  // 3.)  cutNcl  - very week cut on the numebr of clusters
  //              - Length factor 0.6 
  //  
  cutGeom="esdTrack.GetLengthInActiveZone(0,3,236, -5 ,0,0)> 1.*(130-5*abs(esdTrack.fP[4]))"; 
  //  Geomtrical cut in fiducial volume - proper description  in the MC
  //  GetLengthInActiveZone(Int_t mode, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi = 0, TTreeSRedirector* pcstream = 0) const
  //  mode   ==0     - inner param used
  //  deltaY ==3     - cut on edges of sector 
  //                   a.) to avoid dead zone - bias for tracking
  //                   b.) zone with lower Q qualit
  //                   c.) non homogenous Q sample - preferable IROC is skipped -bais in the dEdx
  //  deltaZ ==236   - 
  //  bz = 5 KG      - proper sign should be used ohterwise wrong calculation
  cutNcr="esdTrack.GetTPCClusterInfo(3,1,0,159,1)>0.85*(130-5*abs(esdTrack.fP[4]))";
  //
  // Cut on the number of crossed raws
  // GetTPCClusterInfo(Int_t nNeighbours = 3, Int_t type = 0, Int_t row0 = 0, Int_t row1 = 159, Int_t bitType = 0) 
  // pad row is decalared as crossed if some clusters was found in small neighborhood +- nNeighbours
  //    nNeighbours =3  - +-3 padrows used to define the crossed rows
  //    type = 0        - return number of rows (type==1 returns fraction of clusters)
  //    row0-row1       - upper and lower part of integration
  //  
  cutNcl="esdTrack.fTPCncls>0.7*(130-5*abs(esdTrack.fP[4]))";
  //
  // cut un the number of clusters
  //
  cutFiducial="abs(esdTrack.fP[3])<1&&abs(esdTrack.fD)<1";
}

void makeBBfit(Int_t ntracks=100000){
  //
  // Make BB Bloch fits in bins of the mo
  //
  TFile * ff = TFile::Open("Filtered.root");
  TTreeSRedirector *pcstream = new TTreeSRedirector("dedxFit.root","update");
  TTree * treeMC = (TTree*)ff->Get("highPt");
  //  treeMC->SetCacheSize(6000000000);
  //
  // 1.) Register default ALICE parametrization
  //
  AliTPCParam param;
  param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
  TH3D * hisdEdx3D[4]={0};
  //
  for (Int_t ihis=0; ihis<4; ihis++) {
    massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
    treeMC->SetAlias(TString::Format("cut%s",partName[ihis]), TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]));
    treeMC->SetAlias(TString::Format("dEdxp%s",partName[ihis]), TString::Format("AliTPCParam::BetheBlochAleph(esdTrack.fIp.P()/%f,1)",massPDG[ihis]));
  }
  //
  // 2.) Fill dEdx histograms if not done before
  //

  if (pcstream->GetFile()->Get("RatioP_QMax0Pion3D")==0){    
    //
    TStopwatch timer;
    for (Int_t iDetType=0; iDetType<9; iDetType++){
      for (Int_t ihis=0; ihis<4; ihis++){
	printf("%d\t%d\n",iDetType,ihis);
	timer.Print();
	TString detType="All";
	TString dedx   ="esdTrack.fTPCsignal";
	if (iDetType>0){
	  detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
	  if (iDetType<5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalTot(%d)",(iDetType-1)%4);
	  if (iDetType>5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalMax(%d)",(iDetType-1)%4);
	}  

	TCut cutPDG=TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]).Data();
	TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);	
	TString query= TString::Format("%s/AliTPCParam::BetheBlochAleph(esdTrack.fIp.P()/%f,1):abs(esdTrack.fP[3]):esdTrack.fIp.P()>>%s",dedx.Data(), massPDG[ihis],hname.Data());
	hisdEdx3D[ihis]  = new TH3D(hname.Data(),hname.Data(), 50, 0.2,25, 10,0,1, 100,20,80);
	AliTPCcalibBase::BinLogX(hisdEdx3D[ihis]->GetXaxis());
	treeMC->Draw(query,cutFiducial+cutGeom+cutNcr+cutNcl+cutPDG,"goff",ntracks);
	
	hisdEdx3D[ihis]->GetXaxis()->SetTitle("p (GeV/c)");
	hisdEdx3D[ihis]->GetYaxis()->SetTitle("dEdx/dEdx_{BB} (a.u.)");
	pcstream->GetFile()->cd();
	hisdEdx3D[ihis]->Write(hname.Data());
      }
    }
  }
  delete pcstream;
  //
  // Fit histograms
  //
  pcstream = new TTreeSRedirector("dedxFit.root","update");
  TF1 fg("fg","gaus");
  //
  for (Int_t iDetType=0; iDetType<9; iDetType++){
    TString detType="All";
    if (iDetType>0){
      detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
    }    
    for (Int_t ihis=0; ihis<4; ihis++){ 
      TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);	
      hisdEdx3D[ihis] = (TH3D*)pcstream->GetFile()->Get(hname.Data()); 
      Int_t nbinsP =  hisdEdx3D[0]->GetXaxis()->GetNbins();
      Int_t nbinsTgl =  hisdEdx3D[0]->GetYaxis()->GetNbins();
      //
      for (Int_t ibinP=2; ibinP<nbinsP; ibinP++){
	for (Int_t ibinTgl=2; ibinTgl<nbinsTgl; ibinTgl++){
	  //
	  Double_t pCenter =  hisdEdx3D[0]->GetXaxis()->GetBinCenter(ibinP);
	  Double_t tglCenter =  hisdEdx3D[0]->GetYaxis()->GetBinCenter(ibinTgl);
	  TH1D * hisProj =hisdEdx3D[ihis]->ProjectionZ("xxx", ibinP-1,ibinP+1, ibinTgl-1,ibinTgl+1);
	  Double_t entries = hisProj->GetEntries();
	  if (entries<10) continue;
	  Double_t mean = hisProj->GetMean();
	  Double_t rms = hisProj->GetRMS();
	  hisProj->Fit(&fg,"","");
	  TVectorD vecFit(3, fg.GetParameters());
	  TVectorD vecFitErr(3, fg.GetParErrors());
	  Double_t chi2=fg.GetChisquare();
	  Double_t mass    = massPDG[ihis];
	  Double_t dEdxExp =  AliTPCParam::BetheBlochAleph(pCenter/mass,1);
	  Bool_t isAll=(iDetType==0);
	  Bool_t isMax=(iDetType-1)/4>0;
	  Bool_t isTot=(iDetType-1)/4==0;
	  Int_t dType=(iDetType-1)%4;
	  Double_t dEdx = mean*dEdxExp;
	  //if (
	  (*pcstream)<<"fitdEdxG"<<
	    // dEdx type
	    "iDet="<<iDetType<<      // detector internal number
	    "isAll="<<isAll<<        // full TPC used?
	    "isMax="<<isMax<<        // Qmax charge used ?
	    "isTot="<<isTot<<        // Qtot charge used?
	    "dType="<<dType<<        // Detector region 0-IROC, 1- OROCmedium, 2- OROClong, 3- OROCboth
	    "pType="<<ihis<<         // particle type
	    //
	    "entries="<<entries<<      // entries in histogram
	    "ibinTgl="<<ibinTgl<<      //
	    "ibinP="<<ibinP<<          // 
	    "pCenter="<<pCenter<<      // momentum of center bin
	    "tglCenter="<<tglCenter<<  // tangent lambda
	    //
	    "mass="<<mass<<             // particle mass
	    "dEdxExp="<<dEdxExp<<       // mean expected dEdx in bin 
	    "dEdx="<<dEdx<<             // mean measured dEdx in bin
	    "mean="<<mean<<             // mean measured/expected
	    "rms="<<rms<<               // 
	    "chi2="<<chi2<<             // chi2 of the gausian fit
	    "vecFit.="<<&vecFit<<       // gaus fit param
	    "vecFitErr.="<<&vecFitErr<< // gaus fit error
	    "\n";
	}
      }
    }
  }
  delete pcstream;
  //
  //
  //
}

void FitTransferFunctionAllTheta(){

}


void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream){
  //
  //
  // dEdx_{BB}= f_{tr}(dEdx_{rec}, #phi,#eta)
  //
  // Fit the parematers of transfer function
  // 4 models of transfer function used:
  // 
  //
  // Example usage:
  //    FitTransferFunctionTheta(1,3,8,1,5,0) 
  //    isMax=1; dType=1; tgl=5; skipKaon=0; Double_t maxP=10
  //
  if (!pcstream)  pcstream = new TTreeSRedirector("dedxFit.root","update");
  //
  TTree* treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
  treeFit->SetMarkerStyle(25);
  const char *chfitType[4] = { "dEdx_{BB}=k(#theta)dEdx_{rec}", 
			       "dEdx_{BB}=k(#theta)dEdx_{rec}+#Delta(#theta)", 
			       "dEdx_{BB}=k(#theta,#phi)dEdx_{rec}+#Delta(#theta,#phi)",
			       "dEdx_{BB}=k(#theta,#phi,1/Q)dEdx_{rec}+#Delta(#theta,#phi,1/Q)",
  };
  //
  treeFit->SetAlias("isOK","entries>100&&chi2<80");
  treeFit->SetAlias("dEdxM","(vecFit.fElements[1]*dEdxExp)/50.");
  treeFit->SetAlias("dEdxMErr","(vecFitErr.fElements[1]*dEdxExp)/50.");
  //
  //
  Int_t  npointsMax=30000000;
  TStatToolkit toolkit;
  Double_t chi20=0,chi21=0,chi22=0,chi23=0;
  Int_t    npoints=0;
  TVectorD param0,param1,param2,param3;
  TMatrixD covar0,covar1,covar2,covar3;
  TString fstringFast0="";
  fstringFast0+="dEdxM++";
  //
  TString fstringFast="";
  fstringFast+="dEdxM++";
  fstringFast+="(1/pCenter)++";
  fstringFast+="(1/pCenter)^2++";
  fstringFast+="(1/pCenter)*dEdxM++";
  fstringFast+="((1/pCenter)^2)*dEdxM++";
  //
  //
  TCut cutUse=TString::Format("isOK&&isMax==%d&&dType==%d&&ibinTgl==%d",isMax,dType,tgl).Data();
  TCut cutFit=cutUse+TString::Format("pCenter<%f",maxP).Data();
  if (skipKaon) cutFit+="pType!=3";
  TCut cutDraw =  "(ibinP%2)==0";
  TString *strDeltaFit0 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 1);
  TString *strDeltaFit1 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi21,npoints,param1,covar1,-1,0, npointsMax, 0);
  TString *strDeltaFit2 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi22,npoints,param2,covar2,-1,0, npointsMax, 0);
  //
  fstringFast+="(1/pCenter)/dEdxM++";
  fstringFast+="((1/pCenter)^2)/dEdxM++";
  TString *strDeltaFit3 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi23,npoints,param3,covar3,-1,0, npointsMax, 0);

  treeFit->SetAlias("fitdEdx0",strDeltaFit0->Data());
  treeFit->SetAlias("fitdEdx1",strDeltaFit1->Data());
  treeFit->SetAlias("fitdEdx2",strDeltaFit2->Data());
  treeFit->SetAlias("fitdEdx3",strDeltaFit3->Data());
  
  strDeltaFit0->Tokenize("++")->Print();
  strDeltaFit1->Tokenize("++")->Print();
  strDeltaFit2->Tokenize("++")->Print();
  strDeltaFit3->Tokenize("++")->Print();
  //
  TGraphErrors * graphs[4] ={0};
  //  TCanvas *canvasdEdxFit = new TCanvas(TString::Format("canvasdEdxFit%d",fitType),"canvasdEdxFit",800,800);
  TCanvas *canvasdEdxFit = new TCanvas("canvasdEdxFit","canvasdEdxFit",800,800);
  canvasdEdxFit->SetLeftMargin(0.15);
  canvasdEdxFit->SetRightMargin(0.1);
  canvasdEdxFit->SetBottomMargin(0.15);
  canvasdEdxFit->Divide(2,2,0,0);
 
  for (Int_t fitType=0; fitType<4; fitType++){
    canvasdEdxFit->cd(fitType+1);
    TLegend * legendPart = new TLegend(0.16+0.1*(fitType%2==0),0.16-0.14*(fitType<2),0.89,0.4,TString::Format("%s",chfitType[fitType]));
    legendPart->SetTextSize(0.05);
    legendPart->SetBorderSize(0);
    for (Int_t ihis=3; ihis>=0;ihis--){
      gPad->SetLogx(1);
      TString expr= TString::Format("100*(dEdxExp-fitdEdx%d)/dEdxExp:pCenter:100*dEdxMErr",fitType);
      graphs[ihis]= TStatToolkit::MakeGraphErrors(treeFit,expr,cutUse+cutDraw+TString::Format("pType==%d",ihis).Data(), markers[ihis],colors[ihis],markerSize[ihis]);
      graphs[ihis]->SetMinimum(-5);
      graphs[ihis]->SetMaximum(5);
      graphs[ihis]->GetXaxis()->SetTitle("#it{p} (GeV/c)");
      graphs[ihis]->GetXaxis()->SetTitleSize(0.06);
      graphs[ihis]->GetYaxis()->SetTitleSize(0.06);      
      graphs[ihis]->GetYaxis()->SetTitle("#Delta(d#it{E}dx)/d#it{E}dx_{BB}(#beta#gamma) (%)");
      if (ihis==3) graphs[ihis]->Draw("alp");
      graphs[ihis]->Draw("lp");
      legendPart->AddEntry(graphs[ihis],partName[ihis],"p");
    }
    legendPart->Draw();
  }
  TLatex  latexDraw;
  latexDraw.SetTextSize(0.045);

  const char *chDType[4]={"IROC","OROC medium","OROC long","OROC"};
  {
    if (isMax) latexDraw.DrawLatex(1,4,"Q_{Max}");
    if (!isMax) latexDraw.DrawLatex(1,4,"Q_{tot}");
    latexDraw.DrawLatex(1,3.5,chDType[dType%4]);
    latexDraw.DrawLatex(1,3,TString::Format("#Theta=%1.1f",tgl/10.));
    latexDraw.DrawLatex(1,2.5,TString::Format("Fit: p_{t}<%1.1f (GeV/c)",maxP));
    latexDraw.DrawLatex(1,2.,TString::Format("Fit: Skip Kaon=%d",skipKaon));
  }
  //
  //
  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dTheta%d_UseKaon%d_MaxP%1.0f.pdf",isMax,dType,tgl, skipKaon,maxP));
  (*pcstream)<<"fitTheta"<<
    "isMax="<<isMax<<          // switch is Qmax/Qtot used
    "dType="<<dType<<          // detector Type
    "tgl="<<tgl<<              // tgl number 
    "skipKaon="<<skipKaon<<    // Was kaon dEdx used in the calibration?
    "maxP="<<maxP<<            // Maximal p used in the fit 
    // model 0
    "chi20="<<chi20<<          // chi2
    "param0.="<<&param0<<      // parameters
    "covar0.="<<&covar0<<      // covariance
    // model 1
    "chi21="<<chi21<<          // chi2
    "param1.="<<&param1<<      // parameters
    "covar1.="<<&covar1<<      // covariance
    // model 2
    "chi22="<<chi22<<          // chi2
    "param2.="<<&param2<<      // parameters
    "covar2.="<<&covar2<<      // covariance
    // model 3
    "chi23="<<chi23<<          // chi2
    "param3.="<<&param3<<      // parameters
    "covar3.="<<&covar3<<      // covariance
    "\n";

}
Commit	Line	Data
45fe98dd	1	/*
	2	Macro to make the MC based dEdx fits, and study the influence of thre track selection to the PID
	3
	4	.x $HOME/rootlogon.C
	5	.x $HOME/NimStyle.C
	6	.L $ALICE_ROOT/TPC/macros/data2011/makeBBFit.C+
	7	Init();
	8	MakeESDCutsPID();
	9	makeBBfit(10000000);
	10
	11	*/
	12
	13	#include "TCut.h"
	14	#include "TFile.h"
	15	#include "TTree.h"
	16	#include "TH1.h"
	17	#include "TH2.h"
	18	#include "TH3.h"
	19	#include "TF1.h"
	20	#include "TCanvas.h"
	21	#include "TDatabasePDG.h"
	22	#include "TStatToolkit.h"
	23	#include "TGraphErrors.h"
	24	#include "TStopwatch.h"
	25	#include "TLegend.h"
	26	#include "TLatex.h"
	27	//
	28	#include "TTreeStream.h"
	29	#include "AliTPCParam.h"
	30	#include "AliTPCcalibBase.h"
	31
	32	//
	33	// Global variables
	34	//
	35
	36	TCut cutGeom, cutNcr, cutNcl;
	37	TCut cutFiducial;
	38	Int_t pdgs[4]={11,211,321,2212};
	39	Double_t massPDG[4]={0};
	40	const char *partName[4]={"Electron","Pion","Kaon","Proton"};
	41	const Int_t colors[5]={kGreen,kBlack,kRed,kBlue,5};
	42	const Int_t markers[5]={24,20,21,25,25};
	43	const Float_t markerSize[5]={1.2,1.2,1.0,1.2,1};
	44
	45	void Init(); // init (ALICE) BB parameterization
	46
	47	void FitTransferFunctionAllTheta();
	48	void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream);
	49
	50
	51	void Init(){
	52	//
	53	// 1.) Register default ALICE parametrization
	54	//
	55	AliTPCParam param;
	56	param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
	57	for (Int_t ihis=0; ihis<4; ihis++) massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
	58	}
	59
	60
	61	void MakeESDCutsPID(){
	62	//
	63	// Cuts to be used in the sequenece
	64	// 1.) cutGeom - stronger cut on geometry - perfect agreement data-MC
65	// - Length factor 1
66	// 2.) cutNcr - relativally stong cut on the number of crossed rows to gaurantee tracking performance
67	// relativally good agrement in the MC except of the week decay propability
68	// - Length factor 0.8
69	// 3.) cutNcl - very week cut on the numebr of clusters
70	// - Length factor 0.6
71	//
72	cutGeom="esdTrack.GetLengthInActiveZone(0,3,236, -5 ,0,0)> 1.(130-5abs(esdTrack.fP[4]))";
73	// Geomtrical cut in fiducial volume - proper description in the MC
74	// GetLengthInActiveZone(Int_t mode, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi = 0, TTreeSRedirector* pcstream = 0) const
75	// mode ==0 - inner param used
76	// deltaY ==3 - cut on edges of sector
77	// a.) to avoid dead zone - bias for tracking
78	// b.) zone with lower Q qualit
79	// c.) non homogenous Q sample - preferable IROC is skipped -bais in the dEdx
80	// deltaZ ==236 -
81	// bz = 5 KG - proper sign should be used ohterwise wrong calculation
82	cutNcr="esdTrack.GetTPCClusterInfo(3,1,0,159,1)>0.85(130-5abs(esdTrack.fP[4]))";
83	//
84	// Cut on the number of crossed raws
85	// GetTPCClusterInfo(Int_t nNeighbours = 3, Int_t type = 0, Int_t row0 = 0, Int_t row1 = 159, Int_t bitType = 0)
86	// pad row is decalared as crossed if some clusters was found in small neighborhood +- nNeighbours
87	// nNeighbours =3 - +-3 padrows used to define the crossed rows
88	// type = 0 - return number of rows (type==1 returns fraction of clusters)
89	// row0-row1 - upper and lower part of integration
90	//
91	cutNcl="esdTrack.fTPCncls>0.7(130-5abs(esdTrack.fP[4]))";
92	//
93	// cut un the number of clusters
94	//
95	cutFiducial="abs(esdTrack.fP[3])<1&&abs(esdTrack.fD)<1";
96	}
97
98	void makeBBfit(Int_t ntracks=100000){
99	//
100	// Make BB Bloch fits in bins of the mo
101	//
102	TFile * ff = TFile::Open("Filtered.root");
103	TTreeSRedirector *pcstream = new TTreeSRedirector("dedxFit.root","update");
104	TTree * treeMC = (TTree*)ff->Get("highPt");
105	// treeMC->SetCacheSize(6000000000);
106	//
107	// 1.) Register default ALICE parametrization
108	//
109	AliTPCParam param;
110	param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
111	TH3D * hisdEdx3D[4]={0};
112	//
113	for (Int_t ihis=0; ihis<4; ihis++) {
114	massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
115	treeMC->SetAlias(TString::Format("cut%s",partName[ihis]), TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]));
116	treeMC->SetAlias(TString::Format("dEdxp%s",partName[ihis]), TString::Format("AliTPCParam::BetheBlochAleph(esdTrack.fIp.P()/%f,1)",massPDG[ihis]));
117	}
118	//
119	// 2.) Fill dEdx histograms if not done before
120	//
121
122	if (pcstream->GetFile()->Get("RatioP_QMax0Pion3D")==0){
123	//
124	TStopwatch timer;
125	for (Int_t iDetType=0; iDetType<9; iDetType++){
126	for (Int_t ihis=0; ihis<4; ihis++){
127	printf("%d\t%d\n",iDetType,ihis);
128	timer.Print();
129	TString detType="All";
130	TString dedx ="esdTrack.fTPCsignal";
131	if (iDetType>0){
132	detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
133	if (iDetType<5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalTot(%d)",(iDetType-1)%4);
134	if (iDetType>5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalMax(%d)",(iDetType-1)%4);
135	}
136
137	TCut cutPDG=TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]).Data();
138	TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);
139	TString query= TString::Format("%s/AliTPCParam::BetheBlochAleph(esdTrack.fIp.P()/%f,1):abs(esdTrack.fP[3]):esdTrack.fIp.P()>>%s",dedx.Data(), massPDG[ihis],hname.Data());
140	hisdEdx3D[ihis] = new TH3D(hname.Data(),hname.Data(), 50, 0.2,25, 10,0,1, 100,20,80);
141	AliTPCcalibBase::BinLogX(hisdEdx3D[ihis]->GetXaxis());
142	treeMC->Draw(query,cutFiducial+cutGeom+cutNcr+cutNcl+cutPDG,"goff",ntracks);
143
144	hisdEdx3D[ihis]->GetXaxis()->SetTitle("p (GeV/c)");
145	hisdEdx3D[ihis]->GetYaxis()->SetTitle("dEdx/dEdx_{BB} (a.u.)");
146	pcstream->GetFile()->cd();
147	hisdEdx3D[ihis]->Write(hname.Data());
148	}
149	}
150	}
151	delete pcstream;
152	//
153	// Fit histograms
154	//
155	pcstream = new TTreeSRedirector("dedxFit.root","update");
156	TF1 fg("fg","gaus");
157	//
158	for (Int_t iDetType=0; iDetType<9; iDetType++){
159	TString detType="All";
160	if (iDetType>0){
161	detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
162	}
163	for (Int_t ihis=0; ihis<4; ihis++){
164	TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);
165	hisdEdx3D[ihis] = (TH3D*)pcstream->GetFile()->Get(hname.Data());
166	Int_t nbinsP = hisdEdx3D[0]->GetXaxis()->GetNbins();
167	Int_t nbinsTgl = hisdEdx3D[0]->GetYaxis()->GetNbins();
168	//
169	for (Int_t ibinP=2; ibinP<nbinsP; ibinP++){
170	for (Int_t ibinTgl=2; ibinTgl<nbinsTgl; ibinTgl++){
171	//
172	Double_t pCenter = hisdEdx3D[0]->GetXaxis()->GetBinCenter(ibinP);
173	Double_t tglCenter = hisdEdx3D[0]->GetYaxis()->GetBinCenter(ibinTgl);
174	TH1D * hisProj =hisdEdx3D[ihis]->ProjectionZ("xxx", ibinP-1,ibinP+1, ibinTgl-1,ibinTgl+1);
175	Double_t entries = hisProj->GetEntries();
176	if (entries<10) continue;
177	Double_t mean = hisProj->GetMean();
178	Double_t rms = hisProj->GetRMS();
179	hisProj->Fit(&fg,"","");
180	TVectorD vecFit(3, fg.GetParameters());
181	TVectorD vecFitErr(3, fg.GetParErrors());
182	Double_t chi2=fg.GetChisquare();
183	Double_t mass = massPDG[ihis];
184	Double_t dEdxExp = AliTPCParam::BetheBlochAleph(pCenter/mass,1);
185	Bool_t isAll=(iDetType==0);
186	Bool_t isMax=(iDetType-1)/4>0;
187	Bool_t isTot=(iDetType-1)/4==0;
188	Int_t dType=(iDetType-1)%4;
189	Double_t dEdx = mean*dEdxExp;
190	//if (
191	(*pcstream)<<"fitdEdxG"<<
192	// dEdx type
193	"iDet="<<iDetType<< // detector internal number
194	"isAll="<<isAll<< // full TPC used?
195	"isMax="<<isMax<< // Qmax charge used ?
196	"isTot="<<isTot<< // Qtot charge used?
197	"dType="<<dType<< // Detector region 0-IROC, 1- OROCmedium, 2- OROClong, 3- OROCboth
198	"pType="<<ihis<< // particle type
199	//
200	"entries="<<entries<< // entries in histogram
201	"ibinTgl="<<ibinTgl<< //
202	"ibinP="<<ibinP<< //
203	"pCenter="<<pCenter<< // momentum of center bin
204	"tglCenter="<<tglCenter<< // tangent lambda
205	//
206	"mass="<<mass<< // particle mass
207	"dEdxExp="<<dEdxExp<< // mean expected dEdx in bin
208	"dEdx="<<dEdx<< // mean measured dEdx in bin
209	"mean="<<mean<< // mean measured/expected
210	"rms="<<rms<< //
211	"chi2="<<chi2<< // chi2 of the gausian fit
212	"vecFit.="<<&vecFit<< // gaus fit param
213	"vecFitErr.="<<&vecFitErr<< // gaus fit error
214	"\n";
215	}
216	}
217	}
218	}
219	delete pcstream;
220	//
221	//
222	//
223	}
224
225	void FitTransferFunctionAllTheta(){
226
227	}
228
229
230	void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream){
231	//
232	//
233	// dEdx_{BB}= f_{tr}(dEdx_{rec}, #phi,#eta)
234	//
235	// Fit the parematers of transfer function
236	// 4 models of transfer function used:
237	//
238	//
239	// Example usage:
240	// FitTransferFunctionTheta(1,3,8,1,5,0)
241	// isMax=1; dType=1; tgl=5; skipKaon=0; Double_t maxP=10
242	//
243	if (!pcstream) pcstream = new TTreeSRedirector("dedxFit.root","update");
244	//
245	TTree* treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
246	treeFit->SetMarkerStyle(25);
247	const char *chfitType[4] = { "dEdx_{BB}=k(#theta)dEdx_{rec}",
248	"dEdx_{BB}=k(#theta)dEdx_{rec}+#Delta(#theta)",
249	"dEdx_{BB}=k(#theta,#phi)dEdx_{rec}+#Delta(#theta,#phi)",
250	"dEdx_{BB}=k(#theta,#phi,1/Q)dEdx_{rec}+#Delta(#theta,#phi,1/Q)",
251	};
252	//
253	treeFit->SetAlias("isOK","entries>100&&chi2<80");
254	treeFit->SetAlias("dEdxM","(vecFit.fElements[1]*dEdxExp)/50.");
255	treeFit->SetAlias("dEdxMErr","(vecFitErr.fElements[1]*dEdxExp)/50.");
256	//
257	//
258	Int_t npointsMax=30000000;
259	TStatToolkit toolkit;
260	Double_t chi20=0,chi21=0,chi22=0,chi23=0;
261	Int_t npoints=0;
262	TVectorD param0,param1,param2,param3;
263	TMatrixD covar0,covar1,covar2,covar3;
264	TString fstringFast0="";
265	fstringFast0+="dEdxM++";
266	//
267	TString fstringFast="";
268	fstringFast+="dEdxM++";
269	fstringFast+="(1/pCenter)++";
270	fstringFast+="(1/pCenter)^2++";
271	fstringFast+="(1/pCenter)*dEdxM++";
272	fstringFast+="((1/pCenter)^2)*dEdxM++";
273	//
274	//
275	TCut cutUse=TString::Format("isOK&&isMax==%d&&dType==%d&&ibinTgl==%d",isMax,dType,tgl).Data();
276	TCut cutFit=cutUse+TString::Format("pCenter<%f",maxP).Data();
277	if (skipKaon) cutFit+="pType!=3";
278	TCut cutDraw = "(ibinP%2)==0";
279	TString *strDeltaFit0 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 1);
280	TString *strDeltaFit1 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi21,npoints,param1,covar1,-1,0, npointsMax, 0);
281	TString *strDeltaFit2 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi22,npoints,param2,covar2,-1,0, npointsMax, 0);
282	//
283	fstringFast+="(1/pCenter)/dEdxM++";
284	fstringFast+="((1/pCenter)^2)/dEdxM++";
285	TString *strDeltaFit3 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi23,npoints,param3,covar3,-1,0, npointsMax, 0);
286
287	treeFit->SetAlias("fitdEdx0",strDeltaFit0->Data());
288	treeFit->SetAlias("fitdEdx1",strDeltaFit1->Data());
289	treeFit->SetAlias("fitdEdx2",strDeltaFit2->Data());
290	treeFit->SetAlias("fitdEdx3",strDeltaFit3->Data());
291
292	strDeltaFit0->Tokenize("++")->Print();
293	strDeltaFit1->Tokenize("++")->Print();
294	strDeltaFit2->Tokenize("++")->Print();
295	strDeltaFit3->Tokenize("++")->Print();
296	//
297	TGraphErrors * graphs[4] ={0};
298	// TCanvas *canvasdEdxFit = new TCanvas(TString::Format("canvasdEdxFit%d",fitType),"canvasdEdxFit",800,800);
299	TCanvas *canvasdEdxFit = new TCanvas("canvasdEdxFit","canvasdEdxFit",800,800);
300	canvasdEdxFit->SetLeftMargin(0.15);
301	canvasdEdxFit->SetRightMargin(0.1);
302	canvasdEdxFit->SetBottomMargin(0.15);
303	canvasdEdxFit->Divide(2,2,0,0);
304
305	for (Int_t fitType=0; fitType<4; fitType++){
306	canvasdEdxFit->cd(fitType+1);
307	TLegend * legendPart = new TLegend(0.16+0.1(fitType%2==0),0.16-0.14(fitType<2),0.89,0.4,TString::Format("%s",chfitType[fitType]));
308	legendPart->SetTextSize(0.05);
309	legendPart->SetBorderSize(0);
310	for (Int_t ihis=3; ihis>=0;ihis--){
311	gPad->SetLogx(1);
312	TString expr= TString::Format("100(dEdxExp-fitdEdx%d)/dEdxExp:pCenter:100dEdxMErr",fitType);
313	graphs[ihis]= TStatToolkit::MakeGraphErrors(treeFit,expr,cutUse+cutDraw+TString::Format("pType==%d",ihis).Data(), markers[ihis],colors[ihis],markerSize[ihis]);
314	graphs[ihis]->SetMinimum(-5);
315	graphs[ihis]->SetMaximum(5);
316	graphs[ihis]->GetXaxis()->SetTitle("#it{p} (GeV/c)");
317	graphs[ihis]->GetXaxis()->SetTitleSize(0.06);
318	graphs[ihis]->GetYaxis()->SetTitleSize(0.06);
319	graphs[ihis]->GetYaxis()->SetTitle("#Delta(d#it{E}dx)/d#it{E}dx_{BB}(#beta#gamma) (%)");
320	if (ihis==3) graphs[ihis]->Draw("alp");
321	graphs[ihis]->Draw("lp");
322	legendPart->AddEntry(graphs[ihis],partName[ihis],"p");
323	}
324	legendPart->Draw();
325	}
326	TLatex latexDraw;
327	latexDraw.SetTextSize(0.045);
328
329	const char *chDType[4]={"IROC","OROC medium","OROC long","OROC"};
330	{
331	if (isMax) latexDraw.DrawLatex(1,4,"Q_{Max}");
332	if (!isMax) latexDraw.DrawLatex(1,4,"Q_{tot}");
333	latexDraw.DrawLatex(1,3.5,chDType[dType%4]);
334	latexDraw.DrawLatex(1,3,TString::Format("#Theta=%1.1f",tgl/10.));
335	latexDraw.DrawLatex(1,2.5,TString::Format("Fit: p_{t}<%1.1f (GeV/c)",maxP));
336	latexDraw.DrawLatex(1,2.,TString::Format("Fit: Skip Kaon=%d",skipKaon));
337	}
338	//
339	//
340	canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dTheta%d_UseKaon%d_MaxP%1.0f.pdf",isMax,dType,tgl, skipKaon,maxP));
341	(*pcstream)<<"fitTheta"<<
342	"isMax="<<isMax<< // switch is Qmax/Qtot used
343	"dType="<<dType<< // detector Type
344	"tgl="<<tgl<< // tgl number
345	"skipKaon="<<skipKaon<< // Was kaon dEdx used in the calibration?
346	"maxP="<<maxP<< // Maximal p used in the fit
347	// model 0
348	"chi20="<<chi20<< // chi2
349	"param0.="<<&param0<< // parameters
350	"covar0.="<<&covar0<< // covariance
351	// model 1
352	"chi21="<<chi21<< // chi2
353	"param1.="<<&param1<< // parameters
354	"covar1.="<<&covar1<< // covariance
355	// model 2
356	"chi22="<<chi22<< // chi2
357	"param2.="<<&param2<< // parameters
358	"covar2.="<<&covar2<< // covariance
359	// model 3
360	"chi23="<<chi23<< // chi2
361	"param3.="<<&param3<< // parameters
362	"covar3.="<<&covar3<< // covariance
363	"\n";
364
365	}