TPC/macros/data2011/makeBBFit.C

   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-5*abs(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-5*abs(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-5*abs(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:100*dEdxMErr",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 }