Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
minBiasRecoAna.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file minBiasRecoAna.C
1 #include <iostream>
2 using namespace std;
3 
4 #include "TFile.h"
5 #include "TTree.h"
6 #include "TChain.h"
7 #include "TLegend.h"
8 #include "math.h"
9 #include "TH1.h"
10 #include "TH2.h"
11 #include "TEfficiency.h"
12 
13 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
14  TChain *all = new TChain(treename.c_str());
15  string temp;
16  for (int i = 0; i < filecount; ++i)
17  {
18 
19  ostringstream s;
20  s<<i;
21  temp = name+string(s.str())+extension;
22  all->Add(temp.c_str());
23  }
24  return all;
25 }
26 
27 int reportBackground(TChain* truthTree,TChain* recoTree,TFile* out_file){
28  double signal = truthTree->GetEntries();
29  signal*=(1-.012)*(1-.002)*(1-.011);
30  int background = recoTree->GetEntries();
31  background-=(int) signal;
32  cout<<"For "<<truthTree->GetEntries()<<" truth conversions, I found "<<background<<" background events"<<endl;
33  return (int) signal;
34  /*float pT;
35  float tpT;
36  ttree->SetBranchAddress("photon_pT",&pT);
37  ttree->SetBranchAddress("tphoton_pT",&tpT);
38 
39  TH1F *pTeffPlot = new TH1F("#frac{#Delta#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
40  TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
41  pTeffPlot->Sumw2();
42 
43  for (int event = 0; event < ttree->GetEntries(); ++event)
44  {
45  ttree->GetEvent(event);
46  pTeffPlot->Fill((pT-tpT)/tpT);
47  pTefffuncPlot->Fill(tpT,(pT-tpT)/tpT);
48  }
49  pTeffPlot->Scale(1./ttree->GetEntries(),"width");
50  pTefffuncPlot->Scale(1./ttree->GetEntries(),"width");
51  out_file->Write();*/
52 }
53 
54 void reportCuts(TChain* _treeBackground,int signal){
55  unsigned totalTracks;
56  unsigned passedpTEtaQ;
57  unsigned passedCluster;
58  unsigned passedPair;
59  unsigned passedVtx;
60  unsigned passedPhoton;
61 
62  int sum_totalTracks=-1*signal;
63  int sum_passedpTEtaQ=-1*signal;
64  int sum_passedCluster=-1*signal;
65  int sum_passedPair=-1*signal;
66  int sum_passedVtx=-1*signal;
67  int sum_passedPhoton=-1*signal;
68 
69  _treeBackground->SetBranchAddress("total", &totalTracks);
70  _treeBackground->SetBranchAddress("tracksPEQ", &passedpTEtaQ);
71  _treeBackground->SetBranchAddress("tracks_clus", &passedCluster);
72  _treeBackground->SetBranchAddress("pairs", &passedPair);
73  _treeBackground->SetBranchAddress("vtx", &passedVtx);
74  _treeBackground->SetBranchAddress("photon", &passedPhoton);
75 
76  for(unsigned i=0; i<_treeBackground->GetEntries();i++){
77  _treeBackground->GetEntry(i);
78  sum_totalTracks+=totalTracks;
79  sum_passedpTEtaQ+=passedpTEtaQ;
80  sum_passedCluster+=passedCluster;
81  sum_passedPair+=passedPair;
82  sum_passedVtx+=passedVtx;
83  sum_passedPhoton+=passedPhoton;
84  }
85  cout<<"Did RecoConversionEval with "<<sum_totalTracks<<" total tracks\n\t";
86  cout<<1-(float)sum_passedpTEtaQ/sum_totalTracks<<"+/-"<<sqrt((float)sum_passedpTEtaQ)/sum_totalTracks<<" of tracks were rejected by pTEtaQ\n\t";
87  cout<<1-(float)sum_passedCluster/sum_passedpTEtaQ<<"+/-"<<sqrt((float)sum_passedCluster)/sum_passedpTEtaQ<<" of remaining tracks were rejected by cluster\n\t";
88  cout<<1-(float)sum_passedPair/sum_passedCluster<<"+/-"<<sqrt((float)sum_passedPair)/sum_passedCluster<<" of pairs were rejected by pair cuts\n\t";
89  cout<<1-(float)sum_passedVtx/sum_passedPair<<"+/-"<<sqrt((float)sum_passedVtx)/sum_passedPair<<" of vtx were rejected by vtx cuts\n\t";
90  cout<<1-(float)sum_passedPhoton/sum_passedVtx<<"+/-"<<sqrt((float)sum_passedPhoton)/sum_passedVtx<<" of photons were rejected by photon cuts\n";
91  cout<<sum_passedPhoton<<" background events remain with "<<signal<<" signal events\n";
92 }
93 
94 int analyzeSignal(TChain* _signalCutTree,TFile *outFile){
95  int _b_track_layer ;
96  int _b_track_dlayer ;
97  float _b_track_deta ;
98  float _b_track_pT;
99  float _b_ttrack_pT;
100  float _b_vtx_radius ;
101  float _b_tvtx_radius ;
102  float _b_cluster_prob ;
103  float _b_photon_pT ;
104  float _b_tphoton_pT ;
105  float _b_photon_m ;
106  _signalCutTree->SetBranchAddress("track_deta", &_b_track_deta);
107  _signalCutTree->SetBranchAddress("track_dlayer",&_b_track_dlayer);
108  // _signalCutTree->SetBranchAddress("track_layer", &_b_track_layer);
109  _signalCutTree->SetBranchAddress("track_pT", &_b_track_pT);
110  //_signalCutTree->SetBranchAddress("ttrack_pT", &_b_ttrack_pT);
111  _signalCutTree->SetBranchAddress("vtx_radius", &_b_vtx_radius);
112  //_signalCutTree->SetBranchAddress("tvtx_radius", &_b_tvtx_radius);
113  // _signalCutTree->SetBranchAddress("vtx_chi2", &_b_vtx_chi2);
114  _signalCutTree->SetBranchAddress("cluster_prob", &_b_cluster_prob);
115  _signalCutTree->SetBranchAddress("photon_m", &_b_photon_m);
116  _signalCutTree->SetBranchAddress("photon_pT", &_b_photon_pT);
117  _signalCutTree->SetBranchAddress("tphoton_pT", &_b_tphoton_pT);
118  unsigned pT=0, cluster=0, eta=0, photon=0,rsignal=0;
119 
120  TH1F *rate_plot = new TH1F("rateplot","",60,0,15);
121  rate_plot->Sumw2();
122 
123  for(unsigned i=0; i<_signalCutTree->GetEntries();i++){
124  _signalCutTree->GetEntry(i);
125  if(_b_track_pT<.6) pT++;
126  else if (_b_cluster_prob<0) cluster++;
127  else if (_b_track_deta>.0082) eta++;
128  else if (TMath::Abs(_b_track_dlayer)<=9&&_b_vtx_radius>1.84){
129  if(_b_photon_m<.27||_b_photon_m>8.||_b_photon_pT<.039) photon++;
130  else rsignal++;
131  }
132  rate_plot->Fill(_b_tphoton_pT);
133  }
134  cout<<"pT cut "<<pT<<" events\n";
135  cout<<"cluster cut "<<cluster<<" events\n";
136  cout<<"eta cut "<<eta<<" events\n";
137  cout<<"photon cut "<<photon<<" events\n";
138  rate_plot->Scale(1./6e5);
139  outFile->Wrtie();
140  return rsignal;
141 }
142 
143 void minBiasRecoAna()
144 {
145  string treePath = "/sphenix/user/vassalli/minBiasPythia/conversionembededminBiasanalysis";
146  string truthTreePath = "/sphenix/user/vassalli/minBiasPythia/conversionembededminBiasTruthanalysis";
147  string embedPath = "/sphenix/user/vassalli/RecoConversionTests/conversionembededonlineanalysis";
148  string treeExtension = ".root";
149  unsigned int nFiles=100;
150  TFile *out_file = new TFile("minBiasplots.root","RECREATE");
151  TChain *embed_truth_ttree = handleFile(embedPath,treeExtension,"cutTreeSignal",nFiles);
152  TChain *embed_reco_ttree = handleFile(embedPath,treeExtension,"recoBackground",nFiles);
153  TChain *cut_tree = new TChain("recoBackground");
154  TChain *back_tree = new TChain("recoSignal");
155  TChain *truth_ttree = new TChain("cutTreeSignal");
156  string filename = "/sphenix/user/vassalli/minBiasConversion/conversionanaout.root";
157  back_tree->Add(filename.c_str());
158  cut_tree->Add(filename.c_str());
159  filename = "/sphenix/user/vassalli/minBiasConversion/conversiontruthanaout.root";
160  truth_ttree->Add(filename.c_str());
161  reportCuts(cut_tree,analyzeSignal(truth_ttree));
162 /* cout<<"///////////////////////////////////////////////////////////\n";
163  cout<<"EMBEDED ANALYSIS:\n";
164  reportCuts(embed_reco_ttree,analyzeSignal(embed_truth_ttree));
165 */
166 }