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/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. diff --git a/inc/CData.h b/inc/CData.h index 5551224e..b1cc9a72 100644 --- a/inc/CData.h +++ b/inc/CData.h @@ -260,19 +260,30 @@ class CData TBranch* b_Xoff_intersect; TBranch* b_Yoff_intersect; - TTree* fFriendTree; //! - float Dir_Xoff; //! - float Dir_Yoff; //! - float Dir_Erec; //! + TTree* fStereoFriendTree; //! + float Dir_Xoff; //! + float Dir_Yoff; //! + float Dir_Erec; //! + TTree* fGHFriendTree; //! + float GH_Gamma_Prediction; //! + UChar_t GH_Is_Gamma; //! + vector fXGBFiles; //! - 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 ); + 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 ); - pair get_XYoff_derot( 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* getXGBTree( string file_name, string suffix, string tree_name ); + pair get_XYoff_derot( unsigned int method ); 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/inc/VAnaSumRunParameter.h b/inc/VAnaSumRunParameter.h index aa66c41e..0d971e43 100644 --- a/inc/VAnaSumRunParameter.h +++ b/inc/VAnaSumRunParameter.h @@ -190,11 +190,15 @@ 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; // XGB reconstruction - string fXGB_file_suffix; + string fXGB_stereo_file_suffix; + string fXGB_gh_file_suffix; int f2DAcceptanceMode ; // USE2DACCEPTANCE 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 a19c1c9a..b11a5285 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 }; //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -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,9 +199,10 @@ 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 ); double getArrayCentre_X() { @@ -207,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; @@ -253,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() { @@ -304,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; @@ -317,11 +320,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/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 c0f49312..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; @@ -74,7 +75,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/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/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/inc/VTMVAEvaluator.h b/inc/VTMVAEvaluator.h index 1fdde49f..1e14dd32 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/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 6b1b02f8..7bacdbad 100644 --- a/src/CData.cpp +++ b/src/CData.cpp @@ -10,27 +10,32 @@ #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; + Init( tree ); + + fStereoFriendTree = stereoTree; + fGHFriendTree = ghTree; + initialize_xgb_tree(); +} + - fFriendTree = friendTree; + +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 ); - if( fFriendTree ) - { - fFriendTree->SetBranchAddress( "Dir_Xoff", &Dir_Xoff ); - fFriendTree->SetBranchAddress( "Dir_Yoff", &Dir_Yoff ); - fFriendTree->SetBranchAddress( "Dir_Erec", &Dir_Erec ); - } - else - { - Dir_Xoff = -9999.; - Dir_Yoff = -9999.; - Dir_Erec = -9999.; - } + + fStereoFriendTree = getXGBTree( file_name, stereo_suffix, "StereoAnalysis" ); + fGHFriendTree = getXGBTree( file_name, gamma_hadron_suffix, "Classification" ); + initialize_xgb_tree(); } @@ -40,6 +45,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(); } @@ -53,9 +67,13 @@ Int_t CData::GetEntry( Long64_t entry ) } int a = fChain->GetEntry( entry ); - if( fFriendTree ) + if( fStereoFriendTree ) + { + fStereoFriendTree->GetEntry( entry ); + } + if( fGHFriendTree ) { - fFriendTree->GetEntry( entry ); + fGHFriendTree->GetEntry( entry ); } if( fTelescopeCombination > 0 && fTelescopeCombination != 15 ) @@ -638,35 +656,89 @@ Bool_t CData::Notify() return kTRUE; } +/* + Get Get gamma prediction score. +*/ +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 * 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 && fFriendTree ) + 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 1.e-2; } - else if( iMethod == 3 ) + else if( iMethod == 100 ) { - return Dir_Erec; + return 1.e-2; } - return Erec; + return EChi2; +} + +/* + * Get energy dE + */ +float CData::get_ErecdE( unsigned iMethod ) +{ + if( iMethod == 1 ) + { + return dES; + } + // not defined for XGB methods + else if( iMethod == 2 ) + { + return 1.e-2; + } + else if( iMethod == 100 ) + { + return 1.e-2; + } + return dE; } /* @@ -674,26 +746,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 && fFriendTree ) - { - 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; } @@ -706,26 +769,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 && fFriendTree ) - { - 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; } @@ -1009,3 +1063,66 @@ void CData::initialize_3tel_reconstruction( fTelY = tel_y; fTelZ = tel_z; } + + +/* + 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 ) +{ + if( file_suffix == "" || file_suffix == "None" ) + { + return 0; + } + + 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 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 ); + + TTree* iXGB_tree = ( TTree* )iFile->Get( tree_name.c_str() ); + if(!iXGB_tree ) + { + cout << "CData Error: cannot find " << tree_name << " tree in " << file_name << endl; + exit( EXIT_FAILURE ); + } + cout << "\t Adding " << tree_name << " tree from " << file_name << endl; + return iXGB_tree; +} + +/* + * Initialize XGB trees as kind of friends + * +*/ +void CData::initialize_xgb_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.; + } + if( fGHFriendTree ) + { + fGHFriendTree->SetBranchAddress( "Gamma_Prediction", &GH_Gamma_Prediction ); + fGHFriendTree->SetBranchAddress( "Is_Gamma_70", &GH_Is_Gamma ); + } + else + { + GH_Gamma_Prediction = -9999.; + GH_Is_Gamma = false; + } +} diff --git a/src/VAnaSumRunParameter.cpp b/src/VAnaSumRunParameter.cpp index 151a72d5..ba6ad7c1 100644 --- a/src/VAnaSumRunParameter.cpp +++ b/src/VAnaSumRunParameter.cpp @@ -117,12 +117,14 @@ 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; fEnergyEffectiveAreaSmoothingThreshold = -1.; fDeadTimeCalculationMethod = 0; - fXGB_file_suffix = ""; + fXGB_stereo_file_suffix = ""; + fXGB_gh_file_suffix = ""; // background model fTMPL_fBackgroundModel = 0; @@ -582,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" ) @@ -598,10 +612,15 @@ int VAnaSumRunParameter::readRunParameter( string i_filename ) return 0; } } - else if( temp == "XGBFILESUFFIX" ) + else if( temp == "XGBSTEREOFILESUFFIX" ) + { + fXGB_stereo_file_suffix = temp2; + if( fXGB_stereo_file_suffix == "None" ) fXGB_stereo_file_suffix = ""; + } + else if( temp == "XGBGAMMAHADRONFILESUFFIX" ) { - fXGB_file_suffix = temp2; - if( fXGB_file_suffix == "None" ) fXGB_file_suffix = ""; + fXGB_gh_file_suffix = temp2; + if( fXGB_gh_file_suffix == "None" ) fXGB_gh_file_suffix = ""; } else if( temp == "RATEINTERVALLLENGTH" ) { @@ -1173,22 +1192,48 @@ void VAnaSumRunParameter::printStereoParameter( unsigned int i ) { cout << " (use effective area A_REC)"; } - cout << ", Method " << fEnergyReconstructionMethod; + cout << endl; + cout << "\t Energy reconstruction method " << fEnergyReconstructionMethod; if( fEnergyReconstructionMethod == 0 ) { cout << " (dispBDT energy reconstruction)" << endl; } - else + else if( fEnergyReconstructionMethod == 1 ) { cout << " (lookup table energy reconstruction)" << endl; } - if( fXGB_file_suffix != "" && fXGB_file_suffix != "None" ) + else if( fEnergyReconstructionMethod == 2 ) + { + cout << " (XGB stereo reconstruction)" << endl; + } + cout << "\t Direction reconstruction method " << fDirectionReconstructionMethod; + if( fDirectionReconstructionMethod == 0 ) + { + cout << " (dispBDT direction reconstruction)" << endl; + } + else if( fDirectionReconstructionMethod == 1 ) + { + cout << " (intersection method 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; + } + 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 suffix: " << 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 ) @@ -1466,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..90d02e6d 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 ); @@ -616,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; @@ -691,6 +694,7 @@ bool VDataMCComparision::fillHistograms( string ifile, int iSingleTelescopeCuts ///////////////////////////////////////////////// ///////////////////////////////////////////////// // quality cuts + erec = fData->get_Erec( fEnergyReconstructionMethod ); // nimage cut if( fData->NImages < fNImages_min ) @@ -698,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. ) { @@ -736,13 +741,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 @@ -861,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 ); } @@ -873,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++ ) { @@ -884,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] ); } @@ -892,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 @@ -907,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 @@ -920,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 ) @@ -959,13 +964,13 @@ 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] ) { fHistoArray[EMVA]->fill( - fCuts->getTMVA_EvaluationResult(), weight, log10( fData->get_Erec() ) ); + fCuts->getTMVA_EvaluationResult(), weight, erec_log10 ); } } } @@ -979,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 ); @@ -994,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 ); @@ -1034,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 ); @@ -1282,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/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 816067b7..7966c51c 100644 --- a/src/VGammaHadronCuts.cpp +++ b/src/VGammaHadronCuts.cpp @@ -13,6 +13,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: @@ -26,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 */ @@ -72,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(); } @@ -602,6 +607,14 @@ bool VGammaHadronCuts::readCuts( string i_cutfilename, int iPrint ) { fAnalysisType = MVAAnalysis; } + else if( fGammaHadronCutSelector / 10 == 5 ) + { + fAnalysisType = XGBoostAnalysis; + } + else + { + fAnalysisType = GEO; + } return true; } @@ -691,10 +704,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; @@ -706,12 +724,11 @@ void VGammaHadronCuts::printCutSummary() printTMVA_MVACut(); } } - // other cut parameters - if( fNTel == 2 ) + // XGBoost cuts + if( useXGBoostCuts() ) { - cout << ", size > " << fCut_Size_min; + cout << "XGBoost gamma/hadron separation with fixed 70\% signal efficiency" << endl; } - cout << endl; cout << "Fiducial area (camera) < " << fCut_CameraFiducialSize_max << " deg, "; cout << " stereo reconstruction: " << fCut_Chi2_min << " <= sChi2 <= " << fCut_Chi2_max << endl; cout << "Energy reconstruction: "; @@ -748,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 @@ -786,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. ) ) @@ -811,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 ) { @@ -936,8 +953,8 @@ bool VGammaHadronCuts::applyStereoQualityCuts( unsigned int iEnergyReconstructio // 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 ) @@ -966,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 ) { @@ -1010,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 ) @@ -1035,7 +1050,7 @@ bool VGammaHadronCuts::isGamma( int i, bool bCount, bool fIsOn ) } ///////////////////////////////////////////////////////////////////////////// // apply cut using TMVA reader - else if( useTMVACuts() ) + if( useTMVACuts() ) { if( fDebug ) { @@ -1050,6 +1065,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; } @@ -1080,6 +1110,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) " << ( bool )fData->GH_Is_Gamma; + cout << endl; + } + return ( bool )fData->GH_Is_Gamma; +} + /*! @@ -1303,7 +1350,7 @@ bool VGammaHadronCuts::initTMVAEvaluator( string iTMVAFile, { TDirectory* cDir = gDirectory; - fTMVAEvaluator = new VTMVAEvaluator(); + fTMVAEvaluator = new VTMVAEvaluator( fEnergyReconstructionMethod ); fTMVAEvaluator->setDebug( fDebug ); // constant signal efficiency @@ -1408,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 ) @@ -1487,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.; @@ -1517,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(); @@ -1564,10 +1612,6 @@ void VGammaHadronCuts::terminate() { fTMVAEvaluatorResults->Write(); } - else - { - cout << "No TMVAEvaluator Results." << endl; - } Write(); } @@ -1647,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 906c5993..114f8d24 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; @@ -51,7 +52,8 @@ VInstrumentResponseFunctionRunParameter::VInstrumentResponseFunctionRunParameter fdatafile = ""; fMCdatafile_tree = ""; fMCdatafile_histo = ""; - fXGB_file_suffix = ""; + fXGB_stereo_file_suffix = ""; + fXGB_gh_file_suffix = ""; fze = 0.; fnoise = 0; @@ -157,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" ) { @@ -236,12 +246,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 @@ -341,6 +359,10 @@ bool VInstrumentResponseFunctionRunParameter::readRunParameterFromTextFile( stri } } cout << "========================================" << endl << endl; + if( fIgnoreEnergyReconstructionQuality ) + { + fEnergyReconstructionMethod = 100; + } ////////////////////////////////////////////////////////////////////////////////////// // fill some parameters @@ -610,9 +632,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" ) { @@ -643,6 +669,7 @@ void VInstrumentResponseFunctionRunParameter::print() } 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..6bdf20c6 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 ) ) @@ -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/VStereoAnalysis.cpp b/src/VStereoAnalysis.cpp index de03062c..3147257f 100644 --- a/src/VStereoAnalysis.cpp +++ b/src/VStereoAnalysis.cpp @@ -11,8 +11,6 @@ VStereoAnalysis::VStereoAnalysis( bool ion, string i_hsuffix, VAnaSumRunParamete fDebug = false; fDataFile = 0; - fXGBFile = 0; - fXGB_tree = 0; fInstrumentEpochMinor = "NOT_SET"; fDirTot = iDirTot; fDirTotRun = iDirRun; @@ -111,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 ) { @@ -480,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 @@ -496,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; } @@ -519,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 @@ -576,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 ); } @@ -620,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 @@ -1963,34 +1962,14 @@ 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, + iFileName, + fRunPara->fXGB_stereo_file_suffix, + fRunPara->fXGB_gh_file_suffix + ); // read current (major) epoch from data file VEvndispRunParameter* i_runPara = ( VEvndispRunParameter* )fDataFile->Get( "runparameterV2" ); if( i_runPara ) @@ -2133,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; @@ -2144,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; @@ -2293,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/VTMVAEvaluator.cpp b/src/VTMVAEvaluator.cpp index adeaa657..4d272886 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 ); } 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 0b281939..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 ); @@ -175,39 +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 ); } - // 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->fdatafile, fRunPara->fXGB_stereo_file_suffix, fRunPara->fXGB_gh_file_suffix ); d.initialize_3tel_reconstruction( fRunPara->fRerunStereoReconstruction_3telescopes, fRunPara->fRerunStereoReconstruction_minAngle, @@ -342,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 ) )