StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MiniMcPlots.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "Riostream.h"
3 #include <stdio.h>
4 #include "TSystem.h"
5 #include "TMath.h"
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TH3.h"
9 #include "TProfile.h"
10 #include "TStyle.h"
11 #include "TF1.h"
12 #include "TTree.h"
13 #include "TChain.h"
14 #include "TFile.h"
15 #include "TNtuple.h"
16 #include "TCanvas.h"
17 #include "TMinuit.h"
18 #include "TSpectrum.h"
19 #include "TString.h"
20 #include "TLegend.h"
21 #include "TLine.h"
22 #include "TText.h"
23 #include "TROOT.h"
24 #include "TList.h"
25 #include "TPolyMarker.h"
26 #include "StBichsel/Bichsel.h"
27 #include "BetheBloch.h"
28 #include "TDirIter.h"
29 #include "TTreeIter.h"
30 #else
31 class TSystem;
32 class TMath;
33 class TH1;
34 class TH2;
35 class TH3;
36 class TProfile;
37 class TStyle;
38 class TF1;
39 class TTree;
40 class TChain;
41 class TFile;
42 class TNtuple;
43 class TCanvas;
44 class TMinuit;
45 class TSpectrum;
46 class TString;
47 class TLegend;
48 class TLine;
49 class TText;
50 class TROOT;
51 class TList;
52 class TPolyMarker;
53 class Bichsel;
54 class BetheBloch;
55 class TDirIter;
56 class TTreeIter;
57 #endif
58 TChain *chain = 0;
59 const Char_t *Set[2] = {"tpt","sti"};
60 TFile *fIn[2] = {0,0};
61 const Char_t *charge[2] = {"P","N"};
62 const Char_t *scharge[2] = {"+","-"};
63 const Char_t *PG[2] = {"f","Gl"};
64 static Int_t nPng = 0;
65 //________________________________________________________________________________
66 Int_t accept(Float_t pT, Float_t eta, Int_t nMcHit,Int_t GeantId, Int_t mFitPts,
67  Float_t Qual, Int_t mNAssocGl, Int_t mNAssocPr, Char_t validMc=1, Short_t PrimaryMc=1 ) {
68  if (! validMc || ! PrimaryMc) return 0;
69  // if (pT < 0.2 || TMath::Abs(eta) > 0.5 || nMcHit < 40) return 0;
70 #if 0
71  if (TMath::Abs(eta) > 0.5 || nMcHit < 40) return 0;
72  if (mFitPts > 0 && mFitPts < 25) return 0;
73 #else
74  if (TMath::Abs(eta) > 0.5 || nMcHit < 20) return 0;
75  if (mFitPts > 0 && mFitPts < 15) return 0;
76 #endif
77  if (Qual > 0 && Qual < 90.0) return 0;
78  // if (mNAssocGl > 1 || mNAssocPr > 1) return 0;
79  if (GeantId == 8) return -1;
80  if (GeantId == 9) return 1;
81  return 0;
82 }
83 //________________________________________________________________________________
84 void DrawPng(TCanvas *c) {
85  TString pngName("");
86  if (c) {
87  c->Update(); pngName = c->GetName();
88  pngName.ReplaceAll(" ","_");
89  pngName.ReplaceAll("(","_");
90  pngName.ReplaceAll(")","_");
91  pngName.ReplaceAll("{","_");
92  pngName.ReplaceAll("}","_");
93  pngName.ReplaceAll("<","lt");
94  pngName.ReplaceAll(">","gt");
95  pngName.ReplaceAll("tpt_sti_","");
96  pngName.ReplaceAll("tpt_sti","");
97  pngName.ReplaceAll("GeV/c","");
98  pngName.ReplaceAll(".","_");
99  pngName.ReplaceAll("/","_");
100  pngName.ReplaceAll("^","_");
101  pngName.ReplaceAll("__","_");
102  pngName.ReplaceAll("__","_");
103  pngName += ".png";
104  TVirtualX::Instance()->WritePixmap(c->GetCanvasID(),-1,-1,(Char_t *)pngName.Data());
105  nPng++;
106  cout << "Draw #\t" << nPng << "\t" << pngName << endl;
107  }
108 }
109 //________________________________________________________________________________
110 void DrawpTdif(const Char_t *plot = "pTDiff") {
111  TString Plot(plot);
112  TString c1tit = Plot;
113  Double_t yMin[2] = { 99999., 99999.};
114  Double_t yMax[2] = {-99999., -99999.};
115  Double_t xmin = 0;
116  Double_t xmax = 2.5;
117  TFile *f = 0;
118  for (Int_t c = 0; c < 2; c++) {
119  for (Int_t k = 1; k >= 0; k--) {
120  f = fIn[k];
121  if (! f) continue;
122  cout << "File \t" << f->GetName() << endl;
123  TString name = Plot; name += charge[c];
124  f->cd();
125  TH2D * pTDiff = (TH2D*) f->Get(name);
126  if (! pTDiff) continue;
127  pTDiff->FitSlicesY();
128  TString Mu(name); Mu += "_1";
129  TString Sigma(name); Sigma += "_2";
130  TH1D *mu = (TH1D*)gDirectory->Get(Mu);
131  TH1D *sigma = (TH1D*)gDirectory->Get(Sigma);
132  if (! mu || ! sigma) continue;
133  Int_t nx = mu->GetNbinsX();
134  for (Int_t i = 1; i <= nx; i++) {
135  Double_t x = mu->GetBinCenter(i);
136  if (x < xmin || x > xmax) continue;
137  Double_t y = mu->GetBinContent(i);
138  Double_t dy = TMath::Abs(mu->GetBinError(i));
139  if (dy > 0 && TMath::Abs(y) > dy) {
140  if (y - dy < yMin[0]) yMin[0] = y - dy;
141  if (y + dy > yMax[0]) yMax[0] = y + dy;
142  }
143  // cout << "x0\t" << x << "\ty\t" << y << "+-" << dy << "\tmin,max\t" << yMin[0] << "\t" << yMax[0] << endl;
144  y = sigma->GetBinContent(i);
145  dy = TMath::Abs(sigma->GetBinError(i));
146  if (dy > 0 && TMath::Abs(y) > dy) {
147  if (y - dy < yMin[1]) yMin[1] = y - dy;
148  if (y + dy > yMax[1]) yMax[1] = y + dy;
149  }
150  // cout << "x1\t" << x << "\ty\t" << y << "+-" << dy << "\tmin,max\t" << yMin[1] << "\t" << yMax[1] << endl;
151  }
152  }
153  }
154  // cout << "ymin\t" << yMin[0] << "\t" << yMin[1] << endl;
155  // cout << "ymax\t" << yMax[0] << "\t" << yMax[1] << endl;
156  for (Int_t i = 0; i < 2; i++) {
157  if (yMin[i] < 0) yMin[i] *= 1.1; else yMin[i] *= 0.9;
158  if (yMax[i] > 0) yMax[i] *= 1.1; else yMax[i] *= 0.9;
159  }
160  TCanvas *c1 = new TCanvas(c1tit,c1tit,400,400); //c1->SetLeftMargin(0.24);
161  TString c2tit = c1tit;
162  c2tit += "Sigma";
163  TCanvas *c2 = new TCanvas(c2tit,c2tit,400,400); //c2->SetLeftMargin(0.24);
164  TString same1("e");
165  TString same2("e");
166  TString title1, title2;
167  Double_t scale = 1.e3;
168  TLegend *l1 = 0;
169  TLegend *l2 = 0;
170  if (Plot == "pTDiff" || Plot == "pTDifGl") {
171  l1 = new TLegend(.2,.15,.35,.35);
172  l1->SetBorderSize(0);
173  l1->SetFillColor(0);
174  l1->SetTextSize(0.033);
175  l2 = new TLegend(.2,.65,.45,.85);
176  l2->SetBorderSize(0);
177  l2->SetFillColor(0);
178  l2->SetTextSize(0.033);
179  Double_t ymin = -20; Double_t ymax = 5;
180  title1 = "pT diff. #Delta p_{T} (MeV/c) for ";
181  title2 = "sigma of pT diff.#sigma #Delta p_{T} (MeV/c) for ";
182  if (Plot == "pTDifGl") {
183  title1 += "Globals"; title2 += "Globals";
184  } else {
185  title1 += "Primaries"; title2 += "Primaries";
186  }
187  if (scale*yMin[0] > ymin) ymin = scale*yMin[0];
188  if (scale*yMax[0] < ymax) ymax = scale*yMax[0];
189  c1->cd();
190  if (Plot == "pTDifGl") {ymin = -10; ymax = 2.;}
191  else {ymin = -15.; ymax = 5;}
192  TH1F *dummy = c1->DrawFrame(xmin ,ymin, xmax,ymax, title1);
193  dummy->SetName("hframe1");
194  dummy->SetXTitle("p_{TMc} (GeV/c)");
195  dummy->SetYTitle("Delta p_{T} MeV/c");
196  if (scale*yMin[1] > ymin) ymin = scale*yMin[1];
197  if (scale*yMax[1] < ymax) ymax = scale*yMax[1];
198  ymin = 0; ymax = 50;
199  dummy->Draw();
200  c2->cd();
201  if (Plot == "pTDifGl") {ymin = 0; ymax = 50.;}
202  TH1F *dummy2 = c2->DrawFrame(xmin ,ymin, xmax,ymax,title2);
203  dummy2->SetName("hframe2");
204  dummy2->SetXTitle("p_{TMc} (GeV/c)");
205  dummy2->SetYTitle("sigma Delta p_{T} [MeV/c]");
206  dummy2->Draw();
207  } else if (Plot == "pTRel" || Plot == "pTRelGl") {
208  l1 = new TLegend(.6,.20,.8,.45);
209  l1->SetBorderSize(0);
210  l1->SetFillColor(0);
211  l1->SetTextSize(0.033);
212  l2 = new TLegend(.2,.65,.45,.90);
213  l2->SetBorderSize(0);
214  l2->SetFillColor(0);
215  l2->SetTextSize(0.033);
216  scale = 1e2;
217  Double_t ymin = -5; Double_t ymax = 1;
218  title1 = "Relative pT diff.#delta p_{T} (%) for ";
219  title2 = "sigma of relative pT diff. #sigma #delta p_{T} (%) for ";
220  if (Plot == "pTRelGl") {
221  title1 += "Globals"; title2 += "Globals";
222  } else {
223  title1 += "Primaries"; title2 += "Primaries";
224  }
225  if (scale*yMin[0] > ymin) ymin = scale*yMin[0];
226  if (scale*yMax[0] < ymax) ymax = scale*yMax[0];
227  c1->cd();
228  if (Plot == "pTRelGl") {ymin = -5.; ymax = 1;}
229  else {ymin = -2.; ymax = 0.5;}
230  TH1F *dummy = c1->DrawFrame(xmin ,ymin, xmax,ymax, title1);
231  dummy->SetName(c1->GetName());
232  dummy->SetXTitle("p_{TMc} (GeV/c)");
233  dummy->SetYTitle("#delta p_{T} [%]");
234  dummy->Draw();
235  c2->cd();
236  if (scale*yMin[1] > ymin) ymin = scale*yMin[1];
237  if (scale*yMax[1] < ymax) ymax = scale*yMax[1];
238  ymin = 0; ymax = 5;
239  TH1F *dummy2 = c2->DrawFrame(xmin ,ymin, xmax,ymax, title2);
240  dummy2->SetName(c2->GetName());
241  dummy2->SetXTitle("p_{TMc} (GeV/c)");
242  dummy2->SetYTitle("sigma delta p_{T}) [%]");
243  dummy2->Draw();
244  } else return;
245  TString opt("e");
246  for (Int_t c = 0; c < 2; c++) {
247  for (Int_t k = 1; k >= 0; k--) {
248  f = fIn[k];
249  if (! f) continue;
250  cout << "File \t" << f->GetName() << endl;
251  TString name = Plot; name += charge[c];
252  f->cd();
253  TH2D * pTDiff = (TH2D*) f->Get(name);
254  if (! pTDiff) continue;
255  pTDiff->FitSlicesY();
256  TString Mu(name); Mu += "_1";
257  TString Sigma(name); Sigma += "_2";
258  TH1D *mu = (TH1D*)gDirectory->Get(Mu);
259  TH1D *sigma = (TH1D*)gDirectory->Get(Sigma);
260  if (! mu || ! sigma) continue;
261  mu->SetMarkerStyle(21-k);
262  mu->SetMarkerColor(1+c+2*k);
263  sigma->SetMarkerStyle(21-k);
264  sigma->SetMarkerColor(1+c+2*k);
265  TString caption("pi"); caption += scharge[c];
266  if (! Mu.Contains("Gl",TString::kIgnoreCase)) caption += " Pr.";
267  else caption += " Gl.";
268  caption += Set[k];
269  cout << mu->GetName() << "\t" << caption << endl;
270  l1->AddEntry(mu,caption.Data(),"p");
271  mu->Scale(scale);
272  c1->cd();
273  mu->Draw("same");
274  l2->AddEntry(sigma,caption.Data(),"p");
275  sigma->Scale(scale);
276  c2->cd();
277  sigma->Draw("same");
278  }
279  }
280  c1->cd();
281  l1->Draw();
282  c1->Update();
283  DrawPng(c1);
284  c2->cd();
285  l2->Draw();
286  c2->Update();
287  DrawPng(c2);
288 }
289 //________________________________________________________________________________
290 void DrawEff() {
291 #if 1
292  TCanvas *caneff = new TCanvas("Efficiency","Efficiency");
293  TH1D * pTMcP = (TH1D*) gDirectory->Get("pTMcP");
294  TH1D * pTMcN = (TH1D*) gDirectory->Get("pTMcN");
295  TH1D * pTMcA = (TH1D*) gDirectory->Get("pTMcA");
296 
297  TH1D * pTGlP = (TH1D*) gDirectory->Get("pTGlP");
298  TH1D * pTGlN = (TH1D*) gDirectory->Get("pTGlN");
299  TH1D * pTGlA = (TH1D*) gDirectory->Get("pTGlA");
300 
301  TH1D * pTPrP = (TH1D*) gDirectory->Get("pTPrP");
302  TH1D * pTPrN = (TH1D*) gDirectory->Get("pTPrN");
303  TH1D * pTPrA = (TH1D*) gDirectory->Get("pTPrA");
304 
305  TH1D *effN = (TH1D *) pTPrN->Clone(); effN->SetName("effN"); effN->Reset();
306  TH1D *effP = (TH1D *) pTPrP->Clone(); effP->SetName("effP"); effP->Reset();
307  TH1D *effGlN = (TH1D *) pTGlN->Clone(); effGlN->SetName("effGlN");effGlN->Reset();
308  TH1D *effGlP = (TH1D *) pTGlP->Clone(); effGlP->SetName("effGlP");effGlP->Reset();
309  Int_t nx = effN->GetNbinsX();
310  Double_t val, sum, err;
311  for (Int_t l = 1; l <= nx; l++) {
312  val = pTPrN->GetBinContent(l);
313  sum = pTMcN->GetBinContent(l);
314  if (sum < 1.e-7 || val > sum) {val = 0; err = 0;}
315  else { val /= sum; err = TMath::Sqrt(val*(1.-val)/sum);}
316  effN->SetBinContent(l,100.*val);
317  effN->SetBinError(l,100.*err);
318 
319  val = pTPrP->GetBinContent(l);
320  sum = pTMcP->GetBinContent(l);
321  if (sum < 1.e-7 || val > sum) {val = 0; err = 0;}
322  else { val /= sum; err = TMath::Sqrt(val*(1.-val)/sum);}
323  effP->SetBinContent(l,100.*val);
324  effP->SetBinError(l,100.*err);
325 
326  val = pTGlN->GetBinContent(l);
327  sum = pTMcN->GetBinContent(l);
328  if (sum < 1.e-7 || val > sum) {val = 0; err = 0;}
329  else { val /= sum; err = TMath::Sqrt(val*(1.-val)/sum);}
330  effGlN->SetBinContent(l,100.*val);
331  effGlN->SetBinError(l,100.*err);
332 
333  val = pTGlP->GetBinContent(l);
334  sum = pTMcP->GetBinContent(l);
335  if (sum < 1.e-7 || val > sum) {val = 0; err = 0;}
336  else { val /= sum; err = TMath::Sqrt(val*(1.-val)/sum);}
337  effGlP->SetBinContent(l,100.*val);
338  effGlP->SetBinError(l,100.*err);
339  }
340 
341  effP->SetMarkerStyle(28);
342  effP->SetMarkerColor(4);
343  effN->SetMarkerStyle(23);
344  effN->SetMarkerColor(3);
345  effGlP->SetMarkerStyle(28);
346  effGlP->SetMarkerColor(7);
347  effGlN->SetMarkerStyle(23);
348  effGlN->SetMarkerColor(6);
349  TString title("Efficiency (%)."); title += effP->GetTitle();
350  effP->SetTitle(title);
351  effP->SetXTitle("Mc pT [GeV/c]");
352  effP->SetYTitle("Efficiency (%)");
353  TLegend *l4 = new TLegend(.5,.2,0.9,.4);
354  effP->SetStats(0);
355  effP->SetAxisRange(0,2);
356  effP->Draw(); l4->AddEntry(effP,"#pi+ Pr.","p");
357  effN->Draw("same"); l4->AddEntry(effN,"#pi- Pr.","p");
358 // effGlP->Draw(); l4->AddEntry(effGlP,"#pi+ Gl.","p");
359 // effGlN->Draw("same"); l4->AddEntry(effGlN,"#pi- Gl.","p");
360  l4->Draw();
361  DrawPng(caneff);
362 #endif
363 }
364 //________________________________________________________________________________
365 void DrawPrimVx() {
366  Char_t *XYZ[3] = {"X","Y","Z"};
367  for (Int_t i = 0; i < 3; i++) {
368  TH1* h = (TH1 *) gDirectory->Get(Form("DifPv%s",XYZ[i]));
369  if (! h) continue;
370  TCanvas *c = new TCanvas(Form("PrimaryVertex%sdiff",XYZ[i]),Form("PrimaryVertex%sdiff",XYZ[i]));
371  h->Draw();
372  DrawPng(c);
373  }
374 }
375 //________________________________________________________________________________
376 void Draw(const Char_t *fNameTPT, const Char_t *fNameSTI=0) {
377  if (fNameTPT) fIn[0] = new TFile(fNameTPT);
378  if (fNameSTI) fIn[1] = new TFile(fNameSTI);
379  if (! fIn[0] && ! fIn[1] ) return;
380  TString tFile;
381  if (fNameTPT) {tFile = fNameTPT; tFile.ReplaceAll(".root","");}
382  if (fNameSTI) {tFile = fNameSTI; tFile.ReplaceAll(".root","");}
383  DrawpTdif("pTDiff");
384  DrawpTdif("pTDifGl");
385  DrawpTdif("pTRel");
386  DrawpTdif("pTRelGl");
387  DrawEff();
388  DrawPrimVx();
389 }
390 
391 //________________________________________________________________________________
392 void MiniMcPlots(const Char_t *files="*minimc.root") { // root.exe 'MiniMcPlots.C+("*minimc.root")'
393  if (! files) return;
394  TDirIter Dir(files);
395  TTreeIter iter("StMiniMcTree");
396  TFile *f = 0;
397  Int_t NFiles = 0;
398  Char_t *file = 0;
399  Char_t *file1 = 0;
400  while ((file = (Char_t *) Dir.NextFile())) {iter.AddFile(file); NFiles++; file1 = file;}
401  cout << files << "\twith " << NFiles << " files" << "\t" << file1 << endl;
402  if (! file1 ) return;
403  TString pFile(file1);
404  Int_t i1 = pFile.Index("TbyT/");
405  if (i1 >= 0) i1 += 5;
406  else {
407  i1 = pFile.Index("MC/");
408  if (i1 >= 0) i1 += 3;
409  else {
410  i1 = pFile.Index("./");
411  if (i1 >= 0) i1 += 2;
412  else i1 = 0;
413  }
414  }
415  Int_t i2 = pFile.Index(".root");
416  cout << pFile.Data() << "\t" << i1 << "\t/" << i2 << endl;
417  TString tFile(pFile.Data()+i1,i2-i1);
418  tFile.ReplaceAll("/","_"); cout << "Master file\t" << tFile.Data() << endl;
419  TString Out(tFile);
420  Out += ".15_20.Plots.root";
421  TFile *fOut = new TFile(Out.Data(),"recreate");
422 // const Int_t &mEventId = iter("mEventId");
423 // const Int_t &mRunId = iter("mRunId");
424 // const Int_t &mOriginMult = iter("mOriginMult");
425 // const Int_t &mCentralMult = iter("mCentralMult");
426 // const Int_t &mCentrality = iter("mCentrality");
427 // const Int_t &mNUncorrectedNegativePrimaries = iter("mNUncorrectedNegativePrimaries");
428 // const Int_t &mNUncorrectedPrimaries = iter("mNUncorrectedPrimaries");
429 // const Int_t &mNFtpcWUncorrectedPrimaries = iter("mNFtpcWUncorrectedPrimaries");
430 // const Int_t &mNFtpcEUncorrectedPrimaries = iter("mNFtpcEUncorrectedPrimaries");
431 // const Int_t &mMcMult = iter("mMcMult");
432 // const Int_t &mNMcNch = iter("mNMcNch");
433 // const Int_t &mNMcFtpcWNch = iter("mNMcFtpcWNch");
434 // const Int_t &mNMcFtpcENch = iter("mNMcFtpcENch");
435 // const Int_t &mNMcHminus = iter("mNMcHminus");
436 // const Int_t &mNMcGlobal = iter("mNMcGlobal");
437 // const Int_t &mNMcGoodGlobal20 = iter("mNMcGoodGlobal20");
438 // const Int_t &mNRcGlobal = iter("mNRcGlobal");
439 // const Int_t &mNRcGoodGlobal20 = iter("mNRcGoodGlobal20");
440  const Float_t &mVertexX = iter("mVertexX");
441  const Float_t &mVertexY = iter("mVertexY");
442  const Float_t &mVertexZ = iter("mVertexZ");
443  const Float_t &mMcVertexX = iter("mMcVertexX");
444  const Float_t &mMcVertexY = iter("mMcVertexY");
445  const Float_t &mMcVertexZ = iter("mMcVertexZ");
446 // const Float_t &mMagField = iter("mMagField");
447 // const Float_t &mCenterOfMassEnergy = iter("mCenterOfMassEnergy");
448 // const Float_t &mBackgroundRate = iter("mBackgroundRate");
449 // const Short_t &mBeamMassNumberEast = iter("mBeamMassNumberEast");
450 // const Short_t &mBeamMassNumberWest = iter("mBeamMassNumberWest");
451 // const Float_t &mCtb = iter("mCtb");
452 // const Float_t &mZdcE = iter("mZdcE");
453 // const Float_t &mZdcW = iter("mZdcW");
454  const Int_t &mNMcTrack = iter("mNMcTrack");
455  const Int_t &mNMatchedPair = iter("mNMatchedPair");
456 // const Int_t &mNMergedPair = iter("mNMergedPair");
457 // const Int_t &mNSplitPair = iter("mNSplitPair");
458 // const Int_t &mNGhostPair = iter("mNGhostPair");
459 // const Int_t &mNContamPair = iter("mNContamPair");
460 // const Int_t &mNMatGlobPair = iter("mNMatGlobPair");
461  const Int_t &mMcTracks_ = iter("mMcTracks_");
462  //yf const Char_t *&mMcTracks_mIsValid = iter("mMcTracks.mIsValid"); //[mMcTracks_]
463  const Float_t *&mMcTracks_mPtMc = iter("mMcTracks.mPtMc");
464 // const Float_t *&mMcTracks_mPzMc = iter("mMcTracks.mPzMc");
465  const Float_t *&mMcTracks_mEtaMc = iter("mMcTracks.mEtaMc");
466 // const Float_t *&mMcTracks_mPhiMc = iter("mMcTracks.mPhiMc");
467  const Short_t *&mMcTracks_mNHitMc = iter("mMcTracks.mNHitMc");
468 // const Short_t *&mMcTracks_mNSvtHitMc = iter("mMcTracks.mNSvtHitMc");
469 // const Short_t *&mMcTracks_mNSsdHitMc = iter("mMcTracks.mNSsdHitMc");
470 // const Short_t *&mMcTracks_mNFtpcHitMc = iter("mMcTracks.mNFtpcHitMc");
471  const Short_t *&mMcTracks_mGeantId = iter("mMcTracks.mGeantId");
472 // const Short_t *&mMcTracks_mChargeMc = iter("mMcTracks.mChargeMc");
473 // const Float_t *&mMcTracks_mStopR = iter("mMcTracks.mStopR");
474 // const Short_t *&mMcTracks_mKey = iter("mMcTracks.mKey");
475  const Short_t *&mMcTracks_mNAssocGl = iter("mMcTracks.mNAssocGl");
476  const Short_t *&mMcTracks_mNAssocPr = iter("mMcTracks.mNAssocPr");
477  //yf const Short_t *&mMcTracks_mIsPrimary = iter("mMcTracks.mIsPrimary"); //[mMcTracks_]
478  const Int_t &mMatchedPairs_ = iter("mMatchedPairs_");
479 // const Short_t *&mMatchedPairs_mNCommonHit = iter("mMatchedPairs.mNCommonHit");
480 // const Short_t *&mMatchedPairs_mIsBestContam = iter("mMatchedPairs.mIsBestContam");
481 // const Short_t *&mMatchedPairs_mDominatrack = iter("mMatchedPairs.mDominatrack");
482 // const Short_t *&mMatchedPairs_mDominCommonHit = iter("mMatchedPairs.mDominCommonHit");
483  const Float_t *&mMatchedPairs_mAvgQuality = iter("mMatchedPairs.mAvgQuality");
484  //yf const Char_t *&mMatchedPairs_mIsValid = iter("mMatchedPairs.mIsValid"); //[mMatchedPairs_]
485  const Float_t *&mMatchedPairs_mPtMc = iter("mMatchedPairs.mPtMc");
486 // const Float_t *&mMatchedPairs_mPzMc = iter("mMatchedPairs.mPzMc");
487  const Float_t *&mMatchedPairs_mEtaMc = iter("mMatchedPairs.mEtaMc");
488 // const Float_t *&mMatchedPairs_mPhiMc = iter("mMatchedPairs.mPhiMc");
489  const Short_t *&mMatchedPairs_mNHitMc = iter("mMatchedPairs.mNHitMc");
490 // const Short_t *&mMatchedPairs_mNSvtHitMc = iter("mMatchedPairs.mNSvtHitMc");
491 // const Short_t *&mMatchedPairs_mNSsdHitMc = iter("mMatchedPairs.mNSsdHitMc");
492 // const Short_t *&mMatchedPairs_mNFtpcHitMc = iter("mMatchedPairs.mNFtpcHitMc");
493  const Short_t *&mMatchedPairs_mGeantId = iter("mMatchedPairs.mGeantId");
494 // const Short_t *&mMatchedPairs_mChargeMc = iter("mMatchedPairs.mChargeMc");
495 // const Float_t *&mMatchedPairs_mStopR = iter("mMatchedPairs.mStopR");
496 // const Short_t *&mMatchedPairs_mKey = iter("mMatchedPairs.mKey");
497  const Short_t *&mMatchedPairs_mNAssocGl = iter("mMatchedPairs.mNAssocGl");
498  const Short_t *&mMatchedPairs_mNAssocPr = iter("mMatchedPairs.mNAssocPr");
499  const Float_t *&mMatchedPairs_mPtPr = iter("mMatchedPairs.mPtPr");
500  //yf const Char_t *&mMatchedPairs_mIsValidGl = iter("mMatchedPairs.mIsValidGl"); //[mMatchedPairs_]
501  //yf const Short_t *&mMatchedPairs_mIsPrimary = iter("mMatchedPairs.mIsPrimary");
502 // const Float_t *&mMatchedPairs_mPzPr = iter("mMatchedPairs.mPzPr");
503 // const Float_t *&mMatchedPairs_mEtaPr = iter("mMatchedPairs.mEtaPr");
504 // const Float_t *&mMatchedPairs_mPhiPr = iter("mMatchedPairs.mPhiPr");
505 // const Float_t *&mMatchedPairs_mDcaPr = iter("mMatchedPairs.mDcaPr");
506 // const Float_t *&mMatchedPairs_mDcaXYPr = iter("mMatchedPairs.mDcaXYPr");
507 // const Float_t *&mMatchedPairs_mDcaZPr = iter("mMatchedPairs.mDcaZPr");
508 // const Float_t *&mMatchedPairs_mCurvPr = iter("mMatchedPairs.mCurvPr");
509 // const Float_t *&mMatchedPairs_mTanLPr = iter("mMatchedPairs.mTanLPr");
510  // Float_t **&mMatchedPairs_mErrP = iter("mMatchedPairs.mErrP[5]");
511 // const Float_t *&mMatchedPairs_mChi2Pr = iter("mMatchedPairs.mChi2Pr");
512 // const Short_t *&mMatchedPairs_mFlag = iter("mMatchedPairs.mFlag");
513 // const Float_t *&mMatchedPairs_mDedx = iter("mMatchedPairs.mDedx");
514  const Float_t *&mMatchedPairs_mPtGl = iter("mMatchedPairs.mPtGl");
515 // const Float_t *&mMatchedPairs_mPzGl = iter("mMatchedPairs.mPzGl");
516 // const Float_t *&mMatchedPairs_mEtaGl = iter("mMatchedPairs.mEtaGl");
517 // const Float_t *&mMatchedPairs_mPhiGl = iter("mMatchedPairs.mPhiGl");
518 // const Float_t *&mMatchedPairs_mDcaGl = iter("mMatchedPairs.mDcaGl");
519  const Float_t *&mMatchedPairs_mDcaXYGl = iter("mMatchedPairs.mDcaXYGl");
520 // const Float_t *&mMatchedPairs_mDcaZGl = iter("mMatchedPairs.mDcaZGl");
521 // const Float_t *&mMatchedPairs_mCurvGl = iter("mMatchedPairs.mCurvGl");
522 // const Float_t *&mMatchedPairs_mTanLGl = iter("mMatchedPairs.mTanLGl");
523  // Float_t **&mMatchedPairs_mErrG = iter("mMatchedPairs.mErrG[5]");
524 // const Float_t *&mMatchedPairs_mPidPion = iter("mMatchedPairs.mPidPion");
525 // const Float_t *&mMatchedPairs_mPidProton = iter("mMatchedPairs.mPidProton");
526 // const Float_t *&mMatchedPairs_mPidKaon = iter("mMatchedPairs.mPidKaon");
527 // const Float_t *&mMatchedPairs_mPidElectron = iter("mMatchedPairs.mPidElectron");
528 // const Float_t *&mMatchedPairs_mFirstZ = iter("mMatchedPairs.mFirstZ");
529 // const Float_t *&mMatchedPairs_mLastZ = iter("mMatchedPairs.mLastZ");
530 // const Short_t *&mMatchedPairs_mFirstPadrow = iter("mMatchedPairs.mFirstPadrow");
531 // const Short_t *&mMatchedPairs_mLastPadrow = iter("mMatchedPairs.mLastPadrow");
532 // const Short_t *&mMatchedPairs_mFirstFitPadrow = iter("mMatchedPairs.mFirstFitPadrow");
533 // const Short_t *&mMatchedPairs_mLastFitPadrow = iter("mMatchedPairs.mLastFitPadrow");
534 // const Short_t *&mMatchedPairs_mFirstSector = iter("mMatchedPairs.mFirstSector");
535 // const Short_t *&mMatchedPairs_mLastSector = iter("mMatchedPairs.mLastSector");
536  const Short_t *&mMatchedPairs_mFitPts = iter("mMatchedPairs.mFitPts");
537 // const Short_t *&mMatchedPairs_mFitSvt = iter("mMatchedPairs.mFitSvt");
538 // const Short_t *&mMatchedPairs_mFitSsd = iter("mMatchedPairs.mFitSsd");
539 // const Short_t *&mMatchedPairs_mFitFtpc = iter("mMatchedPairs.mFitFtpc");
540 // const Short_t *&mMatchedPairs_mDedxPts = iter("mMatchedPairs.mDedxPts");
541 // const Short_t *&mMatchedPairs_mAllPts = iter("mMatchedPairs.mAllPts");
542 // const Short_t *&mMatchedPairs_mCharge = iter("mMatchedPairs.mCharge");
543 // const Short_t *&mMatchedPairs_mNAssocMc = iter("mMatchedPairs.mNAssocMc");
544 // const Short_t *&mMatchedPairs_mNPossible = iter("mMatchedPairs.mNPossible");
545 //yf const Char_t *&mMatchedPairs_mIsValidPr = iter("mMatchedPairs.mIsValidPr");
546  //#include "miniMc.h"
547  TH1D *DifPvX = new TH1D("DifPvX","Difference in X for Rc - Mc positions",100,-0.25,0.25);
548  TH1D *DifPvY = new TH1D("DifPvY","Difference in Y for Rc - Mc positions",100,-0.25,0.25);
549  TH1D *DifPvZ = new TH1D("DifPvZ","Difference in Z for Rc - Mc positions",100,-0.25,0.25);
550  TString Title("");
551  if (tFile.Contains("tpt",TString::kIgnoreCase)) Title = "TPT :";
552  else Title = "ITTF:";
553  // Title +=" pT>0.2 | #eta |<0.5 && McHit>39 FitPts>24 Q>90";
554  Title +=" #eta |<0.5 && McHit>39 FitPts>24 Q>90";
555  TH2D * pTDiffP = new TH2D("pTDiffP",Title,40,0.,4.,200,-.2,.2);
556  TH2D * pTDiffN = new TH2D("pTDiffN",Title,40,0.,4.,200,-.2,.2);
557  TH2D * pTRelP = new TH2D("pTRelP",Title,40,0.,4.,200,-.2,.2);
558  TH2D * pTRelN = new TH2D("pTRelN",Title,40,0.,4.,200,-.2,.2);
559 
560  TH2D * pTDifGlP = new TH2D("pTDifGlP",Title,40,0.,4.,200,-.2,.2);
561  TH2D * pTDifGlN = new TH2D("pTDifGlN",Title,40,0.,4.,200,-.2,.2);
562  TH2D * pTRelGlP = new TH2D("pTRelGlP",Title,40,0.,4.,200,-.2,.2);
563  TH2D * pTRelGlN = new TH2D("pTRelGlN",Title,40,0.,4.,200,-.2,.2);
564 
565  TH2D * dcaP = new TH2D("dcaP",Title,20,0.,4.,120,-3.,3.);
566  TH2D * dcaN = new TH2D("dcaN",Title,20,0.,4.,120,-3.,3.);
567 
568  TH1D * pTMcP = new TH1D("pTMcP","pT for MC pion+",100,0.,5.);
569  TH1D * pTMcN = new TH1D("pTMcN","pT for MC pion-",100,0.,5.);
570  TH1D * pTMcA = new TH1D("pTMcA","pT for MC all",100,0.,5.);
571 
572  TH1D * pTPrP = new TH1D("pTPrP","pT for Matched pion+",100,0.,5.);
573  TH1D * pTPrN = new TH1D("pTPrN","pT for Matched pion-",100,0.,5.);
574  TH1D * pTPrA = new TH1D("pTPrA","pT for Matched all",100,0.,5.);
575 
576  TH1D * pTGlP = new TH1D("pTGlP","pT for Global pion+",100,0.,5.);
577  TH1D * pTGlN = new TH1D("pTGlN","pT for Global pion-",100,0.,5.);
578  TH1D * pTGlA = new TH1D("pTGlA","pT for Global all",100,0.,5.);
579 
580  TH1D * EtaMcP = new TH1D("EtaMcP","Eta for MC pion+",100,-2.5,2.5);
581  TH1D * EtaMcN = new TH1D("EtaMcN","Eta for MC pion-",100,-2.5,2.5);
582  TH1D * EtaMcA = new TH1D("EtaMcA","Eta for MC all",100,-2.5,2.5);
583 
584  TH1D * EtaPrP = new TH1D("EtaPrP","Eta for Matched pion+",100,-2.5,2.5);
585  TH1D * EtaPrN = new TH1D("EtaPrN","Eta for Matched pion-",100,-2.5,2.5);
586  TH1D * EtaPrA = new TH1D("EtaPrA","Eta for Matched all",100,-2.5,2.5);
587 
588  TH1D * EtaGlP = new TH1D("EtaGlP","Eta for Global pion+",100,-2.5,2.5);
589  TH1D * EtaGlN = new TH1D("EtaGlN","Eta for Global pion-",100,-2.5,2.5);
590  TH1D * EtaGlA = new TH1D("EtaGlA","Eta for Global all",100,-2.5,2.5);
591 
592  // gROOT->LoadMacro("MiniMc.h");
593  Int_t nread = 0;
594  while (iter.Next()) {
595  DifPvX->Fill(mVertexX - mMcVertexX);
596  DifPvY->Fill(mVertexY - mMcVertexY);
597  DifPvZ->Fill(mVertexZ - mMcVertexZ);
598  if (TMath::Abs(mVertexX - mMcVertexX) > 0.25 ||
599  TMath::Abs(mVertexY - mMcVertexY) > 0.25 ||
600  TMath::Abs(mVertexZ - mMcVertexZ) > 0.25) continue;
601  // cout << "Vertex " << mVertexX << "\t" << mVertexY << "\t" << mVertexZ
602  // << "\tMC\t" << mMcVertexX << "\t" << mMcVertexY << "\t" << mMcVertexZ << endl;
603  // cout << "mNMcTrack " << mNMcTrack
604  // << "\tmNMatchedPair " << mNMatchedPair
605  // << "\tmNMergedPair " << mNMergedPair
606  // << "\tmNSplitPair " << mNSplitPair
607  // << "\tmNGhostPair " << mNGhostPair
608  // << "\tmNContamPair " << mNContamPair
609  // << "\tmNMatGlobPair " << mNMatGlobPair << endl;
610  Int_t k;
611  for (k = 0; k < mNMcTrack; k++) {
612  // cout << "McTrack " << k << "\tpT = " << mMcTracks_mPtMc[k] << "\tpZ = " << mMcTracks_mPzMc[k] << endl;
613  Int_t iok = accept(mMcTracks_mPtMc[k],mMcTracks_mEtaMc[k],mMcTracks_mNHitMc[k],mMcTracks_mGeantId[k],-1,-1.,
614  mMcTracks_mNAssocGl[k],mMcTracks_mNAssocPr[k],1,1);
615  // mMcTracks_mIsValid[k],mMcTracks_mIsPrimary[k]);
616  if (! iok) continue;
617  pTMcA->Fill(mMcTracks_mPtMc[k]);
618  EtaMcA->Fill(mMcTracks_mEtaMc[k]);
619  if (iok < 0) {
620  pTMcN->Fill(mMcTracks_mPtMc[k]);
621  EtaMcN->Fill(mMcTracks_mEtaMc[k]);
622  } else {
623  pTMcP->Fill(mMcTracks_mPtMc[k]);
624  EtaMcP->Fill(mMcTracks_mEtaMc[k]);
625  }
626  }
627  for (k = 0; k < mNMatchedPair; k++) {
628  // cout << "mMatchedPairs " << k << "\tpT = " << mMatchedPairs_mPtMc[k] << "\tpZ = " << mMatchedPairs_mPzMc[k] << endl;
629  if (mMatchedPairs_mFitPts[k] < 15) continue;
630  Int_t iok = accept(mMatchedPairs_mPtMc[k],mMatchedPairs_mEtaMc[k],mMatchedPairs_mNHitMc[k],
631  mMatchedPairs_mGeantId[k],mMatchedPairs_mFitPts[k],mMatchedPairs_mAvgQuality[k],
632  mMatchedPairs_mNAssocGl[k],mMatchedPairs_mNAssocPr[k]);
633  //yf mMatchedPairs_mIsValid[k],mMatchedPairs_mIsPrimary[k]);
634  if (! iok) continue;
635  //yf if (! mMatchedPairs_mIsValidGl[k]) continue;
636  Double_t EtaMc = mMatchedPairs_mEtaMc[k];
637  Double_t pTMc = mMatchedPairs_mPtMc[k];
638  Double_t dif = mMatchedPairs_mPtGl[k]-pTMc;
639  Double_t rel = dif/pTMc;
640  EtaGlA->Fill(EtaMc);
641  pTGlA->Fill(pTMc);
642  if (iok < 0) {
643  pTGlN->Fill(pTMc);
644  EtaGlN->Fill(EtaMc);
645  pTDifGlN->Fill(pTMc,dif);
646  pTRelGlN->Fill(pTMc,rel);
647  } else {
648  pTGlP->Fill(pTMc);
649  EtaGlP->Fill(EtaMc);
650  pTDifGlP->Fill(pTMc,dif);
651  pTRelGlP->Fill(pTMc,rel);
652  }
653  //yf if (! mMatchedPairs_mIsValidPr[k]) continue;
654  dif = mMatchedPairs_mPtPr[k]-pTMc;
655  rel = dif/pTMc;
656  EtaPrA->Fill(EtaMc);
657  pTPrA->Fill(pTMc);
658  if (iok < 0) {
659  pTPrN->Fill(pTMc);
660  EtaPrN->Fill(EtaMc);
661  pTDiffN->Fill(pTMc,dif);
662  pTRelN->Fill(pTMc,rel);
663  dcaN->Fill(pTMc,mMatchedPairs_mDcaXYGl[k]);
664  } else {
665  pTPrP->Fill(pTMc);
666  EtaPrP->Fill(EtaMc);
667  pTDiffP->Fill(pTMc,dif);
668  pTRelP->Fill(pTMc,rel);
669  dcaP->Fill(pTMc,mMatchedPairs_mDcaXYGl[k]);
670  }
671  }
672  nread++;
673  // if (nread>100) break;
674  if (nread%1000 == 1) cout << "read " << nread << " events" << endl;
675  }
676  fOut->Write();
677  delete fOut;
678  if (Title.BeginsWith("TPT",TString::kIgnoreCase)) Draw(Out, 0);
679  else Draw(0, Out);
680 }
Generated by   doxygen 1.8.5