//Jordan Seneca //Usage: PlotRecoTrks filename //#include "/sps/km3net/users/jseneca/docopt.cpp/docopt.h" #include "PlotRecoTrks.hh" #include "Det.hh" using namespace std; int main(int argc, char *argv[] ){ cout << "Input file: " << argv[1] << endl ; string DETECTOR_FILE = "/pbs/throng/km3net/detectors/km3net_jul13_90m.det"; const uint rebin_size_E = 100; const uint rebin_size_nhits = 10; const uint region_n = 3; const Double_t logE_bins[11] = { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 }; TFile* file = TFile::Open(argv[1], "READ"); TFile infile(argv[1], "READ", "infile"); TH3D* h3dAng[region_n]; TH3D* h4dAng[region_n]; TH3D* h4dAngOld[region_n]; TH3D* h3dErec[region_n]; TH3D* h4dErec[region_n]; TH3D* h4dErecOld[region_n]; TH2D* h3dE[region_n]; TH2D* h4dE[region_n]; TH2D* hMcE[region_n]; TH3D* hPos[region_n]; for( uint iregion = 0; iregion < region_n; iregion++ ){ cout << "loading: " << (TString)"h3dAngENhits" + to_string(iregion) << endl; h3dAng[iregion] = (TH3D*)infile.Get( (TString)"h3dAngENhits" + to_string(iregion) ); h4dAng[iregion] = (TH3D*)infile.Get( (TString)"h4dAngENhits" + to_string(iregion) ); h4dAngOld[iregion] = (TH3D*)infile.Get( (TString)"h4dAngOldENhits" + to_string(iregion) ); h3dErec[iregion] = (TH3D*)infile.Get( (TString)"h3dErecENhits" + to_string(iregion) ); h4dErec[iregion] = (TH3D*)infile.Get( (TString)"h4dErecENhits" + to_string(iregion) ); h4dErecOld[iregion] = (TH3D*)infile.Get( (TString)"h4dErecOldENhits" + to_string(iregion) ); h3dE[iregion] = (TH2D*)infile.Get( (TString)"h3dENhits" + to_string( iregion ) ); h4dE[iregion] = (TH2D*)infile.Get( (TString)"h4dENhits" + to_string( iregion ) ); hMcE[iregion] = (TH2D*)infile.Get( (TString)"hMcENhits" + to_string( iregion ) ); hPos[iregion] = (TH3D*)infile.Get( (TString)"hPosition" + to_string( iregion ) ); } TH2D* h3dAngDist = (TH2D*)infile.Get("h3dAngDist"); TH2D* h4dAngDist = (TH2D*)infile.Get("h4dAngDist"); TH2D* h3dErecDist = (TH2D*)infile.Get("h3dErecDist"); TH2D* h4dErecDist = (TH2D*)infile.Get("h4dErecDist"); TH2D* h3dStopwatch = (TH2D*)infile.Get("h3dStopwatch"); TH2D* h4dStopwatch = (TH2D*)infile.Get("h4dStopwatch"); const uint size_arr_E = h3dAng[0]->GetNbinsY()/rebin_size_E; const uint size_arr_nhits = h3dAng[0]->GetNbinsZ()/rebin_size_nhits; cout << "Size array: " << size_arr_E << endl; cout << "Size array: " << size_arr_nhits << endl; vector< vector< TH1D* > > h3d_dir_E ( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h4d_dir_E ( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h4d_dir_Eold( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h3d_dir_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< vector< TH1D* > > h4d_dir_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< TH1D* > h3d_dirbyE ( region_n, new TH1D( "3ddirbyE", "", 10, logE_bins ) ); vector< TH1D* > h4d_dirbyE ( region_n, new TH1D( "4ddirbyE", "", 10, logE_bins ) ); vector< TH1D* > h4d_dirbyEold( region_n, new TH1D( "4ddirbyEold", "", 10, logE_bins ) ); vector< TH1D* > h3d_dirbynhits( region_n, new TH1D( "3ddirbynhits", "", size_arr_nhits, 1000, 100000 ) ); vector< TH1D* > h4d_dirbynhits( region_n, new TH1D( "4ddirbynhits", "", size_arr_nhits, 1000, 100000 ) ); vector< vector< TH1D* > > h3d_Erec_E ( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h4d_Erec_E ( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h4d_Erec_Eold( region_n, vector< TH1D* >( size_arr_E ) ); vector< vector< TH1D* > > h3d_Erec_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< vector< TH1D* > > h4d_Erec_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< TH1D* > h3d_ErecbyE ( region_n, new TH1D( "3dErecbyE", "", 10, logE_bins ) ); vector< TH1D* > h4d_ErecbyE ( region_n, new TH1D( "4dErecbyE", "", 10, logE_bins ) ); vector< TH1D* > h4d_ErecbyEold( region_n, new TH1D( "4dErecbyEold", "", 10, logE_bins ) ); vector< TH1D* > h3d_Erecbynhits( region_n ); vector< TH1D* > h4d_Erecbynhits( region_n ); vector< vector< TH1D* > > h3d_E_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< vector< TH1D* > > h4d_E_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< vector< TH1D* > > hMc_E_nhits( region_n, vector< TH1D* >( size_arr_nhits ) ); vector< TH2D* > hPosXY( region_n ); vector< TH2D* > hPosXZ( region_n ); vector< TH2D* > hDetXY( region_n ); vector< TH2D* > hDetXZ( region_n ); TCanvas c1("c1", "Comparison", 800, 800); gStyle->SetOptStat(0); for( uint iregion = 0; iregion < region_n; iregion++ ){ cout << "iregion: " << iregion << endl; h3dE[iregion] = new TH2D( (TString)"h3dENhits" + to_string(iregion), "3d E; E [GeV]; hit_N", 1000, 0, 1e6, 100, 0, 1e4); //hPosXY[iregion] = new TH2D( (TString)"posxy" + to_string(iregion), "dom pos; X pos; Y pos", 1000, -500, 500, 1000, -500, 500 ); //hPosXZ[iregion] = new TH2D( (TString)"posxz" + to_string(iregion), "dom pos; X pos; Y pos", 1000, -500, 500, 1000, -500, 500 ); hDetXY[iregion] = new TH2D( (TString)"detector" + to_string(iregion), "dom pos; X pos; Y pos", 300, -600, 600, 300, -600, 600 ); hDetXZ[iregion] = new TH2D( (TString)"detector" + to_string(iregion), "dom pos; X pos; Z pos", 300, -600, 600, 400, 0, 800 ); for(uint ibin = 0; ibin != 10; ibin++ ){ cout << "ibin: " << ibin << endl; const uint low_bin = pow( 2, ibin ); const uint hi_bin = pow( 2, ibin + 1 ); h3d_dir_E[iregion][ibin] = (TH1D*)h3dAng[iregion]->ProjectionX( ( "projection_3dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); h4d_dir_E[iregion][ibin] = (TH1D*)h4dAng[iregion]->ProjectionX( ( "projection_4dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); const string dir_Ename = "dir_E_reg" + to_string( iregion ) + "_" + to_string( (uint)h3dAng[iregion]->GetYaxis()->GetBinLowEdge( hi_bin + 1 ) / 1000 ) + "TeV"; const vector< double > medians = plot_2poisson( &c1, dir_Ename, h3d_dir_E[iregion][ibin], h4d_dir_E[iregion][ibin], "#Delta angle [deg]", ""); h3d_dirbyE[iregion]->SetBinContent( ibin + 1, medians[0] ); h4d_dirbyE[iregion]->SetBinContent( ibin + 1, medians[1] ); h4d_dir_Eold[iregion][ibin] = (TH1D*)h4dAngOld[iregion]->ProjectionX( ( "projection_4dEold" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); const string dir_EnameOld = "dir_Eold_reg" + to_string( iregion ) + "_" + to_string( (uint)h4dAngOld[iregion]->GetYaxis()->GetBinLowEdge( hi_bin + 1 ) / 1000 ) + "TeV"; const vector< double > mediansOld = plot_2poisson( &c1, dir_EnameOld, h3d_dir_E[iregion][ibin], h4d_dir_Eold[iregion][ibin], "#Delta angle [deg]", ""); h4d_dirbyEold[iregion]->SetBinContent( ibin + 1, mediansOld[1] ); h3d_Erec_E[iregion][ibin] = (TH1D*)h3dErec[iregion]->ProjectionX( ( "projection_3dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); h4d_Erec_E[iregion][ibin] = (TH1D*)h4dErec[iregion]->ProjectionX( ( "projection_4dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); const string Erec_Ename = "Erec_E_reg" + to_string( iregion ) + "_" + to_string( (uint)h3dErec[iregion]->GetYaxis()->GetBinLowEdge( hi_bin ) / 1000 ) + "TeV"; const vector< double > gausParams = plot_2gauss( &c1, Erec_Ename, h3d_Erec_E[iregion][ibin], h4d_Erec_E[iregion][ibin], "#frac{E_rec - E_{#nu mc}}{E_true}", "", false, false, true); h3d_ErecbyE[iregion]->SetBinContent( ibin + 1, gausParams[0] ); h3d_ErecbyE[iregion]->SetBinError ( ibin + 1, gausParams[1] ); h4d_ErecbyE[iregion]->SetBinContent( ibin + 1, gausParams[2] ); h4d_ErecbyE[iregion]->SetBinError ( ibin + 1, gausParams[3] ); h4d_Erec_Eold[iregion][ibin] = (TH1D*)h4dErecOld[iregion]->ProjectionX( ( "projection_4dEold" + to_string( iregion ) + to_string( ibin ) ).c_str(), low_bin, hi_bin ); const string Erec_EnameOld = "Erec_Eold_reg" + to_string( iregion ) + "_" + to_string( (uint)h4dErecOld[iregion]->GetYaxis()->GetBinLowEdge( hi_bin ) / 1000 ) + "TeV"; plot_2gauss( &c1, Erec_EnameOld, h3d_Erec_E[iregion][ibin], h4d_Erec_Eold[iregion][ibin], "#frac{E_rec - E_{#nu mc}}{E_true}", "", false, false, true); h4d_ErecbyEold[iregion]->SetBinContent( ibin + 1, gausParams[2] ); h4d_ErecbyEold[iregion]->SetBinError ( ibin + 1, gausParams[3] ); } const string dir_Ename = "dirbyE_reg" + to_string( iregion ); plot_th1( &c1, dir_Ename, h3d_dirbyE[iregion], h4d_dirbyE[iregion], "E [PeV]", "median #Delta angle [deg]", true); const string dir_EnameOld = "dirbyEold_reg" + to_string( iregion ); plot_th1( &c1, dir_EnameOld, h3d_dirbyE[iregion], h4d_dirbyEold[iregion], "E [PeV]", "median #Delta angle [deg]", true); const string Erec_Ename = "ErecbyE_reg" + to_string( iregion ); plot_th1( &c1, Erec_Ename, h3d_ErecbyE[iregion], h4d_ErecbyE[iregion], "E [PeV]", "median #Delta angle [deg]", true, false, false); const string Erec_EnameOld = "ErecbyEold_reg" + to_string( iregion ); plot_th1( &c1, Erec_EnameOld, h3d_ErecbyE[iregion], h4d_ErecbyEold[iregion], "E [PeV]", "median #Delta angle [deg]", true, false, false); for(uint ibin = 0; ibin * rebin_size_nhits != h3dAng[iregion]->GetNbinsZ(); ibin ++ ){ cout << "ibin: " << ibin << endl; const uint bin = ibin * rebin_size_nhits; h3d_dir_nhits[iregion][ibin] = (TH1D*)h3dAng[iregion]->ProjectionX( ( "projection_3dNhits" + to_string( iregion ) + to_string( ibin ) ).c_str(), 0, -1, bin, bin + rebin_size_nhits ); h4d_dir_nhits[iregion][ibin] = (TH1D*)h4dAng[iregion]->ProjectionX( ( "projection_4dNhits" + to_string( iregion ) + to_string( ibin ) ).c_str(), 0, -1, bin, bin + rebin_size_nhits ); const string dir_Nhitsname = "dir_Nhits_reg" + to_string( iregion ) + "_nh" + to_string( (uint)h3dAng[iregion]->GetZaxis()->GetBinLowEdge( bin + rebin_size_nhits + 1 ) ); const vector< double > medians = plot_2poisson( &c1, dir_Nhitsname, h3d_dir_nhits[iregion][ibin], h4d_dir_nhits[iregion][ibin], "#Delta angle [deg]", ""); h3d_dirbynhits[iregion]->SetBinContent( ibin + 1, medians[0] ); h4d_dirbynhits[iregion]->SetBinContent( ibin + 1, medians[1] ); h3d_Erec_nhits[iregion][ibin] = (TH1D*)h3dErec[iregion]->ProjectionX( ( "projection_3dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), 0, -1, bin, bin + rebin_size_nhits ); h4d_Erec_nhits[iregion][ibin] = (TH1D*)h4dErec[iregion]->ProjectionX( ( "projection_4dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), 0, -1, bin, bin + rebin_size_nhits ); const string Erec_nhitsname = "Erec_nhits_reg" + to_string( iregion ) + "_nh" + to_string( (uint)h3dErec[iregion]->GetZaxis()->GetBinLowEdge( bin + rebin_size_nhits + 1 ) ); plot_2gauss( &c1, Erec_nhitsname, h3d_Erec_nhits[iregion][ibin], h4d_Erec_nhits[iregion][ibin], "E_rec", "nhits"); h3d_E_nhits[iregion][ibin] = (TH1D*)h3dE[iregion]->ProjectionY( ( "projection_3dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), bin, bin + rebin_size_nhits ); h4d_E_nhits[iregion][ibin] = (TH1D*)h4dE[iregion]->ProjectionY( ( "projection_4dE" + to_string( iregion ) + to_string( ibin ) ).c_str(), bin, bin + rebin_size_nhits ); hMc_E_nhits[iregion][ibin] = (TH1D*)hMcE[iregion]->ProjectionY( ( "projection_McE" + to_string( iregion ) + to_string( ibin ) ).c_str(), bin, bin + rebin_size_nhits ); const string E_nhitsname = "E_nhits_reg" + to_string( iregion ) + "_nh" + to_string( (uint)h3dE[iregion]->GetZaxis()->GetBinLowEdge( bin + rebin_size_nhits + 1 ) ); plot_2gauss( &c1, E_nhitsname, h3d_E_nhits[iregion][ibin], h4d_E_nhits[iregion][ibin], "E", "nhits"); } hPosXY[iregion] = (TH2D*)hPos[iregion]->Project3D("yx"); hPosXZ[iregion] = (TH2D*)hPos[iregion]->Project3D("zx"); Det detector( DETECTOR_FILE ); DetDim detdim( detector ); uint idom = 0; while( idom < detdim.get_dom_count() ){ Vec dom_pos = detdim.get_dom_pos( idom++ ); hDetXY[iregion]->Fill( dom_pos.x, dom_pos.y ); hDetXZ[iregion]->Fill( dom_pos.x, dom_pos.z ); //hDetXY[iregion]->Fill( detdim.get_origin().x, detdim.get_origin().y ); //hDetXZ[iregion]->Fill( detdim.get_origin().x, detdim.get_origin().z ); } hPosXY[iregion]->Draw(); hDetXY[iregion]->Draw("same"); hDetXY[iregion]->SetMarkerColor(2); c1.SaveAs( (TString) "rec_plots/hPosXY" + to_string(iregion) + ".root"); c1.SaveAs( (TString) "rec_plots/hPosXY" + to_string(iregion) + ".png"); hPosXZ[iregion]->Draw(); hDetXZ[iregion]->Draw("same"); hDetXZ[iregion]->SetMarkerColor(2); c1.SaveAs( (TString) "rec_plots/hPosXZ" + to_string(iregion) + ".root"); c1.SaveAs( (TString) "rec_plots/hPosXZ" + to_string(iregion) + ".png"); const string dir_nhitsname = "dirbynhits_reg" + to_string( iregion ); plot_th1( &c1, dir_nhitsname, h3d_dirbynhits[iregion], h4d_dirbynhits[iregion], "N_{hits}", "median #Delta angle [deg]", true); } TProfile* p3d_dir_dist = (TProfile*)h3dAngDist->ProfileY("p3ddir"); TProfile* p4d_dir_dist = (TProfile*)h4dAngDist->ProfileY("p4ddir"); plot_th1(&c1, "dir_dist", p3d_dir_dist, p4d_dir_dist, "Shortest distance to can [m]", "Angular reconstruction [deg]", false); plot_th2(&c1, "dir_dist3d", h3dAngDist, "Angular reconstruction", "Shortest distance to can [m]", false, false); plot_th2(&c1, "dir_dist4d", h4dAngDist, "Angular reconstruction", "Shortest distance to can [m]", false, false); TProfile* p3d_Erec_dist = (TProfile*)h3dErecDist->ProfileY("p3dE", 1, -1, "o"); TProfile* p4d_Erec_dist = (TProfile*)h4dErecDist->ProfileY("p4dE", 1, -1, "o"); plot_th1(&c1, "Erec_dist", p3d_Erec_dist, p4d_Erec_dist, "Shortest distance to can [m]", "Energy reconstruction", false, false, false); plot_th2(&c1, "Erec_dist3d", h3dErecDist, "Energy reconstruction", "Shortest distance to can [m]", false, false); plot_th2(&c1, "Erec_dist4d", h4dErecDist, "Energy reconstruction", "Shortest distance to can [m]", false, false); TProfile* p3dStopwatch = (TProfile*)h3dStopwatch->ProfileY("p3dstop"); TProfile* p4dStopwatch = (TProfile*)h4dStopwatch->ProfileY("p4dstop"); plot_th1(&c1, "stopwatch", p3dStopwatch, p4dStopwatch, "Time spent reconstructing [s]", "", false, true); cout << "Plots done." << endl; };