EMCal-QA/macro/QA.C

0001 // c++ includes --
0002 #include <string>
0003 #include <vector>
0004 #include <iostream>
0005
0006 // -- root includes --
0007 #include <TH2F.h>
0008 #include <TFile.h>
0009 #include <TTree.h>
0010
0011 using std::cout;
0012 using std::endl;
0013 using std::string;
0014 using std::vector;
0015 using std::min;
0016 using std::max;
0017 using std::to_string;
0018
0019 namespace myAnalysis {
0020
0021     void init(const string& inputFile);
0022     void analyze(UInt_t nevents=0);
0023     void finalize(const string &outputFile="test.root");
0024
0025     TFile* input;
0026     TTree* data;
0027
0028     // QA
0029     TH1F* hclusterEta;
0030     TH1F* hclusterPhi;
0031     TH1F* hclusterE;
0032     TH1F* hclusterE_2;
0033     TH1F* hclusterPt;
0034     TH1F* hclusterChisq;
0035     TH1F* hclusterChisq_2;
0036     TH1F* hclusterChisq_3;
0037
0038     // 2D correlations
0039     TH2F* h2clusterPhiVsEta;
0040     TH2F* h2clusterPhiVsEta_2;
0041     TH2F* h2clusterPhiVsEta_3;
0042     TH2F* h2clusterChisqVsE;
0043     TH2F* h2clusterChisqVsE_2;
0044     TH2F* h2clusterE_CalibVsE;
0045
0046     UInt_t  eta_bins = 48;
0047     Float_t eta_min  = -1.152;
0048     Float_t eta_max  = 1.152;
0049
0050     UInt_t  phi_bins = 126;
0051     Float_t phi_min  = -3.15;
0052     Float_t phi_max  = 3.15;
0053
0054     UInt_t  e_bins  = 500;
0055     Float_t e_min   = 0;
0056     Float_t e_max   = 5e4;
0057
0058     UInt_t  pt_bins  = 500;
0059     Float_t pt_min   = 0;
0060     Float_t pt_max   = 5e4;
0061
0062     UInt_t  chi2_bins  = 500;
0063     Float_t chi2_min   = 0;
0064     Float_t chi2_max   = 5e5;
0065
0066     vector<Float_t>* clusterPhi     = 0;
0067     vector<Float_t>* clusterEta     = 0;
0068     vector<Float_t>* clusterE       = 0;
0069     vector<Float_t>* clusterPt      = 0;
0070     vector<Float_t>* clusterChisq   = 0;
0071 }
0072
0073 void myAnalysis::init(const string& inputFile) {
0074     cout << "Start Init" << endl;
0075     input = TFile::Open(inputFile.c_str());
0076     data  = (TTree*)input->Get("T");
0077
0078     data->SetBranchAddress("clusterPhi",&clusterPhi);
0079     data->SetBranchAddress("clusterEta",&clusterEta);
0080     data->SetBranchAddress("clusterE",&clusterE);
0081     data->SetBranchAddress("clustrPt",&clusterPt);
0082     data->SetBranchAddress("clusterChisq",&clusterChisq);
0083
0084     // QA
0085     hclusterEta       = new TH1F("hclusterEta", "Cluster #eta; #eta; Counts", eta_bins, eta_min, eta_max);
0086     hclusterPhi       = new TH1F("hclusterPhi", "Cluster #phi; #phi; Counts", phi_bins, phi_min, phi_max);
0087     hclusterE         = new TH1F("hclusterE", "Cluster Energy; ADC; Counts", e_bins, e_min, e_max);
0088     hclusterPt        = new TH1F("hclusterPt", "Cluster p_{T}; ADC; Counts", pt_bins, pt_min, pt_max);
0089     hclusterChisq     = new TH1F("hclusterChisq", "Cluster #chi^2; #chi^2; Counts", chi2_bins, chi2_min, chi2_max);
0090     h2clusterPhiVsEta = new TH2F("h2clusterPhiVsEta", "Cluster; #eta; #phi", eta_bins, eta_min, eta_max, phi_bins, phi_min, phi_max);
0091     h2clusterChisqVsE = new TH2F("h2clusterChisqVsE", "Cluster #chi^2 vs E; Energy; #chi^2", e_bins, e_min, e_max, chi2_bins, chi2_min, chi2_max);
0092     h2clusterChisqVsE_2 = new TH2F("h2clusterChisqVsE_2", "Cluster #chi^2 vs E; Energy; #chi^2", e_bins, e_min, e_max, 100, 0, 100);
0093
0094     // QA with cuts
0095     hclusterE_2         = new TH1F("hclusterE_2", "Cluster Energy (#chi^{2} < 10); ADC; Counts", e_bins, e_min, e_max);
0096     hclusterChisq_2     = new TH1F("hclusterChisq_2", "Cluster #chi^{2}; #chi^{2}; Counts", 100, 0, 100);
0097     hclusterChisq_3     = new TH1F("hclusterChisq_3", "Cluster #chi^{2} (500 #leq ADC < 16000); #chi^{2}; Counts", 100, 0, 100);
0098     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);
0099     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);
0100
0101     cout << "Finish Init" << endl;
0102 }
0103
0104 void myAnalysis::analyze(UInt_t nevents) {
0105     cout << "Start Analyze" << endl;
0106
0107     // if nevents is 0 then use all events otherwise use the set number of events
0108     nevents = (nevents) ? nevents : data->GetEntries();
0109
0110     Float_t clusterEta_min     = 0;
0111     Float_t clusterEta_max     = 0;
0112     Float_t clusterPhi_min     = 0;
0113     Float_t clusterPhi_max     = 0;
0114     Float_t clusterE_min       = 0;
0115     Float_t clusterE_max       = 0;
0116     Float_t clusterPt_min      = 0;
0117     Float_t clusterPt_max      = 0;
0118     Float_t clusterChisq_min   = 0;
0119     Float_t clusterChisq_max   = 0;
0120     cout << "Starting Event Loop" << endl;
0121     for(UInt_t i = 0; i < nevents; ++i) {
0122
0123         if(i%500 == 0) cout << "Progress: " << i*100./nevents << " %" << endl;
0124
0125         data->GetEntry(i);
0126
0127         UInt_t nclusters = clusterPhi->size();
0128
0129         for (UInt_t j = 0; j < nclusters; ++j) {
0130             Float_t clusterEta_val     = clusterEta->at(j);
0131             Float_t clusterPhi_val     = clusterPhi->at(j);
0132             Float_t clusterE_val       = clusterE->at(j);
0133             Float_t clusterPt_val      = clusterPt->at(j);
0134             Float_t clusterChisq_val   = clusterChisq->at(j);
0135
0136             clusterEta_min     = min(clusterEta_min, clusterEta_val);
0137             clusterEta_max     = max(clusterEta_max, clusterEta_val);
0138             clusterPhi_min     = min(clusterPhi_min, clusterPhi_val);
0139             clusterPhi_max     = max(clusterPhi_max, clusterPhi_val);
0140             clusterE_min       = min(clusterE_min, clusterE_val);
0141             clusterE_max       = max(clusterE_max, clusterE_val);
0142             clusterPt_min      = min(clusterPt_min, clusterPt_val);
0143             clusterPt_max      = max(clusterPt_max, clusterPt_val);
0144             clusterChisq_min   = min(clusterChisq_min, clusterChisq_val);
0145             clusterChisq_max   = max(clusterChisq_max, clusterChisq_val);
0146
0147             hclusterEta->Fill(clusterEta_val);
0148             hclusterPhi->Fill(clusterPhi_val);
0149
0150             h2clusterPhiVsEta->Fill(clusterEta_val, clusterPhi_val);
0151
0152             hclusterE->Fill(clusterE_val);
0153             hclusterPt->Fill(clusterPt_val);
0154             hclusterChisq->Fill(clusterChisq_val);
0155             hclusterChisq_2->Fill(clusterChisq_val);
0156
0157             h2clusterChisqVsE->Fill(clusterE_val, clusterChisq_val);
0158             h2clusterChisqVsE_2->Fill(clusterE_val, clusterChisq_val);
0159
0160             if(clusterChisq_val < 10) {
0161                 hclusterE_2->Fill(clusterE_val);
0162             }
0163
0164             if(clusterE_val >= 500 && clusterE_val < 16000) {
0165                 hclusterChisq_3->Fill(clusterChisq_val);
0166             }
0167
0168             if(clusterChisq_val < 10 && clusterE_val >= 500 && clusterE_val < 16000) {
0169                 h2clusterPhiVsEta_2->Fill(clusterEta_val, clusterPhi_val);
0170             }
0171
0172             if(clusterChisq_val < 10 && clusterE_val >= 8.5e3 && clusterE_val < 9.5e3) {
0173                 h2clusterPhiVsEta_3->Fill(clusterEta_val, clusterPhi_val);
0174             }
0175         }
0176     }
0177
0178     cout << "events processed: " << nevents << endl;
0179     cout << "clusterEta_min: " << clusterEta_min << " clusterEta_max: " << clusterEta_max << endl;
0180     cout << "clusterPhi_min: " << clusterPhi_min << " clusterPhi_max: " << clusterPhi_max << endl;
0181     cout << "clusterE_min: " << clusterE_min << " clusterE_max: " << clusterE_max << endl;
0182     cout << "clusterPt_min: " << clusterPt_min << " clusterPt_max: " << clusterPt_max << endl;
0183     cout << "clusterChisq_min: " << clusterChisq_min << " clusterChisq_max: " << clusterChisq_max << endl;
0184     cout << "Finish Analyze" << endl;
0185 }
0186
0187 void myAnalysis::finalize(const string& outputFile) {
0188     TFile output(outputFile.c_str(),"recreate");
0189     output.cd();
0190
0191     hclusterEta->Write();
0192     hclusterPhi->Write();
0193     h2clusterPhiVsEta->Write();
0194     hclusterE->Write();
0195     hclusterPt->Write();
0196     hclusterChisq->Write();
0197     h2clusterChisqVsE->Write();
0198     h2clusterChisqVsE_2->Write();
0199
0200     h2clusterPhiVsEta_2->Write();
0201     h2clusterPhiVsEta_3->Write();
0202     hclusterE_2->Write();
0203     hclusterChisq_2->Write();
0204     hclusterChisq_3->Write();
0205
0206     input->Close();
0207     output.Close();
0208 }
0209
0210 void QA(const string &inputFile,
0211         const string &outputFile="test.root",
0212         UInt_t nevents=0) {
0213
0214     myAnalysis::init(inputFile);
0215     myAnalysis::analyze(nevents);
0216     myAnalysis::finalize(outputFile);
0217 }
0218
0219 # ifndef __CINT__
0220 int main(int argc, char* argv[]) {
0221     if(argc < 2 || argc > 4){
0222         cout << "usage: ./bin/QA inputFile outputFile events" << endl;
0223         return 1;
0224     }
0225
0226     string input;
0227     string output = "test.root";
0228     UInt_t events = 0;
0229
0230     if(argc >= 2) {
0231         input = argv[1];
0232     }
0233     if(argc >= 3) {
0234         output = argv[2];
0235     }
0236     if(argc >= 4) {
0237         events = atoi(argv[3]);
0238     }
0239
0240     QA(input, output, events);
0241
0242     cout << "done" << endl;
0243     return 0;
0244 }
0245 # endif