From 1793e12ff228687b41d2618538f53bf8c8000400 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 24 Dec 2025 11:46:14 +0100
Subject: [PATCH 01/16] obsolete code

---
 inc/VGammaHadronCuts.h   |   8 ---
 inc/VTMVAEvaluator.h     |   2 +-
 src/VGammaHadronCuts.cpp | 120 +++------------------------------------
 src/VTMVAEvaluator.cpp   |   2 +-
 4 files changed, 10 insertions(+), 122 deletions(-)

diff --git a/inc/VGammaHadronCuts.h b/inc/VGammaHadronCuts.h
index 4b900d2a..a19c1c9a 100644
--- a/inc/VGammaHadronCuts.h
+++ b/inc/VGammaHadronCuts.h
@@ -104,14 +104,6 @@ class VGammaHadronCuts : public VAnalysisUtilities
         unsigned int    fTMVAWeightFileIndex_Zmax;
         map< unsigned int, double > fTMVASignalEfficiency;
         map< unsigned int, double > fTMVA_MVACut;
-        string          fTMVAOptimizeSignalEfficiencyParticleNumberFile;
-        double          fTMVAParticleNumberFile_Conversion_Rate_to_seconds;
-        double          fTMVAOptimizeSignalEfficiencySignificance_Min;
-        double          fTMVAOptimizeSignalEfficiencySignalEvents_Min;
-        double          fTMVAOptimizeSignalEfficiencyObservationTime_h;
-        double          fTMVAFixedSignalEfficiencyMax;
-        double          fTMVAMinSourceStrength;
-        double          fTMVAFixedThetaCutMin;
         double          fTMVA_EvaluationResult;
         VTMVAEvaluatorResults* fTMVAEvaluatorResults;
 
diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h
index cda5f993..34e5af1a 100644
--- a/inc/VTMVAEvaluator.h
+++ b/inc/VTMVAEvaluator.h
@@ -237,7 +237,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         }
         void   setSignalEfficiency( double iSignalEfficiency = -99. );
         void   setSignalEfficiency( map< unsigned int, double > iMSignalEfficiency );
-        void   setSmoothAndInterPolateMVAValues( bool iS = true )
+        void   setSmoothAndInterpolateMVAValues( bool iS = true )
         {
             fSmoothAndInterpolateMVAValues = iS;
         }
diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp
index f1f87a2b..37632803 100644
--- a/src/VGammaHadronCuts.cpp
+++ b/src/VGammaHadronCuts.cpp
@@ -1,7 +1,6 @@
 /*! \class VGammaHadronCuts
   \brief class containing parameter cut definitions
 
-
   List of cut selectors:
 
   gamma/hadron cut selector values consist of two digits: ID1+ID2*10
@@ -28,19 +27,6 @@
     cut selector = 22 : apply event probability cuts
     cut selector = 10 : apply cuts from a tree AND apply MSCW/MSCL cuts
 
-  direction cut selector value consist of one digit:
-
-     0: energy independent direction cut (theta2 cut; set fArrayTheta2_max and fArrayTheta2_min)
-        (in cut file: *theta2cut <theta2cut_min> <theta2cut_max>
-     1: energy dependent TF1 function
-        (in cut file: *theta2file ...)
-     2: graph from instrument response function calculation
-        (in cut file: *theta2file ... IRF)
-     3: TMVA: use gamma/hadron part evaluation (no direction cuts applied in VGammaHadronCuts)
-     4: TMVA: get direction cut from TMVA evaluator for each event
-     5: TMVA: get direction cut from theta2 graph obtained during initialization of TMVA evaluator
-
-
 */
 
 #include "VGammaHadronCuts.h"
@@ -83,13 +69,6 @@ VGammaHadronCuts::VGammaHadronCuts()
     fTMVAWeightFile = "";
     fTMVASignalEfficiency.clear();
     fTMVA_MVACut.clear();
-    fTMVAOptimizeSignalEfficiencyParticleNumberFile = "";
-    fTMVAParticleNumberFile_Conversion_Rate_to_seconds = 60.;
-    fTMVAOptimizeSignalEfficiencySignificance_Min = 5.;
-    fTMVAOptimizeSignalEfficiencySignalEvents_Min = 10.;
-    fTMVAOptimizeSignalEfficiencyObservationTime_h = 50.;
-    fTMVAFixedSignalEfficiencyMax = 1.;
-    fTMVAFixedThetaCutMin = 0.;
     fTMVA_EvaluationResult = -99.;
     fTMVAEvaluatorResults = 0;
 
@@ -568,60 +547,6 @@ bool VGammaHadronCuts::readCuts( string i_cutfilename, int iPrint )
                     }
                 }
             }
-            else if( iCutVariable == "TMVACUTS" )
-            {
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAOptimizeSignalEfficiencySignificance_Min;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAOptimizeSignalEfficiencySignalEvents_Min;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> temp;
-                    // observing time is given as "50h", or "5m", "5s"
-                    if( temp.find( "h" ) != temp.npos )
-                    {
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h = atof( temp.substr( 0, temp.find( "h" ) ).c_str() );
-                    }
-                    else if( temp.find( "m" ) != temp.npos )
-                    {
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h = atof( temp.substr( 0, temp.find( "m" ) ).c_str() );
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h /= 60.;
-                    }
-                    else if( temp.find( "s" ) != temp.npos )
-                    {
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h = atof( temp.substr( 0, temp.find( "s" ) ).c_str() );
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h /= 3600.;
-                    }
-                    else
-                    {
-                        fTMVAOptimizeSignalEfficiencyObservationTime_h = atof( temp.c_str() );
-                    }
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAOptimizeSignalEfficiencyParticleNumberFile;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAFixedSignalEfficiencyMax;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAMinSourceStrength;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAFixedThetaCutMin;
-                }
-                if(!( is_stream >> std::ws ).eof() )
-                {
-                    is_stream >> fTMVAParticleNumberFile_Conversion_Rate_to_seconds;
-                }
-            }
             else if( iCutVariable == "TMVASignalEfficiency" )
             {
                 while(!( is_stream >> std::ws ).eof() )
@@ -775,24 +700,10 @@ void VGammaHadronCuts::printCutSummary()
         cout << "weight files: " << fTMVAWeightFile;
         cout << " (" << fTMVAWeightFileIndex_Emin << "," << fTMVAWeightFileIndex_Emax << ")";
         cout << " (" << fTMVAWeightFileIndex_Zmin << "," << fTMVAWeightFileIndex_Zmax << ")" << endl;
-        if( fTMVAOptimizeSignalEfficiencyParticleNumberFile.size() > 0. )
-        {
-            cout << "using optimal signal efficiency with a requirement of ";
-            cout << fTMVAOptimizeSignalEfficiencySignificance_Min << " sigma and at least ";
-            cout << fTMVAOptimizeSignalEfficiencySignalEvents_Min << " signal events in ";
-            cout << fTMVAOptimizeSignalEfficiencyObservationTime_h << " h observing time" << endl;
-            cout << "reading particle counts from " << fTMVAOptimizeSignalEfficiencyParticleNumberFile;
-            cout << " (unit of rate graphs: " << fTMVAParticleNumberFile_Conversion_Rate_to_seconds << ")" << endl;
-            cout << "   (max signal efficiency: " << fTMVAFixedSignalEfficiencyMax << ")" << endl;
-            cout << "   (min source strength: " << fTMVAMinSourceStrength << ")" << endl;
-        }
-        else
+        if( fDebug )
         {
-            if( fDebug )
-            {
-                printSignalEfficiency();
-                printTMVA_MVACut();
-            }
+            printSignalEfficiency();
+            printTMVA_MVACut();
         }
     }
     // other cut parameters
@@ -1010,7 +921,7 @@ bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructio
     }
 
     /////////////////////////////////////////////////////////////////////////////////
-    // apply cuts on second max
+    // apply cuts on size second max
     if( fData->SizeSecondMax < fCut_SizeSecondMax_min || fData->SizeSecondMax > fCut_SizeSecondMax_max )
     {
         if( bCount && fStats )
@@ -1396,26 +1307,13 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile,
 
     fTMVAEvaluator->setDebug( fDebug );
     // smoothing of MVA values
-    fTMVAEvaluator->setSmoothAndInterPolateMVAValues( true );
-    // set parameters for optimal MVA cut value search
-    // (always assume an alpha value of 0.2)
-    if( fTMVAOptimizeSignalEfficiencyParticleNumberFile.size() > 0. )
-    {
-        fTMVAEvaluator->setSensitivityOptimizationParameters( fTMVAOptimizeSignalEfficiencySignificance_Min,
-                fTMVAOptimizeSignalEfficiencySignalEvents_Min,
-                fTMVAOptimizeSignalEfficiencyObservationTime_h,
-                1. / 5. );
-        fTMVAEvaluator->setSensitivityOptimizationFixedSignalEfficiency( fTMVAFixedSignalEfficiencyMax );
-        fTMVAEvaluator->setParticleNumberFile(
-            fTMVAOptimizeSignalEfficiencyParticleNumberFile, fTMVAParticleNumberFile_Conversion_Rate_to_seconds );
-        fTMVAEvaluator->setSensitivityOptimizationMinSourceStrength( fTMVAMinSourceStrength );
-    }
-    // set a constant signal efficiency
-    else if( fTMVASignalEfficiency.size() > 0 )
+    fTMVAEvaluator->setSmoothAndInterpolateMVAValues( true );
+    // constant signal efficiency
+    if( fTMVASignalEfficiency.size() > 0 )
     {
         fTMVAEvaluator->setSignalEfficiency( fTMVASignalEfficiency );
     }
-    // set a fixed probability threshold or (for TMVA) a fixed MVA cut value
+    // fixed threshold or (for TMVA) a fixed MVA cut value
     else if( fTMVA_MVACut.size() > 0 )
     {
         fTMVAEvaluator->setTMVACutValue( fTMVA_MVACut );
@@ -1819,7 +1717,6 @@ double VGammaHadronCuts::getReconstructedXcore()
     {
         return -9999.;
     }
-
     return fData->Xcore;
 }
 
@@ -1829,6 +1726,5 @@ double VGammaHadronCuts::getReconstructedYcore()
     {
         return -9999.;
     }
-
     return fData->Ycore;
 }
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 32a9bc1e..a4028922 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -52,7 +52,7 @@ void VTMVAEvaluator::reset()
     fTMVA_EvaluationResult = -99.;
     fTMVACutValueNoVec = -99.;
 
-    setSmoothAndInterPolateMVAValues();
+    setSmoothAndInterpolateMVAValues();
 }
 
 /*

From 1fa12c65547cc511972a33e68a29d09d0b4d943d Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 24 Dec 2025 13:47:19 +0100
Subject: [PATCH 02/16] more cleanup

---
 inc/VInstrumentResponseFunctionRunParameter.h |  2 +-
 inc/VTMVAEvaluator.h                          | 14 ++++++-------
 src/VGammaHadronCuts.cpp                      |  2 +-
 src/VTMVAEvaluator.cpp                        | 21 +++++++++----------
 4 files changed, 19 insertions(+), 20 deletions(-)

diff --git a/inc/VInstrumentResponseFunctionRunParameter.h b/inc/VInstrumentResponseFunctionRunParameter.h
index 486ba1a7..c0f49312 100644
--- a/inc/VInstrumentResponseFunctionRunParameter.h
+++ b/inc/VInstrumentResponseFunctionRunParameter.h
@@ -114,7 +114,7 @@ class VInstrumentResponseFunctionRunParameter : public TNamed
         bool                  readRunParameterFromTextFile( string iFile );
         bool                  testRunparameters();
 
-        ClassDef( VInstrumentResponseFunctionRunParameter, 19 );
+        ClassDef( VInstrumentResponseFunctionRunParameter, 20 );
 };
 
 #endif
diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h
index 34e5af1a..589c2cac 100644
--- a/inc/VTMVAEvaluator.h
+++ b/inc/VTMVAEvaluator.h
@@ -117,8 +117,8 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         double                  fOptimizationBackgroundAlpha;
         double                  fOptimizationObservingTime_h;
 
-        double   fTMVA_EvaluationResult;         // result from TVMA evaluator
-        bool     fSmoothAndInterpolateMVAValues;
+        double   fTMVA_EvaluationResult;         // result from TMVA evaluator
+        bool     fsmoothAndInterpolateMVAValues;
         double   fAverageZenithPerRun;           // (rough) average zenith angle of run
 
         string   fTMVAMethodName;
@@ -169,10 +169,10 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         double           interpolate_mva_evaluation();
         TGraph*          readNonNoffGraphsFromFile( TFile* iF, double i_ze_min, double i_ze_max, bool bIsOn = true );
         void             reset();
-        void             smoothAndInterPolateMVAValue( TH1F*, TH1F*,
+        void             smoothAndInterpolateMVAValue( TH1F*, TH1F*,
                 unsigned int iE_min, unsigned int iE_max, unsigned int iZ_min, unsigned int iZ_max );
-        void             smoothAndInterPolateMVAValue_EnergyOnly( TH1F*, TH1F* );
-        void             smoothAndInterPolateMVAValue_Energy_and_Zenith( TH1F*, TH1F*,
+        void             smoothAndInterpolateMVAValue_EnergyOnly( TH1F*, TH1F* );
+        void             smoothAndInterpolateMVAValue_Energy_and_Zenith( TH1F*, TH1F*,
                 unsigned int iE_min, unsigned int iE_max, unsigned int iZ_min, unsigned int iZ_max );
 
     public:
@@ -237,9 +237,9 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         }
         void   setSignalEfficiency( double iSignalEfficiency = -99. );
         void   setSignalEfficiency( map< unsigned int, double > iMSignalEfficiency );
-        void   setSmoothAndInterpolateMVAValues( bool iS = true )
+        void   setsmoothAndInterpolateMVAValues( bool iS = true )
         {
-            fSmoothAndInterpolateMVAValues = iS;
+            fsmoothAndInterpolateMVAValues = iS;
         }
         void   setSpectralIndexForEnergyWeighting( double iS = -2. )
         {
diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp
index 37632803..9c1f325a 100644
--- a/src/VGammaHadronCuts.cpp
+++ b/src/VGammaHadronCuts.cpp
@@ -1307,7 +1307,7 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile,
 
     fTMVAEvaluator->setDebug( fDebug );
     // smoothing of MVA values
-    fTMVAEvaluator->setSmoothAndInterpolateMVAValues( true );
+    fTMVAEvaluator->setsmoothAndInterpolateMVAValues( true );
     // constant signal efficiency
     if( fTMVASignalEfficiency.size() > 0 )
     {
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index a4028922..3de92838 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -52,7 +52,7 @@ void VTMVAEvaluator::reset()
     fTMVA_EvaluationResult = -99.;
     fTMVACutValueNoVec = -99.;
 
-    setSmoothAndInterpolateMVAValues();
+    setsmoothAndInterpolateMVAValues();
 }
 
 /*
@@ -84,7 +84,6 @@ vector< string > VTMVAEvaluator::getTrainingVariables( string iXMLFile, vector<
 
     while( getline( is, is_line ) )
     {
-
         // number of variables
         if( is_line.find( "NVar=\"" ) != string::npos )
         {
@@ -368,7 +367,7 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
     //////////////////////////////////////////////////////////////////////////////////////
     // create and initialize TMVA readers
     // loop over all  energy bins: open one weight (XML) file per energy bin
-    //looping over spectral energy and zenith angle bins
+    // looping over spectral energy and zenith angle bins
     for( unsigned int b = 0; b < fTMVAData.size(); b++ )
     {
         fTMVAData[b]->fTMVAReader = new TMVA::Reader();
@@ -510,9 +509,9 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
     }
 
     // smooth and interpolate
-    if( fParticleNumberFileName.size() > 0 && fSmoothAndInterpolateMVAValues )
+    if( fParticleNumberFileName.size() > 0 && fsmoothAndInterpolateMVAValues )
     {
-        smoothAndInterPolateMVAValue( 0, 0,
+        smoothAndInterpolateMVAValue( 0, 0,
                                       iWeightFileIndex_Emin, iWeightFileIndex_Emax,
                                       iWeightFileIndex_Zmin, iWeightFileIndex_Zmax );
     }
@@ -1295,7 +1294,7 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
 
     //////////////////////////////////////////////////////
     // loop over different source strengths (in Crab Units)
-    // (hardwired: start at 0.001 CU to 30 CU)
+    // (up to 30 CU)
     unsigned int iSourceStrengthStepSizeN =
         ( unsigned int )(( log10( 30. ) - log10( fOptimizationMinSourceStrength ) ) / 0.005 );
     cout << "VTVMAEvaluator::optimizeSensitivity(), source strength steps: " << iSourceStrengthStepSizeN << endl;
@@ -1582,7 +1581,7 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
  (energy axis only)
 
 */
