makeBBFit.C

Go to the documentation of this file.

#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};
TTree * treeFit=0;

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

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

void Init(){

  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(Double_t fgeom=1, Double_t fcr=0.85, Double_t fcl=0.7 ){

  cutGeom=TString::Format("esdTrack.GetLengthInActiveZone(0,3,236, -5 ,0,0)> %f*(130-5*abs(esdTrack.fP[4]))",fgeom); 
  //  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=TString::Format("esdTrack.GetTPCClusterInfo(3,1,0,159,1)>%f*(130-5*abs(esdTrack.fP[4]))",fcr);
  //
  // 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=TString::Format("esdTrack.fTPCncls>%f*(130-5*abs(esdTrack.fP[4]))",fcl);
  //
  // cut un the number of clusters
  //
  cutFiducial="abs(esdTrack.fP[3])<1&&abs(esdTrack.fD)<1";
}

void makeBBfit(Int_t ntracks=100000){

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

  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
  treeFit=(TTree*)pcstream->GetFile()->Get("fitdEdxG"); 
  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
  treeFit->SetMarkerStyle(25);
  treeFit->SetCacheSize(200000000);
  for (Int_t isMax=0; isMax<=1; isMax++){
    for (Int_t dType=0; dType<4; dType++){      
      for (Int_t skipKaon=0; skipKaon<=1; skipKaon++){
  for (Float_t maxP=3; maxP<21; maxP+=3){
    printf("%d\t%d\t%d\t%f\n",isMax,dType,skipKaon,maxP);
    FitTransferFunctionAll(isMax,dType,skipKaon,maxP,pcstream,1);
  }
      }
    }
  }
  delete pcstream;
}


void FitTransferFunctionScanTheta(){

  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
  treeFit=(TTree*)pcstream->GetFile()->Get("fitdEdxG"); 
  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
  treeFit->SetMarkerStyle(25);
  treeFit->SetCacheSize(200000000);
  for (Int_t isMax=0; isMax<=1; isMax++){
    for (Int_t dType=0; dType<4; dType++){
      
      for (Int_t tgl=2; tgl<10; tgl++){
  for (Int_t skipKaon=0; skipKaon<=1; skipKaon++){
    for (Float_t maxP=3; maxP<21; maxP+=3){
      printf("%d\t%d\t%d\t%d\t%f\n",isMax,dType,tgl,skipKaon,maxP);
      FitTransferFunctionTheta(isMax,dType,tgl,skipKaon,maxP,pcstream,1);
    }
  }
      }
    }
  }
  delete pcstream;
}


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

  if (!pcstream)  pcstream = new TTreeSRedirector("dedxFit.root","update");
  //
  if (!treeFit){
    treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
    treeFit->SetMarkerStyle(25);
    treeFit->SetCacheSize(100000000);
  }
  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++";                  // 1
  fstringFast+="(1/pCenter)++";            // 2
  fstringFast+="(1/pCenter)^2++";          // 3
  fstringFast+="(1/pCenter)*dEdxM++";      // 4
  fstringFast+="((1/pCenter)^2)*dEdxM++";  // 5
  //
  //
  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++";       //6
  fstringFast+="((1/pCenter)^2)/dEdxM++";   //7
  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();
  //
  (*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 
    "npoints="<<npoints<<      // number of points for 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";
  //
  if (!doDraw) return;
  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_SkipKaon%d_MaxP%1.0f.pdf",isMax,dType,tgl, skipKaon,maxP));
  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dTheta%d_SkipKaon%d_MaxP%1.0f.png",isMax,dType,tgl, skipKaon,maxP));

}
//
//
//
void FitTransferFunctionAll(Bool_t isMax, Int_t dType, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw){

  if (!pcstream)  pcstream = new TTreeSRedirector("dedxFit.root","update");
  //
  if (!treeFit){
    treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
    treeFit->SetMarkerStyle(25);
    treeFit->SetCacheSize(100000000);
  }
  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++";
  fstringFast0+="(tglCenter-0.5)++";
  fstringFast0+="(tglCenter-0.5)*dEdxM++";
  //
  TString fstringFast1="";
  fstringFast1+="dEdxM++";
  fstringFast1+="(tglCenter-0.5)++";
  fstringFast1+="(tglCenter-0.5)*dEdxM++";
  fstringFast1+="(tglCenter-0.5)^2++";
  fstringFast1+="((tglCenter-0.5)^2)*dEdxM++";
  //
  TString fstringFast2="";
  fstringFast2+="dEdxM++";                  // 1
  fstringFast2+="(1/pCenter)++";            // 2
  fstringFast2+="(1/pCenter)^2++";          // 3
  fstringFast2+="(1/pCenter)*dEdxM++";      // 4
  fstringFast2+="((1/pCenter)^2)*dEdxM++";  // 5
  //
  fstringFast2+="(tglCenter-0.5)*dEdxM++";                  // 8
  fstringFast2+="(tglCenter-0.5)*(1/pCenter)++";            // 9
  fstringFast2+="(tglCenter-0.5)*(1/pCenter)^2++";          // 10
  fstringFast2+="(tglCenter-0.5)*(1/pCenter)*dEdxM++";      // 11
  fstringFast2+="(tglCenter-0.5)*((1/pCenter)^2)*dEdxM++";  // 12
  //
  //
  fstringFast2+="((tglCenter-0.5)^2)*dEdxM++";                  // 15
  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)++";            // 16
  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)^2++";          // 17
  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)*dEdxM++";      // 18
  fstringFast2+="((tglCenter-0.5)^2)*((1/pCenter)^2)*dEdxM++";  // 19
  TString fstringFast3=fstringFast2;
  //
  fstringFast3+="(1/pCenter)/dEdxM++";       // 6
  fstringFast3+="((1/pCenter)^2)/dEdxM++";   // 7
  fstringFast3+="(tglCenter-0.5)*(1/pCenter)/dEdxM++";       // 13
  fstringFast3+="(tglCenter-0.5)*((1/pCenter)^2)/dEdxM++";   // 14
  fstringFast3+="((tglCenter-0.5)^2)*(1/pCenter)/dEdxM++";       // 20
  fstringFast3+="((tglCenter-0.5)^2)*((1/pCenter)^2)/dEdxM++";   // 21
  //
  //
  //
  TCut cutUse=TString::Format("isOK&&isMax==%d&&dType==%d",isMax,dType).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", fstringFast1.Data(),cutFit, chi21,npoints,param1,covar1,-1,0, npointsMax, 0);
  TString *strDeltaFit2 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast2.Data(),cutFit, chi22,npoints,param2,covar2,-1,0, npointsMax, 0);
  //
  TString *strDeltaFit3 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast3.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();
  //
  (*pcstream)<<"fitAll"<<
    "isMax="<<isMax<<          // switch is Qmax/Qtot used
    "dType="<<dType<<          // detector Type
    "skipKaon="<<skipKaon<<    // Was kaon dEdx used in the calibration?
    "maxP="<<maxP<<            // Maximal p used in the fit 
    "npoints="<<npoints<<      // number of points for 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";
  //
  if (!doDraw) return;
  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,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_%dSkipKaon%d_MaxP%1.0f.pdf",isMax,dType, skipKaon,maxP));
  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dSkipKaon%d_MaxP%1.0f.png",isMax,dType, skipKaon,maxP));

}


void DrawFit(){

  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
  TTree * treeTheta=(TTree*)pcstream->GetFile()->Get("fitTheta"); 
  treeTheta->SetMarkerStyle(25);
  TTree * treeAll=(TTree*)pcstream->GetFile()->Get("fitAll"); 
  treeAll->SetMarkerStyle(25);
  //
  //
  //
  //treeFit->Draw("vecFit.fElements[2]/vecFit.fElements[1]*pow(dEdxExp*sqrt(1+tglCenter**2),0.25):dEdxExp:tglCenter","isTot&&dType==1&&entries>500","colz",1000000);

}