AliRoot Core  ee782a0 (ee782a0)
makeBBFit.C
Go to the documentation of this file.
1 
29 #include "TCut.h"
30 #include "TFile.h"
31 #include "TTree.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TH3.h"
35 #include "TF1.h"
36 #include "TCanvas.h"
37 #include "TDatabasePDG.h"
38 #include "TStatToolkit.h"
39 #include "TGraphErrors.h"
40 #include "TStopwatch.h"
41 #include "TLegend.h"
42 #include "TLatex.h"
43 //
44 #include "TTreeStream.h"
45 #include "AliTPCParam.h"
46 #include "AliTPCcalibBase.h"
47 
48 //
49 // Global variables
50 //
51 
54 Int_t pdgs[4]={11,211,321,2212};
55 Double_t massPDG[4]={0};
56 const char *partName[4]={"Electron","Pion","Kaon","Proton"};
57 const Int_t colors[5]={kGreen,kBlack,kRed,kBlue,5};
58 const Int_t markers[5]={24,20,21,25,25};
59 const Float_t markerSize[5]={1.2,1.2,1.0,1.2,1};
60 TTree * treeFit=0;
61 
62 void Init(); // init (ALICE) BB parameterization
63 
64 //
65 void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw);
66 void FitTransferFunctionAll(Bool_t isMax, Int_t dType, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw);
67 
68 void Init(){
70 
71  AliTPCParam param;
73  for (Int_t ihis=0; ihis<4; ihis++) massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
74 }
75 
76 
77 void MakeESDCutsPID(Double_t fgeom=1, Double_t fcr=0.85, Double_t fcl=0.7 ){
92 
93  cutGeom=TString::Format("esdTrack.GetLengthInActiveZone(0,3,236, -5 ,0,0)> %f*(130-5*abs(esdTrack.fP[4]))",fgeom);
94  // Geomtrical cut in fiducial volume - proper description in the MC
95  // GetLengthInActiveZone(Int_t mode, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi = 0, TTreeSRedirector* pcstream = 0) const
96  // mode ==0 - inner param used
97  // deltaY ==3 - cut on edges of sector
98  // a.) to avoid dead zone - bias for tracking
99  // b.) zone with lower Q qualit
100  // c.) non homogenous Q sample - preferable IROC is skipped -bais in the dEdx
101  // deltaZ ==236 -
102  // bz = 5 KG - proper sign should be used ohterwise wrong calculation
103  cutNcr=TString::Format("esdTrack.GetTPCClusterInfo(3,1,0,159,1)>%f*(130-5*abs(esdTrack.fP[4]))",fcr);
104  //
105  // Cut on the number of crossed raws
106  // GetTPCClusterInfo(Int_t nNeighbours = 3, Int_t type = 0, Int_t row0 = 0, Int_t row1 = 159, Int_t bitType = 0)
107  // pad row is decalared as crossed if some clusters was found in small neighborhood +- nNeighbours
108  // nNeighbours =3 - +-3 padrows used to define the crossed rows
109  // type = 0 - return number of rows (type==1 returns fraction of clusters)
110  // row0-row1 - upper and lower part of integration
111  //
112  cutNcl=TString::Format("esdTrack.fTPCncls>%f*(130-5*abs(esdTrack.fP[4]))",fcl);
113  //
114  // cut un the number of clusters
115  //
116  cutFiducial="abs(esdTrack.fP[3])<1&&abs(esdTrack.fD)<1";
117 }
118 
119 void makeBBfit(Int_t ntracks=100000){
125 
126  TFile * ff = TFile::Open("Filtered.root");
127  TTreeSRedirector *pcstream = new TTreeSRedirector("dedxFit.root","update");
128  TTree * treeMC = (TTree*)ff->Get("highPt");
129  // treeMC->SetCacheSize(6000000000);
130  //
131  // 1.) Register default ALICE parametrization
132  //
133  AliTPCParam param;
134  param.RegisterBBParam(param.GetBetheBlochParamAlice(),1);
135  TH3D * hisdEdx3D[4]={0};
136  //
137  for (Int_t ihis=0; ihis<4; ihis++) {
138  massPDG [ihis]= TDatabasePDG::Instance()->GetParticle(pdgs[ihis])->Mass();
139  treeMC->SetAlias(TString::Format("cut%s",partName[ihis]), TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]));
140  treeMC->SetAlias(TString::Format("dEdxp%s",partName[ihis]), TString::Format("AliTPCParam::BetheBlochAleph(esdTrack.fIp.P()/%f,1)",massPDG[ihis]));
141  }
142  //
143  // 2.) Fill dEdx histograms if not done before
144  //
145 
146  if (pcstream->GetFile()->Get("RatioP_QMax0Pion3D")==0){
147  //
148  TStopwatch timer;
149  for (Int_t iDetType=0; iDetType<9; iDetType++){
150  for (Int_t ihis=0; ihis<4; ihis++){
151  printf("%d\t%d\n",iDetType,ihis);
152  timer.Print();
153  TString detType="All";
154  TString dedx ="esdTrack.fTPCsignal";
155  if (iDetType>0){
156  detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
157  if (iDetType<5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalTot(%d)",(iDetType-1)%4);
158  if (iDetType>5) dedx=TString::Format("esdTrack.fTPCdEdxInfo.GetSignalMax(%d)",(iDetType-1)%4);
159  }
160 
161  TCut cutPDG=TString::Format("abs(particle.fPdgCode)==%d",pdgs[ihis]).Data();
162  TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);
163  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());
164  hisdEdx3D[ihis] = new TH3D(hname.Data(),hname.Data(), 50, 0.2,25, 10,0,1, 100,20,80);
165  AliTPCcalibBase::BinLogX(hisdEdx3D[ihis]->GetXaxis());
166  treeMC->Draw(query,cutFiducial+cutGeom+cutNcr+cutNcl+cutPDG,"goff",ntracks);
167 
168  hisdEdx3D[ihis]->GetXaxis()->SetTitle("p (GeV/c)");
169  hisdEdx3D[ihis]->GetYaxis()->SetTitle("dEdx/dEdx_{BB} (a.u.)");
170  pcstream->GetFile()->cd();
171  hisdEdx3D[ihis]->Write(hname.Data());
172  }
173  }
174  }
175  delete pcstream;
176  //
177  // Fit histograms
178  //
179  pcstream = new TTreeSRedirector("dedxFit.root","update");
180  TF1 fg("fg","gaus");
181  //
182  for (Int_t iDetType=0; iDetType<9; iDetType++){
183  TString detType="All";
184  if (iDetType>0){
185  detType=TString::Format("Q%s%d",(iDetType-1)/4>0?"Max":"Tot",(iDetType-1)%4);
186  }
187  for (Int_t ihis=0; ihis<4; ihis++){
188  TString hname= TString::Format("RatioP_%s%s3D",detType.Data(),partName[ihis]);
189  hisdEdx3D[ihis] = (TH3D*)pcstream->GetFile()->Get(hname.Data());
190  Int_t nbinsP = hisdEdx3D[0]->GetXaxis()->GetNbins();
191  Int_t nbinsTgl = hisdEdx3D[0]->GetYaxis()->GetNbins();
192  //
193  for (Int_t ibinP=2; ibinP<nbinsP; ibinP++){
194  for (Int_t ibinTgl=2; ibinTgl<nbinsTgl; ibinTgl++){
195  //
196  Double_t pCenter = hisdEdx3D[0]->GetXaxis()->GetBinCenter(ibinP);
197  Double_t tglCenter = hisdEdx3D[0]->GetYaxis()->GetBinCenter(ibinTgl);
198  TH1D * hisProj =hisdEdx3D[ihis]->ProjectionZ("xxx", ibinP-1,ibinP+1, ibinTgl-1,ibinTgl+1);
199  Double_t entries = hisProj->GetEntries();
200  if (entries<10) continue;
201  Double_t mean = hisProj->GetMean();
202  Double_t rms = hisProj->GetRMS();
203  hisProj->Fit(&fg,"","");
204  TVectorD vecFit(3, fg.GetParameters());
205  TVectorD vecFitErr(3, fg.GetParErrors());
206  Double_t chi2=fg.GetChisquare();
207  Double_t mass = massPDG[ihis];
208  Double_t dEdxExp = AliTPCParam::BetheBlochAleph(pCenter/mass,1);
209  Bool_t isAll=(iDetType==0);
210  Bool_t isMax=(iDetType-1)/4>0;
211  Bool_t isTot=(iDetType-1)/4==0;
212  Int_t dType=(iDetType-1)%4;
213  Double_t dEdx = mean*dEdxExp;
214  //if (
215  (*pcstream)<<"fitdEdxG"<<
216  // dEdx type
217  "iDet="<<iDetType<< // detector internal number
218  "isAll="<<isAll<< // full TPC used?
219  "isMax="<<isMax<< // Qmax charge used ?
220  "isTot="<<isTot<< // Qtot charge used?
221  "dType="<<dType<< // Detector region 0-IROC, 1- OROCmedium, 2- OROClong, 3- OROCboth
222  "pType="<<ihis<< // particle type
223  //
224  "entries="<<entries<< // entries in histogram
225  "ibinTgl="<<ibinTgl<< //
226  "ibinP="<<ibinP<< //
227  "pCenter="<<pCenter<< // momentum of center bin
228  "tglCenter="<<tglCenter<< // tangent lambda
229  //
230  "mass="<<mass<< // particle mass
231  "dEdxExp="<<dEdxExp<< // mean expected dEdx in bin
232  "dEdx="<<dEdx<< // mean measured dEdx in bin
233  "mean="<<mean<< // mean measured/expected
234  "rms="<<rms<< //
235  "chi2="<<chi2<< // chi2 of the gausian fit
236  "vecFit.="<<&vecFit<< // gaus fit param
237  "vecFitErr.="<<&vecFitErr<< // gaus fit error
238  "\n";
239  }
240  }
241  }
242  }
243  delete pcstream;
244  //
245  //
246  //
247 }
248 
249 
252 
253  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
254  treeFit=(TTree*)pcstream->GetFile()->Get("fitdEdxG");
255  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
256  treeFit->SetMarkerStyle(25);
257  treeFit->SetCacheSize(200000000);
258  for (Int_t isMax=0; isMax<=1; isMax++){
259  for (Int_t dType=0; dType<4; dType++){
260  for (Int_t skipKaon=0; skipKaon<=1; skipKaon++){
261  for (Float_t maxP=3; maxP<21; maxP+=3){
262  printf("%d\t%d\t%d\t%f\n",isMax,dType,skipKaon,maxP);
263  FitTransferFunctionAll(isMax,dType,skipKaon,maxP,pcstream,1);
264  }
265  }
266  }
267  }
268  delete pcstream;
269 }
270 
271 
274 
275  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
276  treeFit=(TTree*)pcstream->GetFile()->Get("fitdEdxG");
277  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
278  treeFit->SetMarkerStyle(25);
279  treeFit->SetCacheSize(200000000);
280  for (Int_t isMax=0; isMax<=1; isMax++){
281  for (Int_t dType=0; dType<4; dType++){
282 
283  for (Int_t tgl=2; tgl<10; tgl++){
284  for (Int_t skipKaon=0; skipKaon<=1; skipKaon++){
285  for (Float_t maxP=3; maxP<21; maxP+=3){
286  printf("%d\t%d\t%d\t%d\t%f\n",isMax,dType,tgl,skipKaon,maxP);
287  FitTransferFunctionTheta(isMax,dType,tgl,skipKaon,maxP,pcstream,1);
288  }
289  }
290  }
291  }
292  }
293  delete pcstream;
294 }
295 
296 
297 void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw){
306 
307  if (!pcstream) pcstream = new TTreeSRedirector("dedxFit.root","update");
308  //
309  if (!treeFit){
310  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
311  treeFit->SetMarkerStyle(25);
312  treeFit->SetCacheSize(100000000);
313  }
314  const char *chfitType[4] = { "dEdx_{BB}=k(#theta)dEdx_{rec}",
315  "dEdx_{BB}=k(#theta)dEdx_{rec}+#Delta(#theta)",
316  "dEdx_{BB}=k(#theta,#phi)dEdx_{rec}+#Delta(#theta,#phi)",
317  "dEdx_{BB}=k(#theta,#phi,1/Q)dEdx_{rec}+#Delta(#theta,#phi,1/Q)",
318  };
319  //
320  treeFit->SetAlias("isOK","entries>100&&chi2<80");
321  treeFit->SetAlias("dEdxM","(vecFit.fElements[1]*dEdxExp)/50.");
322  treeFit->SetAlias("dEdxMErr","(vecFitErr.fElements[1]*dEdxExp)/50.");
323  //
324  //
325  Int_t npointsMax=30000000;
327  Double_t chi20=0,chi21=0,chi22=0,chi23=0;
328  Int_t npoints=0;
329  TVectorD param0,param1,param2,param3;
330  TMatrixD covar0,covar1,covar2,covar3;
331  TString fstringFast0="";
332  fstringFast0+="dEdxM++";
333  //
334  TString fstringFast="";
335  fstringFast+="dEdxM++"; // 1
336  fstringFast+="(1/pCenter)++"; // 2
337  fstringFast+="(1/pCenter)^2++"; // 3
338  fstringFast+="(1/pCenter)*dEdxM++"; // 4
339  fstringFast+="((1/pCenter)^2)*dEdxM++"; // 5
340  //
341  //
342  TCut cutUse=TString::Format("isOK&&isMax==%d&&dType==%d&&ibinTgl==%d",isMax,dType,tgl).Data();
343  TCut cutFit=cutUse+TString::Format("pCenter<%f",maxP).Data();
344  if (skipKaon) cutFit+="pType!=3";
345  TCut cutDraw = "(ibinP%2)==0";
346  TString *strDeltaFit0 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 1);
347  TString *strDeltaFit1 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi21,npoints,param1,covar1,-1,0, npointsMax, 0);
348  TString *strDeltaFit2 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi22,npoints,param2,covar2,-1,0, npointsMax, 0);
349  //
350  fstringFast+="(1/pCenter)/dEdxM++"; //6
351  fstringFast+="((1/pCenter)^2)/dEdxM++"; //7
352  TString *strDeltaFit3 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast.Data(),cutFit, chi23,npoints,param3,covar3,-1,0, npointsMax, 0);
353 
354  treeFit->SetAlias("fitdEdx0",strDeltaFit0->Data());
355  treeFit->SetAlias("fitdEdx1",strDeltaFit1->Data());
356  treeFit->SetAlias("fitdEdx2",strDeltaFit2->Data());
357  treeFit->SetAlias("fitdEdx3",strDeltaFit3->Data());
358 
359  strDeltaFit0->Tokenize("++")->Print();
360  strDeltaFit1->Tokenize("++")->Print();
361  strDeltaFit2->Tokenize("++")->Print();
362  strDeltaFit3->Tokenize("++")->Print();
363  //
364  (*pcstream)<<"fitTheta"<<
365  "isMax="<<isMax<< // switch is Qmax/Qtot used
366  "dType="<<dType<< // detector Type
367  "tgl="<<tgl<< // tgl number
368  "skipKaon="<<skipKaon<< // Was kaon dEdx used in the calibration?
369  "maxP="<<maxP<< // Maximal p used in the fit
370  "npoints="<<npoints<< // number of points for fit
371  // model 0
372  "chi20="<<chi20<< // chi2
373  "param0.="<<&param0<< // parameters
374  "covar0.="<<&covar0<< // covariance
375  // model 1
376  "chi21="<<chi21<< // chi2
377  "param1.="<<&param1<< // parameters
378  "covar1.="<<&covar1<< // covariance
379  // model 2
380  "chi22="<<chi22<< // chi2
381  "param2.="<<&param2<< // parameters
382  "covar2.="<<&covar2<< // covariance
383  // model 3
384  "chi23="<<chi23<< // chi2
385  "param3.="<<&param3<< // parameters
386  "covar3.="<<&covar3<< // covariance
387  "\n";
388  //
389  if (!doDraw) return;
390  TGraphErrors * graphs[4] ={0};
391  // TCanvas *canvasdEdxFit = new TCanvas(TString::Format("canvasdEdxFit%d",fitType),"canvasdEdxFit",800,800);
392  TCanvas *canvasdEdxFit = new TCanvas("canvasdEdxFit","canvasdEdxFit",800,800);
393  canvasdEdxFit->SetLeftMargin(0.15);
394  canvasdEdxFit->SetRightMargin(0.1);
395  canvasdEdxFit->SetBottomMargin(0.15);
396  canvasdEdxFit->Divide(2,2,0,0);
397 
398  for (Int_t fitType=0; fitType<4; fitType++){
399  canvasdEdxFit->cd(fitType+1);
400  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]));
401  legendPart->SetTextSize(0.05);
402  legendPart->SetBorderSize(0);
403  for (Int_t ihis=3; ihis>=0;ihis--){
404  gPad->SetLogx(1);
405  TString expr= TString::Format("100*(dEdxExp-fitdEdx%d)/dEdxExp:pCenter:100*dEdxMErr",fitType);
406  graphs[ihis]= TStatToolkit::MakeGraphErrors(treeFit,expr,cutUse+cutDraw+TString::Format("pType==%d",ihis).Data(), markers[ihis],colors[ihis],markerSize[ihis]);
407  graphs[ihis]->SetMinimum(-5);
408  graphs[ihis]->SetMaximum(5);
409  graphs[ihis]->GetXaxis()->SetTitle("#it{p} (GeV/c)");
410  graphs[ihis]->GetXaxis()->SetTitleSize(0.06);
411  graphs[ihis]->GetYaxis()->SetTitleSize(0.06);
412  graphs[ihis]->GetYaxis()->SetTitle("#Delta(d#it{E}dx)/d#it{E}dx_{BB}(#beta#gamma) (%)");
413  if (ihis==3) graphs[ihis]->Draw("alp");
414  graphs[ihis]->Draw("lp");
415  legendPart->AddEntry(graphs[ihis],partName[ihis],"p");
416  }
417  legendPart->Draw();
418  }
419  TLatex latexDraw;
420  latexDraw.SetTextSize(0.045);
421 
422  const char *chDType[4]={"IROC","OROC medium","OROC long","OROC"};
423  {
424  if (isMax) latexDraw.DrawLatex(1,4,"Q_{Max}");
425  if (!isMax) latexDraw.DrawLatex(1,4,"Q_{tot}");
426  latexDraw.DrawLatex(1,3.5,chDType[dType%4]);
427  latexDraw.DrawLatex(1,3,TString::Format("#Theta=%1.1f",tgl/10.));
428  latexDraw.DrawLatex(1,2.5,TString::Format("Fit: p_{t}<%1.1f (GeV/c)",maxP));
429  latexDraw.DrawLatex(1,2.,TString::Format("Fit: Skip Kaon=%d",skipKaon));
430  }
431  //
432  //
433  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dTheta%d_SkipKaon%d_MaxP%1.0f.pdf",isMax,dType,tgl, skipKaon,maxP));
434  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dTheta%d_SkipKaon%d_MaxP%1.0f.png",isMax,dType,tgl, skipKaon,maxP));
435 
436 }
437 //
438 //
439 //
440 void FitTransferFunctionAll(Bool_t isMax, Int_t dType, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw){
446 
447  if (!pcstream) pcstream = new TTreeSRedirector("dedxFit.root","update");
448  //
449  if (!treeFit){
450  treeFit = (TTree*)pcstream->GetFile()->Get("fitdEdxG");
451  treeFit->SetMarkerStyle(25);
452  treeFit->SetCacheSize(100000000);
453  }
454  const char *chfitType[4] = { "dEdx_{BB}=k(#theta)dEdx_{rec}",
455  "dEdx_{BB}=k(#theta)dEdx_{rec}+#Delta(#theta)",
456  "dEdx_{BB}=k(#theta,#phi)dEdx_{rec}+#Delta(#theta,#phi)",
457  "dEdx_{BB}=k(#theta,#phi,1/Q)dEdx_{rec}+#Delta(#theta,#phi,1/Q)",
458  };
459  //
460  treeFit->SetAlias("isOK","entries>100&&chi2<80");
461  treeFit->SetAlias("dEdxM","(vecFit.fElements[1]*dEdxExp)/50.");
462  treeFit->SetAlias("dEdxMErr","(vecFitErr.fElements[1]*dEdxExp)/50.");
463  //
464  //
465  Int_t npointsMax=30000000;
467  Double_t chi20=0,chi21=0,chi22=0,chi23=0;
468  Int_t npoints=0;
469  TVectorD param0,param1,param2,param3;
470  TMatrixD covar0,covar1,covar2,covar3;
471  //
472  //
473  //
474  TString fstringFast0="";
475  fstringFast0+="dEdxM++";
476  fstringFast0+="(tglCenter-0.5)++";
477  fstringFast0+="(tglCenter-0.5)*dEdxM++";
478  //
479  TString fstringFast1="";
480  fstringFast1+="dEdxM++";
481  fstringFast1+="(tglCenter-0.5)++";
482  fstringFast1+="(tglCenter-0.5)*dEdxM++";
483  fstringFast1+="(tglCenter-0.5)^2++";
484  fstringFast1+="((tglCenter-0.5)^2)*dEdxM++";
485  //
486  TString fstringFast2="";
487  fstringFast2+="dEdxM++"; // 1
488  fstringFast2+="(1/pCenter)++"; // 2
489  fstringFast2+="(1/pCenter)^2++"; // 3
490  fstringFast2+="(1/pCenter)*dEdxM++"; // 4
491  fstringFast2+="((1/pCenter)^2)*dEdxM++"; // 5
492  //
493  fstringFast2+="(tglCenter-0.5)*dEdxM++"; // 8
494  fstringFast2+="(tglCenter-0.5)*(1/pCenter)++"; // 9
495  fstringFast2+="(tglCenter-0.5)*(1/pCenter)^2++"; // 10
496  fstringFast2+="(tglCenter-0.5)*(1/pCenter)*dEdxM++"; // 11
497  fstringFast2+="(tglCenter-0.5)*((1/pCenter)^2)*dEdxM++"; // 12
498  //
499  //
500  fstringFast2+="((tglCenter-0.5)^2)*dEdxM++"; // 15
501  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)++"; // 16
502  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)^2++"; // 17
503  fstringFast2+="((tglCenter-0.5)^2)*(1/pCenter)*dEdxM++"; // 18
504  fstringFast2+="((tglCenter-0.5)^2)*((1/pCenter)^2)*dEdxM++"; // 19
505  TString fstringFast3=fstringFast2;
506  //
507  fstringFast3+="(1/pCenter)/dEdxM++"; // 6
508  fstringFast3+="((1/pCenter)^2)/dEdxM++"; // 7
509  fstringFast3+="(tglCenter-0.5)*(1/pCenter)/dEdxM++"; // 13
510  fstringFast3+="(tglCenter-0.5)*((1/pCenter)^2)/dEdxM++"; // 14
511  fstringFast3+="((tglCenter-0.5)^2)*(1/pCenter)/dEdxM++"; // 20
512  fstringFast3+="((tglCenter-0.5)^2)*((1/pCenter)^2)/dEdxM++"; // 21
513  //
514  //
515  //
516  TCut cutUse=TString::Format("isOK&&isMax==%d&&dType==%d",isMax,dType).Data();
517  TCut cutFit=cutUse+TString::Format("pCenter<%f",maxP).Data();
518  if (skipKaon) cutFit+="pType!=3";
519  TCut cutDraw = "(ibinP%2)==0";
520  TString *strDeltaFit0 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast0.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 1);
521  TString *strDeltaFit1 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast1.Data(),cutFit, chi21,npoints,param1,covar1,-1,0, npointsMax, 0);
522  TString *strDeltaFit2 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast2.Data(),cutFit, chi22,npoints,param2,covar2,-1,0, npointsMax, 0);
523  //
524  TString *strDeltaFit3 = TStatToolkit::FitPlane(treeFit,"dEdxExp:dEdxMErr", fstringFast3.Data(),cutFit, chi23,npoints,param3,covar3,-1,0, npointsMax, 0);
525 
526  treeFit->SetAlias("fitdEdx0",strDeltaFit0->Data());
527  treeFit->SetAlias("fitdEdx1",strDeltaFit1->Data());
528  treeFit->SetAlias("fitdEdx2",strDeltaFit2->Data());
529  treeFit->SetAlias("fitdEdx3",strDeltaFit3->Data());
530 
531  strDeltaFit0->Tokenize("++")->Print();
532  strDeltaFit1->Tokenize("++")->Print();
533  strDeltaFit2->Tokenize("++")->Print();
534  strDeltaFit3->Tokenize("++")->Print();
535  //
536  (*pcstream)<<"fitAll"<<
537  "isMax="<<isMax<< // switch is Qmax/Qtot used
538  "dType="<<dType<< // detector Type
539  "skipKaon="<<skipKaon<< // Was kaon dEdx used in the calibration?
540  "maxP="<<maxP<< // Maximal p used in the fit
541  "npoints="<<npoints<< // number of points for fit
542  // model 0
543  "chi20="<<chi20<< // chi2
544  "param0.="<<&param0<< // parameters
545  "covar0.="<<&covar0<< // covariance
546  // model 1
547  "chi21="<<chi21<< // chi2
548  "param1.="<<&param1<< // parameters
549  "covar1.="<<&covar1<< // covariance
550  // model 2
551  "chi22="<<chi22<< // chi2
552  "param2.="<<&param2<< // parameters
553  "covar2.="<<&covar2<< // covariance
554  // model 3
555  "chi23="<<chi23<< // chi2
556  "param3.="<<&param3<< // parameters
557  "covar3.="<<&covar3<< // covariance
558  "\n";
559  //
560  if (!doDraw) return;
561  TGraphErrors * graphs[4] ={0};
562  // TCanvas *canvasdEdxFit = new TCanvas(TString::Format("canvasdEdxFit%d",fitType),"canvasdEdxFit",800,800);
563  TCanvas *canvasdEdxFit = new TCanvas("canvasdEdxFit","canvasdEdxFit",800,800);
564  canvasdEdxFit->SetLeftMargin(0.15);
565  canvasdEdxFit->SetRightMargin(0.1);
566  canvasdEdxFit->SetBottomMargin(0.15);
567  canvasdEdxFit->Divide(2,2,0,0);
568 
569  for (Int_t fitType=0; fitType<4; fitType++){
570  canvasdEdxFit->cd(fitType+1);
571  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]));
572  legendPart->SetTextSize(0.05);
573  legendPart->SetBorderSize(0);
574  for (Int_t ihis=3; ihis>=0;ihis--){
575  gPad->SetLogx(1);
576  TString expr= TString::Format("100*(dEdxExp-fitdEdx%d)/dEdxExp:pCenter:100*dEdxMErr",fitType);
577  graphs[ihis]= TStatToolkit::MakeGraphErrors(treeFit,expr,cutUse+cutDraw+TString::Format("pType==%d",ihis).Data(), markers[ihis],colors[ihis],markerSize[ihis]);
578  graphs[ihis]->SetMinimum(-5);
579  graphs[ihis]->SetMaximum(5);
580  graphs[ihis]->GetXaxis()->SetTitle("#it{p} (GeV/c)");
581  graphs[ihis]->GetXaxis()->SetTitleSize(0.06);
582  graphs[ihis]->GetYaxis()->SetTitleSize(0.06);
583  graphs[ihis]->GetYaxis()->SetTitle("#Delta(d#it{E}dx)/d#it{E}dx_{BB}(#beta#gamma) (%)");
584  if (ihis==3) graphs[ihis]->Draw("alp");
585  graphs[ihis]->Draw("lp");
586  legendPart->AddEntry(graphs[ihis],partName[ihis],"p");
587  }
588  legendPart->Draw();
589  }
590  TLatex latexDraw;
591  latexDraw.SetTextSize(0.045);
592 
593  const char *chDType[4]={"IROC","OROC medium","OROC long","OROC"};
594  {
595  if (isMax) latexDraw.DrawLatex(1,4,"Q_{Max}");
596  if (!isMax) latexDraw.DrawLatex(1,4,"Q_{tot}");
597  latexDraw.DrawLatex(1,3.5,chDType[dType%4]);
598  latexDraw.DrawLatex(1,2.5,TString::Format("Fit: p_{t}<%1.1f (GeV/c)",maxP));
599  latexDraw.DrawLatex(1,2.,TString::Format("Fit: Skip Kaon=%d",skipKaon));
600  }
601  //
602  //
603  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dSkipKaon%d_MaxP%1.0f.pdf",isMax,dType, skipKaon,maxP));
604  canvasdEdxFit->SaveAs(TString::Format("fitTransfer_Max%d_Det_%dSkipKaon%d_MaxP%1.0f.png",isMax,dType, skipKaon,maxP));
605 
606 }
607 
608 
609 void DrawFit(){
611 
612  TTreeSRedirector * pcstream = new TTreeSRedirector("dedxFit.root","update");
613  TTree * treeTheta=(TTree*)pcstream->GetFile()->Get("fitTheta");
614  treeTheta->SetMarkerStyle(25);
615  TTree * treeAll=(TTree*)pcstream->GetFile()->Get("fitAll");
616  treeAll->SetMarkerStyle(25);
617  //
618  //
619  //
620  //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);
621 
622 }
static Double_t BetheBlochAleph(Double_t bb, Int_t type=0)
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void MakeESDCutsPID(Double_t fgeom=1, Double_t fcr=0.85, Double_t fcl=0.7)
Definition: makeBBFit.C:77
TFile * Open(const char *filename, Long64_t &nevents)
TTree * treeFit
Definition: makeBBFit.C:60
Summary of statistics functions.
Double_t massPDG[4]
Definition: makeBBFit.C:55
TFile * GetFile()
Definition: TTreeStream.h:95
TStopwatch timer
Definition: kNNSpeed.C:15
Manager and of geomety classes for set: TPC.
Definition: AliTPCParam.h:18
static TVectorD * GetBetheBlochParamAlice()
void makeBBfit(Int_t ntracks=100000)
Definition: makeBBFit.C:119
TCut cutFiducial
Definition: makeBBFit.C:53
TTreeSRedirector * pcstream
TGraphErrors * MakeGraphErrors(TTree *tree, const char *expr="Entry", const char *cut="1", Int_t mstyle=25, Int_t mcolor=1, Float_t msize=-1, Float_t offset=0.0, Int_t entries=10000000, Int_t firstEntry=0)
const Float_t markerSize[5]
Definition: makeBBFit.C:59
void Init()
Definition: makeBBFit.C:68
TCut cutNcl
Definition: makeBBFit.C:52
npoints
Definition: driftITSTPC.C:85
static void RegisterBBParam(TVectorD *param, Int_t position)
void DrawFit()
Definition: makeBBFit.C:609
Double_t chi2
Definition: AnalyzeLaser.C:7
TString * FitPlane(TTree *tree, const char *drawCommand, const char *formula, const char *cuts, Double_t &chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac=-1, Int_t start=0, Int_t stop=10000000, Bool_t fix0=kFALSE)
void FitTransferFunctionScanAll()
Definition: makeBBFit.C:250
TStatToolkit toolkit
Definition: gainCalib.C:18
TCut cutGeom
Definition: makeBBFit.C:52
void FitTransferFunctionScanTheta()
Definition: makeBBFit.C:272
void FitTransferFunctionTheta(Bool_t isMax, Int_t dType, Int_t tgl, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw)
Definition: makeBBFit.C:297
const char * partName[4]
Definition: makeBBFit.C:56
const Int_t colors[5]
Definition: makeBBFit.C:57
TCut cutNcr
Definition: makeBBFit.C:52
class TVectorT< Double_t > TVectorD
const Int_t markers[5]
Definition: makeBBFit.C:58
Int_t pdgs[4]
Definition: makeBBFit.C:54
class TMatrixT< Double_t > TMatrixD
void FitTransferFunctionAll(Bool_t isMax, Int_t dType, Bool_t skipKaon, Double_t maxP, TTreeSRedirector *pcstream, Bool_t doDraw)
Definition: makeBBFit.C:440