Skip to content

v492 dev5 #312

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
GernotMaier merged 8 commits into from
Oct 26, 2025
Merged

v492 dev5 #312

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
515a38f
Performance mscw_energy
GernotMaier Oct 24, 2025
422bbe2
perf improvements
GernotMaier Oct 24, 2025
16d9850
reserve
GernotMaier Oct 24, 2025
7de7f54
remove obsolete specators
GernotMaier Oct 24, 2025
b872940
fix memory leak
GernotMaier Oct 24, 2025
2133f11
mean energy optimization
GernotMaier Oct 24, 2025
66b9716
default tree number to 100
GernotMaier Oct 25, 2025
445f0df
astyle
GernotMaier Oct 25, 2025
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
83 changes: 35 additions & 48 deletions inc/VStatistics.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,71 +288,58 @@ namespace VStatistics

weighted by cos ze
*/
inline double interpolate( double w1, double ze1, double w2, double ze2, double ze, bool iCos = false,
double iLimitforInterpolation = 0.5, double iMinValidValue = -90. )

template <typename T>
inline T interpolate( T w1, T ze1, T w2, T ze2, T ze,
bool iCos = false,
T limit_low = static_cast<T>( 0.5 ),
T limit_up = static_cast<T>( 1e-5 ) )
{
// don't interpolate if both values are not valid
const T iMinValidValue = static_cast<T>(-98.0 );
const T iLimitforInterpolation = static_cast<T>( 0.05 );

// both values invalid
if( w1 < iMinValidValue && w2 < iMinValidValue )
{
return -99.;
}
return static_cast<T>(-99.0 );

// same x-value, don't interpolate
if( fabs( ze1 - ze2 ) < 1.e-3 )
// same x-value: average or return valid one
if( std::fabs( ze1 - ze2 ) < static_cast<T>( 1e-3 ) )
{
if( w1 < iMinValidValue )
{
return w2;
}
else if( w2 < iMinValidValue )
{
return w1;
}
return ( w1 + w2 ) / 2.;
if( w1 < iMinValidValue ) return w2;
if( w2 < iMinValidValue ) return w1;
return ( w1 + w2 ) / static_cast<T>( 2 );
}

// interpolate
double id = 0.;
double f1 = 0.;
double f2 = 0.;
T id = static_cast<T>( 0 );
T f1 = static_cast<T>( 0 );
T f2 = static_cast<T>( 0 );

if( iCos )
{
id = cos( ze1* TMath::DegToRad() ) - cos( ze2* TMath::DegToRad() );
f1 = 1. - ( cos( ze1* TMath::DegToRad() ) - cos( ze* TMath::DegToRad() ) ) / id;
f2 = 1. - ( cos( ze* TMath::DegToRad() ) - cos( ze2* TMath::DegToRad() ) ) / id;
const T d2r = static_cast<T>( TMath::DegToRad() );
const T cos_ze1 = std::cos( ze1* d2r );
const T cos_ze2 = std::cos( ze2* d2r );
const T cos_ze = std::cos( ze* d2r );

id = cos_ze1 - cos_ze2;
f1 = static_cast<T>( 1 ) - ( cos_ze1 - cos_ze ) / id;
f2 = static_cast<T>( 1 ) - ( cos_ze - cos_ze2 ) / id;
}
else
{
id = ze1 - ze2;
f1 = 1. - ( ze1 - ze ) / id;
f2 = 1. - ( ze - ze2 ) / id;
f1 = static_cast<T>( 1 ) - ( ze1 - ze ) / id;
f2 = static_cast<T>( 1 ) - ( ze - ze2 ) / id;
}

// one of the values is not valid:
// return valid value only when f1 or f2 > iLimitforInterPolation
// one value invalid → conditional return
if( w1 > iMinValidValue && w2 < iMinValidValue )
{
if( f1 > iLimitforInterpolation )
{
return w1;
}
else
{
return -99.;
}
}
else if( w1 < iMinValidValue && w2 > iMinValidValue )
{
if( f2 > iLimitforInterpolation )
{
return w2;
}
else
{
return -99.;
}
}
return ( f1 > iLimitforInterpolation ) ? w1 : static_cast<T>(-99.0 );
if( w1 < iMinValidValue && w2 > iMinValidValue )
return ( f2 > iLimitforInterpolation ) ? w2 : static_cast<T>(-99.0 );

// normal interpolation
return ( w1* f1 + w2* f2 );
}

