#include "Trk.hh" #include "Hit.hh" #include "Vec.hh" #include "Evt.hh" #include "EventFile.hh" #include "JppShowerPdf.hh" #include "JppShowerFit.hh" #include "JppShowerPdf_simple.hh" #include "JppShowerFit_simple.hh" #include "rec2/CascadePdf.hh" #include "Det.hh" #include "rec2/JppGrad.hh" int main( int argc, int** argv ){ EventFile f( "/sps/km3net/users/jseneca/aanet/data/km3_v5/mcv5.1.genhen_nueCC.km3_AAv1.jte.jchain.aashower.94.root" ); f.set_index( 2147 ); Det det( "/sps/km3net/users/jseneca/aanet/data/km3_v5/KM3NeT_-00000001_20171212.detx" ); Trk trk = f.evt.mc_trks[0]; vector pre_hits = f.evt.hits; det.apply( pre_hits ); vector hits; for( auto h = pre_hits.begin(); h != pre_hits.end(); h++ ){ Hit hit = *h; hit.t = hit.t - ( f.evt.mc_t - f.evt.t.AsDouble()*1e9 ); double D = ( trk.pos - hit.pos ).len(); double dt = abs( hit.t - D/v_light ); if( D < 500 && dt < 100 ){ hits.push_back( hit ); } }; cout << "len hits: " << hits.size() << endl; assert( hits.size() ); CascadePdf cascade_pdf( "/pbs/throng/km3net/src/Jpp/master/data/J14p.dat", "/pbs/throng/km3net/src/Jpp/master/data/J13p.dat", 0, 0, //0 : do nothing, 1 = move to shower max, > 1 : sample 5000, //buggy behavior when no background -100, 900, 0.0 ); print( "Made cascadepdf" ); JppShower_simple gradpdf( cascade_pdf, 20 ); print( "Made JppShower_simple" ); gradpdf.init_event( det, hits, trk.pos, 300, true); JppShowerFit_simple fitter_gradpdf( "5d_gradpdf", &gradpdf, "Minuit2", "MIGRAD", 1 ); print( "Made JppShowerFit_simple" ); JppShower pdf( "/pbs/throng/km3net/src/Jpp/master/data/J14p.dat", //Scattered PDF must be first. "/pbs/throng/km3net/src/Jpp/master/data/J13p.dat", 0.0, // blurring of the PDF -100, 900, false, 20, 7e3, 0); pdf.init_event( det, hits, trk.pos, 300, true); JppShowerFit fitter_pdf( "5d_pdf", &pdf, "Minuit2", "MIGRAD", 1 ); print("JppShowerFit ok"); aa_legacy_timer::Timer t1( " GRAD PDF " ); Trk fit_gradpdf = fitter_gradpdf.fit( trk, det , hits, true, true ); aa_legacy_timer::Timer t0( " PDF " ); Trk fit_pdf = fitter_pdf.fit( trk, det , hits, true, true ); print( "MC neutrino: ", trk ); print( "Fit pdf: ", fit_pdf ); print( "Fit gradpdf: ", fit_gradpdf ); t1.report(); };