Skip to content

Commit 84c3737

Browse files
GernotMaierGernotMaier
authored
Merge pull request #312 from VERITAS-Observatory/v492-dev5
v492 dev5
2 parents 81db735 + 445f0df commit 84c3737

File tree

7 files changed

+106
-102
lines changed

7 files changed

+106
-102
lines changed

‎inc/VStatistics.h‎

Lines changed: 35 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -288,71 +288,58 @@ namespace VStatistics
288288
289289
weighted by cos ze
290290
*/
291-
inline double interpolate( double w1, double ze1, double w2, double ze2, double ze, bool iCos = false,
292-
double iLimitforInterpolation = 0.5, double iMinValidValue = -90. )
291+
292+
template <typename T>
293+
inline T interpolate( T w1, T ze1, T w2, T ze2, T ze,
294+
bool iCos = false,
295+
T limit_low = static_cast<T>( 0.5 ),
296+
T limit_up = static_cast<T>( 1e-5 ) )
293297
{
294-
// don't interpolate if both values are not valid
298+
const T iMinValidValue = static_cast<T>(-98.0 );
299+
const T iLimitforInterpolation = static_cast<T>( 0.05 );
300+
301+
// both values invalid
295302
if( w1 < iMinValidValue && w2 < iMinValidValue )
296-
{
297-
return -99.;
298-
}
303+
return static_cast<T>(-99.0 );
299304

300-
// same x-value, don't interpolate
301-
if( fabs( ze1 - ze2 ) < 1.e-3 )
305+
// same x-value: average or return valid one
306+
if( std::fabs( ze1 - ze2 ) < static_cast<T>( 1e-3 ) )
302307
{
303-
if( w1 < iMinValidValue )
304-
{
305-
return w2;
306-
}
307-
else if( w2 < iMinValidValue )
308-
{
309-
return w1;
310-
}
311-
return ( w1 + w2 ) / 2.;
308+
if( w1 < iMinValidValue ) return w2;
309+
if( w2 < iMinValidValue ) return w1;
310+
return ( w1 + w2 ) / static_cast<T>( 2 );
312311
}
313312

314313
// interpolate
315-
double id = 0.;
316-
double f1 = 0.;
317-
double f2 = 0.;
314+
T id = static_cast<T>( 0 );
315+
T f1 = static_cast<T>( 0 );
316+
T f2 = static_cast<T>( 0 );
317+
318318
if( iCos )
319319
{
320-
id = cos( ze1* TMath::DegToRad() ) - cos( ze2* TMath::DegToRad() );
321-
f1 = 1. - ( cos( ze1* TMath::DegToRad() ) - cos( ze* TMath::DegToRad() ) ) / id;
322-
f2 = 1. - ( cos( ze* TMath::DegToRad() ) - cos( ze2* TMath::DegToRad() ) ) / id;
320+
const T d2r = static_cast<T>( TMath::DegToRad() );
321+
const T cos_ze1 = std::cos( ze1* d2r );
322+
const T cos_ze2 = std::cos( ze2* d2r );
323+
const T cos_ze = std::cos( ze* d2r );
324+
325+
id = cos_ze1 - cos_ze2;
326+
f1 = static_cast<T>( 1 ) - ( cos_ze1 - cos_ze ) / id;
327+
f2 = static_cast<T>( 1 ) - ( cos_ze - cos_ze2 ) / id;
323328
}
324329
else
325330
{
326331
id = ze1 - ze2;
327-
f1 = 1. - ( ze1 - ze ) / id;
328-
f2 = 1. - ( ze - ze2 ) / id;
332+
f1 = static_cast<T>( 1 ) - ( ze1 - ze ) / id;
333+
f2 = static_cast<T>( 1 ) - ( ze - ze2 ) / id;
329334
}
330335

331-
// one of the values is not valid:
332-
// return valid value only when f1 or f2 > iLimitforInterPolation
336+
// one value invalid → conditional return
333337
if( w1 > iMinValidValue && w2 < iMinValidValue )
334-
{
335-
if( f1 > iLimitforInterpolation )
336-
{
337-
return w1;
338-
}
339-
else
340-
{
341-
return -99.;
342-
}
343-
}
344-
else if( w1 < iMinValidValue && w2 > iMinValidValue )
345-
{
346-
if( f2 > iLimitforInterpolation )
347-
{
348-
return w2;
349-
}
350-
else
351-
{
352-
return -99.;
353-
}
354-
}
338+
return ( f1 > iLimitforInterpolation ) ? w1 : static_cast<T>(-99.0 );
339+
if( w1 < iMinValidValue && w2 > iMinValidValue )
340+
return ( f2 > iLimitforInterpolation ) ? w2 : static_cast<T>(-99.0 );
355341

342+
// normal interpolation
356343
return ( w1* f1 + w2* f2 );
357344
}
358345

