Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QA.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QA.C
1 // c++ includes --
2 #include <string>
3 #include <vector>
4 #include <iostream>
5 
6 // -- root includes --
7 #include <TH2F.h>
8 #include <TFile.h>
9 #include <TTree.h>
10 
11 using std::cout;
12 using std::endl;
13 using std::string;
14 using std::vector;
15 using std::min;
16 using std::max;
17 using std::to_string;
18 
19 namespace myAnalysis {
20 
21  void init(const string& inputFile);
22  void analyze(UInt_t nevents=0);
23  void finalize(const string &outputFile="test.root");
24 
25  TFile* input;
26  TTree* data;
27 
28  // QA
29  TH1F* hclusterEta;
30  TH1F* hclusterPhi;
31  TH1F* hclusterE;
32  TH1F* hclusterE_2;
33  TH1F* hclusterPt;
34  TH1F* hclusterChisq;
35  TH1F* hclusterChisq_2;
36  TH1F* hclusterChisq_3;
37 
38  // 2D correlations
39  TH2F* h2clusterPhiVsEta;
40  TH2F* h2clusterPhiVsEta_2;
41  TH2F* h2clusterPhiVsEta_3;
42  TH2F* h2clusterChisqVsE;
43  TH2F* h2clusterChisqVsE_2;
44  TH2F* h2clusterE_CalibVsE;
45 
46  UInt_t eta_bins = 48;
47  Float_t eta_min = -1.152;
48  Float_t eta_max = 1.152;
49 
50  UInt_t phi_bins = 126;
51  Float_t phi_min = -3.15;
52  Float_t phi_max = 3.15;
53 
54  UInt_t e_bins = 500;
55  Float_t e_min = 0;
56  Float_t e_max = 5e4;
57 
58  UInt_t pt_bins = 500;
59  Float_t pt_min = 0;
60  Float_t pt_max = 5e4;
61 
62  UInt_t chi2_bins = 500;
63  Float_t chi2_min = 0;
64  Float_t chi2_max = 5e5;
65 
66  vector<Float_t>* clusterPhi = 0;
67  vector<Float_t>* clusterEta = 0;
68  vector<Float_t>* clusterE = 0;
69  vector<Float_t>* clusterPt = 0;
70  vector<Float_t>* clusterChisq = 0;
71 }
72 
73 void myAnalysis::init(const string& inputFile) {
74  cout << "Start Init" << endl;
75  input = TFile::Open(inputFile.c_str());
76  data = (TTree*)input->Get("T");
77 
78  data->SetBranchAddress("clusterPhi",&clusterPhi);
79  data->SetBranchAddress("clusterEta",&clusterEta);
80  data->SetBranchAddress("clusterE",&clusterE);
81  data->SetBranchAddress("clustrPt",&clusterPt);
82  data->SetBranchAddress("clusterChisq",&clusterChisq);
83 
84  // QA
85  hclusterEta = new TH1F("hclusterEta", "Cluster #eta; #eta; Counts", eta_bins, eta_min, eta_max);
86  hclusterPhi = new TH1F("hclusterPhi", "Cluster #phi; #phi; Counts", phi_bins, phi_min, phi_max);
87  hclusterE = new TH1F("hclusterE", "Cluster Energy; ADC; Counts", e_bins, e_min, e_max);
88  hclusterPt = new TH1F("hclusterPt", "Cluster p_{T}; ADC; Counts", pt_bins, pt_min, pt_max);
89  hclusterChisq = new TH1F("hclusterChisq", "Cluster #chi^2; #chi^2; Counts", chi2_bins, chi2_min, chi2_max);
90  h2clusterPhiVsEta = new TH2F("h2clusterPhiVsEta", "Cluster; #eta; #phi", eta_bins, eta_min, eta_max, phi_bins, phi_min, phi_max);
91  h2clusterChisqVsE = new TH2F("h2clusterChisqVsE", "Cluster #chi^2 vs E; Energy; #chi^2", e_bins, e_min, e_max, chi2_bins, chi2_min, chi2_max);
92  h2clusterChisqVsE_2 = new TH2F("h2clusterChisqVsE_2", "Cluster #chi^2 vs E; Energy; #chi^2", e_bins, e_min, e_max, 100, 0, 100);
93 
94  // QA with cuts
95  hclusterE_2 = new TH1F("hclusterE_2", "Cluster Energy (#chi^{2} < 10); ADC; Counts", e_bins, e_min, e_max);
96  hclusterChisq_2 = new TH1F("hclusterChisq_2", "Cluster #chi^{2}; #chi^{2}; Counts", 100, 0, 100);
97  hclusterChisq_3 = new TH1F("hclusterChisq_3", "Cluster #chi^{2} (500 #leq ADC < 16000); #chi^{2}; Counts", 100, 0, 100);
98  h2clusterPhiVsEta_2 = new TH2F("h2clusterPhiVsEta_2", "Cluster (500 #leq ADC < 16000 and #chi^{2} < 10); #eta; #phi", eta_bins, eta_min, eta_max, phi_bins, phi_min, phi_max);
99  h2clusterPhiVsEta_3 = new TH2F("h2clusterPhiVsEta_3", "Cluster (8500 #leq ADC < 9500 and #chi^{2} < 10); #eta; #phi", 96, eta_min, eta_max, 256, phi_min, phi_max);
100 
101  cout << "Finish Init" << endl;
102 }
103 
104 void myAnalysis::analyze(UInt_t nevents) {
105  cout << "Start Analyze" << endl;
106 
107  // if nevents is 0 then use all events otherwise use the set number of events
108  nevents = (nevents) ? nevents : data->GetEntries();
109 
110  Float_t clusterEta_min = 0;
111  Float_t clusterEta_max = 0;
112  Float_t clusterPhi_min = 0;
113  Float_t clusterPhi_max = 0;
114  Float_t clusterE_min = 0;
115  Float_t clusterE_max = 0;
116  Float_t clusterPt_min = 0;
117  Float_t clusterPt_max = 0;
118  Float_t clusterChisq_min = 0;
119  Float_t clusterChisq_max = 0;
120  cout << "Starting Event Loop" << endl;
121  for(UInt_t i = 0; i < nevents; ++i) {
122 
123  if(i%500 == 0) cout << "Progress: " << i*100./nevents << " %" << endl;
124 
125  data->GetEntry(i);
126 
127  UInt_t nclusters = clusterPhi->size();
128 
129  for (UInt_t j = 0; j < nclusters; ++j) {
130  Float_t clusterEta_val = clusterEta->at(j);
131  Float_t clusterPhi_val = clusterPhi->at(j);
132  Float_t clusterE_val = clusterE->at(j);
133  Float_t clusterPt_val = clusterPt->at(j);
134  Float_t clusterChisq_val = clusterChisq->at(j);
135 
136  clusterEta_min = min(clusterEta_min, clusterEta_val);
137  clusterEta_max = max(clusterEta_max, clusterEta_val);
138  clusterPhi_min = min(clusterPhi_min, clusterPhi_val);
139  clusterPhi_max = max(clusterPhi_max, clusterPhi_val);
140  clusterE_min = min(clusterE_min, clusterE_val);
141  clusterE_max = max(clusterE_max, clusterE_val);
142  clusterPt_min = min(clusterPt_min, clusterPt_val);
143  clusterPt_max = max(clusterPt_max, clusterPt_val);
144  clusterChisq_min = min(clusterChisq_min, clusterChisq_val);
145  clusterChisq_max = max(clusterChisq_max, clusterChisq_val);
146 
147  hclusterEta->Fill(clusterEta_val);
148  hclusterPhi->Fill(clusterPhi_val);
149 
150  h2clusterPhiVsEta->Fill(clusterEta_val, clusterPhi_val);
151 
152  hclusterE->Fill(clusterE_val);
153  hclusterPt->Fill(clusterPt_val);
154  hclusterChisq->Fill(clusterChisq_val);
155  hclusterChisq_2->Fill(clusterChisq_val);
156 
157  h2clusterChisqVsE->Fill(clusterE_val, clusterChisq_val);
158  h2clusterChisqVsE_2->Fill(clusterE_val, clusterChisq_val);
159 
160  if(clusterChisq_val < 10) {
161  hclusterE_2->Fill(clusterE_val);
162  }
163 
164  if(clusterE_val >= 500 && clusterE_val < 16000) {
165  hclusterChisq_3->Fill(clusterChisq_val);
166  }
167 
168  if(clusterChisq_val < 10 && clusterE_val >= 500 && clusterE_val < 16000) {
169  h2clusterPhiVsEta_2->Fill(clusterEta_val, clusterPhi_val);
170  }
171 
172  if(clusterChisq_val < 10 && clusterE_val >= 8.5e3 && clusterE_val < 9.5e3) {
173  h2clusterPhiVsEta_3->Fill(clusterEta_val, clusterPhi_val);
174  }
175  }
176  }
177 
178  cout << "events processed: " << nevents << endl;
179  cout << "clusterEta_min: " << clusterEta_min << " clusterEta_max: " << clusterEta_max << endl;
180  cout << "clusterPhi_min: " << clusterPhi_min << " clusterPhi_max: " << clusterPhi_max << endl;
181  cout << "clusterE_min: " << clusterE_min << " clusterE_max: " << clusterE_max << endl;
182  cout << "clusterPt_min: " << clusterPt_min << " clusterPt_max: " << clusterPt_max << endl;
183  cout << "clusterChisq_min: " << clusterChisq_min << " clusterChisq_max: " << clusterChisq_max << endl;
184  cout << "Finish Analyze" << endl;
185 }
186 
187 void myAnalysis::finalize(const string& outputFile) {
188  TFile output(outputFile.c_str(),"recreate");
189  output.cd();
190 
191  hclusterEta->Write();
192  hclusterPhi->Write();
193  h2clusterPhiVsEta->Write();
194  hclusterE->Write();
195  hclusterPt->Write();
196  hclusterChisq->Write();
197  h2clusterChisqVsE->Write();
198  h2clusterChisqVsE_2->Write();
199 
200  h2clusterPhiVsEta_2->Write();
201  h2clusterPhiVsEta_3->Write();
202  hclusterE_2->Write();
203  hclusterChisq_2->Write();
204  hclusterChisq_3->Write();
205 
206  input->Close();
207  output.Close();
208 }
209 
210 void QA(const string &inputFile,
211  const string &outputFile="test.root",
212  UInt_t nevents=0) {
213 
214  myAnalysis::init(inputFile);
215  myAnalysis::analyze(nevents);
216  myAnalysis::finalize(outputFile);
217 }
218 
219 # ifndef __CINT__
220 int main(int argc, char* argv[]) {
221  if(argc < 2 || argc > 4){
222  cout << "usage: ./bin/QA inputFile outputFile events" << endl;
223  return 1;
224  }
225 
226  string input;
227  string output = "test.root";
228  UInt_t events = 0;
229 
230  if(argc >= 2) {
231  input = argv[1];
232  }
233  if(argc >= 3) {
234  output = argv[2];
235  }
236  if(argc >= 4) {
237  events = atoi(argv[3]);
238  }
239 
240  QA(input, output, events);
241 
242  cout << "done" << endl;
243  return 0;
244 }
245 # endif