#ifndef PARAMSHH_INCLUDED #define PARAMSHH_INCLUDED #include "constants.hh" #include "../utilJpp.hh" #include "Hit.hh" #include "Vec.hh" #include "io_util.hh" //! Paramaters describing shower-pmt relation struct Params { double E, d, z, theta, phi; Params(): E(0.0), d(0.0), z(0.0), theta(0.0), phi(0.0) {} Params(double E, double d, double z, double theta, double phi) : E(E), d(d), z(z), theta(theta), phi(phi) {} void init( const Trk& trk, const Vec& pos, const Vec& dir ) { // trk.dir.normalize(); // dir.normalize(); const Vec v = pos - trk.pos; const double y = v.dot( trk.dir ); double u = dir.dot( trk.dir ); if ( u > 1.0 ) u = 1.0; if ( u < -1.0 ) u = -1.0; if ( trk.E > 1e14 ) E = 1e14; else if ( trk.E < 1 ) E = 1; else E = trk.E; d = v.len(); if ( d == 0.0 ) z = 0.0; //Pick arbitrary direction else z = y/d; if ( z > 1.0 ) z = 1.0; if ( z < -1.0 ) z = -1.0; theta = acos(u); phi = angle_between ( v - y * trk.dir, dir - u * trk.dir ); if ( phi != phi ) phi = 0.0; assert( E == E ); assert( d == d ); assert( z == z ); assert( phi == phi ); if( theta != theta ){ cout << "u: " << u << endl; cout << trk << endl; cout << pos << endl; cout << dir << endl; } assert( theta == theta ); } Params( const Trk& trk, const Vec& pos, const Vec& dir ) { init( trk, pos, dir ); } Params( const Trk& trk, const Hit& hit ) { init( trk, hit.pos, hit.dir ); } Params( const Trk& trk, const Pmt& pmt ) { init( trk, pmt.pos, pmt.dir ); } // copy Params( const Params& p2 ) { E = p2.E; d = p2.d; z = p2.z; theta = p2.theta; phi = p2.phi; } //! compute parameters corresponding to moving 'length' along the shower direction Params propagate( double length ) const { Params r = *this; const double Z0 = z * d; const double X2 = d*d - Z0*Z0; // X remains unchanged const double Z1 = Z0 - length; r.d = sqrt(Z1*Z1 + X2); r.z = Z1 / r.d; if ( r.z > 1.0 ) r.z = 1.0; if ( r.z < -1.0 ) r.z = -1.0; return r; } void print( std::ostream& out = std::cout ) const { out << " E=" << E << " d=" << d << " z=" << z << " theta="<< theta << " phi=" << phi << endl; } }; //template Params::Params(const Trk&, const Pmt&); // template instantiation ADD_PRINTING( Params ); //! Parameters, including time (residual) struct Paramst : public Params { double t; Paramst( const Params& p, double t ) : Params(p), t(t) {} Paramst() : Params(), t(0.0) {} Paramst(double E, double d, double z, double theta, double phi , double t) : Params( E,d,z,theta,phi), t(t) {} Paramst( const Trk& trk, const Hit& hit ) : Params( trk, hit ), t( hit.t - ( trk.t + d / aa::v_light ) ) {} void init_t( const Trk &trk, const Hit &hit ){ init( trk, hit.pos, hit.dir ); t = hit.t - ( trk.t + d / aa::v_light ); }; // copy Paramst( const Paramst& p2 ) : Paramst( p2.E, p2.d, p2.z, p2.theta, p2.phi, p2.t ) {} Paramst propagate( double length ) const { Params r = Params::propagate( length ); // slicing return Paramst( r, t - length / aa::c_light - ( r.d - d ) / aa::v_light ) ; } void print( std::ostream& out = std::cout ) const { Params::print(); out << " t=" << t ; } }; ADD_PRINTING( Paramst ); Trk propagate( Trk& trk, double length ); Trk get_elong_trk( const Trk &trk, const double fraction, const double precision = 1e-6); #endif