//Author: Jordan Seneca //Use: ./HistRecoTrks inpath* outfile #include #include #include #include #include //#include "/sps/km3net/users/jseneca/docopt.cpp/docopt.h" #include "util.hh" #include "Vec.hh" #include "Mat.hh" #include "Trk.hh" #include "Evt.hh" #include "EventFile.hh" #include "TH1D.h" #include "TH2D.h" #include "TH3D.h" #include "TString.h" #include "TTree.h" #include "TTreeReader.h" #include "TTreeReaderValue.h" #include "TFile.h" #include "TCanvas.h" #include "HistRecoTrks.hh" static const char USAGE[] = R"( HistRecoTrks. Usage: HistRecoTrks )"; using namespace std; int main(int argc, char *argv[] ){ /* map args = docopt::docopt(USAGE, { argv + 1, argv + argc }, true, "HistRecoTrks 0.1"); for(auto const& arg: args){ cout<< arg.first << arg.second << endl; } */ cout << "File path to be searched: " << argv[1] << endl ; TFile outfile(argv[2], "RECREATE", "outfile"); string DETECTOR_FILE( argv[3] ); std::cout << "Still alive here?" << std::endl; const uint region_n = 3; //All compared against Energy for TH2D and Energy + NHits for TH3D 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++ ){ vector< Double_t > angbin_vec = get_log_vector(1e-3, 180, 1000); for (Int_t i=0;i<=1000;i++) { } const vector< Double_t > Ebin_vec = get_lin_vector(1e5, 1e8, 1000); const vector< Double_t > hitbin_vec = get_lin_vector(1e3, 1e5, 100); h3dAng[iregion] = new TH3D( (TString)"h3dAngENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]); h4dAng[iregion] = new TH3D( (TString)"h4dAngENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]); h4dAngOld[iregion] = new TH3D( (TString)"h4dAngOldENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]); h3dErec[iregion] = new TH3D( (TString)"h3dErecENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5); h4dErec[iregion] = new TH3D( (TString)"h4dErecENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5); h4dErecOld[iregion] = new TH3D( (TString)"h4dErecOldENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5); h3dE[iregion] = new TH2D( (TString)"h3dENhits" + to_string(iregion), "3d E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5); h4dE[iregion] = new TH2D( (TString)"h4dENhits" + to_string(iregion), "4d E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5); hMcE[iregion] = new TH2D( (TString)"hMcENhits" + to_string(iregion), "mc E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5); hPos[iregion] = new TH3D( (TString)"hPosition" + to_string(iregion), "pos mapping; pos x [m]; pos y [m]; pos z [m]", 300, -600, 600, 300, -600, 600, 400, 0, 800); //Don't make this histo too big or the file will reach the byte limit } TH2D* h3dAngDist = new TH2D("h3dAngDist", " ;delta ang[deg]; d_can [m]", 1000, 0, 180, 200, -300, 550); TH2D* h4dAngDist = new TH2D("h4dAngDist", " ;delta ang[deg]; d_can [m]", 1000, 0, 180, 200, -300, 550); TH2D* h3dErecDist = new TH2D("h3dErecDist", "E res comparison; delta E[deg]; d_can [m]", 200, -2, 2, 200, -300, 550); TH2D* h4dErecDist = new TH2D("h4dErecDist", "E res comparison; delta E[deg]; d_can [m]", 200, -2, 2, 200, -300, 550); //This against HitN TH2D* h3dStopwatch = new TH2D("h3dStopwatch", " ;reco_T [s]; hit_N", 1000, 0, 1e3, 1000, 0, 1e5); TH2D* h4dStopwatch = new TH2D("h4dStopwatch", " ;reco_T [s]; hit_N", 30, 0, 1e3, 1000, 0, 1e5); vector filelist = {}; glob_t globbuf; glob(argv[1], 0, NULL, &globbuf); for(uint i = 0; i < globbuf.gl_pathc; i++){ filelist.push_back( globbuf.gl_pathv[i] ); cout << "Added file: " << globbuf.gl_pathv[i] << endl; } cout << "Filelist length: " << filelist.size() << endl; Vec coord_origin; for(vector::const_iterator filename = filelist.begin(); filename != filelist.end(); filename++){ EventFile file( (*filename).c_str() ); cout << "Found event file " << *filename << endl; Head header = file.header; Det detector( DETECTOR_FILE ); DetDim detdim( detector ); cout << "Height: " << detdim.height() << endl; cout << "Radius: " << detdim.radius() << endl; try{ coord_origin = header.coord_origin(); } catch(std::invalid_argument err){ //cout << "Error, :" << err << endl; coord_origin = detector.get_cg(); } uint n_events = 0; uint n_contained = 0; uint n_inside = 0; uint n_outside = 0; uint iregion = 0; //0 contained, 1 inside, 2 outside for( EventFile::iterator ievt = file.begin(); ievt != file.end(); ievt++ ){ Evt evt = *ievt; n_events += 1; Trk mc_neut = evt.mc_trks[0]; mc_neut.pos += coord_origin; const double dist_edge = get_dist_to_edge( mc_neut, detector, detdim ); if( dist_edge < -100 ){ //cout << "Event contained: " << dist_edge << endl; iregion = 0; n_contained++; } else if( dist_edge < 0 ){ //cout << "Event inside: " << dist_edge << endl; iregion = 1; n_inside++; } else { //cout << "Event outside: " << dist_edge << endl; iregion = 2; n_outside++; }; //cout << "evt.hits.size(): " << evt.hits.size() << endl; hMcE[iregion]->Fill( mc_neut.E, evt.hits.size() ); hPos[iregion]->Fill( mc_neut.pos.x, mc_neut.pos.y, mc_neut.pos.z ); for( vector::iterator trk =evt.trks.begin(); trk != evt.trks.end(); trk ++ ){ if(trk->rec_type == 103){ Reso myres( *trk, mc_neut); h3dAng[iregion] ->Fill( myres.ang_res, mc_neut.E, evt.hits.size() ); h3dErec[iregion]->Fill( myres.E_res , mc_neut.E, evt.hits.size() ); h3dE[iregion] ->Fill( trk->E, evt.hits.size() ); //h3dStopwatch ->Fill( 0, evt.hits.size() ); h3dAngDist ->Fill( myres.ang_res, dist_edge); h3dErecDist ->Fill( myres.E_res , dist_edge); continue; } if(trk->rec_type == 104){ Reso myres(*trk, mc_neut); h4dAng[iregion] ->Fill( myres.ang_res, mc_neut.E, evt.hits.size() ); h4dErec[iregion]->Fill( myres.E_res , mc_neut.E, evt.hits.size() ); h4dE[iregion] ->Fill( trk->E, evt.hits.size() ); h4dStopwatch ->Fill( evt.getusr("newshowerfit_time"), evt.hits.size() ); h4dAngDist ->Fill( myres.ang_res, dist_edge); h4dErecDist ->Fill( myres.E_res , dist_edge); continue; } if(trk->rec_type == 114){ Reso myres(*trk, mc_neut); h4dAngOld[iregion] ->Fill( myres.ang_res, mc_neut.E, evt.hits.size() ); h4dErecOld[iregion]->Fill( myres.E_res , mc_neut.E, evt.hits.size() ); continue; } } } cout << "analysed " << n_events << " events" << endl; cout << "contained: " << n_contained << " events" << endl; cout << "inside : " << n_inside << " events" << endl; cout << "outside : " << n_outside << " events" << endl; } cout << "Writing file to " << argv[2] << "... " << endl; outfile.Write(); cout << "Closing file... " << endl; outfile.Close(); cout << "Hists done." << endl; };