-void VTMVAEvaluator::smoothAndInterPolateMVAValue_EnergyOnly(
+void VTMVAEvaluator::smoothAndInterpolateMVAValue_EnergyOnly(
     TH1F* effS, TH1F* effB )
 {
     int z = 0;
@@ -1659,7 +1658,7 @@ void VTMVAEvaluator::smoothAndInterPolateMVAValue_EnergyOnly(
  (energy and zenith angle dependent)
 
 */
-void VTMVAEvaluator::smoothAndInterPolateMVAValue_Energy_and_Zenith(
+void VTMVAEvaluator::smoothAndInterpolateMVAValue_Energy_and_Zenith(
     TH1F* effS, TH1F* effB,
     unsigned int iWeightFileIndex_Emin, unsigned int iWeightFileIndex_Emax,
     unsigned int iWeightFileIndex_Zmin, unsigned int iWeightFileIndex_Zmax )
@@ -1770,7 +1769,7 @@ void VTMVAEvaluator::smoothAndInterPolateMVAValue_Energy_and_Zenith(
  note: signal and background efficiencies are not updated
 
 */
-void VTMVAEvaluator::smoothAndInterPolateMVAValue(
+void VTMVAEvaluator::smoothAndInterpolateMVAValue(
     TH1F* effS, TH1F* effB,
     unsigned int iWeightFileIndex_Emin, unsigned int iWeightFileIndex_Emax,
     unsigned int iWeightFileIndex_Zmin, unsigned int iWeightFileIndex_Zmax )
@@ -1786,12 +1785,12 @@ void VTMVAEvaluator::smoothAndInterPolateMVAValue(
     // energy dependent TMVA cut optimization only
     if( iWeightFileIndex_Zmax == iWeightFileIndex_Zmin )
     {
-        smoothAndInterPolateMVAValue_EnergyOnly( effS, effB );
+        smoothAndInterpolateMVAValue_EnergyOnly( effS, effB );
     }
     // energy and zenith angle dependent
     else
     {
-        smoothAndInterPolateMVAValue_Energy_and_Zenith(
+        smoothAndInterpolateMVAValue_Energy_and_Zenith(
             effS, effB,
             iWeightFileIndex_Emin, iWeightFileIndex_Emax,
             iWeightFileIndex_Zmin, iWeightFileIndex_Zmax );

From cd2e3aa66b8d8f65a5abbcd83f433e76719c0ae7 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 24 Dec 2025 13:55:47 +0100
Subject: [PATCH 03/16] remove obsolete code

---
 Makefile                              |  55 ---------
 src/writeParticleRateFilesForTMVA.cpp | 164 --------------------------
 2 files changed, 219 deletions(-)
 delete mode 100644 src/writeParticleRateFilesForTMVA.cpp

diff --git a/Makefile b/Makefile
index a6ec7d4e..a3725f8c 100644
--- a/Makefile
+++ b/Makefile
@@ -226,7 +226,6 @@ all VTS:	evndisp \
 	VTS.getRunListFromDB \
 	VTS.getLaserRunFromDB \
 	VTS.getRun_TimeElevAzim \
-	writeParticleRateFilesForTMVA \
 	writelaserinDB \
 	logFile \
 	printCrabSensitivity \
@@ -1060,60 +1059,6 @@ writeVTSWPPhysSensitivityFiles:	$(WRITEVTSPHYSOBJ)
 	$(LD) $(LDFLAGS) $^ $(GLIBS) $(OutPutOpt) ./bin/$@
 	@echo "$@ done"
 
-########################################################
-# writeParticleRateFilesForTMVA
-########################################################
-WRITEPARTPHYSOBJ=	./obj/writeParticleRateFilesForTMVA.o \
-			./obj/VGlobalRunParameter.o ./obj/VGlobalRunParameter_Dict.o \
-			./obj/CRunSummary.o ./obj/CRunSummary_Dict.o \
-			./obj/VAstronometry.o ./obj/VAstronometry_Dict.o \
-			./obj/VInstrumentResponseFunctionReader.o ./obj/VInstrumentResponseFunctionReader_Dict.o \
-			./obj/VSensitivityCalculator.o ./obj/VSensitivityCalculator_Dict.o \
-			./obj/CEffArea.o ./obj/CEffArea_Dict.o \
-			./obj/VAnalysisUtilities.o ./obj/VAnalysisUtilities_Dict.o \
-			./obj/VHistogramUtilities.o ./obj/VHistogramUtilities_Dict.o \
-			./obj/VInstrumentResponseFunctionData.o ./obj/VInstrumentResponseFunctionData_Dict.o \
-			./obj/VPlotUtilities.o ./obj/VPlotUtilities_Dict.o \
-			./obj/VGammaHadronCuts.o ./obj/VGammaHadronCuts_Dict.o \
-			./obj/VGammaHadronCutsStatistics.o ./obj/VGammaHadronCutsStatistics_Dict.o \
-			./obj/VTMVAEvaluator.o ./obj/VTMVAEvaluator_Dict.o \
-			./obj/VTMVARunDataEnergyCut.o ./obj/VTMVARunDataEnergyCut_Dict.o \
-			./obj/VTMVARunDataZenithCut.o ./obj/VTMVARunDataZenithCut_Dict.o \
-			./obj/VInstrumentResponseFunctionRunParameter.o ./obj/VInstrumentResponseFunctionRunParameter_Dict.o \
-			./obj/Ctelconfig.o ./obj/CData.o  \
-			./obj/VTMVADispAnalyzer.o ./obj/VMeanScaledVariables.o \
-			./obj/VDispAnalyzer.o ./obj/VDispTableReader.o ./obj/VDispTableReader_Dict.o ./obj/VDispTableAnalyzer.o \
-			./obj/VEmissionHeightCalculator.o ./obj/VSimpleStereoReconstructor.o ./obj/VGrIsuAnalyzer.o \
-			./obj/VSpectralFitter.o ./obj/VSpectralFitter_Dict.o \
-			./obj/VEnergyThreshold.o ./obj/VEnergyThreshold_Dict.o \
-			./obj/VRunList.o ./obj/VRunList_Dict.o \
-			./obj/VEnergySpectrumfromLiterature.o ./obj/VEnergySpectrumfromLiterature_Dict.o \
-			./obj/VEnergySpectrum.o ./obj/VEnergySpectrum_Dict.o \
-			./obj/VMathsandFunctions.o ./obj/VMathsandFunctions_Dict.o  \
-			./obj/VDifferentialFlux.o ./obj/VDifferentialFlux_Dict.o \
-			./obj/VMonteCarloRateCalculator.o ./obj/VMonteCarloRateCalculator_Dict.o \
-			./obj/VMonteCarloRunHeader.o ./obj/VMonteCarloRunHeader_Dict.o \
-			./obj/VStatistics_Dict.o \
-			./obj/VEvndispRunParameter.o ./obj/VEvndispRunParameter_Dict.o \
-			./obj/VSkyCoordinatesUtilities.o \
-			./obj/VTimeMask.o ./obj/VTimeMask_Dict.o \
-			./obj/VAnaSumRunParameter.o ./obj/VAnaSumRunParameter_Dict.o \
-			./obj/VImageCleaningRunParameter.o ./obj/VImageCleaningRunParameter_Dict.o \
-			./obj/VDispAnalyzer.o ./obj/VDispTableReader.o ./obj/VDispTableReader_Dict.o ./obj/VDispTableAnalyzer.o \
-			./obj/VTMVADispAnalyzer.o \
-			./obj/VUtilities.o
-
-ifeq ($(ASTRONMETRY),-DASTROSLALIB)
-    WRITECTAPHYSOBJ += ./obj/VASlalib.o
-endif
-
-./obj/writeParticleRateFilesForTMVA.o: 	./src/writeParticleRateFilesForTMVA.cpp
-	$(CXX) $(CXXFLAGS) -c -o $@ $<
-
-writeParticleRateFilesForTMVA:	$(WRITEPARTPHYSOBJ)
-	$(LD) $(LDFLAGS) $^ $(GLIBS) $(OutPutOpt) ./bin/$@
-	@echo "$@ done"
-
 ########################################################
 # combineLookupTables
 ########################################################
diff --git a/src/writeParticleRateFilesForTMVA.cpp b/src/writeParticleRateFilesForTMVA.cpp
deleted file mode 100644
index 020c7640..00000000
--- a/src/writeParticleRateFilesForTMVA.cpp
+++ /dev/null
@@ -1,164 +0,0 @@
-/* writeParticleRateFilesForTMVA
-
-   write files with particle number spectra for on (gamma) and off (protons+electrons) counts
-   (differential counting rates)
-
-   files are needed e.g. for setting the optimal cut value for TMVA cuts
-
-   Input:
-       * Combined anasum files for each zenith angle bin
-
-   Output:
-       * root file with signal and background rates per zenith angle and energy bi
-
-*/
-
-
-#include <iostream>
-#include <string>
-#include <vector>
-
-#include "TFile.h"
-#include "TTree.h"
-#include "TH1D.h"
-#include "TH2D.h"
-#include "TMath.h"
-#include "TGraph2DErrors.h"
-
-using namespace std;
-
-/*
-
-    fill graphs with particle numbers for on (gamma) and off (protons+electrons) counts
-
-    files are needed e.g. for setting the optimal cut value for TMVA cuts
-
-*/
-int main( int argc, char* argv[] )
-{
-    if( argc < 2 )
-    {
-        cout << endl;
-        cout << "writeParticleRateFilesForTMVA <combined anasum root file> <outputfile>";
-        cout << endl;
-        exit( EXIT_FAILURE );
-    }
-
-    cout << endl;
-    cout << "writeParticleRateFilesForTMVA" << endl;
-    cout << "========================================" << endl;
-    cout << endl ;
-
-    string iDataDir = argv[1];
-    string iOutFil = argv[2];
-
-    /* Name and directory to the COMBINED anasum root file */
-    TFile* f1 = new TFile( iDataDir.c_str(), "OPEN" );
-    if( f1->IsZombie() )
-    {
-        cout << "Error reading file " << iDataDir << endl;
-        cout << "exiting..." << endl;
-        exit( EXIT_FAILURE );
-    }
-
-    char reddataname[256], histnameOn[256], histnameOff[256];
-    TTree* RunSum = ( TTree* )f1->Get( "total_1/stereo/tRunSummary" );
-    if(!RunSum )
-    {
-        cout << "Error reading run summary tree" << endl;
-        cout << "exiting..." << endl;
-        exit( EXIT_FAILURE );
-    }
-    if( RunSum->GetEntries() == 0 )
-    {
-        cout << "Error reading run summary tree - empty tree" << endl;
-        cout << "exiting..." << endl;
-        exit( EXIT_FAILURE );
-    }
-
-    Int_t runOn = 0;
-    Double_t OffNorm = 0.;
-    Double_t elevationOff = 0.;
-    Double_t  alphaNorm = 0.;
-    RunSum->SetBranchAddress( "runOn", &runOn );
-    RunSum->SetBranchAddress( "OffNorm", &OffNorm );
-    RunSum->SetBranchAddress( "elevationOff", &elevationOff );
-
-    // assume that all histograms for all runs have the same binning
-    RunSum->GetEntry( 0 );
-    sprintf( histnameOn, "run_%d/stereo/energyHistograms/herecCounts_on", runOn );
-    sprintf( histnameOff, "run_%d/stereo/energyHistograms/herecCounts_off", runOn );
-    TH1D* t2 = ( TH1D* )f1->Get( histnameOff );
-    unsigned int iBin = t2->GetNbinsX();
-
-    TGraph2DErrors* tRatePerEnergyON  = new TGraph2DErrors();
-    tRatePerEnergyON->SetMarkerStyle( 20 );
-    TGraph2DErrors* tRatePerEnergyOFF = new TGraph2DErrors();
-    tRatePerEnergyOFF->SetMarkerStyle( 20 );
-    TGraph2DErrors* tRatePerEnergySignal = new TGraph2DErrors();
-    tRatePerEnergySignal->SetMarkerStyle( 20 );
-
-    int m = 0;
-    for( unsigned int i = 0; i < ( RunSum->GetEntries() - 1 ); i++ )
-    {
-
-        RunSum->GetEntry( i );
-        alphaNorm = OffNorm;
-        double iZenithperrun = 90.0 - elevationOff;
-        sprintf( histnameOn, "run_%d/stereo/energyHistograms/herecCounts_on", runOn );
-        sprintf( histnameOff, "run_%d/stereo/energyHistograms/herecCounts_off", runOn );
-        TH1D* t1 = ( TH1D* )f1->Get( histnameOn );
-        TH1D* t2 = ( TH1D* )f1->Get( histnameOff );
-
-        // get the observation time of each run
-        sprintf( reddataname, "run_%d/stereo/pointingDataReduced", runOn );
-        TTree* tReddata = ( TTree* )f1->Get( reddataname );
-        if(!tReddata )
-        {
-            cout << "Error reading pointingDataReduced tree for run " << runOn << endl;
-            cout << "exiting..." << endl;
-            exit( EXIT_FAILURE );
-        }
-        Double_t tObs = tReddata->GetEntries();
-
-        cout << "Run Number: " << runOn ;
-        cout << ", Zenith angle = " << iZenithperrun ;
-        cout << ", Time for Run = " << tObs << "s" << endl;
-
-        // get the bin content of the ON and OFF events per second and scale the off-counts with the alpha factor
-        for( unsigned int k = 0; k < iBin; k++ )
-        {
-            double n_on        = ( t1->GetBinContent( k ) ) * 1.0 / tObs;
-            double n_on_error  = ( t1->GetBinError( k ) ) * 1.0 / tObs;
-            double n_off       = ( t2->GetBinContent( k ) ) * alphaNorm / tObs;
-            double n_off_error = ( t2->GetBinError( k ) ) * alphaNorm / tObs;
-
-            /* ON rate (signal and background) */
-            tRatePerEnergyON->SetPoint( m, t1->GetBinCenter( k ), iZenithperrun, n_on );
-            tRatePerEnergyON->SetPointError( m, 0., 0., n_on_error );
-
-            /* Signal rate (ON minus OFF) */
-            tRatePerEnergySignal->SetPoint( m, t1->GetBinCenter( k ), iZenithperrun, TMath::Abs( n_on - n_off ) );
-            tRatePerEnergySignal->SetPointError( m, 0., 0., ( TMath::Abs( n_on_error + n_off_error ) / 2.0 ) );
-
-            /* OFF rate */
-            tRatePerEnergyOFF->SetPoint( m, t1->GetBinCenter( k ), iZenithperrun, n_off );
-            tRatePerEnergyOFF->SetPointError( m, 0., 0., n_off_error );
-
-            m++;
-        }
-    }
-
-    // write particle number file
-    TFile* SignalRateFile = new TFile( iOutFil.c_str(), "RECREATE" );
-
-    // write graphs of ON rate to file
-    tRatePerEnergyON->Write( "gONRate" );
-    // write graphs of OFF rate to file
-    tRatePerEnergyOFF->Write( "gBGRate" );
-    // write graphs of signal rate to file
-    tRatePerEnergySignal->Write( "gSignalRate" );
-    SignalRateFile->Close();
-
-    return 0;
-}

From 0d52ca14e3d6feeaa752c6b87668be94ee999c0e Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 24 Dec 2025 14:48:49 +0100
Subject: [PATCH 04/16] improved docstrings

---
 src/VTMVAEvaluator.cpp          |  2 +-
 src/calculateCrabRateFromMC.cpp | 31 +++++++++++++++++++------------
 2 files changed, 20 insertions(+), 13 deletions(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 3de92838..e282c57f 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -493,7 +493,7 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
         /////////////////////////////////////////////////////////
         // get optimal signal efficiency (from maximum signal/noise ratio)
         /////////////////////////////////////////////////////////
-        if( fParticleNumberFileName.size() > 0 )
+        if( readNonNoffGraphsFromFile.size() > 0 )
         {
             cout << endl;
             cout << "======================= optimize sensitivity =======================" << endl;
diff --git a/src/calculateCrabRateFromMC.cpp b/src/calculateCrabRateFromMC.cpp
index 369b6d8e..4d20d9ac 100644
--- a/src/calculateCrabRateFromMC.cpp
+++ b/src/calculateCrabRateFromMC.cpp
@@ -1,7 +1,7 @@
 /*! \file VTS.calculateCrabRateFromMC.cpp
-    \brief calculate gamma ray rate from Crab Nebula with effective areas from MC
-
+    \brief calculate gamma-ray and background rates equivalent to Crab Nebula observations
 
+    Use effective areas for signal rates and observations for background rates.
 
 */
 
@@ -26,7 +26,7 @@ using namespace std;
  * fill profile histogram with a new rate entry
  *
  */
-void fill_profilehistogram_for_TMVA(
+void fill_profile_histogram_for_TMVA(
     TProfile2D* fMCRates_TMVA,
     double i_E, double i_Ze, double iRate )
 {
@@ -94,6 +94,12 @@ vector< double > read_zenith_bins( string fEffAreaFile )
     return tmp_zebin_edges;
 }
 
+/*
+    * read energy bin pairs from TMVA run parameter file
+    *
+    * ENERGYBINEDGES -1.5 -0.5 -1. 0. -0.5 0.5 0.0 1. 0.5 2.0
+    * ENERGYBINS -1.50 -0.75 -0.25 0.5 2.0
+*/
 vector< pair< double, double > > read_energy_minmax_pairs(
     vector< double > tmp_ebins_histo, string iEnergyKeyWord )
 {
@@ -125,9 +131,7 @@ vector< pair< double, double > > read_energy_minmax_pairs(
  * ENERGYBINEDGES -1.5 -0.5 -1. 0. -0.5 0.5 0.0 1. 0.5 2.0
  * ENERGYBINS -1.50 -0.75 -0.25 0.5 2.0
  */
-vector< double > read_energy_bins(
-    string iTMVAParameterFile,
-    string iEnergyKeyWord )
+vector< double > read_energy_bins(string iTMVAParameterFile, string iEnergyKeyWord )
 {
     vector< double > tmp_e;
     ifstream is( iTMVAParameterFile.c_str(), ifstream::in );
@@ -244,7 +248,7 @@ TTree* fillMCRates(
     float pedvar = 0.;
     float MCrate;
 
-    TTree* fMC = new TTree( "fMCRate", "MC rate predictions" );
+    TTree* fMC = new TTree( "fMCRate", "MC rate predictions (1/min)" );
     fMC->Branch( "ze", &ze, "ze/F" );
     fMC->Branch( "az", &az, "az/I" );
     fMC->Branch( "Woff", &Woff, "Woff/F" );
@@ -279,7 +283,7 @@ TTree* fillMCRates(
     t->SetBranchAddress( "eff", t_eff );
 
     // rate calculator
-    VMonteCarloRateCalculator* fMCR = new VMonteCarloRateCalculator();
+    VMonteCarloRateCalculator fMCR;
     // hardwired Whipple spectrum
     double fWhippleNorm = 3.250e-11;
     double fWhippleIndex = -2.440;
@@ -305,7 +309,7 @@ TTree* fillMCRates(
             feffectivearea.push_back( t_eff[e] );
         }
 
-        MCrate = fMCR->getMonteCarloRate(
+        MCrate = fMCR.getMonteCarloRate(
                      fenergy, feffectivearea,
                      fWhippleIndex, fWhippleNorm, 1., fEnergyThreshold, 1.e7, false );
         MCrate *= ( 1. - fDeadTime / 100. );
@@ -318,11 +322,11 @@ TTree* fillMCRates(
 
         for( unsigned int e = 0; e < ebins.size(); e++ )
         {
-            fill_profilehistogram_for_TMVA(
+            fill_profile_histogram_for_TMVA(
                 h,
                 0.5 * ( ebins[e].first + ebins[e].second ),
                 ze,
-                fMCR->getMonteCarloRate(
+                fMCR.getMonteCarloRate(
                     fenergy, feffectivearea,
                     fWhippleIndex, fWhippleNorm, 1.,
                     TMath::Power( 10., ebins[e].first ),
@@ -447,6 +451,8 @@ void fillBackgroundRates( string iRunList, TProfile2D* h, vector< pair< double,
 /*
  * convert 2D histograms to 2D graphs
  *
+ * convert rate from 1/min to 1/s
+ *
  */
 TGraph2DErrors* getTGraph2D( vector< TProfile2D* > h, string iGraphName )
 {
@@ -470,7 +476,7 @@ TGraph2DErrors* getTGraph2D( vector< TProfile2D* > h, string iGraphName )
             {
                 i_z += h[0]->GetBinContent( b_x, b_y );
             }
-            // expect rate graphs in 1./s
+            // convert expect rate graphs in 1./s
             i_g->SetPoint(
                 z,
                 h[1]->GetXaxis()->GetBinCenter( b_x ),
@@ -531,6 +537,7 @@ int main( int argc, char* argv[] )
     {
         tmp_ebins = read_energy_minmax_pairs( tmp_ebins_histo, "ENERGYBINS" );
     }
+    // ENERGYBINEDGES for VTS analysis
     else
     {
         tmp_ebins = read_energy_minmax_pairs(

From 7c14c32e67fe7ed4b47f6bc7f5a9d37ffeb3d568 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 24 Dec 2025 14:56:32 +0100
Subject: [PATCH 05/16] airmass interpolation

---
 docs/changes/336.feature.md |  1 +
 src/VTMVAEvaluator.cpp      | 12 +++++++++---
 2 files changed, 10 insertions(+), 3 deletions(-)
 create mode 100644 docs/changes/336.feature.md

diff --git a/docs/changes/336.feature.md b/docs/changes/336.feature.md
new file mode 100644
index 00000000..067eaec9
--- /dev/null
+++ b/docs/changes/336.feature.md
@@ -0,0 +1 @@
+Improve TMVA interpolation from average zenith to average airmass.
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index e282c57f..35763acc 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -2136,8 +2136,10 @@ TGraph* VTMVAEvaluator::fillfromGraph2D( TObject* i_G, double i_ze_min, double i
         TGraph2D* iG2D = ( TGraph2D* )i_G;
         TGraph* iG1D = new TGraph( iG2D->GetN() );
         Double_t* x = iG2D->GetX();
-        // TEMP: weight by cos?
-        double ze_mean = 0.5 * ( i_ze_min + i_ze_max );
+        double z1 = i_ze_min * TMath::DegToRad();
+        double z2 = i_ze_max * TMath::DegToRad();
+        double avg_airmass = 0.5 * ( (1.0 / TMath::Cos(z1)) + (1.0 / TMath::Cos(z2)) );
+        double ze_mean = TMath::ACos(1.0 / avg_airmass) * TMath::RadToDeg();
 
         for( int i = 0; i < iG2D->GetN(); i++ )
         {
@@ -2180,7 +2182,11 @@ void VTMVAEvaluator::calculate_average_zenith_angle()
     {
         fAverageZenithPerRun = 90. - i_ze / i_n;
     }
-    cout << "VTMVAEvaluator: average zenith " << fAverageZenithPerRun << endl;
+    else
+    {
+        fAverageZenithPerRun = -1.;
+    }
+    cout << "VTMVAEvaluator: average zenith estimate " << fAverageZenithPerRun << " deg" << endl;
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////

From f4bd09d72c5ffd21b0d77ac23aea567d888a2db4 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 12:08:52 +0100
Subject: [PATCH 06/16] cleanup

---
 inc/VAnaSumRunParameter.h                 |   5 +-
 inc/VTMVAEvaluator.h                      |  14 +-
 inc/VTMVARunData.h                        |   4 +-
 src/VAnaSumRunParameter.cpp               |   5 -
 src/VGammaHadronCuts.cpp                  |   8 +-
 src/VStereoMaps.cpp                       |  48 +---
 src/VTMVAEvaluator.cpp                    | 274 +---------------------
 src/VTMVARunData.cpp                      |   2 +-
 src/trainTMVAforGammaHadronSeparation.cpp |   4 +-
 9 files changed, 25 insertions(+), 339 deletions(-)

diff --git a/inc/VAnaSumRunParameter.h b/inc/VAnaSumRunParameter.h
index fd2f11e7..aa66c41e 100644
--- a/inc/VAnaSumRunParameter.h
+++ b/inc/VAnaSumRunParameter.h
@@ -110,9 +110,6 @@ class VAnaSumRunParameterDataClass : public TNamed
 
         string fEffectiveAreaFile;                // file with effective areas, use NOFILE if not available
 
-        // all models
-        unsigned int fNBoxSmooth;
-
         // ON/OFF MODEL
         double fOO_alpha;
 
@@ -139,7 +136,7 @@ class VAnaSumRunParameterDataClass : public TNamed
         {
             return fRunOn < x.fRunOn;
         }
-        ClassDef( VAnaSumRunParameterDataClass, 3 );
+        ClassDef( VAnaSumRunParameterDataClass, 4 );
 };
 
 class VAnaSumRunParameter : public TNamed, public VGlobalRunParameter
diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h
index 589c2cac..3099b4a4 100644
--- a/inc/VTMVAEvaluator.h
+++ b/inc/VTMVAEvaluator.h
@@ -118,7 +118,6 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         double                  fOptimizationObservingTime_h;
 
         double   fTMVA_EvaluationResult;         // result from TMVA evaluator
-        bool     fsmoothAndInterpolateMVAValues;
         double   fAverageZenithPerRun;           // (rough) average zenith angle of run
 
         string   fTMVAMethodName;
@@ -169,11 +168,6 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         double           interpolate_mva_evaluation();
         TGraph*          readNonNoffGraphsFromFile( TFile* iF, double i_ze_min, double i_ze_max, bool bIsOn = true );
         void             reset();
-        void             smoothAndInterpolateMVAValue( TH1F*, TH1F*,
-                unsigned int iE_min, unsigned int iE_max, unsigned int iZ_min, unsigned int iZ_max );
-        void             smoothAndInterpolateMVAValue_EnergyOnly( TH1F*, TH1F* );
-        void             smoothAndInterpolateMVAValue_Energy_and_Zenith( TH1F*, TH1F*,
-                unsigned int iE_min, unsigned int iE_max, unsigned int iZ_min, unsigned int iZ_max );
 
     public:
 
@@ -181,8 +175,6 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         ~VTMVAEvaluator() {};
 
         bool    evaluate( bool interpolate_mva = false, bool use_average_zenith_angle = true );
-        vector< double > getBackgroundEfficiency();
-        vector< bool >   getOptimumCutValueFound();
         vector< double > getSignalEfficiency();
         double  getOptimalTheta2Cut( double iEnergy_log10TeV, double iZe = -9999 );
         vector< double > getTMVACutValue();
@@ -237,10 +229,6 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         }
         void   setSignalEfficiency( double iSignalEfficiency = -99. );
         void   setSignalEfficiency( map< unsigned int, double > iMSignalEfficiency );
-        void   setsmoothAndInterpolateMVAValues( bool iS = true )
-        {
-            fsmoothAndInterpolateMVAValues = iS;
-        }
         void   setSpectralIndexForEnergyWeighting( double iS = -2. )
         {
             fSpectralIndexForEnergyWeighting = iS;
@@ -253,7 +241,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         }
         void   setTMVAMethod( string iMethodName = "BDT" );
 
-        ClassDef( VTMVAEvaluator, 36 );
+        ClassDef( VTMVAEvaluator, 37 );
 };
 
 #endif
diff --git a/inc/VTMVARunData.h b/inc/VTMVARunData.h
index 613e615b..cc5a09ec 100644
--- a/inc/VTMVARunData.h
+++ b/inc/VTMVARunData.h
@@ -99,7 +99,7 @@ class VTMVARunData : public TNamed
         VTMVARunData();
         ~VTMVARunData() {}
         void print();
-        VTableLookupRunParameter* getTLRunParameter();
+        VTableLookupRunParameter* getTableLookupRunParameters();
         bool readConfigurationFile( char* );
         bool openDataFiles();
         void setDebug( bool iB = true )
@@ -112,7 +112,7 @@ class VTMVARunData : public TNamed
         }
         void updateTrainingEvents( string iVarName, unsigned int iNEvents );
 
-        ClassDef( VTMVARunData, 12 );
+        ClassDef( VTMVARunData, 13 );
 };
 
 #endif
diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp
index 84ec1a27..151a72d5 100644
--- a/src/VAnaSumRunParameter.cpp
+++ b/src/VAnaSumRunParameter.cpp
@@ -62,9 +62,6 @@ VAnaSumRunParameterDataClass::VAnaSumRunParameterDataClass()
 
     fEffectiveAreaFile = "";                      // file with effective areas, use NOFILE if not available
 
-    // smoothing algorithm (don't use it if you don't know it)
-    fNBoxSmooth = 0;
-
     // ON/OFF MODEL
     fOO_alpha = 0.;
 
@@ -1243,7 +1240,6 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i )
             cout << "TEMPLATE BACKGROUND MODEL" << endl;
         }
     }
-    cout << "\t NBoxsmooth: " << fRunList[i].fNBoxSmooth << endl;
     cout << endl;
 }
 
@@ -1364,7 +1360,6 @@ void VAnaSumRunParameter::reset( VAnaSumRunParameterDataClass it )
 
     it.fAcceptanceFile = "";
 
-    it.fNBoxSmooth = 0;
     it.fOO_alpha = 0.;
 
     it.fRM_RingRadius = 0.;
diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp
index 9c1f325a..816067b7 100644
--- a/src/VGammaHadronCuts.cpp
+++ b/src/VGammaHadronCuts.cpp
@@ -1065,7 +1065,7 @@ bool VGammaHadronCuts::applyTMVACut( int i )
     {
         cout << "VGammaHadronCuts::applyTMVACut event " << i;
         cout << ", signal efficiency " << fTMVASignalEfficiency.size();
-        cout << ", probability threshold/MVA cut " << fTMVA_MVACut.size();
+        cout << ", MVA cut " << fTMVA_MVACut.size();
         cout << " (" << fTMVAEvaluator << ")";
         cout << endl;
     }
@@ -1306,8 +1306,6 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile,
     fTMVAEvaluator = new VTMVAEvaluator();
 
     fTMVAEvaluator->setDebug( fDebug );
-    // smoothing of MVA values
-    fTMVAEvaluator->setsmoothAndInterpolateMVAValues( true );
     // constant signal efficiency
     if( fTMVASignalEfficiency.size() > 0 )
     {
@@ -1321,8 +1319,8 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile,
     else
     {
         cout << "VGammaHadronCuts::initTMVAEvaluator error: unclear TMVA cut settings" << endl;
-        cout << "\t TMVASignalEfficiency: " << fTMVASignalEfficiency.size() << endl;
-        cout << "\t TMVAProbabilityThreshold: " << fTMVA_MVACut.size() << endl;
+        cout << "\t MVA Signal Efficiency: " << fTMVASignalEfficiency.size() << endl;
+        cout << "\t MVA Threshold: " << fTMVA_MVACut.size() << endl;
         cout << "exiting... " << endl;
         exit( EXIT_FAILURE );
     }
diff --git a/src/VStereoMaps.cpp b/src/VStereoMaps.cpp
index f8588a36..d08c6c0c 100644
--- a/src/VStereoMaps.cpp
+++ b/src/VStereoMaps.cpp
@@ -388,50 +388,18 @@ void VStereoMaps::makeTwoDStereo_BoxSmooth( double i_xderot, double i_yderot, do
                 continue;
             }
 
-            // additional smoothing to get ride of edge feature (not a default feature!)
-            // get filling factor (fraction of bin in theta2 circle)
-            if( fRunList.fNBoxSmooth > 0 )
+            // theta2 cut
+            i_r = sqrt(( i_xderot - i_xbin ) * ( i_xderot - i_xbin ) + ( i_yderot - i_ybin ) * ( i_yderot - i_ybin ) );
+            if( i_r <= thetaCutMax )
             {
-                double p = 0.;
-                for( unsigned int f = 0; f < fRunList.fNBoxSmooth; f++ )
+                hmap_stereo->Fill( i_xbin, i_ybin );
+                hmap_alpha->Fill( i_xbin, i_ybin, i_weight );
+                if( hmap_ratio )
                 {
-                    i_xbin = fRandom->Uniform( hmap_stereo->GetXaxis()->GetBinLowEdge( i + 1 ), hmap_stereo->GetXaxis()->GetBinUpEdge( i + 1 ) );
-                    i_ybin = fRandom->Uniform( hmap_stereo->GetYaxis()->GetBinLowEdge( j + 1 ), hmap_stereo->GetYaxis()->GetBinUpEdge( j + 1 ) );
-                    i_r = sqrt(( i_xderot - i_xbin ) * ( i_xderot - i_xbin ) + ( i_yderot - i_ybin ) * ( i_yderot - i_ybin ) );
-                    if( i_r <= thetaCutMax )
-                    {
-                        p++;
-                    }
-                }
-                p /= ( double )( fRunList.fNBoxSmooth );
-                i_xbin = hmap_stereo->GetXaxis()->GetBinCenter( i + 1 );
-                i_ybin = hmap_stereo->GetYaxis()->GetBinCenter( j + 1 );
-                if( fRandom->Uniform() < p )
-                {
-                    hmap_stereo->Fill( i_xbin, i_ybin );
-                    hmap_alpha->Fill( i_xbin, i_ybin, i_weight );
-                    if( hmap_ratio )
-                    {
-                        hmap_ratio->Fill( i_MeanSignalBackgroundAreaRatio );
-                    }
+                    hmap_ratio->Fill( i_MeanSignalBackgroundAreaRatio );
                 }
             }
-            // default: filling of the sky map at the corresponding bin
-            else
-            {
-                // theta2 cut
-                i_r = sqrt(( i_xderot - i_xbin ) * ( i_xderot - i_xbin ) + ( i_yderot - i_ybin ) * ( i_yderot - i_ybin ) );
-                if( i_r <= thetaCutMax )
-                {
-                    hmap_stereo->Fill( i_xbin, i_ybin );
-                    hmap_alpha->Fill( i_xbin, i_ybin, i_weight );
-                    if( hmap_ratio )
-                    {
-                        hmap_ratio->Fill( i_MeanSignalBackgroundAreaRatio );
-                    }
-                }
-            }
-        }
+    }
     }
 }
 
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 35763acc..fc9c1a01 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -51,8 +51,6 @@ void VTMVAEvaluator::reset()
     setTMVAErrorFraction();
     fTMVA_EvaluationResult = -99.;
     fTMVACutValueNoVec = -99.;
-
-    setsmoothAndInterpolateMVAValues();
 }
 
 /*
@@ -235,8 +233,8 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
                         bGoodRun = false;
                     }
                 }
-                // allow that first files are missing
-                // (this happens when there are no training events in the first energy or zenith bins)
+                // allow that first or last files are missing
+                // (this happens when there are no training events in these energy or zenith bins)
                 if(!bGoodRun )
                 {
                     if( i == iMinMissingBin || j == jMinMissingBin )
@@ -349,8 +347,8 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
 
                 iF.Close();
             }
-        }//end loop on zenith bins
-    }//end loop on energy bins
+        } // end loop on zenith bins
+    } // end loop on energy bins
 
     // after this stage, there should be no energy/zenith bins (both of them are combined)
 
@@ -366,18 +364,16 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
 
     //////////////////////////////////////////////////////////////////////////////////////
     // create and initialize TMVA readers
-    // loop over all  energy bins: open one weight (XML) file per energy bin
-    // looping over spectral energy and zenith angle bins
+    // loop over all energy / zenith bins: open one weight (XML) file per bin
     for( unsigned int b = 0; b < fTMVAData.size(); b++ )
     {
         fTMVAData[b]->fTMVAReader = new TMVA::Reader();
         if( fDebug )
         {
-            cout << "INITIALIZE TMVA file: " << fTMVAData[b]->fTMVAFileName << endl;
+            cout << "Initialize TMVA file: " << fTMVAData[b]->fTMVAFileName << endl;
         }
         //////////////////////////////////////////
         // set TMVA cut value
-        // (optimization later)
 
         // fixed signal efficiency
         if( fTMVACutValueNoVec < -1. && fSignalEfficiencyNoVec > 0. )
@@ -393,11 +389,9 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
         // no optimization took place
         fTMVAData[b]->fTMVAOptimumCutValueFound = false;
 
-        // weight file for this energy bin
-
         if( fDebug )
         {
-            cout << "reading TMVA XML weight file: " << fTMVAData[b]->fTMVAFileNameXML << endl;
+            cout << "Reading TMVA XML weight file: " << fTMVAData[b]->fTMVAFileNameXML << endl;
         }
 
         // get list of training variables
@@ -508,14 +502,6 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
 
     }
 
-    // smooth and interpolate
-    if( fParticleNumberFileName.size() > 0 && fsmoothAndInterpolateMVAValues )
-    {
-        smoothAndInterpolateMVAValue( 0, 0,
-                                      iWeightFileIndex_Emin, iWeightFileIndex_Emax,
-                                      iWeightFileIndex_Zmin, iWeightFileIndex_Zmax );
-    }
-
     // print some info to screen
     cout << "VTMVAEvaluator: Initialized " << fTMVAData.size() << " MVA readers " << endl;
 
@@ -1575,227 +1561,6 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
     return true;
 }
 
-/*
- smoothing of optimal cut value vs energy curves
-
- (energy axis only)
-
-*/
-void VTMVAEvaluator::smoothAndInterpolateMVAValue_EnergyOnly(
-    TH1F* effS, TH1F* effB )
-{
-    int z = 0;
-    // fill graph to be smoothed
-    TGraph* iG = new TGraph( 1 );
-    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
-    {
-        if( fTMVAData[i] )
-        {
-            if( fTMVAData[i]->fTMVAOptimumCutValueFound )
-            {
-                iG->SetPoint( z, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fTMVACutValue );
-                z++;
-            }
-        }
-    }
-    // smooth graph
-    TGraph* iGout = new TGraph( 1 );
-    TGraphSmooth* iGSmooth = new TGraphSmooth( "t" );
-    iGout = ( TGraph* )iGSmooth->SmoothKern( iG, "normal", 0.5, 100 );
-
-    // fill smoothed and interpolated values into MVA vector
-    // set all points to 'optimized'
-    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
-    {
-        if( fTMVAData[i] )
-        {
-            cout << "\t TMVA values: unsmoothed at ";
-            cout << TMath::Power( 10., fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
-            cout << " TeV, \t" << fTMVAData[i]->fTMVACutValue;
-            fTMVAData[i]->fTMVACutValue = iGout->Eval( fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
-            cout << ", smoothed " << fTMVAData[i]->fTMVACutValue;
-            if(!fTMVAData[i]->fTMVAOptimumCutValueFound )
-            {
-                cout << " (interpolated non-optimal value)";
-            }
-            cout << " (" << i << ")" << endl;
-            fTMVAData[i]->fTMVAOptimumCutValueFound = true;
-
-            // get efficiency histograms
-            TFile iTMVAFile( fTMVAData[i]->fTMVAFileName.c_str() );
-            TH1F* effS = getEfficiencyHistogram( "effS", &iTMVAFile, fTMVAData[i]->fTMVAMethodTag_2 );
-            TH1F* effB = getEfficiencyHistogram( "effB", &iTMVAFile, fTMVAData[i]->fTMVAMethodTag_2 );
-
-            if( effS )
-            {
-                fTMVAData[i]->fSignalEfficiency = effS->GetBinContent( effS->FindBin( fTMVAData[i]->fTMVACutValue ) );
-            }
-            if( effB )
-            {
-                fTMVAData[i]->fBackgroundEfficiency = effB->GetBinContent( effB->FindBin( fTMVAData[i]->fTMVACutValue ) );
-                // background efficiency might be zero -> fill it with first non-zero value
-                if( fTMVAData[i]->fBackgroundEfficiency < 1.e-8 )
-                {
-                    int iS = effB->FindBin( fTMVAData[i]->fTMVACutValue );
-                    for( int j = iS; j > 0; j-- )
-                    {
-                        if( effB->GetBinContent( j ) > 0. )
-                        {
-                            fTMVAData[i]->fBackgroundEfficiency = effB->GetBinContent( j );
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-/*
-
- smoothing of optimal cut value vs energy curves
-
- (energy and zenith angle dependent)
-
-*/
-void VTMVAEvaluator::smoothAndInterpolateMVAValue_Energy_and_Zenith(
-    TH1F* effS, TH1F* effB,
-    unsigned int iWeightFileIndex_Emin, unsigned int iWeightFileIndex_Emax,
-    unsigned int iWeightFileIndex_Zmin, unsigned int iWeightFileIndex_Zmax )
-{
-
-    unsigned int Ebins = iWeightFileIndex_Emax - iWeightFileIndex_Emin + 1;
-    unsigned int ZEbins = iWeightFileIndex_Zmax - iWeightFileIndex_Zmin + 1;
-
-    // 2D histogram to be smoothed (energy and zenith angle dependence)
-    TH2D* iH2 = new TH2D( "h1", "Smooth cut values",
-                          Ebins, iWeightFileIndex_Emin, iWeightFileIndex_Emax,
-                          ZEbins, iWeightFileIndex_Zmin, iWeightFileIndex_Zmax );
-    Double_t effS_value[Ebins * ZEbins];
-    for( unsigned int l = 0; l < fTMVAData.size(); l++ )
-    {
-        if( fTMVAData[l] )
-        {
-            TFile iTMVAFile( fTMVAData[l]->fTMVAFileName.c_str() );
-            if( iTMVAFile.IsZombie() )
-            {
-                cout << "VTMVAEvaluator: error reading file with cut efficiencies: ";
-                cout << fTMVAData[l]->fTMVAFileName << endl;
-                cout << "Exiting..." << endl;
-                exit( EXIT_FAILURE );
-            }
-            TH1F* effS = getEfficiencyHistogram( "effS", &iTMVAFile, fTMVAData[l]->fTMVAMethodTag_2 );
-            effS_value[l] = effS->GetBinContent( effS->FindBin( fTMVAData[l]->fTMVACutValue ) );
-        }
-    }
-    for( unsigned int z = 0; z < fTMVAData.size(); z++ )
-    {
-        if(!fTMVAData[z] )
-        {
-            continue;
-            cout << "VTMVAEvaluator: error reading file with cut efficiencies for bins: ";
-            cout << ", entry " << z << endl;
-            cout << "Exiting..." << endl;
-            exit( EXIT_FAILURE );
-        }
-        // get signal efficiency histograms and cut values
-        TFile iTMVAFile( fTMVAData[z]->fTMVAFileName.c_str() );
-        TH1F* effS = getEfficiencyHistogram( "effS", &iTMVAFile, fTMVAData[z]->fTMVAMethodTag_2 );
-        if( fTMVAData[z]->fTMVAOptimumCutValueFound )
-        {
-            iH2->SetBinContent( fTMVAData[z]->fEnergyBin, fTMVAData[z]->fZenithBin, fTMVAData[z]->fTMVACutValue );
-        }
-        // bins without optimal cut value and not in highest energy bin
-        else if(!fTMVAData[z]->fTMVAOptimumCutValueFound
-                && ( fTMVAData[z]->fZenithBin < ZEbins )
-                && ( fTMVAData[z]->fEnergyBin != Ebins - 1 ) )
-        {
-            for( int k = 0; k < effS->GetNbinsX(); k++ )
-            {
-                // search for similar cut efficiency in neighbouring bin
-                unsigned int i_alt_index = z + 1;
-                if( z > 0 )
-                {
-                    i_alt_index = z - 1;
-                }
-                if( i_alt_index < fTMVAData.size()
-                        && TMath::Abs( effS->GetBinContent( k ) - effS_value[i_alt_index] ) < 0.001 )
-                {
-                    fTMVAData[z]->fTMVACutValue = effS->GetBinCenter( k );
-                    effS_value[z] = effS->GetBinContent( k );
-                }
-            }
-        }
-        // bins without optimal cut value and in highest energy bin
-        else if(!fTMVAData[z]->fTMVAOptimumCutValueFound
-                && ( fTMVAData[z]->fZenithBin < ZEbins ) && ( fTMVAData[z]->fEnergyBin == Ebins - 1 ) )
-        {
-            for( unsigned int l = 0; l < ZEbins; l++ )
-            {
-                if( effS_value[fTMVAData.size() - l] > 1.e-10 )
-                {
-                    for( int k = 0; k < effS->GetNbinsX(); k++ )
-                    {
-                        if( TMath::Abs( effS->GetBinContent( k ) - effS_value[fTMVAData.size() - l] ) < 0.0001 )
-                        {
-                            fTMVAData[z]->fTMVACutValue = effS->GetBinCenter( k );
-                            effS_value[z] = effS->GetBinContent( k );
-                        }
-                    }
-                }
-            }
-            if( fTMVAData[z]->fTMVACutValue == -99 )
-            {
-                cout << "Error: no optimal cut value found for this bin" << endl;
-            }
-            cout << "\t TMVA values: at " << TMath::Power( 10., fTMVAData[z]->fSpectralWeightedMeanEnergy_Log10TeV );
-            cout << " TeV, \t" << fTMVAData[z]->fTMVACutValue;
-            if(!fTMVAData[z]->fTMVAOptimumCutValueFound )
-            {
-                cout << " (interpolated non-optimal value)";
-            }
-            cout << ", signal efficiency: " << effS_value[z];
-            cout << " (" << z << ")" << endl;
-        }
-    }
-}
-
-/*
-
- smoothing of optimal cut value vs energy curves
-
- missing (non-optimized) are interpolated
-
- note: signal and background efficiencies are not updated
-
-*/
-void VTMVAEvaluator::smoothAndInterpolateMVAValue(
-    TH1F* effS, TH1F* effB,
-    unsigned int iWeightFileIndex_Emin, unsigned int iWeightFileIndex_Emax,
-    unsigned int iWeightFileIndex_Zmin, unsigned int iWeightFileIndex_Zmax )
-{
-    if( fTMVAData.size() == 0 )
-    {
-        return;
-    }
-
-    cout << "Smooth and interpolate MVA cut values" << endl;
-
-    //////////////////////////////////////////////
-    // energy dependent TMVA cut optimization only
-    if( iWeightFileIndex_Zmax == iWeightFileIndex_Zmin )
-    {
-        smoothAndInterpolateMVAValue_EnergyOnly( effS, effB );
-    }
-    // energy and zenith angle dependent
-    else
-    {
-        smoothAndInterpolateMVAValue_Energy_and_Zenith(
-            effS, effB,
-            iWeightFileIndex_Emin, iWeightFileIndex_Emax,
-            iWeightFileIndex_Zmin, iWeightFileIndex_Zmax );
-    }
-}
 
 void VTMVAEvaluator::plotEfficiencyPlotsPerBin( unsigned int iBin, TGraph* iGSignal_to_sqrtNoise,
         TGraph* iGSignal_to_sqrtNoise_Smooth, TH1F* hEffS, TH1F* hEffB,
@@ -2030,7 +1795,6 @@ double VTMVAEvaluator::getTMVACutValue( unsigned int iEnergyBin, double iE_min_l
 }
 
 
-
 void VTMVAEvaluator::printSensitivityOptimizationParameters()
 {
     cout << "VTMAEvaluator: MVA cut parameter is optimized for: " << endl;
@@ -2042,18 +1806,6 @@ void VTMVAEvaluator::printSensitivityOptimizationParameters()
     cout << "\t" << fOptimizationMinSourceStrength << " minimum source strength" << endl;
 }
 
-vector< double > VTMVAEvaluator::getBackgroundEfficiency()
-{
-    vector< double > iA;
-    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
-    {
-        if( fTMVAData[i] )
-        {
-            iA.push_back( fTMVAData[i]->fBackgroundEfficiency );
-        }
-    }
-    return iA;
-}
 
 vector< double > VTMVAEvaluator::getSignalEfficiency()
 {
@@ -2081,18 +1833,6 @@ vector< double > VTMVAEvaluator::getTMVACutValue()
     return iA;
 }
 
-vector< bool > VTMVAEvaluator::getOptimumCutValueFound()
-{
-    vector< bool > iA;
-    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
-    {
-        if( fTMVAData[i] )
-        {
-            iA.push_back( fTMVAData[i]->fTMVAOptimumCutValueFound );
-        }
-    }
-    return iA;
-}
 
 /*
  * read graph with on/off numbers as a function of energy and possibly zenith angle
diff --git a/src/VTMVARunData.cpp b/src/VTMVARunData.cpp
index 745cc15f..11877d2e 100644
--- a/src/VTMVARunData.cpp
+++ b/src/VTMVARunData.cpp
@@ -769,7 +769,7 @@ bool VTMVARunData::readConfigurationFile( char* iC )
     return true;
 }
 
-VTableLookupRunParameter* VTMVARunData::getTLRunParameter()
+VTableLookupRunParameter* VTMVARunData::getTableLookupRunParameters()
 {
     TDirectory* iG_CurrentDirectory = gDirectory;
     if( fSignalFileName.size() > 0 )
diff --git a/src/trainTMVAforGammaHadronSeparation.cpp b/src/trainTMVAforGammaHadronSeparation.cpp
index 466aadda..588ad960 100644
--- a/src/trainTMVAforGammaHadronSeparation.cpp
+++ b/src/trainTMVAforGammaHadronSeparation.cpp
@@ -295,9 +295,9 @@ bool train( VTMVARunData* iRun,
         }
         iSignalTree_reduced->Write();
         iBackgroundTree_reduced->Write();
-        if( iRun->getTLRunParameter() )
+        if( iRun->getTableLookupRunParameters() )
         {
-            iRun->getTLRunParameter()->Write();
+            iRun->getTableLookupRunParameters()->Write();
         }
         cout << "Writing reduced event lists for training: ";
         cout << gDirectory->GetName() << endl;

From e33b0b4f9767f467ecf7a57c1c831c5d60a37f7c Mon Sep 17 00:00:00 2001
From: GernotMaier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 12:13:47 +0100
Subject: [PATCH 07/16] typo

---
 src/VTMVAEvaluator.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index fc9c1a01..510f5757 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -487,7 +487,7 @@ bool VTMVAEvaluator::initializeWeightFiles( string iWeightFileName,
         /////////////////////////////////////////////////////////
         // get optimal signal efficiency (from maximum signal/noise ratio)
         /////////////////////////////////////////////////////////
-        if( readNonNoffGraphsFromFile.size() > 0 )
+        if( fParticleNumberFileName.size() > 0 )
         {
             cout << endl;
             cout << "======================= optimize sensitivity =======================" << endl;

From fb5d36ff4124cdbfa693f4768200e1fa8a2ec125 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 16:14:11 +0100
Subject: [PATCH 08/16] readability

---
 src/VTMVAEvaluator.cpp | 106 +++++++++++++++++++++--------------------
 1 file changed, 54 insertions(+), 52 deletions(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 510f5757..c739ae8a 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -1166,9 +1166,10 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
     }
 
     //////////////////////////////////////////////////////
-    // get mean energy of the considered bins
+    // get spectral weighted mean energy of the considered bins
     // interval [fTMVAData[iDataBin]->fEnergyCut_Log10TeV_min,fTMVAData[iDataBin]->fEnergyCut_Log10TeV_max]
     // make sure that energy is not lower or higher then minimum/maximum bins in the rate graphs
+    // This is the energy at which the event numbers are calculated.
 
     double x = 0.;
     double p = 0.;
@@ -1215,7 +1216,6 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
     // Interpolate between the values of the TGraph2D
     //
     // Convert the observing time in seconds as the particle rate is given in 1/seconds
-    // Get the value of the middle of the energy and zenith angle bin
     Non = i_on->Eval( fTMVAData[iDataBin]->fSpectralWeightedMeanEnergy_Log10TeV )
           * fOptimizationObservingTime_h * fParticleNumberFile_Conversion_Rate_to_seconds;
     Nof = i_of->Eval( fTMVAData[iDataBin]->fSpectralWeightedMeanEnergy_Log10TeV )
@@ -1280,16 +1280,18 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
 
     //////////////////////////////////////////////////////
     // loop over different source strengths (in Crab Units)
-    // (up to 30 CU)
+    // start from low values to high values; optimal cut value
+    // is found when we reach fir time the required significance
+    // source strength steps on log scale (up to 30 CU)
     unsigned int iSourceStrengthStepSizeN =
         ( unsigned int )(( log10( 30. ) - log10( fOptimizationMinSourceStrength ) ) / 0.005 );
-    cout << "VTVMAEvaluator::optimizeSensitivity(), source strength steps: " << iSourceStrengthStepSizeN << endl;
+    cout << "VTVMAEvaluator::optimizeSensitivity source strength steps: " << iSourceStrengthStepSizeN << endl;
     for( unsigned int s = 0; s < iSourceStrengthStepSizeN; s++ )
     {
         double iSourceStrength = log10( fOptimizationMinSourceStrength ) + s * 0.005;
         iSourceStrength = TMath::Power( 10., iSourceStrength );
 
-        // source (excess) events
+        // excess events
         Ndif = ( Non - Nof ) * iSourceStrength;
 
         // first quick pass to see if there is a change of reaching the required fOptimizationSourceSignificance
@@ -1343,14 +1345,17 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
         // loop over all signal efficiency bins
         for( int i = 1; i < effS->GetNbinsX(); i++ )
         {
-            if( effB->GetBinContent( i ) > 0. && Nof > 0. )
+            double signalEff = effS->GetBinContent( i );
+            double signalEff_mva = effS->GetBinCenter( i );
+            double backEff   = effB->GetBinContent( i );
+            if( backEff > 0. && Nof > 0. )
             {
                 if( fOptimizationBackgroundAlpha > 0. )
                 {
-                    // optimize signal/sqrt(noise)
+                    // optimize signal/sqrt(noise) - Li & Ma significance
                     i_Signal_to_sqrtNoise = VStatistics::calcSignificance(
-                                                effS->GetBinContent( i ) * Ndif + effB->GetBinContent( i ) * Nof,
-                                                effB->GetBinContent( i ) * Nof / fOptimizationBackgroundAlpha,
+                                                signalEff * Ndif + backEff * Nof,
+                                                backEff * Nof / fOptimizationBackgroundAlpha,
                                                 fOptimizationBackgroundAlpha );
                 }
                 else
@@ -1360,22 +1365,22 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
                 if( fDebug )
                 {
                     cout << "___________________________________________________________" << endl;
-                    cout << i << "\t" << Non << "\t" << effS->GetBinContent( i )  << "\t";
-                    cout << Nof << "\t" << effB->GetBinContent( i ) << "\t";
+                    cout << i << "\t" << Non << "\t" << signalEff  << "\t";
+                    cout << Nof << "\t" << backEff << "\t";
                     cout << Ndif << endl;
-                    cout << "\t" << effS->GetBinContent( i ) * Ndif;
-                    cout << "\t" << effS->GetBinContent( i ) * Ndif + effB->GetBinContent( i ) * Nof;
-                    cout << "\t" << effS->GetBinContent( i ) * Non + effB->GetBinContent( i ) * Nof;
-                    cout << "\t" << effB->GetBinContent( i ) * Nof << endl;
+                    cout << "\t" << signalEff * Ndif;
+                    cout << "\t" << signalEff * Ndif + backEff * Nof;
+                    cout << "\t" << signalEff * Non + backEff * Nof;
+                    cout << "\t" << backEff * Nof << endl;
                 }
-                if( effS->GetBinContent( i ) * Ndif > 0. )
+                if( signalEff * Ndif > 0. )
                 {
-                    iGSignalEvents->SetPoint( z_SB, effS->GetBinCenter( i ),  effS->GetBinContent( i ) * Ndif );
-                    iGBackgroundEvents->SetPoint( z_SB, effS->GetBinCenter( i ), effB->GetBinContent( i ) * Nof );
+                    iGSignalEvents->SetPoint( z_SB, signalEff_mva,  signalEff * Ndif );
+                    iGBackgroundEvents->SetPoint( z_SB, signalEff_mva, backEff * Nof );
                     z_SB++;
                 }
                 // check that a minimum number of off events is available
-                if( effB->GetBinContent( i ) * Nof < fOptimizationMinBackGroundEvents )
+                if( backEff * Nof < fOptimizationMinBackGroundEvents )
                 {
                     if( fDebug )
                     {
@@ -1387,10 +1392,10 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
                 // add results to a graph
                 if( iGSignal_to_sqrtNoise && i_Signal_to_sqrtNoise > 1.e-2 )
                 {
-                    iGSignal_to_sqrtNoise->SetPoint( z, effS->GetBinCenter( i ), i_Signal_to_sqrtNoise );
+                    iGSignal_to_sqrtNoise->SetPoint( z, signalEff_mva, i_Signal_to_sqrtNoise );
                     if( fDebug )
                     {
-                        cout << "\t SET " << z << "\t" << effS->GetBinCenter( i ) << "\t" << i_Signal_to_sqrtNoise << endl;
+                        cout << "\t SET " << z << "\t" << signalEff_mva << "\t" << i_Signal_to_sqrtNoise << endl;
                     }
                     z++;
                 }
@@ -1400,40 +1405,37 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
                     cout << "___________________________________________________________" << endl;
                 }
             }
-        } // END loop over all signal efficiency bins
+        } // END loop over all signal efficiency bins for a given source strength
         /////////////////////////
         // determine position of maximum significance
         // fill a histogram from these values, smooth it, and determine position of maximum significance
         double i_xmax = -99.;
-        if( iGSignal_to_sqrtNoise )
+        TGraphSmooth* iGSmooth = new TGraphSmooth( "s" );
+        iGSignal_to_sqrtNoise_Smooth = iGSmooth->SmoothKern( iGSignal_to_sqrtNoise, "normal", 0.05, 100 );
+        // get maximum values
+        double x = 0.;
+        double y = 0.;
+        double i_ymax = -99.;
+        for( int i = 0; i < iGSignal_to_sqrtNoise_Smooth->GetN(); i++ )
         {
-            TGraphSmooth* iGSmooth = new TGraphSmooth( "s" );
-            iGSignal_to_sqrtNoise_Smooth = iGSmooth->SmoothKern( iGSignal_to_sqrtNoise, "normal", 0.05, 100 );
-            // get maximum values
-            double x = 0.;
-            double y = 0.;
-            double i_ymax = -99.;
-            for( int i = 0; i < iGSignal_to_sqrtNoise_Smooth->GetN(); i++ )
+            iGSignal_to_sqrtNoise_Smooth->GetPoint( i, x, y );
+            if( y > i_ymax )
             {
-                iGSignal_to_sqrtNoise_Smooth->GetPoint( i, x, y );
-                if( y > i_ymax )
-                {
-                    i_ymax = y;
-                    i_xmax = x;
-                }
-                // stop after first maximum (makes maximum research more robust to fluctuations of the
-                // background efficiency)
-                else if( y < i_ymax )
-                {
-                    break;
-                }
+                i_ymax = y;
+                i_xmax = x;
+            }
+            // stop after first maximum (makes maximum research more robust to fluctuations of the
+            // background efficiency)
+            else if( y < i_ymax )
+            {
+                break;
             }
-            i_SignalEfficiency_AtMaximum     = effS->GetBinContent( effS->FindBin( i_xmax ) );
-            i_BackgroundEfficiency_AtMaximum = effB->GetBinContent( effB->FindBin( i_xmax ) );
-            i_TMVACutValue_AtMaximum         = i_xmax;
-            i_Signal_to_sqrtNoise_atMaximum  = i_ymax;
-            i_SourceStrength_atMaximum       = iSourceStrength;
         }
+        i_SignalEfficiency_AtMaximum     = effS->GetBinContent( effS->FindBin( i_xmax ) );
+        i_BackgroundEfficiency_AtMaximum = effB->GetBinContent( effB->FindBin( i_xmax ) );
+        i_TMVACutValue_AtMaximum         = i_xmax;
+        i_Signal_to_sqrtNoise_atMaximum  = i_ymax;
+        i_SourceStrength_atMaximum       = iSourceStrength;
         ///////////////////////////////////////////////////////
         // check if value if really at the optimum or if information is missing from background efficiency curve
         // (check if maximum was find in the last bin or if next bin content is zero)
@@ -1469,7 +1471,7 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
         }
 
         // check detection criteria
-        if( i_Signal_to_sqrtNoise_atMaximum > fOptimizationSourceSignificance
+        if( i_Signal_to_sqrtNoise_atMaximum >= fOptimizationSourceSignificance
                 && Ndif < fOptimizationMinSignalEvents )
         {
             cout << "\t passed significance but not signal events criterium";
@@ -1477,8 +1479,8 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
             cout << "sig " << i_Signal_to_sqrtNoise_atMaximum;
             cout << ", Ndif " << Ndif << endl;
         }
-        if( i_Signal_to_sqrtNoise_atMaximum > fOptimizationSourceSignificance
-                && Ndif > fOptimizationMinSignalEvents )
+        if( i_Signal_to_sqrtNoise_atMaximum >= fOptimizationSourceSignificance
+                && Ndif >= fOptimizationMinSignalEvents )
         {
             break;
         }
@@ -1723,7 +1725,6 @@ void VTMVAEvaluator::plotEfficiencyPlotsPerBin( unsigned int iBin, TGraph* iGSig
             iLMinBack->Draw();
         }
     }
-
 }
 
 void VTMVAEvaluator::setTMVAMethod( string iMethodName )
@@ -1799,8 +1800,9 @@ void VTMVAEvaluator::printSensitivityOptimizationParameters()
 {
     cout << "VTMAEvaluator: MVA cut parameter is optimized for: " << endl;
     cout << "\t" << fOptimizationObservingTime_h << " hours of observing time" << endl;
-    cout << "\t" << fOptimizationSourceSignificance << " minimum significance" << endl;
+    cout << "\t" << fOptimizationSourceSignificance << " sigma minimum significance" << endl;
     cout << "\t" << fOptimizationMinSignalEvents << " minimum number of on events" << endl;
+    cout << "\t" << fOptimizationMinBackGroundEvents << " minimum number of off events" << endl;
     cout << "\t" << fOptimizationBackgroundAlpha << " signal to background area ratio" << endl;
     cout << "\t" << fOptimizationFixedSignalEfficiency << " maximum signal efficiency" << endl;
     cout << "\t" << fOptimizationMinSourceStrength << " minimum source strength" << endl;

From cab09dfa8324ddf6a995fc9f9dc8b740f66ecc25 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 16:36:09 +0100
Subject: [PATCH 09/16] improved plotting

---
 src/VTMVAEvaluator.cpp | 91 ++++++++++++++++++++++++++----------------
 1 file changed, 57 insertions(+), 34 deletions(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index c739ae8a..220a8587 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -863,14 +863,33 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         cout << "TMVAEvaluator::plotSignalAndBackgroundEfficiencies error: signal efficiency vector with size 0" << endl;
         return 0;
     }
+    unsigned int i_ze_min = 9999;
+    unsigned int i_ze_max = 0;
+    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
+    {
+        if( fTMVAData[i]->fZenithBin < i_ze_min )
+        {
+            i_ze_min = fTMVAData[i]->fZenithBin;
+        }
+        if( fTMVAData[i]->fZenithBin > i_ze_max )
+        {
+            i_ze_max = fTMVAData[i]->fZenithBin;
+        }
+    }
+    unsigned int i_n_ze_bins = i_ze_max - i_ze_min + 1;
 
-    // fill graphs
-    TGraphAsymmErrors* igSignal = new TGraphAsymmErrors( 1 );
-    TGraphAsymmErrors* igSignalOpt = new TGraphAsymmErrors( 1 );
-    TGraphAsymmErrors* igBck = new TGraphAsymmErrors( 1 );
-    TGraphAsymmErrors* igBckOpt = new TGraphAsymmErrors( 1 );
-    TGraphAsymmErrors* igCVa = new TGraphAsymmErrors( 1 );
-    TGraphAsymmErrors* igCVaOpt = new TGraphAsymmErrors( 1 );
+    vector<TGraphAsymmErrors*> igSignal_per_ze;
+    vector<TGraphAsymmErrors*> igBck_per_ze;
+    vector<TGraphAsymmErrors*> igCVa_per_ze;
+    for( unsigned int j = 0; j < i_n_ze_bins; j++ )
+    {
+        igSignal_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+        igSignalOpt_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+        igBck_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+        igBckOpt_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+        igCVa_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+        igCVaOpt_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
+    }
 
     unsigned int z_opt = 0;
     unsigned int z_noOpt = 0;
@@ -883,6 +902,12 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         {
             continue;
         }
+        TGraphAsymmErrors* igSignal = igSignal_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+        TGraphAsymmErrors* igSignalOpt = igSignalOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+        TGraphAsymmErrors* igBck = igBck_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+        TGraphAsymmErrors* igBckOpt = igBckOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+        TGraphAsymmErrors* igCVa = igCVa_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+        TGraphAsymmErrors* igCVaOpt = igCVaOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
 
         if( fTMVAData[i]->fSignalEfficiency < 0. || fTMVAData[i]->fBackgroundEfficiency < 0. )
         {
@@ -967,36 +992,30 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
     hnull->SetMaximum( 1. );
     plot_nullHistogram( iCanvas, hnull, false, false, 1.5, iE_min, iE_max );
 
-    setGraphPlottingStyle( igSignal, 1, 1., 20 );
-    setGraphPlottingStyle( igSignalOpt, 1, 1., 24 );
-    if( igBck )
+    for( unsigned int j = 0; j < i_n_ze_bins; j++ )
     {
-        setGraphPlottingStyle( igBck, 2, 1., 21 );
-    }
-    if( igBckOpt )
-    {
-        setGraphPlottingStyle( igBckOpt, 2, 1., 25 );
-    }
+        setGraphPlottingStyle( igSignal_per_ze[j], 1, 1., 20 + j );
+        setGraphPlottingStyle( igSignalOpt_per_ze[j], 1, 1., 24 + j );
+        setGraphPlottingStyle( igBck_per_ze[j], 2, 1., 21 + j );
+        setGraphPlottingStyle( igBckOpt_per_ze[j], 2, 1., 25 + j );
 
-    igSignal->Draw( "pl" );
-    if( z_noOpt )
-    {
-        igSignalOpt->Draw( "pl" );
-    }
-    if( igBck )
-    {
-        igBck->Draw( "pl" );
-    }
-    if( igBckOpt && z_noOpt > 0 )
-    {
-        igBckOpt->Draw( "pl" );
+        igSignal_per_ze[j]->Draw( "pl" );
+        igBck_per_ze[j]->Draw( "pl" );
+        if( z_noOpt )
+        {
+            igSignalOpt_per_ze[j]->Draw( "pl" );
+        }
+        if( igBckOpt_per_ze[j] && z_noOpt > 0 )
+        {
+            igBckOpt_per_ze[j]->Draw( "pl" );
+        }
     }
     if( fPrintPlotting )
     {
         iCanvas->Print( "MVA-SignalBackgroundEfficiency.pdf" );
     }
 
-    // plot MVA cut value
+    // plot MVA cut values
     if( igCVa )
     {
         TCanvas* iCVACanvas = new TCanvas( "iCVACanvas", "MVA cut value", 500, 10, 400, 400 );
@@ -1010,12 +1029,16 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         hnull->SetMinimum( iMVA_min );
         hnull->SetMaximum( iMVA_max );
         plot_nullHistogram( iCanvas, hnull, false, false, 1.3, iE_min, iE_max );
-        setGraphPlottingStyle( igCVa, 1, 1., 20 );
-        igCVa->Draw( "p" );
-        if( igCVaOpt && z_noOpt > 0 )
+        for( unsigned int j = 0; j < i_n_ze_bins; j++ )
         {
-            setGraphPlottingStyle( igCVaOpt, 1, 1., 24 );
-            igCVaOpt->Draw( "p" );
+            setGraphPlottingStyle( igCVa_per_ze[j], 1, 1., 20 + j );
+            setGraphPlottingStyle( igCVaOpt_per_ze[j], 1, 1., 24 + j );
+
+            igCVa_per_ze[j]->Draw( "pl" );
+            if( z_noOpt )
+            {
+                igCVaOpt_per_ze[j]->Draw( "pl" );
+            }
         }
         if( fPrintPlotting )
         {

From d52bf5a2e9f1a68b5bd12fb87b947fd0f35c2688 Mon Sep 17 00:00:00 2001
From: GernotMaier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 17:19:31 +0100
Subject: [PATCH 10/16] improved plotting

---
 inc/VTMVAEvaluator.h     |   2 +-
 macros/optimizeBDTcuts.C |   2 +-
 src/VTMVAEvaluator.cpp   | 208 ++++++++++++++++++++-------------------
 3 files changed, 107 insertions(+), 105 deletions(-)

diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h
index 3099b4a4..1fdde49f 100644
--- a/inc/VTMVAEvaluator.h
+++ b/inc/VTMVAEvaluator.h
@@ -193,7 +193,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities
         {
             return fIsZombie;
         }
-        TGraphAsymmErrors* plotSignalAndBackgroundEfficiencies( bool iLogY = true, double iYmin = 1.e-4, double iMVA_min = -1., double iMVA_max = 1. );
+        void plotSignalAndBackgroundEfficiencies( bool iLogY = true, double iYmin = 1.e-4, double iMVA_min = -1., double iMVA_max = 1. );
         void   printOptimizedMVACutValues( string iEpoch = "V6" );
         void   printSensitivityOptimizationParameters();
         void   printSignalEfficiency();
diff --git a/macros/optimizeBDTcuts.C b/macros/optimizeBDTcuts.C
index fc1ad1d5..2b65c62e 100644
--- a/macros/optimizeBDTcuts.C
+++ b/macros/optimizeBDTcuts.C
@@ -57,7 +57,7 @@ void optimizeBDTcuts(
         weightFileIndex_Zmin, weightFileIndex_Zmax );
 
     // plotting
-    a.plotSignalAndBackgroundEfficiencies();
+    a.plotSignalAndBackgroundEfficiencies(true, 7.e-3, -0.2, 0.7);
 
     // print results to screen
     a.printSignalEfficiency();
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 220a8587..590e1734 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -855,13 +855,12 @@ bool VTMVAEvaluator::initializeDataStructures( CData* iC )
  * plot signal and background efficiencies
  *
  */
-TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
-    bool iLogY, double iYmin, double iMVA_min, double iMVA_max )
+void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmin, double iMVA_min, double iMVA_max )
 {
     if( fTMVAData.size() == 0 )
     {
         cout << "TMVAEvaluator::plotSignalAndBackgroundEfficiencies error: signal efficiency vector with size 0" << endl;
-        return 0;
+        return;
     }
     unsigned int i_ze_min = 9999;
     unsigned int i_ze_max = 0;
@@ -879,8 +878,11 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
     unsigned int i_n_ze_bins = i_ze_max - i_ze_min + 1;
 
     vector<TGraphAsymmErrors*> igSignal_per_ze;
+    vector<TGraphAsymmErrors*> igSignalOpt_per_ze;
     vector<TGraphAsymmErrors*> igBck_per_ze;
+    vector<TGraphAsymmErrors*> igBckOpt_per_ze;
     vector<TGraphAsymmErrors*> igCVa_per_ze;
+    vector<TGraphAsymmErrors*> igCVaOpt_per_ze;
     for( unsigned int j = 0; j < i_n_ze_bins; j++ )
     {
         igSignal_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
@@ -891,10 +893,12 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         igCVaOpt_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
     }
 
-    unsigned int z_opt = 0;
-    unsigned int z_noOpt = 0;
+    vector< unsigned int > z_opt(i_n_ze_bins, 0 );
+    vector< unsigned int > z_noOpt(i_n_ze_bins, 0 );
 
     double iMinBck = 1.;
+    double iE_min =  1.e99;
+    double iE_max = -1.e99;
 
     for( unsigned int i = 0; i < fTMVAData.size(); i++ )
     {
@@ -902,12 +906,24 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         {
             continue;
         }
-        TGraphAsymmErrors* igSignal = igSignal_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
-        TGraphAsymmErrors* igSignalOpt = igSignalOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
-        TGraphAsymmErrors* igBck = igBck_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
-        TGraphAsymmErrors* igBckOpt = igBckOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
-        TGraphAsymmErrors* igCVa = igCVa_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
-        TGraphAsymmErrors* igCVaOpt = igCVaOpt_per_ze[fTMVAData[i]->fZenithBin - i_ze_min];
+
+        // max energy range
+        if( fTMVAData[i] && fTMVAData[i]->fEnergyCut_Log10TeV_min < iE_min )
+        {
+            iE_min = fTMVAData[i]->fEnergyCut_Log10TeV_min;
+        }
+        if( fTMVAData[i] && fTMVAData[i]->fEnergyCut_Log10TeV_max > iE_max )
+        {
+            iE_max = fTMVAData[i]->fEnergyCut_Log10TeV_max;
+        }
+        unsigned int ze_bin = fTMVAData[i]->fZenithBin - i_ze_min;
+
+        TGraphAsymmErrors* igSignal = igSignal_per_ze[ze_bin];
+        TGraphAsymmErrors* igSignalOpt = igSignalOpt_per_ze[ze_bin];
+        TGraphAsymmErrors* igBck = igBck_per_ze[ze_bin];
+        TGraphAsymmErrors* igBckOpt = igBckOpt_per_ze[ze_bin];
+        TGraphAsymmErrors* igCVa = igCVa_per_ze[ze_bin];
+        TGraphAsymmErrors* igCVaOpt = igCVaOpt_per_ze[ze_bin];
 
         if( fTMVAData[i]->fSignalEfficiency < 0. || fTMVAData[i]->fBackgroundEfficiency < 0. )
         {
@@ -918,35 +934,35 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
         }
         if( fTMVAData[i]->fTMVAOptimumCutValueFound )
         {
-            igSignal->SetPoint( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fSignalEfficiency );
-            igSignal->SetPointEXlow( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igSignal->SetPointEXhigh( z_opt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igSignal->SetPoint( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fSignalEfficiency );
+            igSignal->SetPointEXlow( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igSignal->SetPointEXhigh( z_opt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            igBck->SetPoint( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fBackgroundEfficiency );
-            igBck->SetPointEXlow( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igBck->SetPointEXhigh( z_opt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igBck->SetPoint( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fBackgroundEfficiency );
+            igBck->SetPointEXlow( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igBck->SetPointEXhigh( z_opt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            igCVa->SetPoint( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fTMVACutValue );
-            igCVa->SetPointEXlow( z_opt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igCVa->SetPointEXhigh( z_opt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igCVa->SetPoint( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fTMVACutValue );
+            igCVa->SetPointEXlow( z_opt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igCVa->SetPointEXhigh( z_opt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            z_opt++;
+            z_opt[ze_bin]++;
         }
         else if( fTMVAData[i]->fTMVACutValue > -90. )
         {
-            igSignalOpt->SetPoint( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fSignalEfficiency );
-            igSignalOpt->SetPointEXlow( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igSignalOpt->SetPointEXhigh( z_noOpt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igSignalOpt->SetPoint( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fSignalEfficiency );
+            igSignalOpt->SetPointEXlow( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igSignalOpt->SetPointEXhigh( z_noOpt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            igBckOpt->SetPoint( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fBackgroundEfficiency );
-            igBckOpt->SetPointEXlow( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igBckOpt->SetPointEXhigh( z_noOpt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igBckOpt->SetPoint( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fBackgroundEfficiency );
+            igBckOpt->SetPointEXlow( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igBckOpt->SetPointEXhigh( z_noOpt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            igCVaOpt->SetPoint( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fTMVACutValue );
-            igCVaOpt->SetPointEXlow( z_noOpt, fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
-            igCVaOpt->SetPointEXhigh( z_noOpt, fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
+            igCVaOpt->SetPoint( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV, fTMVAData[i]->fTMVACutValue );
+            igCVaOpt->SetPointEXlow( z_noOpt[ze_bin], fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV - fTMVAData[i]->fEnergyCut_Log10TeV_min );
+            igCVaOpt->SetPointEXhigh( z_noOpt[ze_bin], fTMVAData[i]->fEnergyCut_Log10TeV_max - fTMVAData[i]->fSpectralWeightedMeanEnergy_Log10TeV );
 
-            z_noOpt++;
+            z_noOpt[ze_bin]++;
         }
         if( fTMVAData[i]->fBackgroundEfficiency < iMinBck )
         {
@@ -956,97 +972,83 @@ TGraphAsymmErrors* VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(
 
     // plot everything
     TCanvas* iCanvas = new TCanvas( "cSignalAndBackgroundEfficiencies", "signal and background efficiencies",
-                                    10, 10, 400, 400 );
-    iCanvas->SetGridx( 0 );
-    iCanvas->SetGridy( 0 );
-    iCanvas->SetLeftMargin( 0.13 );
-    if( iLogY )
-    {
-        iCanvas->SetLogy();
-    }
-    else
-    {
-        iCanvas->SetLogy( 0 );
-    }
+                                    10, 10, 1200, 800 );
+    iCanvas->Divide(i_n_ze_bins, 3);
     iCanvas->Draw();
-
-    double iE_min =  1.e99;
-    double iE_max = -1.e99;
-    for( unsigned int i = 0; i < fTMVAData.size(); i++ )
-    {
-        if( fTMVAData[i] && fTMVAData[i]->fEnergyCut_Log10TeV_min < iE_min )
-        {
-            iE_min = fTMVAData[i]->fEnergyCut_Log10TeV_min;
-        }
-        if( fTMVAData[i] && fTMVAData[i]->fEnergyCut_Log10TeV_max > iE_max )
-        {
-            iE_max = fTMVAData[i]->fEnergyCut_Log10TeV_max;
-        }
-    }
-
-    TH1D* hnull = new TH1D( "hnullcSignalAndBackgroundEfficiencies", "", 100, iE_min, iE_max );
-    hnull->SetStats( 0 );
-    hnull->SetXTitle( "energy [TeV]" );
-    hnull->SetYTitle( "signal/background efficiency" );
-    hnull->SetMinimum( iYmin );
-    hnull->SetMaximum( 1. );
-    plot_nullHistogram( iCanvas, hnull, false, false, 1.5, iE_min, iE_max );
-
+    char hname[200];
+    char htitle[200];
     for( unsigned int j = 0; j < i_n_ze_bins; j++ )
     {
+        TPad *iPad = (TPad*)iCanvas->cd(j+1);
+        iPad->SetLeftMargin( 0.13 );
+        sprintf( hname, "hnullcSignalEfficiencies_%d", j );
+        sprintf( htitle, "signal efficiency (ze %d)", j );
+        TH1D* hnull = new TH1D( hname, htitle, 100, iE_min, iE_max );
+        hnull->SetStats( 0 );
+        hnull->SetXTitle( "energy [TeV]" );
+        hnull->SetYTitle( "signal efficiency" );
+        hnull->SetMinimum( 0. );
+        hnull->SetMaximum( 1. );
+        plot_nullHistogram( iPad, hnull, false, false, 1.5, iE_min, iE_max );
         setGraphPlottingStyle( igSignal_per_ze[j], 1, 1., 20 + j );
         setGraphPlottingStyle( igSignalOpt_per_ze[j], 1, 1., 24 + j );
-        setGraphPlottingStyle( igBck_per_ze[j], 2, 1., 21 + j );
-        setGraphPlottingStyle( igBckOpt_per_ze[j], 2, 1., 25 + j );
 
         igSignal_per_ze[j]->Draw( "pl" );
-        igBck_per_ze[j]->Draw( "pl" );
-        if( z_noOpt )
+        if( z_noOpt[j] )
         {
             igSignalOpt_per_ze[j]->Draw( "pl" );
         }
-        if( igBckOpt_per_ze[j] && z_noOpt > 0 )
+
+        // background efficiency
+        iPad = (TPad*)iCanvas->cd( j+1+i_n_ze_bins );
+        iPad->SetLeftMargin( 0.13 );
+        if( iLogY ) iPad->SetLogy();
+        else iPad->SetLogy( 0 );
+        sprintf( hname, "hnullcBackgroundEfficiencies_%d", j );
+        sprintf( htitle, "background efficiency (ze %d)", j );
+        TH1D* hnull_bck = new TH1D( hname, htitle, 100, iE_min, iE_max );
+        hnull_bck->SetStats( 0 );
+        hnull_bck->SetXTitle( "energy [TeV]" );
+        hnull_bck->SetYTitle( "background efficiency" );
+        hnull_bck->SetMinimum( iYmin );
+        hnull_bck->SetMaximum( 1. );
+        plot_nullHistogram( iPad, hnull_bck, false, false, 1.5, iE_min, iE_max );
+        setGraphPlottingStyle( igBck_per_ze[j], 2, 1., 21 + j );
+        setGraphPlottingStyle( igBckOpt_per_ze[j], 2, 1., 25 + j );
+
+        igBck_per_ze[j]->Draw( "pl" );
+        if( igBckOpt_per_ze[j] && z_noOpt[j] > 0 )
         {
             igBckOpt_per_ze[j]->Draw( "pl" );
         }
+
+        // MVA cut
+        iPad = (TPad*)iCanvas->cd( j+1+i_n_ze_bins*2 );
+
+        sprintf(hname, "hnull_mvacMVACuts_%d", j );
+        sprintf( htitle, "MVA cut variable (ze %d)", j );
+        TH1D* hnull_mva = new TH1D( hname, htitle, 100, iE_min, iE_max );
+        hnull_mva->SetStats( 0 );
+        hnull_mva->SetXTitle( "energy [TeV]" );
+        hnull_mva->SetYTitle( "MVA cut variable" );
+        hnull_mva->SetMinimum( iMVA_min );
+        hnull_mva->SetMaximum( iMVA_max );
+        plot_nullHistogram( iPad, hnull_mva, false, false, 1.3, iE_min, iE_max );
+        setGraphPlottingStyle( igCVa_per_ze[j], 1, 1., 20 + j );
+        setGraphPlottingStyle( igCVaOpt_per_ze[j], 1, 1., 24 + j );
+
+        igCVa_per_ze[j]->Draw( "pl" );
+        if( z_noOpt[j] )
+        {
+            igCVaOpt_per_ze[j]->Draw( "pl" );
+        }
     }
+
     if( fPrintPlotting )
     {
         iCanvas->Print( "MVA-SignalBackgroundEfficiency.pdf" );
     }
 
-    // plot MVA cut values
-    if( igCVa )
-    {
-        TCanvas* iCVACanvas = new TCanvas( "iCVACanvas", "MVA cut value", 500, 10, 400, 400 );
-        iCVACanvas->SetGridx( 0 );
-        iCVACanvas->SetGridy( 0 );
-
-        TH1D* hnull = new TH1D( "hnullcMVACuts", "", 100, iE_min, iE_max );
-        hnull->SetStats( 0 );
-        hnull->SetXTitle( "energy [TeV]" );
-        hnull->SetYTitle( "MVA cut variable" );
-        hnull->SetMinimum( iMVA_min );
-        hnull->SetMaximum( iMVA_max );
-        plot_nullHistogram( iCanvas, hnull, false, false, 1.3, iE_min, iE_max );
-        for( unsigned int j = 0; j < i_n_ze_bins; j++ )
-        {
-            setGraphPlottingStyle( igCVa_per_ze[j], 1, 1., 20 + j );
-            setGraphPlottingStyle( igCVaOpt_per_ze[j], 1, 1., 24 + j );
-
-            igCVa_per_ze[j]->Draw( "pl" );
-            if( z_noOpt )
-            {
-                igCVaOpt_per_ze[j]->Draw( "pl" );
-            }
-        }
-        if( fPrintPlotting )
-        {
-            iCVACanvas->Print( "MVA-MVACut.pdf" );
-        }
-    }
-
-    return igCVa;
 }
 
 void VTMVAEvaluator::setSignalEfficiency( double iSignalEfficiency )

From 15c53c677c059135fd17339f1cf97b808afd33ba Mon Sep 17 00:00:00 2001
From: GernotMaier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 17:35:09 +0100
Subject: [PATCH 11/16] ze range

---
 src/VTMVAEvaluator.cpp | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 590e1734..ba408f43 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -1903,11 +1903,20 @@ TGraph* VTMVAEvaluator::fillfromGraph2D( TObject* i_G, double i_ze_min, double i
         TGraph2D* iG2D = ( TGraph2D* )i_G;
         TGraph* iG1D = new TGraph( iG2D->GetN() );
         Double_t* x = iG2D->GetX();
+        Double_t* y = iG2D->GetY();
+
+        // reset ze max to max of Graph
+        double ze_max_graph = TMath::MaxElement(iG2D->GetN(), y);
+        if( ze_max_graph < i_ze_max ) i_ze_max = ze_max_graph;
+
         double z1 = i_ze_min * TMath::DegToRad();
         double z2 = i_ze_max * TMath::DegToRad();
         double avg_airmass = 0.5 * ( (1.0 / TMath::Cos(z1)) + (1.0 / TMath::Cos(z2)) );
         double ze_mean = TMath::ACos(1.0 / avg_airmass) * TMath::RadToDeg();
 
+        cout << "Graph filling from " << i_G->GetName() << ": average (airmass) ze ";
+        cout << ze_mean << " (" << i_ze_min << ", " << i_ze_max << ") deg " << endl;
+
         for( int i = 0; i < iG2D->GetN(); i++ )
         {
             Double_t y = iG2D->Interpolate( x[i], ze_mean );

From d7c5fd96676d59f4a3cf3a3898cdb3222a69738e Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Thu, 25 Dec 2025 17:54:55 +0100
Subject: [PATCH 12/16] comment

---
 src/VTMVAEvaluator.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index ba408f43..51c44269 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -744,7 +744,7 @@ bool VTMVAEvaluator::evaluate( bool interpolate_mva, bool use_average_zenith_ang
         fTMVA_EvaluationResult = fTMVAData[iDataBin]->fTMVAReader->EvaluateMVA( fTMVAData[iDataBin]->fTMVAMethodTag_2 );
     }
 
-    // apply MVA cut (cut value not interpolated in zenith)
+    // apply MVA cut
     if( fTMVA_EvaluationResult < fTMVAData[iDataBin]->fTMVACutValue )
     {
         return false;

From 396c5c0584aab47c6328a8810aca59275b081146 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Sat, 3 Jan 2026 16:01:08 +0100
Subject: [PATCH 13/16] typos

---
 src/VTMVAEvaluator.cpp          | 2 +-
 src/calculateCrabRateFromMC.cpp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index 51c44269..a97fce46 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -1306,7 +1306,7 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
     //////////////////////////////////////////////////////
     // loop over different source strengths (in Crab Units)
     // start from low values to high values; optimal cut value
-    // is found when we reach fir time the required significance
+    // is found when we reach first time the required significance
     // source strength steps on log scale (up to 30 CU)
     unsigned int iSourceStrengthStepSizeN =
         ( unsigned int )(( log10( 30. ) - log10( fOptimizationMinSourceStrength ) ) / 0.005 );
diff --git a/src/calculateCrabRateFromMC.cpp b/src/calculateCrabRateFromMC.cpp
index 4d20d9ac..cc04350e 100644
--- a/src/calculateCrabRateFromMC.cpp
+++ b/src/calculateCrabRateFromMC.cpp
@@ -476,7 +476,7 @@ TGraph2DErrors* getTGraph2D( vector< TProfile2D* > h, string iGraphName )
             {
                 i_z += h[0]->GetBinContent( b_x, b_y );
             }
-            // convert expect rate graphs in 1./s
+            // convert rate graphs to 1./s
             i_g->SetPoint(
                 z,
                 h[1]->GetXaxis()->GetBinCenter( b_x ),

From bd6c57e07835751f8826a922d5fa62dbdb52814a Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Sat, 3 Jan 2026 16:02:00 +0100
Subject: [PATCH 14/16] changelog

---
 docs/changes/336.maintenance.md | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 docs/changes/336.maintenance.md

diff --git a/docs/changes/336.maintenance.md b/docs/changes/336.maintenance.md
new file mode 100644
index 00000000..2f351cd5
--- /dev/null
+++ b/docs/changes/336.maintenance.md
@@ -0,0 +1,2 @@
+Removed deprecated smoothing and interpolation functionality for MVA cut values.
+Cleaned up unused parameters and obsolete code related to box smoothing and TMVA optimization.

From bad13b935ccad39aa20661e327830143bd0f2dea Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Sat, 3 Jan 2026 16:07:47 +0100
Subject: [PATCH 15/16] date

---
 LICENSE | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/LICENSE b/LICENSE
index 4d2a0cd4..0aff9286 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2022-2025 Eventdisplay Developers
+Copyright (c) 2022-2026 Eventdisplay Developers
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

From 8a668adc43c197b483b88f1142045e3d01559d76 Mon Sep 17 00:00:00 2001
From: GernotMaier <gernot.maier@desy.de>
Date: Sat, 3 Jan 2026 16:19:36 +0100
Subject: [PATCH 16/16] astyle

---
 src/VStereoMaps.cpp             |  2 +-
 src/VTMVAEvaluator.cpp          | 30 +++++++++++++++---------------
 src/calculateCrabRateFromMC.cpp |  2 +-
 3 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/VStereoMaps.cpp b/src/VStereoMaps.cpp
index d08c6c0c..7a3e379c 100644
--- a/src/VStereoMaps.cpp
+++ b/src/VStereoMaps.cpp
@@ -399,7 +399,7 @@ void VStereoMaps::makeTwoDStereo_BoxSmooth( double i_xderot, double i_yderot, do
                     hmap_ratio->Fill( i_MeanSignalBackgroundAreaRatio );
                 }
             }
-    }
+        }
     }
 }
 
diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp
index a97fce46..adeaa657 100644
--- a/src/VTMVAEvaluator.cpp
+++ b/src/VTMVAEvaluator.cpp
@@ -855,7 +855,7 @@ bool VTMVAEvaluator::initializeDataStructures( CData* iC )
  * plot signal and background efficiencies
  *
  */
-void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmin, double iMVA_min, double iMVA_max )
+void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies( bool iLogY, double iYmin, double iMVA_min, double iMVA_max )
 {
     if( fTMVAData.size() == 0 )
     {
@@ -893,8 +893,8 @@ void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmi
         igCVaOpt_per_ze.push_back( new TGraphAsymmErrors( 1 ) );
     }
 
-    vector< unsigned int > z_opt(i_n_ze_bins, 0 );
-    vector< unsigned int > z_noOpt(i_n_ze_bins, 0 );
+    vector< unsigned int > z_opt( i_n_ze_bins, 0 );
+    vector< unsigned int > z_noOpt( i_n_ze_bins, 0 );
 
     double iMinBck = 1.;
     double iE_min =  1.e99;
@@ -973,13 +973,13 @@ void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmi
     // plot everything
     TCanvas* iCanvas = new TCanvas( "cSignalAndBackgroundEfficiencies", "signal and background efficiencies",
                                     10, 10, 1200, 800 );
-    iCanvas->Divide(i_n_ze_bins, 3);
+    iCanvas->Divide( i_n_ze_bins, 3 );
     iCanvas->Draw();
     char hname[200];
     char htitle[200];
     for( unsigned int j = 0; j < i_n_ze_bins; j++ )
     {
-        TPad *iPad = (TPad*)iCanvas->cd(j+1);
+        TPad *iPad = ( TPad* )iCanvas->cd( j + 1 );
         iPad->SetLeftMargin( 0.13 );
         sprintf( hname, "hnullcSignalEfficiencies_%d", j );
         sprintf( htitle, "signal efficiency (ze %d)", j );
@@ -1000,7 +1000,7 @@ void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmi
         }
 
         // background efficiency
-        iPad = (TPad*)iCanvas->cd( j+1+i_n_ze_bins );
+        iPad = ( TPad* )iCanvas->cd( j + 1 + i_n_ze_bins );
         iPad->SetLeftMargin( 0.13 );
         if( iLogY ) iPad->SetLogy();
         else iPad->SetLogy( 0 );
@@ -1023,9 +1023,9 @@ void VTMVAEvaluator::plotSignalAndBackgroundEfficiencies(bool iLogY, double iYmi
         }
 
         // MVA cut
-        iPad = (TPad*)iCanvas->cd( j+1+i_n_ze_bins*2 );
+        iPad = ( TPad* )iCanvas->cd( j + 1 + i_n_ze_bins * 2 );
 
-        sprintf(hname, "hnull_mvacMVACuts_%d", j );
+        sprintf( hname, "hnull_mvacMVACuts_%d", j );
         sprintf( htitle, "MVA cut variable (ze %d)", j );
         TH1D* hnull_mva = new TH1D( hname, htitle, 100, iE_min, iE_max );
         hnull_mva->SetStats( 0 );
@@ -1379,8 +1379,8 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
                 {
                     // optimize signal/sqrt(noise) - Li & Ma significance
                     i_Signal_to_sqrtNoise = VStatistics::calcSignificance(
-                                                signalEff * Ndif + backEff * Nof,
-                                                backEff * Nof / fOptimizationBackgroundAlpha,
+                                                signalEff* Ndif + backEff* Nof,
+                                                backEff* Nof / fOptimizationBackgroundAlpha,
                                                 fOptimizationBackgroundAlpha );
                 }
                 else
@@ -1400,8 +1400,8 @@ bool VTMVAEvaluator::optimizeSensitivity( unsigned int iDataBin )
                 }
                 if( signalEff * Ndif > 0. )
                 {
-                    iGSignalEvents->SetPoint( z_SB, signalEff_mva,  signalEff * Ndif );
-                    iGBackgroundEvents->SetPoint( z_SB, signalEff_mva, backEff * Nof );
+                    iGSignalEvents->SetPoint( z_SB, signalEff_mva,  signalEff* Ndif );
+                    iGBackgroundEvents->SetPoint( z_SB, signalEff_mva, backEff* Nof );
                     z_SB++;
                 }
                 // check that a minimum number of off events is available
@@ -1906,13 +1906,13 @@ TGraph* VTMVAEvaluator::fillfromGraph2D( TObject* i_G, double i_ze_min, double i
         Double_t* y = iG2D->GetY();
 
         // reset ze max to max of Graph
-        double ze_max_graph = TMath::MaxElement(iG2D->GetN(), y);
+        double ze_max_graph = TMath::MaxElement( iG2D->GetN(), y );
         if( ze_max_graph < i_ze_max ) i_ze_max = ze_max_graph;
 
         double z1 = i_ze_min * TMath::DegToRad();
         double z2 = i_ze_max * TMath::DegToRad();
-        double avg_airmass = 0.5 * ( (1.0 / TMath::Cos(z1)) + (1.0 / TMath::Cos(z2)) );
-        double ze_mean = TMath::ACos(1.0 / avg_airmass) * TMath::RadToDeg();
+        double avg_airmass = 0.5 * (( 1.0 / TMath::Cos( z1 ) ) + ( 1.0 / TMath::Cos( z2 ) ) );
+        double ze_mean = TMath::ACos( 1.0 / avg_airmass ) * TMath::RadToDeg();
 
         cout << "Graph filling from " << i_G->GetName() << ": average (airmass) ze ";
         cout << ze_mean << " (" << i_ze_min << ", " << i_ze_max << ") deg " << endl;
diff --git a/src/calculateCrabRateFromMC.cpp b/src/calculateCrabRateFromMC.cpp
index cc04350e..a02be0a8 100644
--- a/src/calculateCrabRateFromMC.cpp
+++ b/src/calculateCrabRateFromMC.cpp
@@ -131,7 +131,7 @@ vector< pair< double, double > > read_energy_minmax_pairs(
  * ENERGYBINEDGES -1.5 -0.5 -1. 0. -0.5 0.5 0.0 1. 0.5 2.0
  * ENERGYBINS -1.50 -0.75 -0.25 0.5 2.0
  */
-vector< double > read_energy_bins(string iTMVAParameterFile, string iEnergyKeyWord )
+vector< double > read_energy_bins( string iTMVAParameterFile, string iEnergyKeyWord )
 {
     vector< double > tmp_e;
     ifstream is( iTMVAParameterFile.c_str(), ifstream::in );