Skip to content
Merged
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion DL3-IRFs/inc/VDL3IRFs.h
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ class VDL3IRFs
string name,
char* col_name,
char* col_unit,
bool MEV_BACKGROUND_UNIT = false );
bool MEV_BACKGROUND_UNIT = false,
bool include_uncertainty = false );


public:
Expand Down
109 changes: 79 additions & 30 deletions DL3-IRFs/src/VDL3IRFs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,7 @@ bool VDL3IRFs::write_fits_table_header( string irftype, char *instrument )
(char*)"RESPONSE",
(char*)"HDUCLAS1" );

if( irftype == "PSF_3GAUSS" )
if( irftype == "PSF_3GAUSS" || irftype == "PSF_TABLE" )
{
write_fits_keyword( (char*)"HDUCLAS2",
(char*)"PSF",
Expand Down Expand Up @@ -363,10 +363,10 @@ bool VDL3IRFs::write_psf_table( TH3F *h, char *instrument )
nRows = 0;
char* tType[nCol] = { (char*)"ENERG_LO",
(char*)"ENERG_HI",
(char*)"RAD_LO",
(char*)"RAD_HI",
(char*)"THETA_LO",
(char*)"THETA_HI",
(char*)"MIGRA_LO",
(char*)"MIGRA_HI",
(char*)"RPSF" };
char* tUnit[nCol] = { (char*)"TeV",
(char*)"TeV",
Expand All @@ -378,10 +378,10 @@ bool VDL3IRFs::write_psf_table( TH3F *h, char *instrument )
// e_true axis is x-axis
char x_form[10];
sprintf( x_form, "%dE", h->GetNbinsX() );
// Offset angle from source position
// Offset angle from source position (y-axis, RAD)
char y_form[10];
sprintf( y_form, "%dE", h->GetNbinsY() );
// offset angle axis is z-axis
// Field of view offset angle (z-axis, THETA)
char z_form[10];
sprintf( z_form, "%dE", h->GetNbinsZ() );
// mig
Expand All @@ -406,7 +406,7 @@ bool VDL3IRFs::write_psf_table( TH3F *h, char *instrument )
tType,
tForm,
tUnit,
"POINT SPREAD FUNCTION",
"PSF_TABLE",
&status ) )
Comment thread
orelgueta marked this conversation as resolved.
{
return printerror( status );
Expand Down Expand Up @@ -440,10 +440,12 @@ bool VDL3IRFs::write_psf_table( TH3F *h, char *instrument )
}
table.push_back( data );

write_table( table );
write_fits_table_header( "PSF_TABLE", instrument );

return true;
bool writing_success = write_table( table );
if( writing_success )
{
write_fits_table_header( "PSF_TABLE", instrument );
}
return writing_success;
}