Expand Down
16 changes: 16 additions & 0 deletions inc/VTMVADispAnalyzer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,22 @@ class VTMVADispAnalyzer
float fcross;
float fRcore;
float fEHeight;
// spectators (nowhere used)
float iMCe0;
float iMCxoff;
float iMCyoff;
float iMCxcore;
float iMCycore;
float iMCrcore;
float iNImages;
float cen_x;
float cen_y;
float cosphi;
float sinphi;
float temp1;
float temp2;
float temp3;
float temp4;

public:

Expand Down
2 changes: 1 addition & 1 deletion inc/VTableCalculator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ class VTableCalculator
bool createMedianApprox( int i, int j );
double getWeightMeanBinContent( TH2F*, int, int, double, double );
void fillMPV( TH2F*, int, int, TH1F*, double, double );
double interpolate( TH2F* h, double x, double y, bool iError );
double interpolate( TH2F* h, float x, float y, bool iError );
bool readHistograms();
void setBinning();
void setConstants( bool iPE = false );
Expand Down
33 changes: 16 additions & 17 deletions src/VDispAnalyzer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -785,32 +785,33 @@ void VDispAnalyzer::calculateMeanEnergy(
double* img_weight )
{
vector< double > energy_tel;
vector< double > energy_weight;
energy_tel.reserve( disp_energy_T.size() );
double energy_weight = 0.;
double w = 0.;
fdisp_energy = 0.;
fdisp_energy_NT = 0.;
fdisp_energy_chi = 0.;
fdisp_energy_dEs = 0.;

for( unsigned int i = 0; i < disp_energy_T.size(); i++ )
{
if( disp_energy_T[i] > 0. && img_size[i] > 0. )
{
energy_tel.push_back( disp_energy_T[i] );
if( img_weight )
{
energy_weight.push_back( img_size[i] * img_weight[i] * img_size[i] * img_weight[i] );
energy_weight = img_size[i] * img_weight[i] * img_size[i] * img_weight[i];
}
else
{
energy_weight.push_back( img_size[i] * img_size[i] );
energy_weight = img_size[i] * img_size[i];
}
fdisp_energy += ( double )disp_energy_T[i] * energy_weight;
w += energy_weight;
fdisp_energy_NT++;
}
}

double w = 0.;
fdisp_energy = 0.;
for( unsigned int j = 0; j < energy_tel.size(); j++ )
{
fdisp_energy += energy_tel[j] * energy_weight[j];
w += energy_weight[j];
}

fdisp_energy_NT = energy_tel.size();
if( w > 0. )
{
fdisp_energy /= w;
Expand All @@ -834,19 +835,17 @@ void VDispAnalyzer::calculateMeanEnergy(
// calculate chi2 and dE
// (note: different definition for dE
// than in lookup table code)
fdisp_energy_chi = 0.;
fdisp_energy_dEs = 0.;
for( unsigned int i = 0; i < energy_tel.size(); i++ )
{
fdisp_energy_chi += ( energy_tel[i] - fdisp_energy ) *
( energy_tel[i] - fdisp_energy ) /
fdisp_energy / fdisp_energy;
fdisp_energy_dEs += TMath::Abs( energy_tel[i] - fdisp_energy ) / fdisp_energy;
}
if( energy_tel.size() > 1 )
if( fdisp_energy_NT > 1 )
{
fdisp_energy_chi /= (( float )energy_tel.size() - 1. );
fdisp_energy_dEs /= ( float )energy_tel.size();
fdisp_energy_chi /= ( fdisp_energy_NT - 1 );
fdisp_energy_dEs /= ( float )fdisp_energy_NT;
}
else
{
Expand Down
30 changes: 15 additions & 15 deletions src/VTMVADispAnalyzer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,21 +37,21 @@ VTMVADispAnalyzer::VTMVADispAnalyzer( string iFile, vector<ULong64_t> iTelTypeLi
fEHeight = 0.;

// spectators (nowhere used)
float iMCe0 = 0.;
float iMCxoff = 0.;
float iMCyoff = 0.;
float iMCxcore = 0.;
float iMCycore = 0.;
float iMCrcore = 0.;
float iNImages = 0.;
float cen_x = 0;
float cen_y = 0;
float cosphi = 0.;
float sinphi = 0.;
float temp1 = 0.;
float temp2 = 0.;
float temp3 = 0.;
float temp4 = 0.;
iMCe0 = 0.;
iMCxoff = 0.;
iMCyoff = 0.;
iMCxcore = 0.;
iMCycore = 0.;
iMCrcore = 0.;
iNImages = 0.;
cen_x = 0;
cen_y = 0;
cosphi = 0.;
sinphi = 0.;
temp1 = 0.;
temp2 = 0.;
temp3 = 0.;
temp4 = 0.;

// list of telescope types: required to selected correct BDT weight file
fTelescopeTypeList = iTelTypeList;
Expand Down
42 changes: 22 additions & 20 deletions src/VTableCalculator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -600,7 +600,6 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
// fill width/length/energy into a 1D and 2D histogram
// (chi2 is here an external weight (from e.g. spectral weighting))
//============================================================================================================
// PRELIMINARY_START
if( fEnergy )
{
if( w[tel] < Oh[is][ir]->GetXaxis()->GetXmin() || w[tel] > Oh[is][ir]->GetXaxis()->GetXmax() )
Expand All @@ -609,7 +608,6 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
cout << Oh[is][ir]->GetXaxis()->GetXmin() << "\t" << Oh[is][ir]->GetXaxis()->GetXmax() << endl;
}
}
// PRELIMINARY_END
//============================================================================================================
Oh[is][ir]->Fill( w[tel], chi2 );
}
Expand Down Expand Up @@ -666,6 +664,10 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
vector< double > sigma2_tel;
vector< double > sigma2_tel_noRadiusWeigth;
vector< double > sigma_tel;
energy_tel.reserve( ntel );
sigma2_tel.reserve( ntel );
sigma2_tel_noRadiusWeigth.reserve( ntel );
sigma_tel.reserve( ntel );

// reset everything
for( tel = 0; tel < ntel; tel++ )
Expand Down Expand Up @@ -1016,7 +1018,7 @@ bool VTableCalculator::readHistograms()
return false;
}

double VTableCalculator::interpolate( TH2F* h, double x, double y, bool iError )
double VTableCalculator::interpolate( TH2F* h, float x, float y, bool iError )
{
if(!h )
{
Expand Down Expand Up @@ -1046,33 +1048,33 @@ double VTableCalculator::interpolate( TH2F* h, double x, double y, bool iError )
i_y--;
}

double e1 = 0.;
double e2 = 0.;
double v = 0.;
float e1 = 0.;
float e2 = 0.;
float v = 0.;

// first interpolate on distance axis, then on size axis
if(!iError )
{
e1 = VStatistics::interpolate( h->GetBinContent( i_x, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
h->GetBinContent( i_x, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
y, false, 0.5, 1.e-5 );
e2 = VStatistics::interpolate( h->GetBinContent( i_x + 1, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
h->GetBinContent( i_x + 1, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
y, false, 0.5, 1.e-5 );
v = VStatistics::interpolate( e1, h->GetXaxis()->GetBinCenter( i_x ),
e2, h->GetXaxis()->GetBinCenter( i_x + 1 ),
x, false, 0.5, 1.e-5 );
e1 = VStatistics::interpolate( static_cast<float>( h->GetBinContent( i_x, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
static_cast<float>( h->GetBinContent( i_x, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
y, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
e2 = VStatistics::interpolate( static_cast<float>( h->GetBinContent( i_x + 1, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
static_cast<float>( h->GetBinContent( i_x + 1, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
y, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
v = VStatistics::interpolate( e1, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x ) ),
e2, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x + 1 ) ),
x, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
}
else
{
e1 = VStatistics::interpolate( h->GetBinError( i_x, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
h->GetBinError( i_x, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
e1 = VStatistics::interpolate( static_cast<float>( h->GetBinError( i_x, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
static_cast<float>( h->GetBinError( i_x, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
y, false );
e2 = VStatistics::interpolate( h->GetBinError( i_x + 1, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
h->GetBinError( i_x + 1, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
e2 = VStatistics::interpolate( static_cast<float>( h->GetBinError( i_x + 1, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
static_cast<float>( h->GetBinError( i_x + 1, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
y, false );

v = VStatistics::interpolate( e1, h->GetXaxis()->GetBinCenter( i_x ), e2, h->GetXaxis()->GetBinCenter( i_x + 1 ), x, false );
v = VStatistics::interpolate( e1, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x ) ), e2, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x + 1 ) ), x, false );
}
// final check on consistency of results
// (don't expect to reconstruct anything below 1 GeV)
Expand Down
2 changes: 1 addition & 1 deletion src/trainTMVAforAngularReconstruction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -876,7 +876,7 @@ int main( int argc, char* argv[] )
iQualityCut = argv[7];
}
// TMVA options (default options derived from hyperparameter optimisation on CTAO prod3 simulations)
string iTMVAOptions = "NTrees=300:BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:GradBaggingFraction=0.5:nCuts=20:MaxDepth=10:";
string iTMVAOptions = "NTrees=100:BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:GradBaggingFraction=0.5:nCuts=20:MaxDepth=10:";
iTMVAOptions += "PruneMethod=ExpectedError:RegressionLossFunctionBDTG=Huber:MinNodeSize=0.02:VarTransform=N";
if( argc >= 9 )
{
Expand Down
Loading