StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TbyTPlots.C
1 /*
2  .L TbyTPlots.C+
3  Init();
4  */
5 #if !defined(__CINT__) || defined(__MAKECINT__)
6 #include "Riostream.h"
7 #include <stdio.h>
8 #include "TSystem.h"
9 #include "TMath.h"
10 #include "TH1.h"
11 #include "TH2.h"
12 #include "TH3.h"
13 #include "TProfile.h"
14 #include "TStyle.h"
15 #include "TF1.h"
16 #include "TTree.h"
17 #include "TChain.h"
18 #include "TFile.h"
19 #include "TNtuple.h"
20 #include "TCanvas.h"
21 #include "TMinuit.h"
22 #include "TSpectrum.h"
23 #include "TString.h"
24 #include "TLegend.h"
25 #include "TLine.h"
26 #include "TText.h"
27 #include "TROOT.h"
28 #include "TList.h"
29 #include "TPolyMarker.h"
30 #include "TFileSet.h"
31 #include "TDataSetIter.h"
32 #else
33 class TSystem;
34 class TMath;
35 class TH1;
36 class TH2;
37 class TH3;
38 class TProfile;
39 class TStyle;
40 class TF1;
41 class TTree;
42 class TChain;
43 class TFile;
44 class TNtuple;
45 class TCanvas;
46 class TMinuit;
47 class TSpectrum;
48 class TString;
49 class TLegend;
50 class TLine;
51 class TText;
52 class TROOT;
53 class TList;
54 class TPolyMarker;
55 #endif
56 TFile *fOut = 0;
57 TCanvas *c1 = 0;
58 TChain *fChain = 0;
59 const Int_t NfitPtsHist = 7;
60 struct mc_data_t {
61  float oldPtGl;
62  float newPtGl;
63  float oldEtaGl;
64  float newEtaGl;
65  float oldPhiGl;
66  float newPhiGl;
67  float oldPGl;
68  float newPGl;
69  float oldFitPtsGl;
70  float newFitPtsGl;
71  float oldPtPr;
72  float newPtPr;
73  float oldEtaPr;
74  float newEtaPr;
75  float oldPhiPr;
76  float newPhiPr;
77  float oldPPr;
78  float newPPr;
79  float oldFitPtsPr;
80  float newFitPtsPr;
81  float oldDedx;
82  float newDedx;
83  float oldCharge;
84  float newCharge;
85  float maxPing;
86  float Prim;
87  float oldChi2Gl0;
88  float newChi2Gl0;
89  float oldChi2Gl1;
90  float newChi2Gl1;
91  float oldChi2Pr0;
92  float newChi2Pr0;
93  float oldChi2Pr1;
94  float newChi2Pr1;
95  float firstHitsDist;
96  float lastHitsDist;
97  float oldPrimX;
98  float oldPrimY;
99  float oldPrimZ;
100  float newPrimX;
101  float newPrimY;
102  float newPrimZ;
103 };
104 static mc_data_t data;
105 static Float_t refMult;
106 const Char_t *GP[2] = {"Gl","Pr"};
107 const Char_t *GPTitle[2] = {"Global","Primary"};
108 static TString NameEffpT[4] = {"pTnew","pTold","pTOnew","pTOold"};
109 static TString NameEffPhi[4] = {"Phinew","Phiold","PhiOnew","PhiOold"};
110 const Int_t minNFP = 10; // 25; // 10;
111 const Int_t effNFP = 15; // 25; // 15;
112 static TString Old("Old");
113 static TString New("New");
114 const Char_t *charge[2] = {"pos","neg"};
115 TFile *fIn = (TFile *) gDirectory;
116 TString gTitle;
117 TString cTitle;
118 //________________________________________________________________________________
119 void DrawPng(TCanvas *c) {
120  static Int_t nPng = 0;
121  if (! c) return;
122  TString pngName("");
123  c->Update(); pngName = c->GetName();
124  pngName.ReplaceAll(" ","_");
125  pngName.ReplaceAll("(","_");
126  pngName.ReplaceAll(")","_");
127  pngName.ReplaceAll("{","_");
128  pngName.ReplaceAll("}","_");
129  pngName.ReplaceAll("<","lt");
130  pngName.ReplaceAll(">","gt");
131  pngName.ReplaceAll("old_new_","");
132  pngName.ReplaceAll("old_new","");
133  pngName.ReplaceAll("Old_New_","");
134  pngName.ReplaceAll("Old_New","");
135  pngName.ReplaceAll("GeV/c","");
136  pngName.ReplaceAll(".","_");
137  pngName.ReplaceAll("/","_");
138  pngName.ReplaceAll("^","_");
139  pngName.ReplaceAll("__","_");
140  pngName.ReplaceAll("__","_");
141  pngName.ReplaceAll("#","");
142  pngName += ".png";
143  c->SaveAs(pngName);
144  nPng++;
145  cout << "Draw #\t" << nPng << "\t" << pngName << endl;
146 }
147 //________________________________________________________________________________
148 void DrawFitPnts(Int_t ping=0) {
149  // gStyle->SetPalette(1,0);
150  //________________________________________________________________________________
151  // Fit pts correlation
152  //________________________________________________________________________________
153  TString Name("fitPtsHist");
154  if (ping) Name += "Ping";
155  TH2F* fitPtsHist = (TH2F*) fIn->Get(Name);
156  if (! fitPtsHist) return;
157  cTitle += gTitle;
158  cTitle += "Fit Points Corr";
159  if (ping) cTitle += " Matched";
160  c1 = new TCanvas(cTitle,cTitle,400,400); c1->SetLeftMargin(0.14);
161  c1->SetLogz();
162  fitPtsHist->SetStats(0);
163  fitPtsHist->GetXaxis()->SetRange(0,46);
164  // fitPtsHist->GetYaxis()->SetRange(0,45);
165  fitPtsHist->Draw("colz");
166  Int_t nx = fitPtsHist->GetNbinsX();
167  Int_t ny = fitPtsHist->GetNbinsY();
168  Double_t meanX = 0, meanY = 0;
169  Double_t X = 0, Y = 0;
170  for (Int_t ix = 1; ix <= nx; ix++) {
171  Double_t x = fitPtsHist->GetXaxis()->GetBinCenter(ix);
172  for (Int_t iy = 1; iy <= ny; iy++) {
173  Double_t y = fitPtsHist->GetYaxis()->GetBinCenter(iy);
174  Double_t cont = fitPtsHist->GetBinContent(ix,iy);
175  if (x >= effNFP) {X += cont; meanX += x*cont;}
176  if (y >= effNFP) {Y += cont; meanY += y*cont;}
177  }
178  }
179  if (X > 0) meanX /= X;
180  if (Y > 0) meanY /= Y;
181  TLegend* leg2 = new TLegend(0.12,.85,.95,.95);
182  leg2->SetTextSize(0.033);
183  fitPtsHist->SetMarkerStyle(20);
184  TProfile *pfx = fitPtsHist->ProfileX("_pfx",2,-1,"e"); pfx->Draw("same");
185  fitPtsHist->SetMarkerColor(0);
186  TProfile *pfy = fitPtsHist->ProfileY("_pfy",2,-1,"e"); pfy->Draw("same");
187  leg2->AddEntry(pfy,Form("profY, <fit points. %s> = %7.2f for NFP >= %i",New.Data(),meanX,effNFP));
188  leg2->AddEntry(pfx,Form("profX, <fit points. %s> = %7.2f for NFP >= %i",Old.Data(),meanY,effNFP));
189  leg2->Draw();
190  fitPtsHist->SetMarkerColor(1);
191  DrawPng(c1);
192 }
193 //________________________________________________________________________________
194 void DrawRelMomDifNft(Int_t k=0, Int_t kase=0) {
195  //________________________________________________________________________________
196  // Momentum difference, vs fit pts
197  //________________________________________________________________________________
198  TH2F* pTDifNFP[2];
199  if (kase == 0) for (Int_t c = 0; c < 2; c++) pTDifNFP[c] = (TH2F *) fIn->Get(Form("pTDifNFP%s%s",GP[k],charge[c]));
200  else for (Int_t c = 0; c < 2; c++) pTDifNFP[c] = (TH2F *) fIn->Get(Form("pTDifNFP5%s%s",GP[k],charge[c]));
201  if (! pTDifNFP[0] || ! pTDifNFP[1]) return;
202  if (pTDifNFP[0]->GetEntries() <100 || pTDifNFP[1]->GetEntries() <100) return;
203  cTitle = gTitle;
204  cTitle += "P diff vs NoFitPnts for "; cTitle += GPTitle[k];
205  if (kase) cTitle += " with pT > 0.5 GeV/c";
206  c1 = new TCanvas(cTitle,cTitle,400,400); c1->SetLeftMargin(0.14);
207  TString cTitle2(cTitle);
208  cTitle2 += " Shift";
209  TCanvas* c2 = new TCanvas(cTitle2,cTitle2,400,400); c2->SetLeftMargin(0.14);
210  TString cTitle3(cTitle);
211  cTitle3 += " Sigma";
212  TCanvas* c3 = new TCanvas(cTitle3,cTitle3,400,400); c3->SetLeftMargin(0.14);
213  Double_t ymax = pTDifNFP[0]->GetMaximum();
214  if (pTDifNFP[1]->GetMaximum() > ymax) ymax = pTDifNFP[1]->GetMaximum();
215  c1->cd();
216  c1->SetLogy();
217  c1->SetTicks(1,1);
218  TString opt("e");
219  TLegend* leg2 = new TLegend(0.72,.4,1,1);
220  leg2->SetBorderSize(0);
221  leg2->SetFillColor(0);
222  leg2->SetTextSize(0.033);
223  Int_t nx = pTDifNFP[0]->GetNbinsX();
224  // Int_t ny = pTDifNFP->GetNbinsY();
225  TH1 *proj = 0;
226  for (Int_t i = 1; i <= nx; i++) {
227  for (Int_t c = 0; c < 2; c++) {
228  proj = pTDifNFP[c]->ProjectionY(Form("%s_prj_%i",pTDifNFP[c]->GetName(),i),i,i,"e");
229  proj->SetStats(0);
230  proj->SetMaximum(1.1*ymax);
231  proj->SetMinimum(10);
232  proj->SetMarkerColor(i+1);
233  proj->SetMarkerStyle(20+c);
234  proj->SetTitle(pTDifNFP[c]->GetYaxis()->GetTitle());
235  proj->Draw(opt); opt = "esame";
236  Int_t xmin = (Int_t) pTDifNFP[c]->GetXaxis()->GetBinLowEdge(i);
237  Int_t xmax = (Int_t) pTDifNFP[c]->GetXaxis()->GetBinUpEdge(i);
238  leg2->AddEntry(proj,Form("%s %i<Nfp <%i",charge[c],xmin,xmax));
239  cout << i << ")" << charge[c] << " = [" << xmin << "," << xmax << "] Mean, RMS "
240  << proj->GetMean() << ", " << proj->GetRMS() << endl;
241  }
242  }
243  leg2->Draw();
244  c2->cd();
245  TString name;
246  TLegend* leg7 = new TLegend(0.62,.65,.85,.85);
247  leg7->SetBorderSize(0);
248  leg7->SetFillColor(0);
249  leg7->SetTextSize(0.033);
250  TLegend* leg8 = new TLegend(0.62,.65,.85,.85);
251  leg8->SetBorderSize(0);
252  leg8->SetFillColor(0);
253  leg8->SetTextSize(0.033);
254  ymax = 0;
255  Double_t ymin = 0;
256  Double_t smax = 0;
257  for (Int_t c = 0; c < 2; c++) {
258  pTDifNFP[c]->FitSlicesY();
259  name = pTDifNFP[c]->GetName(); name += "_1";
260  TH1 *mu = (TH1 *) fIn->Get(name);
261  TH1 *sigma = (TH1 *) fIn->Get(name);
262  if (mu->GetMaximum() > ymax) ymax = mu->GetMaximum();
263  if (mu->GetMinimum() < ymin) ymin = mu->GetMinimum();
264  if (sigma->GetMaximum() > smax) smax = sigma->GetMaximum();
265  }
266  for (Int_t c = 0; c < 2; c++) {
267  c2->cd();
268  pTDifNFP[c]->FitSlicesY();
269  name = pTDifNFP[c]->GetName(); name += "_1";
270  TH1 *mu = (TH1 *) fIn->Get(name); mu->SetMaximum(1.1*ymax); if (ymin < 0) mu->SetMinimum(1.1*ymin);
271  mu->SetMarkerStyle(20);
272  mu->SetMarkerColor(1+c);
273  // mu->SetTitle("Relative Shift versus no. of New fit points");
274  mu->SetStats(0);
275  leg7->AddEntry(mu,charge[c],"p");
276  TString X(pTDifNFP[c]->GetYaxis()->GetTitle());
277  TString T("Shift:");
278  T += X; T += " "; T += GP[k]; if (kase) T += " with pT > 0.5 GeV/c";
279  mu->SetTitle(T.Data());
280  if (! c) mu->Draw();
281  else mu->Draw("same");
282  c3->cd();
283  name = pTDifNFP[c]->GetName(); name += "_2";
284  TH1 *sigma = (TH1 *) fIn->Get(name); //sigma->SetMaximum(0.1);
285  sigma->SetMarkerStyle(20);
286  sigma->SetMarkerColor(1+c);
287  // sigma->SetTitle("Relative Resolution versus no. of New fit points");
288  sigma->SetStats(0);
289  leg8->AddEntry(sigma,charge[c],"p");
290  T = "Sigma: ";
291  T += X; T += " "; T += GP[k]; if (kase) T += " with pT > 0.5 GeV/c";
292  sigma->SetTitle(T.Data());
293  if (! c) sigma->Draw();
294  else sigma->Draw("same");
295  }
296  c2->cd(); leg7->Draw();
297  c3->cd(); leg8->Draw();
298  DrawPng(c1);
299  DrawPng(c2);
300  DrawPng(c3);
301 }
302 //________________________________________________________________________________
303 void EffRefMult(Int_t k = 0) {
304  //________________________________________________________________________________
305  // Efficiency versus refMult and pT
306  //________________________________________________________________________________
307  TH2F *pTEf[4];
308  Double_t ymax = 0;
309  for (Int_t i = 0; i < 4; i++) {
310  pTEf[i] = (TH2F *) fIn->Get(Form("%s%s",NameEffpT[i].Data(),GP[k]));
311  if (! pTEf[i]) continue;
312  Double_t y = pTEf[i]->GetMaximum();
313  if (y > ymax) ymax = y;
314  }
315  // cout << "ymax\t" << ymax << endl;
316  if (ymax < 100) return;
317  cTitle = gTitle; cTitle += GP[k];
318  TString hTitle("Momentum Dist"); hTitle += " for "; hTitle += GPTitle[k];
319  cTitle += " "; cTitle += hTitle;
320  c1 = new TCanvas(cTitle,cTitle,400,400); c1->SetLeftMargin(0.14);
321  c1->SetLogy();
322  c1->SetTicks(1,1);
323  // gStyle->SetOptTitle(0);
324  // gStyle->SetOptStat(0);
325  TH1F *dummy = c1->DrawFrame(0.0,1,5.0, 10*ymax);
326  dummy->SetTitle(hTitle);
327  dummy->SetXTitle("pT (GeV/c)");
328  c1->Update();
329  TLegend* leg3 = new TLegend(.35,.6,.9,.9);
330  leg3->SetBorderSize(0);
331  leg3->SetFillColor(0);
332  leg3->SetTextSize(0.023);
333  TH1D *proj;
334  for (Int_t i = 0; i < 4; i++) {
335  if (pTEf[i]) {
336  TAxis *yax = pTEf[i]->GetYaxis();
337  Int_t ny = yax->GetNbins();
338  for (Int_t j = 1; j < ny; j++) {
339  proj = pTEf[i]->ProjectionX(Form("%sHi%i",pTEf[i]->GetName(),j),j,j,"e");
340  // cout << proj->GetName() << "\t" << proj->GetEntries() << endl;
341  if (proj->GetEntries() <= 0) continue;
342  proj->SetMarkerStyle(20+i);
343  proj->SetMarkerColor(j);
344  proj->Draw("same");
345  Int_t xmin = (Int_t) yax->GetBinLowEdge(j);
346  Int_t xmax = (Int_t) yax->GetBinUpEdge(j);
347  TString title(proj->GetTitle());
348  title.ReplaceAll(" versus pTnew","");
349  title.ReplaceAll(" versus pTold","");
350  title.ReplaceAll(" for Global","");
351  title.ReplaceAll(" for Primaries","");
352  leg3->AddEntry(proj,Form("%s for %i < Mult < %i",title.Data(),xmin,xmax),"P");
353  }
354  }
355  }
356  leg3->Draw();
357  DrawPng(c1);
358 }
359 //________________________________________________________________________________
360 void DrawEfficiency(Int_t k = 0, Double_t pmax = 2.5) {
361  //________________________________________________________________________________
362  // Efficiency itself
363  //________________________________________________________________________________
364  TH2F *pTEf[4];
365  TH1D *pTEfP[4];
366  for (Int_t i = 0; i < 4; i++) {
367  pTEf[i] = (TH2F *) fIn->Get(Form("%s%s",NameEffpT[i].Data(),GP[k]));
368  if (pTEf[i]) pTEfP[i] = pTEf[i]->ProjectionX(Form("%s_all",pTEf[i]->GetName()),-1,-1,"e");
369  else pTEfP[i] = 0;
370  }
371  Double_t emin = 0.6;
372  if (k == 1) emin = 0.9;
373  cTitle = gTitle;
374  TString hTitle("");
375  hTitle += GPTitle[k]; hTitle +=" track efficiencies vs pT";
376  cTitle += hTitle;
377  TCanvas* c1 = new TCanvas(cTitle,cTitle,400,400);// c1->SetLeftMargin(0.14);
378  // c1->SetTicks(1,1);
379  TH1F* dummyeff= c1->DrawFrame(0.0,emin,pmax, 1.01);
380  dummyeff->SetTitle(hTitle);
381  dummyeff->SetXTitle("pT (GeV/c)");
382  dummyeff->SetYTitle(Form("fit points>=%i rel. efficiency",effNFP));
383  TLegend* leg4 = new TLegend(.2,.15,.5,.4);
384 #if 0
385  leg4->SetBorderSize(0);
386  leg4->SetFillColor(0);
387  leg4->SetTextSize(0.033);
388 #endif
389  TAxis *y = pTEf[0]->GetYaxis();
390  Int_t ny = y->GetNbins();
391  Double_t yMin = y->GetXmin();
392  Double_t yMax = y->GetXmax();
393  Int_t mstep = 2;
394  Int_t m = ny + 1;
395  for (; m >= 0; m -= mstep) {
396  if (m < 0) m = 0;
397  Int_t m1 = -1, m2 = -1;
398  if (m == 0) {}
399  else {
400  m1 = m-mstep+1;
401  m2 = m;
402  }
403  for (Int_t i = 0; i < 4; i++) {
404  pTEfP[i] = pTEf[i]->ProjectionX(Form("%s_all_%i",pTEf[i]->GetName(),m),m1,m2,"e");
405  }
406  if (! pTEfP[0] || pTEfP[0]->GetEntries() < 1e3) continue;
407  if (! pTEfP[1] || pTEfP[1]->GetEntries() < 1e3) continue;
408  TH1D *effNewP = new TH1D(*pTEfP[3]); effNewP->SetName(Form("eff%s%i",New.Data(),m)); effNewP->SetTitle(hTitle);
409  effNewP->Divide(pTEfP[3],pTEfP[1],1.,1.,"b");
410  TH1D *effOldP = new TH1D(*pTEfP[2]); effOldP->SetName(Form("eff%s%i",Old.Data(),m)); effOldP->SetTitle(hTitle);
411  effOldP->Divide(pTEfP[2],pTEfP[0],1.,1.,"b");
412  effOldP->SetMarkerStyle(21);
413  effOldP->SetMarkerColor(1);
414  effNewP->SetMarkerStyle(20);
415  effNewP->SetMarkerColor(1);
416  if (m) {
417  effOldP->SetMarkerColor(m/mstep+2);
418  effNewP->SetMarkerColor(m/mstep+2);
419  }
420  Int_t ymin = (Int_t) pTEf[0]->GetYaxis()->GetBinLowEdge(m-mstep+1);
421  Int_t ymax = (Int_t) pTEf[0]->GetYaxis()->GetBinUpEdge(m);
422  c1->cd();
423  TString NewName;
424  if (! m) NewName = " for All";
425  else {
426  if (ymin < yMin) NewName = Form(" for Mult < %i", ymax);
427  else if (ymax > yMax) NewName = Form(" for Mult > %i", ymin);
428  else NewName = Form(" for %i < Mult < %i", ymin, ymax);
429  }
430  TString NewcName = c1->GetName(); NewcName += NewName;
431  cout << NewcName.Data() << endl;
432  effOldP->Draw("same");
433  effNewP->Draw("same");
434  leg4->AddEntry(effOldP,Form("%s efficiency %s",Old.Data(),NewName.Data()));
435  leg4->AddEntry(effNewP,Form("%s efficiency %s",New.Data(),NewName.Data()));
436  TCanvas* c3 = new TCanvas(NewcName,NewcName,400,400);// c3->SetLeftMargin(0.14);
437  // c3->SetTicks(1,1);
438  TH1F* dummyeff3= c3->DrawFrame(0.0,emin,pmax, 1.01);
439  TString ttitle(hTitle);
440  ttitle += NewName;
441  dummyeff3->SetTitle(ttitle);
442  dummyeff3->SetXTitle("pT (GeV/c)");
443  dummyeff3->SetYTitle(Form("fit points>=%i rel. efficiency",effNFP));
444  effOldP->Draw("same");
445  effNewP->Draw("same");
446  TLegend* leg5 = new TLegend(.2,.15,.8,.35);
447  leg5->AddEntry(effOldP,Form("%s efficiency %s",Old.Data(),NewName.Data()));
448  leg5->AddEntry(effNewP,Form("%s efficiency %s",New.Data(),NewName.Data()));
449  leg5->Draw();
450  DrawPng(c3);
451  }
452  c1->cd();
453  leg4->Draw();
454  DrawPng(c1);
455 }
456 //________________________________________________________________________________
457 void DrawPhiEfficiency(Int_t k = 0) {
458  //________________________________________________________________________________
459  // Efficiency itself
460  //________________________________________________________________________________
461  TH2F *PhiEf[4];
462  TH1D *PhiEfP[4];
463  for (Int_t i = 0; i < 4; i++) {
464  PhiEf[i] = (TH2F *) fIn->Get(Form("%s%s",NameEffPhi[i].Data(),GP[k]));
465  if (PhiEf[i]) PhiEfP[i] = PhiEf[i]->ProjectionX(Form("%s_all",PhiEf[i]->GetName()),-1,-1,"e");
466  else PhiEfP[i] = 0;
467  }
468  Double_t emin = 0.7;
469  if (k == 1) emin = 0.9;
470  cTitle = gTitle;
471  TString hTitle("");
472  hTitle += GPTitle[k]; hTitle +=" track efficiencies vs #phi";
473  cTitle += hTitle;
474  TCanvas* c1 = new TCanvas(cTitle,cTitle,400,400);// c1->SetLeftMargin(0.14);
475  // c1->SetTicks(1,1);
476  TH1F* dummyeff= c1->DrawFrame(-180,emin,180, 1.01);
477  dummyeff->SetTitle(hTitle);
478  dummyeff->SetXTitle("#phi (degree)");
479  dummyeff->SetYTitle(Form("fit points>=%i rel. efficiency",effNFP));
480  TLegend* leg4 = new TLegend(.2,.15,.5,.4);
481 #if 0
482  leg4->SetBorderSize(0);
483  leg4->SetFillColor(0);
484  leg4->SetTextSize(0.033);
485 #endif
486  TAxis *y = PhiEf[0]->GetYaxis();
487  Int_t ny = y->GetNbins();
488  Double_t yMin = y->GetXmin();
489  Double_t yMax = y->GetXmax();
490  Int_t mstep = 2;
491  Int_t m = ny + 1;
492  for (; m >= 0; m -= mstep) {
493  if (m < 0) m = 0;
494  Int_t m1 = -1, m2 = -1;
495  if (m == 0) {}
496  else {
497  m1 = m-mstep+1;
498  m2 = m;
499  }
500  for (Int_t i = 0; i < 4; i++) {
501  PhiEfP[i] = PhiEf[i]->ProjectionX(Form("%s_all_%i",PhiEf[i]->GetName(),m),m1,m2,"e");
502  }
503  if (! PhiEfP[0] || PhiEfP[0]->GetEntries() < 1e3) continue;
504  if (! PhiEfP[1] || PhiEfP[1]->GetEntries() < 1e3) continue;
505  TH1D *effNewP = new TH1D(*PhiEfP[3]); effNewP->SetName(Form("eff%s%i",New.Data(),m)); effNewP->SetTitle(hTitle);
506  effNewP->Divide(PhiEfP[3],PhiEfP[1],1.,1.,"b");
507  TH1D *effOldP = new TH1D(*PhiEfP[2]); effOldP->SetName(Form("eff%s%i",Old.Data(),m)); effOldP->SetTitle(hTitle);
508  effOldP->Divide(PhiEfP[2],PhiEfP[0],1.,1.,"b");
509  effOldP->SetMarkerStyle(21);
510  effOldP->SetMarkerColor(1);
511  effNewP->SetMarkerStyle(20);
512  effNewP->SetMarkerColor(1);
513  if (m) {
514  effOldP->SetMarkerColor(m/mstep+2);
515  effNewP->SetMarkerColor(m/mstep+2);
516  }
517  Int_t ymin = (Int_t) PhiEf[0]->GetYaxis()->GetBinLowEdge(m-mstep+1);
518  Int_t ymax = (Int_t) PhiEf[0]->GetYaxis()->GetBinUpEdge(m);
519  c1->cd();
520  TString NewName;
521  if (! m) NewName = " for All";
522  else {
523  if (ymin < yMin) NewName = Form(" for Mult < %i", ymax);
524  else if (ymax > yMax) NewName = Form(" for Mult > %i", ymin);
525  else NewName = Form(" for %i < Mult < %i", ymin, ymax);
526  }
527  TString NewcName = c1->GetName(); NewcName += NewName;
528  cout << NewcName.Data() << endl;
529  effOldP->Draw("same");
530  effNewP->Draw("same");
531  leg4->AddEntry(effOldP,Form("%s efficiency %s",Old.Data(),NewName.Data()));
532  leg4->AddEntry(effNewP,Form("%s efficiency %s",New.Data(),NewName.Data()));
533  TCanvas* c3 = new TCanvas(NewcName,NewcName,400,400);// c3->SetLeftMargin(0.14);
534  // c3->SetTicks(1,1);
535  TH1F* dummyeff3= c3->DrawFrame(-180,emin,180, 1.01);
536  TString ttitle(hTitle);
537  ttitle += NewName;
538  dummyeff3->SetTitle(ttitle);
539  dummyeff3->SetXTitle("#phi (degrees)");
540  dummyeff3->SetYTitle(Form("fit points>=%i rel. efficiency",effNFP));
541  effOldP->Draw("same");
542  effNewP->Draw("same");
543  TLegend* leg5 = new TLegend(.2,.15,.8,.35);
544  leg5->AddEntry(effOldP,Form("%s efficiency %s",Old.Data(),NewName.Data()));
545  leg5->AddEntry(effNewP,Form("%s efficiency %s",New.Data(),NewName.Data()));
546  leg5->Draw();
547  DrawPng(c3);
548  }
549  c1->cd();
550  leg4->Draw();
551  DrawPng(c1);
552 }
553 //________________________________________________________________________________
554 void DrawEffVsMult(Int_t k = 0) {
555  //________________________________________________________________________________
556  // Efficiency itself
557  //________________________________________________________________________________
558  cTitle = gTitle;
559  TString hTitle("");
560  hTitle += " Efficiencies"; hTitle += " versus Multiplicity for "; hTitle += GPTitle[k];
561  cTitle += hTitle;
562  TCanvas* cnv4 = new TCanvas(cTitle,cTitle,400,400); c1 = cnv4; c1->SetLeftMargin(0.14);
563  c1->SetTicks(1,1);
564  TH2F *pTEf[4];
565  TH1D *pTEfP[4];
566  for (Int_t i = 0; i < 4; i++) {
567  pTEf[i] = (TH2F *) fIn->Get(Form("%s%s",NameEffpT[i].Data(),GP[k]));
568  if (pTEf[i]) pTEfP[i] = pTEf[i]->ProjectionY(Form("%s_all",pTEf[i]->GetName()),-1,-1,"e");
569  else pTEfP[i] = 0;
570  }
571  TH1F* dummyeff= cnv4->DrawFrame(0.0,0.6,700, 1.05);
572  dummyeff->SetTitle(hTitle);
573  dummyeff->SetXTitle("(uncorrected) Multiplicity");
574  dummyeff->SetYTitle(Form("fit points>=%i rel. efficiency",effNFP));
575  TLegend* leg4 = new TLegend(.2,.15,.5,.4);
576 #if 0
577  leg4->SetBorderSize(0);
578  leg4->SetFillColor(0);
579  leg4->SetTextSize(0.033);
580 #endif
581  TAxis *x = pTEf[0]->GetXaxis();
582  Double_t xMin = x->GetXmin();
583  Double_t xMax = x->GetXmax();
584  Double_t pTs[4] = {0.1, 0.2, 0.5, 2.0};
585  for (Int_t m = 0; m < 5; m++) {
586  if (m < 0) m = 0;
587  Int_t m1 = -1, m2 = -1;
588  if (m == 0) {}
589  else {
590  m1 = x->FindBin(pTs[m-1]);
591  if (m < 4) m2 = x->FindBin(pTs[m])-1;
592  }
593  for (Int_t i = 0; i < 4; i++) {
594  pTEfP[i] = pTEf[i]->ProjectionY(Form("%s_all_%i",pTEf[i]->GetName(),m),m1,m2,"e");
595  }
596  if (! pTEfP[0] || pTEfP[0]->GetEntries() < 1e3) continue;
597  if (! pTEfP[1] || pTEfP[1]->GetEntries() < 1e3) continue;
598  TH1D *effNewP = new TH1D(*pTEfP[3]); effNewP->SetName(Form("eff%s%i",New.Data(),m)); effNewP->SetTitle(hTitle);
599  effNewP->Divide(pTEfP[3],pTEfP[1],1.,1.,"b");
600  TH1D *effOldP = new TH1D(*pTEfP[2]); effOldP->SetName(Form("eff%s%i",Old.Data(),m)); effOldP->SetTitle(hTitle);
601  effOldP->Divide(pTEfP[2],pTEfP[0],1.,1.,"b");
602  effOldP->SetMarkerStyle(21);
603  effOldP->SetMarkerColor(1);
604  effNewP->SetMarkerStyle(20);
605  effNewP->SetMarkerColor(1);
606  if (m) {
607  effOldP->SetMarkerColor(m+1);
608  effNewP->SetMarkerColor(m+1);
609  }
610  Double_t xmin = xMin;
611  Double_t xmax = xMax;
612  if (m1 >= 0) xmin = x->GetBinLowEdge(m1);
613  if (m2 >= 0) xmax = x->GetBinUpEdge(m2);
614  effOldP->Draw("same");
615  effNewP->Draw("same");
616  if (! m) {
617  leg4->AddEntry(effOldP,Form("%s efficiency for ALL",Old.Data()));
618  leg4->AddEntry(effNewP,Form("%s efficiency for ALL",New.Data()));
619  } else {
620  if (xmin < xMin) {
621  leg4->AddEntry(effOldP,Form("%s efficiency for pT < %4.1f",Old.Data(), xmax));
622  leg4->AddEntry(effNewP,Form("%s efficiency for pT < %4.1f",New.Data(), xmax));
623  } else if (xmax > xMax) {
624  leg4->AddEntry(effOldP,Form("%s efficiency for pT > %4.1f",Old.Data(), xmin));
625  leg4->AddEntry(effNewP,Form("%s efficiency for pT > %4.1f",New.Data(), xmin));
626  } else {
627  leg4->AddEntry(effOldP,Form("%s efficiency for %4.1f < pT < %4.1f", Old.Data(), xmin, xmax));
628  leg4->AddEntry(effNewP,Form("%s efficiency for %4.1f < pT < %4.1f", New.Data(), xmin, xmax));
629  }
630  }
631  }
632  leg4->Draw();
633  DrawPng(c1);
634 }
635 //________________________________________________________________________________
636 void DrawpTDiff(Int_t k=0, Double_t pmax = 2.5) {
637  //________________________________________________________________________________
638  // Pt difference, vs pt
639  //________________________________________________________________________________
640  cTitle = gTitle;
641  TString hTitle("");
642  hTitle += Form("<p_{T}^{%s} - p_{T}^{%s}> (GeV/c) versus p_{T}",Old.Data(),New.Data()); hTitle += " for "; hTitle += GPTitle[k];
643  cTitle += "<p_{T}^{Old} - p_{T}^{New}> (GeV/c) versus p_{T}"; cTitle += " for "; cTitle += GPTitle[k];
644  TCanvas* ptdiffgrcnv = new TCanvas(cTitle,cTitle,400,400); c1 = ptdiffgrcnv; c1->SetLeftMargin(0.14);
645  TH1F* dummy2 = c1->DrawFrame(0,-.05,pmax, 0.05);
646  dummy2->SetTitle(hTitle);
647  dummy2->SetXTitle("p_{T} (GeV/c)");
648  dummy2->SetYTitle(Form("p_{T}^{%s} - p_{T}^{%s} GeV/c", Old.Data(), New.Data()));
649  dummy2->GetYaxis()->SetTitleOffset(1.75);
650  // c1->SetLeftMargin(0.15);
651  // c1->SetTicks(1,1);
652  TLegend* leg7 = new TLegend(0.22,.65,.45,.85);
653  leg7->SetBorderSize(0);
654  leg7->SetFillColor(0);
655  leg7->SetTextSize(0.033);
656  for (Int_t i = 0; i < 2; i++) {
657  TH2F *pTdiff = (TH2F *) fIn->Get(Form("pTdiff%s%s",charge[i],GP[k]));
658  if (! pTdiff) continue;
659  TProfile *profx = pTdiff->ProfileX();
660  profx->SetMarkerStyle(20+i);
661  profx->SetMarkerColor(i+1);
662  profx->Draw("same");
663  if (i == 0) leg7->AddEntry(profx,"Positives","P");
664  else leg7->AddEntry(profx,"Negatives","P");
665  }
666  leg7->Draw();
667  DrawPng(c1);
668 }
669 //________________________________________________________________________________
670 void DrawPhiDiff(Int_t k=0, const Char_t *opt = "Phi") {
671  //________________________________________________________________________________
672  // Pt difference, vs pt
673  //________________________________________________________________________________
674  cTitle = gTitle;
675  TString hTitle;
676  hTitle += Form("(p_{T}^{%s} - p_{T}^{%s})/p_{T}^{%s} versus #phi_{%s} ", Old.Data(), New.Data(),New.Data(),New.Data()); hTitle += " for "; hTitle += GPTitle[k];
677  if (TString(opt).Contains("Phi5")) hTitle += " with pT > 0.5 GeV/c";
678  if (TString(opt).Contains("EP")) hTitle += " and #eta > 0";
679  if (TString(opt).Contains("EN")) hTitle += " and #eta < 0";
680  cTitle += hTitle;
681  TCanvas* ptdiffgrcnv = new TCanvas(cTitle,cTitle,400,400); c1 = ptdiffgrcnv; c1->SetLeftMargin(0.14);
682  TH1F* dummy2 = c1->DrawFrame(-TMath::Pi(),-.01,TMath::Pi(),.02);
683  dummy2->SetTitle(hTitle);
684  dummy2->SetXTitle("#phi");
685  dummy2->SetYTitle(Form("(p_{T}^{%s} - p_{T}^{%s})/p_{T}^{%s}", Old.Data(), New.Data(), New.Data()));
686  dummy2->GetYaxis()->SetTitleOffset(1.75);
687  // c1->SetLeftMargin(0.15);
688  // c1->SetTicks(1,1);
689  TLegend* leg7 = new TLegend(0.22,.65,.45,.85);
690  leg7->SetBorderSize(0);
691  leg7->SetFillColor(0);
692  leg7->SetTextSize(0.033);
693  TF1 *gPhi = (TF1 *) gROOT->GetFunction("gPhi");
694  if (! gPhi) {
695  gPhi = new TF1("gPhi","([0]+[1]*x)*sin(x+[2])+[3]",-4,4);
696  gPhi->SetParameters(-1.34684e-03,-2.05069e-04,4.76465e-01,4.30224e-03);
697  }
698  for (Int_t i = 0; i < 2; i++) {
699  TH2F *pTdiff = (TH2F *) fIn->Get(Form("%sdiffR%s%s",opt,charge[i],GP[k]));
700  if (! pTdiff) continue;
701  pTdiff->FitSlicesY(0,0,0,100,"qnr");
702  TH1 *profx = (TH1 *) fIn->Get(Form("%sdiffR%s%s_1",opt,charge[i],GP[k]));
703  if (! profx) continue;
704  profx->Fit(gPhi,"0");
705  profx->SetMarkerStyle(20);
706  profx->SetMarkerColor(i+1);
707  TF1 *fu = profx->GetFunction("gPhi");
708  if (fu) {
709  fu->SetLineColor(i+1);
710  fu->Draw("same");
711  }
712  profx->Draw("same");
713  if (i == 0) leg7->AddEntry(profx,"Positives","P");
714  else leg7->AddEntry(profx,"Negatives","P");
715  }
716  leg7->Draw();
717  DrawPng(c1);
718 }
719 //________________________________________________________________________________
720 void DrawRpTDiff(Int_t k = 0, const Char_t *opt="pTdiffR", Double_t pmax = 2.5) {
721  //________________________________________________________________________________
722  // relative Pt difference, vs pt
723  //________________________________________________________________________________
724 
725  cTitle = gTitle;
726  TString yTitle, sTitle;
727  if (TString(opt) == "pTInvdiffR") {
728  cTitle += "Inv pT diff (perc.) vs pT";
729  yTitle = Form("(1/p_{T}^{%s} - 1/p_{T}^{%s})/1/p_{T}^{%s}", Old.Data(), New.Data(), New.Data());
730  sTitle = Form("Sigma of Relative difference 1/pT (%s - %s)/%s vs pT", Old.Data(), New.Data(), New.Data());
731  } else {
732  cTitle += "pT diff (perc.) vs pT";
733  yTitle = Form("(p_{T}^{%s} - p_{T}^{%s})/p_{T}^{%s}", Old.Data(), New.Data(), New.Data());
734  sTitle = Form("Sigma of Relative difference pT (%s - %s)/%s vs pT", Old.Data(), New.Data(), New.Data());
735  }
736  cTitle += " for "; cTitle += GPTitle[k];
737  c1 = new TCanvas(cTitle,cTitle,400,400); c1->SetLeftMargin(0.14);
738  c1->SetTicks(1,1);
739  c1->SetLeftMargin(0.15);
740  TString cTitle2(cTitle);
741  cTitle2 += "Sigma";
742  TCanvas *c2 = new TCanvas(cTitle2,cTitle2,400,400); c2->SetLeftMargin(0.14);
743  c2->SetTicks(1,1);
744  c2->SetLeftMargin(0.15);
745  TString cTitle3(cTitle);
746  c1->cd();
747  TLegend* leg8 = new TLegend(0.22,.65,.45,.85);
748  leg8->SetBorderSize(0);
749  leg8->SetFillColor(0);
750  leg8->SetTextSize(0.033);
751  TLegend* leg9 = new TLegend(0.22,.65,.45,.85);
752  leg9->SetBorderSize(0);
753  leg9->SetFillColor(0);
754  leg9->SetTextSize(0.033);
755  for (Int_t i = 0; i < 2; i++) {
756  TH2F *pTdiffR = (TH2F *) fIn->Get(Form("%s%s%s",opt,charge[i],GP[k]));
757  if (! pTdiffR) {cout << "Histogram " << Form("%s%s%s",opt,charge[i],GP[k]) << " is not found" << endl; continue;}
758  c1->cd(1);
759  pTdiffR->FitSlicesY();
760  TH1 *mu = (TH1 *) fIn->Get(Form("%s%s%s_1",opt,charge[i],GP[k]));
761  if (mu) {
762  c1->cd();
763  mu->SetAxisRange(0,pmax);
764  mu->SetTitle(cTitle);
765  mu->SetXTitle("pT (GeV/c)");
766  mu->SetYTitle(yTitle);
767  mu->GetYaxis()->SetTitleOffset(1.75);
768  mu->SetMinimum(-0.02);
769  mu->SetMaximum(0.03);
770  mu->SetStats(0);
771  mu->SetMarkerStyle(20); mu->SetMarkerColor(i+1);
772  if (i == 0) mu->Draw();
773  else mu->Draw("same");
774  if (i == 0) leg8->AddEntry(mu,"Positives (fit)","P");
775  else leg8->AddEntry(mu,"Negatives (fit)","P");
776  } else {cout << "Histogram " << Form("%s%s%s_1",opt,charge[i],GP[k]) << " is not found" << endl;}
777  TH1 *sigma = (TH1 *) fIn->Get(Form("%s%s%s_2",opt,charge[i],GP[k]));
778  if (sigma) {
779  c2->cd();
780  sigma->SetMarkerStyle(20); sigma->SetMarkerColor(i+1);
781  sigma->SetStats(0);
782  if (i == 0) {
783  sigma->SetAxisRange(0,pmax);
784  sigma->SetTitle(sTitle);
785  sigma->SetXTitle("pT (GeV/c)");
786  TString s("#sigma(");
787  s += yTitle; s+= ")";
788  sigma->SetYTitle(s);
789  sigma->GetYaxis()->SetTitleOffset(1.75);
790  sigma->SetMinimum(0);
791  sigma->SetMaximum(0.10);
792  sigma->Draw();
793  } else
794  sigma->Draw("same");
795  if (i == 0) leg9->AddEntry(mu,"Positives (fit)","P");
796  else leg9->AddEntry(mu,"Negatives (fit)","P");
797 
798  } else {cout << "Histogram " << Form("%s%s%s_2",opt,charge[i],GP[k]) << " is not found" << endl;}
799  }
800  c1->cd();
801  leg8->Draw();
802  c2->cd();
803  leg9->Draw();
804  DrawPng(c1);
805  DrawPng(c2);
806 }
807 //________________________________________________________________________________
808 void Init(const Char_t *file="Plots.root") {
809  TString pwd(gSystem->BaseName( gSystem->WorkingDirectory()));
810  TObjArray *obj = pwd.Tokenize("_");
811  Int_t nParsed = obj->GetEntries();
812  if (nParsed == 2) {
813  Old = ((TObjString *) obj->At(0))->GetName();
814  New = ((TObjString *) obj->At(1))->GetName();
815  NameEffpT[0] = Form("pT%s",New.Data());
816  NameEffpT[1] = Form("pT%s",Old.Data());
817  NameEffpT[2] = Form("pTO%s",New.Data());
818  NameEffpT[3] = Form("pTO%s",Old.Data());
819  NameEffPhi[0] = Form("Phi%s",New.Data());
820  NameEffPhi[1] = Form("Phi%s",Old.Data());
821  NameEffPhi[2] = Form("PhiO%s",New.Data());
822  NameEffPhi[3] = Form("PhiO%s",Old.Data());
823  }
824  gTitle = file;
825  if (gTitle != "") {
826  fIn = new TFile(gTitle);
827  // gTitle.ReplaceAll(".Plots.root","");
828  gTitle = "";
829  }
830 }
831 //________________________________________________________________________________
832 void DrawPrimVx() {
833  const Char_t *XYZ[3] = {"X","Y","Z"};
834  for (Int_t i = 0; i < 3; i++) {
835  TH1* h = (TH1 *) gDirectory->Get(Form("DifPv%s",XYZ[i]));
836  if (! h) continue;
837  TCanvas *c = new TCanvas(Form("PrimaryVertex%sdiff",XYZ[i]),Form("PrimaryVertex%sdiff",XYZ[i]));
838  h->Draw();
839  DrawPng(c);
840  }
841 }
842 //________________________________________________________________________________
843 void DrawCharge() {
844  const Char_t *hNames[2] = {"Charge","Charge15"};
845  // const Char_t *fNames[2] = {Form("%sVs%sChargeAllPr",New.Data(),Old.Data()),Form("%sVs%sChargePrNFPGE15",New.Data(),Old.Data())};
846  const Char_t *fNames[2] = {"NewVsOldChargeAllPr","NewVsOldChargePrNFPGE15"};
847  for (Int_t i = 0; i < 2; i++) {
848  TH2 *h = (TH2 *) gDirectory->Get(hNames[i]);
849  if (! h ) continue;
850  TCanvas *c = new TCanvas(fNames[i],fNames[i]);
851  h->Draw("text");
852  DrawPng(c);
853  }
854 }
855 //________________________________________________________________________________
856 void Draw(const Char_t *file="Plots.root") {
857  Init(file);
858  DrawFitPnts(); DrawFitPnts(1);
859  DrawRelMomDifNft(0,0); DrawRelMomDifNft(1,0);
860  DrawRelMomDifNft(0,5); DrawRelMomDifNft(1,5);
861  EffRefMult(0); EffRefMult(1);
862  DrawEfficiency(0); DrawEfficiency(1);
863  DrawPhiEfficiency(0); DrawPhiEfficiency(1);
864  DrawpTDiff(0); DrawpTDiff(1);
865  DrawRpTDiff(0);DrawRpTDiff(1);
866  DrawRpTDiff(0,"pTInvdiffR"); DrawRpTDiff(1,"pTInvdiffR");
867  DrawPrimVx();
868  DrawCharge();
869 #if 0
870  DrawPhiDiff(0);DrawPhiDiff(1);
871  DrawPhiDiff(0,"Phi5");DrawPhiDiff(1,"Phi5");
872  DrawPhiDiff(0,"PhiEP");DrawPhiDiff(1,"PhiEP");
873  DrawPhiDiff(0,"Phi5EP");DrawPhiDiff(1,"Phi5EP");
874  DrawPhiDiff(0,"PhiEN");DrawPhiDiff(1,"PhiEN");
875  DrawPhiDiff(0,"Phi5EN");DrawPhiDiff(1,"Phi5EN");
876 #endif
877 }
878 //________________________________________________________________________________
879 void TbyTPlots(const Char_t *file = 0, Int_t Nentries=0) {
880  TString Out(file);
881  Init(Out.Data());
882  if (Out == "") {
883  TString TreeName("trackMateComp");
884  fChain = new TChain(TreeName);
885  TFileSet dir(".");
886  TDataSetIter next(&dir,0);
887  TDataSet *set = 0;
888  TFile *f = 0;
889  Int_t FileNo = 0;
890  Int_t Ntotal = 0;
891  TString tFile("");
892  while ((set = next())) {
893  TString Title(set->GetTitle());
894  if (Title != "file") continue;
895  TString Name(set->GetName());
896  if (! Name.EndsWith(".root")) continue;
897  if (! (Name.BeginsWith("trackMateFile") || Name.BeginsWith("TbyT"))) continue;
898  f = new TFile(Name);
899  TChain *tree = (TChain *) f->Get(TreeName);
900  if (tree) {
901  const Char_t *pFile = f->GetName();
902  ULong64_t nEvents = tree->GetEntries();
903  if (nEvents > 0) {
904  FileNo++;
905  Ntotal += nEvents;
906  cout << "\tAdd "<< FileNo << "\t" << nEvents << "\tEvents, Total = " << Ntotal << endl;
907  fChain->Add(pFile);
908  }
909  delete f;
910  }
911  }
912  fChain->SetBranchAddress("data_array",&data.oldPtGl);
913  fChain->SetBranchAddress("ev_array",&refMult);
914  tFile.ReplaceAll(".root","");
915  tFile.ReplaceAll("trackMate_physics_","");
916  tFile.ReplaceAll("trackMate_","");
917  tFile.ReplaceAll("5043073_raw_1010010_","");
918  tFile.ReplaceAll("5046043_raw_1010010_","");
919  Out = tFile;
920  Out += "Plots.root";
921  TFile *fOut = new TFile(Out.Data(),"recreate");
922  cout << "Opened " << Out << endl;
923  // book
924  TString Title;
925  TString Name;
926  TH1D *DifPvX = new TH1D("DifPvX",Form("Difference in X for %s - %s positions",New.Data(),Old.Data()),100,-0.25,0.25);
927  TH1D *DifPvY = new TH1D("DifPvY",Form("Difference in Y for %s - %s positions",New.Data(),Old.Data()),100,-0.25,0.25);
928  TH1D *DifPvZ = new TH1D("DifPvZ",Form("Difference in Z for %s - %s positions",New.Data(),Old.Data()),100,-0.25,0.25);
929  TH2D *Charge = new TH2D("Charge",Form("Charge flip between %s and %s for all primaries",New.Data(),Old.Data()),3,-1.5,1.5,3,-1.5,1.5);
930  Charge->SetXTitle(Form("%s charge",Old.Data()));
931  Charge->SetYTitle(Form("%s charge",New.Data()));
932  TH2D *Charge15 = new TH2D("Charge15",Form("Charge flip between %s and %s for all primaries with NFP >= 15",New.Data(),Old.Data()),3,-1.5,1.5,3,-1.5,1.5);
933  Charge15->SetXTitle(Form("%s charge",Old.Data()));
934  Charge15->SetYTitle(Form("%s charge",New.Data()));
935  //________________________________________________________________________________
936  // Fit pts correlation
937  //________________________________________________________________________________
938  Name = "fitPtsHist", Title = Form("Fit Pts %s vs %s",Old.Data(), New.Data());
939  TH2F* fitPtsHist = new TH2F(Name,Title,50,0,50,50,0,50);
940  fitPtsHist->SetXTitle(Form("fit points, %s",New.Data()));
941  fitPtsHist->SetYTitle(Form("fit points, %s",Old.Data()));
942  Name = "fitPtsHistPing", Title = Form("Fit Pts %s vs %s matched", Old.Data(), New.Data());
943  TH2F* fitPtsHistPing = new TH2F(Name,Title,50,0,50,50,0,50);
944  fitPtsHist->SetXTitle(Form("fit points, %s",New.Data()));
945  fitPtsHist->SetYTitle(Form("fit points, %s",Old.Data()));
946  //________________________________________________________________________________
947  // Relative Momentum difference, vs fit pts
948  //________________________________________________________________________________
949  TH2F* pTDifNFP[2][2];
950  TH2F* pTDifNFP5[2][2];
951  for (Int_t k = 0; k < 2; k++) {
952  for (Int_t c = 0; c < 2; c++) {
953  Name = Form("pTDifNFP%s%s",GP[k],charge[c]);
954  Title = Form("Relative Momentum difference vs no. of %s fit pts for %s tracks",New.Data(),GPTitle[k]);
955  pTDifNFP[k][c] = new TH2F(Name,Title, 7,10.,45.,200,-.2,.2);
956  pTDifNFP[k][c]->SetYTitle(Form("Momentum Difference (pT_%s - pT_%s)/pT_%s",Old.Data(),New.Data(),New.Data()));
957  pTDifNFP[k][c]->SetXTitle(Form("No. %s fit points",New.Data()));
958  pTDifNFP[k][c]->SetMarkerStyle(20);
959  Name = Form("pTDifNFP5%s%s",GP[k],charge[c]);
960  Title = Form("Relative Momentum difference vs no. of %s fit pts for %s tracks with pT > 0.5 GeV/c",New.Data(),GPTitle[k]);
961  pTDifNFP5[k][c] = new TH2F(Name,Title, 7,10.,45.,200,-.2,.2);
962  pTDifNFP5[k][c]->SetYTitle(Form("Momentum Difference (pT_%s - pT_%s)/pT_%s", Old.Data(), New.Data(), New.Data()));
963  pTDifNFP5[k][c]->SetXTitle(Form("No. %s fit points",New.Data()));
964  pTDifNFP5[k][c]->SetMarkerStyle(20);
965  }
966  }
967  //________________________________________________________________________________
968  // Efficiency versus Mult and pT
969  //________________________________________________________________________________
970  //
971  TH2F *pTEf[2][4];
972  const Char_t *TitleEffPt[4] = {Form("%s Fit Pts>%i %s versus pT_%s |#eta| < 0.5",New.Data(),effNFP,New.Data(),New.Data()),
973  Form("%s Fit Pts>%i %s versus pT_%s |#eta| < 0.5",Old.Data(),effNFP,Old.Data(),Old.Data()),
974  Form("%s & %s Fit Pts>%i %s and %s versus pT_%s |#eta| < 0.5",New.Data(),Old.Data(),effNFP,Old.Data(),New.Data(),New.Data()),
975  Form("%s & %s Fit Pts>%i %s and %s versus pT_%s |#eta| < 0.5",New.Data(),Old.Data(),effNFP,Old.Data(),New.Data(),Old.Data())};
976  for (Int_t k = 0; k < 2; k++) {
977  for (Int_t i = 0; i < 4; i++) {
978  Name = Form("%s%s",NameEffpT[i].Data(),GP[k]);
979  Title = Form("%s for %s",TitleEffPt[i],GPTitle[k]);
980  pTEf[k][i] = new TH2F(Name,Title,70,0,7,7,0,700);
981  pTEf[k][i]->Sumw2();
982  pTEf[k][i]->SetMarkerStyle(20+i);
983  pTEf[k][i]->SetMarkerColor(1+i);
984  }
985  }
986  // Efficiency versus Mult and Phi
987  //________________________________________________________________________________
988  //
989  const Char_t *TitleEffPhi[4] = {Form("%s Fit Pts>%i %s versus #phi_%s |#eta| < 0.5",New.Data(),effNFP,New.Data(),New.Data()),
990  Form("%s Fit Pts>%i %s versus #phi_%s |#eta| < 0.5",Old.Data(),effNFP,Old.Data(),Old.Data()),
991  Form("%s & %s Fit Pts>%i %s and %s versus #phi_%s |#eta| < 0.5",New.Data(),Old.Data(),effNFP,Old.Data(),New.Data(),New.Data()),
992  Form("%s & %s Fit Pts>%i %s and %s versus #phi_%s |#eta| < 0.5",New.Data(),Old.Data(),effNFP,Old.Data(),New.Data(),Old.Data())};
993  TH2F *PhiEf[2][4];
994  for (Int_t k = 0; k < 2; k++) {
995  for (Int_t i = 0; i < 4; i++) {
996  Name = Form("%s%s",NameEffPhi[i].Data(),GP[k]);
997  Title = Form("%s for %s",TitleEffPhi[i],GPTitle[k]);
998  PhiEf[k][i] = new TH2F(Name,Title,180,-180,180,7,0,700);
999  PhiEf[k][i]->Sumw2();
1000  PhiEf[k][i]->SetMarkerStyle(20+i);
1001  PhiEf[k][i]->SetMarkerColor(1+i);
1002  }
1003  }
1004  //________________________________________________________________________________
1005  // Pt difference, vs pt
1006  //________________________________________________________________________________
1007  TH2F* pTdiff[2][2];
1008  TH2F* pTdiffR[2][2];
1009  TH2F* pTInvdiffR[2][2];
1010  TH2F* PhidiffR[3][2][2];
1011  TH2F* Phi5diffR[3][2][2];
1012  for (Int_t k = 0; k < 2; k++) {
1013  for (Int_t i = 0; i < 2; i++) {
1014  Name = Form("pTdiff%s%s",charge[i],GP[k]);
1015  Title = Form("pT_%s - pT_%s for %s",Old.Data(),New.Data(),GPTitle[k]);
1016  pTdiff[k][i] = new TH2F(Name,Title,40,0,4,200,-.8,.8);
1017  pTdiff[k][i]->SetXTitle(Form("p_{T} %s (GeV/c)",New.Data()));
1018  TString YTitle = Form("p_{T} Difference (pT_%s - pT_%s),",Old.Data(),New.Data());
1019  TString title("");
1020  if (i == 0)
1021  title += " +";
1022  else
1023  title += " -";
1024  title += "charge";
1025  if (k == 0) title += " for Globals";
1026  else title += " for Primaries";
1027  TString yTitle = YTitle; yTitle += title;
1028  pTdiff[k][i]->SetYTitle(yTitle);
1029  pTdiff[k][i]->Sumw2();
1030  pTdiff[k][i]->SetMarkerStyle(20);
1031  pTdiff[k][i]->SetMarkerColor(i+2);
1032  Name = Form("pTdiffR%s%s",charge[i],GP[k]);
1033  pTdiffR[k][i] = new TH2F(Name,Title,40,0,4,600,-0.15,0.15);
1034  pTdiffR[k][i]->SetXTitle(Form("p_{T} %s (GeV/c)",New.Data()));
1035  YTitle = Form("(pT_%s - pT_%s)/pT_%s,",Old.Data(),New.Data(),New.Data());
1036  yTitle = YTitle; yTitle += title;
1037  pTdiffR[k][i]->SetYTitle(yTitle);
1038  pTdiffR[k][i]->Sumw2();
1039  pTdiffR[k][i]->SetMarkerStyle(20);
1040  pTdiffR[k][i]->SetMarkerColor(i+2);
1041 
1042  Name = Form("pTInvdiffR%s%s",charge[i],GP[k]);
1043  pTInvdiffR[k][i] = new TH2F(Name,Title,40,0.,4.,600,-0.15,0.15);
1044  pTInvdiffR[k][i]->SetXTitle(Form("pT_%s",New.Data()));
1045  YTitle = Form("(1/pT_%s - 1/pT_%s)/(1/pT_%s), versus pT_%s",Old.Data(),New.Data(),New.Data(),New.Data());
1046  yTitle = YTitle; yTitle += title;
1047  pTInvdiffR[k][i]->SetYTitle(yTitle);
1048  pTInvdiffR[k][i]->Sumw2();
1049  pTInvdiffR[k][i]->SetMarkerStyle(20);
1050  pTInvdiffR[k][i]->SetMarkerColor(i+2);
1051  const Char_t *NameEta[3] = {"","EP","EN"};
1052  const Char_t *TitleEta[3] = {"all #eta", "|#eta| > 0","|#eta| < 0"};
1053  for (Int_t l = 0; l < 3; l++) {
1054  Name = Form("Phi%sdiffR%s%s",NameEta[l],charge[i],GP[k]);
1055  Title = Form("pT_%s - pT_%s for %s and %s",Old.Data(),New.Data(),GPTitle[k],TitleEta[l]);
1056  PhidiffR[l][k][i] = new TH2F(Name,Title,72,-TMath::Pi(),TMath::Pi(),600,-0.15,0.15);
1057  PhidiffR[l][k][i]->SetXTitle(Form("#phi_%s",New.Data()));
1058  YTitle = Form("(pT_%s - pT_%s)/(pT_%s), versus #phi_%s",Old.Data(),New.Data(),New.Data(),New.Data());
1059  yTitle = YTitle; yTitle += title;
1060  PhidiffR[l][k][i]->SetYTitle(yTitle);
1061  PhidiffR[l][k][i]->Sumw2();
1062  PhidiffR[l][k][i]->SetMarkerStyle(20);
1063  PhidiffR[l][k][i]->SetMarkerColor(i+2);
1064 
1065  Name = Form("Phi5%sdiffR%s%s",NameEta[l],charge[i],GP[k]);
1066  Phi5diffR[l][k][i] = new TH2F(Name,Title,72,-TMath::Pi(),TMath::Pi(),600,-0.15,0.15);
1067  Phi5diffR[l][k][i]->SetXTitle(Form("#phi_%s",New.Data()));
1068  YTitle = Form("(pT_%s - pT_%s)/(pT_%s), versus #phi_{%s} for pT > 0.5 GeV/c",Old.Data(),New.Data(),New.Data(),New.Data());
1069  yTitle = YTitle; yTitle += title;
1070  Phi5diffR[l][k][i]->SetYTitle(yTitle);
1071  Phi5diffR[l][k][i]->Sumw2();
1072  Phi5diffR[l][k][i]->SetMarkerStyle(20);
1073  Phi5diffR[l][k][i]->SetMarkerColor(i+2);
1074  }
1075  }
1076  }
1077  // sector 20 and sector 04 pT distributins
1078  TH1D *secpT[2][6];
1079  const Char_t *secNames[4] = {"20","04","17","07"};
1080  const Char_t *pTNames[6] = {Old.Data(), New.Data(),
1081  Form("%sNot%s",Old.Data(),New.Data()), Form("%sNot%s", New.Data(),Old.Data()),
1082  Form("%sAnd%s",Old.Data(),New.Data()), Form("%sAnd%s", New.Data(),Old.Data())};
1083  const Char_t *pTitles[6] = {Old.Data(), New.Data(),
1084  Form("Reconstructed by %s and Not by %s",Old.Data(),New.Data()), Form("reconstructed by %s and Not by %s", New.Data(),Old.Data()),
1085  Form("Reconstructed by both %s And %s",Old.Data(),New.Data()), Form("reconstructed by both %s And %s", New.Data(),Old.Data())};
1086  for (Int_t i = 0; i < 4; i++) { // s = 0 => Sector 20; s = 1 => Sector 04, s = 2 => Sector 17; s = 3 => Sector 07
1087  for (Int_t t = 0; t < 6; t++) {
1088  Name = pTNames[t]; Name += secNames[i];
1089  Title = "Sector "; Title += secNames[i]; Title += " "; Title += pTitles[i];
1090  secpT[i][t] = new TH1D(Name,Title,100,0,10);
1091  }
1092  }
1093  TH2F *etaphi[2];
1094  etaphi[0] = new TH2F(Form("EtaPhi%s",Old.Data()),Form("#phi versus #eta for %s",Old.Data()),
1095  200,-2,2,90,-180,180);
1096  etaphi[1] = new TH2F(Form("EtaPhi%s",New.Data()),Form("#phi versus #eta for %s",New.Data()),
1097  200,-2,2,90,-180,180);
1098  // Loop
1099  Long64_t nentries = fChain->GetEntriesFast();
1100  if (Nentries > 0 && nentries > Nentries) nentries = Nentries;
1101  Int_t nbytes = 0, nb = 0;
1102  Float_t RefMultOld = -1;
1103  for (Long64_t jentry=0; jentry<nentries;jentry++) {
1104  Long64_t ientry = fChain->LoadTree(jentry);
1105  if (ientry < 0) break;
1106  if (jentry%10000 == 0) cout << "Read entry " << jentry << endl;
1107  nb = fChain->GetEntry(jentry); nbytes += nb;
1108  if (data.newFitPtsGl <= 0 && data.oldFitPtsGl <= 0) continue;
1109  if ((data.newFitPtsGl > 0 && TMath::Abs(data.newEtaGl) > 0.5) ||
1110  (data.oldFitPtsGl > 0 && TMath::Abs(data.oldEtaGl) > 0.5)) continue;
1111  fitPtsHist->Fill(data.newFitPtsGl,data.oldFitPtsGl);
1112  Double_t pTdiffGl = -9999;
1113  Double_t pTdiffGlR = -9999;
1114  Double_t pTInvdiffGlR = -9999;
1115  Double_t pTdiffPr = -9999;
1116  Double_t pTdiffPrR = -9999;
1117  Double_t pTInvdiffPrR = -9999;
1118  Int_t oldSec = -1;
1119  if (data.oldFitPtsGl >= effNFP) {
1120  etaphi[0]->Fill(data.oldEtaGl,TMath::RadToDeg()*data.oldPhiGl);
1121  if (data.oldEtaGl < 0) {
1122  if ( -45 < TMath::RadToDeg()*data.oldPhiGl && TMath::RadToDeg()*data.oldPhiGl < -15) oldSec = 0;
1123  if (-135 < TMath::RadToDeg()*data.oldPhiGl && TMath::RadToDeg()*data.oldPhiGl < -105) oldSec = 2;
1124  } else {
1125  if ( -45 < TMath::RadToDeg()*data.oldPhiGl && TMath::RadToDeg()*data.oldPhiGl < -15) oldSec = 1;
1126  if (-135 < TMath::RadToDeg()*data.oldPhiGl && TMath::RadToDeg()*data.oldPhiGl < -105) oldSec = 3;
1127  }
1128  }
1129  Int_t newSec = -1;
1130  if (data.newFitPtsGl >= effNFP) {
1131  etaphi[1]->Fill(data.newEtaGl,TMath::RadToDeg()*data.newPhiGl);
1132  if (data.newEtaGl < 0) {
1133  if ( -45 < TMath::RadToDeg()*data.newPhiGl && TMath::RadToDeg()*data.newPhiGl < -15) newSec = 0;
1134  if (-135 < TMath::RadToDeg()*data.newPhiGl && TMath::RadToDeg()*data.newPhiGl < -105) newSec = 2;
1135  } else {
1136  if ( -45 < TMath::RadToDeg()*data.newPhiGl && TMath::RadToDeg()*data.newPhiGl < -15) newSec = 1;
1137  if (-135 < TMath::RadToDeg()*data.newPhiGl && TMath::RadToDeg()*data.newPhiGl < -105) newSec = 3;
1138  }
1139  }
1140  if (data.newFitPtsGl >= effNFP) {
1141  pTEf[0][0]->Fill(data.newPtGl,refMult);
1142  PhiEf[0][0]->Fill(TMath::RadToDeg()*data.newPhiGl,refMult);
1143  if (newSec >= 0) {
1144  secpT[newSec][1]->Fill(data.newPtGl); // StiCA
1145  if (oldSec < 0)
1146  secpT[newSec][3]->Fill(data.newPtGl); // StiCA but not Sti
1147  else
1148  secpT[newSec][5]->Fill(data.newPtGl); // StiCA and Sti
1149  }
1150  }
1151  if (data.oldFitPtsGl >= effNFP) {
1152  pTEf[0][1]->Fill(data.oldPtGl,refMult);
1153  PhiEf[0][1]->Fill(TMath::RadToDeg()*data.oldPhiGl,refMult);
1154  if (oldSec >= 0) {
1155  secpT[oldSec][0]->Fill(data.oldPtGl); // Sti
1156  if (newSec < 0)
1157  secpT[oldSec][2]->Fill(data.oldPtGl); // Sti but not StiCA
1158  else
1159  secpT[oldSec][4]->Fill(data.oldPtGl); // Sti and StiCA
1160  }
1161  }
1162  Int_t Prim = (Int_t) data.Prim;
1163  if (Prim%10 && data.newFitPtsPr >= effNFP) {pTEf[1][0]->Fill(data.newPtPr,refMult); PhiEf[1][0]->Fill(TMath::RadToDeg()*data.newPhiPr,refMult);}
1164  if (Prim/10 && data.oldFitPtsPr >= effNFP) {pTEf[1][1]->Fill(data.oldPtPr,refMult); PhiEf[1][1]->Fill(TMath::RadToDeg()*data.oldPhiPr,refMult);}
1165  // Matched
1166 #if 0
1167  if (data.newFitPtsGl - data.maxPing >= 3) continue;
1168 #else
1169  if (data.newFitPtsGl > 2*data.maxPing ||
1170  data.oldFitPtsGl > 2*data.maxPing ) continue;
1171 #endif
1172  if (data.newFitPtsGl < minNFP || data.oldFitPtsGl < minNFP) continue;
1173  // if (data.firstHitsDist < 0 || data.firstHitsDist> 1) continue;
1174  fitPtsHistPing->Fill(data.newFitPtsGl,data.oldFitPtsGl);
1175  pTdiffGl = data.oldPtGl - data.newPtGl;
1176  pTdiffGlR = data.oldPtGl/data.newPtGl - 1;
1177  pTInvdiffGlR = data.newPtGl/data.oldPtGl - 1;
1178  if (Prim == 11) {
1179  if (RefMultOld != refMult) {
1180  RefMultOld = refMult;
1181  if (! (data.newPrimX == 0 && data.newPrimY == 0 && data.newPrimZ == 0) &&
1182  ! (data.oldPrimX == 0 && data.oldPrimY == 0 && data.oldPrimZ == 0)) {
1183  DifPvX->Fill(data.newPrimX - data.oldPrimX);
1184  DifPvY->Fill(data.newPrimY - data.oldPrimY);
1185  DifPvZ->Fill(data.newPrimZ - data.oldPrimZ);
1186  }
1187  }
1188  Charge->Fill(data.oldCharge,data.newCharge);
1189  if (data.newFitPtsGl >= 15 && data.oldFitPtsGl >= 15) Charge15->Fill(data.oldCharge,data.newCharge);
1190  Charge15->Fill(data.oldCharge,data.newCharge);
1191  pTdiffPr = data.oldPtPr - data.newPtPr;
1192  pTdiffPrR = data.oldPtPr/data.newPtPr - 1;
1193  pTInvdiffPrR = data.newPtPr/data.oldPtPr - 1;
1194 
1195 
1196  }
1197  Int_t charge = 0;
1198  if ((data.newCharge && data.newCharge < 0) ||
1199  (data.oldCharge && data.oldCharge < 0)) charge = 1;
1200  pTDifNFP[0][charge]->Fill(data.newFitPtsGl,pTdiffGlR);
1201  if (data.newPtGl > 0.5) pTDifNFP5[0][charge]->Fill(data.newFitPtsGl,pTdiffGlR);
1202  if (Prim == 11) {
1203  pTDifNFP[1][charge]->Fill(data.newFitPtsPr,pTdiffPrR);
1204  if (data.newPtPr > 0.5) pTDifNFP5[1][charge]->Fill(data.newFitPtsPr,pTdiffPrR);
1205  }
1206  // Efficiency for more than effNFP fit points
1207  if (data.newFitPtsGl < effNFP || data.oldFitPtsGl < effNFP) continue;
1208  pTdiff[0][charge]->Fill(data.newPtGl,pTdiffGl);
1209  pTdiffR[0][charge]->Fill(data.newPtGl,pTdiffGlR);
1210  pTInvdiffR[0][charge]->Fill(data.newPtGl,pTInvdiffGlR);
1211  PhidiffR[0][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1212  if (data.newEtaGl > 0) PhidiffR[1][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1213  else PhidiffR[2][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1214  if (data.newPtGl > 0.5) {
1215  Phi5diffR[0][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1216  if (data.newEtaGl > 0) Phi5diffR[1][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1217  else Phi5diffR[2][0][charge]->Fill(data.newPhiGl,pTdiffGlR);
1218  }
1219  pTEf[0][2]->Fill(data.newPtGl,refMult); PhiEf[0][2]->Fill(TMath::RadToDeg()*data.newPhiGl,refMult);
1220  pTEf[0][3]->Fill(data.oldPtGl,refMult); PhiEf[0][3]->Fill(TMath::RadToDeg()*data.oldPhiGl,refMult);
1221  if (Prim == 11) {
1222  pTdiff[1][charge]->Fill(data.newPtPr,pTdiffPr);
1223  pTdiffR[1][charge]->Fill(data.newPtPr,pTdiffPrR);
1224  pTInvdiffR[1][charge]->Fill(data.newPtPr,pTInvdiffPrR);
1225  PhidiffR[0][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1226  if (data.newEtaPr > 0) PhidiffR[1][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1227  else PhidiffR[2][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1228  if (data.newPtPr > 0.5) {
1229  Phi5diffR[0][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1230  if (data.newEtaPr > 0) Phi5diffR[1][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1231  else Phi5diffR[2][1][charge]->Fill(data.newPhiPr,pTdiffPrR);
1232  }
1233  pTEf[1][2]->Fill(data.newPtPr,refMult); PhiEf[1][2]->Fill(TMath::RadToDeg()*data.newPhiPr,refMult);
1234  pTEf[1][3]->Fill(data.oldPtPr,refMult); PhiEf[1][3]->Fill(TMath::RadToDeg()*data.oldPhiPr,refMult);
1235  }
1236  }
1237  if (fOut) fOut->Write();
1238  delete fOut;
1239  }
1240  Draw(Out);
1241 }
1242 
mc_data_t
Definition: TbyTPlots.C:60
StarPicoDst::Charge
Charge
Define charge.
Definition: StPicoCommon.h:20
TVolumeView::Draw
virtual void Draw(Option_t *depth="3")
Draw Referenced node with current parameters.
Definition: TVolumeView.cxx:511
Generated by   doxygen 1.8.5