#ifndef __JPP_SHOWER_FIT__ #define __JPP_SHOWER_FIT__ #include "Evt.hh" #include "showrec_util.hh" #include "JppShowerPdf.hh" #include "Math/GSLMinimizer.h" #include "Math/Minimizer.h" #include "Math/Factory.h" #include "Math/Functor.h" #include "TRandom2.h" #include "TError.h" #include "TStopwatch.h" #include "Fit/FitConfig.h" #include "Math/GSLMinimizer.h" template< class pdftype > class JppShowerFit { protected: uint verbose; double total_time; int nfits; uint call_count; vector fixed_vars; ROOT::Math::Functor _score_functor; //------------------------------------------- // minimizer and recording control //------------------------------------------- string _minimizer_name; string _minimizer_algo; uint _strategy; bool record; double tolerance; double eval_offset; //------------------------------------------- // Hit and DOM selection //------------------------------------------- double max_dist; bool incl_unhits; // -------------------------- // The Pdf used in the fit //--------------------------- JppShower< pdftype > *pdf; JppUtil ju; public: ~JppShowerFit(){}; JppShowerFit(){}; JppShowerFit( JppShower< pdftype > *pdf_, string minimizer_name = "GSLMultiMin", string minimizer_algo = "BFGS2", uint strategy = 1, uint verb = 0): pdf(pdf_){ _init(minimizer_name, minimizer_algo, strategy, verb); }; void _init_event( Trk& start_track , Det& det, vector& hits ) { pdf->init_event( det, hits, start_track.pos, max_dist, incl_unhits ); _score_functor = ROOT::Math::Functor( this, &JppShowerFit::score_function, 7 ); }; void _init(string minimizer_name, string minimizer_algo, uint strategy, uint verb ) { nfits = 0; call_count = 0; total_time = 0; verbose = verb; record = false; fixed_vars = pack( 7 ); tolerance = 1000; _minimizer_name = minimizer_name; //Minuitx2 _minimizer_algo = minimizer_algo; _strategy = 1; max_dist = 300; incl_unhits = true; }; void init_offset( const double* pars ) { static Trk trk_off; ju.to_trk( pars, trk_off); try { eval_offset = pdf->eval_lik(trk_off); } catch ( JLANG::JException err ){ cout << "Error in score function: " << err << ", setting result to 0." << endl; eval_offset = 0; } }; double score_function( const double* pars) // delegate to pdf->eval { call_count += 1; static Trk trk; ju.to_trk( pars, trk ); double res = pdf->eval_lik(trk) - eval_offset; if( verbose ) cout << trk << ", L: " << res << endl; if( std::isnan(res) ) return 1e10; return res; }; // -------------------------- // which parameters to fix //--------------------------- void fix_vars( int v1, int v2=-1, int v3=-1, int v4=-1, int v5=-1, int v6=-1, int v7=-1 ) { fixed_vars = pack( v1, v2, v3, v4, v5, v6, v7 ); }; void set_incl_unhits( bool setting = true ){ incl_unhits = setting; }; void set_verbose( uint setting = 1 ){ verbose = setting; }; void set_max_dist( double dist = 0.0 ){ max_dist = dist; }; double get_max_dist(){ return max_dist; }; void set_strategy( uint strategy = 1 ){ _strategy = strategy; }; Trk fit( Trk& start_track , Det& det, vector& hits, int do_init_event /*=1*/ ) { TStopwatch W; if( do_init_event ) _init_event( start_track, det, hits ); //cout << " name: " << _minimizer_name.c_str() << endl; //cout << " algo: " << _minimizer_algo.c_str() << endl; ROOT::Math::Minimizer* M = ROOT::Math::Factory::CreateMinimizer( _minimizer_name.c_str(), _minimizer_algo.c_str()); M->SetMaxFunctionCalls(10000); // for Minuit/Minuit2 M->SetMaxIterations(30); // for GSL M->SetTolerance( tolerance ); //Determines the point at which minimization should stop M->SetPrintLevel( verbose * 999 - 1 ); // -1 seems to be the lowser setting in Gslminimizer JppGradFunc< pdftype > gf; gf.pdf = pdf; //Variables are: [ log(energy), posx, posy, posz, theta, phi, time] //See showrec_util.hh -> to_trk double vstart[7]; ju.to_pars( start_track, vstart ); vector step = pack( 0.5, 1.0, 1.0, 1.0, 0.2, 0.2, 5.0 ); gf.init_offset( vstart ); //Needs to be called before the clone next line //M->SetFunction( gf ); M->SetFunction( static_cast &>(gf)); if (verbose) print ( "fixing ", fixed_vars ); for (int i=0; i<7; i++ ) { if ( contains( fixed_vars, i ) ) { M -> SetFixedVariable( i, ju.parnames()[i], vstart[i] ); } else { M -> SetVariable( i, ju.parnames()[i], vstart[i], step[i] ); if (verbose) cout << "Start " << ju.parnames()[i] << ": " << vstart[i] << endl; if( _minimizer_name != "GSLMultiMin" ) M -> SetVariableLimits( i, ju.parlimits()[i].first, ju.parlimits()[i].second ); if( _minimizer_name != "GSLMultiMin" && verbose) cout << "Limits " << ju.parnames()[i] << ": " << ju.parlimits()[i].first << ", " << ju.parlimits()[i].second << endl; } } // ---------------------------------------- // do the fit // ---------------------------------------- //ErrorHandlerFunc_t oldhandler = SetErrorHandler( fit_errorhandler ); cout << "Starting minimization here" << endl; M -> SetStrategy( _strategy ); //0 faster, 2 more accurate M -> Minimize(); //SetErrorHandler( oldhandler ); Trk r; ju.to_trk( M -> X(), r ); r.lik = M -> MinValue(); r.fitinf.resize(1); r.fitinf[0] = W.RealTime(); total_time += r.fitinf[0]; nfits++; r.error_matrix.resize(49); M -> GetCovMatrix( &r.error_matrix[0] ); M->PrintResults(); cout << "strategy: " << M->Strategy() << endl; if(verbose) cout << "True " << start_track << ", L: " << pdf->eval_lik(start_track) << endl; if(verbose) cout << "Fit " << r << ", L: " << pdf->eval_lik(r) << endl; if(verbose) cout << "Fit lik" << r.lik << ", Time seconds: " << r.fitinf[0] << endl; if(verbose) cout << "Iterations: " << M->NIterations() << ", Function calls: " << M->NCalls() << endl; delete M; return r; } }; inline void use_jppshowerfitmerged(){ //Needed to make templates JppShowerFit c1; JppShowerFit c2; }; #endif