diff --git a/inc/VPlotAnasumHistograms.h b/inc/VPlotAnasumHistograms.h index ed53dbcd..fa9ce5e4 100644 --- a/inc/VPlotAnasumHistograms.h +++ b/inc/VPlotAnasumHistograms.h @@ -101,7 +101,8 @@ class VPlotAnasumHistograms : public VAnalysisUtilities, public VPlotUtilities, void plot_mscPlots( int irebin = 2, double xmin = -2., double xmax = 4., string mscwfile = "" ); void plot_qualityHistograms( double iSourceStrength = 1., bool bUpper = true, int iMethod = 0 ); TCanvas* plot_theta2( double t2min = 0., double t2max = 0.3, int irbin = 5, - double setYMax = -9999., double setYMin = -9999. ); + double setYMax = -9999., double setYMin = -9999., + bool containment_lines = false ); void plot_theta2Correction(); void plot_UncorrelatedSkyPlots(); void plot_CorrelatedSkyPlots(); @@ -144,7 +145,7 @@ class VPlotAnasumHistograms : public VAnalysisUtilities, public VPlotUtilities, } bool setRunNumber( int iRun ); // select run for plotting - ClassDef( VPlotAnasumHistograms, 16 ); + ClassDef( VPlotAnasumHistograms, 17 ); }; #endif diff --git a/inc/VPlotInstrumentResponseFunction.h b/inc/VPlotInstrumentResponseFunction.h index 9df8b45d..4f4246c6 100644 --- a/inc/VPlotInstrumentResponseFunction.h +++ b/inc/VPlotInstrumentResponseFunction.h @@ -19,6 +19,7 @@ #include "VHistogramUtilities.h" #include "VInstrumentResponseFunctionReader.h" #include "VPlotUtilities.h" +#include "VUtilities.h" using namespace std; diff --git a/inc/VUtilities.h b/inc/VUtilities.h index 2eac76bd..a9986727 100644 --- a/inc/VUtilities.h +++ b/inc/VUtilities.h @@ -48,5 +48,20 @@ namespace VUtilities return (*p == 0 ) ; } + + // friendlier colors + inline int color_id(unsigned int index) + { + // Array containing the specified ROOT color ID + int colors[] = {12, 633, 9, 418, 801, 881, 900, 32, 423, 393, 798, 616, 409}; + unsigned int size = sizeof(colors) / sizeof(colors[0]); + + if (index >= size) + { + return index; + } + return colors[index]; + } + } #endif diff --git a/src/VPlotAnasumHistograms.cpp b/src/VPlotAnasumHistograms.cpp index 24ad1c20..9c38fdf9 100644 --- a/src/VPlotAnasumHistograms.cpp +++ b/src/VPlotAnasumHistograms.cpp @@ -645,7 +645,7 @@ void VPlotAnasumHistograms::plot_CorrelatedSkyPlots() * * */ -TCanvas* VPlotAnasumHistograms::plot_theta2( double t2min, double t2max, int irbin, double setYMax, double setYMin ) +TCanvas* VPlotAnasumHistograms::plot_theta2( double t2min, double t2max, int irbin, double setYMax, double setYMin, bool containment_lines ) { // int iPlotPSF = 0; @@ -690,33 +690,43 @@ TCanvas* VPlotAnasumHistograms::plot_theta2( double t2min, double t2max, int irb htheta2_off->Draw( "hist e" ); htheta2_on->Draw( "hist e same" ); - // get 68% containment radius (up to theta2 = 0.05deg2) + // get 68% containment radius (up to theta2 = 0.2deg2) double nt2 = 0.; - for( int i = 1; i < htheta2_diff->GetXaxis()->FindBin( 0.05 ); i++ ) + double max_integration_radius = 0.1; + for( int i = 1; i < htheta2_diff->GetXaxis()->FindBin( max_integration_radius ); i++ ) { nt2 += htheta2_diff->GetBinContent( i ); } double nt2_68 = 0.; + double t_68 = -1.; + double t_95 = 0.; if( nt2 > 0. ) { - for( int i = 1; i < htheta2_diff->GetXaxis()->FindBin( 0.05 ); i++ ) + for( int i = 1; i < htheta2_diff->GetXaxis()->FindBin( max_integration_radius ); i++ ) { nt2_68 += htheta2_diff->GetBinContent( i ); - if( nt2_68 / nt2 > 0.68 ) + if( nt2_68 / nt2 > 0.68 && t_68 < 0. ) { - cout << "Theta2 containment radius (68%, binning dependent): " << htheta2_diff->GetXaxis()->GetBinLowEdge( i ) << " deg2" << endl; + t_68 = htheta2_diff->GetXaxis()->GetBinLowEdge( i ); + } + if( nt2_68 / nt2 > 0.95 ) + { + t_95 = htheta2_diff->GetXaxis()->GetBinLowEdge( i ); break; } } } + cout << "Theta2 containment radius (68%, binning dependent): " << t_68 << " deg2" << endl; + cout << "Theta2 containment radius (95%, binning dependent): " << t_95 << " deg2" << endl; c_t2_diff->cd(); + c_t2_diff->SetLeftMargin( 0.14 ); htheta2_diff->SetFillColor( 418 ); htheta2_diff->SetXTitle( "#Theta^{2} [deg^{2}]" ); htheta2_diff->SetTitle( "" ); sprintf( hname, "Number of events / %.3f deg^{2}", htheta2_diff->GetXaxis()->GetBinWidth( 2 ) ); htheta2_diff->SetYTitle( hname ); - htheta2_diff->GetYaxis()->SetTitleOffset( 1.6 ); + htheta2_diff->GetYaxis()->SetTitleOffset( 3.0 ); setHistogramPlottingStyle( htheta2_diff, 418, 2, 1, 1, irbin, 1001 ); htheta2_diff->SetAxisRange( t2min, t2max ); if( setYMax > -999. ) @@ -730,6 +740,18 @@ TCanvas* VPlotAnasumHistograms::plot_theta2( double t2min, double t2max, int irb htheta2_diff->Draw( "hist e" ); + if( containment_lines ) + { + TLine *iL_68 = new TLine( t_68, 0., t_68, htheta2_diff->GetMaximum()*0.5); + iL_68->SetLineStyle(2); + iL_68->SetLineColor(2); + iL_68->Draw(); + TLine *iL_95 = new TLine( t_95, 0., t_95, htheta2_diff->GetMaximum()*0.5); + iL_95->SetLineStyle(2); + iL_95->SetLineColor(2); + iL_95->Draw(); + } + if( fDebug ) { @@ -1132,7 +1154,7 @@ TCanvas* VPlotAnasumHistograms::plot_radec( int sPlot, double rmax, double zmin, { sprintf( hname, "c_skysig_%d_%d", fRunNumber, sPlot ); sprintf( htitle, "sky map, run %d", fRunNumber ); - c_skysig = new TCanvas( hname, htitle, 10, 10, 400, 400 ); + c_skysig = new TCanvas( hname, htitle, 10, 10, 800, 800 ); c_skysig->Draw(); c_skysig->SetRightMargin( 0.14 ); c_skysig->SetLeftMargin( 0.11 ); diff --git a/src/VPlotInstrumentResponseFunction.cpp b/src/VPlotInstrumentResponseFunction.cpp index d2d01a90..ad43c92b 100644 --- a/src/VPlotInstrumentResponseFunction.cpp +++ b/src/VPlotInstrumentResponseFunction.cpp @@ -102,7 +102,7 @@ bool VPlotInstrumentResponseFunction::addInstrumentResponseData( string iFile, d iTempIRFReader->setDebug( fDebug ); if( iColor < 0 ) { - iColor = fData.size() + 1; + iColor = VUtilities::color_id(fData.size() + 1); } if( iLineStyle < 0 ) { @@ -1574,7 +1574,9 @@ TCanvas* VPlotInstrumentResponseFunction::plotPSF( vector< double > i_Energy_TeV sprintf( hname, "Theta_ID_%d", iCumulative ); } TCanvas* c = new TCanvas( hname, hname, 10, 10, 600, 600 ); - c->Divide( TMath::Nint( sqrt( i_Energy_TeV_lin.size() ) ), TMath::Nint( sqrt( i_Energy_TeV_lin.size() ) ) ); + cout << "PSF plotting: " << i_Energy_TeV_lin.size() << " energies ("; + cout << TMath::Ceil( sqrt( i_Energy_TeV_lin.size() ) ) << ")" << endl; + c->Divide( TMath::Ceil( sqrt( i_Energy_TeV_lin.size() ) ), TMath::Ceil( sqrt( i_Energy_TeV_lin.size() ) ) ); for( unsigned int j = 0; j < i_Energy_TeV_lin.size(); j++ ) { c->cd( j + 1 ); @@ -1632,10 +1634,12 @@ TCanvas* VPlotInstrumentResponseFunction::plotPSF( vector< double > i_Energy_TeV TH1D* hCumu = get_Cumulative_Histogram( h, true, true ); if( iCumulative ) { + hCumu->SetLineColor( VUtilities::color_id(i+1) ); hCumu->Draw( "same" ); } else { + h->SetLineColor( VUtilities::color_id(i+1) ); // rebin if( iPlotTheta2 ) { diff --git a/src/printCrabSensitivity.cpp b/src/printCrabSensitivity.cpp index e786c9cd..78157bf8 100644 --- a/src/printCrabSensitivity.cpp +++ b/src/printCrabSensitivity.cpp @@ -44,7 +44,7 @@ void print_sensitivity( string anasum_file, double alpha = 1. / 6. ) t->GetEntry( t->GetEntries() - 1 ); - cout << "Rates (on/Off): " << Rate << "\t" << RateOff << endl; + cout << "Rates (on/off): " << Rate << "\t" << RateOff << endl; vector< double > fSourceStrength; fSourceStrength.push_back( 1. );