#include "JppTrackPdf.hh" #include #include #include #include #include "util.hh" #include "Jeep/JMessage.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JPhysics/JPDFToolkit.hh" #include "Vec.hh" #include "Mat.hh" #include "JGeometry3D/JPosition3D.hh" #include "JGeometry3D/JRotation3D.hh" #include "JGeometry3D/JDirection3D.hh" #include "JGeometry3D/JAxis3D.hh" #include "rootutil.hh" using namespace JPP; template void load_jpp_pdf( T& pdf , string fn ) { string f = fn; pdf.load(f.c_str() ); pdf.compile(); pdf.setExceptionHandler(new Function1D_t::JDefaultResult(JMATH::zero)); //?? cout << "loaded " << f << endl; } JppTrackPdf::JppTrackPdf( string muon_direct_pdf, string muon_scattered_pdf, string em_direct_pdf, string em_scattered_pdf, string delta_direct_pdf, string delta_scattered_pdf, const double blur, const double time_window_start, const double time_window_end, const double k40_rate, bool v ): _k40_rate( k40_rate ), verbose( v ); { pdfs.resize(6); load_jpp_track_pdf( pdfs[0], muon_direct_pdf ); load_jpp_track_pdf( pdfs[1], muon_scattered_pdf ); load_jpp_track_pdf( pdfs[2], em_direct_pdf ); load_jpp_track_pdf( pdfs[3], em_scattered_pdf ); load_jpp_track_pdf( pdfs[4], delta_direct_pdf ); load_jpp_track_pdf( pdfs[5], delta_scattered_pdf ); } double JppTrackPdf::npe(double E, double D, double cos_pmt) { const uint steps = 100; Vec dir = Vec(-1, 0, 0); const double pmt_a = acos(cos_pmt); std::vector dirs; for(uint i = 0; i < steps; i++){ const double phi = 2 * pi * i / steps; const Vec v = Vec(0, 0, 1).rotate_x(pmt_a).rotate_z(phi); dirs.push_back( v ); } print(dirs); Mat Rot = Mat( dir ); std::vector rot_dirs; for(uint idir = 0; idir < dirs.size(); idir ++){ rot_dirs[idir] = Rot * dirs[idir]; } print(rot_dirs); return 1.0; } double JppTrackPdf::npe(double E, double R, double theta, double phi, double dt ) { const double deltas = getDeltaRaysFromMuon(E); terms = { pdfs[0](R, theta, phi, dt).V , pdfs[1](R, theta, phi, dt).V , pdfs[2](R, theta, phi, dt).V * E , pdfs[3](R, theta, phi, dt).V * E , pdfs[4](R, theta, phi, dt).V * deltas , pdfs[5](R, theta, phi, dt).V * deltas}; return sum(terms); } double JppTrackPdf::dnpe_dt(double E, double R, double theta, double phi, double dt) { const double deltas = getDeltaRaysFromMuon(E); terms = { pdfs[0](R, theta, phi, dt).f , pdfs[1](R, theta, phi, dt).f , pdfs[2](R, theta, phi, dt).f * E , pdfs[3](R, theta, phi, dt).f * E , pdfs[4](R, theta, phi, dt).f * deltas , pdfs[5](R, theta, phi, dt).f * deltas }; return sum(terms); } double JppTrackPdf::d2npe_dt2(double E, double R, double theta, double phi, double dt) { const double deltas = getDeltaRaysFromMuon(E); terms = { pdfs[0](R, theta, phi, dt).fp , pdfs[1](R, theta, phi, dt).fp , pdfs[2](R, theta, phi, dt).fp * E , pdfs[3](R, theta, phi, dt).fp * E , pdfs[4](R, theta, phi, dt).fp * deltas, pdfs[5](R, theta, phi, dt).fp * deltas}; return sum(terms); } // double JppTrackPdf::npe( const Trk& trk, // const Hit& hit ) // { // HitParams p( trk, hit); // return npe( trk.E, p.R, p.theta, p.phi, p.dt ); // } // double JppTrackPdf::dnpe_dt( const Trk& trk, // const Hit& hit ) // { // HitParams p( trk, hit); // return dnpe_dt( trk.E, p.R, p.theta, p.phi, p.dt ); // } // double JppTrackPdf::d2npe_dt2( const Trk& trk, // const Hit& hit ) // { // HitParams p( trk, hit); // return d2npe_dt2( trk.E, p.R, p.theta, p.phi, p.dt ); // } // double JppTrackPdf::eval_L(Trk& trk){ // { // if (trk.E > 1e14 ) trk.E = 1e14; // //const double alpha = trk.E * 1e-6; // energy scale factor ( H3 is for 1 PeV ). // double L=0; // double L0=0; // double L1=0; // int nhit =0; // foreach( dom, domsums ) // { // { // TIMER( hit-pmt-loop ) // vector::const_iterator ipmtdir = dom.hit_pmts.begin(); // vector< vector >::const_iterator ihit_times = dom.hit_times.begin(); // for(; // ipmtdir != dom.hit_pmts.end(); // ipmtdir++, ihit_times++){ // nhit ++; // Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop // double mu1 = npe(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi); // // foreach(hit_time, *ihit_times){ //Need to take into account same pmt being hit several times // //Hit_params this_hit(trk, dom.pos, *ipmtdir, hit_time);//much of this can be pulled out of loop // //double mu1 = dnpe_dt(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi, this_hit.dt); // const double mu = mu1 + k40prob; // double logp1 = log( 1 - exp( - mu ) ); // //if (do_clamp) clamp( logp1, minlogp, maxlogp ); // L -= logp1; // L1 -= logp1; // } // } // { // TIMER( unhit-pmt-loop ) // for(vector::const_iterator ipmtdir = dom.hit_pmts.begin(); // ipmtdir != dom.hit_pmts.end(); ipmtdir++){ // Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop // double mu1 = npe(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi); // const double mu = mu1 + k40prob; // double logp0 = -mu; // //if (do_clamp) clamp( logp1, minlogp, maxlogp ); // L -= logp0; // L0 -= logp0; // } // } // if ( L != L ) // { // print("likelihood is nan @ r,z,a ="); // print("trk Energy: ", trk.E); // print("Likelihood hits:", L1); // print("Likelihood no hits:", L0); // } // } // dom // //test_params(trk); //A-ok // return L; // } // }; // void JppTrackPdf::test_params(Trk& trk){ // { // foreach( dom, domsums ) // { // { // TIMER( hit-pmt-loop ) // vector::const_iterator ipmtdir = dom.hit_pmts.begin(); // vector< vector >::const_iterator ihit_times = dom.hit_times.begin(); // for(; // ipmtdir != dom.hit_pmts.end(); // ipmtdir++, ihit_times++){ // foreach(hit_time, *ihit_times){ //Need to take into account same pmt being hit several times // Hit_params this_hit(trk, dom.pos, *ipmtdir, hit_time);//much of this can be pulled out of loop // if(round(this_hit.D * 1000) != round(calc_D(trk.pos, dom.pos) * 1000)) // print("test_params: bad D, ", this_hit.D , " vs. ", calc_D(trk.pos, dom.pos)); // if(round(this_hit.R * 1000) != round(calc_R(trk.pos, dom.pos, trk.dir) * 1000)) // print("test_params: bad R, ", this_hit.R , " vs. ", calc_R(trk.pos, dom.pos, trk.dir)); // if(round(this_hit.cd * 1000) != round(calc_cd(trk.pos, dom.pos, trk.dir) * 1000)) // print("test_params: bad cd, ", this_hit.cd , " vs. ", calc_cd(trk.pos, dom.pos, trk.dir)); // if(round(this_hit.dt * 1000) != round(calc_dt(trk.pos, dom.pos, trk.t, hit_time) * 1000)) // print("test_params: bad dt, ", this_hit.dt , " vs. ", calc_dt(trk.pos, dom.pos, trk.t, hit_time)); // } // } // } // { // TIMER( unhit-pmt-loop ) // for(vector::const_iterator ipmtdir = dom.hit_pmts.begin(); // ipmtdir != dom.hit_pmts.end(); ipmtdir++){ // Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop // if(round(this_hit.D * 1000) != round(calc_D(trk.pos, dom.pos) * 1000)) // print("test_params: bad D, ", this_hit.D , " vs. ", calc_D(trk.pos, dom.pos)); // if(round(this_hit.R * 1000) != round(calc_R(trk.pos, dom.pos, trk.dir) * 1000)) // print("test_params: bad R, ", this_hit.R , " vs. ", calc_R(trk.pos, dom.pos, trk.dir)); // if(round(this_hit.cd * 1000) != round(calc_cd(trk.pos, dom.pos, trk.dir) * 1000)) // print("test_params: bad cd, ", this_hit.cd , " vs. ", calc_cd(trk.pos, dom.pos, trk.dir)); // } // } // } // dom // } // };