#ifndef EFFECTIVEAREA_INCLUDED_HH #define EFFECTIVEAREA_INCLUDED_HH #include "io_util.hh" #include "TH2D.h" #include "TMath.h" #include "Flux.hh" #include "Astro.hh" #include "definitions.hh" /*! Effective area vs energy */ struct EffectiveArea1D { TH1D H; flavor flav; channel chan; EffectiveArea1D() {} EffectiveArea1D( TH1D h , flavor flav, channel chan ) : H(h), flav(flav), chan(chan) {} /*! Return the amount of GeV's in a logarithmic bin */ double GetWidthE (int i) { double Emin = pow(10.,H.GetXaxis()->GetBinLowEdge(i)); double Emax = pow(10.,H.GetXaxis()->GetBinUpEdge(i)); return Emax-Emin; } void prnt() { print ( " H has ", H.GetNbinsX(), H.Integral() ); } /*! Integrate a function(logE) (e.g. a flux) over the effective area. Func must have units GeV-1 */ template< typename F > double integrate ( const F& func ) { double c = 0; for (int bx = 1; bx <= H.GetNbinsX() ; bx++ ) { double aeff = H.GetBinContent(bx); double loge = H.GetBinCenter(bx); double width = GetWidthE( bx ); double f = func( loge ); //print ( aeff, loge, width, f ); c += aeff * width * f; } return c; } /*! Integrate a Aeff* flux to get the total event rate. The flux should provide the method dNdE(flavour,log_E) */ double GetTotalRate( Flux& F, double livetime = 365*24*3600 ) { double r = integrate( [&](double loge) { return F.dNdE( flav, loge ); } ) * livetime ; return r; } /*! Return a histogratm of the Rate per energy bin. */ TH1D GetRate( Flux& flux , double livetime = 365*24*3600 ) { TH1D r = H; r.SetName("r"); for (int bx = 1; bx <= r.GetNbinsX() ; bx++ ) { double aeff = H.GetBinContent(bx); double loge = r.GetXaxis()->GetBinCenter(bx); double f = flux.dNdE( flav, loge ); //print( "xxx", f, aeff, GetWidthE(bx) , livetime); r.SetBinContent(bx, aeff * f * GetWidthE(bx) * livetime ); } return r; } /*! Version of GetRate that returns 4 vectors: bin-center Energy, Flux, Aeff, Rate */ vector< vector< double > > GetRateInfo( Flux& flux, int nseconds = 365*24*3600 ) { vector< vector< double >> r; r.resize(4); TH1D h = GetRate( flux, nseconds ); for (int bx = 1; bx <= h.GetNbinsX() ; bx++ ) { r[0].push_back( h.GetBinCenter(bx) ); r[1].push_back( flux.dNdE( flav, h.GetBinCenter(bx) )); r[2].push_back( H.GetBinContent(bx) ); // big H = Aeff r[3].push_back( h.GetBinContent(bx) ); // rate } return r; } ~EffectiveArea1D() {} }; /*! Effective Area vs E_true and costheta_true */ struct EffectiveArea { TH2D Haeff_costheta; TH2D Haeff_sindecl; bool Haeff_sindecl_init = false; vector< EffectiveArea1D > Haeff_slices_sindecl; // the 1-d slices of Haeff_sindecl bool GetAreaAtDecl_init = false ; Det det; // ugly but.... todo: make one global detector flavor flav; channel chan; EffectiveArea () { //print("warning default Aeff ctor" ); det.latitude = -100; } // EffectiveArea(TH2D H_costheta, TH2D H_sindecl) : // Haeff_costheta(H_costheta),Haeff_sindecl(H_sindecl) {} EffectiveArea(TH2D H_costheta, Det& det): Haeff_costheta(H_costheta), det(det) { if (det.latitude < 0.4 || det.latitude > 0.7 ) fatal ("wrong lattitude"); } /* compute Aeff vs sin(declination) */ void ensure_Haeff_sindecl(); // in cc //functions to extract information from Haeff and Haeff_sindec double GetPhaseSpace_costheta (int i) { return GetWidthE_costheta(i) * 2*pi * Haeff_costheta.GetYaxis()->GetBinWidth(1); } double GetPhaseSpace_sindecl (int i) { ensure_Haeff_sindecl(); return GetWidthE_sindecl(i) * 2*pi * Haeff_sindecl.GetYaxis()->GetBinWidth(1); } double GetWidthE_costheta (int i) { double Emin = TMath::Power(10.,Haeff_costheta.GetXaxis()->GetBinLowEdge(i)); double Emax = TMath::Power(10.,Haeff_costheta.GetXaxis()->GetBinUpEdge(i)); return Emax-Emin; } double GetWidthE_sindecl (int i) { ensure_Haeff_sindecl(); double Emin = TMath::Power(10.,Haeff_sindecl.GetXaxis()->GetBinLowEdge(i)); double Emax = TMath::Power(10.,Haeff_sindecl.GetXaxis()->GetBinUpEdge(i)); return Emax-Emin; } /*! Integrate over energy at the specified cos(zenith). Func must have units GeV-1 */ template< typename F> double integrate ( double cos_zenith, F func ) { double c = 0; int by = Haeff_costheta.GetYaxis()->FindBin( cos_zenith ); for (int bx = 1; bx <= Haeff_costheta.GetNbinsX() ; bx++ ) { double aeff = Haeff_costheta.GetBinContent(bx,by); double loge = Haeff_costheta.GetXaxis()->GetBinCenter(bx); double width = GetWidthE_costheta( bx ); c += aeff * width * func( loge ); } return c; } /*! return the total rate. The function will call the flux object with the correct flavor */ double rate( double cos_zenith, Flux* PF ) { return integrate( cos_zenith, [&] (double logEtrue ) { return PF->dNdE( flav , logEtrue); } ); } /*! return the rate (th2d) for the given flux and the given time-span */ TH2D GetRate( DiffuseFlux* F , int nseconds = 365*24*3600 ) { TH2D r = Haeff_costheta; r.SetName("rate"); for (int bx = 1; bx <= Haeff_costheta.GetNbinsX(); bx++ ) { double loge = Haeff_costheta.GetXaxis()->GetBinCenter(bx); for (int by = 1; by <= Haeff_costheta.GetNbinsY(); by++ ) { double costheta = Haeff_costheta.GetYaxis()->GetBinCenter(by); double flux = F->dNdEdOmega(flav, loge, costheta ); double aeff = Haeff_costheta.GetBinContent(bx,by); r.SetBinContent( bx, by, aeff * flux * GetPhaseSpace_costheta(bx) ); } } r.Scale(nseconds); return r; } /*! return the total rate for the given flux and the given time-span */ double GetTotalRate( DiffuseFlux* F , int nseconds = 365*24*3600 ) { double rate = 0; for (int bx = 1; bx <= Haeff_costheta.GetNbinsX(); bx++ ) { double loge = Haeff_costheta.GetXaxis()->GetBinCenter(bx); for (int by = 1; by <= Haeff_costheta.GetNbinsY(); by++ ) { double costheta = Haeff_costheta.GetYaxis()->GetBinCenter(by); double aeff = Haeff_costheta.GetBinContent(bx,by); double flux = F->dNdEdOmega(flav, loge, costheta ); rate += aeff * flux * GetPhaseSpace_costheta(bx); } } return rate * nseconds; } EffectiveArea1D GetAreaAtDecl( double decl ) { if (!GetAreaAtDecl_init) { ensure_Haeff_sindecl(); int n = Haeff_sindecl.GetYaxis()->GetNbins(); Haeff_slices_sindecl.resize( n + 1 ); print ("initializing GetAreaAtDecl for",n,"slices", flav, chan ); for (int b = 1 ; b <= n ; b++ ) { TH1D* h = Haeff_sindecl.ProjectionX( ("hsindeclslice"+str(b)+str(chan)+str(flav)).c_str(),b,b); h->SetDirectory(0); Haeff_slices_sindecl[b] = EffectiveArea1D( *h, flav, chan ); } GetAreaAtDecl_init = true; } int b = Haeff_sindecl.GetYaxis()->FindBin( sin(decl) ); return Haeff_slices_sindecl[b]; } EffectiveArea1D GetAverage() { print ("EffectiveArea1D GetAverage() -- warning: this leaks "); // don't remove print before fixing the leak, please EffectiveArea1D r( *Haeff_costheta.ProjectionX("aeff_average"), flav, chan ); r.H.Scale( 1.0 / Haeff_costheta.GetNbinsY() ); return r; } //! return total rate of point source at declination= decl double GetRateDecl( double decl, Flux* F, int nseconds = 365*24*3600 ) { ensure_Haeff_sindecl(); double sd = sin(decl); double r =0; int by = Haeff_sindecl.GetYaxis()->FindBin( sd ); for (int bx = 1; bx <= Haeff_sindecl.GetNbinsX() ; bx++ ) { double loge = Haeff_sindecl.GetXaxis()->GetBinCenter(bx); double flux = F->dNdE( flav , loge ); double aeff = Haeff_sindecl.GetBinContent(bx,by); r += flux * aeff * GetWidthE_sindecl( bx ); } return r * nseconds; } virtual EffectiveArea *clone() const { return new EffectiveArea(*this); } virtual ~EffectiveArea() {} ClassDef( EffectiveArea, 1) }; #endif