#include "JppGrad.hh" #include "stringutil.hh" JppGrad::JppGrad( vector< string > shower_pdfs, double sigma_blur ,// = 0, double R_bg ,// background rate in [Hz] double t_start ,// start of hit time residual window double t_end, uint verbose, uint debug ): _background_rate( R_bg ), _window_start( t_start ), _window_end( t_end ), _verbose( verbose ), _debug( debug ), _dx_dtes(), _dXs( 20, vector< vector < double > >( 3, vector< double >( 6, 0.0 ) ) ), _sum_dnpes_a( 7, vector( 3, 0 ) ), _sum_dXs_a ( 7, vector( 3, 0 ) ), _sum_dTs_a ( 7, vector( 6, 0 ) ), _sum_dTes_a ( 7, vector( 7, 0 ) ), _sum_dnpes_n( 7, vector( 3, 0 ) ), _sum_dXs_n ( 7, vector( 3, 0 ) ), _sum_dTs_n ( 7, vector( 6, 0 ) ), _sum_dTes_n ( 7, vector( 7, 0 ) ), _dL_a( 7, 0.0 ), _dL_n( 7, 0.0 ){ for( uint ii = 0; ii != 7; ii++ ){ uint iii = ii >= 4 ? ii + 3 : ii; //Skip over the cart direction dims returned by jacobian objects kk_kki[ii] = iii; jj_jji[ii] = iii; }; bool build = !file_exists("totalpdf.dat"); if (build) { cout << "Building JppGrad pdf from direct and scattered constituents..." << endl; cout << shower_pdfs[0] << endl; cout << shower_pdfs[1] << endl; cout << sigma_blur << endl; cout << R_bg << endl; cout << t_start << endl; cout << t_end << endl; Timer tim("constructPDF", true); _pdf_gradfuncs.resize(2); _npe_gradfuncs.resize(1); print ("loading scattered photon table."); _pdf_gradfuncs[0].load( shower_pdfs[0].c_str() ); print ("loading direct photon table."); _pdf_gradfuncs[1].load( shower_pdfs[1].c_str() ); print ("computing pdfs"); for( vector< JppShowerGradPdf_t >::iterator ipdf = _pdf_gradfuncs.begin(); ipdf != _pdf_gradfuncs.end(); ++ipdf ) { std::cout << "x" << std::flush; ipdf->transform(JppShowerGradPdf_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function ipdf->setExceptionHandler( new JppShowerGradPdf_t::JDefaultResult( JMATH::zero ) ); auto x = new PDG_Function1D_t::JDefaultResult( JMATH::zero ) ; for (JppShowerGradPdf_t::super_iterator p = ipdf->super_begin(); p != ipdf->super_end(); ++p) { p.getValue().setExceptionHandler(x); } }; print ("adding pdfs"); add_jpp_shower_pdfs( _pdf_gradfuncs[0], _pdf_gradfuncs[1] ); print("Making npe table."); _npe_gradfuncs[0] = JppShowerGradNpe_t( _pdf_gradfuncs[0] ); //Make NPE table for total only print("smearing combined table with sigma=", sigma_blur ); if ( sigma_blur > 0 ){ _pdf_gradfuncs[0].blur( sigma_blur ); // only blur total : save some time } print("background rate:", R_bg, "Hz"); print("time window:", t_start, t_end ); print("JppGrad constructed in", tim.value(), "sec"); _pdf_gradfunc = &_pdf_gradfuncs[ 0 ]; //Could remove pdf_funcs vector if we really don't need components _npe_gradfunc = &_npe_gradfuncs[ 0 ]; print ("saving (?)"); _pdf_gradfunc->store("totalpdf.dat"); } else // don't build, but load { print("loading totalpdf.dat"); _pdf_gradfunc = new JppShowerGradPdf_t; _pdf_gradfunc->load("totalpdf.dat"); _pdf_gradfunc->transform(JppShowerGradPdf_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function _pdf_gradfunc->setExceptionHandler( new JppShowerGradPdf_t::JDefaultResult( JMATH::zero ) ); print("building npe table"); _npe_gradfunc = new JppShowerGradNpe_t( *_pdf_gradfunc ); } print ("done with JppGrad ctor."); }; void JppGrad::eval_lik_gradient(double &r, vector &ds, const Trk &trk, vector hits, vector unhits) { TIMEFUNC // Timer time0("time0"); // out debug : 1.66632e-07 // set elong : 0.000437399 // fill ds : 3.44014e-07 // hit loop : 0.00382506 // unhit loop: 1.06737e-07 //----------inner---------- // thit norm : 1.72007e-07 // thit paramst : 3.31191e-06 // thit jacob : 0.000140122 // thit A : 0.00121552 // thit dnpedtdts : 0.00123817 // thit dnpedts : 0.00122572 // thit dte : 1.07735e-06 // thit dhitL : 2.48028e-07 // thit sanity : 2.61082e-07 // thit debug : 1.55114e-07 // tunhit norm : 0 // tunhit params : 0 // tunhit jacob : 0 // tunhit dnpedts : 0 // tunhit dunhitL : 0 // double t0 = time0.time_since_last_start(); if (_debug) { bad_dXs = 0; bad_int_dXs = 0; bad_dtdXs = 0; call_dXs = 0; call_int_dXs = 0; call_dtdXs = 0; pdf_call_count = 0; npe_call_count = 0; bad_pdf_call_count = 0; bad_npe_call_count = 0; if (_debug > 1) { _sum_dnpes_a.assign(7, vector(3, 0)); _sum_dXs_a.assign(7, vector(3, 0)); _sum_dTs_a.assign(7, vector(6, 0)); _sum_dTes_a.assign(7, vector(7, 0)); _sum_dnpes_n.assign(7, vector(3, 0)); _sum_dXs_n.assign(7, vector(3, 0)); _sum_dTs_n.assign(7, vector(6, 0)); _sum_dTes_n.assign(7, vector(7, 0)); _dL_a.assign(7, 0.0); _dL_n.assign(7, 0.0); }; }; // double t1 = time0.time_since_last_start(); { Timer _ ("set_elon"); set_elongation(trk); } ds.assign(_N, 0.0); r = 0.0; foreach (hit, hits) { TIMER(hitloop) // Here we access the pdf for the first and only time set_paramst(hit); // 0.96%, 1.74e-7 if (hit.pattern_flags == 2) { continue; // Hit flagged as bad } set_jacobians(hit); // 39.1%, 7.1e-6 // double t7 = time0.time_since_last_start(); // thit_jacob += t7 - t6; const vector dnpe_dtdTs = this->dnpe_dtdTs(hit); // double t8 = time0.time_since_last_start(); // thit_dnpedtdts += t8 - t7; vector dnpe_dTs; if (npe(hit) < _fudge2) dnpe_dTs = this->dnpe_dTs(hit); // 26.0%, 4.7e-6 else dnpe_dTs.assign(7, 0.0); // double t9 = time0.time_since_last_start(); // thit_dnpedts += t9 - t8; _A = dnpe_dt(hit); // double t10 = time0.time_since_last_start(); // thit_A += t10 - t9; r -= get_hit_L(_A, npe(hit), true); // double t10b = time0.time_since_last_start(); // thit_hitL += t10b - t10; for (uint ii = 0; ii != _N; ii++) // gradient { // double t11 = time0.time_since_last_start(); double dt_e = 0; for (vector>>::iterator idTs = _dX_dTes.begin(); idTs != _dX_dTes.end(); idTs++) { dt_e += (*idTs)[5][ii] / (float)_elong_trks.size(); }; // double t12 = time0.time_since_last_start(); // thit_dte += t12 - t11; ds[ii] -= dhit_L(dnpe_dtdTs[ii], dnpe_dTs[ii], _A, dt_e); // double t13 = time0.time_since_last_start(); // thit_dhitL += t13 - t12; if (ds[ii] < -1e6 or ds[ii] != ds[ii]) { cout << "ds[ii]: " << ds[ii] << endl; cout << "bad LL grad ii: " << ii << endl; cout << "dnpe_dtdTs: " << dnpe_dtdTs[ii] << endl; cout << "dnpe_dTs: " << dnpe_dTs[ii] << endl; cout << "dt_e: " << dt_e << endl; cout << "_A: " << _A << endl; cout << "trk: " << trk << endl; cout << "hit: " << hit << endl; if (_debug) cin.ignore(); }; // double t14 = time0.time_since_last_start(); // thit_sanity += t14 - t13; if (_debug > 1) { print("debug"); fill_components(trk, hit); // Must come first as the trk is set double dnpe_dtdT = this->num_dnpe_dtdTs(trk, hit, 1e-5)[ii]; double dnpe_dT; if (npe(hit) < _fudge2) dnpe_dT = this->num_dnpe_dTs(trk, hit, 1e-5)[ii]; else dnpe_dT = 0.0; double num_dt_e = 0.0; for (vector::iterator itrk = _elong_trks.begin(); itrk != _elong_trks.end(); ++itrk) { num_dt_e += num_dX_dTs(*itrk, hit, 1e-5)[5][ii] / (float)_elong_trks.size(); }; _dL_n[ii] -= dhit_L(dnpe_dtdT, dnpe_dT, _A, num_dt_e); } // double t15 = time0.time_since_last_start(); // thit_debug += t15 - t14; } } // double t16 = time0.time_since_last_start(); /* double tunhit_norm = 0.0; double tunhit_params = 0.0; double tunhit_jacob= 0.0; double tunhit_dnpedts = 0.0; double tunhit_dunhitL = 0.0; double tunhit_unhitL = 0.0; */ foreach (unhit, unhits) { TIMER(unhit) set_params(unhit); // has to be first if (unhit.flag == 2) continue; // Hit flagged as bad // double t19 = time0.time_since_last_start(); // tunhit_params += t19 - t18; set_jacobians(unhit); // double t20 = time0.time_since_last_start(); // tunhit_jacob += t20 - t19; const vector dnpe_int_dTs = this->dnpe_int_dTs(unhit); // double t21 = time0.time_since_last_start(); // tunhit_dnpedts += t21 - t20; r -= -npe(unhit); // double t21c = time0.time_since_last_start(); // tunhit_unhitL += t21c - t21; for (uint ii = 0; ii != _N; ii++) { // double t21a = time0.time_since_last_start(); double dL_unhit = -dnpe_int_dTs[ii]; ds[ii] -= dL_unhit; // double t21b = time0.time_since_last_start(); // tunhit_dunhitL += t21b - t21a; if (ds[ii] != ds[ii]) { cout << "nan encountered in ll grad ii: " << ii << endl; cout << "trk: " << trk << endl; cout << "unhit: " << unhit << endl; if (_debug) cin.ignore(); }; if (_debug > 1) _dL_n[ii] -= -this->num_dnpe_dTs(trk, unhit, 1e-5)[ii]; // Yes, 2 negative signs }; if (_debug > 1) fill_components(trk, unhit); } if (_debug) { cout << "Good pdf calls: " << pdf_call_count - bad_npe_call_count << ", bad: " << bad_npe_call_count << endl; cout << "Good npe calls: " << npe_call_count - bad_pdf_call_count << ", bad: " << bad_pdf_call_count << endl; // cout << "Good calls to dnpe_int_dXs: " << call_int_dXs - bad_int_dXs << ", bad: " << bad_int_dXs << endl; } // double t23 = time0.time_since_last_start(); /* cout << "out debug : " << t1 - t0 << endl; cout << "set elong : " << t2 - t1 << endl; cout << "fill ds : " << t3 - t2 << endl; cout << "hit loop : " << t16 - t3 << endl; cout << "unhit loop: " << t23 - t16 << endl; cout << "grand tot : " << t23 - t0 << endl; cout << "----------inner----------" << endl; cout << "thit norm : " << thit_norm << endl; cout << "thit paramst : " << thit_paramst << endl; cout << "thit jacob : " << thit_jacob << endl; cout << "thit A : " << thit_A << endl; cout << "thit dnpedtdts : " << thit_dnpedtdts << endl; cout << "thit dnpedts : " << thit_dnpedts << endl; cout << "thit dte : " << thit_dte << endl; cout << "thit dhitL : " << thit_dhitL << endl; cout << "thit dhitL : " << thit_hitL << endl; cout << "thit sanity : " << thit_sanity << endl; cout << "thit debug : " << thit_debug << endl; cout << "tunhit norm : " << tunhit_norm << endl; cout << "tunhit params : " << tunhit_params << endl; cout << "tunhit jacob : " << tunhit_jacob << endl; cout << "tunhit dnpedts : " << tunhit_dnpedts << endl; cout << "tunhit dunhitL : " << tunhit_dunhitL << endl; cout << "tunhit unhitL : " << tunhit_unhitL << endl; cin.ignore(); */ };