#ifndef HISTOCONV_INCD #define HISTOCONV_INCD #include "Det.hh" #include "foreach.hh" #include "rootutil.hh" #include "TH2D.h" struct convert_decl { TH2D pivot_matrix; }; /*! for every sin-decl, you have a slice of cos-theta with the bin-content * indication how much time the source spends there */ inline TH2D conversion_table_costheta_sindecl( Det& det, int nbinsxy, int Nsamples ) { static TH2D R; if ( R.GetNbinsY() == nbinsxy ) return R; // todo: check the rest of the parameters print ("generating new conversion table for nbins=", nbinsxy ); R = TH2D("h_conversion_table","conversion from costheta_to_sindecl", nbinsxy, -1,1, nbinsxy, -1,1 ); TTimeStamp T(2020, 1, 1, 2, 2, 2, 2 ); for (int bx = 1; bx <= R.GetNbinsX() ; bx++ ) { double sindecl = R.GetXaxis()->GetBinCenter(bx); double decl = asin( sindecl ); for ( int i = 0; i < Nsamples ; i++ ) { double ra = i * 2 * pi / Nsamples; EquatorialCoords coords( ra, decl ); Vec v = -coords.track_direction( det, T ); double costheta = v.z; R.Fill( sindecl,costheta, 1.0/Nsamples); } } return R; } /*! take a 2d histogram with cos-zenith on the y-axis and turn it into a histogram * with sin-decl on the y-axis * */ inline TH2D* h_declzen_conversion( TH2D* p = nullptr ) { static TH2D* x; if ( p ) x = p; return x; } inline TH2D convert_to_sindecl( Det& det, TH2D& H ) { static TH2D C; if ( H.GetNbinsY() != C.GetNbinsY() ) { print ( "gen new conv table", H.GetNbinsY() , C.GetNbinsY()); C = conversion_table_costheta_sindecl( det, H.GetNbinsY(), 1000 ); h_declzen_conversion(&C); } TH2D R = H; // copy R.Reset(); R.SetName( ( string( H.GetName()) + "__sindecl" ).c_str() ); R.GetYaxis()->SetTitle("sin(decl)"); for (int hy = 1; hy <= H.GetNbinsY(); hy++ ) { for (int ry = 1; ry <= R.GetNbinsY(); ry++ ) { double w = C.GetBinContent( ry, hy ); for (int hx=1; hx <= H.GetNbinsX(); hx++ ) { double u = R.GetBinContent( hx, ry ); R.SetBinContent( hx, ry, u + w * H.GetBinContent( hx, hy ) ); } } } return R; } /*! Given a 3d histogram who's y-axis contains cos-zenith, Return one * with y-axis representing sin(decl) */ inline TH3D convert_to_sindecl( Det& det, TH3D& H ) { const TH2D& C = conversion_table_costheta_sindecl( det, H.GetNbinsY(), 1000 ); TH3D R = H; // copy R.Reset(); R.SetName( ( string( H.GetName()) + "__sindecl" ).c_str() ); R.GetYaxis()->SetTitle("sin(decl)"); if ( H.Integral() == 0 ) fatal("empty histogram in convert_to_sindecl"); for (int hy = 1; hy <= H.GetNbinsY(); hy++ ) { for (int ry = 1; ry <= R.GetNbinsY(); ry++ ) { double w = C.GetBinContent( ry, hy ); if ( w == 0 ) { continue; } for (int hz=1; hz <= H.GetNbinsZ(); hz++ ) { for (int hx=1; hx <= H.GetNbinsX(); hx++ ) { double u = R.GetBinContent( hx, ry , hz ); R.SetBinContent( hx, ry, hz, u + w * H.GetBinContent( hx, hy, hz ) ); } } } } return R; } #endif