‎inc/VTMVADispAnalyzer.h‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,22 @@ class VTMVADispAnalyzer
4848
float fcross;
4949
float fRcore;
5050
float fEHeight;
51+
// spectators (nowhere used)
52+
float iMCe0;
53+
float iMCxoff;
54+
float iMCyoff;
55+
float iMCxcore;
56+
float iMCycore;
57+
float iMCrcore;
58+
float iNImages;
59+
float cen_x;
60+
float cen_y;
61+
float cosphi;
62+
float sinphi;
63+
float temp1;
64+
float temp2;
65+
float temp3;
66+
float temp4;
5167

5268
public:
5369

‎inc/VTableCalculator.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -126,7 +126,7 @@ class VTableCalculator
126126
bool createMedianApprox( int i, int j );
127127
double getWeightMeanBinContent( TH2F*, int, int, double, double );
128128
void fillMPV( TH2F*, int, int, TH1F*, double, double );
129-
double interpolate( TH2F* h, double x, double y, bool iError );
129+
double interpolate( TH2F* h, float x, float y, bool iError );
130130
bool readHistograms();
131131
void setBinning();
132132
void setConstants( bool iPE = false );

‎src/VDispAnalyzer.cpp‎

Lines changed: 16 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -785,32 +785,33 @@ void VDispAnalyzer::calculateMeanEnergy(
785785
double* img_weight )
786786
{
787787
vector< double > energy_tel;
788-
vector< double > energy_weight;
788+
energy_tel.reserve( disp_energy_T.size() );
789+
double energy_weight = 0.;
790+
double w = 0.;
791+
fdisp_energy = 0.;
792+
fdisp_energy_NT = 0.;
793+
fdisp_energy_chi = 0.;
794+
fdisp_energy_dEs = 0.;
795+
789796
for( unsigned int i = 0; i < disp_energy_T.size(); i++ )
790797
{
791798
if( disp_energy_T[i] > 0. && img_size[i] > 0. )
792799
{
793800
energy_tel.push_back( disp_energy_T[i] );
794801
if( img_weight )
795802
{
796-
energy_weight.push_back( img_size[i] * img_weight[i] * img_size[i] * img_weight[i] );
803+
energy_weight = img_size[i] * img_weight[i] * img_size[i] * img_weight[i];
797804
}
798805
else
799806
{
800-
energy_weight.push_back( img_size[i] * img_size[i] );
807+
energy_weight = img_size[i] * img_size[i];
801808
}
809+
fdisp_energy += ( double )disp_energy_T[i] * energy_weight;
810+
w += energy_weight;
811+
fdisp_energy_NT++;
802812
}
803813
}
804814

805-
double w = 0.;
806-
fdisp_energy = 0.;
807-
for( unsigned int j = 0; j < energy_tel.size(); j++ )
808-
{
809-
fdisp_energy += energy_tel[j] * energy_weight[j];
810-
w += energy_weight[j];
811-
}
812-
813-
fdisp_energy_NT = energy_tel.size();
814815
if( w > 0. )
815816
{
816817
fdisp_energy /= w;
@@ -834,19 +835,17 @@ void VDispAnalyzer::calculateMeanEnergy(
834835
// calculate chi2 and dE
835836
// (note: different definition for dE
836837
// than in lookup table code)
837-
fdisp_energy_chi = 0.;
838-
fdisp_energy_dEs = 0.;
839838
for( unsigned int i = 0; i < energy_tel.size(); i++ )
840839
{
841840
fdisp_energy_chi += ( energy_tel[i] - fdisp_energy ) *
842841
( energy_tel[i] - fdisp_energy ) /
843842
fdisp_energy / fdisp_energy;
844843
fdisp_energy_dEs += TMath::Abs( energy_tel[i] - fdisp_energy ) / fdisp_energy;
845844
}
846-
if( energy_tel.size() > 1 )
845+
if( fdisp_energy_NT > 1 )
847846
{
848-
fdisp_energy_chi /= (( float )energy_tel.size() - 1. );
849-
fdisp_energy_dEs /= ( float )energy_tel.size();
847+
fdisp_energy_chi /= ( fdisp_energy_NT - 1 );
848+
fdisp_energy_dEs /= ( float )fdisp_energy_NT;
850849
}
851850
else
852851
{

‎src/VTMVADispAnalyzer.cpp‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -37,21 +37,21 @@ VTMVADispAnalyzer::VTMVADispAnalyzer( string iFile, vector<ULong64_t> iTelTypeLi
3737
fEHeight = 0.;
3838

3939
// spectators (nowhere used)
40-
float iMCe0 = 0.;
41-
float iMCxoff = 0.;
42-
float iMCyoff = 0.;
43-
float iMCxcore = 0.;
44-
float iMCycore = 0.;
45-
float iMCrcore = 0.;
46-
float iNImages = 0.;
47-
float cen_x = 0;
48-
float cen_y = 0;
49-
float cosphi = 0.;
50-
float sinphi = 0.;
51-
float temp1 = 0.;
52-
float temp2 = 0.;
53-
float temp3 = 0.;
54-
float temp4 = 0.;
40+
iMCe0 = 0.;
41+
iMCxoff = 0.;
42+
iMCyoff = 0.;
43+
iMCxcore = 0.;
44+
iMCycore = 0.;
45+
iMCrcore = 0.;
46+
iNImages = 0.;
47+
cen_x = 0;
48+
cen_y = 0;
49+
cosphi = 0.;
50+
sinphi = 0.;
51+
temp1 = 0.;
52+
temp2 = 0.;
53+
temp3 = 0.;
54+
temp4 = 0.;
5555

5656
// list of telescope types: required to selected correct BDT weight file
5757
fTelescopeTypeList = iTelTypeList;

‎src/VTableCalculator.cpp‎

Lines changed: 22 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -600,7 +600,6 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
600600
// fill width/length/energy into a 1D and 2D histogram
601601
// (chi2 is here an external weight (from e.g. spectral weighting))
602602
//============================================================================================================
603-
// PRELIMINARY_START
604603
if( fEnergy )
605604
{
606605
if( w[tel] < Oh[is][ir]->GetXaxis()->GetXmin() || w[tel] > Oh[is][ir]->GetXaxis()->GetXmax() )
@@ -609,7 +608,6 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
609608
cout << Oh[is][ir]->GetXaxis()->GetXmin() << "\t" << Oh[is][ir]->GetXaxis()->GetXmax() << endl;
610609
}
611610
}
612-
// PRELIMINARY_END
613611
//============================================================================================================
614612
Oh[is][ir]->Fill( w[tel], chi2 );
615613
}
@@ -666,6 +664,10 @@ double VTableCalculator::calc( int ntel, float* r, float* s, float* w, double* m
666664
vector< double > sigma2_tel;
667665
vector< double > sigma2_tel_noRadiusWeigth;
668666
vector< double > sigma_tel;
667+
energy_tel.reserve( ntel );
668+
sigma2_tel.reserve( ntel );
669+
sigma2_tel_noRadiusWeigth.reserve( ntel );
670+
sigma_tel.reserve( ntel );
669671

670672
// reset everything
671673
for( tel = 0; tel < ntel; tel++ )
@@ -1016,7 +1018,7 @@ bool VTableCalculator::readHistograms()
10161018
return false;
10171019
}
10181020

1019-
double VTableCalculator::interpolate( TH2F* h, double x, double y, bool iError )
1021+
double VTableCalculator::interpolate( TH2F* h, float x, float y, bool iError )
10201022
{
10211023
if(!h )
10221024
{
@@ -1046,33 +1048,33 @@ double VTableCalculator::interpolate( TH2F* h, double x, double y, bool iError )
10461048
i_y--;
10471049
}
10481050

1049-
double e1 = 0.;
1050-
double e2 = 0.;
1051-
double v = 0.;
1051+
float e1 = 0.;
1052+
float e2 = 0.;
1053+
float v = 0.;
10521054

10531055
// first interpolate on distance axis, then on size axis
10541056
if(!iError )
10551057
{
1056-
e1 = VStatistics::interpolate( h->GetBinContent( i_x, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
1057-
h->GetBinContent( i_x, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
1058-
y, false, 0.5, 1.e-5 );
1059-
e2 = VStatistics::interpolate( h->GetBinContent( i_x + 1, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
1060-
h->GetBinContent( i_x + 1, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
1061-
y, false, 0.5, 1.e-5 );
1062-
v = VStatistics::interpolate( e1, h->GetXaxis()->GetBinCenter( i_x ),
1063-
e2, h->GetXaxis()->GetBinCenter( i_x + 1 ),
1064-
x, false, 0.5, 1.e-5 );
1058+
e1 = VStatistics::interpolate( static_cast<float>( h->GetBinContent( i_x, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
1059+
static_cast<float>( h->GetBinContent( i_x, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
1060+
y, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
1061+
e2 = VStatistics::interpolate( static_cast<float>( h->GetBinContent( i_x + 1, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
1062+
static_cast<float>( h->GetBinContent( i_x + 1, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
1063+
y, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
1064+
v = VStatistics::interpolate( e1, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x ) ),
1065+
e2, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x + 1 ) ),
1066+
x, false, static_cast<float>( 0.5 ), static_cast<float>( 1.e-5 ) );
10651067
}
10661068
else
10671069
{
1068-
e1 = VStatistics::interpolate( h->GetBinError( i_x, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
1069-
h->GetBinError( i_x, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
1070+
e1 = VStatistics::interpolate( static_cast<float>( h->GetBinError( i_x, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
1071+
static_cast<float>( h->GetBinError( i_x, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
10701072
y, false );
1071-
e2 = VStatistics::interpolate( h->GetBinError( i_x + 1, i_y ), h->GetYaxis()->GetBinCenter( i_y ),
1072-
h->GetBinError( i_x + 1, i_y + 1 ), h->GetYaxis()->GetBinCenter( i_y + 1 ),
1073+
e2 = VStatistics::interpolate( static_cast<float>( h->GetBinError( i_x + 1, i_y ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y ) ),
1074+
static_cast<float>( h->GetBinError( i_x + 1, i_y + 1 ) ), static_cast<float>( h->GetYaxis()->GetBinCenter( i_y + 1 ) ),
10731075
y, false );
10741076

1075-
v = VStatistics::interpolate( e1, h->GetXaxis()->GetBinCenter( i_x ), e2, h->GetXaxis()->GetBinCenter( i_x + 1 ), x, false );
1077+
v = VStatistics::interpolate( e1, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x ) ), e2, static_cast<float>( h->GetXaxis()->GetBinCenter( i_x + 1 ) ), x, false );
10761078
}
10771079
// final check on consistency of results
10781080
// (don't expect to reconstruct anything below 1 GeV)

‎src/trainTMVAforAngularReconstruction.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -876,7 +876,7 @@ int main( int argc, char* argv[] )
876876
iQualityCut = argv[7];
877877
}
878878
// TMVA options (default options derived from hyperparameter optimisation on CTAO prod3 simulations)
879-
string iTMVAOptions = "NTrees=300:BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:GradBaggingFraction=0.5:nCuts=20:MaxDepth=10:";
879+
string iTMVAOptions = "NTrees=100:BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:GradBaggingFraction=0.5:nCuts=20:MaxDepth=10:";
880880
iTMVAOptions += "PruneMethod=ExpectedError:RegressionLossFunctionBDTG=Huber:MinNodeSize=0.02:VarTransform=N";
881881
if( argc >= 9 )
882882
{

0 commit comments

Comments
 (0)