Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
observHister.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file observHister.C
1 #include <iostream>
2 using namespace std;
3 
4 #include "TFile.h"
5 #include "Scalar.h"
6 #include "TTree.h"
7 #include "TChain.h"
8 #include "TLegend.h"
9 #include "math.h"
10 #include "TH1.h"
11 #include "TH2.h"
12 #include "TEfficiency.h"
13 #include "TLine.h"
14 #include "TGraphAsymmErrors.h"
15 
16 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
17  TChain *all = new TChain(treename.c_str());
18  string temp;
19  for (int i = 0; i < filecount; ++i)
20  {
21 
22  ostringstream s;
23  s<<i;
24  temp = name+string(s.str())+extension;
25  all->Add(temp.c_str());
26  }
27  return all;
28 }
29 
30 void matchPlot(TChain* ttree,TFile* out_file){
31  int b_match;
32  int b_unmatch;
33  ttree->SetBranchAddress("nMatched", & b_match );
34  ttree->SetBranchAddress("nUnmatched", &b_unmatch );
35  unsigned totalmatch=0;
36  unsigned totaltrack=0;
37  for (int event = 0; event < ttree->GetEntries(); ++event)
38  {
39  ttree->GetEvent(event);
40 
41  totalmatch+=b_match;
42  totaltrack+=b_match;
43  totaltrack+=b_unmatch;
44  }
45  cout<<"Match "<<totalmatch<<" tracks\n";
46  cout<<"Matched "<<(float)totalmatch/totaltrack*100.<<"\% of tracks"<<endl;
47 }
48 
49 void makeVtxR(TChain* ttree,TFile* out_file){
50  float vtxr;
51  float tvtxr;
52  ttree->SetBranchAddress("vtx_radius",&vtxr);
53  ttree->SetBranchAddress("tvtx_radius",&tvtxr);
54 
55  std::vector<TH1F*> plots;
56  plots.push_back(new TH1F("vtx_reco","",40,0,30));
57  plots.push_back(new TH1F("vtx_truth","",40,0,30));
58 
59  plots[0]->Sumw2();
60  plots[1]->Sumw2();
61 
62  double calc=0;
63  for (int event = 0; event < ttree->GetEntries(); ++event)
64  {
65  ttree->GetEvent(event);
66  plots[0]->Fill(vtxr);
67  plots[1]->Fill(tvtxr);
68  calc+=TMath::Abs(vtxr-tvtxr);
69  }
70  calc/=ttree->GetEntries();
71  for (int i = 0; i < 2; ++i)
72  {
73  plots[i]->Scale(1./ttree->GetEntries(),"width");
74  }
75  out_file->Write();
76  std::cout<<"mean deviation="<<calc<<std::endl;
77 }
78 
79 void makepTEff(TChain* ttree,TFile* out_file){
80  float pT;
81  float tpT;
82  float track_pT;
83  ttree->SetBranchAddress("photon_pT",&pT);
84  ttree->SetBranchAddress("tphoton_pT",&tpT);
85  ttree->SetBranchAddress("track_pT",&track_pT);
86 
87  TH1F *pTeffPlot = new TH1F("#frac{#Delta#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
88  TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
89  TH1F *trackpTDist = new TH1F("truthpt","",40,0,35);
90  pTeffPlot->Sumw2();
91  pTefffuncPlot->Sumw2();
92  trackpTDist->Sumw2();
93  unsigned lowpTCount=0;
94 
95  for (int event = 0; event < ttree->GetEntries(); ++event)
96  {
97  ttree->GetEvent(event);
98  pTeffPlot->Fill((pT-tpT)/tpT);
99  pTefffuncPlot->Fill(tpT,(pT-tpT)/tpT);
100  trackpTDist->Fill(track_pT);
101  if(tpT<.6)lowpTCount++;
102  }
103  pTeffPlot->Scale(1./ttree->GetEntries(),"width");
104  pTefffuncPlot->Scale(1./ttree->GetEntries(),"width");
105  trackpTDist->Scale(1./ttree->GetEntries(),"width");
106  cout<<"Signal rejection through pT cut= "<<(float)lowpTCount/ttree->GetEntries()<<endl;
107  cout<<"Signal rejection through pT cut= "<<sqrt((float)lowpTCount)/ttree->GetEntries()<<endl;
108  out_file->Write();
109 }
110 void testCuts(TChain* ttree,TFile* out_file){
111  float dphi;
112  float prob;
113  float track_pT;
114  float deta;
115  float radius;
116  int layer;
117  int dlayer;
118  ttree->SetBranchAddress("cluster_dphi",&dphi);
119  ttree->SetBranchAddress("cluster_prob",&prob);
120  ttree->SetBranchAddress("track_layer",&layer);
121  ttree->SetBranchAddress("track_dlayer",&dlayer);
122  ttree->SetBranchAddress("track_pT",&track_pT);
123  ttree->SetBranchAddress("track_deta",&deta);
124  ttree->SetBranchAddress("vtx_radius",&radius);
125 
126  TH1F *layerDist = new TH1F("layer","",16,-.5,15.5);
127  TH1F *probDist = new TH1F("clust_prob","",30,-.5,1.);
128  TH1F *deta_plot = new TH1F("deta","",30,-.001,.01);
129  TH1F *dlayer_plot = new TH1F("dlayer","",11,-.5,10.5);
130  TH1F *r_plot = new TH1F("signal_vtx_radius_dist","",26,-.5,25.5);
131  layerDist->Sumw2();
132  probDist->Sumw2();
133  deta_plot->Sumw2();
134  dlayer_plot->Sumw2();
135  r_plot->Sumw2();
136  unsigned badLayCount=0;
137  unsigned badClusCount=0;
138  unsigned bigDetaCount=0;
139  unsigned shortRadiusCount=0;
140 
141  for (int event = 0; event < ttree->GetEntries(); ++event)
142  {
143  ttree->GetEvent(event);
144  if(layer==0)badLayCount++;
145  if(dphi<0&&track_pT>.6){
146  badClusCount++;
147  }
148  if(track_pT>.6&&dphi>=0){
149  layerDist->Fill(layer);
150  probDist->Fill(prob);
151  deta_plot->Fill(deta);
152  dlayer_plot->Fill(TMath::Abs(dlayer));
153  if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
154  else{
155  r_plot->Fill(radius);
156  if(radius<1.84) shortRadiusCount++;
157  }
158  }
159 
160  }
161  layerDist->Scale(1./ttree->GetEntries());
162  deta_plot->Scale(1./ttree->GetEntries());
163  dlayer_plot->Scale(1./ttree->GetEntries());
164  r_plot->Scale(1./ttree->GetEntries());
165  cout<<"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
166  cout<<"error= "<<sqrt((float)badLayCount)/ttree->GetEntries()<<endl;
167  cout<<"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
168  cout<<"error= "<<sqrt((float)badClusCount)/ttree->GetEntries()<<endl;
169  cout<<"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
170  cout<<"error= "<<sqrt((float)bigDetaCount)/ttree->GetEntries()<<endl;
171  cout<<"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
172  cout<<"error= "<<sqrt((float)shortRadiusCount)/ttree->GetEntries()<<endl;
173  out_file->Write();
174 }
175 void makeRefitDist(TChain* ttree, TFile *out_file){
176  float diffx;
177  float diffy;
178  float diffz;
179  ttree->SetBranchAddress("refitdiffx",&diffx);
180  ttree->SetBranchAddress("refitdiffy",&diffy);
181  ttree->SetBranchAddress("refitdiffz",&diffz);
182  TH1F *diffplotx= new TH1F("x_{0}-x_{refit}","",40,-10,10);
183  TH1F *diffploty= new TH1F("y_{0}-y_{refit}","",40,-10,10);
184  TH1F *diffplotz= new TH1F("z_{0}-z_{refit}","",40,-10,10);
185 
186  diffplotx->Sumw2();
187  diffploty->Sumw2();
188  diffplotz->Sumw2();
189 
190  for (int event = 0; event < ttree->GetEntries(); ++event)
191  {
192  ttree->GetEvent(event);
193  diffplotx->Fill(diffx);
194  diffploty->Fill(diffy);
195  diffplotz->Fill(diffz);
196  }
197 
198  diffplotx->Scale(1./ttree->GetEntries(),"width");
199  diffploty->Scale(1./ttree->GetEntries(),"width");
200  diffplotz->Scale(1./ttree->GetEntries(),"width");
201 
202  out_file->Write();
203 }
204 
205 void observHister()
206 {
207  string treePath = "/sphenix/user/vassalli/minBiasConversion/conversiontruthanaout.root";
208  string treeExtension = ".root";
209  unsigned int nFiles=100;
210  TFile *out_file = new TFile("obsrvplots.root","RECREATE");
211  TChain *ttree = new TChain("observTree");
212  ttree->Add(treePath.c_str());
213  cout<<"Total events= "<<ttree->GetEntries()<<'\n';
214  matchPlot(ttree,out_file);
215 }