StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
standardPlots.C
1 /*
2  *
3  * $Id: standardPlots.C,v 1.6 2002/08/26 22:17:15 andrewar Exp $
4  * A. Rose, WSU
5  *
6  *
7  * $Log: standardPlots.C,v $
8  * Revision 1.6 2002/08/26 22:17:15 andrewar
9  * Added multi-file input functionality (see doRun() for examples of use). Added
10  * cut string to provided easy equivalence of cut function in on-the-fly hists.
11  *
12  * Revision 1.5 2002/08/19 19:33:20 pruneau
13  * eliminated cout when unnecessary, made helix member of the EventFiller
14  *
15  * Revision 1.4 2002/07/08 14:32:48 pruneau
16  * added efficiency plots
17  *
18  * Revision 1.3 2002/07/02 18:58:52 andrewar
19  * fixed bug with Pt plot
20  *
21  * Revision 1.2 2002/06/26 14:27:54 andrewar
22  * Added cut function and track efficiency hists.
23  *
24  * Revision 1.1 2002/06/12 20:23:35 andrewar
25  * Initial commit. Methods file for standard evaluation plot package.
26  * Few plots are defined; plots will be added as cuts and characteristics
27  * of the tracker are explored.
28  *
29  *
30  *
31  *
32  *
33  */
34 
35 #define standardPlots_cxx
36 #define MAX_TRACKS 7000
37 
38 #include <iostream.h>
39 
40 #include "TSystem.h"
41 #include "TH2.h"
42 #include "TH3.h"
43 #include "TProfile.h"
44 #include "TStyle.h"
45 #include "TCanvas.h"
46 #include "TF1.h"
47 #include "standardPlots.h"
48 
49 void doRun(char* outfile)
50 {
51  standardPlots old(gSystem, "/star/data22/ITTF/EvalData/MCNtuple","*tpt.minimc.root");
52  standardPlots closed(gSystem);
53  standardPlots open(gSystem, "/star/data22/ITTF/EvalData/TempHist/Window20","*evts.minimc.root");
54 
55 
56  old.SetTrackCutNHit(1,55);
57  closed.SetTrackCutNHit(1,55);
58  open.SetTrackCutNHit(1,55);
59  old.SetTrackCutNAHit(1,55);
60  closed.SetTrackCutNAHit(1,55);
61  open.SetTrackCutNAHit(1,55);
62 
63  //do lowest mult first
64  int lowMult=0, highMult=2000;
65  old.SetEventCutMult(lowMult,highMult);
66  closed.SetEventCutMult(lowMult,highMult);
67  open.SetEventCutMult(lowMult,highMult);
68 
69  old.makeTrackEffPlots(100);
70  closed.makeTrackEffPlots(100);
71  open.makeTrackEffPlots(100);
72 
73  TH1D *oldLowMultLowHit=new TH1D("oldLowMultLowHit","",100,0,5);
74  oldLowMultLowHit->Divide(old.primaryTrackPt,old.mcTrackEffPt);
75  TH1D *closedLowMultLowHit=new TH1D("closedLowMultLowHit","",100,0,5);
76  closedLowMultLowHit->Divide(closed.primaryTrackPt,closed.mcTrackEffPt);
77  TH1D *openLowMultLowHit=new TH1D("openLowMultLowHit","",100,0,5);
78  openLowMultLowHit->Divide(open.primaryTrackPt,open.mcTrackEffPt);
79 
80  lowMult=200;
81  highMult=4000;
82  old.SetEventCutMult(lowMult,highMult);
83  closed.SetEventCutMult(lowMult,highMult);
84  open.SetEventCutMult(lowMult,highMult);
85 
86  old.makeTrackEffPlots(100);
87  closed.makeTrackEffPlots(100);
88  open.makeTrackEffPlots(100);
89 
90  TH1D *oldMedMultLowHit=new TH1D("oldMedMultLowHit","",100,0,5);
91  oldMedMultLowHit->Divide(old.primaryTrackPt,old.mcTrackEffPt);
92  TH1D *closedMedMultLowHit=new TH1D("closedMedMultLowHit","",100,0,5);
93  closedMedMultLowHit->Divide(closed.primaryTrackPt,closed.mcTrackEffPt);
94  TH1D *openMedMultLowHit=new TH1D("openMedMultLowHit","",100,0,5);
95  openMedMultLowHit->Divide(open.primaryTrackPt,open.mcTrackEffPt);
96 
97  lowMult=4000;
98  highMult=6000;
99  old.SetEventCutMult(lowMult,highMult);
100  closed.SetEventCutMult(lowMult,highMult);
101  open.SetEventCutMult(lowMult,highMult);
102 
103  old.makeTrackEffPlots(100);
104  closed.makeTrackEffPlots(100);
105  open.makeTrackEffPlots(100);
106 
107  TH1D *oldHiMultLowHit=new TH1D("oldHiMultLowHit","",100,0,5);
108  oldHiMultLowHit->Divide(old.primaryTrackPt,old.mcTrackEffPt);
109  TH1D *closedHiMultLowHit=new TH1D("closedHiMultLowHit","",100,0,5);
110  closedHiMultLowHit->Divide(closed.primaryTrackPt,closed.mcTrackEffPt);
111  TH1D *openHiMultLowHit=new TH1D("openHiMultLowHit","",100,0,5);
112  openHiMultLowHit->Divide(open.primaryTrackPt,open.mcTrackEffPt);
113 
114 
115 
116 
117 
118  TFile *oFile = new TFile(outfile,"RECREATE");
119  oldHiMultLowHit->Write();
120  closedHiMultLowHit->Write();
121  openHiMultLowHit->Write();
122  oldMedMultLowHit->Write();
123  closedMedMultLowHit->Write();
124  openMedMultLowHit->Write();
125  oldLowMultLowHit->Write();
126  closedLowMultLowHit->Write();
127  openLowMultLowHit->Write();
128 
129  oFile->Close();
130 
131 
132 }
133 
134 void standardPlots::Loop()
135 {
136 // In a ROOT session, you can do:
137 // Root > .L standardPlots.C
138 // Root > standardPlots t
139 // Root > t.GetEntry(12); // Fill t data members with entry number 12
140 // Root > t.Show(); // Show values of entry 12
141 // Root > t.Show(16); // Read and show values of entry 16
142 // Root > t.Loop(); // Loop on all entries
143 //
144 
145 // This is the loop skeleton
146 // To read only selected branches, Insert statements like:
147 // METHOD1:
148 // fChain->SetBranchStatus("*",0); // disable all branches
149 // fChain->SetBranchStatus("branchname",1); // activate branchname
150 // METHOD2: replace line
151 // fChain->GetEntry(i); // read all branches
152 //by b_branchname->GetEntry(i); //read only this branch
153  if (fChain == 0) return;
154  nentries = Int_t(fChain->GetEntries());
155  nbytes = 0;
156  nb = 0;
157  for (Int_t jentry=0; jentry<nentries;jentry++) {
158  Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
159  if(!ientry%100) cout <<"Local file entry "<<ientry<<endl;
160  nb = fChain->GetEntry(jentry); nbytes += nb;
161  }
162 }
163 
164 void standardPlots::makeTrackEffPlots(int nEvents)
165 {
166 
167  if (fChain == 0) return;
168 
169  showCuts();
170 
171  //Hists are now initialized in InitHists()
172 
173  nentries = Int_t(fChain->GetEntries());
174  nbytes = 0;
175  nb = 0;
176  for (Int_t jentry=0; jentry<nentries;jentry++) {
177  int ientry=LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
178  nb = fChain->GetEntry(jentry); nbytes += nb;
179 
180  if (!Cut(ientry)) continue;
181 
182  //do some global event characteristic stuff
183  float mPrimaryTrackEff = (float)mMatchedPairs_/(float)mMcTracks_;
184  primaryTrackEffMult->Fill((float)mMcTracks_,mPrimaryTrackEff);
185 
186 
187  double pt;
188  double eta;
189  int nHits;
190  int gId;
191  double q;
192  int iTrack;
193 
194  //make Pt and Eta spectra of MC tracks
195  for(iTrack=0;iTrack<mMcTracks_;iTrack++)
196  {
197 
198 
199  if(!mcTrackCut(ientry,iTrack)) continue; //next track if mcTrackCut doesn't pass
200  pt = mMcTracks_mPtMc[iTrack];
201  eta = mMcTracks_mEtaMc[iTrack];
202  nHits = mMcTracks_mNHitMc[iTrack];
203  gId = mMcTracks_mGeantId[iTrack];
204  q = mMcTracks_mChargeMc[iTrack];
205 
206  mcTrackEffEta->Fill(eta);
207  mcTrackEffPt->Fill(pt);
208  }
209 
210  //make Pt and Eta spectra of matched tracks
211  for(iTrack=0;iTrack<mMatchedPairs_;iTrack++)
212  {
213  if (!trackCut(ientry,iTrack)) continue; //next track if !trackCut
214  pt = mMatchedPairs_mPtMc[iTrack];
215  eta = mMatchedPairs_mEtaMc[iTrack];
216  nHits = mMatchedPairs_mNHitMc[iTrack];
217  gId = mMatchedPairs_mGeantId[iTrack];
218  q = mMatchedPairs_mChargeMc[iTrack];
219 
220  primaryTrackEta->Fill(eta);
221  primaryTrackPt->Fill(pt);
222  }
223 
224 
225  }//end loop over events
226 
227 }
228 
229 
230 void standardPlots::makeMomentumPlots()
231 {
232  cout <<"standardPlots::makeMomentumPlots | Creating Standard Momentum Plots"<<endl;
233  if (fChain == 0) return;
234 
235  nentries = int(fChain->GetEntries());
236  nbytes = 0;
237  nb = 0;
238 
239  TH3D* resolpteta = new TH3D("resolpteta","3d dpt/pt vs pt vs eta",20,-1,1,20,0,2,50,-.5,.5);
240  TH3D* resolpeta = new TH3D("resolpeta","3d dp/p vs p vs eta",20,-1,1,20,0,2,50,-.5,.5);
241  resolpeta->SetXTitle("eta");
242  resolpeta->SetYTitle("p");
243  resolpeta->SetZTitle("#delta p/p");
244 
245  cout <<"standardPlots::makeMomentumPlots +-+Looping over events:"<<endl;
246  cout <<"standardPlots::makeMomentumPlots +-There are "<<nentries<<" events."<<endl;
247 
248  for (int jentry=0; jentry<nentries;jentry++)
249  {
250  int ientry = LoadTree(jentry); //in case of a TChain, ientry
251  //is the entry number in the
252  //current file
253  nb = fChain->GetEntry(jentry);
254  nbytes += nb;
255  if (Cut(ientry) < 0) continue;
256 
257  //now loop over all Matched tracks in event
258  for(int tMatched=0;tMatched<mMatchedPairs_;tMatched++)
259  {
260  if(mMatchedPairs_mFitPts[tMatched]>=24 && mMatchedPairs_mDcaXYGl[tMatched]<1)
261  {
262  float dptp = (mMatchedPairs_mPtMc[tMatched]-mMatchedPairs_mPtPr[tMatched])
263  /mMatchedPairs_mPtPr[tMatched];
264  resolpteta->Fill(mMatchedPairs_mEtaPr[tMatched],mMatchedPairs_mPtMc[tMatched],dptp);
265  dptp = (sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
266  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched])
267  -sqrt(mMatchedPairs_mPtPr[tMatched]*mMatchedPairs_mPtPr[tMatched]
268  +mMatchedPairs_mPzPr[tMatched]*mMatchedPairs_mPzPr[tMatched]))
269  /sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
270  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched]);
271  resolpeta->Fill(mMatchedPairs_mEtaPr[tMatched],
272  sqrt(mMatchedPairs_mPtMc[tMatched]*mMatchedPairs_mPtMc[tMatched]
273  +mMatchedPairs_mPzMc[tMatched]*mMatchedPairs_mPzMc[tMatched]),
274  dptp);
275  }//end cut
276  }//end loop over tracks
277  }//end loop over events
278  cout <<"standardPlots::makeMomentumPlots +-Done."<<endl;
279 
280  //-----------Slice and fit 3D hists-----------
281 
282  cout <<"standardPlots::makeMomentumPlots +-Slicing 3D Hists:"<<endl;
283  return;
284  //Transverse Momentum
285  TString projTit("resolslice");
286  char* extraTitle = "ptYY";
287  TF1* gaus2 = new TF1("gaus2","gaus(0)+gaus(3)",-.25,.25);
288  gaus2->SetParameter(0,500);
289  gaus2->SetParameter(1,0);
290  gaus2->SetParameter(2,.08);
291  gaus2->SetParameter(3,100);
292  gaus2->SetParameter(4,0);
293  gaus2->SetParameter(5,.04);
294  gaus2->SetParName(0,"area");
295  gaus2->SetParName(1,"mean");
296  gaus2->SetParName(2,"width");
297  gaus2->SetParName(3,"area2");
298  gaus2->SetParName(4,"mean2");
299  gaus2->SetParName(5,"width2");
300 
301 
302  TH1D* resolutionPtAtEtaZero = new TH1D("resolutionPtAtEtaZero","RMS (Pt_rec - Pt_mc) vs Pt (eta=0)",20,0,2);
303  TH1D* energylossPtAtEtaZero = new TH1D("energylossPtAtEtaZero","energy loss vs Pt (|eta| < 0.1)",20,0,2);
304  for (int i=2;i<=20;++i)
305  {
306  sprintf(extraTitle,"pt%02i",i);
307  TString currentTit = projTit + extraTitle;
308 
309  TH1D* resolslice = resolpteta->ProjectionZ(currentTit.Data(),10,11,i,i,"e");
310 
311  resolslice->Fit("gaus");
312 // //etaResol->ProjectionZ("etaresolPt",11,11,i,i)->Fit("gaus");
313  //canvas->Modified();
314  //canvas->Update();
315  energylossPtAtEtaZero->SetBinContent(i,resolslice->GetMean());
316  energylossPtAtEtaZero->SetBinError(i,gaus2->GetParError(1));
317  resolutionPtAtEtaZero->SetBinContent(i,gaus2->GetParameter(2));
318  resolutionPtAtEtaZero->SetBinError(i,gaus2->GetParError(2));
319  }
320  energylossPtAtEtaZero->SetMarkerStyle(20);
321  energylossPtAtEtaZero->SetMarkerColor(4);
322  energylossPtAtEtaZero->SetXTitle("pt (GeV/c)");
323  resolutionPtAtEtaZero->SetMarkerStyle(20);
324  resolutionPtAtEtaZero->SetMarkerColor(4);
325  resolutionPtAtEtaZero->SetXTitle("pt (GeV/c)");
326  resolutionPtAtEtaZero->Draw();
327 
328  TH1D* resolutionPAtEtaZero = new TH1D("resolutionPAtEtaZero","RMS (P_rec - P_mc) vs P (eta=0)",20,0,2);
329  TH1D* energylossPAtEtaZero = new TH1D("energylossPAtEtaZero","energy loss vs P (|eta| < 0.1)",20,0,2);
330  for (int i=2;i<=20;++i)
331  {
332  sprintf(extraTitle,"pt%02i",i);
333  TString currentTit = projTit + extraTitle;
334 
335  TH1D* resolslice = resolpeta->ProjectionZ(currentTit.Data(),10,11,i,i,"e");
336 
337  resolslice->Fit("gaus");
338  energylossPAtEtaZero->SetBinContent(i,resolslice->GetMean());
339  energylossPAtEtaZero->SetBinError(i,gaus2->GetParError(1));
340  resolutionPAtEtaZero->SetBinContent(i,gaus2->GetParameter(2));
341  resolutionPAtEtaZero->SetBinError(i,gaus2->GetParError(2));
342  }
343  energylossPAtEtaZero->SetMarkerStyle(20);
344  energylossPAtEtaZero->SetMarkerColor(4);
345  energylossPAtEtaZero->SetXTitle("pt (GeV/c)");
346  resolutionPAtEtaZero->SetMarkerStyle(20);
347  resolutionPAtEtaZero->SetMarkerColor(4);
348  resolutionPAtEtaZero->SetXTitle("pt (GeV/c)");
349  resolutionPAtEtaZero->Draw();
350 
351  //make momentum resolution plot
352 
353  TProfile *myMomResPlot=new TProfile("myMomResPlot","dp/p vs p.",20,0,4);
354  TProfile *mySigmaMomResPlot=new TProfile("mySigmaMomResPlot","RMS(dp/p) vs p.",20,0,4);
355  float totalM=0; float totalMmc=0;
356  nbytes = 0;
357  nb = 0;
358  for (Int_t jentry=0; jentry<nentries;jentry++) {
359  Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
360  nb = fChain->GetEntry(jentry); nbytes += nb;
361  for(int iTrack=0; iTrack<mMatchedPairs_; iTrack++)
362  {
363  if(!trackCut(ientry,iTrack)) continue;
364  totalM=mMatchedPairs_mPtPr[iTrack]*mMatchedPairs_mPtPr[iTrack]
365  +mMatchedPairs_mPzPr[iTrack]*mMatchedPairs_mPzPr[iTrack];
366  totalM=sqrt(totalM);
367 
368  totalMmc=mMatchedPairs_mPtMc[iTrack]*mMatchedPairs_mPtMc[iTrack]
369  +mMatchedPairs_mPzMc[iTrack]*mMatchedPairs_mPzMc[iTrack];
370  totalMmc=sqrt(totalMmc);
371  myMomResPlot->Fill(totalMmc, (totalMmc-totalM)/totalMmc);
372  }
373 
374  }
375  for(int iBin=0; iBin<20;iBin++)
376  {
377  float mom = 4.*((float)(iBin+1))/20.-.1;
378  mySigmaMomResPlot->Fill(mom,myMomResPlot->GetBinError(iBin));
379  }
380  return;
381 }
382 
383 void standardPlots::makeFitPointsPlots()
384 {
385 
386 
387 
388  if (fChain == 0) return;
389 
390  nentries = Int_t(fChain->GetEntries());
391 
392  nbytes = 0;
393  nb = 0;
394  cout <<"Making Plots..."<<endl;
395  for (Int_t jentry=0; jentry<nentries;jentry++) {
396  Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
397  nb = fChain->GetEntry(jentry);
398  nbytes += nb;
399  if (Cut(ientry) < 0) continue;
400 
401  for(Int_t iTrack = 0; iTrack<mMatchedPairs_;iTrack++)
402  {
403  if(trackCut(ientry, iTrack))
404  {
405  fitPointsPlot->Fill((float)mMatchedPairs_mFitPts[iTrack]/(float)mMatchedPairs_mNPossible[iTrack]);
406  fitPointsUsed->Fill((float)mMatchedPairs_mFitPts[iTrack]);
407  fitPointsPossible->Fill(mMatchedPairs_mNPossible[iTrack]);
408  fitPointsPt->Fill(mMatchedPairs_mPtMc[iTrack],
409  (float)mMatchedPairs_mFitPts[iTrack]);
410  fitPointsEta->Fill(mMatchedPairs_mEtaMc[iTrack],
411  (float)mMatchedPairs_mFitPts[iTrack]);
412  }
413  }
414 
415  }
416 
417 }
418 
419 void standardPlots::makeHitEffPlots()
420 {
421  cout <<"standardPlots::hitPlots | Plotting Hit Efficiencies"<<endl;
422 
423  if (fChain == 0) return;
424 
425  nentries = Int_t(fChain->GetEntries());
426  nbytes = 0;
427  nb = 0;
428  for (Int_t jentry=0; jentry<nentries;jentry++)
429  {
430  Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
431  nb = fChain->GetEntry(jentry);
432  nbytes += nb;
433  if (Cut(ientry) < 0) continue;
434 
435  for(Int_t nMatched = 0; nMatched<mMatchedPairs_;nMatched++)
436  {
437 
438  if(!trackCut(ientry, nMatched)) continue;
439 
440  float hitEffRatio=(float)mMatchedPairs_mNCommonHit[nMatched]
441  /(float)mMatchedPairs_mNHitMc[nMatched];
442  hitEffMult->Fill(mNMcTrack, hitEffRatio);
443  //cout <<mNMcTrack <<" "<<hitEffRatio<<endl;
444  hitEffPt->Fill(mMatchedPairs_mPtMc[nMatched],
445  hitEffRatio);
446  hitEffEta->Fill(mMatchedPairs_mEtaMc[nMatched],
447  hitEffRatio);
448 
449  hitN->Fill((float)mMatchedPairs_mNCommonHit[nMatched]);
450 
451  lastHit->Fill(mMatchedPairs_mLastFitPadrow[nMatched]);
452  firstHit->Fill(mMatchedPairs_mFirstFitPadrow[nMatched]);
453  lastHitPt->Fill(mMatchedPairs_mPtMc[nMatched],mMatchedPairs_mLastFitPadrow[nMatched]);
454  firstHitPt->Fill(mMatchedPairs_mPtMc[nMatched],mMatchedPairs_mFirstFitPadrow[nMatched]);
455  lastHitEta->Fill(mMatchedPairs_mEtaMc[nMatched],mMatchedPairs_mLastFitPadrow[nMatched]);
456  firstHitEta->Fill(mMatchedPairs_mEtaMc[nMatched],mMatchedPairs_mFirstFitPadrow[nMatched]);
457  }
458 
459 
460  }//end loop over events
461  hitEffMult->Draw();
462 }
463 
464 void standardPlots::makeGeantIdPlots()
465 {
466  if (fChain == 0) return;
467 
468  nentries = Int_t(fChain->GetEntries());
469 
470  nbytes = 0;
471  nb = 0;
472  cout <<"Making Plots..."<<endl;
473  for (Int_t jentry=0; jentry<nentries;jentry++) {
474  Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
475  nb = fChain->GetEntry(jentry);
476  nbytes += nb;
477  if (Cut(ientry) < 0) continue;
478 
479  for(Int_t iTrack = 0; iTrack<mMatchedPairs_;iTrack++)
480  {
481  if(trackCut(ientry, iTrack))
482  {
483  mGeantIdPlot->Fill(mMatchedPairs_mGeantId[iTrack]);
484  }
485  }
486 
487  }
488 }
Cut
Definition: Cut.h:18
Generated by   doxygen 1.8.5