StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
standardPlotsClaude.C
1 /*
2  *
3  * $Id: standardPlotsClaude.C,v 1.2 2002/09/05 21:27:18 pruneau Exp $
4  * A. Rose, WSU
5  *
6  *
7  * $Log: standardPlotsClaude.C,v $
8  * Revision 1.2 2002/09/05 21:27:18 pruneau
9  * Fixed problem with StiRootSimpleTrackFilter::makeNewObject
10  *
11  * Revision 1.1 2002/07/11 18:21:31 pruneau
12  * new plots added
13  *
14  * Revision 1.4 2002/07/08 14:32:48 pruneau
15  * added efficiency plots
16  *
17  * Revision 1.3 2002/07/02 18:58:52 andrewar
18  * fixed bug with Pt plot
19  *
20  * Revision 1.2 2002/06/26 14:27:54 andrewar
21  * Added cut function and track efficiency hists.
22  *
23  * Revision 1.1 2002/06/12 20:23:35 andrewar
24  * Initial commit. Methods file for standard evaluation plot package.
25  * Few plots are defined; plots will be added as cuts and characteristics
26  * of the tracker are explored.
27  *
28  *
29  *
30  *
31  *
32  */
33 
34 #define standardPlots_cxx
35 #include "standardPlots.h"
36 #include "TH2.h"
37 #include "TH3.h"
38 #include "TProfile.h"
39 #include "TStyle.h"
40 #include "TCanvas.h"
41 #include "TFile.h"
42 #include "TF1.h"
43 #include "TCanvas.h"
44 #include <iostream.h>
45 
46 #define MAX_TRACKS 7000
47 
48 void run()
49 {
50  TChain * c = new TChain("StMiniMcTree");
51  standardPlots t(c);
52  //t.makeTrackEffPlots("/star/u/pruneau/ittf5/","trackEff.root");
53  t.makeTrackEffPlots("/star/u/pruneau/ittf5/","rcf0183_20_300evts.minimc.root");
54 }
55 
56 void standardPlots::Loop()
57 {
58 // In a ROOT session, you can do:
59 // Root > .L standardPlots.C
60 // Root > standardPlots t
61 // Root > t.GetEntry(12); // Fill t data members with entry number 12
62 // Root > t.Show(); // Show values of entry 12
63 // Root > t.Show(16); // Read and show values of entry 16
64 // Root > t.Loop(); // Loop on all entries
65 //
66 
67 // This is the loop skeleton
68 // To read only selected branches, Insert statements like:
69 // METHOD1:
70 // fChain->SetBranchStatus("*",0); // disable all branches
71 // fChain->SetBranchStatus("branchname",1); // activate branchname
72 // METHOD2: replace line
73 // fChain->GetEntry(i); // read all branches
74 //by b_branchname->GetEntry(i); //read only this branch
75  if (fChain == 0) return;
76  nentries = Int_t(fChain->GetEntries());
77  nbytes = 0;
78  nb = 0;
79  for (Int_t jentry=0; jentry<nentries;jentry++) {
80  ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
81  nb = fChain->GetEntry(jentry); nbytes += nb;
82  }
83 }
84 
85 void standardPlots::makeTrackEffPlots(const TString & path,const TString &fileName)
86 {
87  if (fChain == 0) return;
88 
89 
90  TFile outputFile(path+"/"+fileName,"RECREATE");
91  long j=0;
92 
93  TString name;
94  int ptBins = 16;
95  double ptMin = 0.;
96  double ptMax = 4.;
97  int pzBins = 20;
98  double pzMin = -4.;
99  double pzMax = 4.;
100  int etaBins = 15;
101  double etaMin = -1.5;
102  double etaMax = 1.5;
103  // Book histograms 0:inclusive 1/2/3: centrality bins
104  TString ptString = "pt";
105  TString dptString = "Dpt";
106  TString dpt2String = "Dpt2";
107  TString pzString = "pz";
108  TString dpzString = "Dpz";
109  TString dpz2String = "Dpz2";
110  TString etaString = "eta";
111  TString detaString = "Deta";
112  TString deta2String = "Deta2";
113 
114  double ptBinSizes[] = {0.,0.2,0.4,0.6,0.8,1.0,1.5,2.,3.,4.,5.};
115 
116  TString hitEffString = "hitEff";
117  TString hitEffMultString = "hitEffMult";
118  TString hitEffPtString = "hitEffPt";
119  TString hitEffEtaString = "hitEffEta";
120  name = hitEffMultString+j;
121  hitEffMult = new TProfile(name,name, 25, 0, MAX_TRACKS);
122 
123  for (j=0;j<4;j++)
124  {
125  name = ptString+"Mc"+j;
126  hMcPt[j] = new TH1D(name,name, ptBins, ptMin, ptMax);
127  name = ptString+j;
128  hPt[j] = new TH1D(name,name, ptBins, ptMin, ptMax);
129  name = ptString+"R"+j;
130  hPtR[j] = new TH1D(name,name, ptBins, ptMin, ptMax);
131  name = etaString+"Mc"+j;
132  hMcEta[j] = new TH1D(name,name, etaBins, etaMin, etaMax);
133  name = etaString+j;
134  hEta[j] = new TH1D(name,name, etaBins, etaMin, etaMax);
135  name = etaString+"R"+j;
136  hEtaR[j] = new TH1D(name,name, etaBins, etaMin, etaMax);
137 
138  name = dptString+"_"+j;
139  pDPt[j] = new TProfile(name,name, 10, ptBinSizes, "");
140  name = dpzString+"_"+j;
141  pDPz[j] = new TProfile(name,name, pzBins, pzMin, pzMax);
142  name = detaString+"_"+j;
143  pDEta[j] = new TProfile(name,name, etaBins, etaMin, etaMax);
144  name = dpt2String+"_"+j;
145  pDPt2[j] = new TProfile(name,name, 10, ptBinSizes, "");
146  name = dpz2String+"_"+j;
147  pDPz2[j] = new TProfile(name,name, pzBins, pzMin, pzMax);
148  name = deta2String+"_"+j;
149  pDEta2[j] = new TProfile(name,name, etaBins, etaMin, etaMax);
150 
151  name = dptString+"R_"+j;
152  pDPtR[j] = new TH1D(name,name, 10, ptBinSizes);
153  name = dpzString+"R_"+j;
154  pDPzR[j] = new TH1D(name,name, pzBins, pzMin, pzMax);
155  name = detaString+"R_"+j;
156  pDEtaR[j] = new TH1D(name,name, etaBins, etaMin, etaMax);
157  name = dpt2String+"R_"+j;
158  pDPt2R[j] = new TH1D(name,name, 10, ptBinSizes);
159  name = dpz2String+"R_"+j;
160  pDPz2R[j] = new TH1D(name,name, pzBins, pzMin, pzMax);
161  name = deta2String+"R_"+j;
162  pDEta2R[j] = new TH1D(name,name, etaBins, etaMin, etaMax);
163 
164  name = hitEffString+j;
165  hitEff[j] = new TH1D(name,name, 25, 0.,1.);
166  name = hitEffPtString+j;
167  hitEffPt[j] = new TProfile(name,name, 25, 0.,5.);
168  name = hitEffEtaString+j;
169  hitEffEta[j] = new TProfile(name,name, 25, -1.5,1.5);
170  }
171 
172  nentries = Int_t(fChain->GetEntries());
173  nbytes = 0;
174  nb = 0;
175 
176  double pt=0,mcPt=0,pz=0,mcPz=0,eta=0,mcEta=0,dpt=0,dpz=0,deta=0,hitRatio=0;
177  double nCommon;
178  int nHits=0;
179  int gId=0;
180  double q=0;
181  int iTrack=0;
182  int iC=0;
183 
184  // Event Loop
185  for (Int_t jentry=0; jentry<nentries;jentry++)
186  {
187  //in case of a TChain, ientry is the entry number in the current file
188  ientry = LoadTree(jentry);
189  nb = fChain->GetEntry(jentry); nbytes += nb;
190  //if (!Cut(ientry)) continue;
191 
192  if (mMcMult<2000)
193  iC=1;
194  else if (mMcMult<5000)
195  iC=2;
196  else
197  iC=3;
198 
199  for(Int_t nMatched = 0; nMatched<mMatchedPairs_;nMatched++)
200  {
201  pt = mMatchedPairs_mPtMc[nMatched];
202  eta = mMatchedPairs_mEtaMc[nMatched];
203  nCommon = mMatchedPairs_mNCommonHit[nMatched];
204  nHits = mMatchedPairs_mNHitMc[nMatched];
205  gId = mMatchedPairs_mGeantId[nMatched];
206  q = mMatchedPairs_mChargeMc[nMatched];
207  hitRatio = nCommon/nHits;
208  hitEffPt[0]->Fill(pt,hitRatio);
209  hitEffEta[0]->Fill(eta,hitRatio);
210  if (pt>0.2 && fabs(eta)<=0.5 && gId>3 && fabs(q)>0 && nHits>25)
211  {
212  hitEff[0]->Fill(hitRatio);
213  hitEff[iC]->Fill(hitRatio);
214  hitEffMult->Fill(mNMcTrack, hitRatio);
215  }
216  if (fabs(eta)<=0.5 && gId>3 && fabs(q)>0 && nHits>25)
217  {
218  hitEffPt[0]->Fill(pt,hitRatio);
219  hitEffPt[iC]->Fill(pt,hitRatio);
220  }
221  if (pt>0.2 && gId>3 && fabs(q)>0 && nHits>25)
222  {
223  hitEffEta[0]->Fill(eta,hitRatio);
224  hitEffEta[iC]->Fill(eta,hitRatio);
225  }
226  }
227 
228  //make Pt and Eta spectra of MC tracks
229  for(iTrack=0;iTrack<mMcTracks_;iTrack++)
230  {
231  //if(!mcTrackCut(ientry,iTrack)) continue; //next track if mcTrackCut doesn't pass
232  pt = mMcTracks_mPtMc[iTrack];
233  eta = mMcTracks_mEtaMc[iTrack];
234  nHits = mMcTracks_mNHitMc[iTrack];
235  gId = mMcTracks_mGeantId[iTrack];
236  q = mMcTracks_mChargeMc[iTrack];
237 
238  if (fabs(eta)<=0.5 && gId>3 && fabs(q)>0 && nHits>25)
239  {
240  hMcPt[0]->Fill(pt);
241  hMcPt[iC]->Fill(pt);
242  }
243  if (pt>0.2 && gId>3 && fabs(q)>0 && nHits>25)
244  {
245  hMcEta[0]->Fill(eta);
246  hMcEta[iC]->Fill(eta);
247  }
248  }
249 
250  //make Pt and Eta spectra of matched tracks
251  for(iTrack=0;iTrack<mMatchedPairs_;iTrack++)
252  {
253  //if (!trackCut(ientry,iTrack)) continue; //next track if trackCut
254 
255  mcPt = mMatchedPairs_mPtMc[iTrack];
256  mcPz = mMatchedPairs_mPzMc[iTrack];
257  mcEta = mMatchedPairs_mEtaMc[iTrack];
258  pt = mMatchedPairs_mPtPr[iTrack];
259  pz = mMatchedPairs_mPzPr[iTrack];
260  eta = mMatchedPairs_mEtaPr[iTrack];
261  nHits = mMatchedPairs_mNHitMc[iTrack];
262  gId = mMatchedPairs_mGeantId[iTrack];
263  q = mMatchedPairs_mChargeMc[iTrack];
264 
265  if (fabs(eta)<=0.5 && gId>3 && fabs(q)>0 && nHits>25)
266  {
267  hPt[0]->Fill(mcPt);
268  hPt[iC]->Fill(mcPt);
269  dpt = pt-mcPt;
270  pDPt[0]->Fill(mcPt,dpt);
271  pDPt[iC]->Fill(mcPt,dpt);
272  pDPt2[0]->Fill(mcPt,dpt*dpt);
273  pDPt2[iC]->Fill(mcPt,dpt*dpt);
274  }
275  if (pt>0.2 && gId>3 && fabs(q)>0 && nHits>25)
276  {
277  hEta[0]->Fill(mcEta);
278  hEta[iC]->Fill(mcEta);
279  dpz = pz-mcPz;
280  pDPz[0]->Fill(mcPz,dpz);
281  pDPz[iC]->Fill(mcPz,dpz);
282  pDPz2[0]->Fill(mcPz,dpz*dpz);
283  pDPz2[iC]->Fill(mcPz,dpz*dpz);
284  deta = eta-mcEta;
285  pDEta[0]->Fill(mcEta,deta);
286  pDEta[iC]->Fill(mcEta,deta);
287  pDEta2[0]->Fill(mcEta,deta*deta);
288  pDEta2[iC]->Fill(mcEta,deta*deta);
289  }
290  }
291  }
292  // Save
293  hitEffMult->Write();
294  for (j=0;j<4;j++)
295  {
296 
297  hitEff[j]->Write();
298  hitEffPt[j]->Write();
299  hitEffEta[j]->Write();
300 
301  hPtR[j]->Divide(hPt[j],hMcPt[j]);
302  hEtaR[j]->Divide(hEta[j],hMcEta[j]);
303 
304  hPt[j]->Write();
305  hMcPt[j]->Write();
306  hPtR[j]->Write();
307  hEta[j]->Write();
308  hMcEta[j]->Write();
309  hEtaR[j]->Write();
310 
311  pDPt[j]->Write();
312  pDPz[j]->Write();
313  pDEta[j]->Write();
314  pDPt2[j]->Write();
315  pDPz2[j]->Write();
316  pDEta2[j]->Write();
317 
318  calculateRel(pDPt[j],pDPt2[j],pDPtR[j], pDPt2R[j]);
319  calculateRel(pDPz[j],pDPz2[j],pDPzR[j], pDPz2R[j]);
320  calculateRel(pDEta[j],pDEta2[j],pDEtaR[j], pDEta2R[j]);
321 
322  pDPtR[j]->Write();
323  pDPzR[j]->Write();
324  pDEtaR[j]->Write();
325  pDPt2R[j]->Write();
326  pDPz2R[j]->Write();
327  pDEta2R[j]->Write();
328 
329  }
330  print(path,fileName);
331  outputFile.Close();
332 }
333 
334 void standardPlots::calculateRel(const TProfile *p1, const TProfile *p2,
335  TH1D *h1, TH1D *h2)
336 {
337  double x=0, mean=0, eMean=0, relMean=0, eRelMean=0, var=0,eVar=0,s=0,eS=0;
338  int n = p1->GetNbinsX();
339  for (int k=1;k<=n;k++)
340  {
341  mean = p1->GetBinContent(k);
342  eMean = p1->GetBinError(k);
343  x = fabs(p1->GetBinCenter(k));
344  if (x==0) continue;
345  relMean = mean/x; // fractional bias
346  eRelMean = eMean/x;
347  h1->SetBinContent(k,relMean);
348  h1->SetBinError(k,eRelMean);
349 
350  var = p2->GetBinContent(k);
351  eVar = p2->GetBinError(k);
352  s = sqrt(var)/x;
353  if (var==0.)
354  eS=0;
355  else
356  eS = s*eVar/var;
357  h2->SetBinContent(k,s);
358  h2->SetBinError(k,eS);
359  }
360 }
361 
362 void standardPlots::print(const TString & path,const TString &fileName)
363 {
364  TCanvas * c = new TCanvas();
365  hMcPt[0]->Draw();
366  hPt[0]->Draw("SAME");
367  c->Print(path+"nTrackVsPt.gif");
368  cout<<"set the options"<<endl;
369  int done;
370  cin >> done;
371  for (int j=0;j<4;j++)
372  {
373  hPt[j]->SetLineColor(j+1);
374  hMcPt[j]->SetLineColor(j+1);
375  hPtR[j]->SetLineColor(j+1);
376  hEta[j]->SetLineColor(j+1);
377  hMcEta[j]->SetLineColor(j+1);
378  hEtaR[j]->SetLineColor(j+1);
379  pDPt[j]->SetLineColor(1+j);
380  pDPt2[j]->SetLineColor(1+j);
381  pDPtR[j]->SetLineColor(1+j);
382  pDPt2R[j]->SetLineColor(1+j);
383  pDPz[j]->SetLineColor(1+j);
384  pDPz2[j]->SetLineColor(1+j);
385  pDPzR[j]->SetLineColor(1+j);
386  pDPz2R[j]->SetLineColor(1+j);
387  pDEta[j]->SetLineColor(1+j);
388  pDEta2[j]->SetLineColor(1+j);
389  pDEtaR[j]->SetLineColor(1+j);
390  pDEta2R[j]->SetLineColor(1+j);
391 
392  hPtR[j]->SetMinimum(0.);
393  hPtR[j]->SetMaximum(1.1);
394  hEtaR[j]->SetMinimum(0.);
395  hEtaR[j]->SetMaximum(1.1);
396  }
397 
398  hPtR[0]->Draw();
399  hPtR[1]->Draw("SAME");
400  hPtR[2]->Draw("SAME");
401  hPtR[3]->Draw("SAME");
402  c->Print(path+"trackEffVsPt.gif");
403 
404  hMcEta[0]->Draw();
405  hEta[0]->Draw("SAME");
406  c->Print(path+"nTrackVsEta.gif");
407 
408  hEtaR[0]->Draw();
409  hEtaR[1]->Draw("SAME");
410  hEtaR[2]->Draw("SAME");
411  hEtaR[3]->Draw("SAME");
412  c->Print(path+"trackEffVsEta.gif");
413 
414 
415  pDPt[0]->SetMinimum(-1.);
416  pDPt[0]->SetMaximum(1.);
417  pDPt[0]->Draw();
418  pDPt[1]->Draw("SAME");
419  pDPt[2]->Draw("SAME");
420  pDPt[3]->Draw("SAME");
421  c->Print(path+"dPtVsPt.gif");
422 
423  pDPt2[0]->SetMinimum(-1.);
424  pDPt2[0]->SetMaximum(1.);
425  pDPt2[0]->Draw();
426  pDPt2[1]->Draw("SAME");
427  pDPt2[2]->Draw("SAME");
428  pDPt2[3]->Draw("SAME");
429  c->Print(path+"dPt2VsPt.gif");
430 
431  pDPtR[0]->SetMinimum(-0.3);
432  pDPtR[0]->SetMaximum(0.3);
433  pDPtR[0]->Draw();
434  pDPtR[1]->Draw("SAME");
435  pDPtR[2]->Draw("SAME");
436  pDPtR[3]->Draw("SAME");
437  c->Print(path+"dPtRelVsPt.gif");
438 
439  pDPt2R[3]->SetMinimum(0.);
440  pDPt2R[3]->SetMaximum(0.3);
441  //pDPt2R[0]->Draw();
442  //pDPt2R[1]->Draw("SAME");
443  //pDPt2R[2]->Draw("SAME");
444  pDPt2R[3]->Draw();
445  c->Print(path+"dPt2RelVsPt.gif");
446 
447  pDPz[0]->SetMinimum(-0.2);
448  pDPz[0]->SetMaximum(0.2);
449  pDPz[0]->Draw();
450  pDPz[1]->Draw("SAME");
451  pDPz[2]->Draw("SAME");
452  pDPz[3]->Draw("SAME");
453  c->Print(path+"dPzVsPz.gif");
454 
455  pDPz2[0]->SetMinimum(0.);
456  pDPz2[0]->SetMaximum(0.5);
457  pDPz2[0]->Draw();
458  //pDPz2[1]->Draw("SAME");
459  //pDPz2[2]->Draw("SAME");
460  //pDPz2[3]->Draw("SAME");
461  c->Print(path+"dPz2VsPz.gif");
462 
463  pDPzR[0]->SetMinimum(-0.2);
464  pDPzR[0]->SetMaximum(0.2);
465  pDPzR[0]->Draw();
466  pDPzR[1]->Draw("SAME");
467  pDPzR[2]->Draw("SAME");
468  pDPzR[3]->Draw("SAME");
469  c->Print(path+"dPzRelVsPz.gif");
470 
471  pDPz2R[0]->SetMinimum(0.);
472  pDPz2R[0]->SetMaximum(0.5);
473  pDPz2R[0]->Draw();
474  //pDPz2R[1]->Draw("SAME");
475  //pDPz2R[2]->Draw("SAME");
476  //pDPz2R[3]->Draw("SAME");
477  c->Print(path+"dPz2RelVsPz.gif");
478 
479  pDEta[0]->SetMinimum(-0.5);
480  pDEta[0]->SetMaximum(0.5);
481  pDEta[0]->Draw();
482  pDEta[1]->Draw("SAME");
483  pDEta[2]->Draw("SAME");
484  pDEta[3]->Draw("SAME");
485  c->Print(path+"dEtaVsEta.gif");
486 
487  pDEta2[0]->SetMinimum(0.);
488  pDEta2[0]->SetMaximum(0.5);
489  pDEta2[0]->Draw();
490  pDEta2[1]->Draw("SAME");
491  pDEta2[2]->Draw("SAME");
492  pDEta2[3]->Draw("SAME");
493  c->Print(path+"dEta2VsEta.gif");
494 
495  pDEtaR[0]->SetMinimum(-0.5);
496  pDEtaR[0]->SetMaximum(0.5);
497  pDEtaR[0]->Draw();
498  pDEtaR[1]->Draw("SAME");
499  pDEtaR[2]->Draw("SAME");
500  pDEtaR[3]->Draw("SAME");
501  c->Print(path+"dEtaRelVsEta.gif");
502 
503  pDEta2R[0]->SetMinimum(0.);
504  pDEta2R[0]->SetMaximum(0.5);
505  pDEta2R[0]->Draw();
506  //pDEta2R[1]->Draw("SAME");
507  //pDEta2R[2]->Draw("SAME");
508  //pDEta2R[3]->Draw("SAME");
509  c->Print(path+"dEta2RelVsEta.gif");
510 }
511 
512 void standardPlots::makeMomentumPlots()
513 {
514  cout <<"standardPlots::makeMomentumPlots | Creating Standard Momentum Plots"<<endl;
515  if (fChain == 0) return;
516 
517  nentries = int(fChain->GetEntries());
518  nbytes = 0;
519  nb = 0;
520 
521  TH3D* resolpteta = new TH3D("resolpteta","3d dpt/pt vs pt vs eta",20,-1,1,20,0,2,50,-.5,.5);
522  TH3D* resolpeta = new TH3D("resolpeta","3d dp/p vs p vs eta",20,-1,1,20,0,2,50,-.5,.5);
523  resolpeta->SetXTitle("eta");
524  resolpeta->SetYTitle("p");
525  resolpeta->SetZTitle("#delta p/p");
526 
527  cout <<"standardPlots::makeMomentumPlots +-+Looping over events:"<<endl;
528  cout <<"standardPlots::makeMomentumPlots +-There are "<<nentries<<" events."<<endl;
529 
530  for (int jentry=0; jentry<nentries;jentry++)
531  {
532  ientry = LoadTree(jentry); //in case of a TChain, ientry
533  //is the entry number in the
534  //current file
535  nb = fChain->GetEntry(jentry);
536  nbytes += nb;
537  // if (Cut(ientry) < 0) continue;
538 
539  //now loop over all Matched tracks in event
540  for(int tMatched=0;tMatched<mMatchedPairs_;tMatched++)
541  {
542  if(mMatchedPairs_mFitPts[tMatched]>=24 && mMatchedPairs_mDcaXYGl[tMatched]<1)
543  {
544  float dptp = (mMatchedPairs_mPtMc[tMatched]-mMatchedPairs_mPtPr[tMatched])
545  /mMatchedPairs_mPtPr[tMatched];
546  resolpteta->Fill(mMatchedPairs_mEtaPr[tMatched],mMatchedPairs_mPtMc[tMatched],dptp);
547  dptp = (sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
548  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched])
549  -sqrt(mMatchedPairs_mPtPr[tMatched]*mMatchedPairs_mPtPr[tMatched]
550  +mMatchedPairs_mPzPr[tMatched]*mMatchedPairs_mPzPr[tMatched]))
551  /sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
552  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched]);
553  resolpeta->Fill(mMatchedPairs_mEtaPr[tMatched],
554  sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
555  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched]),
556  dptp);
557  }//end cut
558  }//end loop over tracks
559  }//end loop over events
560  cout <<"standardPlots::makeMomentumPlots +-Done."<<endl;
561 
562  //-----------Slice and fit 3D hists-----------
563 
564  cout <<"standardPlots::makeMomentumPlots +-Slicing 3D Hists:"<<endl;
565 
566  //Transverse Momentum
567  TString projTit("resolslice");
568  char* extraTitle = "ptYY";
569  TF1* gaus2 = new TF1("gaus2","gaus(0)+gaus(3)",-.25,.25);
570  gaus2->SetParameter(0,500);
571  gaus2->SetParameter(1,0);
572  gaus2->SetParameter(2,.08);
573  gaus2->SetParameter(3,100);
574  gaus2->SetParameter(4,0);
575  gaus2->SetParameter(5,.04);
576  gaus2->SetParName(0,"area");
577  gaus2->SetParName(1,"mean");
578  gaus2->SetParName(2,"width");
579  gaus2->SetParName(3,"area2");
580  gaus2->SetParName(4,"mean2");
581  gaus2->SetParName(5,"width2");
582 
583 
584  TH1D* resolutionPtAtEtaZero = new TH1D("resolutionPtAtEtaZero","RMS (Pt_rec - Pt_mc) vs Pt (eta=0)",20,0,2);
585  TH1D* energylossPtAtEtaZero = new TH1D("energylossPtAtEtaZero","energy loss vs Pt (|eta| < 0.1)",20,0,2);
586  for (int i=2;i<=20;++i)
587  {
588  sprintf(extraTitle,"pt%02i",i);
589  TString currentTit = projTit + extraTitle;
590 
591  TH1D* resolslice = resolpteta->ProjectionZ(currentTit.Data(),10,11,i,i,"e");
592 
593  resolslice->Fit("gaus");
594 // //etaResol->ProjectionZ("etaresolPt",11,11,i,i)->Fit("gaus");
595  //canvas->Modified();
596  //canvas->Update();
597  energylossPtAtEtaZero->SetBinContent(i,resolslice->GetMean());
598  energylossPtAtEtaZero->SetBinError(i,gaus2->GetParError(1));
599  resolutionPtAtEtaZero->SetBinContent(i,gaus2->GetParameter(2));
600  resolutionPtAtEtaZero->SetBinError(i,gaus2->GetParError(2));
601  }
602  energylossPtAtEtaZero->SetMarkerStyle(20);
603  energylossPtAtEtaZero->SetMarkerColor(4);
604  energylossPtAtEtaZero->SetXTitle("pt (GeV/c)");
605  resolutionPtAtEtaZero->SetMarkerStyle(20);
606  resolutionPtAtEtaZero->SetMarkerColor(4);
607  resolutionPtAtEtaZero->SetXTitle("pt (GeV/c)");
608  resolutionPtAtEtaZero->Draw();
609 
610  TH1D* resolutionPAtEtaZero = new TH1D("resolutionPAtEtaZero","RMS (P_rec - P_mc) vs P (eta=0)",20,0,2);
611  TH1D* energylossPAtEtaZero = new TH1D("energylossPAtEtaZero","energy loss vs P (|eta| < 0.1)",20,0,2);
612  for (int i=2;i<=20;++i)
613  {
614  sprintf(extraTitle,"pt%02i",i);
615  TString currentTit = projTit + extraTitle;
616 
617  TH1D* resolslice = resolpeta->ProjectionZ(currentTit.Data(),10,11,i,i,"e");
618 
619  resolslice->Fit("gaus");
620  energylossPAtEtaZero->SetBinContent(i,resolslice->GetMean());
621  energylossPAtEtaZero->SetBinError(i,gaus2->GetParError(1));
622  resolutionPAtEtaZero->SetBinContent(i,gaus2->GetParameter(2));
623  resolutionPAtEtaZero->SetBinError(i,gaus2->GetParError(2));
624  }
625  energylossPAtEtaZero->SetMarkerStyle(20);
626  energylossPAtEtaZero->SetMarkerColor(4);
627  energylossPAtEtaZero->SetXTitle("pt (GeV/c)");
628  resolutionPAtEtaZero->SetMarkerStyle(20);
629  resolutionPAtEtaZero->SetMarkerColor(4);
630  resolutionPAtEtaZero->SetXTitle("pt (GeV/c)");
631  resolutionPAtEtaZero->Draw();
632 
633  //make momentum resolution plot
634 
635  TProfile *myMomResPlot=new TProfile("myMomResPlot","dp/p vs p.",20,0,4);
636  TProfile *mySigmaMomResPlot=new TProfile("mySigmaMomResPlot","RMS(dp/p) vs p.",20,0,4);
637  float totalM=0; float totalMmc=0;
638  nbytes = 0;
639  nb = 0;
640  for (Int_t jentry=0; jentry<nentries;jentry++) {
641  ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
642  nb = fChain->GetEntry(jentry); nbytes += nb;
643  for(int iTrack=0; iTrack<mMatchedPairs_; iTrack++)
644  {
645  if(!trackCut(ientry,iTrack)) continue;
646  totalM=mMatchedPairs_mPtPr[iTrack]*mMatchedPairs_mPtPr[iTrack]
647  +mMatchedPairs_mPzPr[iTrack]*mMatchedPairs_mPzPr[iTrack];
648  totalM=sqrt(totalM);
649 
650  totalMmc=mMatchedPairs_mPtMc[iTrack]*mMatchedPairs_mPtMc[iTrack]
651  +mMatchedPairs_mPzMc[iTrack]*mMatchedPairs_mPzMc[iTrack];
652  totalMmc=sqrt(totalMmc);
653  myMomResPlot->Fill(totalMmc, (totalMmc-totalM)/totalMmc);
654  }
655 
656  }
657  for(int iBin=0; iBin<20;iBin++)
658  {
659  float mom = 4.*((float)(iBin+1))/20.-.1;
660  mySigmaMomResPlot->Fill(mom,myMomResPlot->GetBinError(iBin));
661  }
662  return;
663 }
664 
665 void standardPlots::makeFitPointsPlots()
666 {
667 
668  TH1D* fitPointsPlot=new TH1D("fitPointsPlot","Number of Fit Points / Number Possible", 50, 0,1);
669  TH1D* fitPointsUsed=new TH1D("fitPointsUsed","Number of Fit Points Used", 50, 0, 50);
670  TH1D* fitPointsPossible=new TH1D("fitPointsPossible","Number of Fit Points Possible", 50, 0,50);
671 
672 
673 
674  if (fChain == 0) return;
675 
676  nentries = Int_t(fChain->GetEntries());
677 
678  nbytes = 0;
679  nb = 0;
680  for (Int_t jentry=0; jentry<nentries;jentry++) {
681  ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
682  nb = fChain->GetEntry(jentry);
683  nbytes += nb;
684  // if (Cut(ientry) < 0) continue;
685 
686  for(Int_t nMatched = 0; nMatched<mMatchedPairs_;nMatched++)
687  {
688  if(mMatchedPairs_mFitPts[nMatched]>=24 && mMatchedPairs_mDcaXYGl[nMatched]<1
689  && mMatchedPairs_mPtMc[nMatched] >.4)
690  {
691  fitPointsPlot->Fill(mMatchedPairs_mFitPts[nMatched]/mMatchedPairs_mNPossible[nMatched]);
692  fitPointsUsed->Fill(mMatchedPairs_mFitPts[nMatched]);
693  fitPointsPossible->Fill(mMatchedPairs_mNPossible[nMatched]);
694  }
695  }
696 
697  }
698  fitPointsPlot->Draw();
699 
700 
701  fitPointsUsed->SetMarkerColor(3);
702  fitPointsUsed->SetLineColor(3);
703  fitPointsPossible->SetMarkerColor(4);
704  fitPointsPossible->SetLineColor(4);
705  fitPointsPossible->Draw();
706  fitPointsUsed->Draw("same");
707 }
708 
Generated by   doxygen 1.8.5