v492-dev17

CHANGELOG.md

@@ -3,7 +3,7 @@
 All notable changes to the Eventdisplay_v4 project will be documented in this file.
 Changes for upcoming releases can be found in the [docs/changes](docs/changes) directory.
 Note that changes before release v492.0 are not documented here, but can be found in the
-[GitHub repository](https://github.com/VERITAS-Observatory/VTS-SimPipe/releases).
+[GitHub repository](https://github.com/VERITAS-Observatory/EventDisplay_v4/releases).
 
 This changelog is generated using [Towncrier](https://towncrier.readthedocs.io/).
 

README.md

@@ -8,12 +8,12 @@
 ## Overview
 
 Eventdisplay is a reconstruction and analysis pipeline for data of
-Imaging Atmospheric Cherenkov Telescopes.
+Imaging Atmospheric Cherenkov Telescopes (IACT).
 It has been primarily developed for the analysis of VERITAS and CTAO data to generate results for scientific publications.
 This repository contains the Eventdisplay version used for VERITAS analysis ([CTAO version](https://github.com/Eventdisplay/Eventdisplay))
 
 * Original developers: Gernot Maier and Jamie Holder
-* Authors and contributors: [CITATION.cff](CITATION.cff)
+* Authors and contributors: [CITATION.cff](CITATION.cff) and [CONTRIBUTORS.md](CONTRIBUTORS.md)
 * License: [LICENSE](LICENSE)
 * Contributing: [CONTRIBUTING.md](CONTRIBUTING.md)
 * Changelog: [CHANGELOG.md](CHANGELOG.md)

docs/changes/351.bugfix.md

@@ -0,0 +1,2 @@
+Fixes long list of bugs listed in [Issue #352](https://github.com/VERITAS-Observatory/EventDisplay_v4/issues/352).
+Fix A2: `VTraceHandler::calculateTraceSum_slidingWindow()` fRaw branch returned `FADC[1]` (second FADC sample, ped-subtracted) instead of the accumulated raw charge sum. Affected pedestal calculation path.

docs/changes/352.bugfix.md

@@ -0,0 +1,5 @@
+# Bug fix: non-XGB jobs no longer probe XGB sidecar files (E4)
+
+XGB sidecar files (stereo and gamma-hadron) are now opened only when the active
+reconstruction method and gamma-hadron cut type actually require them, eliminating
+spurious ROOT file-open errors in non-XGB jobs.

docs/changes/353.bugfix.md

@@ -0,0 +1,7 @@
+# Bug fix: explicit NegWeightTreatment for TMVA Gradient BDT angular reconstruction
+
+Added `NegWeightTreatment=Pray` to the default TMVA options in
+`trainTMVAforAngularReconstruction` to match actual TMVA behaviour:
+`BoostType=Grad` does not support `InverseBoostNegWeights` (TMVA's global
+default) and silently replaces it with `Pray`; the options string now
+declares this explicitly so the training configuration matches what is used.

inc/CData.h

@@ -30,6 +30,7 @@ class CData
         vector< double > fTelX;
         vector< double > fTelY;
         vector< double > fTelZ;
+        string fDataFileName;
 
         void reconstruct_3tel_images( long unsigned int );
         void reconstruct_3tel_images_direction();
@@ -280,6 +281,7 @@ class CData
         float get_Yoff( unsigned int method );
         void   initialize_xgb_tree();
         TTree* getXGBTree( string file_name, string suffix, string tree_name );
+        bool   loadGHXGBTree( string gh_suffix );
         pair<float, float> get_XYoff_derot( unsigned int method );
         virtual Long64_t LoadTree( Long64_t entry );
         float            get_GH_Gamma_Prediction();

src/CData.cpp

@@ -31,6 +31,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, string file_name, string stere
     fShort = bShort;
     fVersion = 6;
     fTelescopeCombination = 0;
+    fDataFileName = file_name;
     Init( tree );
 
     fStereoFriendTree = getXGBTree( file_name, stereo_suffix, "StereoAnalysis" );
@@ -1077,7 +1078,15 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name
         return 0;
     }
 
-    file_name = file_name.replace( file_name.find( ".root" ), 5, "." + file_suffix + ".root" );
+    size_t iRootPos = file_name.rfind( ".root" );
+    if( iRootPos != string::npos && iRootPos + 5 == file_name.size() )
+    {
+        file_name.replace( iRootPos, 5, "." + file_suffix + ".root" );
+    }
+    else
+    {
+        file_name += "." + file_suffix + ".root";
+    }
     TFile *iFile = TFile::Open( file_name.c_str() );
     if(!iFile || iFile->IsZombie() )
     {
@@ -1126,3 +1135,25 @@ void CData::initialize_xgb_tree()
         GH_Is_Gamma = false;
     }
 }
+
+/*
+ * Deferred loading of the XGB gamma-hadron friend tree.
+ * Call this after the analysis type is known (e.g. after readCuts()),
+ * when the caller has confirmed that XGB gamma-hadron cuts are active.
+ * Idempotent: a second call is a no-op if the tree is already loaded.
+ */
+bool CData::loadGHXGBTree( string gh_suffix )
+{
+    if( fGHFriendTree )
+    {
+        return true;
+    }
+    if( fDataFileName.empty() )
+    {
+        cout << "CData::loadGHXGBTree error: no data file name stored" << endl;
+        return false;
+    }
+    fGHFriendTree = getXGBTree( fDataFileName, gh_suffix, "Classification" );
+    initialize_xgb_tree();
+    return ( fGHFriendTree != 0 );
+}

src/VArrayAnalyzer.cpp

@@ -944,7 +944,7 @@ int VArrayAnalyzer::rcs_method_0( unsigned int iMethod )
 bool VArrayAnalyzer::fillShowerCore( unsigned int iMeth, float ximp, float yimp )
 {
     // check validity
-    if(!isnormal( ximp ) || !isnormal( ximp ) )
+    if(!isnormal( ximp ) || !isnormal( yimp ) )
     {
         ximp = -99999.;
         yimp = -99999.;
@@ -1243,7 +1243,7 @@ int VArrayAnalyzer::rcs_method_3( unsigned int iMethod )
             iangdiff = sin( fabs( atan( m[jj] ) - atan( m[ii] ) ) );
 
             // discard all pairs with almost parallel lines
-            float i_diff =  fabs( atan( m[0] ) - atan( m[1] ) );
+            float i_diff =  fabs( atan( mm[0] ) - atan( mm[1] ) );
             if( i_diff < fEvndispReconstructionParameter->fAxesAngles_min[iMethod] / TMath::RadToDeg() ||
                     fabs( 180. / TMath::RadToDeg() - i_diff ) < fEvndispReconstructionParameter->fAxesAngles_min[iMethod] / TMath::RadToDeg() )
             {

src/VCalibrationData.cpp

@@ -367,6 +367,7 @@ void  VCalibrationData::setAverageTZero( double iAverageTzero, bool iLowGain )
     if( iLowGain )
     {
         fAverageTZero_lowgain = iAverageTzero;
+        return;
     }
 
     fAverageTZero_highgain = iAverageTzero;

src/VEmissionHeightCalculator.cpp

@@ -53,12 +53,18 @@ double VEmissionHeightCalculator::getEmissionHeight( float* cen_x, float* cen_y,
                     fImageDistance = TMath::Tan( imageDistance( cen_x[i], cen_x[j], cen_y[i], cen_y[j] ) / TMath::RadToDeg() );
                     if( fImageDistance > 0. )
                     {
+                        double iLogSizeI = log10( size[i] );
+                        double iLogSizeJ = log10( size[j] );
+                        if( iLogSizeI <= 0. || iLogSizeJ <= 0. )
+                        {
+                            continue;
+                        }
                         // get distance between the two telescopes in shower coordinates
                         fTelescopeDistanceSC = getTelescopeDistanceSC( i, j, az, el );
                         // calculate emission height [km]
                         iEmissionHeightTemp = fTelescopeDistanceSC / fImageDistance / 1.e3;
                         // weight for pairwise emission height calculation
-                        iEmissionHeightWeightTemp = 1. / (( 1. / log10( size[i] ) ) + ( 1. / log10( size[j] ) ) );
+                        iEmissionHeightWeightTemp = 1. / (( 1. / iLogSizeI ) + ( 1. / iLogSizeJ ) );
                         iEmissionHeightWeight    += iEmissionHeightWeightTemp;
                         iEmissionHeight          += iEmissionHeightTemp * iEmissionHeightWeightTemp;
                         iEmissionHeight2         += iEmissionHeightTemp * iEmissionHeightTemp * iEmissionHeightWeightTemp;
@@ -105,6 +111,9 @@ void VEmissionHeightCalculator::setTelescopePositions( vector< float > x, vector
     {
         return;
     }
+    fTelX.clear();
+    fTelY.clear();
+    fTelZ.clear();
     for( unsigned int i = 0; i < fNTel; i++ )
     {
         fTelX.push_back( x[i] );
@@ -121,6 +130,9 @@ void VEmissionHeightCalculator::setTelescopePositions( vector< float > x, vector
 void VEmissionHeightCalculator::setTelescopePositions( unsigned int ntel, double* x, double* y, double* z )
 {
     fNTel = ntel;
+    fTelX.clear();
+    fTelY.clear();
+    fTelZ.clear();
     for( unsigned int i = 0; i < ntel; i++ )
     {
         fTelX.push_back( x[i] );

src/VImageParameterCalculation.cpp

@@ -739,6 +739,10 @@ float VImageParameterCalculation::correctSizeInMuonRing()
     {
 
         float xi_tmp = ( float )i / ( float )numSteps;
+        if( xi_tmp <= 0. )
+        {
+            continue;
+        }
         float kTest_tmp = 2.0 * ngExi[i] * kRatio / xi_tmp;
 
         if( fabs( 1.0 - kTest_tmp ) < fabs( 1.0 - kTest ) )
@@ -1544,11 +1548,13 @@ vector<bool> VImageParameterCalculation::calcLL( bool iUseSums2 )
         // image centroid should not be outside of the FOV (by more than 10%)
         if( fData->getDetectorGeometry() &&  fData->getTelID() < fData->getDetectorGeometry()->getFieldofView().size() )
         {
+            fdistXmin = -1.1 * 0.5 * fData->getDetectorGeometry()->getFieldofView()[fData->getTelID()];
             fdistXmax = 1.1 * 0.5 * fData->getDetectorGeometry()->getFieldofView()[fData->getTelID()];
         }
         // should never land here
         else
         {
+            fdistXmin = -5.;
             fdistXmax = 5.;
         }
     }
@@ -1590,11 +1596,13 @@ vector<bool> VImageParameterCalculation::calcLL( bool iUseSums2 )
         // image centroid should not be outside of the FOV (by more than 10%)
         if( fData->getDetectorGeometry() &&  fData->getTelID() < fData->getDetectorGeometry()->getFieldofView().size() )
         {
+            fdistYmin = -1.1 * 0.5 * fData->getDetectorGeometry()->getFieldofView()[fData->getTelID()];
             fdistYmax = 1.1 * 0.5 * fData->getDetectorGeometry()->getFieldofView()[fData->getTelID()];
         }
         else
         {
             // should never land here
+            fdistYmin = -5.;
             fdistYmax = 5.;
         }
     }

src/VInstrumentResponseFunctionData.cpp

@@ -476,7 +476,7 @@ void VInstrumentResponseFunctionData::fill( double iWeight )
         iDiff = sqrt(( fData->Xcore - fData->MCxcore ) * ( fData->Xcore - fData->MCxcore ) +
                      ( fData->Ycore - fData->MCycore ) * ( fData->Ycore - fData->MCycore ) );
         // core error
-        iError = sqrt( fData->Xcore* fData->Xcore + fData->Xcore* fData->Xcore ) -
+        iError = sqrt( fData->Xcore* fData->Xcore + fData->Ycore* fData->Ycore ) -
                  sqrt( fData->MCxcore* fData->MCxcore + fData->MCycore* fData->MCycore );
         // relative error
         if( sqrt( fData->MCxcore * fData->MCxcore + fData->MCycore * fData->MCycore ) > 0. )

src/VInstrumentResponseFunctionRunParameter.cpp

@@ -356,6 +356,12 @@ bool VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile( stri
                     is_stream >> fRerunStereoReconstruction_3telescopes;
                 }
             }
+            else
+            {
+                cout << "VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile error: unknown run parameter key " << temp << endl;
+                cout << "line: " << is_line << endl;
+                return false;
+            }
         }
     }
     cout << "========================================" << endl << endl;

src/VStereoAnalysis.cpp

@@ -348,6 +348,20 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d
     // initialize cuts
     setCuts( fRunPara->fRunList[fHisCounter], irun );
 
+    // load XGB gamma-hadron friend tree now that the analysis type is known
+    if( fCuts && fCuts->useXGBoostCuts()
+            && fDataRun
+            && !fRunPara->fXGB_gh_file_suffix.empty() )
+    {
+        if(!fDataRun->loadGHXGBTree( fRunPara->fXGB_gh_file_suffix ) )
+        {
+            cout << "VStereoAnalysis::fillHistograms error: failed to load XGB gamma-hadron friend tree" << endl;
+            cout << "suffix: " << fRunPara->fXGB_gh_file_suffix << endl;
+            cout << "exiting..." << endl;
+            exit( EXIT_FAILURE );
+        }
+    }
+
     // define histograms
     fDirTotRun[fHisCounter]->cd();
     fHisto[fHisCounter]->setRunNumber( irun );
@@ -1967,8 +1981,9 @@ CData* VStereoAnalysis::getDataFromFile( int i_runNumber )
             false,
             false,
             iFileName,
-            fRunPara->fXGB_stereo_file_suffix,
-            fRunPara->fXGB_gh_file_suffix
+            ( fRunPara->fEnergyReconstructionMethod == 2 || fRunPara->fDirectionReconstructionMethod == 2 )
+                ? fRunPara->fXGB_stereo_file_suffix : "",
+            ""  // GH XGB tree loaded later, after analysis type is known from cut file
         );
         // read current (major) epoch from data file
         VEvndispRunParameter* i_runPara = ( VEvndispRunParameter* )fDataFile->Get( "runparameterV2" );

src/VTMVARunData.cpp

@@ -720,7 +720,7 @@ bool VTMVARunData::readConfigurationFile( char* iC )
                 // check sanity
                 if( iZenithCut_min.size() < 2 )
                 {
-                    cout << "VTMVARunData::readConfigurationFile error: need at least one zenith bin ";
+                    cout << "VTMVARunData::readConfigurationFile error: insufficient number of zenith bins ";
                     cout << iZenithCut_min.size() << endl;
                     return false;
                 }
@@ -795,7 +795,7 @@ bool VTMVARunData::fillEnergyCutData(
     // check sanity
     if( iEnergyCut_Log10TeV_min.size() < 1 )
     {
-        cout << "VTMVARunData::readConfigurationFile error: need at least two energy bins ";
+        cout << "VTMVARunData::readConfigurationFile error: insufficient number of energy bins ";
         cout << iEnergyCut_Log10TeV_min.size() << endl;
         return false;
     }

src/VTableLookupRunParameter.cpp

@@ -172,7 +172,7 @@ bool VTableLookupRunParameter::fillParameters( int argc, char* argv[] )
             }
             else
             {
-                cout << "unknown parameter, choose 1=fill or 2=read lookup tables" << endl;
+                cout << "unknown parameter, choose 1=fill or 0=read lookup tables" << endl;
                 return false;
             }
         }

src/VTraceHandler.cpp

@@ -322,6 +322,19 @@ vector< float >& VTraceHandler::getPulseTiming( int fFirst, int fLast, int fTFir
     {
         fFirst = 0;
     }
+    int nTrace = ( int )fpTrace.size();
+    if( nTrace < 2 )
+    {
+        return fpulsetiming;
+    }
+    if( fTFirst < 0 )
+    {
+        fTFirst = 0;
+    }
+    if( fTLast > nTrace - 1 )
+    {
+        fTLast = nTrace - 1;
+    }
     // reset pulse timing vector
     for( unsigned int i = 0; i < fpulsetiming.size(); i++ )
     {
@@ -343,7 +356,12 @@ vector< float >& VTraceHandler::getPulseTiming( int fFirst, int fLast, int fTFir
     // (loop backwards over pulse)
     bool bBreak = false;
     fFindPulseTiming = false;
-    for( int i = maxpos; i >= fTFirst ; i-- )
+    int i_start = maxpos;
+    if( i_start > nTrace - 2 )
+    {
+        i_start = nTrace - 2;
+    }
+    for( int i = i_start; i >= fTFirst ; i-- )
     {
         i_trace = fpTrace[i] - fPed;
         // loop over all pulse level
@@ -753,7 +771,7 @@ double VTraceHandler::calculateTraceSum_slidingWindow( unsigned int iSearchStart
         fSumWindowFirst = n - iIntegrationWindow;
         fSumWindowLast  = n;
 
-        return FADC[1];
+        return charge;
     }
 
     ////////////////////////////////////////
@@ -780,7 +798,10 @@ double VTraceHandler::calculateTraceSum_slidingWindow( unsigned int iSearchStart
                 fSumWindowLast = n;
             }
         }
-        xmax = xmax - FADC[i] + FADC[i + iIntegrationWindow];
+        if( i + iIntegrationWindow < n )
+        {
+            xmax = xmax - FADC[i] + FADC[i + iIntegrationWindow];
+        }
     }
     // arrival times (weighted average)
     float tcharge = 0.;

src/anasum.cpp

@@ -105,7 +105,7 @@ int main( int argc, char* argv[] )
     // stereo analysis (default)
     if( analysisType == 3 || analysisType == 4 )
     {
-        anasum->doStereoAnalysis(( analysisType == 3 ) || ( analysisType == 5 ) );
+        anasum->doStereoAnalysis( analysisType == 3 );
     }
     else
     {