#ifndef DETRESPONSE_INCLUDED_HH #define DETRESPONSE_INCLUDED_HH #include "rootutil.hh" #include "ana/search/EffectiveArea.hh" #include "TRandom.h" #include "PsfParam.hh" #include "definitions.hh" #include "histogram_conversion.hh" #include "statutil.hh" /*! 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("\n\nfailed to find object called", name, "\nin TFile", f.GetName(), "\n\n"); } dst = *((T *)p); dst.SetDirectory(0); } struct DetResponse { string name; channel chan; flavor flav; EffectiveArea Aeff; PsfParam psf; // smooth parametrisation. Det det; DetResponse() {} // 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); }; struct DefaultDetResponse : public DetResponse { string name; // detector resolution TH3D h_COSZENREC; // x = etrue, y=costhetatrue, z =costheta_rec TH3D h_EREC; // x = etrue, y=costheattrue, z =E_rec TH3D h_ANGERR; // x = etrue, y=costhetatrue, z =log(angle-muon-rec-track) TH3D h_EREC_sindecl; // like h_EREC bug with sindecl on the y-axis // expected number of events map h_BKGREC_3D; // x = erec y=costhetarec (Background contributions) z = azimuth TH3D h_BKGREC_azi_sindecl; map h_BKGREC; // x = erec y=costhetarec (Background contributions) TH2D h_BKGREC_sindecl; TH3D h_EFFAREA_3D; // Aeff vs energy vs zenith vs azimuth TH2D h_EFFAREA; // Aeff vs energy and zenith angle // maps to store histograms for NC and CC interactions ["CC/NC"] map h_COSZENREC_ALL; map h_EREC_ALL; map h_ANGERR_ALL; map h_BKGREC_3D_ALL; map h_BKGREC_ALL; map h_EFFAREA_3D_ALL; map h_EFFAREA_ALL; DefaultDetResponse() { init_det(); } DefaultDetResponse(string ingredients_file, flavor fl = numu, vector fl_int = {numuCC}, // you can do NC, CC or ALL here string rec_type = "track", double src_width = 0., int psf_mergebins = 2, bool gaus = true, double cos_zenith_cut = 999 ) { print ("DefaultDetResponse ctor"); init_det(); TFile f(ingredients_file.c_str()); // check if we have 2D or 3D histograms for EFFAREA and BKGREC TObject *tmp = f.Get(("h_EFFAREA_numuCC_" + rec_type).c_str()); if (!tmp) fatal("failed to get ","h_EFFAREA_numuCC_" + rec_type ) bool use_3D_hist = tmp->InheritsFrom("TH3"); print ("DefaultDetResponse 1"); if (fl != muon) { name = "DefaultDetResponse for " + flavor_name(fl) + "/" + rec_type + ", using 3D histograms " + use_3D_hist; // extract histograms from file and store for (auto fl_int_tmp : fl_int) { cout << "before read_histograms " << flavor_name(fl) << " " << flavor_interaction_name(fl_int_tmp) << endl; read_histograms(f, flavor_name(fl), fl_int_tmp, rec_type, use_3D_hist); } // add CC and NC for ALL, or copy histograms if only 1 interaction type add_histograms(flavor_name(fl), fl_int, rec_type); f.Close(); // zero out all bins in h_EFFARA wtih cos_zenith > cos_zenith_cut. // this to ensure misreconstructed events are not counted in ///the signal efficiency for (int bx = 1; bx <= h_EFFAREA.GetNbinsX(); bx++) { for (int by = 1; by <= h_EFFAREA.GetNbinsY(); by++) { double costheta = h_EFFAREA.GetYaxis()->GetBinCenter(by); int bin = h_EFFAREA.GetBin( bx, by ); if ( costheta > cos_zenith_cut ) h_EFFAREA.SetBinContent( bin , 0 ); } } // now make the Aeff object Aeff = EffectiveArea(h_EFFAREA, det); // 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_EREC, 1.0, true); // compute the energy response as function on sindecl print("computing h_EREC_sindecl", flavor_name(fl), rec_type); h_EREC_sindecl = convert_to_sindecl(det, h_EREC); print(h_EREC_sindecl.Integral()); check(h_COSZENREC); check(h_EREC); check(h_ANGERR); check(h_BKGREC["ATMNU_conv"]); check(h_BKGREC["ATMNU_prompt"]); if (psf_mergebins != 0) { std::cout << "initializing PSF" << std::flush; psf.name = "PSF for " + name; psf.init(h_ANGERR, psf_mergebins, src_width, gaus); } h_BKGREC["ATM_all"] = h_BKGREC["ATMNU_conv"]; h_BKGREC["ATM_all"].Add(&h_BKGREC["ATMNU_prompt"]); // h_BKGREC["ATM_all"].Add(&h_BKGREC["ATMMUON"]); } else // muon flavor { cout << "getting " << "h_ATMMU_reco_" << rec_type << endl; TObject *tmp = f.Get(("h_ATMMU_reco_" + rec_type).c_str()); bool use_3D_hist = tmp->InheritsFrom("TH3"); for (auto fl_int_tmp : fl_int) { read_histograms(f, flavor_name(fl), fl_int_tmp, rec_type, use_3D_hist); } add_histograms(flavor_name(fl), fl_int, rec_type); f.Close(); check(h_BKGREC["ATMMUON"]); h_BKGREC["ATM_all"] = h_BKGREC["ATMMUON"]; } // compute_background vs sindecl h_BKGREC_sindecl = convert_to_sindecl(det, h_BKGREC["ATM_all"]); cout << "sum background: " << h_BKGREC["ATM_all"].Integral() << endl; } void init_det() { det.longitude = 0.278819; det.latitude = 0.633407; det.meridian_convergence_angle = 0.0100733; } /*! Function reads the histograms from the ingredients file 3D histograms are projected to 2D histograms */ void read_histograms(TFile &f, string flavor, flavor_interaction fl_int, string rec_type, bool use_3D_hist) { string idstr = flavor_interaction_name(fl_int) + "_" + rec_type; print("Reading IRFs for flavor=", flavor, " channel/rectype = ", rec_type, flavor_interaction_name(fl_int)); if (fl_int != muonmuon) { string idstr = flavor_interaction_name(fl_int) + "_" + rec_type; read(f, h_COSZENREC_ALL[flavor_interaction_name(fl_int)], "h_COSZENREC_" + idstr); read(f, h_ANGERR_ALL[flavor_interaction_name(fl_int)], "h_ANGERR_" + idstr); read(f, h_EREC_ALL[flavor_interaction_name(fl_int)], "h_EREC_" + idstr); // extract 3D hists and store scaled 2D histograms if (use_3D_hist) { read(f, h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)], "h_EFFAREA_" + idstr); read(f, h_BKGREC_3D_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_conv_" + idstr); read(f, h_BKGREC_3D_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_prompt_" + idstr); // read( f, h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] , "h_ATMMU_reco_" + rec_type ); // convert to 2d histograms h_EFFAREA_ALL[flavor_interaction_name(fl_int)] = *((TH2D *)h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)].Project3D("yx")); h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)].Project3D("yx")); h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)].Project3D("yx")); // h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] = *((TH2D*) h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)].Project3D("yx") ); h_EFFAREA_ALL[flavor_interaction_name(fl_int)].Scale(1.0 / h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)].GetNbinsZ()); } // assume 2D histograms else { read(f, h_EFFAREA_ALL[flavor_interaction_name(fl_int)], "h_EFFAREA_" + idstr); read(f, h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_conv_" + idstr); read(f, h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_prompt_" + idstr); // read( f, h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] , "h_ATMMU_reco_" + rec_type ); } } else { if (use_3D_hist) { read(f, h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)], "h_ATMMU_reco_" + rec_type); h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)].Project3D("yx")); } else { read(f, h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)], "h_ATMMU_reco_" + rec_type); } } } /*! Function adds interaction contributions to one histogram */ void add_histograms(string flavor, vector interactions, string rec_type) { cout << "adding histograms" << endl; cout << "interaction " << flavor_interaction_name(interactions[0]) << endl; string idstr = flavor + "_" + rec_type; if (startswith(flavor, "a") || startswith(flavor, "n")) { h_COSZENREC = *((TH3D *)h_COSZENREC_ALL[flavor_interaction_name(interactions[0])].Clone(("h_COSZENREC_ALL_" + idstr).c_str())); h_EREC = *((TH3D *)h_EREC_ALL[flavor_interaction_name(interactions[0])].Clone(("h_EREC_ALL_" + idstr).c_str())); h_ANGERR = *((TH3D *)h_ANGERR_ALL[flavor_interaction_name(interactions[0])].Clone(("h_ANGERR_ALL_" + idstr).c_str())); h_EFFAREA = *((TH2D *)h_EFFAREA_ALL[flavor_interaction_name(interactions[0])].Clone(("h_EFFAREA_ALL_" + idstr).c_str())); h_BKGREC["ATMNU_conv"] = *((TH2D *)h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMNU_reco_conv_ALL_" + idstr).c_str())); h_BKGREC["ATMNU_prompt"] = *((TH2D *)h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMNU_reco_prompt_ALL_" + idstr).c_str())); if (interactions.size() > 1) { for (unsigned i = 1; i < interactions.size(); i++) { cout << "interaction " << flavor_interaction_name(interactions[i]) << endl; h_COSZENREC.Add(&h_COSZENREC_ALL[flavor_interaction_name(interactions[i])]); h_EREC.Add(&h_EREC_ALL[flavor_interaction_name(interactions[i])]); h_ANGERR.Add(&h_ANGERR_ALL[flavor_interaction_name(interactions[i])]); h_EFFAREA.Add(&h_EFFAREA_ALL[flavor_interaction_name(interactions[i])]); h_BKGREC["ATMNU_conv"].Add(&h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(interactions[i])]); h_BKGREC["ATMNU_prompt"].Add(&h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(interactions[i])]); } } } else if (startswith(flavor, "m")) { h_BKGREC["ATMMUON"] = *((TH2D *)h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMMUON_" + idstr).c_str())); // make sure to not double add the muons! } else { fatal(" Error: flavour ", flavor, " not implemented"); } } /*! 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) { dbg(""); double a = interpolate(&h_EREC, E_true, cos_theta_true, E_rec, true); double b = interpolate(&h_COSZENREC, E_true, cos_theta_true, cos_theta_rec, true); return a * b / (2 * pi); } /*! P(E_rec, cos_theta_rec | E_true, cos_theta_true ) */ 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_EREC.GetBin(bin_E_true, bin_cos_theta_true, bin_E_rec); double a = h_EREC.GetBinContent(bb); if (!isfinite(a)) fatal("h_errr failed: ", &h_EREC, 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); } /*! P ( E_rec, rec_angle_from_true | E_true, cos_theta ) */ // virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error_deg, double E_true, double cos_theta_true) // { // return dP_dlogErec_dOmega_reco_error_log( E_rec, log10(reco_error_deg), E_true, cos_theta_true ); // } /*! P ( E_rec, log10(rec_angle_from_true) | E_true, cos_theta ) */ virtual double dP_dlogErec_dOmega_reco_error_log(double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true) { // MEMCHECK // print( E_rec, log_reco_error_deg, E_true, cos_theta_true ); double dP_dOmega = psf.dP_dOmega(E_true, log_reco_error_deg); double r = interpolate(&h_EREC, E_true, cos_theta_true, E_rec, true) * dP_dOmega; if (::isnan(r)) { fatal("nan detected in dP_dlogErec_dOmega_reco_error_log", E_true, cos_theta_true, E_rec, dP_dOmega); } return r; } /*! for a given true energy and declination, return the fraction of events with Erec > cutoff_erec */ double erec_fraction_above_sindecl(double loge, double sindecl, double cutoff_logerec) { int bx = h_EREC_sindecl.GetXaxis()->FindBin(loge); int by = h_EREC_sindecl.GetYaxis()->FindBin(sindecl); double n(0), m(0); for (int bz = 1; bz <= h_EREC_sindecl.GetNbinsZ(); bz++) { // double logerec = h_EREC_sindecl.GetZaxis()->GetBinCenter( bz ); double xlow = h_EREC_sindecl.GetZaxis()->GetBinLowEdge(bz); double xup = h_EREC_sindecl.GetZaxis()->GetBinUpEdge(bz); double width = h_EREC_sindecl.GetZaxis()->GetBinWidth(bz); int bin = h_EREC_sindecl.GetBin(bx, by, bz); double c = h_EREC_sindecl.GetBinContent(bin); if (xlow > cutoff_logerec) n += c; // the whole bin is above the cut else if (xup > cutoff_logerec) n += c * (xup - cutoff_logerec) / width; m += c; // print ("erec/n/m", bx,by,bz, bin, logerec, n, m, sindecl, h_EREC_sindecl.Integral() ); } if (m == 0) return 0; return n / m; } /*! for a given true energy and declination, return the fraction of events with Erec within an energy bin */ double erec_fraction_range_above_sindecl(double loge, double sindecl, double cutoff_low_logerec, double cutoff_high_logerec) { int bx = h_EREC_sindecl.GetXaxis()->FindBin(loge); int by = h_EREC_sindecl.GetYaxis()->FindBin(sindecl); double n(0), m(0); for (int bz = 1; bz <= h_EREC_sindecl.GetNbinsZ(); bz++) { double xlow = h_EREC_sindecl.GetZaxis()->GetBinLowEdge(bz); double xup = h_EREC_sindecl.GetZaxis()->GetBinUpEdge(bz); int bin = h_EREC_sindecl.GetBin(bx, by, bz); double c = h_EREC_sindecl.GetBinContent(bin); // which fraction of our bin overlaps with the requested enery interval? Interval ival( xlow, xup ); // from statutil.hh double weight = ival.fraction_in_range( cutoff_low_logerec, cutoff_high_logerec ); if (weight <0 || weight > 1 ) fatal ("programming error weight", weight); n += weight * c; m += c; } if (m == 0) return 0; return n / m; } /*! The number of background events in a cone of size cutoff_angle and with * Erec > cutoff */ double background_rate(double livetime, double declination, double cutoff_angle, double cutoff_logerec) { int by = h_BKGREC_sindecl.GetYaxis()->FindBin(sin(declination)); double r(0), f(0); int firstbin = h_BKGREC_sindecl.GetXaxis()->FindBin(cutoff_logerec); for (int bx = firstbin; bx <= h_BKGREC_sindecl.GetNbinsX(); bx++) { if (bx == firstbin) f = 1 - (cutoff_logerec - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx); else f = 1; // print("f=",f); r += f * h_BKGREC_sindecl.GetBinContent(bx, by); } double g = 4 * pi * pow(sin(cutoff_angle / 2), 2) / (2 * pi * h_BKGREC_sindecl.GetYaxis()->GetBinWidth(by)); return r * g * livetime / (365.0 * 24 * 60 * 60); } double background_rate_range(double livetime, double declination, double angle_min, double angle_max, double logerec_min, double logerec_max) { int by = h_BKGREC_sindecl.GetYaxis()->FindBin(sin(declination)); double r(0), f(0); int firstbin = h_BKGREC_sindecl.GetXaxis()->FindBin(logerec_min); int lastbin = h_BKGREC_sindecl.GetXaxis()->FindBin(logerec_max); if (firstbin == lastbin) { // print("oops"); f = (logerec_max - logerec_min) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(firstbin); r += f * h_BKGREC_sindecl.GetBinContent(firstbin, by); } else { for (int bx = firstbin; bx <= lastbin; bx++) { if (bx == firstbin) f = 1 - (logerec_min - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx); else if (bx == lastbin) f = (logerec_max - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx); else f = 1; // print("bx", bx ,"f=",f, "bincont", h_BKGREC_sindecl.GetBinContent( bx, by ), "r", r); r += f * h_BKGREC_sindecl.GetBinContent(bx, by); } } double g = (4 * pi * pow(sin(angle_max / 2), 2) - 4 * pi * pow(sin(angle_min / 2), 2)) / (2 * pi * h_BKGREC_sindecl.GetYaxis()->GetBinWidth(by)); return r * g * livetime / (365.0 * 24 * 60 * 60); } /*! return a human-readable table with rate information for flux */ Table compute_rates_table(Flux &flux, double livetime, // seconds bool point, double declination, double cutoff_angle, double cutoff_logerec); // in cc /*! Return total signal and background, optionally filling table */ vector compute_rates(Flux &flux, double livetime, // seconds bool point, double declination, double cutoff_angle, double cutoff_logerec, Table *T = nullptr); // in cc /*! Return total signal and background in energy/angle bin */ vector compute_rates_angle_energy_range(Flux &flux, double livetime, // seconds bool point, double declination, double angle_min, double angle_max, double logerec_min, double logerec_max); // in cc virtual ~DefaultDetResponse() { // print ("dtor of DefaultDetResponse"); } ClassDef(DefaultDetResponse, 1); }; /*! A collection of DetResponse's for each channel and neutrino flavor */ struct FullDetResponse { string name; // vector flavors={ muon, nue, anue, numu, anumu, nutau, anutau}; vector flavors; // let's build it from fl_int // we decide which ingredients are actually used by the list provided by user in fl_int vector channels; vector fl_int; string filename; map> response_map; // by pdg-type of incomming neutrino and by observed channeld (track, shower, ...) FullDetResponse() {} /*! Initialize from File. \param cos_zenith_cut : for bins above this cut (downgong = 1, upgoing = -1), the Aeff will be set to zero */ FullDetResponse(string name, string filename_, vector fl_int = {numuCC, numuNC, nueCC, nueNC, nutauNC, nutauCCshowerdecay, nutauCCmuondecay, anumuCC, anumuNC, anueCC, anueGLRES, anueNC, anutauNC, anutauCCshowerdecay, anutauCCmuondecay, muonmuon}, vector channels = {track, shower}, double src_width = 0., bool gaus = true, double cos_zenith_cut = 999 ) : channels(channels) { flavors = get_flavor_from_flavor_interaction(fl_int); filename = filename_; int psf_mergebins = 2; for (auto chan : channels) for (auto flav : flavors) { if (get_flavor_interaction(flav, fl_int).size() > 0) { cout << "DEBUG in Fulldetresp " << flavor_name(flav) << " " << channel_name(chan) << endl; response_map[flav][chan] = DefaultDetResponse(filename, flav, get_flavor_interaction(flav, fl_int), channel_name(chan), src_width, psf_mergebins, gaus, cos_zenith_cut ); response_map[flav][chan].flav = flav; response_map[flav][chan].chan = chan; response_map[flav][chan].Aeff.flav = flav; response_map[flav][chan].Aeff.chan = chan; } } } Det &det() { return (response_map.begin()->second.begin()->second.det); } DefaultDetResponse &get(flavor nutype, channel chan) { if (response_map.find(nutype) == response_map.end()) { fatal("FullDetResponse object has no DetResponse for flavor", nutype); } if (response_map[nutype].find(chan) == response_map[nutype].end()) { fatal("FullDetResponse for flavor", nutype, "has no channel", chan); } return response_map[nutype][chan]; } double dP_dlogErec_dOmega_zenith(flavor nutype, channel chan, double E_rec, double cos_theta_rec, double E_true, double cos_theta_true) { return get(nutype, chan).dP_dlogErec_dOmega_zenith(E_rec, cos_theta_rec, E_true, cos_theta_true); } double dP_dlogErec_dOmega_reco_error(flavor nutype, channel chan, double E_rec, double reco_error, double E_true, double cos_theta_true) { return get(nutype, chan).dP_dlogErec_dOmega_reco_error(E_rec, reco_error, E_true, cos_theta_true); } double dP_dlogErec_dOmega_zenith_bins(flavor nutype, channel chan, int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true) { return get(nutype, chan).dP_dlogErec_dOmega_zenith_bins(bin_E_rec, bin_cos_theta_rec, bin_E_true, bin_cos_theta_true); } double dP_dlogErec_dOmega_reco_error_log(flavor nutype, channel chan, double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true) { return get(nutype, chan).dP_dlogErec_dOmega_reco_error_log(E_rec, log_reco_error_deg, E_true, cos_theta_true); } EffectiveArea &Aeff(flavor nutype, channel chan) { return get(nutype, chan).Aeff; } void show_rates(Flux &flux, double livetime, bool point, double declination, double cutoff_angle, double cutoff_logerec); // void compute_rates( Flux& flux, bool point, double declination = -99 , double cutoff_angle = pi , double cutoff_erec = 0 ); // in cc void compute_rates_summary(Flux &flux, bool point, double declination = -99, double max_angle = pi); // in cc virtual ~FullDetResponse(){}; ClassDef(FullDetResponse, 1); }; #endif