From 79f61b1b40a18789b81cba8a2ffcce6463797ed6 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Sat, 3 Jan 2026 16:57:24 +0100 Subject: [PATCH 01/21] XGboost cuts --- inc/CData.h | 7 ++++-- inc/VGammaHadronCuts.h | 8 +++++-- src/CData.cpp | 49 +++++++++++++++++++++++++++++++--------- src/VGammaHadronCuts.cpp | 30 ++++++++++++++++++++++++ 4 files changed, 79 insertions(+), 15 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index 5551224e..7da27a1c 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -260,12 +260,15 @@ class CData TBranch* b_Xoff_intersect; TBranch* b_Yoff_intersect; - TTree* fFriendTree; //! + TTree* fStereoFriendTree; //! float Dir_Xoff; //! float Dir_Yoff; //! float Dir_Erec; //! + TTree* fGHFriendTree; //! + float GH_Gamma_Prediction; //! + bool GH_Is_Gamma; //! - CData( TTree* tree = 0, bool bMC = false, bool bShort = false, TTree* friendTree = 0 ); + CData( TTree* tree = 0, bool bMC = false, bool bShort = false, TTree* stereoTree = 0, TTree* ghTree = 0 ); virtual ~CData(); virtual Int_t GetEntry( Long64_t entry ); float get_Erec( unsigned int method = 0 ); diff --git a/inc/VGammaHadronCuts.h b/inc/VGammaHadronCuts.h index 4b900d2a..bb3ce118 100644 --- a/inc/VGammaHadronCuts.h +++ b/inc/VGammaHadronCuts.h @@ -32,7 +32,7 @@ using namespace std; //////////////////////////////////////////////////////////////////////////////// // analysis types //////////////////////////////////////////////////////////////////////////////// -enum E_AnalysisType { GEO = 0, MVAAnalysis = 1 }; +enum E_AnalysisType { GEO = 0, MVAAnalysis = 1, XGBoostAnalysis = 2 }; //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -325,11 +325,15 @@ class VGammaHadronCuts : public VAnalysisUtilities { return ( fAnalysisType == MVAAnalysis ); } + bool useXGBoostCuts() + { + return ( fAnalysisType == XGBoostAnalysis ); + } bool useOrbitalPhaseCuts() { return fUseOrbitalPhaseCuts; } - ClassDef( VGammaHadronCuts, 60 ); + ClassDef( VGammaHadronCuts, 61 ); }; #endif diff --git a/src/CData.cpp b/src/CData.cpp index 6b1b02f8..8fef4cbc 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -10,20 +10,20 @@ #include "CData.h" -CData::CData( TTree* tree, bool bMC, bool bShort, TTree* friendTree ) +CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree *ghTree ) { fMC = bMC; fShort = bShort; fVersion = 6; fTelescopeCombination = 0; - fFriendTree = friendTree; + fStereoFriendTree = friendTree; Init( tree ); - if( fFriendTree ) + if( fStereoFriendTree ) { - fFriendTree->SetBranchAddress( "Dir_Xoff", &Dir_Xoff ); - fFriendTree->SetBranchAddress( "Dir_Yoff", &Dir_Yoff ); - fFriendTree->SetBranchAddress( "Dir_Erec", &Dir_Erec ); + fStereoFriendTree->SetBranchAddress( "Dir_Xoff", &Dir_Xoff ); + fStereoFriendTree->SetBranchAddress( "Dir_Yoff", &Dir_Yoff ); + fStereoFriendTree->SetBranchAddress( "Dir_Erec", &Dir_Erec ); } else { @@ -31,6 +31,17 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* friendTree ) Dir_Yoff = -9999.; Dir_Erec = -9999.; } + fGHFriendTree = ghTree; + if( fGHFriendTree ) + { + fGHFriendTree->SetBranchAddress( "GH_Gamma_Prediction", &GH_Gamma_Prediction ); + fGHFriendTree->SetBranchAddress( "Is_Gamma_70, ", &GH_Is_Gamma ); + } + else + { + GH_Gamma_Prediction = -9999.; + GH_Is_Gamma = false; + } } @@ -53,9 +64,9 @@ Int_t CData::GetEntry( Long64_t entry ) } int a = fChain->GetEntry( entry ); - if( fFriendTree ) + if( fStereoFriendTree ) { - fFriendTree->GetEntry( entry ); + fStereoFriendTree->GetEntry( entry ); } if( fTelescopeCombination > 0 && fTelescopeCombination != 15 ) @@ -638,6 +649,22 @@ Bool_t CData::Notify() return kTRUE; } +/* + Get Gamma/hadron decision from XGB friend tree. +*/ +float CData::get_GH_Gamma_Prediction() +{ + return GH_Gamma_Prediction; +} + +/* + Get Gamma/hadron decision from XGB friend tree. +*/ +bool CData::is_GH_Gamma() +{ + return GH_Is_Gamma; +} + /* * Get energy depending on analysis method @@ -650,7 +677,7 @@ Bool_t CData::Notify() */ float CData::get_Erec( unsigned iMethod ) { - if( iMethod == 0 && fFriendTree ) + if( iMethod == 0 && fStereoFriendTree ) { return Dir_Erec; } @@ -681,7 +708,7 @@ float CData::get_Erec( unsigned iMethod ) */ float CData::get_Xoff( unsigned iMethod ) { - if( iMethod == 0 && fFriendTree ) + if( iMethod == 0 && fStereoFriendTree ) { return Dir_Xoff; } @@ -713,7 +740,7 @@ float CData::get_Xoff( unsigned iMethod ) */ float CData::get_Yoff( unsigned iMethod ) { - if( iMethod == 0 && fFriendTree ) + if( iMethod == 0 && fStereoFriendTree ) { return Dir_Yoff; } diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index f1f87a2b..1b3f3a40 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -14,6 +14,7 @@ 3: apply cuts on probabilities given by a friend to the data tree already at the level of the event quality level 4: TMVA gamma/hadron separation + 5: XGBoost gamma/hadron separation ID1: @@ -677,6 +678,14 @@ bool VGammaHadronCuts::readCuts( string i_cutfilename, int iPrint ) { fAnalysisType = MVAAnalysis; } + else if( fGammaHadronCutSelector / 10 == 5 ) + { + fAnalysisType = XGBoostAnalysis; + } + else + { + fAnalysisType = GEO; + } return true; } @@ -795,6 +804,12 @@ void VGammaHadronCuts::printCutSummary() } } } + // XGBoost cuts + if( useXGBoostCuts() ) + { + // TODO + cout << "XGBoost gamma/hadron separation" << endl; + } // other cut parameters if( fNTel == 2 ) { @@ -1139,6 +1154,21 @@ bool VGammaHadronCuts::isGamma( int i, bool bCount, bool fIsOn ) return false; } } + else if( useXGBoostCuts() ) + { + if( fDebug ) + { + cout << "VGammaHadronCuts::isGamma: applyXGBoostCut" << endl; + } + if(!applyXGBoostCut( i ) ) + { + if( bCount && fStats ) + { + fStats->updateCutCounter( VGammaHadronCutsStatistics::eIsGamma ); + } + return false; + } + } return true; } From 1445bd3da2f08c8a631ace220b53abe0a57e10ac Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 5 Jan 2026 19:41:03 +0100 Subject: [PATCH 02/21] add XGB gamma/hadron reader --- inc/CData.h | 3 ++ inc/VAnaSumRunParameter.h | 3 +- inc/VGammaHadronCuts.h | 1 + inc/VInstrumentResponseFunctionRunParameter.h | 3 +- inc/VStereoAnalysis.h | 2 - src/CData.cpp | 52 +++++++++++++++++++ src/VAnaSumRunParameter.cpp | 28 +++++++--- src/VGammaHadronCuts.cpp | 20 ++++++- ...InstrumentResponseFunctionRunParameter.cpp | 25 ++++++--- src/VStereoAnalysis.cpp | 36 +++---------- src/makeEffectiveArea.cpp | 29 +---------- 11 files changed, 126 insertions(+), 76 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index 7da27a1c..8fd2d36d 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -267,13 +267,16 @@ class CData TTree* fGHFriendTree; //! float GH_Gamma_Prediction; //! bool GH_Is_Gamma; //! + vector fXGBFiles; //! CData( TTree* tree = 0, bool bMC = false, bool bShort = false, TTree* stereoTree = 0, TTree* ghTree = 0 ); + CData( TTree* tree, bool bMC, bool bShort, string stereo_suffix, string gamma_hadron_suffix ); virtual ~CData(); virtual Int_t GetEntry( Long64_t entry ); float get_Erec( unsigned int method = 0 ); float get_Xoff( unsigned int method = 0 ); float get_Yoff( unsigned int method = 0 ); + TTree *getXGBTree( string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method = 0 ); virtual Long64_t LoadTree( Long64_t entry ); virtual void Init( TTree* tree ); diff --git a/inc/VAnaSumRunParameter.h b/inc/VAnaSumRunParameter.h index aa66c41e..96effba7 100644 --- a/inc/VAnaSumRunParameter.h +++ b/inc/VAnaSumRunParameter.h @@ -194,7 +194,8 @@ class VAnaSumRunParameter : public TNamed, public VGlobalRunParameter int fDeadTimeCalculationMethod; // XGB reconstruction - string fXGB_file_suffix; + string fXGB_stereo_file_suffix; + string fXGB_gh_file_suffix; int f2DAcceptanceMode ; // USE2DACCEPTANCE diff --git a/inc/VGammaHadronCuts.h b/inc/VGammaHadronCuts.h index 5d34a821..8acf4362 100644 --- a/inc/VGammaHadronCuts.h +++ b/inc/VGammaHadronCuts.h @@ -198,6 +198,7 @@ class VGammaHadronCuts : public VAnalysisUtilities bool applyStereoQualityCuts( unsigned int iEnergyReconstructionMethod = 0, bool bCount = false, int iEntry = 0, bool fIsOn = false ); bool applyStereoShapeCuts(); bool applyTMVACut( int i ); + bool applyXGBoostCut( int i ); double getArrayCentre_X() { diff --git a/inc/VInstrumentResponseFunctionRunParameter.h b/inc/VInstrumentResponseFunctionRunParameter.h index c0f49312..5d6d2add 100644 --- a/inc/VInstrumentResponseFunctionRunParameter.h +++ b/inc/VInstrumentResponseFunctionRunParameter.h @@ -74,7 +74,8 @@ class VInstrumentResponseFunctionRunParameter : public TNamed string fdatafile; string fMCdatafile_tree; string fMCdatafile_histo; - string fXGB_file_suffix; + string fXGB_stereo_file_suffix; + string fXGB_gh_file_suffix; double fze; int fnoise; diff --git a/inc/VStereoAnalysis.h b/inc/VStereoAnalysis.h index 2d963a3b..186fdcf7 100644 --- a/inc/VStereoAnalysis.h +++ b/inc/VStereoAnalysis.h @@ -220,8 +220,6 @@ class VStereoAnalysis CData* fDataRun; TTree* fDataRunTree; TFile* fDataFile; - TDirectory* fXGBFile; - TTree* fXGB_tree; string fInstrumentEpochMinor; vector< unsigned int > fTelToAnalyze; diff --git a/src/CData.cpp b/src/CData.cpp index 8fef4cbc..ca8bc263 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -44,6 +44,13 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree *ghTr } } +CData::CData( TTree* tree, bool bMC, bool bShort, string stereo_suffix, string gamma_hadron_suffix ) + : CData(tree, bMC, bShort, + getXGBTree(stereo_suffix, "StereoAnalysis"), + getXGBTree(gamma_hadron_suffix, "Classification")) +{ +} + CData::~CData() { @@ -51,6 +58,15 @@ CData::~CData() { return; } + for( unsigned int i = 0; i < fXGBFiles.size(); i++ ) + { + if( fXGBFiles[i] ) + { + fXGBFiles[i]->Close(); + delete fXGBFiles[i]; + } + } + fXGBFiles.clear(); delete fChain->GetCurrentFile(); } @@ -68,6 +84,10 @@ Int_t CData::GetEntry( Long64_t entry ) { fStereoFriendTree->GetEntry( entry ); } + if( fGHFriendTree ) + { + fGHFriendTree->GetEntry( entry ); + } if( fTelescopeCombination > 0 && fTelescopeCombination != 15 ) { @@ -1036,3 +1056,35 @@ void CData::initialize_3tel_reconstruction( fTelY = tel_y; fTelZ = tel_z; } + + +/* + Read XGB friend tree for gamma/hadron separation and stereo reconstruction +*/ +/* + +*/ +TTree *CData::getXGBTree(string file_suffix, string tree_name) +{ + if( file_suffix == "" || file_suffix != "None" ) + { + return 0; + } + + string iFileName = iFileName.replace(iFileName.find( ".root" ), 5, "." + file_suffix + ".root" ); + TFile *iFile = new TFile( iFileName.c_str()); + if( iFile->IsZombie() ) + { + cout << "CData Error: cannot open XGB file " << iFileName << endl; + exit( EXIT_FAILURE ); + } + TTree* iXGB_tree = ( TTree* )iFile->Get( tree_name.c_str() ); + if(!iXGB_tree ) + { + cout << "CData Error: cannot find " << tree_name << " tree in " << iFileName << endl; + exit( EXIT_FAILURE ); + } + fXGBFiles.push_back(iFile); + cout << "Adding " << tree_name << " tree from " << iFileName << endl; + return iXGB_tree; +} diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp index 151a72d5..05689c65 100644 --- a/src/VAnaSumRunParameter.cpp +++ b/src/VAnaSumRunParameter.cpp @@ -122,7 +122,8 @@ VAnaSumRunParameter::VAnaSumRunParameter() fEnergyEffectiveAreaSmoothingIterations = -1; fEnergyEffectiveAreaSmoothingThreshold = -1.; fDeadTimeCalculationMethod = 0; - fXGB_file_suffix = ""; + fXGB_stereo_file_suffix = ""; + fXGB_gh_file_suffix = ""; // background model fTMPL_fBackgroundModel = 0; @@ -598,10 +599,15 @@ int VAnaSumRunParameter::readRunParameter( string i_filename ) return 0; } } - else if( temp == "XGBFILESUFFIX" ) + else if( temp == "XGBSTEREOFILESUFFIX" ) { - fXGB_file_suffix = temp2; - if( fXGB_file_suffix == "None" ) fXGB_file_suffix = ""; + fXGB_stereo_file_suffix = temp2; + if( fXGB_stereo_file_suffix == "None" ) fXGB_stereo_file_suffix = ""; + } + else if( temp == "XGBGAMMAHADRONFILESUFFIX" ) + { + fXGB_gh_file_suffix = temp2; + if( fXGB_gh_file_suffix == "None" ) fXGB_gh_file_suffix = ""; } else if( temp == "RATEINTERVALLLENGTH" ) { @@ -1182,13 +1188,21 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) { cout << " (lookup table energy reconstruction)" << endl; } - if( fXGB_file_suffix != "" && fXGB_file_suffix != "None" ) + if( fXGB_stereo_file_suffix != "" && fXGB_stereo_file_suffix != "None" ) + { + cout << "\t XGB stereo analysis file: " << fXGB_stereo_file_suffix << endl; + } + else + { + cout << "\t no XGB stereo analysis file used" << endl; + } + if( fXGB_gh_file_suffix != "" && fXGB_gh_file_suffix != "None" ) { - cout << "\t XY direction file: " << fXGB_file_suffix << endl; + cout << "\t XGB gamma-hadron separation file: " << fXGB_gh_file_suffix << endl; } else { - cout << "\t no XY direction file used" << endl; + cout << "\t no XGB gamma-hadron separation file used" << endl; } cout << "\t dead time calculation method: "; if( fDeadTimeCalculationMethod == 0 ) diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index cfd822c3..9e7bc0e2 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -718,8 +718,7 @@ void VGammaHadronCuts::printCutSummary() // XGBoost cuts if( useXGBoostCuts() ) { - // TODO - cout << "XGBoost gamma/hadron separation" << endl; + cout << "XGBoost gamma/hadron separation with fixed 70\% signal efficiency" << endl; } // other cut parameters if( fNTel == 2 ) @@ -1110,6 +1109,23 @@ bool VGammaHadronCuts::applyTMVACut( int i ) return false; } +/* + + apply XGBoost cuts + +*/ +bool VGammaHadronCuts::applyXGBoostCut( int i ) +{ + if( fDebug ) + { + cout << "VGammaHadronCuts::applyXGBoostCut event " << i; + cout << ", prediction " << fData->GH_Gamma_Prediction; + cout << ", is gamma (70\% signal efficiency) " << fData->GH_Gamma_IsGamma; + cout << endl; + } + return fData->GH_Gamma_IsGamma; +} + /*! diff --git a/src/VInstrumentResponseFunctionRunParameter.cpp b/src/VInstrumentResponseFunctionRunParameter.cpp index 906c5993..505a105d 100644 --- a/src/VInstrumentResponseFunctionRunParameter.cpp +++ b/src/VInstrumentResponseFunctionRunParameter.cpp @@ -51,7 +51,8 @@ VInstrumentResponseFunctionRunParameter::VInstrumentResponseFunctionRunParameter fdatafile = ""; fMCdatafile_tree = ""; fMCdatafile_histo = ""; - fXGB_file_suffix = ""; + fXGB_stereo_file_suffix = ""; + fXGB_gh_file_suffix = ""; fze = 0.; fnoise = 0; @@ -236,12 +237,20 @@ bool VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile( stri is_stream >> fCutFileName; } } - else if( temp == "XGBFILESUFFIX" ) + else if( temp == "XGBSTEREOFILESUFFIX" ) { if(!( is_stream >> std::ws ).eof() ) { - is_stream >> fXGB_file_suffix; - if( fXGB_file_suffix == "None" ) fXGB_file_suffix = ""; + is_stream >> fXGB_stereo_file_suffix; + if( fXGB_stereo_file_suffix == "None" ) fXGB_stereo_file_suffix = ""; + } + } + else if( temp == "XGBGAMMAHADRONFILESUFFIX" ) + { + if(!( is_stream >> std::ws ).eof() ) + { + is_stream >> fXGB_gh_file_suffix; + if( fXGB_gh_file_suffix == "None" ) fXGB_gh_file_suffix = ""; } } // * SCATTERMODE @@ -610,9 +619,13 @@ void VInstrumentResponseFunctionRunParameter::print() { cout << " MC histograms: " << fMCdatafile_histo << endl; } - if( fXGB_file_suffix.size() > 0 ) + if( fXGB_stereo_file_suffix.size() > 0 ) + { + cout << " XGB file suffix: " << fXGB_stereo_file_suffix << endl; + } + if( fXGB_gh_file_suffix.size() > 0 ) { - cout << " XGB file suffix: " << fXGB_file_suffix << endl; + cout << " XGB Gamma/Hadron file suffix: " << fXGB_gh_file_suffix << endl; } if( fInstrumentEpoch != "NOT_SET" ) { diff --git a/src/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index de03062c..adae0133 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -11,7 +11,6 @@ VStereoAnalysis::VStereoAnalysis( bool ion, string i_hsuffix, VAnaSumRunParamete fDebug = false; fDataFile = 0; - fXGBFile = 0; fXGB_tree = 0; fInstrumentEpochMinor = "NOT_SET"; fDirTot = iDirTot; @@ -1963,34 +1962,13 @@ CData* VStereoAnalysis::getDataFromFile( int i_runNumber ) cout << "exiting..." << endl; exit( EXIT_FAILURE ); } - fXGB_tree = 0; - if( fRunPara->fXGB_file_suffix != "" && fRunPara->fXGB_file_suffix != "None" ) - { - fXGBFile = new TFile( iFileName.replace( - iFileName.find( ".root" ), 5, - "." + fRunPara->fXGB_file_suffix + ".root" ).c_str() - ); - if( fXGBFile->IsZombie() ) - { - cout << "VStereoAnalysis::getDataFromFile() warning: cannot open DispDirection file " - << iFileName << endl; - exit( EXIT_FAILURE ); - } - else - { - fXGB_tree = ( TTree* )fXGBFile->Get( "StereoAnalysis" ); - // backwards compatibility - if(!fXGB_tree ) fXGB_tree = ( TTree* )fXGBFile->Get( "DispDirection" ); - if(!fXGB_tree ) - { - cout << "VStereoAnalysis::getDataFromFile() error: cannot find stereo analysis tree in " - << fXGBFile->GetName() << endl; - exit( EXIT_FAILURE ); - } - cout << "VStereoAnalysis::getDataFromFile(): adding DispDirection from " << fXGBFile->GetName() << endl; - } - } - c = new CData( fDataRunTree, false, false, fXGB_tree ); + c = new CData( + fDataRunTree, + false, + false, + fRunPara->fXGB_stereo_file_suffix, + fRunPara->fXGB_gamma_hadron_file_suffix + ); // read current (major) epoch from data file VEvndispRunParameter* i_runPara = ( VEvndispRunParameter* )fDataFile->Get( "runparameterV2" ); if( i_runPara ) diff --git a/src/makeEffectiveArea.cpp b/src/makeEffectiveArea.cpp index 0b281939..b057e7c0 100644 --- a/src/makeEffectiveArea.cpp +++ b/src/makeEffectiveArea.cpp @@ -180,34 +180,7 @@ int main( int argc, char* argv[] ) exit( EXIT_FAILURE ); } - // XGB file - TFile *fXGBFile = 0; - TTree *fXGB_tree = 0; - if( fRunPara->fXGB_file_suffix != "" && fRunPara->fXGB_file_suffix != "None" ) - { - string xgb_file_name = fRunPara->fdatafile; - xgb_file_name.replace( fRunPara->fdatafile.find( ".root" ), 5, "." + fRunPara->fXGB_file_suffix + ".root" ); - fXGBFile = new TFile( xgb_file_name.c_str() ); - if( fXGBFile->IsZombie() ) - { - cout << "Error: cannot open XGB file " << xgb_file_name << endl; - exit( EXIT_FAILURE ); - } - else - { - fXGB_tree = ( TTree* )fXGBFile->Get( "StereoAnalysis" ); - // backwards compatibility - if(!fXGB_tree ) fXGB_tree = ( TTree* )fXGBFile->Get( "DispDirection" ); - if(!fXGB_tree ) - { - cout << "Error: cannot find stereo analysis tree in " << fXGBFile->GetName() << endl; - exit( EXIT_FAILURE ); - } - cout << "Adding XGB DispDirection from " << fXGBFile->GetName() << endl; - } - } - - CData d( c, true, false, fXGB_tree ); + CData d( c, true, false, fRunPara->fXGB_stereo_file_suffix, fRunPara->fXGB_gh_file_suffix ); d.initialize_3tel_reconstruction( fRunPara->fRerunStereoReconstruction_3telescopes, fRunPara->fRerunStereoReconstruction_minAngle, From 8ede0c5759c082666b04dbeb7c5242347ccce5e6 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 5 Jan 2026 19:46:22 +0100 Subject: [PATCH 03/21] typo --- src/VGammaHadronCuts.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 9e7bc0e2..1722d3ac 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -1123,7 +1123,7 @@ bool VGammaHadronCuts::applyXGBoostCut( int i ) cout << ", is gamma (70\% signal efficiency) " << fData->GH_Gamma_IsGamma; cout << endl; } - return fData->GH_Gamma_IsGamma; + return fData->GH_Is_Gamma; } From 120e9a0a23bec7f0f311fb814f7d4aab69bcfd29 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 5 Jan 2026 19:49:13 +0100 Subject: [PATCH 04/21] typo --- src/VGammaHadronCuts.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 1722d3ac..0f32c1ed 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -1120,7 +1120,7 @@ bool VGammaHadronCuts::applyXGBoostCut( int i ) { cout << "VGammaHadronCuts::applyXGBoostCut event " << i; cout << ", prediction " << fData->GH_Gamma_Prediction; - cout << ", is gamma (70\% signal efficiency) " << fData->GH_Gamma_IsGamma; + cout << ", is gamma (70\% signal efficiency) " << fData->GH_Is_Gamma; cout << endl; } return fData->GH_Is_Gamma; From 021dca67c7bcdcab8aa783db28e83476980f9d6b Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Wed, 7 Jan 2026 08:50:32 +0100 Subject: [PATCH 05/21] compile fixes --- inc/CData.h | 2 ++ src/CData.cpp | 2 +- src/VStereoAnalysis.cpp | 3 +-- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index 8fd2d36d..ee43d7ea 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -279,6 +279,8 @@ class CData TTree *getXGBTree( string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method = 0 ); virtual Long64_t LoadTree( Long64_t entry ); + float get_GH_Gamma_Prediction(); + bool is_GH_Gamma(); virtual void Init( TTree* tree ); virtual Bool_t Notify(); bool isMC() diff --git a/src/CData.cpp b/src/CData.cpp index ca8bc263..48e8a1e7 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -17,7 +17,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree *ghTr fVersion = 6; fTelescopeCombination = 0; - fStereoFriendTree = friendTree; + fStereoFriendTree = stereoTree; Init( tree ); if( fStereoFriendTree ) { diff --git a/src/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index adae0133..d78ccacb 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -11,7 +11,6 @@ VStereoAnalysis::VStereoAnalysis( bool ion, string i_hsuffix, VAnaSumRunParamete fDebug = false; fDataFile = 0; - fXGB_tree = 0; fInstrumentEpochMinor = "NOT_SET"; fDirTot = iDirTot; fDirTotRun = iDirRun; @@ -1967,7 +1966,7 @@ CData* VStereoAnalysis::getDataFromFile( int i_runNumber ) false, false, fRunPara->fXGB_stereo_file_suffix, - fRunPara->fXGB_gamma_hadron_file_suffix + fRunPara->fXGB_gh_file_suffix ); // read current (major) epoch from data file VEvndispRunParameter* i_runPara = ( VEvndispRunParameter* )fDataFile->Get( "runparameterV2" ); From 40a86a4a0069162922fe166d5d08ef9c9db7c08b Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Wed, 7 Jan 2026 08:50:55 +0100 Subject: [PATCH 06/21] astyle --- inc/CData.h | 2 +- src/CData.cpp | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index ee43d7ea..f0d2741a 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -276,7 +276,7 @@ class CData float get_Erec( unsigned int method = 0 ); float get_Xoff( unsigned int method = 0 ); float get_Yoff( unsigned int method = 0 ); - TTree *getXGBTree( string suffix, string tree_name ); + TTree* getXGBTree( string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method = 0 ); virtual Long64_t LoadTree( Long64_t entry ); float get_GH_Gamma_Prediction(); diff --git a/src/CData.cpp b/src/CData.cpp index 48e8a1e7..65af82d1 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -10,7 +10,7 @@ #include "CData.h" -CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree *ghTree ) +CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree* ghTree ) { fMC = bMC; fShort = bShort; @@ -45,9 +45,9 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree *ghTr } CData::CData( TTree* tree, bool bMC, bool bShort, string stereo_suffix, string gamma_hadron_suffix ) - : CData(tree, bMC, bShort, - getXGBTree(stereo_suffix, "StereoAnalysis"), - getXGBTree(gamma_hadron_suffix, "Classification")) + : CData( tree, bMC, bShort, + getXGBTree( stereo_suffix, "StereoAnalysis" ), + getXGBTree( gamma_hadron_suffix, "Classification" ) ) { } @@ -1064,15 +1064,15 @@ void CData::initialize_3tel_reconstruction( /* */ -TTree *CData::getXGBTree(string file_suffix, string tree_name) +TTree* CData::getXGBTree( string file_suffix, string tree_name ) { if( file_suffix == "" || file_suffix != "None" ) { return 0; } - string iFileName = iFileName.replace(iFileName.find( ".root" ), 5, "." + file_suffix + ".root" ); - TFile *iFile = new TFile( iFileName.c_str()); + string iFileName = iFileName.replace( iFileName.find( ".root" ), 5, "." + file_suffix + ".root" ); + TFile *iFile = new TFile( iFileName.c_str() ); if( iFile->IsZombie() ) { cout << "CData Error: cannot open XGB file " << iFileName << endl; @@ -1084,7 +1084,7 @@ TTree *CData::getXGBTree(string file_suffix, string tree_name) cout << "CData Error: cannot find " << tree_name << " tree in " << iFileName << endl; exit( EXIT_FAILURE ); } - fXGBFiles.push_back(iFile); + fXGBFiles.push_back( iFile ); cout << "Adding " << tree_name << " tree from " << iFileName << endl; return iXGB_tree; } From 6a662189cb4f15709ed16c39084924f4c990468b Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Sat, 10 Jan 2026 09:56:56 +0100 Subject: [PATCH 07/21] docstrings --- src/CData.cpp | 3 --- src/VGammaHadronCuts.cpp | 2 ++ 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/CData.cpp b/src/CData.cpp index 65af82d1..28ea35cb 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -1060,9 +1060,6 @@ void CData::initialize_3tel_reconstruction( /* Read XGB friend tree for gamma/hadron separation and stereo reconstruction -*/ -/* - */ TTree* CData::getXGBTree( string file_suffix, string tree_name ) { diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 0f32c1ed..8022135a 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -27,6 +27,8 @@ cut selector = 0 : apply MSCW/MSCL cuts (default) cut selector = 22 : apply event probability cuts cut selector = 10 : apply cuts from a tree AND apply MSCW/MSCL cuts + cut selector = 42: apply TMVA gamma/hadron separation AND apply mean width/length cuts (default TMVA-BDT cut) + cut selector = 52: apply XGBoost gamma/hadron separation AND apply mean width/length cuts */ From c1ea56c789e4f60789117cc0afdf1cec9c1d8c83 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sat, 10 Jan 2026 11:18:32 +0100 Subject: [PATCH 08/21] fix bug in cut selection --- src/VGammaHadronCuts.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 8022135a..8dfe9e5f 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -1026,7 +1026,6 @@ bool VGammaHadronCuts::applyEnergyReconstructionQualityCuts( unsigned int iEnerg */ bool VGammaHadronCuts::isGamma( int i, bool bCount, bool fIsOn ) { - ///////////////////////////////////////////////////////////////////////////// // apply box cuts (e.g. MSCW/MSCL or MWR/MLR) if( fGammaHadronCutSelector % 10 <= 3 ) @@ -1051,7 +1050,7 @@ bool VGammaHadronCuts::isGamma( int i, bool bCount, bool fIsOn ) } ///////////////////////////////////////////////////////////////////////////// // apply cut using TMVA reader - else if( useTMVACuts() ) + if( useTMVACuts() ) { if( fDebug ) { From 4ac2276ce6f2014791ef65032349c74c2db09374 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sat, 10 Jan 2026 12:12:48 +0100 Subject: [PATCH 09/21] fix reading of gh trees --- inc/CData.h | 17 +++---- src/CData.cpp | 91 ++++++++++++++++++++++--------------- src/VAnaSumRunParameter.cpp | 4 +- src/VGammaHadronCuts.cpp | 4 +- src/VStereoAnalysis.cpp | 1 + 5 files changed, 69 insertions(+), 48 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index f0d2741a..3d00f517 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -261,22 +261,23 @@ class CData TBranch* b_Yoff_intersect; TTree* fStereoFriendTree; //! - float Dir_Xoff; //! - float Dir_Yoff; //! - float Dir_Erec; //! + float Dir_Xoff; //! + float Dir_Yoff; //! + float Dir_Erec; //! TTree* fGHFriendTree; //! - float GH_Gamma_Prediction; //! - bool GH_Is_Gamma; //! - vector fXGBFiles; //! + float GH_Gamma_Prediction; //! + UChar_t GH_Is_Gamma; //! + vector fXGBFiles; //! CData( TTree* tree = 0, bool bMC = false, bool bShort = false, TTree* stereoTree = 0, TTree* ghTree = 0 ); - CData( TTree* tree, bool bMC, bool bShort, string stereo_suffix, string gamma_hadron_suffix ); + CData( TTree* tree, bool bMC, bool bShort, string file_name, string stereo_suffix, string gamma_hadron_suffix ); virtual ~CData(); virtual Int_t GetEntry( Long64_t entry ); float get_Erec( unsigned int method = 0 ); float get_Xoff( unsigned int method = 0 ); float get_Yoff( unsigned int method = 0 ); - TTree* getXGBTree( string suffix, string tree_name ); + void initialize_xgb_tree(TTree* stereoTree, TTree* ghTree); + TTree* getXGBTree( string file_name, string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method = 0 ); virtual Long64_t LoadTree( Long64_t entry ); float get_GH_Gamma_Prediction(); diff --git a/src/CData.cpp b/src/CData.cpp index 28ea35cb..f67539e6 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -16,39 +16,27 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree* ghTr fShort = bShort; fVersion = 6; fTelescopeCombination = 0; + Init( tree ); fStereoFriendTree = stereoTree; - Init( tree ); - if( fStereoFriendTree ) - { - fStereoFriendTree->SetBranchAddress( "Dir_Xoff", &Dir_Xoff ); - fStereoFriendTree->SetBranchAddress( "Dir_Yoff", &Dir_Yoff ); - fStereoFriendTree->SetBranchAddress( "Dir_Erec", &Dir_Erec ); - } - else - { - Dir_Xoff = -9999.; - Dir_Yoff = -9999.; - Dir_Erec = -9999.; - } fGHFriendTree = ghTree; - if( fGHFriendTree ) - { - fGHFriendTree->SetBranchAddress( "GH_Gamma_Prediction", &GH_Gamma_Prediction ); - fGHFriendTree->SetBranchAddress( "Is_Gamma_70, ", &GH_Is_Gamma ); - } - else - { - GH_Gamma_Prediction = -9999.; - GH_Is_Gamma = false; - } + + initialize_xgb_tree(fStereoFriendTree, fGHFriendTree); } -CData::CData( TTree* tree, bool bMC, bool bShort, string stereo_suffix, string gamma_hadron_suffix ) - : CData( tree, bMC, bShort, - getXGBTree( stereo_suffix, "StereoAnalysis" ), - getXGBTree( gamma_hadron_suffix, "Classification" ) ) + + +CData::CData( TTree* tree, bool bMC, bool bShort, string file_name, string stereo_suffix, string gamma_hadron_suffix ) { + fMC = bMC; + fShort = bShort; + fVersion = 6; + fTelescopeCombination = 0; + Init( tree ); + + fStereoFriendTree = getXGBTree( file_name, stereo_suffix, "StereoAnalysis" ); + fGHFriendTree = getXGBTree( file_name, gamma_hadron_suffix, "Classification" ); + initialize_xgb_tree(fStereoFriendTree, fGHFriendTree); } @@ -1061,27 +1049,58 @@ void CData::initialize_3tel_reconstruction( /* Read XGB friend tree for gamma/hadron separation and stereo reconstruction */ -TTree* CData::getXGBTree( string file_suffix, string tree_name ) +TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name ) { - if( file_suffix == "" || file_suffix != "None" ) + if( file_suffix == "" || file_suffix == "None" ) { return 0; } - string iFileName = iFileName.replace( iFileName.find( ".root" ), 5, "." + file_suffix + ".root" ); - TFile *iFile = new TFile( iFileName.c_str() ); - if( iFile->IsZombie() ) + file_name = file_name.replace( file_name.find( ".root" ), 5, "." + file_suffix + ".root" ); + TFile *iFile = TFile::Open( file_name.c_str() ); + if( !iFile || iFile->IsZombie() ) { - cout << "CData Error: cannot open XGB file " << iFileName << endl; + cout << "CData Error: cannot open XGB file " << file_name << endl; exit( EXIT_FAILURE ); } + fXGBFiles.push_back( iFile ); + TTree* iXGB_tree = ( TTree* )iFile->Get( tree_name.c_str() ); if(!iXGB_tree ) { - cout << "CData Error: cannot find " << tree_name << " tree in " << iFileName << endl; + cout << "CData Error: cannot find " << tree_name << " tree in " << file_name << endl; exit( EXIT_FAILURE ); } - fXGBFiles.push_back( iFile ); - cout << "Adding " << tree_name << " tree from " << iFileName << endl; + cout << "Adding " << tree_name << " tree from " << file_name << endl; return iXGB_tree; } + +/* + * Initialize XGB trees as kind of friends + * +*/ +void CData::initialize_xgb_tree(TTree* stereoTree, TTree* ghTree ) +{ + if( fStereoFriendTree ) + { + fStereoFriendTree->SetBranchAddress( "Dir_Xoff", &Dir_Xoff ); + fStereoFriendTree->SetBranchAddress( "Dir_Yoff", &Dir_Yoff ); + fStereoFriendTree->SetBranchAddress( "Dir_Erec", &Dir_Erec ); + } + else + { + Dir_Xoff = -9999.; + Dir_Yoff = -9999.; + Dir_Erec = -9999.; + } + if( fGHFriendTree ) + { + fGHFriendTree->SetBranchAddress( "Gamma_Prediction", &GH_Gamma_Prediction ); + fGHFriendTree->SetBranchAddress( "Is_Gamma_60", &GH_Is_Gamma ); + } + else + { + GH_Gamma_Prediction = -9999.; + GH_Is_Gamma = false; + } +} diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp index 05689c65..27c78674 100644 --- a/src/VAnaSumRunParameter.cpp +++ b/src/VAnaSumRunParameter.cpp @@ -1190,7 +1190,7 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) } if( fXGB_stereo_file_suffix != "" && fXGB_stereo_file_suffix != "None" ) { - cout << "\t XGB stereo analysis file: " << fXGB_stereo_file_suffix << endl; + cout << "\t XGB stereo analysis file suffix: " << fXGB_stereo_file_suffix << endl; } else { @@ -1198,7 +1198,7 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) } if( fXGB_gh_file_suffix != "" && fXGB_gh_file_suffix != "None" ) { - cout << "\t XGB gamma-hadron separation file: " << fXGB_gh_file_suffix << endl; + cout << "\t XGB gamma-hadron separation file suffix: " << fXGB_gh_file_suffix << endl; } else { diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 8dfe9e5f..e3cfc1c3 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -1121,10 +1121,10 @@ bool VGammaHadronCuts::applyXGBoostCut( int i ) { cout << "VGammaHadronCuts::applyXGBoostCut event " << i; cout << ", prediction " << fData->GH_Gamma_Prediction; - cout << ", is gamma (70\% signal efficiency) " << fData->GH_Is_Gamma; + cout << ", is gamma (70\% signal efficiency) " << (bool)fData->GH_Is_Gamma; cout << endl; } - return fData->GH_Is_Gamma; + return (bool)fData->GH_Is_Gamma; } diff --git a/src/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index d78ccacb..7507fd34 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -1965,6 +1965,7 @@ CData* VStereoAnalysis::getDataFromFile( int i_runNumber ) fDataRunTree, false, false, + iFileName, fRunPara->fXGB_stereo_file_suffix, fRunPara->fXGB_gh_file_suffix ); From 1d444347603dd379b02d94dfa61ae43b9b103424 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sat, 10 Jan 2026 12:15:00 +0100 Subject: [PATCH 10/21] constructor --- src/makeEffectiveArea.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/makeEffectiveArea.cpp b/src/makeEffectiveArea.cpp index b057e7c0..743ed265 100644 --- a/src/makeEffectiveArea.cpp +++ b/src/makeEffectiveArea.cpp @@ -175,12 +175,12 @@ int main( int argc, char* argv[] ) TChain* c = new TChain( "data" ); if(!c->Add( fRunPara->fdatafile.c_str(), -1 ) ) { - cout << "Error while trying to add mscw data tree from file " << fRunPara->fdatafile << endl; + cout << "Error while trying to add mscw data tree from file " << fRunPara->fdatafile << endl; cout << "exiting..." << endl; exit( EXIT_FAILURE ); } - CData d( c, true, false, fRunPara->fXGB_stereo_file_suffix, fRunPara->fXGB_gh_file_suffix ); + CData d( c, true, false, fRunPara->fdatafile, fRunPara->fXGB_stereo_file_suffix, fRunPara->fXGB_gh_file_suffix ); d.initialize_3tel_reconstruction( fRunPara->fRerunStereoReconstruction_3telescopes, fRunPara->fRerunStereoReconstruction_minAngle, From 9bfd7ac8537588881c3f25281d83570f49bd81fd Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sat, 10 Jan 2026 12:20:08 +0100 Subject: [PATCH 11/21] use 70 --- src/CData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CData.cpp b/src/CData.cpp index f67539e6..dc510664 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -1096,7 +1096,7 @@ void CData::initialize_xgb_tree(TTree* stereoTree, TTree* ghTree ) if( fGHFriendTree ) { fGHFriendTree->SetBranchAddress( "Gamma_Prediction", &GH_Gamma_Prediction ); - fGHFriendTree->SetBranchAddress( "Is_Gamma_60", &GH_Is_Gamma ); + fGHFriendTree->SetBranchAddress( "Is_Gamma_70", &GH_Is_Gamma ); } else { From 75fc4992aea196b3c6f39c69ba93ee363fdc88b9 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 11:38:59 +0100 Subject: [PATCH 12/21] improved printout --- src/CData.cpp | 8 ++++---- src/VGammaHadronCuts.cpp | 17 ++++++++--------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/CData.cpp b/src/CData.cpp index dc510664..5f59fb32 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -21,7 +21,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree* ghTr fStereoFriendTree = stereoTree; fGHFriendTree = ghTree; - initialize_xgb_tree(fStereoFriendTree, fGHFriendTree); + initialize_xgb_tree( fStereoFriendTree, fGHFriendTree ); } @@ -36,7 +36,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, string file_name, string stere fStereoFriendTree = getXGBTree( file_name, stereo_suffix, "StereoAnalysis" ); fGHFriendTree = getXGBTree( file_name, gamma_hadron_suffix, "Classification" ); - initialize_xgb_tree(fStereoFriendTree, fGHFriendTree); + initialize_xgb_tree( fStereoFriendTree, fGHFriendTree ); } @@ -1058,7 +1058,7 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name file_name = file_name.replace( file_name.find( ".root" ), 5, "." + file_suffix + ".root" ); TFile *iFile = TFile::Open( file_name.c_str() ); - if( !iFile || iFile->IsZombie() ) + if(!iFile || iFile->IsZombie() ) { cout << "CData Error: cannot open XGB file " << file_name << endl; exit( EXIT_FAILURE ); @@ -1079,7 +1079,7 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name * Initialize XGB trees as kind of friends * */ -void CData::initialize_xgb_tree(TTree* stereoTree, TTree* ghTree ) +void CData::initialize_xgb_tree( TTree* stereoTree, TTree* ghTree ) { if( fStereoFriendTree ) { diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index e3cfc1c3..1d4e21df 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -702,10 +702,15 @@ void VGammaHadronCuts::printCutSummary() cout << endl; cout << "Orbital Phase bin ( " << fOrbitalPhase_min << ", " << fOrbitalPhase_max << " )"; } + // other cut parameters + if( fNTel == 2 ) + { + cout << ", size > " << fCut_Size_min; + } + cout << endl; // TMVA cuts if( useTMVACuts() ) { - cout << endl; cout << "TMVA gamma/hadron separation with MVA method " << fTMVA_MVAMethod; cout << endl; cout << "weight files: " << fTMVAWeightFile; @@ -722,12 +727,6 @@ void VGammaHadronCuts::printCutSummary() { cout << "XGBoost gamma/hadron separation with fixed 70\% signal efficiency" << endl; } - // other cut parameters - if( fNTel == 2 ) - { - cout << ", size > " << fCut_Size_min; - } - cout << endl; cout << "Fiducial area (camera) < " << fCut_CameraFiducialSize_max << " deg, "; cout << " stereo reconstruction: " << fCut_Chi2_min << " <= sChi2 <= " << fCut_Chi2_max << endl; cout << "Energy reconstruction: "; @@ -1121,10 +1120,10 @@ bool VGammaHadronCuts::applyXGBoostCut( int i ) { cout << "VGammaHadronCuts::applyXGBoostCut event " << i; cout << ", prediction " << fData->GH_Gamma_Prediction; - cout << ", is gamma (70\% signal efficiency) " << (bool)fData->GH_Is_Gamma; + cout << ", is gamma (70\% signal efficiency) " << ( bool )fData->GH_Is_Gamma; cout << endl; } - return (bool)fData->GH_Is_Gamma; + return ( bool )fData->GH_Is_Gamma; } From 43a92a8886ba0983301866832ab92f5b242dac60 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 11:59:17 +0100 Subject: [PATCH 13/21] clean up --- src/VGammaHadronCuts.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 1d4e21df..df5a14f8 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -1610,10 +1610,6 @@ void VGammaHadronCuts::terminate() { fTMVAEvaluatorResults->Write(); } - else - { - cout << "No TMVAEvaluator Results." << endl; - } Write(); } From 798b0f7ba9a4df5a37e8d8de5d078cc4f2a84478 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 16:04:51 +0100 Subject: [PATCH 14/21] Introduce energy and direction reconstruction methods --- inc/CData.h | 4 +- inc/VAnaSumRunParameter.h | 3 + inc/VDataMCComparision.h | 9 ++ inc/VEffectiveAreaCalculator.h | 7 +- inc/VGammaHadronCuts.h | 24 ++-- inc/VInstrumentResponseFunction.h | 4 +- inc/VInstrumentResponseFunctionData.h | 5 + inc/VInstrumentResponseFunctionRunParameter.h | 1 + inc/VRadialAcceptance.h | 7 +- inc/VTMVARunDataEnergyCut.h | 3 +- src/CData.cpp | 103 ++++++++++------- src/VAnaSumRunParameter.cpp | 55 +++++++-- src/VDataMCComparision.cpp | 13 ++- src/VEffectiveAreaCalculator.cpp | 56 ++++----- src/VGammaHadronCuts.cpp | 109 +++--------------- src/VInstrumentResponseFunction.cpp | 14 +-- src/VInstrumentResponseFunctionData.cpp | 10 +- ...InstrumentResponseFunctionRunParameter.cpp | 11 ++ src/VRadialAcceptance.cpp | 2 +- src/VStereoAnalysis.cpp | 23 ++-- src/compareDatawithMC.cpp | 7 ++ src/makeEffectiveArea.cpp | 4 +- src/makeRadialAcceptance.cpp | 2 +- 23 files changed, 240 insertions(+), 236 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index 3d00f517..75944dcd 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -274,9 +274,11 @@ class CData virtual ~CData(); virtual Int_t GetEntry( Long64_t entry ); float get_Erec( unsigned int method = 0 ); + float get_ErecChi2( unsigned int method = 0 ); + float get_ErecdE( unsigned int method = 0 ); float get_Xoff( unsigned int method = 0 ); float get_Yoff( unsigned int method = 0 ); - void initialize_xgb_tree(TTree* stereoTree, TTree* ghTree); + void initialize_xgb_tree( TTree* stereoTree, TTree* ghTree ); TTree* getXGBTree( string file_name, string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method = 0 ); virtual Long64_t LoadTree( Long64_t entry ); diff --git a/inc/VAnaSumRunParameter.h b/inc/VAnaSumRunParameter.h index 96effba7..0d971e43 100644 --- a/inc/VAnaSumRunParameter.h +++ b/inc/VAnaSumRunParameter.h @@ -190,6 +190,9 @@ class VAnaSumRunParameter : public TNamed, public VGlobalRunParameter double fEnergyEffectiveAreaSmoothingThreshold; vector< double > fMCZe; // zenith angle interval for Monte Carlo + // direction reconstruction + unsigned int fDirectionReconstructionMethod; + // dead time calculation method int fDeadTimeCalculationMethod; diff --git a/inc/VDataMCComparision.h b/inc/VDataMCComparision.h index 73b9108c..0db2172f 100644 --- a/inc/VDataMCComparision.h +++ b/inc/VDataMCComparision.h @@ -90,6 +90,10 @@ class VDataMCComparision bool fCalculateMVAValues; string fEpochATM; + // reconstruction methods + unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; + // lists with all histograms TList* hisList; vector hTel; @@ -143,6 +147,11 @@ class VDataMCComparision { fShowerMaxZe_deg = iZe; } + void setStereoReconstructionMethod( unsigned int iEnergyMethod = 0, unsigned int iDirectionMethod = 0 ) + { + fEnergyReconstructionMethod = iEnergyMethod; + fDirectionReconstructionMethod = iDirectionMethod; + } bool setTelescopeCoordinates( double x, double y, double z = 0. ); void setWobbleFromDataTree() { diff --git a/inc/VEffectiveAreaCalculator.h b/inc/VEffectiveAreaCalculator.h index 3aeb8be3..d9d16295 100644 --- a/inc/VEffectiveAreaCalculator.h +++ b/inc/VEffectiveAreaCalculator.h @@ -113,7 +113,6 @@ class VEffectiveAreaCalculator VInstrumentResponseFunctionRunParameter* fRunPara; VGammaHadronCuts* fCuts; - bool fIgnoreEnergyReconstruction; bool fIsotropicArrivalDirections; // effective area calculation @@ -316,7 +315,7 @@ class VEffectiveAreaCalculator ~VEffectiveAreaCalculator(); void cleanup(); - bool fill( TH1D* hE0mc, CData* d, VEffectiveAreaCalculatorMCHistograms* iMC_histo, unsigned int iMethod ); + bool fill( TH1D* hE0mc, CData* d, VEffectiveAreaCalculatorMCHistograms* iMC_histo, unsigned int iEnergyReconstructionMethod, unsigned int iDirectionReconstructionMethod ); TH1D* getHistogramhEmc(); TGraphErrors* getMeanSystematicErrorHistogram(); TTree* getTree() @@ -363,10 +362,6 @@ class VEffectiveAreaCalculator { fEffectiveAreaVsEnergyMC = iMC; } - void setIgnoreEnergyReconstructionCuts( bool iB = false ) - { - fIgnoreEnergyReconstruction = iB; - } void setIsotropicArrivalDirections( bool iB = false ) { fIsotropicArrivalDirections = iB; diff --git a/inc/VGammaHadronCuts.h b/inc/VGammaHadronCuts.h index 8acf4362..b11a5285 100644 --- a/inc/VGammaHadronCuts.h +++ b/inc/VGammaHadronCuts.h @@ -107,6 +107,10 @@ class VGammaHadronCuts : public VAnalysisUtilities double fTMVA_EvaluationResult; VTMVAEvaluatorResults* fTMVAEvaluatorResults; + // reconstruction methods + unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; + // orbital phase analysis TFile* fPhaseCut_File; //! TTree* fPhaseCut_Tree; //! @@ -186,8 +190,8 @@ class VGammaHadronCuts : public VAnalysisUtilities VGammaHadronCuts(); ~VGammaHadronCuts(); - bool applyDirectionCuts( unsigned int iEnergyReconstructionMethod = 0, bool bCount = false, double x0 = -99999., double y0 = -99999. ); - bool applyEnergyReconstructionQualityCuts( unsigned int iEnergyReconstructionMethod = 0, bool bCount = false ); + bool applyDirectionCuts( bool bCount = false, double x0 = -99999., double y0 = -99999. ); + bool applyEnergyReconstructionQualityCuts( bool bCount = false ); bool applyInsideFiducialAreaCut( bool bCount = false ); bool applyInsideFiducialAreaCut( float Xoff, float Yoff, bool bCount = false ); bool applyMCXYoffCut( double x, double y, bool bCount = false ); @@ -195,7 +199,7 @@ class VGammaHadronCuts : public VAnalysisUtilities bool applyMeanStereoShapeCuts(); bool applyMeanScaledStereoShapeCuts(); bool applyPhaseCut( int i ); - bool applyStereoQualityCuts( unsigned int iEnergyReconstructionMethod = 0, bool bCount = false, int iEntry = 0, bool fIsOn = false ); + bool applyStereoQualityCuts( bool bCount = false, int iEntry = 0, bool fIsOn = false ); bool applyStereoShapeCuts(); bool applyTMVACut( int i ); bool applyXGBoostCut( int i ); @@ -208,13 +212,6 @@ class VGammaHadronCuts : public VAnalysisUtilities { return fArrayCentre_Y; } - double getReconstructedEnergy( unsigned int iEnergyReconstructionMethod = 0 ); - double getReconstructedEnergyChi2( unsigned int iEnergyReconstructionMethod = 0 ); - double getReconstructedEnergydE( unsigned int iEnergyReconstructionMethod = 0. ); - double getReconstructedXoff(); - double getReconstructedYoff(); - double getReconstructedXcore(); - double getReconstructedYcore(); int getGammaHadronCutSelector() { return fGammaHadronCutSelector; @@ -254,7 +251,7 @@ class VGammaHadronCuts : public VAnalysisUtilities { return fTMVAEvaluatorResults; } - void initialize(); + void initialize( unsigned int iEnergyMethod, unsigned int iDirectionMethod ); bool isGamma( int i = 0, bool bCount = false, bool fIsOn = true ); bool isMCCuts() { @@ -305,6 +302,11 @@ class VGammaHadronCuts : public VAnalysisUtilities fArrayCentre_X = iX; fArrayCentre_Y = iY; } + void setStereoReconstructionMethod( unsigned int iEnergyMethod = 0, unsigned int iDirectionMethod = 0 ) + { + fEnergyReconstructionMethod = iEnergyMethod; + fDirectionReconstructionMethod = iDirectionMethod; + } void setTelToAnalyze( vector< unsigned int > iTelToAnalyze ) { fTelToAnalyze = iTelToAnalyze; diff --git a/inc/VInstrumentResponseFunction.h b/inc/VInstrumentResponseFunction.h index 796268d4..d94c972e 100644 --- a/inc/VInstrumentResponseFunction.h +++ b/inc/VInstrumentResponseFunction.h @@ -47,8 +47,6 @@ class VInstrumentResponseFunction // histograms are not re-filled but duplicated unsigned int fDuplicationID; - unsigned int fEnergyReconstructionMethod; - // histograms and data vector< vector< VInstrumentResponseFunctionData* > > fIRFData; @@ -115,7 +113,7 @@ class VInstrumentResponseFunction bool initialize( string iName, string iType, unsigned int iNTel, double iMCMaxCoreRadius, double iZe, int iNoise, double iPedvars, double iXoff, double iYoff ); void setDuplicationID( unsigned int iDuplicationID = 9999 ); - void setEnergyReconstructionMethod( unsigned int iMethod ); + void setStereoReconstructionMethod( unsigned int iEnergy, unsigned int iDirection ); void setCuts( VGammaHadronCuts* iCuts ); void setContainmentProbability( double iP = 0.68, double iPError = 0.95 ) { diff --git a/inc/VInstrumentResponseFunctionData.h b/inc/VInstrumentResponseFunctionData.h index a989c133..0884b2e8 100644 --- a/inc/VInstrumentResponseFunctionData.h +++ b/inc/VInstrumentResponseFunctionData.h @@ -74,6 +74,7 @@ class VInstrumentResponseFunctionData : public TObject, public VHistogramUtiliti unsigned int fNTel; unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; // list of histograms enum E_HISTOID { E_DIFF, E_DIFF2, E_LOGDIFF, E_NIMAG, E_DIST, E_ERROR, E_RELA, @@ -106,6 +107,10 @@ class VInstrumentResponseFunctionData : public TObject, public VHistogramUtiliti { fEnergyReconstructionMethod = iMethod; } + void setDirectionReconstructionMethod( unsigned int iMethod = 0 ) + { + fDirectionReconstructionMethod = iMethod; + } void setPlottingStyle( int icolor, double iwidth = 1., int imarker = 20, double isize = 1., int iFillStyle = 0, int iLineStyle = 1 ); void setHistogramEbinning( int iN = 60, double iMin = -2.0, double iMax = 4.0 ) { diff --git a/inc/VInstrumentResponseFunctionRunParameter.h b/inc/VInstrumentResponseFunctionRunParameter.h index 5d6d2add..058fff62 100644 --- a/inc/VInstrumentResponseFunctionRunParameter.h +++ b/inc/VInstrumentResponseFunctionRunParameter.h @@ -43,6 +43,7 @@ class VInstrumentResponseFunctionRunParameter : public TNamed int fGammaHadronCutSelector; unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; unsigned int fEnergyAxisBins_log10; bool fIgnoreEnergyReconstructionQuality; unsigned int fNSpectralIndex; diff --git a/inc/VRadialAcceptance.h b/inc/VRadialAcceptance.h index a786fadb..62965b84 100644 --- a/inc/VRadialAcceptance.h +++ b/inc/VRadialAcceptance.h @@ -42,6 +42,7 @@ class VRadialAcceptance double fCut_CameraFiducialSize_max; unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; // regions excluded from background vector fXE; @@ -154,10 +155,14 @@ class VRadialAcceptance fAzCut_min = iAzMin; //!< cut on Az (shower directory) fAzCut_max = iAzMax; } - void setEnergyReconstructionMethod( unsigned int iEMethod = 1 ) + void setEnergyReconstructionMethod( unsigned int iEMethod = 0 ) { fEnergyReconstructionMethod = iEMethod; } + void setDirectionReconstructionMethod( unsigned int iEMethod = 0 ) + { + fDirectionReconstructionMethod = iEMethod; + } void setSource( double x, double y, double r, double idist, double imaxdist = 5. ); //!< set source position, radius, and minimal distance between source and background void setRegionToExcludeAcceptance( vector x, vector y, vector r ); //set the region to be exclude in the analysis // for ellipsoidal region diff --git a/inc/VTMVARunDataEnergyCut.h b/inc/VTMVARunDataEnergyCut.h index d19ef89b..5eaca2f0 100644 --- a/inc/VTMVARunDataEnergyCut.h +++ b/inc/VTMVARunDataEnergyCut.h @@ -18,13 +18,14 @@ class VTMVARunDataEnergyCut : public TNamed double fEnergyCut_Log10TeV_max; TCut fEnergyCut; unsigned int fEnergyReconstructionMethod; + unsigned int fDirectionReconstructionMethod; VTMVARunDataEnergyCut(); ~VTMVARunDataEnergyCut() {} void print(); - ClassDef( VTMVARunDataEnergyCut, 2 ); + ClassDef( VTMVARunDataEnergyCut, 3 ); }; #endif diff --git a/src/CData.cpp b/src/CData.cpp index 5f59fb32..b5f7b7a0 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -678,30 +678,68 @@ bool CData::is_GH_Gamma() * Get energy depending on analysis method * Methods: * -* 0: return friend tree result (if available) or Erec +* 0: return Erec * 1: return ErecS -* 2: return Erec -* 3: return friend tree result +* 2: return XGB energy (friend tree) +* 100: return MC energy */ float CData::get_Erec( unsigned iMethod ) { - if( iMethod == 0 && fStereoFriendTree ) + if( iMethod == 1 ) + { + return ErecS; + } + else if( iMethod == 2 ) { return Dir_Erec; } - else if( iMethod == 1 ) + else if( iMethod == 100 ) { - return ErecS; + return MCe0; + } + return Erec; +} + +/* + * Get energy chi2 + */ +float CData::get_ErecChi2( unsigned iMethod ) +{ + if( iMethod == 1 ) + { + return EChi2S; } + // not defined for XGB methods else if( iMethod == 2 ) { - return Erec; + return 0.; } - else if( iMethod == 3 ) + else if( iMethod == 100 ) { - return Dir_Erec; + return 0.; } - return Erec; + return EChi2S; +} + +/* + * Get energy dE + */ +float CData::get_ErecdE( unsigned iMethod ) +{ + if( iMethod == 1 ) + { + return dES; + } + // not defined for XGB methods + else if( iMethod == 2 ) + { + return 0.; + } + else if( iMethod == 100 ) + { + return 0.; + } + return dE; } /* @@ -709,26 +747,17 @@ float CData::get_Erec( unsigned iMethod ) * * Methods: * - * 0: return friend tree result (if available or Xoff) - * 1: return Xoff - * 2: return Xoff_intersect - * 3: return friend tree result + * 0: return Xoff + * 1: return Xoff_intersect + * 2: return friend tree result */ float CData::get_Xoff( unsigned iMethod ) { - if( iMethod == 0 && fStereoFriendTree ) - { - return Dir_Xoff; - } - else if( iMethod == 1 ) - { - return Xoff; - } - else if( iMethod == 2 ) + if( iMethod == 1 ) { return Xoff_intersect; } - else if( iMethod == 3 ) + else if( iMethod == 2 ) { return Dir_Xoff; } @@ -741,26 +770,17 @@ float CData::get_Xoff( unsigned iMethod ) * * Methods: * - * 0: return friend tree result (if available or Yoff) - * 1: return Yoff - * 2: return Yoff_intersect - * 3: return friend tree result + * 0: return Yoff + * 1: return Yoff_intersect + * 2: return friend tree result */ float CData::get_Yoff( unsigned iMethod ) { - if( iMethod == 0 && fStereoFriendTree ) - { - return Dir_Yoff; - } - else if( iMethod == 1 ) - { - return Yoff; - } - else if( iMethod == 2 ) + if( iMethod == 1 ) { return Yoff_intersect; } - else if( iMethod == 3 ) + else if( iMethod == 2 ) { return Dir_Yoff; } @@ -1048,6 +1068,8 @@ void CData::initialize_3tel_reconstruction( /* Read XGB friend tree for gamma/hadron separation and stereo reconstruction + + Note that this functions returns 0 in case the XGB file is not found. */ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name ) { @@ -1060,8 +1082,9 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name TFile *iFile = TFile::Open( file_name.c_str() ); if(!iFile || iFile->IsZombie() ) { - cout << "CData Error: cannot open XGB file " << file_name << endl; - exit( EXIT_FAILURE ); + cout << "CData warning: cannot open XGB file " << file_name << endl; + cout << "(this might be ok if no XGB analysis results are requested)" << endl; + return 0; } fXGBFiles.push_back( iFile ); diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp index 27c78674..1d93e534 100644 --- a/src/VAnaSumRunParameter.cpp +++ b/src/VAnaSumRunParameter.cpp @@ -117,6 +117,7 @@ VAnaSumRunParameter::VAnaSumRunParameter() // parameter for energy spectra (in log E) fEnergyReconstructionSpectralIndex = 2.5; fEnergyReconstructionMethod = 0; + fDirectionReconstructionMethod = 0; fEffectiveAreaVsEnergyMC = 1; // default: use effective areas vs reconstructed energy (accurate method) fEnergySpectrumBinSize = 0.05; fEnergyEffectiveAreaSmoothingIterations = -1; @@ -583,11 +584,23 @@ int VAnaSumRunParameter::readRunParameter( string i_filename ) else if( temp == "ENERGYRECONSTRUCTIONMETHOD" ) { fEnergyReconstructionMethod = ( unsigned int )atoi( temp2.c_str() ); - if( fEnergyReconstructionMethod > 1 ) + if( fEnergyReconstructionMethod > 2 ) { - cout << "Unknown parameter for ENERGYRECONSTRUCTIONMETHOD in parameter file " << i_filename << ": " << temp2 << endl; - cout << "allowed values are 0 and 1" << endl; - return 0; + cout << "Invalid value for ENERGYRECONSTRUCTIONMETHOD in parameter file " << i_filename << ": " << temp2 << endl; + cout << "allowed values are 0 (DispBDT), 1 (LT Tables), 2 (XGB stereo)" << endl; + cout << "exiting..." << endl; + exit( EXIT_FAILURE ); + } + } + else if( temp == "DIRECTIONRECONSTRUCTIONMETHOD" ) + { + fDirectionReconstructionMethod = ( unsigned int )atoi( temp2.c_str() ); + if( fDirectionReconstructionMethod > 2 ) + { + cout << "Invalid value for DIRECTIONRECONSTRUCTIONMETHOD in parameter file " << i_filename << ": " << temp2 << endl; + cout << "allowed values are 0 (DispBDT), 1 (Intersection Method), 2 (XGB stereo)" << endl; + cout << "exiting..." << endl; + exit( EXIT_FAILURE ); } } else if( temp == "DEADTIMECALCULATIONMETHOD" ) @@ -1179,15 +1192,33 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) { cout << " (use effective area A_REC)"; } - cout << ", Method " << fEnergyReconstructionMethod; + cout << endl; + cout << "Energy reconstruction method " << fEnergyReconstructionMethod; if( fEnergyReconstructionMethod == 0 ) { cout << " (dispBDT energy reconstruction)" << endl; } - else + else if( fEnergyReconstructionMethod == 1 ) { cout << " (lookup table energy reconstruction)" << endl; } + else if( fEnergyReconstructionMethod == 2 ) + { + cout << " (XGB stereo reconstruction)" << endl; + } + cout << "Direction reconstruction method " << fDirectionReconstructionMethod; + if( fDirectionReconstructionMethod == 0 ) + { + cout << " (dispBDT energy reconstruction)" << endl; + } + else if( fDirectionReconstructionMethod == 1 ) + { + cout << " (lookup table energy reconstruction)" << endl; + } + else if( fDirectionReconstructionMethod == 2 ) + { + cout << " (XGB stereo reconstruction)" << endl; + } if( fXGB_stereo_file_suffix != "" && fXGB_stereo_file_suffix != "None" ) { cout << "\t XGB stereo analysis file suffix: " << fXGB_stereo_file_suffix << endl; @@ -1480,12 +1511,14 @@ bool VAnaSumRunParameter::checkAnasumParameter( string ifile ) } else { - // check energy reconstruction method - if( iIRF->fEnergyReconstructionMethod != fEnergyReconstructionMethod ) + // check stereo reconstruction method + if( iIRF->fEnergyReconstructionMethod != fEnergyReconstructionMethod || iIRF->fDirectionReconstructionMethod != fDirectionReconstructionMethod ) { - cout << "VAnaSumRunParameter::checkAnasumParameter error in energy reconstruction method specified in runparameter file. " << endl; - cout << "\t Effective area file (" << ifile << ") uses energy reconstruction method " << iIRF->fEnergyReconstructionMethod << endl; - cout << "\t but energy reconstruction method " << fEnergyReconstructionMethod << " is requested in the anasum runparameter file. " << endl; + cout << "VAnaSumRunParameter::checkAnasumParameter error in energy or direction reconstruction method specified in runparameter file. " << endl; + cout << "\t Effective area file (" << ifile << ") uses energy methods: " << iIRF->fEnergyReconstructionMethod << " (energy) "; + cout << iIRF->fDirectionReconstructionMethod << " (direction)" << endl; + cout << "\t but reconstruction methods " << fEnergyReconstructionMethod << " (energy) and " << fDirectionReconstructionMethod; + cout << " (direction) are requested in the anasum runparameter file. " << endl; cout << "exiting..." << endl; exit( EXIT_FAILURE ); } diff --git a/src/VDataMCComparision.cpp b/src/VDataMCComparision.cpp index 27486a91..0b1938ac 100644 --- a/src/VDataMCComparision.cpp +++ b/src/VDataMCComparision.cpp @@ -21,7 +21,6 @@ VDataMCComparisionHistogramData::VDataMCComparisionHistogramData( fHis2D_size = 0; fHis2D_sizeHG = 0; fHis2D_sizeLG = 0; - } TH2D* VDataMCComparisionHistogramData::newHistogram( string iName, string iYTitle, string iXTitle, int iNbins, double ix_min, double ix_max, double iy_min, double iy_max ) @@ -141,6 +140,8 @@ VDataMCComparision::VDataMCComparision( string iname, int intel ) fCalculateMVAValues = false; fEpochATM = ""; + setStereoReconstructionMethod(); + // spectral weighting (will be set later correctly, as it is run from MC run header) fSpectralWeight = new VSpectralWeight(); // fSpectralWeight->setMCParameter( 2.0, 0.05, 20. ); @@ -183,7 +184,7 @@ void VDataMCComparision::initialGammaHadronCuts() { VGlobalRunParameter( true ); fCuts = new VGammaHadronCuts(); - fCuts->initialize(); + fCuts->initialize( fEnergyReconstructionMethod, fDirectionReconstructionMethod ); fCuts->resetCutValues(); fCuts->setNTel( 4 ); fCuts->setInstrumentEpoch( fEpochATM ); @@ -736,13 +737,13 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts // MC data if( fInput == 0 ) { - theta2 = ( fData->get_Yoff() - fData->MCyoff ) * ( fData->get_Yoff() - fData->MCyoff ) - + ( fData->get_Xoff() - fData->MCxoff ) * ( fData->get_Xoff() - fData->MCxoff ); + theta2 = ( fData->get_Yoff( fDirectionReconstructionMethod ) - fData->MCyoff ) * ( fData->get_Yoff( fDirectionReconstructionMethod ) - fData->MCyoff ) + + ( fData->get_Xoff( fDirectionReconstructionMethod ) - fData->MCxoff ) * ( fData->get_Xoff( fDirectionReconstructionMethod ) - fData->MCxoff ); } // data runs (on and off runs) else { - pair< float, float > xyoff_derot = fData->get_XYoff_derot(); + pair< float, float > xyoff_derot = fData->get_XYoff_derot( fDirectionReconstructionMethod ); float on_x = xyoff_derot.first + fWobbleEast; float on_y = xyoff_derot.second + fWobbleNorth; // off data @@ -959,7 +960,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts fCuts->newEvent(); // apply stereo quality cuts to ensure good response of TMV if( fCuts->applyInsideFiducialAreaCut( true ) - && fCuts->applyStereoQualityCuts( 1, true, i, true ) ) + && fCuts->applyStereoQualityCuts( true, i, true ) ) { fCuts->applyTMVACut( i ); if( fCuts->getTMVA_EvaluationResult() > -90. && fHistoArray[EMVA] ) diff --git a/src/VEffectiveAreaCalculator.cpp b/src/VEffectiveAreaCalculator.cpp index 30c430c6..1c02ac3f 100644 --- a/src/VEffectiveAreaCalculator.cpp +++ b/src/VEffectiveAreaCalculator.cpp @@ -56,7 +56,6 @@ VEffectiveAreaCalculator::VEffectiveAreaCalculator( VInstrumentResponseFunctionR fNoise.push_back( 0 ); fPedVar.push_back( 0. ); } - setIgnoreEnergyReconstructionCuts( fRunPara->fIgnoreEnergyReconstructionQuality ); setIsotropicArrivalDirections( fRunPara->fIsotropicArrivalDirections ); setWobbleOffset( fRunPara->fXoff, fRunPara->fYoff ); setNoiseLevel( fRunPara->fnoise, fRunPara->fpedvar ); @@ -2017,7 +2016,9 @@ bool VEffectiveAreaCalculator::getMonteCarloSpectra( VEffectiveAreaCalculatorMCH * */ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, - VEffectiveAreaCalculatorMCHistograms* iMC_histo, unsigned int iMethod ) + VEffectiveAreaCalculatorMCHistograms* iMC_histo, + unsigned int iEnergyReconstructionMethod, + unsigned int iDirectionReconstructionMethod ) { bool bDebugCuts = false; // lots of debug output @@ -2030,12 +2031,6 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, } resetHistogramsVectors( ize ); - // do not require successful energy reconstruction - if( fIgnoreEnergyReconstruction ) - { - iMethod = 100; - } - ////////////////////////////////////////////////////////////////// // total Monte Carlo core scatter area (depends on CORSIKA shower core scatter mode) fMC_ScatterArea = 0.; @@ -2066,7 +2061,8 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, cout << "\t noise level: " << fNoise[ize] << " (pedvar: " << fPedVar[ize] << ")" << endl; cout << "\t area (" << fScatterMode[ize] << ") [m^2]: " << fMC_ScatterArea; cout << " (scatter radius " << fAreaRadius[ize] << " [m])" << endl; - cout << "\t energy reconstruction method: " << iMethod << endl; + cout << "\t energy reconstruction method: " << iEnergyReconstructionMethod << endl; + cout << "\t direction reconstruction method: " << iDirectionReconstructionMethod << endl; if( fIsotropicArrivalDirections ) { cout << "\t assuming isotropic arrival directions" << endl; @@ -2130,7 +2126,7 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, { d->GetEntry( i ); - if( d->get_Xoff() < -999. || d->get_Yoff() < -999. ) continue; + if( d->get_Xoff( iDirectionReconstructionMethod ) < -999. || d->get_Yoff( iDirectionReconstructionMethod ) < -999. ) continue; // update cut statistics fCuts->newEvent(); @@ -2164,7 +2160,7 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, // apply reconstruction cuts // apply reconstruction quality cuts - if(!fCuts->applyStereoQualityCuts( iMethod, true, i, true ) ) + if(!fCuts->applyStereoQualityCuts( true, i, true ) ) { continue; } @@ -2175,7 +2171,7 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, { cout << "#1 CUT applyInsideFiducialAreaCut "; cout << fCuts->applyInsideFiducialAreaCut(); - cout << "\t" << fCuts->applyStereoQualityCuts( iMethod, false, i, true ) << endl; + cout << "\t" << fCuts->applyStereoQualityCuts( false, i, true ) << endl; } if(!fCuts->applyInsideFiducialAreaCut( true ) ) { @@ -2191,7 +2187,7 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, bool bDirectionCut = false; if(!fIsotropicArrivalDirections ) { - if(!fCuts->applyDirectionCuts( iMethod, true ) ) + if(!fCuts->applyDirectionCuts( true ) ) { bDirectionCut = true; } @@ -2200,7 +2196,7 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, // (command line option -d) else { - if(!fCuts->applyDirectionCuts( iMethod, true, 0., 0. ) ) + if(!fCuts->applyDirectionCuts( true, 0., 0. ) ) { bDirectionCut = true; } @@ -2212,16 +2208,13 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, ////////////////////////////////////// // apply energy reconstruction quality cut - if(!fIgnoreEnergyReconstruction ) + if( bDebugCuts ) { - if( bDebugCuts ) - { - cout << "#4 EnergyReconstructionQualityCuts " << fCuts->applyEnergyReconstructionQualityCuts( iMethod ) << endl; - } - if(!fCuts->applyEnergyReconstructionQualityCuts( iMethod, true ) ) - { - continue; - } + cout << "#4 EnergyReconstructionQualityCuts " << fCuts->applyEnergyReconstructionQualityCuts() << endl; + } + if(!fCuts->applyEnergyReconstructionQualityCuts( true ) ) + { + continue; } if(!bDirectionCut ) { @@ -2229,23 +2222,14 @@ bool VEffectiveAreaCalculator::fill( TH1D* hE0mc, CData* d, } // skip event if no energy has been reconstructed - // get energy according to which reconstruction method - if( fIgnoreEnergyReconstruction ) + eRecLin = d->get_Erec( iEnergyReconstructionMethod ); + if( eRecLin > 0. ) { - eRec = log10( d->MCe0 ); - eRecLin = d->MCe0; + eRec = log10( eRecLin ); } else { - eRecLin = d->get_Erec( iMethod ); - if( eRecLin > 0. ) - { - eRec = log10( eRecLin ); - } - else - { - continue; - } + continue; } /////////////////////////////////////////// diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index df5a14f8..85cbf19d 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -75,11 +75,13 @@ VGammaHadronCuts::VGammaHadronCuts() fTMVA_EvaluationResult = -99.; fTMVAEvaluatorResults = 0; + setStereoReconstructionMethod(); setArrayCentre(); } -void VGammaHadronCuts::initialize() +void VGammaHadronCuts::initialize( unsigned int iEnergyMethod, unsigned int iDirectionMethod ) { + setStereoReconstructionMethod( iEnergyMethod, iDirectionMethod ); fStats = new VGammaHadronCutsStatistics(); fStats->initialize(); } @@ -763,7 +765,7 @@ void VGammaHadronCuts::printCutSummary() /* ensure event quality, reasonable output from the table variables, etc */ -bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructionMethod, bool bCount, int iEntry, bool fIsOn ) +bool VGammaHadronCuts::applyStereoQualityCuts( bool bCount, int iEntry, bool fIsOn ) { ///////////////////////////////////////////////////////////////////////////////// // apply number of images cut @@ -801,7 +803,7 @@ bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructio } ///////////////////////////////////////////////////////////////////////////////// - if( iEnergyReconstructionMethod != 99 ) + if( fEnergyReconstructionMethod != 99 ) { // quality cut for MSCW/L reconstruction cuts if( fGammaHadronCutSelector % 10 < 1 && ( fData->MSCW < -50. || fData->MSCL < -50. ) ) @@ -826,8 +828,8 @@ bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructio ///////////////////////////////////////////////////////////////////////////////// // check energy reconstruction - double iErec = getReconstructedEnergy( iEnergyReconstructionMethod ); - double iEchi2 = getReconstructedEnergyChi2( iEnergyReconstructionMethod ); + double iErec = fData->get_Erec( fEnergyReconstructionMethod ); + double iEchi2 = fData->get_ErecChi2( fEnergyReconstructionMethod ); if( iErec > 0. && iEchi2 <= fCut_EChi2_min ) { @@ -981,16 +983,15 @@ bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructio iEnergyReconstructionMethod == 100: return always true */ -bool VGammaHadronCuts::applyEnergyReconstructionQualityCuts( unsigned int iEnergyReconstructionMethod, bool bCount ) +bool VGammaHadronCuts::applyEnergyReconstructionQualityCuts( bool bCount ) { - double iErec = getReconstructedEnergy( iEnergyReconstructionMethod ); - double iErecChi2 = getReconstructedEnergyChi2( iEnergyReconstructionMethod ); - double idE = getReconstructedEnergydE( iEnergyReconstructionMethod ); - - if( iEnergyReconstructionMethod == 100 ) + if( fEnergyReconstructionMethod == 100 ) { return true; } + double iErec = fData->get_Erec( fEnergyReconstructionMethod ); + double iErecChi2 = fData->get_ErecChi2( fEnergyReconstructionMethod ); + double idE = fData->get_ErecdE( fEnergyReconstructionMethod ); if( iErecChi2 < fCut_EChi2_min || iErecChi2 > fCut_EChi2_max ) { @@ -1454,7 +1455,7 @@ bool VGammaHadronCuts::initPhaseCuts( string iDir ) bool VGammaHadronCuts::applyInsideFiducialAreaCut( bool bCount ) { - return applyInsideFiducialAreaCut( getReconstructedXoff(), getReconstructedYoff(), bCount ); + return applyInsideFiducialAreaCut( fData->get_Xoff( fDirectionReconstructionMethod ), fData->get_Yoff( fDirectionReconstructionMethod ), bCount ); } bool VGammaHadronCuts::applyInsideFiducialAreaCut( float Xoff, float Yoff, bool bCount ) @@ -1533,7 +1534,7 @@ bool VGammaHadronCuts::applyMCXYoffCut( double xoff, double yoff, bool bCount ) x0, y0: calculate theta2 relative to these points (-99999. if relative to MCx/yoff) */ -bool VGammaHadronCuts::applyDirectionCuts( unsigned int fEnergyReconstructionMethod, bool bCount, double x0, double y0 ) +bool VGammaHadronCuts::applyDirectionCuts( bool bCount, double x0, double y0 ) { double theta2 = 0.; @@ -1563,7 +1564,8 @@ bool VGammaHadronCuts::applyDirectionCuts( unsigned int fEnergyReconstructionMet } // calculate theta2 - theta2 = ( getReconstructedXoff() - x0 ) * ( getReconstructedXoff() - x0 ) + ( getReconstructedYoff() - y0 ) * ( getReconstructedYoff() - y0 ); + theta2 = ( fData->get_Xoff( fDirectionReconstructionMethod ) - x0 ) * ( fData->get_Xoff( fDirectionReconstructionMethod ) - x0 ) + + ( fData->get_Yoff( fDirectionReconstructionMethod ) - y0 ) * ( fData->get_Yoff( fDirectionReconstructionMethod ) - y0 ); // fetch theta2 cut (max) double i_theta2_cut_max = getTheta2Cut_max(); @@ -1689,82 +1691,3 @@ string VGammaHadronCuts::getTelToAnalyzeString() return iTemp.str(); } - -double VGammaHadronCuts::getReconstructedEnergy( unsigned int iEnergyReconstructionMethod ) -{ - if(!fData ) - { - return -99.; - } - return fData->get_Erec( iEnergyReconstructionMethod ); -} - -double VGammaHadronCuts::getReconstructedEnergyChi2( unsigned int iEnergyReconstructionMethod ) -{ - if(!fData ) - { - return -99.; - } - if( iEnergyReconstructionMethod == 0 ) - { - return fData->EChi2; - } - else if( iEnergyReconstructionMethod == 1 ) - { - return fData->EChi2S; - } - return -99.; -} - -double VGammaHadronCuts::getReconstructedEnergydE( unsigned int iEnergyReconstructionMethod ) -{ - if(!fData ) - { - return -99.; - } - if( iEnergyReconstructionMethod == 0 ) - { - return fData->dE; - } - else if( iEnergyReconstructionMethod == 1 ) - { - return fData->dES; - } - return -99.; -} - -double VGammaHadronCuts::getReconstructedXoff() -{ - if(!fData ) - { - return -9999.; - } - return fData->get_Xoff( 0 ); -} - -double VGammaHadronCuts::getReconstructedYoff() -{ - if(!fData ) - { - return -9999.; - } - return fData->get_Yoff( 0 ); -} - -double VGammaHadronCuts::getReconstructedXcore() -{ - if(!fData ) - { - return -9999.; - } - return fData->Xcore; -} - -double VGammaHadronCuts::getReconstructedYcore() -{ - if(!fData ) - { - return -9999.; - } - return fData->Ycore; -} diff --git a/src/VInstrumentResponseFunction.cpp b/src/VInstrumentResponseFunction.cpp index 1410dd58..676f15fd 100644 --- a/src/VInstrumentResponseFunction.cpp +++ b/src/VInstrumentResponseFunction.cpp @@ -15,7 +15,6 @@ VInstrumentResponseFunction::VInstrumentResponseFunction() fData = 0; fAnaCuts = 0; - fEnergyReconstructionMethod = 0; fSpectralWeight = new VSpectralWeight(); @@ -33,8 +32,7 @@ void VInstrumentResponseFunction::setRunParameter( VInstrumentResponseFunctionRu { return; } - fEnergyReconstructionMethod = iRunPara->fEnergyReconstructionMethod; - setEnergyReconstructionMethod( iRunPara->fEnergyReconstructionMethod ); + setStereoReconstructionMethod( iRunPara->fEnergyReconstructionMethod, iRunPara->fDirectionReconstructionMethod ); setMonteCarloEnergyRange( iRunPara->fMCEnergy_min, iRunPara->fMCEnergy_max, TMath::Abs( iRunPara->fMCEnergy_index ) ); fVMinAz = iRunPara->fAzMin; @@ -96,7 +94,8 @@ bool VInstrumentResponseFunction::initialize( string iName, string iType, unsign } i_irf.back()->setData( iZe, ( int )j, fVMinAz[j], fVMaxAz[j], iNoise, iPedvars, fVSpectralIndex[i], iXoff, iYoff ); - i_irf.back()->setEnergyReconstructionMethod( fEnergyReconstructionMethod ); + i_irf.back()->setEnergyReconstructionMethod( fRunPara->fEnergyReconstructionMethod ); + i_irf.back()->setDirectionReconstructionMethod( fRunPara->fDirectionReconstructionMethod ); } fIRFData.push_back( i_irf ); } @@ -174,7 +173,7 @@ bool VInstrumentResponseFunction::fillEventData() } // apply reconstruction quality cuts - if(!fAnaCuts->applyStereoQualityCuts( fEnergyReconstructionMethod, true, i, true ) ) + if(!fAnaCuts->applyStereoQualityCuts( true, i, true ) ) { continue; } @@ -299,7 +298,7 @@ void VInstrumentResponseFunction::setDataTree( CData* iData ) } -void VInstrumentResponseFunction::setEnergyReconstructionMethod( unsigned int iMethod ) +void VInstrumentResponseFunction::setStereoReconstructionMethod( unsigned int iEnergy, unsigned int iDirection ) { for( unsigned int i = 0; i < fIRFData.size(); i++ ) { @@ -307,7 +306,8 @@ void VInstrumentResponseFunction::setEnergyReconstructionMethod( unsigned int iM { if( fIRFData[i][j] ) { - fIRFData[i][j]->setEnergyReconstructionMethod( iMethod ); + fIRFData[i][j]->setEnergyReconstructionMethod( iEnergy ); + fIRFData[i][j]->setDirectionReconstructionMethod( iDirection ); } } } diff --git a/src/VInstrumentResponseFunctionData.cpp b/src/VInstrumentResponseFunctionData.cpp index 25fca1e9..b7d87c74 100644 --- a/src/VInstrumentResponseFunctionData.cpp +++ b/src/VInstrumentResponseFunctionData.cpp @@ -440,7 +440,9 @@ void VInstrumentResponseFunctionData::fill( double iWeight ) bool bPlotResolution_vs_reconstructedEnergy = true; // simple quality check (FOV shouldn't be larger than 50 deg) - if( fData->get_Xoff() < -50. || fData->get_Yoff() < -50. ) + double xoff = fData->get_Xoff( fDirectionReconstructionMethod ); + double yoff = fData->get_Yoff( fDirectionReconstructionMethod ); + if( xoff < -50. || yoff < -50. ) { return; } @@ -460,11 +462,9 @@ void VInstrumentResponseFunctionData::fill( double iWeight ) if( fType_numeric == 0 ) { // angular difference - iDiff = sqrt(( fData->get_Xoff() - fData->MCxoff ) * ( fData->get_Xoff() - fData->MCxoff ) + - ( fData->get_Yoff() - fData->MCyoff ) * ( fData->get_Yoff() - fData->MCyoff ) ); + iDiff = sqrt(( xoff - fData->MCxoff ) * ( xoff - fData->MCxoff ) + ( yoff - fData->MCyoff ) * ( yoff - fData->MCyoff ) ); // error - iError = sqrt( fData->get_Xoff() * fData->get_Xoff() + fData->get_Yoff() * fData->get_Yoff() ) - - sqrt( fData->MCxoff* fData->MCxoff + fData->MCyoff* fData->MCyoff ); + iError = sqrt( xoff* xoff + yoff* yoff ) - sqrt( fData->MCxoff* fData->MCxoff + fData->MCyoff* fData->MCyoff ); // relative error (not sure if it is useful) iErrorRelative = -99.e6; } diff --git a/src/VInstrumentResponseFunctionRunParameter.cpp b/src/VInstrumentResponseFunctionRunParameter.cpp index 505a105d..ffb50628 100644 --- a/src/VInstrumentResponseFunctionRunParameter.cpp +++ b/src/VInstrumentResponseFunctionRunParameter.cpp @@ -17,6 +17,7 @@ VInstrumentResponseFunctionRunParameter::VInstrumentResponseFunctionRunParameter fSpectralIndexStep = 0.1; fEnergyReconstructionMethod = 0; + fDirectionReconstructionMethod = 0; fEnergyAxisBins_log10 = 60; fIgnoreEnergyReconstructionQuality = false; @@ -158,6 +159,14 @@ bool VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile( stri is_stream >> fEnergyReconstructionMethod; } } + // direction reconstruction method + else if( temp == "DIRECTIONRECONSTRUCTIONMETHOD" ) + { + if(!( is_stream >> std::ws ).eof() ) + { + is_stream >> fDirectionReconstructionMethod; + } + } // number of bins on log10 energy axis else if( temp == "ENERGYAXISBINS" ) { @@ -653,9 +662,11 @@ void VInstrumentResponseFunctionRunParameter::print() if( fIgnoreEnergyReconstructionQuality ) { cout << ", ignoring cut on quality of energy reconstruction"; + fEnergyReconstructionMethod = 100; } cout << endl; cout << "energy reconstruction method " << fEnergyReconstructionMethod << endl; + cout << "direction reconstruction method " << fDirectionReconstructionMethod << endl; cout << endl; cout << "input Monte Carlo with following parameters (will be modified later): " << endl; diff --git a/src/VRadialAcceptance.cpp b/src/VRadialAcceptance.cpp index 17f3c16a..7b1fe048 100644 --- a/src/VRadialAcceptance.cpp +++ b/src/VRadialAcceptance.cpp @@ -584,7 +584,7 @@ int VRadialAcceptance::fillAcceptanceFromData( CData* iData, int entry ) bool bPassed = false; // apply some basic quality cuts - if( fCuts->applyInsideFiducialAreaCut() && fCuts->applyStereoQualityCuts( fEnergyReconstructionMethod, false, entry, true ) ) + if( fCuts->applyInsideFiducialAreaCut() && fCuts->applyStereoQualityCuts( false, entry, true ) ) { // gamma/hadron cuts if(!fCuts->isGamma( entry, false ) ) diff --git a/src/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index 7507fd34..6d0705ae 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -109,6 +109,7 @@ VStereoAnalysis::VStereoAnalysis( bool ion, string i_hsuffix, VAnaSumRunParamete // define the cuts fCuts = new VGammaHadronCuts(); + fCuts->setStereoReconstructionMethod( fRunPara->fEnergyReconstructionMethod, fRunPara->fDirectionReconstructionMethod ); char hname[200]; if( fIsOn ) { @@ -478,11 +479,11 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d } // get energy (depending on energy reconstruction method) - iErec = fCuts->getReconstructedEnergy( fRunPara->fEnergyReconstructionMethod ); - iErecChi2 = fCuts->getReconstructedEnergyChi2( fRunPara->fEnergyReconstructionMethod ); + iErec = fDataRun->get_Erec( fRunPara->fEnergyReconstructionMethod ); + iErecChi2 = fDataRun->get_ErecChi2( fRunPara->fEnergyReconstructionMethod ); // get shower direction (depending on shower reconstruction method) - iXoff = fCuts->getReconstructedXoff(); - iYoff = fCuts->getReconstructedYoff(); + iXoff = fDataRun->get_Xoff( fRunPara->fDirectionReconstructionMethod ); + iYoff = fDataRun->get_Yoff( fRunPara->fDirectionReconstructionMethod ); //////////////////////////////////////////////// // apply all quality cuts @@ -494,7 +495,7 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d } // stereo quality cuts (e.g. successful direction, mscw, mscl reconstruction) - if(!fCuts->applyStereoQualityCuts( fRunPara->fEnergyReconstructionMethod, false, i, fIsOn ) ) + if(!fCuts->applyStereoQualityCuts( false, i, fIsOn ) ) { continue; } @@ -517,7 +518,7 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d fDataRun->Ze, iErec, fDataRun->runNumber, bIsGamma, i_theta2 ); // energy reconstruction cut - bEnergyQualityCuts = fCuts->applyEnergyReconstructionQualityCuts( fRunPara->fEnergyReconstructionMethod ); + bEnergyQualityCuts = fCuts->applyEnergyReconstructionQualityCuts(); ///////////////////////////////////////////////////////////////////////////////////////////////////////// // following histograms (theta2, mscw, mscl, core position, etc.) assume source at given target position @@ -618,7 +619,7 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d fHisto[fHisCounter]->hTriggerPatternAfterCuts->Fill( fDataRun->LTrig ); fHisto[fHisCounter]->hImagePatternAfterCuts->Fill( fDataRun->ImgSel ); // make core plots - fHisto[fHisCounter]->hcore->Fill( fCuts->getReconstructedXcore(), fCuts->getReconstructedYcore() ); + fHisto[fHisCounter]->hcore->Fill( fDataRun->Xcore, fDataRun->Ycore ); // ################################## // spectral energy reconstruction // apply energy reconstruction cuts @@ -2111,7 +2112,7 @@ void VStereoAnalysis::fill_TreeWithSelectedEvents( CData* c, double i_xderot, do fTreeSelected_theta2 = i_theta2; fTreeSelected_Xoff = c->Xoff; fTreeSelected_Yoff = c->Yoff; - pair tmp_xy_derot = c->get_XYoff_derot(); + pair tmp_xy_derot = c->get_XYoff_derot( fRunPara->fDirectionReconstructionMethod ); fTreeSelected_Xoff_derot = tmp_xy_derot.first; fTreeSelected_Yoff_derot = tmp_xy_derot.second; fTreeSelected_Xcore = c->Xcore; @@ -2122,7 +2123,7 @@ void VStereoAnalysis::fill_TreeWithSelectedEvents( CData* c, double i_xderot, do fTreeSelected_MLR = c->MLR; fTreeSelected_Erec = c->Erec; fTreeSelected_EChi2 = c->EChi2; - fTreeSelected_ErecS = c->get_Erec(); + fTreeSelected_ErecS = c->get_Erec( fRunPara->fEnergyReconstructionMethod ); fTreeSelected_EChi2S = c->EChi2S; fTreeSelected_EmissionHeight = c->EmissionHeight; fTreeSelected_EmissionHeightChi2 = c->EmissionHeightChi2; @@ -2271,8 +2272,8 @@ void VStereoAnalysis::fill_DL3Tree( CData* c, double i_xderot, double i_yderot, fDL3EventTree_Yderot = i_yderot; // Derotated Gamma Point-Of-Origin (deg, DEC) if( fCuts ) { - fDL3EventTree_Erec = fCuts->getReconstructedEnergy( fRunPara->fEnergyReconstructionMethod ); - fDL3EventTree_Erec_Err = fCuts->getReconstructedEnergydE( fRunPara->fEnergyReconstructionMethod ); // Reconstructed Gamma Energy (TeV) Error + fDL3EventTree_Erec = c->get_Erec( fRunPara->fEnergyReconstructionMethod ); + fDL3EventTree_Erec_Err = c->get_ErecdE( fRunPara->fEnergyReconstructionMethod ); // Reconstructed Gamma Energy (TeV) Error } else { diff --git a/src/compareDatawithMC.cpp b/src/compareDatawithMC.cpp index bb2d371e..add69fe9 100644 --- a/src/compareDatawithMC.cpp +++ b/src/compareDatawithMC.cpp @@ -247,6 +247,11 @@ int main( int argc, char* argv[] ) { fEpochATM = argv[5]; } + unsigned int fEnergyReconstructionMethod = 0; + unsigned int fDirectionReconstructionMethod = 0; + cout << "Stereo reconstruction methods for energy: " << fEnergyReconstructionMethod; + cout << ", direction: " << fDirectionReconstructionMethod << endl; + // test number of telescopes int iNT = 0; @@ -309,6 +314,7 @@ int main( int argc, char* argv[] ) cout << fInputData[i].fType << endl; cout << "----" << endl; fStereoCompare.push_back( new VDataMCComparision( fInputData[i].fType, fInputData[i].fNTelescopes ) ); + fStereoCompare.back()->setStereoReconstructionMethod( fEnergyReconstructionMethod, fDirectionReconstructionMethod ); if( fCalculateMVACut ) { fStereoCompare.back()->setTMVABDTComparision( fEpochATM ); @@ -354,6 +360,7 @@ int main( int argc, char* argv[] ) cout << "DIFF" << endl; cout << "----" << endl; VDataMCComparision* fDiff = new VDataMCComparision( "DIFF", iNT ); + fDiff->setStereoReconstructionMethod( fEnergyReconstructionMethod, fDirectionReconstructionMethod ); // assume 5 background regions fDiff->setOnOffHistograms( fStereoCompareOn, fStereoCompareOff, 1. / 5. ); fDiff->writeHistograms( fOutputfile ); diff --git a/src/makeEffectiveArea.cpp b/src/makeEffectiveArea.cpp index 743ed265..5d1cede1 100644 --- a/src/makeEffectiveArea.cpp +++ b/src/makeEffectiveArea.cpp @@ -101,7 +101,7 @@ int main( int argc, char* argv[] ) ///////////////////////////////////////////////////////////////// // gamma/hadron cuts VGammaHadronCuts* fCuts = new VGammaHadronCuts(); - fCuts->initialize(); + fCuts->initialize( fRunPara->fEnergyReconstructionMethod, fRunPara->fDirectionReconstructionMethod ); fCuts->setNTel( fRunPara->telconfig_ntel, fRunPara->telconfig_arraycentre_X, fRunPara->telconfig_arraycentre_Y ); fCuts->setInstrumentEpoch( fRunPara->getInstrumentATMString() ); fCuts->setTelToAnalyze( fRunPara->fTelToAnalyse ); @@ -315,7 +315,7 @@ int main( int argc, char* argv[] ) } } - fEffectiveAreaCalculator.fill( hE0mc, &d, fMC_histo, fRunPara->fEnergyReconstructionMethod ); + fEffectiveAreaCalculator.fill( hE0mc, &d, fMC_histo, fRunPara->fEnergyReconstructionMethod, fRunPara->fDirectionReconstructionMethod ); fStopWatch.Print(); } diff --git a/src/makeRadialAcceptance.cpp b/src/makeRadialAcceptance.cpp index 54eeea74..66cd39ec 100644 --- a/src/makeRadialAcceptance.cpp +++ b/src/makeRadialAcceptance.cpp @@ -100,7 +100,7 @@ int main( int argc, char* argv[] ) // read gamma/hadron cuts from cut file VGammaHadronCuts* fCuts = new VGammaHadronCuts(); - fCuts->initialize(); + fCuts->initialize( fRunPara->fEnergyReconstructionMethod, fRunPara->fDirectionReconstructionMethod ); fCuts->setNTel( ntel ); fCuts->setTelToAnalyze( teltoana ); if(!fCuts->readCuts( cutfilename, 2 ) ) From 0416d329a4554f0ff5eef79796921d7ddf624f56 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 16:07:29 +0100 Subject: [PATCH 15/21] changelog --- docs/changes/339.feature.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 docs/changes/339.feature.md diff --git a/docs/changes/339.feature.md b/docs/changes/339.feature.md new file mode 100644 index 00000000..d46c4c26 --- /dev/null +++ b/docs/changes/339.feature.md @@ -0,0 +1,12 @@ +Allow to choose energy and direction reconstruction methods: + +``` +Stereo reconstruction +--------------------- +Method ids: 0 (DispBDT), 1 (LT Tables), 2 (XGB stereo) +* ENERGYRECONSTRUCTIONMETHOD 0 +Method ids: 0 (DispBDT), 1 (Intersection Method), 2 (XGB stereo) +* DIRECTIONRECONSTRUCTIONMETHOD 0 +``` + +Consistently use `CData` methods to get variables. From 81df5d65794cfdc57565c0097eeb815f361cb61c Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 16:28:07 +0100 Subject: [PATCH 16/21] require method statement --- inc/CData.h | 12 ++--- inc/VTMVAEvaluator.h | 5 +- src/VDataMCComparision.cpp | 96 ++++++++++++++++++++------------------ src/VGammaHadronCuts.cpp | 6 +-- src/VRadialAcceptance.cpp | 17 ++++--- src/VTMVAEvaluator.cpp | 15 +++--- 6 files changed, 80 insertions(+), 71 deletions(-) diff --git a/inc/CData.h b/inc/CData.h index 75944dcd..2c2101b5 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -273,14 +273,14 @@ class CData CData( TTree* tree, bool bMC, bool bShort, string file_name, string stereo_suffix, string gamma_hadron_suffix ); virtual ~CData(); virtual Int_t GetEntry( Long64_t entry ); - float get_Erec( unsigned int method = 0 ); - float get_ErecChi2( unsigned int method = 0 ); - float get_ErecdE( unsigned int method = 0 ); - float get_Xoff( unsigned int method = 0 ); - float get_Yoff( unsigned int method = 0 ); + float get_Erec( unsigned int method ); + float get_ErecChi2( unsigned int method ); + float get_ErecdE( unsigned int method ); + float get_Xoff( unsigned int method ); + float get_Yoff( unsigned int method ); void initialize_xgb_tree( TTree* stereoTree, TTree* ghTree ); TTree* getXGBTree( string file_name, string suffix, string tree_name ); - pair get_XYoff_derot( unsigned int method = 0 ); + pair get_XYoff_derot( unsigned int method ); virtual Long64_t LoadTree( Long64_t entry ); float get_GH_Gamma_Prediction(); bool is_GH_Gamma(); diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h index 1fdde49f..1246ebe5 100644 --- a/inc/VTMVAEvaluator.h +++ b/inc/VTMVAEvaluator.h @@ -101,6 +101,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities vector< VTMVAEvaluatorData* > fTMVAData; CData* fData; + unsigned int fEnergyReconstructionMethod; map< unsigned int, double > fSignalEfficiencyMap; // from user: energy dependent signal efficiency double fSignalEfficiencyNoVec; @@ -171,7 +172,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities public: - VTMVAEvaluator(); + VTMVAEvaluator(unsigned int iEnergyReconstructionMethod=0); ~VTMVAEvaluator() {}; bool evaluate( bool interpolate_mva = false, bool use_average_zenith_angle = true ); @@ -241,7 +242,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities } void setTMVAMethod( string iMethodName = "BDT" ); - ClassDef( VTMVAEvaluator, 37 ); + ClassDef( VTMVAEvaluator, 38 ); }; #endif diff --git a/src/VDataMCComparision.cpp b/src/VDataMCComparision.cpp index 0b1938ac..952dbc5d 100644 --- a/src/VDataMCComparision.cpp +++ b/src/VDataMCComparision.cpp @@ -617,6 +617,8 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts double rdist1 = 0.; double theta2 = 0.; + double erec = 0.; + double erec_log10 = 0.; int iOldRun = 0; @@ -692,6 +694,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts ///////////////////////////////////////////////// ///////////////////////////////////////////////// // quality cuts + erec = fData->get_Erec(fEnergyReconstructionMethod); // nimage cut if( fData->NImages < fNImages_min ) @@ -699,10 +702,11 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts continue; } // successful energy reconstruction - if( fData->EChi2S < 0 || fData->get_Erec() <= 0 ) + if( fData->EChi2S < 0 || erec <= 0 ) { continue; } + erec_log10 = log10(erec); // successful lookup table reconstruction if( fData->MSCW < -50. || fData->MSCL < -50. ) { @@ -862,11 +866,11 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { if( fHistoArray[ETHETA2] ) { - fHistoArray[ETHETA2]->fill( theta2, weight, log10( fData->get_Erec() ) ); + fHistoArray[ETHETA2]->fill( theta2, weight, erec_log10 ); } if( fHistoArray[ELTHETA2] ) { - fHistoArray[ELTHETA2]->fill( log10( theta2 ), weight, log10( fData->get_Erec() ) ); + fHistoArray[ELTHETA2]->fill( log10( theta2 ), weight, erec_log10 ); } hYt2->Fill( log10( theta2 ), fData->Ycore, weight ); } @@ -874,7 +878,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { if( fHistoArray[EMSCW] ) { - fHistoArray[EMSCW]->fill( fData->MSCW, weight, log10( fData->get_Erec() ) ); + fHistoArray[EMSCW]->fill( fData->MSCW, weight, erec_log10 ); } for( int j = 0; j < fData->NImages; j++ ) { @@ -885,7 +889,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts && iTelImage < fHistoSingleTel[EMWRT].size() && fHistoSingleTel[EMWRT][iTelImage] ) { fHistoSingleTel[EMWRT][iTelImage]->fill( fData->width[iTelImage] / fData->MSCWT[iTelImage], - weight, log10( fData->get_Erec() ), + weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } @@ -893,14 +897,14 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts && fHistoSingleTel[EMSCWT][iTelImage] ) { fHistoSingleTel[EMSCWT][iTelImage]->fill( fData->MSCWT[iTelImage], - weight, log10( fData->get_Erec() ), + weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } } if( fHistoArray[EMWR] ) { - fHistoArray[EMWR]->fill( fData->MWR, weight, log10( fData->get_Erec() ) ); + fHistoArray[EMWR]->fill( fData->MWR, weight, erec_log10 ); } } if( fSingleTelescopeCuts != -1 @@ -908,7 +912,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { if( fHistoArray[EMSCL] ) { - fHistoArray[EMSCL]->fill( fData->MSCL, weight, log10( fData->get_Erec() ) ); + fHistoArray[EMSCL]->fill( fData->MSCL, weight, erec_log10 ); } // Again: this is extremely bad as some index have a range of // NTel and others a range of NImages @@ -921,38 +925,38 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts && iTelImage < fHistoSingleTel[EMLTT].size() && fHistoSingleTel[EMLTT][iTelImage] ) { fHistoSingleTel[EMLTT][iTelImage]->fill( fData->length[iTelImage] / fData->MSCLT[iTelImage], - weight, log10( fData->get_Erec() ), + weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EMSCLT].size() > 0 && iTelImage < fHistoSingleTel[EMSCLT].size() && fHistoSingleTel[EMSCLT][iTelImage] ) { - fHistoSingleTel[EMSCLT][iTelImage]->fill( fData->MSCLT[iTelImage], weight, log10( fData->get_Erec() ), + fHistoSingleTel[EMSCLT][iTelImage]->fill( fData->MSCLT[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } } if( fHistoArray[EMLR] ) { - fHistoArray[EMLR]->fill( fData->MLR, weight, log10( fData->get_Erec() ) ); + fHistoArray[EMLR]->fill( fData->MLR, weight, erec_log10 ); } if( fData->EmissionHeight > 0. && fHistoArray[EEMISSIONHEIGHT] ) { fHistoArray[EEMISSIONHEIGHT]->fill( getCorrectedEmissionHeight( fData->EmissionHeight, fData->Ze ), - weight, log10( fData->get_Erec() ) ); + weight, erec_log10 ); } if( fHistoArray[ENIMAGES] ) { - fHistoArray[ENIMAGES]->fill( fData->NImages, weight, log10( fData->get_Erec() ) ); + fHistoArray[ENIMAGES]->fill( fData->NImages, weight, erec_log10 ); } if( fHistoArray[EPEDVAR] ) { - fHistoArray[EPEDVAR]->fill( fData->meanPedvar_Image, weight, log10( fData->get_Erec() ) ); + fHistoArray[EPEDVAR]->fill( fData->meanPedvar_Image, weight, erec_log10 ); } if( fHistoArray[EIMGSEL] ) { - fHistoArray[EIMGSEL]->fill( fData->ImgSel, weight, log10( fData->get_Erec() ) ); + fHistoArray[EIMGSEL]->fill( fData->ImgSel, weight, erec_log10 ); } if( fCuts ) @@ -966,7 +970,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts if( fCuts->getTMVA_EvaluationResult() > -90. && fHistoArray[EMVA] ) { fHistoArray[EMVA]->fill( - fCuts->getTMVA_EvaluationResult(), weight, log10( fData->get_Erec() ) ); + fCuts->getTMVA_EvaluationResult(), weight, erec_log10 ); } } } @@ -980,13 +984,13 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { if( fHistoArray[EXCORE] ) { - fHistoArray[EXCORE]->fill( fData->Xcore, weight, log10( fData->get_Erec() ) ); + fHistoArray[EXCORE]->fill( fData->Xcore, weight, erec_log10 ); } if( fInput == 0 ) { if( fHistoArray[EYCORE] ) { - fHistoArray[EYCORE]->fill(-1.*fData->Ycore, weight, log10( fData->get_Erec() ) ); + fHistoArray[EYCORE]->fill(-1.*fData->Ycore, weight, erec_log10 ); } hXYcore->Fill( fData->Xcore, -1.*fData->Ycore, weight ); hAzYcore->Fill( fData->Az, -1.*fData->Ycore, weight ); @@ -995,7 +999,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { if( fHistoArray[EYCORE] ) { - fHistoArray[EYCORE]->fill( fData->Ycore, weight, log10( fData->get_Erec() ) ); + fHistoArray[EYCORE]->fill( fData->Ycore, weight, erec_log10 ); } hXYcore->Fill( fData->Xcore, fData->Ycore, weight ); hAzYcore->Fill( fData->Az, fData->Ycore, weight ); @@ -1035,133 +1039,133 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts UInt_t iTelImage = fData->ImgSel_list[t]; // telescope wise quality cuts - if( fData->ntubes[iTelImage] > ntubes_min && fData->size[t] > 0. && fData->get_Erec() > 0. ) + if( fData->ntubes[iTelImage] > ntubes_min && fData->size[t] > 0. && erec > 0. ) { if( fHistoSingleTel[ELENGTH][iTelImage] ) { - fHistoSingleTel[ELENGTH][iTelImage]->fill( fData->length[iTelImage], weight, log10( fData->get_Erec() ), + fHistoSingleTel[ELENGTH][iTelImage]->fill( fData->length[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EWIDTH][iTelImage] ) { - fHistoSingleTel[EWIDTH][iTelImage]->fill( fData->width[iTelImage], weight, log10( fData->get_Erec() ), + fHistoSingleTel[EWIDTH][iTelImage]->fill( fData->width[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ENTUBES][iTelImage] ) { - fHistoSingleTel[ENTUBES][iTelImage]->fill( fData->ntubes[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ENTUBES][iTelImage]->fill( fData->ntubes[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ENLOWGAIN][iTelImage] && fData->nlowgain[iTelImage] ) { - fHistoSingleTel[ENLOWGAIN][iTelImage]->fill( fData->nlowgain[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ENLOWGAIN][iTelImage]->fill( fData->nlowgain[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EDIST][iTelImage] ) { - fHistoSingleTel[EDIST][iTelImage]->fill( fData->dist[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[EDIST][iTelImage]->fill( fData->dist[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ESIZE][iTelImage] ) { - fHistoSingleTel[ESIZE][iTelImage]->fill( log10( fData->size[iTelImage] ), weight, log10( fData->get_Erec() ), + fHistoSingleTel[ESIZE][iTelImage]->fill( log10( fData->size[iTelImage] ), weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ESIZEHG][iTelImage] ) { fHistoSingleTel[ESIZEHG][iTelImage]->fill( log10( fData->size[iTelImage] - fData->size[iTelImage]*fData->fraclow[iTelImage] ), - weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ESIZELG][iTelImage] && fData->fraclow[iTelImage] > 0. ) { fHistoSingleTel[ESIZELG][iTelImage]->fill( log10( fData->size[iTelImage]*fData->fraclow[iTelImage] ), weight, - log10( fData->get_Erec() ), fData->ntubes[iTelImage], + erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EFRACLOW][iTelImage] ) { - fHistoSingleTel[EFRACLOW][iTelImage]->fill( fData->fraclow[iTelImage], weight, log10( fData->get_Erec() ), + fHistoSingleTel[EFRACLOW][iTelImage]->fill( fData->fraclow[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EMAX1][iTelImage] && fData->max1[iTelImage] > 0. ) { - fHistoSingleTel[EMAX1][iTelImage]->fill( log10( fData->max1[iTelImage] ), weight, log10( fData->get_Erec() ), + fHistoSingleTel[EMAX1][iTelImage]->fill( log10( fData->max1[iTelImage] ), weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EMAX2][iTelImage] && fData->max2[iTelImage] > 0. ) { - fHistoSingleTel[EMAX2][iTelImage]->fill( log10( fData->max2[iTelImage] ), weight, log10( fData->get_Erec() ), + fHistoSingleTel[EMAX2][iTelImage]->fill( log10( fData->max2[iTelImage] ), weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EMAX3][iTelImage] && fData->max3[iTelImage] > 0. ) { - fHistoSingleTel[EMAX3][iTelImage]->fill( log10( fData->max3[iTelImage] ), weight, log10( fData->get_Erec() ), + fHistoSingleTel[EMAX3][iTelImage]->fill( log10( fData->max3[iTelImage] ), weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EALPHA][iTelImage] ) { - fHistoSingleTel[EALPHA][iTelImage]->fill( fData->alpha[iTelImage], weight, log10( fData->get_Erec() ), + fHistoSingleTel[EALPHA][iTelImage]->fill( fData->alpha[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ELOS][iTelImage] && fData->size[iTelImage] > 0. ) { fHistoSingleTel[ELOS][iTelImage]->fill( fData->length[iTelImage] / fData->size[iTelImage] * 1000., weight, - log10( fData->get_Erec() ), fData->ntubes[iTelImage], + erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ELOSS][iTelImage] ) { - fHistoSingleTel[ELOSS][iTelImage]->fill( fData->loss[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ELOSS][iTelImage]->fill( fData->loss[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EASYM][iTelImage] ) { - fHistoSingleTel[EASYM][iTelImage]->fill( fData->asym[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[EASYM][iTelImage]->fill( fData->asym[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ECENX][iTelImage] ) { - fHistoSingleTel[ECENX][iTelImage]->fill( fData->cen_x[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ECENX][iTelImage]->fill( fData->cen_x[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ECENY][iTelImage] ) { - fHistoSingleTel[ECENY][iTelImage]->fill( fData->cen_y[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ECENY][iTelImage]->fill( fData->cen_y[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ETGRADX][iTelImage] ) { - fHistoSingleTel[ETGRADX][iTelImage]->fill( fData->tgrad_x[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ETGRADX][iTelImage]->fill( fData->tgrad_x[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[ETELDIST][iTelImage] ) { - fHistoSingleTel[ETELDIST][iTelImage]->fill( fData->R_core[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ETELDIST][iTelImage]->fill( fData->R_core[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } - if( fHistoSingleTel[ERECRAT][iTelImage] && fData->get_Erec() > 0. && fData->ES[iTelImage] > 0. ) + if( fHistoSingleTel[ERECRAT][iTelImage] && erec > 0. && fData->ES[iTelImage] > 0. ) { - fHistoSingleTel[ERECRAT][iTelImage]->fill( fData->ES[iTelImage] / fData->get_Erec(), weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[ERECRAT][iTelImage]->fill( fData->ES[iTelImage] / erec, weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EPEDVART][iTelImage] ) { - fHistoSingleTel[EPEDVART][iTelImage]->fill( fData->meanPedvar_ImageT[iTelImage], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[EPEDVART][iTelImage]->fill( fData->meanPedvar_ImageT[iTelImage], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } if( fHistoSingleTel[EDISPT][iTelImage] ) { - fHistoSingleTel[EDISPT][iTelImage]->fill( fData->Disp_T[t], weight, log10( fData->get_Erec() ), fData->ntubes[iTelImage], + fHistoSingleTel[EDISPT][iTelImage]->fill( fData->Disp_T[t], weight, erec_log10, fData->ntubes[iTelImage], fData->size[iTelImage], fData->fraclow[iTelImage] ); } - if( fHistoSingleTel[ERECRAT][iTelImage] && fData->get_Erec() > 0. ) + if( fHistoSingleTel[ERECRAT][iTelImage] && erec > 0. ) if( hcen_xy[iTelImage] ) { hcen_xy[iTelImage]->Fill( fData->cen_x[iTelImage], fData->cen_y[iTelImage], weight ); @@ -1283,7 +1287,7 @@ TH1D* VDataMCComparision::getAzimuthWeightingHistogram( string ifile ) continue; } // successful energy reconstruction - if( tData->EChi2S < 0 || tData->get_Erec() <= 0 ) + if( tData->EChi2S < 0 || tData->get_Erec(fEnergyReconstructionMethod) <= 0 ) { continue; } diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 85cbf19d..54c07b31 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -953,8 +953,8 @@ bool VGammaHadronCuts::applyStereoQualityCuts( bool bCount, int iEntry, bool fIs // apply cut on difference between stereo intersection and disp method // (for stereo method only: this should always pass) float i_disp_diff = sqrt( - ( fData->get_Xoff() - fData->Xoff_intersect ) * ( fData->get_Xoff() - fData->Xoff_intersect ) + - ( fData->get_Yoff() - fData->Yoff_intersect ) * ( fData->get_Yoff() - fData->Yoff_intersect ) ); + ( fData->get_Xoff(0) - fData->Xoff_intersect ) * ( fData->get_Xoff(0) - fData->Xoff_intersect ) + + ( fData->get_Yoff(0) - fData->Yoff_intersect ) * ( fData->get_Yoff(0) - fData->Yoff_intersect ) ); if( fCut_DispIntersectSuccess && ( fData->Xoff_intersect < -90. || fData->Yoff_intersect < -90. ) ) { if( bCount && fStats ) @@ -1350,7 +1350,7 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile, { TDirectory* cDir = gDirectory; - fTMVAEvaluator = new VTMVAEvaluator(); + fTMVAEvaluator = new VTMVAEvaluator(fEnergyReconstructionMethod); fTMVAEvaluator->setDebug( fDebug ); // constant signal efficiency diff --git a/src/VRadialAcceptance.cpp b/src/VRadialAcceptance.cpp index 7b1fe048..b7e7b7ba 100644 --- a/src/VRadialAcceptance.cpp +++ b/src/VRadialAcceptance.cpp @@ -615,17 +615,20 @@ int VRadialAcceptance::fillAcceptanceFromData( CData* iData, int entry ) // no more cuts after this statement bPassed = true; - idist = sqrt( iData->Xoff* iData->Xoff + iData->Yoff* iData->Yoff ); + double xoff = iData->get_Xoff(fDirectionReconstructionMethod); + double yoff = iData->get_Yoff(fDirectionReconstructionMethod); + + idist = sqrt( xoff*xoff + yoff*yoff); // fill 2D distribution of events - hXYAccTot->Fill( iData->get_Xoff(), iData->get_Yoff() ); - pair xy_derot = iData->get_XYoff_derot(); + hXYAccTot->Fill( xoff, yoff ); + pair xy_derot = iData->get_XYoff_derot(fDirectionReconstructionMethod); hXYAccTotDeRot->Fill( xy_derot.first, xy_derot.second ); hXYAccImgSel[iData->ImgSel]->Fill( xy_derot.first, xy_derot.second ) ; - hXYAccImgSelPreDeRot[iData->ImgSel]->Fill( iData->get_Xoff(), iData->get_Yoff() ) ; + hXYAccImgSelPreDeRot[iData->ImgSel]->Fill( xoff, yoff ); hXYAccNImages[iData->NImages]->Fill( xy_derot.first, xy_derot.second ) ; - hXYAccNImagesPreDeRot[iData->NImages]->Fill( iData->get_Xoff(), iData->get_Yoff() ) ; + hXYAccNImagesPreDeRot[iData->NImages]->Fill( xoff, yoff ); // 1D histograms // Radius Dependent Histograms @@ -666,7 +669,7 @@ int VRadialAcceptance::fillAcceptanceFromData( CData* iData, int entry ) } // fill azimuth angle dependent histograms (camera coordinates) - i_Phi = atan2( iData->get_Yoff(), iData->get_Xoff() ) * TMath::RadToDeg(); + i_Phi = atan2( yoff, xoff ) * TMath::RadToDeg(); hPhiDist->Fill( i_Phi ); for( unsigned int j = 0; j < fPhiMin.size(); j++ ) @@ -735,7 +738,7 @@ int VRadialAcceptance::fillAcceptanceFromData( CData* iData, int entry ) } if( j < hXYAccRun.size() ) { - hXYAccRun[j]->Fill( iData->get_Xoff(), iData->get_Yoff() ); + hXYAccRun[j]->Fill( xoff, yoff ); } break; } diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp index adeaa657..5502c3cd 100644 --- a/src/VTMVAEvaluator.cpp +++ b/src/VTMVAEvaluator.cpp @@ -5,7 +5,7 @@ #include "VTMVAEvaluator.h" -VTMVAEvaluator::VTMVAEvaluator() +VTMVAEvaluator::VTMVAEvaluator(unsigned int iEnergyReconstructionMethod) { fIsZombie = false; @@ -14,6 +14,7 @@ VTMVAEvaluator::VTMVAEvaluator() fData = 0; fTMVAEvaluatorResults = 0; fAverageZenithPerRun = -9999.; + fEnergyReconstructionMethod = iEnergyReconstructionMethod; reset(); } @@ -716,11 +717,11 @@ bool VTMVAEvaluator::evaluate( bool interpolate_mva, bool use_average_zenith_ang { i_ze = fAverageZenithPerRun; } - if( fData->get_Erec() <= 0. ) + if( fData->get_Erec(fEnergyReconstructionMethod) <= 0. ) { return false; } - unsigned int iDataBin = getDataBin( log10( fData->get_Erec() ), i_ze ); + unsigned int iDataBin = getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), i_ze ); if( fDebug ) { cout << "VTMVAEvaluator::evaluate: data bin " << iDataBin; @@ -760,14 +761,14 @@ bool VTMVAEvaluator::evaluate( bool interpolate_mva, bool use_average_zenith_ang */ double VTMVAEvaluator::interpolate_mva_evaluation() { - if(!fData || fData->get_Erec() <= 0 ) + if(!fData || fData->get_Erec(fEnergyReconstructionMethod) <= 0 ) { return 9999.; } set data_bins; for( unsigned int i = 0; i < fTMVAData.size(); i++ ) { - data_bins.insert( getDataBin( log10( fData->get_Erec() ), 0.5 * ( fTMVAData[i]->fZenithCut_max + fTMVAData[i]->fZenithCut_min ) ) ); + data_bins.insert( getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), 0.5 * ( fTMVAData[i]->fZenithCut_max + fTMVAData[i]->fZenithCut_min ) ) ); } TGraph iG(( int )data_bins.size() ); unsigned int i = 0; @@ -790,11 +791,11 @@ double VTMVAEvaluator::interpolate_mva_evaluation() */ unsigned int VTMVAEvaluator::getDataBin() { - if(!fData || fData->get_Erec() <= 0 ) + if(!fData || fData->get_Erec(fEnergyReconstructionMethod) <= 0 ) { return 9999; } - return getDataBin( log10( fData->get_Erec() ), fData->Ze ); + return getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), fData->Ze ); } From 15df7682de685d25c438f0285d4abece0ba293dd Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 16:28:28 +0100 Subject: [PATCH 17/21] astyle --- inc/VTMVAEvaluator.h | 2 +- src/VDataMCComparision.cpp | 6 +++--- src/VGammaHadronCuts.cpp | 6 +++--- src/VRadialAcceptance.cpp | 8 ++++---- src/VTMVAEvaluator.cpp | 14 +++++++------- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h index 1246ebe5..1e14dd32 100644 --- a/inc/VTMVAEvaluator.h +++ b/inc/VTMVAEvaluator.h @@ -172,7 +172,7 @@ class VTMVAEvaluator : public TNamed, public VPlotUtilities public: - VTMVAEvaluator(unsigned int iEnergyReconstructionMethod=0); + VTMVAEvaluator( unsigned int iEnergyReconstructionMethod = 0 ); ~VTMVAEvaluator() {}; bool evaluate( bool interpolate_mva = false, bool use_average_zenith_angle = true ); diff --git a/src/VDataMCComparision.cpp b/src/VDataMCComparision.cpp index 952dbc5d..90d02e6d 100644 --- a/src/VDataMCComparision.cpp +++ b/src/VDataMCComparision.cpp @@ -694,7 +694,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts ///////////////////////////////////////////////// ///////////////////////////////////////////////// // quality cuts - erec = fData->get_Erec(fEnergyReconstructionMethod); + erec = fData->get_Erec( fEnergyReconstructionMethod ); // nimage cut if( fData->NImages < fNImages_min ) @@ -706,7 +706,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts { continue; } - erec_log10 = log10(erec); + erec_log10 = log10( erec ); // successful lookup table reconstruction if( fData->MSCW < -50. || fData->MSCL < -50. ) { @@ -1287,7 +1287,7 @@ TH1D* VDataMCComparision::getAzimuthWeightingHistogram( string ifile ) continue; } // successful energy reconstruction - if( tData->EChi2S < 0 || tData->get_Erec(fEnergyReconstructionMethod) <= 0 ) + if( tData->EChi2S < 0 || tData->get_Erec( fEnergyReconstructionMethod ) <= 0 ) { continue; } diff --git a/src/VGammaHadronCuts.cpp b/src/VGammaHadronCuts.cpp index 54c07b31..7966c51c 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -953,8 +953,8 @@ bool VGammaHadronCuts::applyStereoQualityCuts( bool bCount, int iEntry, bool fIs // apply cut on difference between stereo intersection and disp method // (for stereo method only: this should always pass) float i_disp_diff = sqrt( - ( fData->get_Xoff(0) - fData->Xoff_intersect ) * ( fData->get_Xoff(0) - fData->Xoff_intersect ) + - ( fData->get_Yoff(0) - fData->Yoff_intersect ) * ( fData->get_Yoff(0) - fData->Yoff_intersect ) ); + ( fData->get_Xoff( 0 ) - fData->Xoff_intersect ) * ( fData->get_Xoff( 0 ) - fData->Xoff_intersect ) + + ( fData->get_Yoff( 0 ) - fData->Yoff_intersect ) * ( fData->get_Yoff( 0 ) - fData->Yoff_intersect ) ); if( fCut_DispIntersectSuccess && ( fData->Xoff_intersect < -90. || fData->Yoff_intersect < -90. ) ) { if( bCount && fStats ) @@ -1350,7 +1350,7 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile, { TDirectory* cDir = gDirectory; - fTMVAEvaluator = new VTMVAEvaluator(fEnergyReconstructionMethod); + fTMVAEvaluator = new VTMVAEvaluator( fEnergyReconstructionMethod ); fTMVAEvaluator->setDebug( fDebug ); // constant signal efficiency diff --git a/src/VRadialAcceptance.cpp b/src/VRadialAcceptance.cpp index b7e7b7ba..6bdf20c6 100644 --- a/src/VRadialAcceptance.cpp +++ b/src/VRadialAcceptance.cpp @@ -615,14 +615,14 @@ int VRadialAcceptance::fillAcceptanceFromData( CData* iData, int entry ) // no more cuts after this statement bPassed = true; - double xoff = iData->get_Xoff(fDirectionReconstructionMethod); - double yoff = iData->get_Yoff(fDirectionReconstructionMethod); + double xoff = iData->get_Xoff( fDirectionReconstructionMethod ); + double yoff = iData->get_Yoff( fDirectionReconstructionMethod ); - idist = sqrt( xoff*xoff + yoff*yoff); + idist = sqrt( xoff* xoff + yoff* yoff ); // fill 2D distribution of events hXYAccTot->Fill( xoff, yoff ); - pair xy_derot = iData->get_XYoff_derot(fDirectionReconstructionMethod); + pair xy_derot = iData->get_XYoff_derot( fDirectionReconstructionMethod ); hXYAccTotDeRot->Fill( xy_derot.first, xy_derot.second ); hXYAccImgSel[iData->ImgSel]->Fill( xy_derot.first, xy_derot.second ) ; diff --git a/src/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp index 5502c3cd..4d272886 100644 --- a/src/VTMVAEvaluator.cpp +++ b/src/VTMVAEvaluator.cpp @@ -5,7 +5,7 @@ #include "VTMVAEvaluator.h" -VTMVAEvaluator::VTMVAEvaluator(unsigned int iEnergyReconstructionMethod) +VTMVAEvaluator::VTMVAEvaluator( unsigned int iEnergyReconstructionMethod ) { fIsZombie = false; @@ -717,11 +717,11 @@ bool VTMVAEvaluator::evaluate( bool interpolate_mva, bool use_average_zenith_ang { i_ze = fAverageZenithPerRun; } - if( fData->get_Erec(fEnergyReconstructionMethod) <= 0. ) + if( fData->get_Erec( fEnergyReconstructionMethod ) <= 0. ) { return false; } - unsigned int iDataBin = getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), i_ze ); + unsigned int iDataBin = getDataBin( log10( fData->get_Erec( fEnergyReconstructionMethod ) ), i_ze ); if( fDebug ) { cout << "VTMVAEvaluator::evaluate: data bin " << iDataBin; @@ -761,14 +761,14 @@ bool VTMVAEvaluator::evaluate( bool interpolate_mva, bool use_average_zenith_ang */ double VTMVAEvaluator::interpolate_mva_evaluation() { - if(!fData || fData->get_Erec(fEnergyReconstructionMethod) <= 0 ) + if(!fData || fData->get_Erec( fEnergyReconstructionMethod ) <= 0 ) { return 9999.; } set data_bins; for( unsigned int i = 0; i < fTMVAData.size(); i++ ) { - data_bins.insert( getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), 0.5 * ( fTMVAData[i]->fZenithCut_max + fTMVAData[i]->fZenithCut_min ) ) ); + data_bins.insert( getDataBin( log10( fData->get_Erec( fEnergyReconstructionMethod ) ), 0.5 * ( fTMVAData[i]->fZenithCut_max + fTMVAData[i]->fZenithCut_min ) ) ); } TGraph iG(( int )data_bins.size() ); unsigned int i = 0; @@ -791,11 +791,11 @@ double VTMVAEvaluator::interpolate_mva_evaluation() */ unsigned int VTMVAEvaluator::getDataBin() { - if(!fData || fData->get_Erec(fEnergyReconstructionMethod) <= 0 ) + if(!fData || fData->get_Erec( fEnergyReconstructionMethod ) <= 0 ) { return 9999; } - return getDataBin( log10( fData->get_Erec(fEnergyReconstructionMethod) ), fData->Ze ); + return getDataBin( log10( fData->get_Erec( fEnergyReconstructionMethod ) ), fData->Ze ); } From 67839cf642e616273ad5d8e935860c027a949db8 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 16:43:10 +0100 Subject: [PATCH 18/21] copilot review --- src/CData.cpp | 4 ++-- src/VAnaSumRunParameter.cpp | 8 ++++---- src/VInstrumentResponseFunctionRunParameter.cpp | 5 ++++- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/CData.cpp b/src/CData.cpp index b5f7b7a0..4c6e37fc 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -718,7 +718,7 @@ float CData::get_ErecChi2( unsigned iMethod ) { return 0.; } - return EChi2S; + return EChi2; } /* @@ -1094,7 +1094,7 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name cout << "CData Error: cannot find " << tree_name << " tree in " << file_name << endl; exit( EXIT_FAILURE ); } - cout << "Adding " << tree_name << " tree from " << file_name << endl; + cout << "\t Adding " << tree_name << " tree from " << file_name << endl; return iXGB_tree; } diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp index 1d93e534..ba6ad7c1 100644 --- a/src/VAnaSumRunParameter.cpp +++ b/src/VAnaSumRunParameter.cpp @@ -1193,7 +1193,7 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) cout << " (use effective area A_REC)"; } cout << endl; - cout << "Energy reconstruction method " << fEnergyReconstructionMethod; + cout << "\t Energy reconstruction method " << fEnergyReconstructionMethod; if( fEnergyReconstructionMethod == 0 ) { cout << " (dispBDT energy reconstruction)" << endl; @@ -1206,14 +1206,14 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) { cout << " (XGB stereo reconstruction)" << endl; } - cout << "Direction reconstruction method " << fDirectionReconstructionMethod; + cout << "\t Direction reconstruction method " << fDirectionReconstructionMethod; if( fDirectionReconstructionMethod == 0 ) { - cout << " (dispBDT energy reconstruction)" << endl; + cout << " (dispBDT direction reconstruction)" << endl; } else if( fDirectionReconstructionMethod == 1 ) { - cout << " (lookup table energy reconstruction)" << endl; + cout << " (intersection method reconstruction)" << endl; } else if( fDirectionReconstructionMethod == 2 ) { diff --git a/src/VInstrumentResponseFunctionRunParameter.cpp b/src/VInstrumentResponseFunctionRunParameter.cpp index ffb50628..114f8d24 100644 --- a/src/VInstrumentResponseFunctionRunParameter.cpp +++ b/src/VInstrumentResponseFunctionRunParameter.cpp @@ -359,6 +359,10 @@ bool VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile( stri } } cout << "========================================" << endl << endl; + if( fIgnoreEnergyReconstructionQuality ) + { + fEnergyReconstructionMethod = 100; + } ////////////////////////////////////////////////////////////////////////////////////// // fill some parameters @@ -662,7 +666,6 @@ void VInstrumentResponseFunctionRunParameter::print() if( fIgnoreEnergyReconstructionQuality ) { cout << ", ignoring cut on quality of energy reconstruction"; - fEnergyReconstructionMethod = 100; } cout << endl; cout << "energy reconstruction method " << fEnergyReconstructionMethod << endl; From 35e8a89991bafa7a7c2a8e2ee6eb593ab01b6926 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 17:19:18 +0100 Subject: [PATCH 19/21] chi2 cut --- src/CData.cpp | 8 ++++---- src/VStereoAnalysis.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/CData.cpp b/src/CData.cpp index 4c6e37fc..ed14aa13 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -712,11 +712,11 @@ float CData::get_ErecChi2( unsigned iMethod ) // not defined for XGB methods else if( iMethod == 2 ) { - return 0.; + return 1.e-2; } else if( iMethod == 100 ) { - return 0.; + return 1.e-2; } return EChi2; } @@ -733,11 +733,11 @@ float CData::get_ErecdE( unsigned iMethod ) // not defined for XGB methods else if( iMethod == 2 ) { - return 0.; + return 1.e-2; } else if( iMethod == 100 ) { - return 0.; + return 1.e-2; } return dE; } diff --git a/src/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index 6d0705ae..3147257f 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -575,7 +575,7 @@ double VStereoAnalysis::fillHistograms( int icounter, int irun, double iAzMin, d fHisto[fHisCounter]->hemissC2->Fill( fDataRun->EmissionHeightChi2 ); } // chi2 of energy reconstruction - if( iErecChi2 > 0. ) + if( iErecChi2 >= 0. ) { fHisto[fHisCounter]->herecChi2->Fill( iErecChi2 ); } From 23b023575c69a68dc6e9465312366376d6edd1c2 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 19:38:26 +0100 Subject: [PATCH 20/21] fix initalize tree --- docs/changes/337.feauture.md | 1 + inc/CData.h | 2 +- src/CData.cpp | 4 ++-- 3 files changed, 4 insertions(+), 3 deletions(-) create mode 100644 docs/changes/337.feauture.md diff --git a/docs/changes/337.feauture.md b/docs/changes/337.feauture.md new file mode 100644 index 00000000..94489a35 --- /dev/null +++ b/docs/changes/337.feauture.md @@ -0,0 +1 @@ +Add XGB-based gamma/hadron classification routines. Use with `* cutselection 52 0` in cut file. diff --git a/inc/CData.h b/inc/CData.h index 2c2101b5..b1cc9a72 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -278,7 +278,7 @@ class CData float get_ErecdE( unsigned int method ); float get_Xoff( unsigned int method ); float get_Yoff( unsigned int method ); - void initialize_xgb_tree( TTree* stereoTree, TTree* ghTree ); + void initialize_xgb_tree(); TTree* getXGBTree( string file_name, string suffix, string tree_name ); pair get_XYoff_derot( unsigned int method ); virtual Long64_t LoadTree( Long64_t entry ); diff --git a/src/CData.cpp b/src/CData.cpp index ed14aa13..5647a4cc 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -658,7 +658,7 @@ Bool_t CData::Notify() } /* - Get Gamma/hadron decision from XGB friend tree. + Get Get gamma prediction score. */ float CData::get_GH_Gamma_Prediction() { @@ -1102,7 +1102,7 @@ TTree* CData::getXGBTree( string file_name, string file_suffix, string tree_name * Initialize XGB trees as kind of friends * */ -void CData::initialize_xgb_tree( TTree* stereoTree, TTree* ghTree ) +void CData::initialize_xgb_tree() { if( fStereoFriendTree ) { From b611ca6576f4e08c7c096dc1670f4bb62c6a5e27 Mon Sep 17 00:00:00 2001 From: GernotMaier Date: Sun, 11 Jan 2026 19:40:24 +0100 Subject: [PATCH 21/21] fix init --- src/CData.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/CData.cpp b/src/CData.cpp index 5647a4cc..7bacdbad 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -20,8 +20,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, TTree* stereoTree, TTree* ghTr fStereoFriendTree = stereoTree; fGHFriendTree = ghTree; - - initialize_xgb_tree( fStereoFriendTree, fGHFriendTree ); + initialize_xgb_tree(); } @@ -36,7 +35,7 @@ CData::CData( TTree* tree, bool bMC, bool bShort, string file_name, string stere fStereoFriendTree = getXGBTree( file_name, stereo_suffix, "StereoAnalysis" ); fGHFriendTree = getXGBTree( file_name, gamma_hadron_suffix, "Classification" ); - initialize_xgb_tree( fStereoFriendTree, fGHFriendTree ); + initialize_xgb_tree(); }