DST_MC/DeltaR_check/CheckDeltaR.C

0001 #include "CheckDeltaR.h"
0002 #include <iostream>
0003 using namespace std;
0004
0005 CheckDeltaR::CheckDeltaR(string run_type, string input_folder_directory, string input_file_name, string input_tree_name, string out_folder_directory)
0006 :
0007 run_type(run_type),
0008 input_folder_directory(input_folder_directory),
0009 input_file_name(input_file_name),
0010 input_tree_name(input_tree_name),
0011 out_folder_directory(out_folder_directory)
0012 {
0013     track_pair.clear();
0014     track_pair_deltaR.clear();
0015     true_track_used_tag.clear();
0016     reco_track_used_tag.clear();
0017
0018     plot_text = (run_type == "MC") ? "Simulation" : "Work-in-progress";
0019
0020     cout << "CheckDeltaR::CheckDeltaR()" << endl;
0021     file_in = TFile::Open(Form("%s/%s", input_folder_directory.c_str(), input_file_name.c_str()));
0022     tree_in = (TTree *)file_in->Get(input_tree_name.c_str());
0023     readeta = new ReadEta_v3(tree_in);
0024
0025     cout<<"In CheckDeltaR::CheckDeltaR(), tree entry : "<<tree_in -> GetEntries()<<endl;
0026
0027     system(Form("mkdir %s", out_folder_directory.c_str()));
0028
0029     Init();
0030 }
0031
0032 void CheckDeltaR::Init()
0033 {
0034     c1 = new TCanvas("","",950,800);
0035     c1 -> cd();
0036
0037     ltx = new TLatex();
0038     ltx->SetNDC();
0039     ltx->SetTextSize(0.045);
0040     ltx->SetTextAlign(31);
0041
0042     DeltaPhi_All_1D = new TH1F("","DeltaPhi_All_1D",  100, -5, 5);
0043     DeltaPhi_All_1D -> GetXaxis() -> SetTitle("#Delta#phi [degree]");
0044     DeltaPhi_All_1D -> GetYaxis() -> SetTitle("Entry");
0045     DeltaPhi_All_1D -> GetXaxis() -> SetNdivisions(505);
0046
0047     DeltaEta_All_1D = new TH1F("","DeltaEta_All_1D", 100, -5, 5);
0048     DeltaEta_All_1D -> GetXaxis() -> SetTitle("#Delta#eta");
0049     DeltaEta_All_1D -> GetYaxis() -> SetTitle("Entry");
0050     DeltaEta_All_1D -> GetXaxis() -> SetNdivisions(505);
0051
0052     DeltaPhiDeltaEta_All_2D = new TH2F("","DeltaPhiDeltaEta_All_2D",  100, -5, 5, 100, -5, 5);
0053     DeltaPhiDeltaEta_All_2D -> GetXaxis() -> SetTitle("#Delta#phi [degree]");
0054     DeltaPhiDeltaEta_All_2D -> GetYaxis() -> SetTitle("#Delta#eta");
0055     DeltaPhiDeltaEta_All_2D -> GetXaxis() -> SetNdivisions(505);
0056
0057
0058     DeltaPhi_link_1D = new TH1F("","DeltaPhi_link_1D",  100, -5, 5);
0059     DeltaPhi_link_1D -> GetXaxis() -> SetTitle("#Delta#phi [degree]");
0060     DeltaPhi_link_1D -> GetYaxis() -> SetTitle("Entry");
0061     DeltaPhi_link_1D -> GetXaxis() -> SetNdivisions(505);
0062
0063     DeltaEta_link_1D = new TH1F("","DeltaEta_link_1D", 100, -0.15, 0.15);
0064     DeltaEta_link_1D -> GetXaxis() -> SetTitle("#Delta#eta");
0065     DeltaEta_link_1D -> GetYaxis() -> SetTitle("Entry");
0066     DeltaEta_link_1D -> GetXaxis() -> SetNdivisions(505);
0067
0068     DeltaPhiDeltaEta_link_2D = new TH2F("","DeltaPhiDeltaEta_link_2D",  100, -5, 5, 100, -5, 5);
0069     DeltaPhiDeltaEta_link_2D -> GetXaxis() -> SetTitle("#Delta#phi [degree]");
0070     DeltaPhiDeltaEta_link_2D -> GetYaxis() -> SetTitle("#Delta#eta");
0071     DeltaPhiDeltaEta_link_2D -> GetXaxis() -> SetNdivisions(505);
0072
0073     DeltaR_best_1D = new TH1F("","DeltaR_best_1D", 100, -10, 10);
0074     DeltaR_best_1D -> GetXaxis() -> SetTitle("#DeltaR, with weighted #Delta#eta");
0075     DeltaR_best_1D -> GetYaxis() -> SetTitle("Entry");
0076     DeltaR_best_1D -> GetXaxis() -> SetNdivisions(505);
0077
0078     DeltaPhiDeltaEta_best_2D = new TH2F("","DeltaPhiDeltaEta_best_2D",  100, -5, 5, 100, -5, 5);
0079     DeltaPhiDeltaEta_best_2D -> GetXaxis() -> SetTitle("#Delta#phi [degree]");
0080     DeltaPhiDeltaEta_best_2D -> GetYaxis() -> SetTitle("Weighted #Delta#eta");
0081     DeltaPhiDeltaEta_best_2D -> GetXaxis() -> SetNdivisions(505);
0082
0083     DeltaEtaTrueEta_best_2D = new TH2F("","DeltaEtaTrueEta_best_2D",  100, -3, 3, 100, -0.15, 0.15);
0084     DeltaEtaTrueEta_best_2D -> GetXaxis() -> SetTitle("True track #eta");
0085     DeltaEtaTrueEta_best_2D -> GetYaxis() -> SetTitle("#Delta#eta [reco - true]");
0086     DeltaEtaTrueEta_best_2D -> GetXaxis() -> SetNdivisions(505);
0087
0088     DeltaEtaRecoZ_best_2D = new TH2F("","DeltaEtaRecoZ_best_2D",  100, -400, 100, 100, -0.15, 0.15);
0089     DeltaEtaRecoZ_best_2D -> GetXaxis() -> SetTitle("Reco zvtx [mm]");
0090     DeltaEtaRecoZ_best_2D -> GetYaxis() -> SetTitle("#Delta#eta [reco - true]");
0091     DeltaEtaRecoZ_best_2D -> GetXaxis() -> SetNdivisions(505);
0092
0093     DeltaEtaTrueZ_best_2D = new TH2F("","DeltaEtaTrueZ_best_2D",  100, -400, 100, 100, -0.15, 0.15);
0094     DeltaEtaTrueZ_best_2D -> GetXaxis() -> SetTitle("True zvtx [mm]");
0095     DeltaEtaTrueZ_best_2D -> GetYaxis() -> SetTitle("#Delta#eta [reco - true]");
0096     DeltaEtaTrueZ_best_2D -> GetXaxis() -> SetNdivisions(505);
0097
0098 }
0099
0100 void CheckDeltaR::Processing(int run_Nevent)
0101 {
0102     cout << "CheckDeltaR::Processing()" << endl;
0103
0104     double best_delta_phi;
0105     double best_delta_eta;
0106     double best_delta_R;
0107
0108     int final_run_Nevent = (run_Nevent == -1 || run_Nevent > tree_in -> GetEntries()) ? tree_in -> GetEntries() : run_Nevent;
0109
0110     for (int i = 0; i < final_run_Nevent/*tree_in -> GetEntries()*/; i++)
0111     {
0112         readeta -> LoadTree(i);
0113         readeta -> GetEntry(i);
0114         if (i % 1 == 0) cout << "Processing entry " << i << endl;
0115
0116         if (readeta -> NTrueTrack > 500 || readeta -> NTrueTrack < 100) {continue;}
0117         if (readeta -> True_zvtx > -180 || readeta -> True_zvtx < -220) {continue;}
0118
0119         for (int true_i = 0; true_i < readeta -> NTrueTrack; true_i++)
0120         {
0121             for (int reco_i = 0; reco_i < readeta -> N2Clu_track; reco_i++)
0122             {
0123                 double delta_phi = readeta -> RecoTrack_phi_d -> at(reco_i) - readeta -> TrueTrack_phi_d -> at(true_i);
0124                 double delta_eta = readeta -> RecoTrack_eta_d -> at(reco_i) - readeta -> TrueTrack_eta_d -> at(true_i);
0125                 DeltaPhi_All_1D -> Fill(delta_phi);
0126                 DeltaEta_All_1D -> Fill(delta_eta);
0127                 DeltaPhiDeltaEta_All_2D -> Fill(delta_phi, delta_eta);
0128
0129                 if (reco_i == 0)
0130                 {
0131                     best_delta_phi = delta_phi;
0132                     best_delta_eta = delta_eta;
0133                 }
0134                 else
0135                 {
0136                     if (abs(delta_phi) < abs(best_delta_phi))
0137                     {
0138                         best_delta_phi = delta_phi;
0139                     }
0140                     if (abs(delta_eta) < abs(best_delta_eta))
0141                     {
0142                         best_delta_eta = delta_eta;
0143                     }
0144
0145                 }
0146
0147                 track_pair.push_back({double(true_i), double(reco_i), delta_phi, delta_eta, readeta -> TrueTrack_eta_d -> at(true_i)});
0148                 track_pair_deltaR.push_back(sqrt(pow(delta_phi,2) + pow(delta_eta * eta_weight,2)));
0149             }
0150
0151             DeltaPhi_link_1D -> Fill(best_delta_phi);
0152             DeltaEta_link_1D -> Fill(best_delta_eta);
0153         }
0154
0155         long long vec_size = track_pair_deltaR.size();
0156         long long ind[track_pair_deltaR.size()];
0157         TMath::Sort(vec_size, &track_pair_deltaR[0], ind, false);
0158
0159         for (int index = 0; index < track_pair_deltaR.size(); index++)
0160         {
0161             int better_index = ind[index];
0162             if (true_track_used_tag[int(track_pair[better_index][0])] != 0) {continue;}
0163             if (reco_track_used_tag[int(track_pair[better_index][1])] != 0) {continue;}
0164
0165             DeltaR_best_1D -> Fill(track_pair_deltaR[better_index]);
0166             DeltaPhiDeltaEta_best_2D -> Fill(track_pair[better_index][2], track_pair[better_index][3]*eta_weight);
0167
0168             if (track_pair_deltaR[better_index] < 2)
0169             {
0170                 DeltaEtaTrueEta_best_2D -> Fill(track_pair[better_index][4], track_pair[better_index][3]);
0171                 DeltaEtaRecoZ_best_2D -> Fill(readeta -> Evt_zvtx, track_pair[better_index][3]);
0172                 DeltaEtaTrueZ_best_2D -> Fill(readeta -> True_zvtx, track_pair[better_index][3]);
0173             }
0174
0175             true_track_used_tag[int(track_pair[better_index][0])] += 1;
0176             reco_track_used_tag[int(track_pair[better_index][1])] += 1;
0177         }
0178
0179         ClearEvt();
0180     } // note : event loop
0181
0182     return;
0183 }
0184
0185 void CheckDeltaR::ClearEvt()
0186 {
0187     track_pair.clear();
0188     track_pair_deltaR.clear();
0189     true_track_used_tag.clear();
0190     reco_track_used_tag.clear();
0191 }
0192
0193 void CheckDeltaR::PrintPlots()
0194 {
0195     c1 -> cd();
0196     DeltaPhi_All_1D -> Draw("hist");
0197     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0198     c1 -> Print(Form("%s/DeltaPhi_All_1D.pdf",out_folder_directory.c_str()));
0199     c1 -> Clear();
0200
0201     c1 -> cd();
0202     DeltaEta_All_1D -> Draw("hist");
0203     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0204     c1 -> Print(Form("%s/DeltaEta_All_1D.pdf",out_folder_directory.c_str()));
0205     c1 -> Clear();
0206
0207     c1 -> cd();
0208     DeltaPhiDeltaEta_All_2D -> Draw("colz0");
0209     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0210     c1 -> Print(Form("%s/DeltaPhiDeltaEta_All_2D.pdf",out_folder_directory.c_str()));
0211     c1 -> Clear();
0212
0213     c1 -> cd();
0214     DeltaPhi_link_1D -> Draw("hist");
0215     cout<<"DeltaPhi_link_1D, stddev : "<<DeltaPhi_link_1D -> GetStdDev()<<endl;
0216     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0217     c1 -> Print(Form("%s/DeltaPhi_link_1D.pdf",out_folder_directory.c_str()));
0218     c1 -> Clear();
0219
0220     c1 -> cd();
0221     DeltaEta_link_1D -> Draw("hist");
0222     cout<<"DeltaEta_link_1D, stddev : "<<DeltaEta_link_1D -> GetStdDev()<<endl;
0223     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0224     c1 -> Print(Form("%s/DeltaEta_link_1D.pdf",out_folder_directory.c_str()));
0225     c1 -> Clear();
0226
0227     c1 -> cd();
0228     DeltaR_best_1D -> Draw("hist");
0229     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0230     c1 -> Print(Form("%s/DeltaR_best_1D.pdf",out_folder_directory.c_str()));
0231     c1 -> Clear();
0232
0233     c1 -> cd();
0234     DeltaPhiDeltaEta_best_2D -> Draw("colz0");
0235     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0236     c1 -> Print(Form("%s/DeltaPhiDeltaEta_best_2D.pdf",out_folder_directory.c_str()));
0237     c1 -> Clear();
0238
0239     DeltaEtaTrueEta_best_2D -> Draw("colz0");
0240     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0241     c1 -> Print(Form("%s/DeltaEtaTrueEta_best_2D.pdf",out_folder_directory.c_str()));
0242     c1 -> Clear();
0243
0244     DeltaEtaRecoZ_best_2D -> Draw("colz0");
0245     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0246     c1 -> Print(Form("%s/DeltaEtaRecoZ_best_2D.pdf",out_folder_directory.c_str()));
0247     c1 -> Clear();
0248
0249     DeltaEtaTrueZ_best_2D -> Draw("colz0");
0250     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX INTT}} %s", plot_text.c_str()));
0251     c1 -> Print(Form("%s/DeltaEtaTrueZ_best_2D.pdf",out_folder_directory.c_str()));
0252     c1 -> Clear();
0253
0254     return;
0255 }