HcalSingleHadronCalib/macros/TrackCaloDist.C

0001 #include <iostream>
0002 #include <TH2D.h>
0003 #include <TH1D.h>
0004 #include <TChain.h>
0005 #include <TMath.h>
0006 #include <TTree.h>
0007 #include <TFile.h>
0008 #include <sstream> //std::ostringstsream
0009 #include <fstream> //std::ifstream
0010 #include <iostream> //std::cout, std::endl
0011 #include <cmath>
0012 #include <TGraphErrors.h>
0013 #include <TGraph.h>
0014 #include <TSpectrum.h>
0015 #include <TF1.h>
0016 #include <TTreeReader.h>
0017 #include <TTreeReaderValue.h>
0018 #include <TTreeReaderArray.h>
0019 #include <string>
0020 #include <set>
0021
0022 using namespace std;
0023
0024 void TrackCaloDist()
0025 {
0026     TFile *out = new TFile("TrackCaloDist.root","RECREATE");
0027
0028     TH1D* h_dist[6]; // distance between track (eta,phi) and shower vertex (eta,phi)
0029     TH1D* h_emcal[6]; // ratio of energy deposition in emcal to track momentum
0030     TH1D* h_ihcal[6]; // ratio of energy deposition in ihcal to track momentum
0031     TH1D* h_ohcal[6]; // ratio of energy deposition in ohcal to track momentum
0032     TH2D* h_totale_dist[6]; // total energy deposition from the shower vs distance between track and shower
0033
0034     // only filled in the case of an ohcal shower
0035     TH1D* ohcal_top = new TH1D("ohcal_top","",300,0,3); // ratio of energy in leading tower of ohcal to track momentum
0036     TH1D* ohcal_towers = new TH1D("ohcal_towers","",300,0,3); // ratio of energy in ohcal tower region to track mom
0037     TH1D* ohcal_total = new TH1D("ohcal_total","",300,0,3); // ratio of energy in all calo tower region to track mom
0038     TH1D* ohcal_ihcal_mip = new TH1D("ohcal_ihcal_mip","",300,0,1); // energy in ihcal in case of ohcal shower
0039     TH1D* ohcal_emcal_mip = new TH1D("ohcal_emcal_mip","",300,0,1); // energy in emcal in case of ohcal shower
0040     TH2D* h_4towers = new TH2D("h_4towers","",5,-2.5,2.5,5,-2.5,2.5); // 5x5 tower energy deposition in ohcal
0041
0042     // classifications of shower location
0043     string detector[6] = { "mip", "ohcal", "magnet", "ihcal", "emcal", "tracker"};
0044     for (int i = 0; i < 6; i++) {
0045         h_dist[i] = new TH1D(TString::Format("h_dist_%s",detector[i].c_str()),"",200,0,5);
0046         h_emcal[i] = new TH1D(TString::Format("h_emcale_%s",detector[i].c_str()),"",100,0,3);
0047         h_ihcal[i] = new TH1D(TString::Format("h_ihcale_%s",detector[i].c_str()),"",100,0,3);
0048         h_ohcal[i] = new TH1D(TString::Format("h_ohcale_%s",detector[i].c_str()),"",100,0,3);
0049         h_totale_dist[i] = new TH2D(TString::Format("h_totale_dist_%s",detector[i].c_str()),"",100,0,3,200,0,5);
0050     }
0051
0052     TH2D* h_ohcal_dist = new TH2D("h_ohcale_dist","",100,0,3,200,0,5);
0053     TH1D* h_class = new TH1D("h_class","",7,-0.5,6.5);
0054
0055     // distance between the track (eta,phi) and the shower vertex for each classification
0056     // energy distribution for each calorimeter for each classification
0057
0058     TFile *file = new TFile("TruthCaloTree.root");
0059     TTreeReader reader("T",file);
0060
0061     TTreeReaderValue<int> n(reader, "n_truth");
0062     TTreeReaderArray<float> E(reader, "truthenergy");
0063     TTreeReaderArray<int> pid(reader, "truthpid");
0064     TTreeReaderArray<float> eta(reader, "trutheta");
0065     TTreeReaderArray<float> phi(reader, "truthphi");
0066     TTreeReaderArray<float> p(reader, "truthp");
0067     TTreeReaderArray<float> pt(reader, "truthpt");
0068     TTreeReaderArray<int> track_id(reader, "truth_track_id");
0069
0070     TTreeReaderValue<int> n_child(reader, "n_child");
0071     TTreeReaderArray<int> child_pid(reader, "child_pid");
0072     TTreeReaderArray<int> child_parent_id(reader, "child_parent_id");
0073     TTreeReaderArray<int> child_vertex_id(reader, "child_vertex_id");
0074
0075     TTreeReaderValue<int> BH_n(reader, "BH_n");
0076     TTreeReaderArray<int> BH_track_id(reader, "BH_track_id");
0077     TTreeReaderArray<float> BH_e(reader, "BH_e");
0078
0079     TTreeReaderValue<int> n_vertex(reader, "n_vertex");
0080     TTreeReaderArray<int> vertex_id(reader, "vertex_id");
0081     TTreeReaderArray<float> x(reader, "vertex_x");
0082     TTreeReaderArray<float> y(reader, "vertex_y");
0083     TTreeReaderArray<float> z(reader, "vertex_z");
0084
0085     TTreeReaderArray<float> tower_sim_E(reader, "tower_sim_E");
0086     TTreeReaderArray<int> tower_sim_ieta(reader, "tower_sim_ieta");
0087     TTreeReaderArray<int> tower_sim_iphi(reader, "tower_sim_iphi");
0088     TTreeReaderArray<float> tower_sim_eta(reader, "tower_sim_eta");
0089     TTreeReaderArray<float> tower_sim_phi(reader, "tower_sim_phi");
0090     TTreeReaderValue<int> tower_sim_n(reader, "tower_sim_n");
0091
0092     TTreeReaderArray<float> EMcal_sim_E(reader, "EMcal_sim_E");
0093     TTreeReaderArray<int> EMcal_sim_ieta(reader, "EMcal_sim_ieta");
0094     TTreeReaderArray<int> EMcal_sim_iphi(reader, "EMcal_sim_iphi");
0095     TTreeReaderArray<float> EMcal_sim_eta(reader, "EMcal_sim_eta");
0096     TTreeReaderArray<float> EMcal_sim_phi(reader, "EMcal_sim_phi");
0097     TTreeReaderValue<int> EMcal_sim_n(reader, "EMcal_sim_n");
0098
0099     TTreeReaderArray<float> hcalIN_sim_E(reader, "hcalIN_sim_E");
0100     TTreeReaderArray<int> hcalIN_sim_ieta(reader, "hcalIN_sim_ieta");
0101     TTreeReaderArray<int> hcalIN_sim_iphi(reader, "hcalIN_sim_iphi");
0102     TTreeReaderArray<float> hcalIN_sim_eta(reader, "hcalIN_sim_eta");
0103     TTreeReaderArray<float> hcalIN_sim_phi(reader, "hcalIN_sim_phi");
0104     TTreeReaderValue<int> hcalIN_sim_n(reader, "hcalIN_sim_n");
0105
0106     const int nphi_hcal = 64; // number of hcal towers in phi direction
0107     const int neta_hcal = 24;
0108     const int nphi_emcal = 256;
0109     const int neta_emcal = 96;
0110     const float max_eta = 0.8; // maximum eta value for isolated track
0111     const float neigh_max_eta = 1.0;
0112     const float track_pt_cut = 0.5; // pt cut for neighboring particles
0113     const float high_pt_cut = 4.0;
0114     const float de = 0.045833; // distance from center of one calorimeter tower to edge in eta
0115     const float dp = 0.049087; // distacne from center of one calorimeter tower to edge in phi
0116     const float ihcal_sf = 0.162166; // sampling fractions
0117     const float ohcal_sf = 0.0338021;
0118     const float emcal_sf = 0.02;
0119     const float emcal_inner = 92; // inner radius of emcal
0120     const float emcal_outer = 116;
0121     const float ihcal_inner = 117.27;
0122     const float ihcal_outer = 134.42;
0123     const float ohcal_inner = 183.3;
0124     const float ohcal_outer = 274.317;
0125
0126     double dr, delta_r, delta_eta, delta_phi, dr_charge, dphi, deta, tower_dr;
0127     float mipEM, mipIH, mipOH;
0128     double vr, v_phi, v_eta;
0129     float avgtracks[10] = {0.0};
0130     int ntracks_centrality[10] = {0};
0131     int centrality_bins = 10;
0132     set<int> charged_pid = {-3334,-3312,-3222,-3112,-2212,-431,-411,-321,-211,-13,-11,11,
0133                          13,211,321,411,431,2212,3112,3222,3312,3334};
0134     set<int> neutral_pid = {-3322,-3212,-3122,-2112,-421,-14,-12,12,14,22,111,130,310,421,
0135                          2112,3122,3212,3322};
0136     int cent[] = {40,50,60,70,80};
0137
0138     float five_by_five[5][5] = {{0.0}};
0139     int child;
0140     set<int> vertex;
0141     float v_radius, v_rad;
0142     int classify;
0143
0144     float eta_map[24],eta_mapIH[24],eta_mapEM[96];
0145     float phi_map[64],phi_mapIH[64],phi_mapEM[256];
0146
0147
0148     // first find the mapping of (ieta,iphi) calorimeter towers to the actual (eta,phi) coordinates of the towers
0149     while (reader.Next()) {
0150         for (int i = 0; i < *tower_sim_n; i++) {
0151             if (tower_sim_ieta[i] >= 0 && tower_sim_ieta[i] < neta_hcal
0152                 && tower_sim_iphi[i] >= 0 && tower_sim_iphi[i] < nphi_hcal) {
0153                 phi_map[tower_sim_iphi[i]] = tower_sim_phi[i];
0154                 eta_map[tower_sim_ieta[i]] = tower_sim_eta[i];
0155             }
0156         }
0157
0158         for (int i = 0; i < *hcalIN_sim_n; i++) {
0159             if (hcalIN_sim_ieta[i] >= 0 && hcalIN_sim_ieta[i] < neta_hcal
0160                 && hcalIN_sim_iphi[i] >= 0 && hcalIN_sim_iphi[i] < nphi_hcal) {
0161                 phi_mapIH[hcalIN_sim_iphi[i]] = hcalIN_sim_phi[i];
0162                 eta_mapIH[hcalIN_sim_ieta[i]] = hcalIN_sim_eta[i];
0163             }
0164         }
0165
0166         for (int i = 0; i < nphi_hcal; i++) {
0167             if (phi_map[i] < 0.0) phi_map[i] += 2*M_PI;
0168             if (phi_mapIH[i] < 0.0) phi_mapIH[i] += 2*M_PI;
0169         }
0170
0171         for (int i = 0; i < *EMcal_sim_n; i++) {
0172             if (EMcal_sim_ieta[i] >= 0 && EMcal_sim_ieta[i] < neta_emcal
0173                 && EMcal_sim_iphi[i] >= 0 && EMcal_sim_iphi[i] < nphi_emcal) {
0174                 phi_mapEM[EMcal_sim_iphi[i]] = EMcal_sim_phi[i];
0175                 eta_mapEM[EMcal_sim_ieta[i]] = EMcal_sim_eta[i];
0176             }
0177         }
0178
0179         for (int i = 0; i < nphi_emcal; i++) {
0180             if (phi_mapEM[i] < 0.0) phi_mapEM[i] += 2*M_PI;
0181         }
0182     }
0183
0184     reader.Restart();
0185     int f = 0;
0186     while (reader.Next()) {
0187
0188         if ( f % 1000 == 0) cout << f << endl;
0189         if (!vertex.empty()) vertex.clear();
0190         classify = -1;
0191
0192         float e_sim[24][64] = {{0.0}};
0193         float e_simIH[24][64] = {{0.0}};
0194         float e_simEM[96][256] = {{0.0}};
0195         float total_energy = 0.0;
0196
0197         // find the energy distribution in the calo towers for the event
0198         for (int i = 0; i < *tower_sim_n; i++) {
0199             if (tower_sim_ieta[i] >= 0 && tower_sim_ieta[i] < neta_hcal
0200                 && tower_sim_iphi[i] >= 0 && tower_sim_iphi[i] < nphi_hcal) {
0201                 e_sim[tower_sim_ieta[i]][tower_sim_iphi[i]] += tower_sim_E[i]/ohcal_sf;
0202             }
0203         }
0204         for (int i = 0; i < *hcalIN_sim_n; i++) {
0205             if (hcalIN_sim_ieta[i] >= 0 && hcalIN_sim_ieta[i] < neta_hcal
0206                 && hcalIN_sim_iphi[i] >= 0 && hcalIN_sim_iphi[i] < nphi_hcal) {
0207                 e_simIH[hcalIN_sim_ieta[i]][hcalIN_sim_iphi[i]] += hcalIN_sim_E[i]/ihcal_sf;
0208             }
0209         }
0210         for (int i = 0; i < *EMcal_sim_n; i++) {
0211             if (EMcal_sim_ieta[i] >= 0 && EMcal_sim_ieta[i] < neta_emcal
0212                 && EMcal_sim_iphi[i] >= 0 && EMcal_sim_iphi[i] < nphi_emcal) {
0213                 e_simEM[EMcal_sim_ieta[i]][EMcal_sim_iphi[i]] += EMcal_sim_E[i]/emcal_sf;
0214             }
0215         }
0216
0217         // loop over truth level particles
0218         for (int i = 0; i < *n; i++) {
0219
0220             // find charged isolated high pt particle and check that it is isolated
0221             if (pt[i] >= high_pt_cut && charged_pid.find(pid[i])!=charged_pid.end() && abs(double(eta[i])) < max_eta) {
0222                 if (phi[i] < 0.0) phi[i] += 2*M_PI;
0223                 dr = 100.0;
0224                 for (int j = 0; j < *n; j++) {
0225                     if (i != j && pt[j] > track_pt_cut && abs(double(eta[j])) < neigh_max_eta) {
0226                         if (phi[j] < 0.0) phi[j] += 2*M_PI;
0227                         delta_phi = double(phi[i]) - double(phi[j]);
0228                         if (delta_phi > M_PI) delta_phi -= M_PI;
0229                         delta_eta = double(eta[i]) - double(eta[j]);
0230                         delta_r = sqrt(delta_phi*delta_phi + delta_eta*delta_eta);
0231                         if (delta_r < dr) dr = delta_r;
0232                     }
0233                 }
0234
0235                 // check where particle begins showering
0236                 if (dr >= 0.2) {
0237                     for (int j = 0; j < *n_child; j++) {
0238                         if (child_parent_id[j] == track_id[i] && abs(child_pid[j]) != 11 && child_pid[j] != 22) vertex.insert(child_vertex_id[j]);
0239                     }
0240                     v_radius = ohcal_outer + 1;
0241                     for (set<int>::iterator it = vertex.begin(); it != vertex.end(); it++) {
0242                         child = 0;
0243                         for (int j = 0 ; j < *n_child; j++) {
0244                             if (child_vertex_id[j] == *it && abs(child_pid[j]) != 11 && child_pid[j] != 22) {
0245                                 child++;
0246                             }
0247                         }
0248                         if (child > 1) {
0249                             // now find the location of the shower vertex
0250                             for (int v = 0; v < *n_vertex; v++) {
0251                                 if (vertex_id[v] == *it) {
0252                                     v_rad = sqrt(x[v]*x[v] + y[v]*y[v]);
0253                                     if (v_rad < v_radius) {
0254                                         v_radius = v_rad;
0255                                         v_phi = atan2(y[v],x[v]);
0256                                         v_eta = atanh(z[v] / sqrt(x[v]*x[v] + y[v]*y[v] + z[v]*z[v]));
0257                                         dphi = abs(double(v_eta - eta[i]));
0258                                         deta = abs(double(v_phi - phi[i]));
0259                                         if (dphi >= M_PI) dphi -= M_PI;
0260                                         vr = sqrt(dphi*dphi + deta*deta);
0261                                     }
0262                                 }
0263                             }
0264                         }
0265                     }
0266
0267
0268                     // find the energy deposition and number of towers in tower region
0269                     int towers[200][2];
0270                     float top_energy = 0.0;
0271                     int ti, tj;
0272                     int numOH = 0;
0273                     int tn = 0;
0274                     mipEM = 0.0;
0275                     mipIH = 0.0;
0276                     mipOH = 0.0;
0277                     for (int j = 0; j < neta_hcal; j++) {
0278                         for (int k = 0; k < nphi_hcal; k++) {
0279                             dphi = abs(double(eta_map[j] - eta[i]));
0280                             deta = abs(double(phi_map[k] - phi[i]));
0281                             if (dphi >= M_PI) dphi -= M_PI;
0282                             tower_dr = sqrt(dphi*dphi + deta*deta);
0283                             // find towers in (eta,phi) range of track
0284                             if (tower_dr <= 0.2) {
0285                                 towers[tn][0] = j;
0286                                 towers[tn][1] = k;
0287                                 mipIH += e_simIH[j][k];
0288                                 mipOH += e_sim[j][k];
0289                                 tn++;
0290                             }
0291                         }
0292                     }
0293
0294                     for (int j = 0; j < neta_emcal; j++) {
0295                         for(int k = 0; k < nphi_emcal; k++) {
0296                             dphi = abs(double(eta_mapEM[j] - eta[i]));
0297                             deta = abs(double(phi_mapEM[k] - phi[i]));
0298                             if (dphi >= M_PI) dphi -= M_PI;
0299                             tower_dr = sqrt(dphi*dphi + deta*deta);
0300                             if (tower_dr <= 0.2) {
0301                                 mipEM += e_simEM[j][k];
0302                             }
0303                         }
0304                     }
0305
0306                     int j;
0307                     int k;
0308                     double mipTotal = mipEM + mipIH + mipOH;
0309
0310                     for (int t = 0; t < tn; t++) {
0311                         j = towers[t][0];
0312                         k = towers[t][1];
0313                         if (e_sim[j][k] > top_energy) {
0314                             top_energy = e_sim[j][k];
0315                             ti = j;
0316                             tj = k;
0317                         }
0318                     }
0319
0320                     // based on shower location classification fill histograms
0321                     if (v_radius < emcal_inner) {
0322                         h_emcal[5]->Fill(mipEM/p[i]);
0323                         h_ihcal[5]->Fill(mipIH/p[i]);
0324                         h_ohcal[5]->Fill(mipOH/p[i]);
0325                         h_dist[5]->Fill(vr);
0326                         h_class->Fill(5);
0327                         h_totale_dist[5]->Fill(mipTotal/p[i],vr);
0328                     } else if (v_radius >= emcal_inner && v_radius <= emcal_outer) {
0329                         h_emcal[4]->Fill(mipEM/p[i]);
0330                         h_ihcal[4]->Fill(mipIH/p[i]);
0331                         h_ohcal[4]->Fill(mipOH/p[i]);
0332                         h_dist[4]->Fill(vr);
0333                         h_class->Fill(4);
0334                         h_totale_dist[4]->Fill(mipTotal/p[i],vr);
0335                     } else if (v_radius >= ihcal_inner && v_radius <= ihcal_outer) {
0336                         h_emcal[3]->Fill(mipEM/p[i]);
0337                         h_ihcal[3]->Fill(mipIH/p[i]);
0338                         h_ohcal[3]->Fill(mipOH/p[i]);
0339                         h_dist[3]->Fill(vr);
0340                         h_class->Fill(3);
0341                         h_totale_dist[3]->Fill(mipTotal/p[i],vr);
0342                     } else if (v_radius > ihcal_outer && v_radius < ohcal_inner) {
0343                         h_emcal[2]->Fill(mipEM/p[i]);
0344                         h_ihcal[2]->Fill(mipIH/p[i]);
0345                         h_ohcal[2]->Fill(mipOH/p[i]);
0346                         h_dist[2]->Fill(vr);
0347                         h_class->Fill(2);
0348                         h_totale_dist[2]->Fill(mipTotal/p[i],vr);
0349                     } else if (v_radius >= ohcal_inner && v_radius <= ohcal_outer) {
0350                         h_emcal[1]->Fill(mipEM/p[i]);
0351                         h_ihcal[1]->Fill(mipIH/p[i]);
0352                         h_ohcal[1]->Fill(mipOH/p[i]);
0353                         h_dist[1]->Fill(vr);
0354                         h_class->Fill(1);
0355                         h_ohcal_dist->Fill(mipOH/p[i],vr);
0356                         h_totale_dist[1]->Fill(mipTotal/p[i],vr);
0357                         ohcal_top->Fill(top_energy/p[i]);
0358                         ohcal_towers->Fill(mipOH/p[i]);
0359                         ohcal_total->Fill(mipTotal/p[i]);
0360                         ohcal_ihcal_mip->Fill(mipIH);
0361                         ohcal_emcal_mip->Fill(mipEM);
0362
0363                         // fill 5x5 tower energy distribution for ohcal showers
0364                         for (int fi = -2; fi < 3; fi++) {
0365                             for (int fj = -2; fj < 3; fj++) {
0366                                 five_by_five[fi+2][fj+2] += e_sim[ti+fi][tj+fj];
0367                                 h_4towers->SetBinContent(h_4towers->FindBin(fi,fj),five_by_five4[fi+2][fj+2]);
0368                             }
0369                         }
0370                     } else if (v_radius > ohcal_outer) {
0371                         // check if the event mips through the entire system
0372                         for (int j = 0; j < *BH_n; j++) {
0373                             if (BH_track_id[j] == track_id[i]) {
0374                                 h_emcal[0]->Fill(mipEM/p[i]);
0375                                 h_ihcal[0]->Fill(mipIH/p[i]);
0376                                 h_ohcal[0]->Fill(mipOH/p[i]);
0377                                 h_class->Fill(0);
0378                                 h_totale_dist[0]->Fill(mipTotal/p[i],vr);
0379                             }
0380                         }
0381                     }
0382                 }
0383             }
0384         }
0385         f++;
0386     }
0387     out->Write();
0388     out->Close();
0389 }