#ifndef BACKGROUNDCOMPONENT_HH #define BACKGROUNDCOMPONENT_HH #include "Component.hh" #include "dayfraction.hh" #include "Mat.hh" /*! A background-component is meant for e.g. atmspheric neutrino background. Such a contribution is typically taking from data, with maybe the normalisation a free parameter. */ struct BackgroundComponent : public Component { FullDetResponse& detres; map< channel, TH2D> hbkg; // histogram of all the background for this channel (sum over flavors) map< channel, TH1D> hbkg_erec; // 1d pdf of erec only. map< channel, TH3D > _pdf_erec_dayfraction_decl; // pdf of Erec and time-of-day vs decl map< channel, vector > _pdf_erec_dayfraction_decl_slices; map< channel, TH1D > _dP_dOmega_decl; // number of events found at decl (modulo norm) map< channel, TH2D > _dP_dOmega_dlogerec_decl; // integrated over time. vector _dp_domega_dloge; // search-event cache vector< channel> channels; BackgroundComponent( string name, FullDetResponse& detr ) : Component(name), detres(detr) { // add up all the background channels norm.value = 0; for ( auto chan : detr.channels ) { for (auto flav : detr.flavors ) { if( startswith(flavor_name(flav),"a") || startswith(flavor_name(flav),"n" )){ if ( hbkg.find( chan ) == hbkg.end() ) { hbkg[chan] = detr.get( flav, chan ).h_BKGREC["ATMNU_conv"]; // copy } else { hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMNU_conv"] ); } hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMNU_prompt"] ); } else { if ( hbkg.find( chan ) == hbkg.end() ) { hbkg[chan] = detr.get( flav, chan ).h_BKGREC["ATMMUON"]; } else { hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMMUON"] ); } } } norm.value += hbkg[chan].Integral(); normalize( hbkg[chan] ); init_cos_zenith_decl(chan); } norm.fit_max = norm.value * 1.3; norm.fit_min = norm.value * 0.7; info(); } void init_cos_zenith_decl(channel chan) { _pdf_erec_dayfraction_decl[chan] = TH3D("_pdf_erec_dayfraction_decl", "_pdf_erec_dayfraction_decl", hbkg[chan].GetNbinsX(), hbkg[chan].GetXaxis()->GetXmin(), hbkg[chan].GetXaxis()->GetXmax(), 100,0,1, 40,-1,1); //_rate_at_sindecl = TH1D( "_rate_at_sindecl","background rate_vs_decl", 40,-1,1 ); print ("initializing background vs declination histogram"); TH3D& H = _pdf_erec_dayfraction_decl[chan]; for (int bx =1 ; bx <= H.GetXaxis()->GetNbins() ; bx ++ ) { double logerec = H.GetXaxis()->GetBinCenter(bx); print (bx, "/", H.GetXaxis()->GetNbins() ); for (int bz =1 ; bz <= H.GetZaxis()->GetNbins() ; bz ++ ) { double sindecl = H.GetZaxis()->GetBinCenter( bz ); for (int by =1 ; by <= H.GetYaxis()->GetNbins() ; by ++ ) { double dayf = H.GetYaxis()->GetBinCenter( by ); double costheta = dayfraction_to_coszenith( 0, asin(sindecl), detres.det(), dayf ); // ra = 0 on purpose double v = dP_dOmega_dlogE( chan, logerec, costheta ); int bin = H.GetBin( bx, by, bz ); H.SetBinContent( bin , v / H.GetYaxis()->GetNbins() ); } // by }//bz }//bx _pdf_erec_dayfraction_decl_slices[chan] = split_z( _pdf_erec_dayfraction_decl[chan] ); _dP_dOmega_dlogerec_decl[chan] = * (( TH2D*) H.Project3D("ZX")); //_dP_dOmega_dlogerec_decl[chan].Scale( H.GetXaxis()->GetBinWidth(1) ); _dP_dOmega_decl[chan] = *(H.ProjectionZ()); _dP_dOmega_decl[chan].Scale( H.GetXaxis()->GetBinWidth(1) ); } TH2D& get_dayf_logerec_hist( channel chan, double sindecl ) { int b = _pdf_erec_dayfraction_decl[chan].GetZaxis()->FindBin( sindecl); return _pdf_erec_dayfraction_decl_slices[chan][b]; } double dP_dOmega_dlogErec_at_declination( channel chan, double logerec, double sindecl ) { return _dP_dOmega_dlogerec_decl[chan].Interpolate( logerec, sindecl ); } double dN_dOmega_dlogErec_at_declination( channel chan, double logerec, double sindecl ) { return norm.value * dP_dOmega_dlogErec_at_declination( chan, logerec, sindecl ); } /*! fill the vector _dp_domega_dloge with dP_dOmega_dlogE, which does not depend on any of the Params */ virtual void init( DataSet& dset ) { _dp_domega_dloge.resize( dset.data.size() ); for (unsigned i = 0 ; i < dset.data.size(); i++ ) { auto& sev = dset.data[i]; if ( sev.index != i ) fatal ("dataset indices not correct"); _dp_domega_dloge[i] = dP_dOmega_dlogE( (channel) sev.chan, sev.Eproxy, cos( sev.zenith ) ); } } virtual void deinit() { _dp_domega_dloge.resize(0); } virtual void update() { // nothing to do: only norm can change. } virtual double dN_dOmega_dlogE(SearchEvent &sev ) { TIMEFUNC double p = 0; if (_dp_domega_dloge.size() == 0 ) { p = dP_dOmega_dlogE( (channel) sev.chan, sev.Eproxy, cos( sev.zenith ) ); } else p = _dp_domega_dloge[ sev.index ]; double r = norm.value * p; if (r<0) fatal( "illegal value", r, sev.str() ); return r; } double dP_dOmega_dlogE( channel chan, double E_rec, double cos_zenith_rec ) { return get_hist(chan).Interpolate( E_rec, cos_zenith_rec ) / (2*pi); } double check_integral( channel chan ) { TH2D& h = get_hist(chan); double n = 0; for ( int bx = 1 ; bx <= h.GetXaxis()->GetNbins() ; bx++ ) { double dloge = h.GetXaxis()->GetBinWidth(bx); for ( int by = 1 ; by <= h.GetYaxis()->GetNbins() ; by++ ) { double logerc = h.GetXaxis()->GetBinCenter(bx); double coszen = h.GetYaxis()->GetBinCenter(by); double sr = h.GetYaxis()->GetBinWidth(by) * 2*pi; // solid angle in this bin n += sr * dloge * dP_dOmega_dlogE( chan, logerc, coszen ); } } return n; } double check_integral_declination( channel chan ) { TH2D& h = get_hist(chan); double n = 0; for ( int bx = 1 ; bx <= h.GetXaxis()->GetNbins() ; bx++ ) { double dloge = h.GetXaxis()->GetBinWidth(bx); for ( int by = 1 ; by <= h.GetYaxis()->GetNbins() ; by++ ) { double logerc = h.GetXaxis()->GetBinCenter(bx); double sindecl = h.GetYaxis()->GetBinCenter(by); double sr = h.GetYaxis()->GetBinWidth(by) * 2*pi; // solid angle in this bin n += sr * dloge * dP_dOmega_dlogErec_at_declination( chan, logerc, sindecl ); } } return n; } TH2D& get_hist( channel chan ) { auto c = hbkg.find( chan ); if ( c== hbkg.end() ) { fatal ("background called for a channel it does not have ", chan ); } return (*c).second; } virtual SearchEvent generate_event() // nb : only track for the moment { TIMEFUNC double z, E; channel chan = track; hbkg[chan].GetRandom2( E, z ); SearchEvent r( acos(z), 2*pi*gRandom->Rndm(), period.random(), E ); r.compute_coordinates( detres.det() ); r.chan = chan; //dbg ( r.chan , r.__str__() ); return r; } /*! version that generates background events in a small region of the * celestial sky. */ vector generate_events_region( unsigned howmany, channel chan, EquatorialCoords center, double radius ) { vector v(howmany); double sindecl = sin ( center.dec() ); TH2D& h = get_dayf_logerec_hist( chan, sindecl ); double dayf; Mat R( center.asVec() ); for (unsigned i =0; iRndm() ) *radius; double phi = 2* pi * gRandom->Rndm(); Vec v; v.set_angles( x, phi ); r.coordinates.fromVec( R* v ); r.t = dayfraction_to_time( dayf ); r.compute_source_direction( detres.det() ); // lets just check the pdf is not zero double p = dN_dOmega_dlogE( r ); if ( p <= 0 ) fatal(" background generation produced an impossible event "); } return v; } void info() { Table T ("channel","binsx","binsy","integral"); T.title = "background-component info"; for( auto p : hbkg ) { auto& h = p.second; auto x = p.first; T << int(x) << h.GetNbinsX() << h.GetNbinsY() << h.Integral(); } print (T); } // -----boilerplate----- virtual BackgroundComponent *clone() const { return new BackgroundComponent(*this); } virtual ~BackgroundComponent() {} ClassDef( BackgroundComponent, 1) }; #endif