EMCal-commissioning/packetCorrelation/draw_corr.C

0001 #include "TFile.h"
0002 #include "TTree.h"
0003
0004 #include "TStyle.h"
0005 #include "TLatex.h"
0006 #include "TH1D.h"
0007 #include "TH2D.h"
0008 #include "TCanvas.h"
0009
0010 #include <iostream>
0011 #include <iomanip>
0012 #include <vector>
0013 #include <tuple>
0014 #include <string>
0015 #include <cmath>
0016 #include <algorithm>
0017
0018
0019 int round8(float num) {
0020       int inum = static_cast<int>(num);
0021         return inum - inum % 8;
0022 }
0023
0024 void draw_corr() {
0025     const char * run = "../9613-0000.root";
0026     const char * runNum = "EMCal run 9613";
0027
0028     // Are your tree variables named correctly??
0029
0030     TFile *f = new TFile(run,"READ");
0031
0032     f->ls();
0033
0034     TTree *t = (TTree*) f->Get("towerntup");
0035     const int entries = t->GetEntries();
0036
0037
0038     t->Show(0);
0039
0040     std::vector<float> *emcTowE = 0;
0041     std::vector<float> *emciEta = 0;
0042     std::vector<float> *emciPhi = 0;
0043
0044     t->SetBranchAddress("energy",&emcTowE);
0045     t->SetBranchAddress("etabin",&emciEta);
0046     t->SetBranchAddress("phibin",&emciPhi);
0047
0048     TH2D *EDist = new TH2D("EDist",";EMCal tower phi;EMCal tower eta",256,-.5,255.5,96,-.5,95.5);
0049     TH2D *EDist_gap = new TH2D("EDist_gap",";EMCal tower phi;EMCal tower eta",256,-.5,255.5,96,-.5,95.5);
0050     TH2D* packet_corrs[128][128];
0051     for (int i = 0; i < 128; i++) {
0052         for (int j = 0; j < 128; j++) {
0053             std::string corrIndex = std::to_string(i) + "-" + std::to_string(j);
0054             const char * ccorrIndex = corrIndex.c_str();
0055             packet_corrs[i][j] = new TH2D(ccorrIndex,";Energy;Energy",256,0,20000,256,0,20000);
0056         }
0057     }
0058
0059     auto packetE = new float[entries][128];
0060     TH2D *correlations = new TH2D("correlations",";packets",128,0,128,128,0,128);
0061
0062     float tower_E[256 * 96] = { 0 }; // variables used for average energies
0063     float tower_ET[256 * 96] = { 0 };
0064
0065     TCanvas *tc = new TCanvas();
0066     gStyle->SetOptStat(0);
0067     TLatex l;
0068     l.SetNDC();
0069     l.SetTextSize(0.03);
0070
0071
0072
0073     std::cout << "I see " << entries << " events!" << std::endl;
0074
0075     for (int e = 0 ; e < entries; e++ ) {
0076       t->GetEntry( e );
0077
0078       for (unsigned i = 0; i < emcTowE->size(); i++) {
0079            tower_E[i] += emcTowE->at(i);
0080       }
0081
0082       if ( e % 1000 == 0 ) {
0083         std::cout << "\r" << static_cast<int>(static_cast<float>(e) / entries * 100) << "%" << std::flush;
0084       }
0085     }
0086     std::cout << endl;
0087
0088     // Find the various hot/cold spots
0089     std::cout << endl << "Finding hot/cold spots..." << endl;
0090
0091     TH1D *h_tower_E = new TH1D("h_tower_E",";Tower energy;counts",200,0,200000);
0092     TH1D *h_tower_E_top = new TH1D("h_tower_E_top",";Tower energy;counts",200,0,200000);
0093
0094     float sumE = 0.0;
0095     float hottest = 0.0;
0096     for (int i = 0; i < 256 * 96; i++) {
0097         if (emciEta->at(i) >= 80) {
0098            h_tower_E_top->Fill(tower_E[i]);
0099         }
0100         h_tower_E->Fill(tower_E[i]);
0101         sumE += tower_E[i];
0102         if (tower_E[i] > hottest) {
0103             hottest = tower_E[i];
0104         }
0105     }
0106     float aveE = sumE / (256 * 96);
0107     float stdDev = h_tower_E->GetStdDev();
0108
0109     std::vector<int> hot_spots;
0110     std::vector<int> warm_spots;
0111     std::vector<int> cold_spots;
0112     std::vector<int> cool_spots;
0113     std::vector<int> normal_spots;
0114
0115     float cool_tol = 0.65;
0116
0117     for (int i = 0; i < 256 * 96; i++) {
0118
0119         if (tower_E[i] > (aveE+stdDev*3) && tower_E[i + 1] < (aveE+stdDev*3) && tower_E[i - 1] < (aveE+stdDev*3)) {
0120            hot_spots.push_back(i);
0121         }
0122         else if (tower_E[i] > (aveE + stdDev*3)) {
0123            warm_spots.push_back(i);
0124         }
0125         else if (tower_E[i] < aveE * cool_tol && tower_E[i + 1] > aveE * cool_tol && tower_E[i - 1] > aveE * cool_tol) {
0126            cold_spots.push_back(i);
0127         }
0128         else if (tower_E[i] < aveE * cool_tol) {
0129            cool_spots.push_back(i);
0130         }
0131         else {
0132            normal_spots.push_back(i);
0133         }
0134     }
0135
0136     gPad->SetLogy(1);
0137     h_tower_E->Draw();
0138     l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0139     l.DrawLatex(0.15,0.91,runNum);
0140     tc->Print("tower_E_1D.pdf");
0141
0142     h_tower_E_top->Draw();
0143     l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0144     l.DrawLatex(0.15,0.91,runNum);
0145     tc->Print("tower_E_top_1D.pdf");
0146     gPad->SetLogy(0);
0147
0148
0149     // Print out the locations of the very hot/cold spots
0150
0151      std::cout << "hot spots" << endl;
0152      for (unsigned i = 0; i < hot_spots.size(); i++) {
0153        std::cout << "(ieta, iphi) = (";
0154        std::cout << emciEta->at(hot_spots.at(i)) << ", " << emciPhi->at(hot_spots.at(i)) << ") " << endl;
0155      }
0156      std::cout << endl << "warm spots" << endl;
0157      for (unsigned i = 0; i < warm_spots.size(); i++ ) {
0158          std::cout << "(ieta, iphi) = (";
0159          std::cout << emciEta->at(warm_spots.at(i)) << ", " << emciPhi->at(warm_spots.at(i)) << ") " << endl;
0160      }
0161      std::cout << endl << "cold spots" << endl;
0162      for (unsigned i = 0; i < cold_spots.size(); i++) {
0163        std::cout << "(ieta, iphi) = (";
0164        std::cout << emciEta->at(cold_spots.at(i)) << ", " << emciPhi->at(cold_spots.at(i)) << ") " << endl;
0165      }
0166      std::cout << endl << "cool spots" << endl;
0167      std::vector<std::tuple<int,int>> cool_points;
0168
0169      for (unsigned i = 0; i < cool_spots.size(); i++) {
0170        int rEta = round8(emciEta->at(cool_spots.at(i)));
0171        int rPhi = round8(emciPhi->at(cool_spots.at(i)));
0172
0173        if (rEta != 0 &&
0174            std::find(cool_points.begin(), cool_points.end(), std::make_tuple(rEta,rPhi)) == cool_points.end()) {
0175          cool_points.push_back(std::make_tuple(rEta, rPhi));
0176          std::cout << "(ieta, iphi) = (" << rEta << ", " << rPhi << ")" << endl;
0177        }
0178      }
0179
0180      bool skip = false;
0181      for (unsigned i = 0; i < emcTowE->size(); i++ ) {
0182
0183        EDist->Fill( emciPhi->at(i), emciEta->at(i), tower_E[i] );
0184
0185        for (unsigned j = 0; j < hot_spots.size(); j++ ) {
0186           if (i == static_cast<unsigned>(hot_spots.at(j))) {
0187              skip = true;
0188            }
0189         }
0190         for (unsigned j = 0; j < warm_spots.size(); j++ ) {
0191           if (i == static_cast<unsigned>(warm_spots.at(j))) {
0192              skip = true;
0193           }
0194         }
0195
0196         if (!skip) {
0197           EDist_gap->Fill(emciPhi->at(i), emciEta->at(i), tower_E[i]);
0198         }
0199         skip = false;
0200     }
0201
0202     gPad->SetLogz(1);
0203     EDist->Draw("colz");
0204     l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0205     l.DrawLatex(0.15,0.91,runNum);
0206     l.SetTextSize(0.035);
0207     l.DrawLatex(0.65,0.92,"#Sigma_{evt} ADC");
0208     l.SetTextSize(0.03);
0209     tc->Print("EDist.pdf");
0210     gPad->SetLogz(0);
0211
0212     EDist_gap->Draw("colz");
0213     l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0214     l.DrawLatex(0.15,0.91,runNum);
0215     l.SetTextSize(0.035);
0216     l.DrawLatex(0.65,0.92,"#Sigma_{evt} ADC");
0217     l.SetTextSize(0.03);
0218     tc->Print("EDistGap.pdf");
0219
0220
0221     // packet correlations
0222
0223     std::cout << endl << "Finding packet energies..." << endl;
0224     bool hot = false;
0225     int ieta = 0;
0226     int iphi = 0;
0227     int packet = 0;
0228     for (int e = 0; e < entries; e++) {
0229
0230         t->GetEntry(e);
0231
0232         for (unsigned i = 0; i < emcTowE->size(); i++) {
0233             ieta = emciEta->at(i);
0234             iphi = emciPhi->at(i);
0235             for (unsigned j = 0; j < hot_spots.size(); j++ ){
0236                 if (ieta == emciEta->at(hot_spots.at(j)) && iphi == emciPhi->at(hot_spots.at(j))) {
0237                     hot = true;
0238                 }
0239             }
0240             for (unsigned j = 0; j < warm_spots.size(); j++ ) {
0241                 if (ieta == emciEta->at(warm_spots.at(j)) && iphi == emciPhi->at(warm_spots.at(j))) {
0242                     hot = true;
0243                 }
0244             }
0245
0246             packet += iphi / 8 * 4;
0247             if (ieta / 24 == 0) { packet += ieta / 24 + 1; }
0248             else if (ieta / 24 == 1) { packet += ieta / 24 - 1; }
0249             else { packet += ieta / 24; }
0250
0251             if (!hot) { packetE[e][packet] += emcTowE->at(i); } // 24 = 96/4 and 8 = 256/32
0252             hot = false;
0253             packet = 0;
0254         }
0255         if ( e % 1000 == 0 ) {
0256            std::cout << "\r" << static_cast<int>(static_cast<float>(e) / entries * 100) << "%" << std::flush;
0257         }
0258     }
0259     std::cout << endl;
0260
0261     std::cout << endl << "Finding correlations..." << endl;
0262
0263     for (int e = 0; e < entries; e++ ) {
0264         for (int i = 0; i < 128; i++) {
0265             for (int j = 0; j < 128; j++) {
0266                 packet_corrs[i][j]->Fill(packetE[e][i],packetE[e][j]);
0267                 correlations->Fill(i,j,packet_corrs[i][j]->GetCorrelationFactor() / entries);
0268             }
0269         }
0270         if ( e % 1000 == 0 ) {
0271            std::cout << "\r" << static_cast<int>(static_cast<float>(e) / entries * 100) << "%" << std::flush;
0272         }
0273     }
0274     std::cout << endl;
0275
0276
0277     correlations->SetMinimum(.7);
0278     correlations->Draw("colz");
0279     l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0280     l.DrawLatex(0.15,0.91,runNum);
0281     l.SetTextSize(0.035);
0282     l.DrawLatex(0.65,0.92,"Packet Correlations");
0283     l.SetTextSize(0.03);
0284     tc->Print("packet_correlations.pdf");
0285
0286     bool stop = false;
0287     std::string quit = "";
0288     std::cout << "Would you like to see a specific correlation? (yes/no)" << endl;
0289     std::cin >> quit;
0290     if (quit == "no") { stop = true; }
0291     while (!stop) {
0292         std::cout << "Which correlation would you like to see? (i j)" << endl;
0293         int ans1 = 0;
0294         int ans2 = 0;
0295         std::cin >> ans1 >> ans2;
0296
0297         std::string title = "Packet Correlations " + std::to_string(ans1) + "-" + std::to_string(ans2);
0298         std::cout << title << endl;
0299         const char * ctitle = title.c_str();
0300
0301         packet_corrs[ans1][ans2]->Draw("colz");
0302         l.DrawLatex(0.15,0.95,"sPHENIX Internal");
0303         l.DrawLatex(0.15,0.91,runNum);
0304         l.SetTextSize(0.035);
0305         l.DrawLatex(0.65,0.92,ctitle);
0306         l.SetTextSize(0.03);
0307         tc->Print("specific_correlation.pdf");
0308
0309         std::cout << "Would you like to keep going? (yes/no) " << endl;
0310         std::cin >> quit;
0311         if (quit == "no") {
0312             stop = true;
0313         }
0314     }
0315 }