/*
Expand Down Expand Up @@ -555,8 +557,10 @@ bool VDL3IRFs::write_psf_gauss( TH2F *h, char *instrument )
}

bool writing_success = write_table( table );
write_fits_table_header( "PSF_3GAUSS", instrument );

if( writing_success )
{
write_fits_table_header( "PSF_3GAUSS", instrument );
}
return writing_success;
}

Expand Down Expand Up @@ -764,7 +768,8 @@ bool VDL3IRFs::write_effarea( TH2F *h, char *instrument )
"EFFECTIVE AREA",
(char*)"EFFAREA",
(char*)"m**2",
false );
false,
true );
write_fits_table_header( "AEFF_2D", instrument );
Comment thread
orelgueta marked this conversation as resolved.
Outdated
return writing_success;
}
Expand All @@ -791,35 +796,50 @@ bool VDL3IRFs::write_histo2D( TH2F *h,
string name,
char* col_name,
char* col_unit,
bool MEV_BACKGROUND_UNIT )
bool MEV_BACKGROUND_UNIT,
bool include_uncertainty )
{
if( !h ) return false;

int status = 0;
const int nCol = 5;
const int nCol = include_uncertainty ? 6 : 5;
const int maxCol = 6; // Maximum possible columns (with uncertainty)
long nRows = h->GetNbinsX() * h->GetNbinsY();
nRows = 0;
char* tType[nCol] = { (char*)"ENERG_LO",
(char*)"ENERG_HI",
(char*)"THETA_LO",
(char*)"THETA_HI",
(char*)col_name };
char* tUnit[nCol] = { (char*)"TeV",
(char*)"TeV",
(char*)"deg",
(char*)"deg",
(char*)col_unit };

// Column names - create error column name using std::string
string err_col_name;
if( include_uncertainty )
{
err_col_name = string(col_name) + "_ERR";
}

// Arrays sized for maximum columns; fits_create_tbl uses only first nCol entries
char* tType[maxCol] = { (char*)"ENERG_LO",
(char*)"ENERG_HI",
(char*)"THETA_LO",
(char*)"THETA_HI",
(char*)col_name,
include_uncertainty ? (char*)err_col_name.c_str() : (char*)"" };

Comment on lines +824 to +838
char* tUnit[maxCol] = { (char*)"TeV",
(char*)"TeV",
(char*)"deg",
(char*)"deg",
(char*)col_unit,
(char*)col_unit };
char x_form[10];
sprintf( x_form, "%dE", h->GetNbinsX() );
char y_form[10];
sprintf( y_form, "%dE", h->GetNbinsY() );
char z_form[10];
sprintf( z_form, "%dE", h->GetNbinsX()*h->GetNbinsY() );
char* tForm[nCol] = { &x_form[0],
&x_form[0],
&y_form[0],
&y_form[0],
&z_form[0] };
char* tForm[maxCol] = { &x_form[0],
&x_form[0],
&y_form[0],
&y_form[0],
&z_form[0],
&z_form[0] };

///////////////
// create empty table
Expand All @@ -845,6 +865,18 @@ bool VDL3IRFs::write_histo2D( TH2F *h,
{
return printerror( status );
}
// set dimensions for uncertainty column if present
if( include_uncertainty )
{
if( fits_write_tdim( fptr,
6,
2,
naxes,
&status ) )
{
return printerror( status );
}
}
///////////////
// write data
vector< vector< float > > table = get_baseline_axes( h );
Expand All @@ -858,6 +890,7 @@ bool VDL3IRFs::write_histo2D( TH2F *h,

// data
vector< float > data;
vector< float > errors;
for( int j = 0; j < h->GetNbinsY(); j++ )
{
for( int i = 0; i < h->GetNbinsX(); i++ )
Expand All @@ -868,15 +901,31 @@ bool VDL3IRFs::write_histo2D( TH2F *h,
{
data.push_back( h->GetBinContent( i+1, j+1 )
* norm_mev_background[i] );
if( include_uncertainty )
{
// Get the symmetric bin error and scale it with the background normalisation
float error = h->GetBinError( i+1, j+1 );
errors.push_back( error * norm_mev_background[i] );
}
}
else
{
data.push_back( h->GetBinContent( i+1, j+1 ) );
if( include_uncertainty )
{
errors.push_back( h->GetBinError( i+1, j+1 ) );
}
}

}
}
table.push_back( data );

// Add uncertainty column if requested
if( include_uncertainty )
{
table.push_back( errors );
}

return write_table( table );
}
Expand Down
13 changes: 8 additions & 5 deletions DL3-IRFs/src/convertSensitivityFilesToFITS.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -87,12 +87,15 @@ int main( int argc, char* argv[] )
(TH2F*)fData->Get( "AngResEtrue_offaxis" ),
(char*)getArrayName(fData) );

/*
cout << "Writing gamma-ray point-spread function (3D table)" << endl;
a.write_psf_table(
(TH3F*)fData->Get( "AngularPSF2DEtrue_offaxis" ),
(char*)getArrayName(fData) );
*/
TH3F* psf_table_hist = (TH3F*)fData->Get( "AngularPSF2DEtrue_offaxis" );
if( psf_table_hist ) {
if( !a.write_psf_table( psf_table_hist, (char*)getArrayName(fData) ) ) {
cerr << "Error: Failed to write PSF_TABLE" << endl;
}
Comment thread
orelgueta marked this conversation as resolved.
} else {
Comment thread
orelgueta marked this conversation as resolved.
Outdated
cout << "Warning: AngularPSF2DEtrue_offaxis histogram not found, skipping PSF_TABLE" << endl;
}

// edisp
// note: note using migration matrix MigMatrixNoTheta2cut_offaxis
Expand Down
Loading