#include "PointSourceComponent.hh" TH3D& PointSourceComponent::_pdf_erec_dayfraction_loga ( flavor flav, channel chan ) { static TTimeStamp T0 (2020, 1, 1, 0, 0, 0 ); { TH3D& r = HFC[flav][chan]; if ( r.GetEntries() > 0 ) return r; } print(name, "building lookup histogram.", flav, chan ); string nam = "HDR" + flavor_name(flav) + channel_name(chan); string tit = "pdf of erec, dayfraction, loga for " + flavor_name(flav) + channel_name(chan); // we use h_ANGERR to derive the binnin of the nes histogram DefaultDetResponse& detres = detector_response.get( flav, chan ); TAxis* XX = detres.h_ANGERR.GetXaxis(); TAxis* ZZ = detres.h_ANGERR.GetZaxis(); HFC[flav][chan] = TH3D( nam.c_str(), tit.c_str(), XX->GetNbins(), XX->GetXmin(), XX->GetXmax(), 48 , 0 , 1 , // dayfraction ZZ->GetNbins(), ZZ->GetXmin(), ZZ->GetXmax() ); TH3D& H = HFC[flav][chan]; H.SetXTitle("Erec (GeV)"); H.SetYTitle("Dayfraction"); H.SetZTitle("alpha (rad)"); // fill it. TTimeStamp T; for ( int bx = 1; bx <= H.GetNbinsX() ; bx++ ) // Erec slice { print ( " erec slice : ", bx, "/", H.GetNbinsX() ); for ( int by = 1; by <= H.GetNbinsY() ; by++ ) // time of day from 0 to 1 { double erec = H.GetXaxis()->GetBinCenter(bx); double dayf = H.GetYaxis()->GetBinCenter(by); // compute the true zenith angle for this time of day T.SetSec( T0.GetSec() + dayf * siderial_day() ) ; // todo : use period EquatorialCoords eq( 0 , dec.value ); // nb: ra set to zero on purpose Vec dtrue = -eq.track_direction( detres.det , T ); double theta_true = dtrue.theta(); if ( bx == 1 ) { print ("DAYF", dayf, theta_true ); } for ( int bz = 1; bz <= H.GetNbinsZ() ; bz++ ) // log angular error { double loga = H.GetZaxis()->GetBinCenter(bz); double J = detres.psf.d_loga_d_omega ( loga ); double alpha = pow(10, loga) * pi / 180 ; double v; if ( alpha > pi ) v = 0; else v = dN_dOmega_dlogErec_alpha( track, alpha , erec, theta_true , false ) / J ; if ( v < 0 ) { print ("warning : negative value", v, "for", flav, chan , "in N_dOmega_dlogErec_alpha at", bz, "apha=" , alpha, "erec=",erec, "costhetatrue=",cos(theta_true) , "J=", J ); v = 0; } int b = H.GetBin(bx, by, bz); H.SetBinContent( b, v ); } } } return H; } Table PointSourceComponent::rate_info( flavor flav, channel chan ) { Table T("Ebin","logE","Flux (/GeV m2 s)","Aeff@decl (m)","Rate (/yr)"); auto v = detector_response.Aeff(flav, chan ).GetAreaAtDecl( dec.value ).GetRateInfo( *flux ); double c =0; for(unsigned i = 0; i