#ifndef DETRESPONSE_INCLUDED_HH #define DETRESPONSE_INCLUDED_HH #include "rootutil.hh" #include "ana/search/EffectiveArea.hh" #include "TRandom.h" struct DetResponse { string name; EffectiveArea Aeff; DetResponse() {} // for background virtual double bg_dP_dlogErec_dOmega( double logE_rec, double cos_theta_rec ) {return 0;} // the following are for signal virtual double dP_dlogErec_dOmega_zenith( double E_rec, double cos_theta_rec, double E_true, double cos_theta_true) {return 0;}; virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error, double E_true, double cos_theta_true) {return 0;}; virtual double dP_dlogErec_dOmega_zenith_bins(int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true) {return 0;} virtual ~DetResponse() {} void dus(); ClassDef( DetResponse, 1); }; /*! read object from root file... reasonably */ template void read( TFile& f , T& dst , string name ) { TObject* p = f.Get(name.c_str()); if (!p) fatal("failed to find object called", name, "in TFile", f.GetName()); dst = *((T*) p); } struct DefaultDetResponse : public DetResponse { TH3D h_COSZENREC; // x = etrue, y=costhetatrue, z =costheta_rec TH3D h_EERR; // x = etrue, y=costheattrue, z =E_rec TH3D h_ANGERR; // x = etrue, y=costhetatrue, z =log(angle-muon-rec-track) TH2D h_ATMNU; // x = etrue, y=costhetatrue TH2D h_ATMNUREC; // x = erec y=costhetarec double n_background; // number of background events from h_ATMMU TH2D A; // Aeff DefaultDetResponse() {} DefaultDetResponse( string ingredients_file, string flavour = "numuCC") { TFile f( ingredients_file.c_str() ); read( f, A , "h_EFFAREA" + flavour ); Aeff.Haeff = &A; read( f, h_COSZENREC ,"h_COSZENREC"+flavour ); read( f, h_ANGERR ,"h_ANGERR"+flavour ); read( f, h_EERR ,"h_EERR"+flavour ); read( f, h_ATMNU ,"h_ATMNU"+flavour ); // these are pdfs of the quantity on the z-axis. In other words // they are dP/ ( dlogE dcostheta). // So they must be normalized to that sum ( content * binsize ) = 1 normalize_z( h_COSZENREC ,1.0, true ); normalize_z( h_ANGERR ,1.0, true ); normalize_z( h_EERR ,1.0, true ); check( h_COSZENREC ); check( h_ANGERR ); check( h_EERR ); check( h_ATMNU ); h_ATMNUREC = smeared_bin( h_ATMNU, [this](int a, int b, int c, int d) { double x = this->dP_dlogErec_dOmega_zenith_bins(a,b,c,d); //print( a,b,c,d, x ); return x; }); n_background = h_ATMNU.Integral(); h_ATMNUREC.SetName("h_ATMNUREC"); normalize( h_ATMNUREC ); } /*! the probability P( Erec, cos(theta)_rec | Etrue , cos(theta)_true ) is factorised as P( cos(theta)_rec | Etrue , cos(theta)_true ) * P( Erec | Etrue , cos(theta)_true ) */ virtual double dP_dlogErec_dOmega_zenith(double E_rec, double cos_theta_rec, double E_true, double cos_theta_true) { double a = interpolate( &h_EERR , E_true, cos_theta_true, E_rec ) ; double b = interpolate( &h_COSZENREC, E_true, cos_theta_true, cos_theta_rec ); return a* b / (2*pi); } virtual double dP_dlogErec_dOmega_zenith_bins(int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true) { double bb = h_EERR.GetBin(bin_E_true, bin_cos_theta_true, bin_E_rec ); double a = h_EERR .GetBinContent(bb ) ; if (!isfinite(a)) fatal("h_errr failed: ", &h_EERR, bin_E_true, bin_cos_theta_true, bin_E_rec , a ); double b = h_COSZENREC.GetBinContent( bin_E_true, bin_cos_theta_true, bin_cos_theta_rec ); if (!isfinite(b)) fatal("h_cos failed: ", bin_E_true, bin_cos_theta_true, bin_cos_theta_rec); return a* b / (2*pi); } virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error_deg, double E_true, double cos_theta_true) { const double log10_reco_error_deg = log10( reco_error_deg ); double dP_dlog10alpha = interpolate( &h_ANGERR, E_true, cos_theta_true, log10_reco_error_deg ); double dP_dalpha = dP_dlog10alpha / ( reco_error_deg * log(10) ); double dP_dOmega = dP_dalpha / (2*pi * sin( reco_error_deg*pi/180) ); return interpolate( &h_EERR, E_true, cos_theta_true, E_rec ) * dP_dOmega; } virtual double dP_dlogErec_dOmega_reco_error_log( double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true) { return dP_dlogErec_dOmega_reco_error( E_rec, pow(10, log_reco_error_deg), E_true, cos_theta_true ); } // dP/dOmega = dP/(dcostheta * 2*pi) virtual double bg_dP_dlogErec_dOmega( double logE_rec, double cos_theta_rec ) { int bin = h_ATMNUREC.FindBin( logE_rec, cos_theta_rec ); return h_ATMNUREC.GetBinContent( bin ) / (2*pi); } virtual ~DefaultDetResponse() {} ClassDef( DefaultDetResponse, 1); }; #endif