diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/CMakeLists.txt b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..c2c6a30568d525bb749ef237540bca297c679296 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/CMakeLists.txt @@ -0,0 +1,53 @@ +# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + +# Declare the package name: +atlas_subdir( EGammaVariableCorrection ) + +# External dependencies: +find_package( ROOT COMPONENTS Core MathCore Hist Graph RIO ) + +# Component(s) in the package: +atlas_add_library( EGammaVariableCorrectionLib + EGammaVariableCorrection/*.h + Root/*.cxx + PUBLIC_HEADERS EGammaVariableCorrection + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools EgammaAnalysisInterfacesLib EgammaAnalysisHelpersLib + PRIVATE_LINK_LIBRARIES PathResolver xAODEgamma xAODEventShape) + +if( NOT XAOD_STANDALONE ) + atlas_add_component( EGammaVariableCorrection + src/components/*.cxx + LINK_LIBRARIES EGammaVariableCorrectionLib ) +endif() + +atlas_add_dictionary( ElectronPhotonVariableCorrectionDict + EGammaVariableCorrection/ElectronPhotonVariableCorrectionDict.h + EGammaVariableCorrection/selection.xml + LINK_LIBRARIES EGammaVariableCorrectionLib ) + +# Add tests +if( XAOD_STANDALONE ) + atlas_add_executable( Test_ElectronPhotonVariableCorrectionBase + util/testElectronPhotonVariableCorrectionBase.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} EGammaVariableCorrectionLib) + + atlas_add_executable( Test_IsoCorrection + util/testIsoCorrection.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} IsolationCorrectionsLib EGammaVariableCorrectionLib) + + atlas_add_executable( Test_ElectronPhotonVariableCorrectionTool + util/testElectronPhotonVariableCorrectionTool.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} EGammaVariableCorrectionLib) +endif() + +atlas_add_executable( Test_ElectronPhotonVariableCorrectionTool_DictionaryToolHandle + util/testElectronPhotonVariableCorrectionTool_DictionaryToolHandle.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} EGammaVariableCorrectionLib) + +# Install files from the package: +atlas_install_data( data/*.conf ) diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h new file mode 100644 index 0000000000000000000000000000000000000000..cf0fdb16fbfd02230f032aa20af5bff283c43dcd --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h @@ -0,0 +1,287 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ElectronPhotonVariableCorrectionBase_H +#define ElectronPhotonVariableCorrectionBase_H + +/** @class ElectronPhotonVariableCorrectionBase + * @brief Class to correct electron and photon MC variables. + * @details For a detailed documentation, please refer to the [gitlab readme](https://gitlab.cern.ch/atlas/athena/-/blob/21.2/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/README.md) + * + * @author Nils Gillwald (DESY) nils.gillwald@desy.de + * @date February 2020 +*/ + +//ATLAS includes +#include "AsgTools/AsgTool.h" +#include "PATInterfaces/CorrectionCode.h" + +//EDM includes +#include "xAODEgamma/Electron.h" +#include "xAODEgamma/Photon.h" + +//Root includes +#include "TF1.h" + +// forward declarations +class TGraph; +class TEnv; + +// =========================================================================== +// Class ElectronPhotonVariableCorrectionBase +// =========================================================================== + +class ElectronPhotonVariableCorrectionBase : public asg::AsgTool +{ + +public: + /** @brief Standard constructor + * @param myname Internal name of the class instance, so they can be distinguished + */ + ElectronPhotonVariableCorrectionBase(const std::string& myname); + //! @brief Standard destructor + ~ElectronPhotonVariableCorrectionBase() {}; + + /** @brief Initialize the class instance + * @details Reads out the configuration file set via the setProperty(ConfigFile, "/path/to/conf/file.conf") function + * and sets up the class instance accordingly + */ + virtual StatusCode initialize() override; + + /** @brief Apply the correction given in the conf file to the passed photon + * @param photon The photon which should be corrected + * @details Pass a photon which should be corrected by the tool to this function. The tool gets the original photon variable value, + * corrects it according to the configuration file given, and overwrites the original value of the variable. The original variable + * is saved as <variable_name>_original. Note that the variable to be corrected should be an auxiliary variable, and that the photon + * or photon container, respectively, must not be const if you want to use this function. If you are running on const objects, please + * make a deep / shallow copy of the container or use the correctedCopy function of this class. + */ + const CP::CorrectionCode applyCorrection( xAOD::Photon& photon ) const; + + /** @brief Apply the correction given in the conf file to the passed electron + * @param electron The electron which should be corrected + * @details Pass an electron which should be corrected by the tool to this function. The tool gets the original electron variable value, + * corrects it according to the configuration file given, and overwrites the original value of the variable. The original variable + * is saved as <variable_name>_original. Note that the variable to be corrected should be an auxiliary variable, and that the electron + * or electron container, respectively, must not be const if you want to use this function. If you are running on const objects, please + * make a deep / shallow copy of the container or use the correctedCopy function of this class. + */ + const CP::CorrectionCode applyCorrection( xAOD::Electron& electron ) const; + + /** @brief Make a corrected copy of the passed photon according to the given conf file + * @param in_photon The original photon of which a corrected copy should be made + * @param out_photon Empty new photon where the corrected copy will be stored in + * @details Makes a corrected copy of the passed photon object according to the correction provided in the configuration file. The in_photon + * is copied to the out_photon, which is then passed to the applyCorrection function. + */ + const CP::CorrectionCode correctedCopy( const xAOD::Photon& in_photon, xAOD::Photon*& out_photon ) const; + + /** @brief Make a corrected copy of the passed electron according to the given conf file + * @param in_electron The original electron of which a corrected copy should be made + * @param out_electron Empty new electron where the corrected copy will be stored in + * @details Makes a corrected copy of the passed electron object according to the correction provided in the configuration file. The in_electron + * is copied to the out_electron, which is then passed to the applyCorrection function. + */ + const CP::CorrectionCode correctedCopy( const xAOD::Electron& in_electron, xAOD::Electron*& out_electron) const; + + //! @brief Returns the variable which should be corrected according to the passed configuration file + const std::string getCorrectionVariable() { return m_correctionVariable; }; + + /** @brief Define the categories of EGamma objects tool can be applied to + * @details The tool can be applied to electrons and photons, but the latter could also have different corrections according to whether they are + * converted or unconverted photons. The categories are designed in such a way, that the analyzer can freely choose whether he_she_* wants all objects + * (i.e. e and y) to be corrected in the same way, or wheter he_she_* wants to apply different corrections to different types of objects. + */ + enum class EGammaObjects{ + Failure = 0, + unconvertedPhotons = 1, + convertedPhotons, + allPhotons, + allElectrons, + allEGammaObjects + }; //end enum EGammaObjects + + //! @brief Returns the type of EGamma object to which the correction should be applied. + ElectronPhotonVariableCorrectionBase::EGammaObjects isAppliedTo() { return m_applyToObjects; }; + + /** @brief Check if the ApplyTo flag passed in the conf file is compatible with converted photons + * @return True if the ApplyTo flag is compatible, False if not + */ + bool applyToConvertedPhotons() const; + + /** @brief Check if the ApplyTo flag passed in the conf file is compatible with uconverted photons + * @return True if the ApplyTo flag is compatible, False if not + */ + bool applyToUnconvertedPhotons() const; + + /** @brief Check if the ApplyTo flag passed in the conf file is compatible with electrons + * @return True if the ApplyTo flag is compatible, False if not + */ + bool applyToElectrons() const; //check if passed flag is compatible with electron + +private: + //! @brief Use enum and not string for type of function parameter in order to do faster comparisons + enum class parameterType{ + Failure = 0, + EtaDependentTGraph = 1, + PtDependentTGraph, + EtaBinned, + PtBinned, + EtaTimesPtBinned, + EventDensity + }; // end enum ParameterType + //! @brief The name of the configuration file + std::string m_configFile; + //! @brief The name of the variable to correct + std::string m_correctionVariable; + //! @brief Values of discontinuities in the variable which should not be corrected + std::vector<float> m_uncorrectedDiscontinuities; + //! @brief Function to use for the variable correction, TF1 style + std::string m_correctionFunctionString; + //! @brief The actual TF1 correction function + std::unique_ptr<TF1> m_correctionFunctionTF1; + //! @brief Number of parameters of the variable correction function + int m_numberOfFunctionParameters; + //! @brief Map of the correction function parameter number to the parameter type + std::vector<parameterType> m_ParameterTypeVector; + //! @brief Copy of the TGraph from the root file, stored if needed by the respective correction function parameter + std::vector<TGraph*> m_graphCopies; + //! @brief List of eta/pt dependent values, stored if needed by the respective correction function parameter + std::vector<std::vector<float>> m_binValues; + //! @brief List of bin boundaries in eta, stored if needed by any correction function parameter + std::vector<float> m_etaBins; + //! @brief List of bin boundaries in pT, stored if needed by any correction function parameter + std::vector<float> m_ptBins; + //! @brief List of bools whether a parameter should use linear interpolation in pT if it's some kind of pT binned parameter + std::vector<bool> m_interpolatePtFlags; + //! @brief The type of objects to which the specific conf file settings are allowed to be applied to + ElectronPhotonVariableCorrectionBase::EGammaObjects m_applyToObjects; + //! @brief Store if already retrieved eta binning + bool m_retrievedEtaBinning = false; + //! @brief Store if already retrieved pt binning + bool m_retrievedPtBinning = false; + //! @brief Accessor for the variable to be corrected + std::unique_ptr<SG::AuxElement::Accessor<float>> m_variableToCorrect; + //! @brief Accessor to store the original value of the corrected variable + std::unique_ptr<SG::AuxElement::Accessor<float>> m_originalVariable; + + /** @brief Convert input string to a parameter function type + * @param input The string to convert to ElectronPhotonVariableCorrectionBase::parameterType + * @return ElectronPhotonVariableCorrectionBase::parameterType, returns ElectronPhotonVariableCorrectionBase::parameterType::Failure if no match found + */ + ElectronPhotonVariableCorrectionBase::parameterType stringToParameterType( const std::string& input ) const; + + /** @brief Convert input string to egamma object type + * @param input The string to convert to ElectronPhotonVariableCorrectionBase::EGammaObjects + * @return ElectronPhotonVariableCorrectionBase::EGammaObjects, returns ElectronPhotonVariableCorrectionBase::EGammaObjects::Failure if no match found + */ + ElectronPhotonVariableCorrectionBase::EGammaObjects stringToEGammaObject( const std::string& input ) const; + + /** @brief Check if the photon which was passed to the tool has the correct type, if only (un)converted photons were requested to be corrected + * @param photon The photon which was passed to the tool to be corrected + */ + bool passedCorrectPhotonType(const xAOD::Photon& photon) const; + + /** @brief Check if the value passed is equal to one of the values passed via the UncorrectedDiscontinuites flag and should thus not be corrected + * @param value The value which should be checked + */ + bool isEqualToUncorrectedDiscontinuity(const float value) const; + + /** @brief Get the relevant information for a correction function parameter from the given conf file + * @param env The given TEnv, which is used to read out the current conf file + * @param parameter_number The parameter number with respect to the m_correctionFunctionTF1 + * @param type The type of the respective parameter with this parameter_number + * @details The relevant information for the parameter with the given parameter number and given parameter type is retrieved from the current configuration + * file, and saved in the according member variables of the class, i.e. m_graphCopies, m_binValues, m_etaBins, m_ptBins + */ + const StatusCode getParameterInformationFromConf(TEnv& env, const int parameter_number, const ElectronPhotonVariableCorrectionBase::parameterType type); + + /** @brief Get the actual parameters entering the correction TF1 for the current e/y object to be corrected + * @param properties The vector where the values of the correction TF1 parameters will be saved in + * @param pt The pT of the current e/y object to be corrected + * @param absEta The absolute eta of the current e/y object to be corrected + * @details As every electron/photon has different values of pT/eta, the correction function must be adapted accordingly for every e/y. The according values of + * each of the correction function parameters are updated with this function. + */ + const StatusCode getCorrectionParameters(std::vector<float>& properties, const float pt, const float absEta) const; + + /** @brief Get the correction function parameter value if its type is eta- or pT-binned + * @param return_parameter_value The respective correction function parameter value is saved in this parameter + * @param evalPoint pT or eta evaluation point - i.e., the eta or pT value of the current e/y object to be corrected. Used to find the correct eta/pT bin to use for the correction + * @param binning The eta or pT binning + * @param parameter_number The number of the parameter with respect to the correction TF1. Needed in order to retrieve the correct values matching this parameter. + */ + const StatusCode get1DBinnedParameter(float& return_parameter_value, const float evalPoint, const std::vector<float>& binning, const int parameter_number) const; + + /** @brief Get the correction function parameter value if its type is eta- and pT-binned + * @param return_parameter_value The respective correction function parameter value is saved in this parameter + * @param etaEvalPoint eta evaluation point - i.e., the eta value of the current e/y object to be corrected. Used to find the correct eta bin to use for the correction + * @param ptEvalPoint pT evaluation point - i.e., the pT value of the current e/y object to be corrected. Used to find the correct pT bin to use for the correction + * @param parameter_number The number of the parameter with respect to the correction TF1. Needed in order to retrieve the correct values matching this parameter. + */ + const StatusCode get2DBinnedParameter(float& return_parameter_value, const float etaEvalPoint, const float ptEvalPoint, const int parameter_number) const; + + /** @brief Find the bin number in which the evalPoint is located in the binning binning. + * @param return_bin The respective bin number is saved in this parameter + * @param evalPoint The evaluation point for which the bin number should be found + * @param binning The binning which should be evaluated + */ + const StatusCode findBin(int& return_bin, const float evalPoint, const std::vector<float>& binning) const; + + /** @brief Return the interpolation flag of parameter parameter_number as a boolean. + * @param parameter_number Number of the parameter for which the interpolation flag should be checked + * @return The interpolation flag of parameter parameter_number (boolean) + */ + bool getInterpolationFlag(const int parameter_number) const; + + /** @brief Given a point x, approximates the value via linear interpolation based on the two nearest bin centers. Re-implementation of Double_t TH1::Interpolate( Double_t x) const. + * @param return_parameter_value The interpolated parameter value is saved in this parameter + * @param evalPoint The point for which the interpolation should be done + * @param bin The bin in which the evalPoint is located with respect to binning + * @param binning The binning based on which the interpolation should be done + * @param binValues The bin values according to the binning given in binning + */ + const StatusCode interpolate(float& return_parameter_value, const float evalPoint, const int bin, const std::vector<float>& binning, const std::vector<float>& binValues) const; + + /** @brief Get the bin center of a bin bin_int using the binning binning + * @param return_bin_center The bin center is saved in this parameter + * @param binning The binning which should be used to find the bin centers + * @param bin_int The bin for which the bin center should be found + */ + const StatusCode getBinCenter(float& return_bin_center, const std::vector<float>& binning, const int bin_int) const; + + /** @brief Returns the linearly intrpolated value of value given the bin centers and bin values + * @param value The x-value at which the interpolation should be done + * @param left_bin_center The x-value of the left bin center used for the interpolation + * @param left_bin_value The y-value of the left bin at the left bin center + * @param right_bin_center The x-value of the right bin center used for the interpolation + * @param right_bin_value The y-value of the right bin at the right bin center + */ + float interpolate_function(const float value, const float left_bin_center, const float left_bin_value, const float right_bin_center, const float right_bin_value) const; + + /** @brief Get the events energy density from the eventShape + * @param value The respective correction function parameter value is saved in this parameter + * @param eventShapeContainer The name of the respective event shape container to use + */ + const StatusCode getDensity(float& value, const std::string& eventShapeContainer) const; + + /** @brief Get the e/y kinematic properties + * @param egamma_object The e/y object to get the kinemativ properties of + * @param pt The pT value is saved in this parameter + * @param absEta The absolute eta value is saved in this parameter + * @details As every electron/photon has different values of pT/eta, the correction function must be adapted accordingly for every e/y. The according values of + * eta and pt are updated with this function. + */ + const StatusCode getKinematicProperties(const xAOD::Egamma& egamma_object, float& pt, float& absEta) const; + + /** @brief Actual application of the correction to the variable + * @param return_corrected_variable The corrected variable value is saved in this parameter + * @param original_variable The original value of the corrected variable + * @param properties The vector containing the correction TF1 parameters so the correction TF1 can be set for the respective e/y object + */ + const StatusCode correct(float& return_corrected_variable, const float original_variable, std::vector<float>& properties) const; + +}; //end class ElectronPhotonVariableCorrectionBase + +#endif diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionDict.h b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionDict.h new file mode 100644 index 0000000000000000000000000000000000000000..154b6df66e155098ad2a1da8531df7f023af2403 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionDict.h @@ -0,0 +1,15 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +// Dear emacs, this is -*-c++-*- + +#ifndef ELECTRONPHOTONVARIABLECORRECTIONDICT_H +#define ELECTRONPHOTONVARIABLECORRECTIONDICT_H + + +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h" +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h" + + +#endif diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h new file mode 100644 index 0000000000000000000000000000000000000000..82551565881201042df97ac675597c2c8c7f7972 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h @@ -0,0 +1,131 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ElectronPhotonVariableCorrectionTool_H +#define ElectronPhotonVariableCorrectionTool_H + +/** @class ElectronPhotonVariableCorrectionTool + * @brief Tool to correct electron and photon MC variables using the ElectronPhotonVariableCorrectionBase class. + * @details For a detailed documentation, please refer to the [gitlab readme](https://gitlab.cern.ch/atlas/athena/-/blob/21.2/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/README.md) + * + * @author Nils Gillwald (DESY) nils.gillwald@desy.de + * @date February 2020 +*/ + +//ATLAS includes +#include "AsgTools/AsgTool.h" +#include "EgammaAnalysisInterfaces/IElectronPhotonShowerShapeFudgeTool.h" +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h" + +//EDM includes +#include "xAODEgamma/Electron.h" +#include "xAODEgamma/Photon.h" + +// =========================================================================== +// Class ElectronPhotonVariableCorrectionTool +// =========================================================================== + +class ElectronPhotonVariableCorrectionTool : public asg::AsgTool, virtual public IElectronPhotonShowerShapeFudgeTool +{ + +/// Declare the interface that the class provides +ASG_TOOL_CLASS(ElectronPhotonVariableCorrectionTool, IElectronPhotonShowerShapeFudgeTool) + +public: + /** @brief Standard constructor + * @param myname Internal name of the tool instance, so they can be distinguished + */ + ElectronPhotonVariableCorrectionTool(const std::string& myname); + //! @brief Standard destructor + ~ElectronPhotonVariableCorrectionTool() {}; + + /** @brief Initialize the tool instance + * @details Reads out the configuration file set via the setProperty(ConfigFile, "/path/to/conf/file.conf") function + * and sets up the tool instance accordingly + */ + virtual StatusCode initialize() override; + + /** @brief Apply the correction given in the conf file to the passed photon + * @param photon The photon which should be corrected + * @details Pass a photon which should be corrected by the tool to this function. The tool gets the original photon variable value, + * corrects it according to the configuration file given, and overwrites the original value of the variable. The original variable + * is saved as <variable_name>_original. Note that the variable to be corrected should be an auxiliary variable, and that the photon + * or photon container, respectively, must not be const if you want to use this function. If you are running on const objects, please + * make a deep / shallow copy of the container or use the correctedCopy function of this class. + */ + virtual const CP::CorrectionCode applyCorrection( xAOD::Photon& photon ) const override; + + /** @brief Apply the correction given in the conf file to the passed electron + * @param electron The electron which should be corrected + * @details Pass an electron which should be corrected by the tool to this function. The tool gets the original electron variable value, + * corrects it according to the configuration file given, and overwrites the original value of the variable. The original variable + * is saved as <variable_name>_original. Note that the variable to be corrected should be an auxiliary variable, and that the electron + * or electron container, respectively, must not be const if you want to use this function. If you are running on const objects, please + * make a deep / shallow copy of the container or use the correctedCopy function of this class. + */ + virtual const CP::CorrectionCode applyCorrection( xAOD::Electron& electron ) const override; + + /** @brief Make a corrected copy of the passed photon according to the given conf file + * @param in_photon The original photon of which a corrected copy should be made + * @param out_photon Empty new photon where the corrected copy will be stored in + * @details Makes a corrected copy of the passed photon object according to the correction provided in the configuration file. The in_photon + * is copied to the out_photon, which is then passed to the applyCorrection function. + */ + virtual const CP::CorrectionCode correctedCopy( const xAOD::Photon& in_photon, xAOD::Photon*& out_photon ) const override; + + /** @brief Make a corrected copy of the passed electron according to the given conf file + * @param in_electron The original electron of which a corrected copy should be made + * @param out_electron Empty new electron where the corrected copy will be stored in + * @details Makes a corrected copy of the passed electron object according to the correction provided in the configuration file. The in_electron + * is copied to the out_electron, which is then passed to the applyCorrection function. + */ + virtual const CP::CorrectionCode correctedCopy( const xAOD::Electron& in_electron, xAOD::Electron*& out_electron) const override; + +private: + //! @brief The configuration file for the tool + std::string m_configFile; + //! @brief The configuration files for the base class - converted photons + std::vector<std::string> m_convertedPhotonConfFiles; + //! @brief The configuration files for the base class - unconverted photons + std::vector<std::string> m_unconvertedPhotonConfFiles; + //! @brief The configuration files for the base class - electrons + std::vector<std::string> m_electronConfFiles; + //! @brief The base class instances for single variable correction - converted photons + std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>> m_convertedPhotonTools; + //! @brief The base class instances for single variable correction - unconverted photons + std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>> m_unconvertedPhotonTools; + //! @brief The base class instances for single variable correction - electrons + std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>> m_electronTools; + + //! @brief Configure, initialize, and check the base class instances saved in the m_*Tools members + const StatusCode initializeCorrectionTools(); + + /** @brief Initialize the base clase instances saved in the m_*Tools members + * @param name The name of the tool holder which should be included in the base class instances names + * @param confFiles The conf files for the base class instances which should be created and initialized + * @param toolHolder The tool holder where the base class instances are created + */ + const StatusCode initializeTools( const std::string& name, const std::vector<std::string>& confFiles, std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>>& toolHolder ); + + /** @brief Get current name of the variable to be corrected by the current base class instance, in order to coherently name the base class instances + * @param variableName The current correction variable name is saved in this parameter + * @param confFile The respective conf file of the current base class instance to be read-out + */ + const StatusCode getCorrectionVariableName( std::string &variableName, const std::string& confFile ) const; + + /** @brief Locate all config files for the base class instances passed to the tool in m_configFile in the ATLAS file system + * @param confFiles The conf files which shoud be used for the setup and initialization of the base class instances + */ + const StatusCode findAllConfigFiles( std::vector<std::string>& confFiles ); + + /** @brief Check if the ApplyTo flag from the base class conf file does match with the container type spcified in the tool conf file + * @param confFiles The conf files used to initialze the base class instances now stored in the respective toolHolder + * @param toolHolder The base class instances which were initialized using the respective confFiles + * @param toolHolderType The type of base class instances which should be stored in the respective toolHolder with respect to the tespectif confFiles + */ + const StatusCode applyToFlagMatchesToolHolder( const std::vector<std::string>& confFiles, const std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>>& toolHolder, ElectronPhotonVariableCorrectionBase::EGammaObjects toolHolderType ); + +}; //end class ElectronPhotonVariableCorrectionTool + +#endif diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/selection.xml b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/selection.xml new file mode 100644 index 0000000000000000000000000000000000000000..6d1a76b00e023a5ee9bd5da216bd80b4338264aa --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/EGammaVariableCorrection/selection.xml @@ -0,0 +1,3 @@ +<lcgdict> + <class name="ElectronPhotonVariableCorrectionTool" /> +</lcgdict> diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionBase.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionBase.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b528b91d0c0586109bc8f7e4a034216edd604ca7 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionBase.cxx @@ -0,0 +1,899 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + @class ElectronPhotonVariableCorrectionBase + @brief Tool to correct electron and photon MC variables. + + @author Nils Gillwald (DESY) nils.gillwald@desy.de + @date February 2020 +**/ + +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h" +#include "EgammaAnalysisHelpers/AsgEGammaConfigHelper.h" +#include "xAODEventShape/EventShape.h" + +// EDM includes +#include "xAODEgamma/EgammaEnums.h" +#include "TEnv.h" + +//Root includes +#include "TF1.h" +#include "TGraph.h" +#include "TFile.h" + +//allow advanced math for the TF1 +#include "TMath.h" + +#include "PathResolver/PathResolver.h" + +// =========================================================================== +// Standard Constructor +// =========================================================================== +ElectronPhotonVariableCorrectionBase::ElectronPhotonVariableCorrectionBase(const std::string& myname) : + AsgTool(myname) +{ + //declare the needed properties + declareProperty("ConfigFile",m_configFile="","The configuration file to use"); +} + +// =========================================================================== +// Initialize +// =========================================================================== +StatusCode ElectronPhotonVariableCorrectionBase::initialize() +{ + // Locate configuration file, abort if not found + std::string configFile; + if (!m_configFile.empty()) + { + configFile = PathResolverFindCalibFile(m_configFile); + if (configFile == "") + { + ATH_MSG_ERROR("Could not locate configuration file " << m_configFile); + return StatusCode::FAILURE; + } + ATH_MSG_DEBUG("Use configuration file " << m_configFile << " " << configFile); + } + else + { + ATH_MSG_ERROR("Config file string is empty. Please provide a config file to the tool."); + return StatusCode::FAILURE; + } + + // retrieve properties from configuration file, using TEnv class + TEnv env(configFile.c_str()); + // Send warning if duplicates found in conf file + env.IgnoreDuplicates(false); + + //retrieve variable to correct + if (env.Lookup("Variable")) + { + m_correctionVariable = env.GetValue("Variable",""); + } + else + { + ATH_MSG_ERROR("Correction variable is empty or not in configuration file."); + return StatusCode::FAILURE; + } + + // initialize the variable aux element accessors + // variable to be corrected + m_variableToCorrect = std::make_unique<SG::AuxElement::Accessor<float>>(m_correctionVariable); + // save original value under different name + m_originalVariable = std::make_unique<SG::AuxElement::Accessor<float>>(m_correctionVariable + "_original"); + + // Get the function used to correct the variable + if (env.Lookup("Function")) + { + m_correctionFunctionString = env.GetValue("Function","failure"); + // initialize the actual correction function TF1 + m_correctionFunctionTF1 = std::make_unique<TF1>(TF1("m_correctionFunctionTF1",m_correctionFunctionString.c_str())); + } + else + { + ATH_MSG_ERROR("Correction function is empty or not in configuration file."); + return StatusCode::FAILURE; + } + + // Get the number of parameters used in the correction function + if (env.Lookup("nFunctionParameters")) + { + m_numberOfFunctionParameters = env.GetValue("nFunctionParameters",-1); + } + else + { + ATH_MSG_ERROR("You did not specify the number of parameters in the correction function."); + return StatusCode::FAILURE; + } + + // resize all vectors used for getting parameters to m_numberOfFunctionParameters + m_ParameterTypeVector.resize(m_numberOfFunctionParameters); + m_binValues.resize(m_numberOfFunctionParameters); + m_graphCopies.resize(m_numberOfFunctionParameters); + m_interpolatePtFlags.resize(m_numberOfFunctionParameters); + + // Save the type of all parameters in the correction function (assuming m_numberOfFunctionParameters parameters) + for (int parameter_itr = 0; parameter_itr < m_numberOfFunctionParameters; parameter_itr++) + { + // if the parameter #parameter_itr is in the conf file, save its type + TString parameterType = TString::Format("Parameter%dType",parameter_itr); + if (env.Lookup(parameterType)) + { + // convert string to ParameterType, fail if non-existing type + ElectronPhotonVariableCorrectionBase::parameterType type = stringToParameterType(env.GetValue(parameterType.Data(),"")); + if( type == ElectronPhotonVariableCorrectionBase::parameterType::Failure ) + { + ATH_MSG_ERROR("Parameter " << parameter_itr << " read-in failed, not an allowed parameter type."); + return StatusCode::FAILURE; + } + // save type, get according type information and save it + m_ParameterTypeVector.at(parameter_itr) = type; + ATH_CHECK(getParameterInformationFromConf(env, parameter_itr, type)); + } + // else fail + else + { + ATH_MSG_ERROR("Did not find Parameter" << parameter_itr << ", although you specified there were " << m_numberOfFunctionParameters << " parameters for the correction function."); + return StatusCode::FAILURE; + } + } // end loop over all function parameters + + // check to which EGamma object the conf file should be applied to, if flag is set + if (env.Lookup("ApplyTo")) + { + std::string applyToObjectsFlag = env.GetValue("ApplyTo","Failure"); + m_applyToObjects = stringToEGammaObject(applyToObjectsFlag); + // fail if not passed a proper type + if (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::Failure) + { + ATH_MSG_ERROR("You did not correctly specify the object type in the ApplyTo flag."); + return StatusCode::FAILURE; + } + } + // else fail + else + { + ATH_MSG_ERROR("You did not specify to which objects this conf file should be applied to (ApplyTo)."); + return StatusCode::FAILURE; + } + + // check if there are any (discrete) values which should be left uncorrected + if (env.Lookup("UncorrectedDiscontinuities")) + { + m_uncorrectedDiscontinuities = AsgConfigHelper::HelperFloat("UncorrectedDiscontinuities", env); + // if flag is given, but no values, fail + if (m_uncorrectedDiscontinuities.size() < 1) + { + ATH_MSG_ERROR("Did not find any discontinuities to not correct, despite finding the flag UncorrectedDiscontinuities."); + return StatusCode::FAILURE; + } + } + + //everything worked out, so + return StatusCode::SUCCESS; +} + +// =========================================================================== +// Application of correction +// =========================================================================== +const CP::CorrectionCode ElectronPhotonVariableCorrectionBase::applyCorrection(xAOD::Photon& photon) const +{ + // check if we should only deal with converted / unconverted photons + if (!passedCorrectPhotonType(photon)) + { + ATH_MSG_ERROR("You specified in the conf file that the tool should only be used for (un-)converted photons, but passed the other conversion type."); + return CP::CorrectionCode::Error; + } + + // From the object, get the variable value according to the variable from the conf file + // if variable not found, fail + float original_variable = 0.; + if( m_variableToCorrect->isAvailable(photon) ) + { + original_variable = (*m_variableToCorrect)(photon); + //Save the original value to the photon under different name + (*m_originalVariable)(photon) = original_variable; + // check if tool should skip correcting this variable, as it's from some discontinuity + if (isEqualToUncorrectedDiscontinuity(original_variable)) + { + return CP::CorrectionCode::Ok; + } + } + else + { + ATH_MSG_ERROR("The correction variable \"" << m_correctionVariable << "\" provided in the conf file is not available."); + return CP::CorrectionCode::Error; + } + + //declare objects needed to retrieve photon properties + std::vector<float> properties; //safe value of function parameter i at place i + properties.resize(m_numberOfFunctionParameters); + float absEta; //safe absolute value of eta of event + float pt; //safe pt of event + + //Get all the properties from the photon and fill them to properties, eta, pt + if (getKinematicProperties(photon, pt, absEta) != StatusCode::SUCCESS) + { + ATH_MSG_ERROR("Could not retrieve kinematic properties of this photon object."); + return CP::CorrectionCode::Error; + } + ATH_MSG_VERBOSE("Got the photon kinematics , pT = " << pt << " |eta| = " << absEta); + if (getCorrectionParameters(properties, pt, absEta) != StatusCode::SUCCESS) + { + ATH_MSG_ERROR("Could not get the correction parameters for this photon object."); + return CP::CorrectionCode::Error; + } + for (auto p : properties) + ATH_MSG_VERBOSE("prop " << p); + + // Apply the correction, write to the corrected AuxElement + correct((*m_variableToCorrect)(photon),original_variable, properties).ignore(); // ignore as will always return SUCCESS + + // everything worked out, so + return CP::CorrectionCode::Ok; +} + +const CP::CorrectionCode ElectronPhotonVariableCorrectionBase::applyCorrection(xAOD::Electron& electron) const +{ + if (!(m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allElectrons || m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects)) + { + ATH_MSG_ERROR("You want to correct electrons, but passed a conf file with ApplyTo flag not set for electrons. Are you using the correct conf file?"); + return CP::CorrectionCode::Error; + } + + // From the object, get the variable value according to the variable from the conf file + // if variable not found, fail + float original_variable = 0.; + if (m_variableToCorrect->isAvailable(electron)) + { + original_variable = (*m_variableToCorrect)(electron); + //Save the original value to the photon under different name + (*m_originalVariable)(electron) = original_variable; + // check if tool should skip correcting this variable, as it's from some discontinuity + if (isEqualToUncorrectedDiscontinuity(original_variable)) + { + return CP::CorrectionCode::Ok; + } + } + else + { + ATH_MSG_ERROR("The correction variable \"" << m_correctionVariable << "\" provided in the conf file is not available."); + return CP::CorrectionCode::Error; + } + + //declare objects needed to retrieve electron properties + std::vector<float> properties; //safe value of function parameter i at place i + properties.resize(m_numberOfFunctionParameters); + float absEta; //safe absolute value of eta of event + float pt; //safe pt of event + + //Get all the properties from the electron and fill them to properties, eta, pt + if (getKinematicProperties(electron, pt, absEta) != StatusCode::SUCCESS) + { + ATH_MSG_ERROR("Could not retrieve kinematic properties of this electron object."); + return CP::CorrectionCode::Error; + } + ATH_MSG_VERBOSE("Got the electron kinematics , pT = " << pt << " |eta| = " << absEta); + if (getCorrectionParameters(properties, pt, absEta) != StatusCode::SUCCESS) + { + ATH_MSG_ERROR("Could not get the correction parameters for this electron object."); + return CP::CorrectionCode::Error; + } + + for (auto p : properties) + ATH_MSG_VERBOSE("prop " << p); + + // Apply the correction, write to the corrected AuxElement + correct((*m_variableToCorrect)(electron),original_variable, properties).ignore(); // ignore as will always return SUCCESS + + // everything worked out, so + return CP::CorrectionCode::Ok; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::correct(float& return_corrected_variable, const float original_variable, std::vector<float>& properties) const +{ + // set the parameters of the correction function + for (unsigned int parameter_itr = 0; parameter_itr < properties.size(); parameter_itr++) + { + m_correctionFunctionTF1->SetParameter(parameter_itr,properties.at(parameter_itr)); + } + + ATH_MSG_VERBOSE("original var " << original_variable); + + // Calculate corrected value + return_corrected_variable = m_correctionFunctionTF1->Eval(original_variable); + + ATH_MSG_VERBOSE("corrected var " << return_corrected_variable); + + // everything worked out, so + return StatusCode::SUCCESS; +} + +// =========================================================================== +// Corrected Copies +// =========================================================================== +const CP::CorrectionCode ElectronPhotonVariableCorrectionBase::correctedCopy( const xAOD::Photon& in_photon, xAOD::Photon*& out_photon ) const +{ + out_photon = new xAOD::Photon(in_photon); + return applyCorrection(*out_photon); +} + +const CP::CorrectionCode ElectronPhotonVariableCorrectionBase::correctedCopy( const xAOD::Electron& in_electron, xAOD::Electron*& out_electron) const +{ + out_electron = new xAOD::Electron(in_electron); + return applyCorrection(*out_electron); +} + +// =========================================================================== +// Helper Functions +// =========================================================================== + +bool ElectronPhotonVariableCorrectionBase::isEqualToUncorrectedDiscontinuity(const float value) const +{ + // if no values set, return false as there is nothing to check + if (m_uncorrectedDiscontinuities.size() < 1) + { + return false; + } + + // check all discontinuities which where passed + for (unsigned int value_itr = 0; value_itr < m_uncorrectedDiscontinuities.size(); value_itr++) + { + if (value == m_uncorrectedDiscontinuities.at(value_itr)) + { + // if the value is equal to one of the discontinuities, no need to check further + return true; + } + } + // if we eer get here, the value was never equal to a discontinuity + return false; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::getKinematicProperties(const xAOD::Egamma& egamma_object, float& pt, float& absEta) const +{ + // just reteriving eta and pt is probably less expensive then checking if I need it and + // then retrieve it only if I actually need it + + // protection against bad clusters + const xAOD::CaloCluster* cluster = egamma_object.caloCluster(); + if ( cluster == nullptr ) + { + ATH_MSG_ERROR("EGamma object calorimeter cluster is NULL: Cluster " << cluster); + return StatusCode::FAILURE; + } + + // Fill variables + // eta position in second sampling + absEta = std::abs(cluster->etaBE(2)); + // transverse energy in calorimeter (using eta position in second sampling) + const double energy = cluster->e(); + double et = 0.; + if (absEta<999.) { + const double cosheta = cosh(absEta); + et = (cosheta != 0.) ? energy/cosheta : 0.; + } + pt = et; + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::getParameterInformationFromConf(TEnv& env, const int parameter_number, const ElectronPhotonVariableCorrectionBase::parameterType type) +{ + // don't want to write the same code multiple times, so set flags when to retrieve eta/pt bins + bool getEtaBins = false; + bool getPtBins = false; + // form strings according to which parameter to retrieve + TString filePathKey = TString::Format("Parameter%dFile",parameter_number); + TString graphNameKey = TString::Format("Parameter%dGraphName",parameter_number); + TString binValues = TString::Format("Parameter%dValues",parameter_number); + TString interpolate = TString::Format("Parameter%dInterpolate",parameter_number); + // helpers + TString filePath = ""; + TString graphName = ""; + + // according to the parameter type, retrieve the information from conf + if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaDependentTGraph || type == ElectronPhotonVariableCorrectionBase::parameterType::PtDependentTGraph) + { + // check if necessary information is in conf, else fail + if (env.Lookup(filePathKey)) + { + //get the path to the root file where the graph is saved + filePath = PathResolverFindCalibFile(env.GetValue(filePathKey.Data(),"")); + // fail if file not found + if (filePath == "") + { + ATH_MSG_ERROR("Could not locate Parameter" << parameter_number << " TGraph file."); + return StatusCode::FAILURE; + } + } + else + { + ATH_MSG_ERROR("Could not retrieve Parameter" << parameter_number << " file path."); + return StatusCode::FAILURE; + } + // check if necessary information is in conf, else fail + if (env.Lookup(graphNameKey)) + { + //get the name of the TGraph + graphName = env.GetValue(graphNameKey.Data(),""); + } + else + { + ATH_MSG_ERROR("Could not retrieve Parameter" << parameter_number << " graph name."); + return StatusCode::FAILURE; + } + // open file, if it works, try to find graph, get graph, store a copy, else warning + fail + std::unique_ptr<TFile> file (new TFile(filePath.Data(),"READ")); + // check if file is open - if open, get graph, if not, fail + if (file->IsOpen()) + { + // if graph exists, get it, else fail + if (file->Get(graphName)) + { + std::unique_ptr<TGraph> graph ((TGraph*)file->Get(graphName.Data())); + m_graphCopies.at(parameter_number) = (TGraph*)graph->Clone(); //Or use copy constructor? + file->Close(); + } + else + { + ATH_MSG_ERROR("Could not find TGraph " << graphName << " in file " << filePath); + return StatusCode::FAILURE; + } + } + else + { + ATH_MSG_ERROR("Could not open Parameter" << parameter_number << " TGraph file " << filePath.Data()); + return StatusCode::FAILURE; + } + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaBinned ) + { + //get eta binning later + getEtaBins = true; + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::PtBinned ) + { + //get pt binning later + getPtBins = true; + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaTimesPtBinned ) + { + //get eta and pt binning later + getEtaBins = true; + getPtBins = true; + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EventDensity ) + { + // nothing has to be retrieved, no additional parameters for EventDensity currently + return StatusCode::SUCCESS; + } + + // if needed and not already retrieved, get eta binning + if (getEtaBins && !m_retrievedEtaBinning) + { + // check if necessary information is in conf, else fail + if (env.Lookup("EtaBins")) + { + //get the eta binning (global!) + m_etaBins = AsgConfigHelper::HelperFloat("EtaBins",env); + //force that the low bin edges are given by the conf file, starting with 0 + if (m_etaBins.at(0) != 0.) + { + ATH_MSG_ERROR("Lowest bin edge given for parameter " << parameter_number << " is not 0. Please provide the lower bin edges of your correction binning in the conf file, starting with 0."); + return StatusCode::FAILURE; + } + // don't want to retrieve the same thing twice from conf + m_retrievedEtaBinning = true; + } + else + { + ATH_MSG_ERROR("Could not retrieve eta binning."); + return StatusCode::FAILURE; + } + } + // if needed and not already retrieved, get pt binning + if (getPtBins && !m_retrievedPtBinning) + { + // check if necessary information is in conf, else fail + if (env.Lookup("PtBins")) + { + //get the pt binning (global!) + m_ptBins = AsgConfigHelper::HelperFloat("PtBins",env); + //force that the low bin edges are given by the conf file, starting with 0 + if (m_ptBins.at(0) != 0.) + { + ATH_MSG_ERROR("Lowest bin edge given for parameter " << parameter_number << " is not 0. Please provide the lower bin edges of your correction binning in the conf file, starting with 0."); + return StatusCode::FAILURE; + } + // don't want to retrieve the same thing twice from conf + m_retrievedPtBinning = true; + } + else + { + ATH_MSG_ERROR("Could not retrieve pt binning."); + return StatusCode::FAILURE; + } + } + //check if interpolation should be done + if (getPtBins) + { + // default is false + m_interpolatePtFlags.at(parameter_number) = env.GetValue(interpolate.Data(),false); + } + if ( getEtaBins || getPtBins) + { + // check if necessary information is in conf, else fail + if (env.Lookup(binValues)) + { + //get the binned values of the eta/pt dependent parameter + m_binValues.at(parameter_number) = AsgConfigHelper::HelperFloat(binValues.Data(),env); + } + else + { + ATH_MSG_ERROR("Could not retrieve binned values."); + return StatusCode::FAILURE; + } + } + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::getCorrectionParameters(std::vector<float>& properties, const float pt, const float absEta) const +{ + // according to the parameter type, get the actual parameter going to the correction function + // for this, loop over the parameter type vector + for (unsigned int parameter_itr = 0; parameter_itr < m_ParameterTypeVector.size(); parameter_itr++) + { + ElectronPhotonVariableCorrectionBase::parameterType type = m_ParameterTypeVector.at(parameter_itr); + if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaDependentTGraph) + { + // evaluate TGraph at abs(eta) + properties.at(parameter_itr) = m_graphCopies.at(parameter_itr)->Eval(absEta); + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::PtDependentTGraph) + { + // evaluate TGraph at pt + properties.at(parameter_itr) = m_graphCopies.at(parameter_itr)->Eval(pt); + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaBinned) + { + // get value of correct eta bin + ATH_CHECK(get1DBinnedParameter(properties.at(parameter_itr),absEta,m_etaBins,parameter_itr)); + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::PtBinned) + { + // get value of correct pt bin + ATH_CHECK(get1DBinnedParameter(properties.at(parameter_itr),pt,m_ptBins,parameter_itr)); + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EtaTimesPtBinned) + { + // get value of correct eta x pt bin + ATH_CHECK(get2DBinnedParameter(properties.at(parameter_itr),absEta,pt,parameter_itr)); + } + else if (type == ElectronPhotonVariableCorrectionBase::parameterType::EventDensity) + { + // get event density + ATH_CHECK(getDensity(properties.at(parameter_itr), "TopoClusterIsoCentralEventShape")); + } + } + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::get1DBinnedParameter(float& return_parameter_value, const float evalPoint, const std::vector<float>& binning, const int parameter_number) const +{ + ANA_MSG_VERBOSE("EvalPoint: " << evalPoint); + // need to find the bin in which the evalPoint is + int bin = -1; + ATH_CHECK(findBin(bin, evalPoint, binning)); + + // calculate return value + // if interpolation flag is true, interpolate + if (getInterpolationFlag(parameter_number)) + { + ATH_CHECK(interpolate(return_parameter_value, evalPoint, bin, binning, m_binValues.at(parameter_number))); + } + else + { + return_parameter_value = m_binValues.at(parameter_number).at(bin); + } + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::get2DBinnedParameter(float& return_parameter_value, const float etaEvalPoint, const float ptEvalPoint, const int parameter_number) const +{ + //need to find eta bin, and need to find pt bin + //from this, calculate which parameter of the list is needed to be returned. + int etaBin = -1; + int ptBin = -1; + + ATH_MSG_VERBOSE("eta = " << etaEvalPoint << " pt = " << ptEvalPoint); + + ATH_CHECK(findBin(etaBin, etaEvalPoint, m_etaBins)); + ATH_CHECK(findBin(ptBin, ptEvalPoint, m_ptBins)); + + // get the corresponding pt x eta bin found + /* Note: Assuming that the values are binned in pt x eta in the conf file: + * eta bin 0 | eta bin 1 | eta bin 2 | eta bin 3 | eta bin 4 | etc. + * pt bin 0 0 1 2 3 4 + * pt bin 1 5 6 7 8 9 + * pt bin 2 10 11 12 13 14 + * pt bin 3 15 16 17 18 19 + * pt bin 4 20 21 22 23 24 + * etc. + * the correct parameter is saved in the vector at (number of eta bins) * ptBinNumber + etaBinNumber + * */ + + ATH_MSG_VERBOSE(" eta bin " << etaBin << " pT bin " << ptBin); + + int bin_number_in_bin_values = m_etaBins.size() * ptBin + etaBin; + + // calculate return value + // if interpolation flag is true, interpolate + if (getInterpolationFlag(parameter_number)) + { + // create the vector of correction values binned in pT at the found eta bin + std::vector<float> tmp_binValuesAtEtaSlice; + // from the full binning vector, need to cut one of the columns + for (unsigned int binValue_itr = etaBin; binValue_itr < m_binValues.at(parameter_number).size(); binValue_itr += m_etaBins.size()) + { + tmp_binValuesAtEtaSlice.push_back(m_binValues.at(parameter_number).at(binValue_itr)); + } + // interpolate using the above slice of the correction values + ATH_CHECK(interpolate(return_parameter_value, ptEvalPoint, ptBin, m_ptBins, tmp_binValuesAtEtaSlice)); + } + else + { + return_parameter_value = m_binValues.at(parameter_number).at(bin_number_in_bin_values); + } + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::findBin(int& return_bin, const float evalPoint, const std::vector<float>& binning) const +{ + // need to find the bin in which the evalPoint is + return_bin = -1; + // if the evalPoint is < 0, something is very wrong + if (evalPoint < binning.at(0)) + { + // I use it for eta and pT, so error message is tailored to that... + ATH_MSG_ERROR("Abs(Eta) or pT of object is smaller than 0."); + return StatusCode::FAILURE; + } + // loop over bin boundaries + //run only up to binning.size()-1, as running to binning.size() will get a seg fault for the boundary check + for (unsigned int bin_itr = 0; bin_itr < binning.size()-1; bin_itr++) + { + // if the evaluation point is within the checked bins boundaries, this is the value we want + if (evalPoint > binning.at(bin_itr) && evalPoint < binning.at(bin_itr + 1)) + { + // we found the according bin, so return the according value + return_bin = bin_itr; + // we can stop now + break; + } + } + //if this point is ever reached and bin == -1, the evaluation point is larger then the largest lowest bin edge + //if that's the case, return the value for the last bin + //the -1 is because the parameter numbering in a vector starts at 0 + if (return_bin == -1) + { + return_bin = m_binValues.size()-1; + } + + // everythin went fine, so + return StatusCode::SUCCESS; +} + +bool ElectronPhotonVariableCorrectionBase::getInterpolationFlag(const int parameter_number) const +{ + bool do_interpolation = false; + // get parameter number type + ElectronPhotonVariableCorrectionBase::parameterType type = m_ParameterTypeVector.at(parameter_number); + // check if parameter has a type for which we need to check if interpolation is wanted + if (type == ElectronPhotonVariableCorrectionBase::parameterType::PtBinned || type == ElectronPhotonVariableCorrectionBase::parameterType::EtaTimesPtBinned) + { + // get the flag + do_interpolation = m_interpolatePtFlags.at(parameter_number); + } + return do_interpolation; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::interpolate(float& return_parameter_value, const float evalPoint, const int bin, const std::vector<float>& binning, const std::vector<float>& binValues) const +{ + // check if passed binning is consistent + if (binning.size() != binValues.size()) + { + ATH_MSG_ERROR("Binning and bin values have different sizes."); + return StatusCode::FAILURE; + } + + // if evalPoint is left to the leftmost bin center, return the leftmost bin center without interpolation + float leftmost_bin_center = 0; + ANA_CHECK(getBinCenter(leftmost_bin_center, binning, 0)); + if (evalPoint <= leftmost_bin_center) + { + return_parameter_value = binValues.at(0); + return StatusCode::SUCCESS; + } + + // if evalPoint is right to the rightmost bin center, return the rightmost bin center without interpolation + float rightmost_bin_center = 0; + ANA_CHECK(getBinCenter(rightmost_bin_center, binning, binning.size()-1)) + if (evalPoint >= rightmost_bin_center) + { + return_parameter_value = binValues.at(binValues.size()-1); + return StatusCode::SUCCESS; + } + + float left_bin_center = 0; + float right_bin_center = 0; + float left_bin_value = 0; + float right_bin_value = 0; + float current_bin_center = 0; + ANA_CHECK(getBinCenter(current_bin_center, binning, bin)); + + // else interpolate using next left or next right bin + if (evalPoint <= current_bin_center) + { + //interpolate left + ANA_CHECK(getBinCenter(left_bin_center, binning, bin-1)); + ANA_CHECK(getBinCenter(right_bin_center, binning, bin)); + left_bin_value = binValues.at(bin-1); + right_bin_value = binValues.at(bin); + } + else // evalPoint is right from bin center + { + //interpolate right + ANA_CHECK(getBinCenter(left_bin_center, binning, bin)); + ANA_CHECK(getBinCenter(right_bin_center, binning, bin+1)); + left_bin_value = binValues.at(bin); + right_bin_value = binValues.at(bin+1); + } + ATH_MSG_VERBOSE("bin centers : " << left_bin_center << " " << right_bin_center << " current : " << current_bin_center << " values : " << left_bin_value << " " << right_bin_value); + + // calculate return value + return_parameter_value = interpolate_function(evalPoint, left_bin_center, left_bin_value, right_bin_center, right_bin_value); + + // everything went fine, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionBase::getBinCenter(float& return_bin_center, const std::vector<float>& binning, const int bin_int) const +{ + if (bin_int < 0) + { + ATH_MSG_ERROR("Bin number must be a non-negative integer."); + return StatusCode::FAILURE; + } + + // implicitly convert to loong unsigend int for comparisons, to get rid of compiler warnings resulting from comparisons of int and unsigned int + long unsigned int bin = bin_int; + + if (bin >= binning.size()) + { + ATH_MSG_ERROR("The requested bin is out of range of the passed binning."); + return StatusCode::FAILURE; + } + + // need special treatment for rightmost bin center: + // it goes up to infinity...assume it's as broad as the + // next to the rightmost bin for the interpolation + if (bin == binning.size()-1) + { + //calculate the width of the next to rightmost bin + float bin_width = binning.at(bin) - binning.at(bin-1); + return_bin_center = binning.at(bin) + 0.5 * bin_width; + } + // normal calculation + else + { + return_bin_center = 0.5 * (binning.at(bin) + binning.at(bin+1)); + } + + //everything went fine, so + return StatusCode::SUCCESS; +} + +float ElectronPhotonVariableCorrectionBase::interpolate_function(const float value, const float left_bin_center, const float left_bin_value, const float right_bin_center, const float right_bin_value) const +{ + return left_bin_value + (value - left_bin_center) * (right_bin_value - left_bin_value) / (right_bin_center - left_bin_center); +} + +const StatusCode ElectronPhotonVariableCorrectionBase::getDensity(float& value, const std::string& eventShapeContainer) const +{ + // retrieve the event shape container + const xAOD::EventShape* evtShape = nullptr; + if(evtStore()->retrieve(evtShape, eventShapeContainer).isFailure()){ + ATH_MSG_ERROR("Cannot retrieve density container " << eventShapeContainer); + return StatusCode::FAILURE; + } + // get the density from the container + value = evtShape->getDensity(xAOD::EventShape::Density); + return StatusCode::SUCCESS; +} + +ElectronPhotonVariableCorrectionBase::parameterType ElectronPhotonVariableCorrectionBase::stringToParameterType( const std::string& input ) const +{ + // return parameter type according to string given in conf file + if( input == "EtaDependentTGraph") + { return ElectronPhotonVariableCorrectionBase::parameterType::EtaDependentTGraph; } + else if( input == "PtDependentTGraph") + { return ElectronPhotonVariableCorrectionBase::parameterType::PtDependentTGraph; } + else if( input == "EtaBinned") + { return ElectronPhotonVariableCorrectionBase::parameterType::EtaBinned; } + else if( input == "PtBinned") + { return ElectronPhotonVariableCorrectionBase::parameterType::PtBinned; } + else if( input == "EtaTimesPtBinned") + { return ElectronPhotonVariableCorrectionBase::parameterType::EtaTimesPtBinned; } + else if( input == "EventDensity") + { return ElectronPhotonVariableCorrectionBase::parameterType::EventDensity; } + else + { + // if not a proper type, return failure type - check and fail on this! + ATH_MSG_ERROR(input.c_str() << " is not an allowed parameter type."); + return ElectronPhotonVariableCorrectionBase::parameterType::Failure; + } +} + +ElectronPhotonVariableCorrectionBase::EGammaObjects ElectronPhotonVariableCorrectionBase::stringToEGammaObject( const std::string& input ) const +{ + // return object type which correction should be applied to + if( input == "unconvertedPhotons" ) + { return ElectronPhotonVariableCorrectionBase::EGammaObjects::unconvertedPhotons; } + else if( input == "convertedPhotons" ) + { return ElectronPhotonVariableCorrectionBase::EGammaObjects::convertedPhotons; } + else if( input == "allPhotons" ) + { return ElectronPhotonVariableCorrectionBase::EGammaObjects::allPhotons; } + else if( input == "allElectrons" ) + { return ElectronPhotonVariableCorrectionBase::EGammaObjects::allElectrons; } + else if( input == "allEGammaObjects" ) + { return ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects; } + else + { + // if not a proper object type, return failure type - check and fail on this! + ATH_MSG_ERROR(input.c_str() << " is not an allowed EGamma object type to apply corrections to."); + return ElectronPhotonVariableCorrectionBase::EGammaObjects::Failure; + } +} + +bool ElectronPhotonVariableCorrectionBase::passedCorrectPhotonType(const xAOD::Photon& photon) const +{ + // retrieve if photon is converted or unconverted + bool isConvertedPhoton = xAOD::EgammaHelpers::isConvertedPhoton(&photon); + bool isUnconvertedPhoton = !isConvertedPhoton; + + // check if conf file ApplyTo flag matches photon conversion type + return ((applyToConvertedPhotons() && isConvertedPhoton) || (applyToUnconvertedPhotons() && isUnconvertedPhoton)); +} + +bool ElectronPhotonVariableCorrectionBase::applyToConvertedPhotons() const +{ + bool applyToAllEGamma = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects); + bool applyToAllPhotons = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allPhotons); + bool applyToConvertedPhotons = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::convertedPhotons); + return (applyToAllEGamma || applyToAllPhotons || applyToConvertedPhotons); +} + +bool ElectronPhotonVariableCorrectionBase::applyToUnconvertedPhotons() const +{ + bool applyToAllEGamma = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects); + bool applyToAllPhotons = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allPhotons); + bool applyToUnconvertedPhotons = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::unconvertedPhotons); + return (applyToAllEGamma || applyToAllPhotons || applyToUnconvertedPhotons); +} + +bool ElectronPhotonVariableCorrectionBase::applyToElectrons() const +{ + bool applyToAllEGamma = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects); + bool applyToAllElectrons = (m_applyToObjects == ElectronPhotonVariableCorrectionBase::EGammaObjects::allElectrons); + return (applyToAllEGamma || applyToAllElectrons); +} diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionTool.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c57b4f12c4fb7ba04190702f30b624218c66aaba --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/Root/ElectronPhotonVariableCorrectionTool.cxx @@ -0,0 +1,257 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + @class ElectronPhotonVariableCorrectionTool + @brief Tool to correct electron and photon MC variables. + + @author Nils Gillwald (DESY) nils.gillwald@desy.de + @date February 2020 +**/ + +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h" +#include "EgammaAnalysisHelpers/AsgEGammaConfigHelper.h" +#include "PathResolver/PathResolver.h" +#include "TEnv.h" + +// =========================================================================== +// Standard Constructor +// =========================================================================== +ElectronPhotonVariableCorrectionTool::ElectronPhotonVariableCorrectionTool(const std::string& myname) : + AsgTool(myname) +{ + //declare the needed properties + declareProperty("ConfigFile",m_configFile="", "The configuration file to use"); +} + +StatusCode ElectronPhotonVariableCorrectionTool::initialize() +{ + // Locate configuration file, abort if not found + std::string configFile; + if (!m_configFile.empty()) + { + configFile = PathResolverFindCalibFile(m_configFile); + if (configFile == "") + { + ATH_MSG_ERROR("Could not locate configuration file " << m_configFile); + return StatusCode::FAILURE; + } + ATH_MSG_DEBUG("Use configuration file " << m_configFile); + } + else + { + ATH_MSG_ERROR("Config file string is empty. Please provide a config file to the tool."); + return StatusCode::FAILURE; + } + + // Retreive properties from configuration file, using TEnv class + TEnv env(configFile.c_str()); + // Send warning if duplicates found in conf file + env.IgnoreDuplicates(false); + + // retrieve different object type conf files + if (env.Lookup("ElectronConfigs")) + { + m_electronConfFiles = AsgConfigHelper::HelperString("ElectronConfigs",env); + } + if (env.Lookup("ConvertedPhotonConfigs")) + { + m_convertedPhotonConfFiles = AsgConfigHelper::HelperString("ConvertedPhotonConfigs",env); + } + if (env.Lookup("UnconvertedPhotonConfigs")) + { + m_unconvertedPhotonConfFiles = AsgConfigHelper::HelperString("UnconvertedPhotonConfigs",env); + } + + // check if any conf files were received + if (m_electronConfFiles.size() + m_convertedPhotonConfFiles.size() + m_unconvertedPhotonConfFiles.size() < 1) + { + ANA_MSG_ERROR("You did not provide any config files for the ElectronPhotonVariableCorrectionBase to the ElectronPhotonVariableCorrectionTool."); + return StatusCode::FAILURE; + } + + ATH_MSG_VERBOSE("number of files for the electron case : " << m_electronConfFiles.size()); + for (auto fN : m_electronConfFiles) + ATH_MSG_VERBOSE("file " << fN); + + // initialize the ElectronPhotonVariableCorrectionTools + ANA_CHECK(initializeCorrectionTools()); + + ANA_MSG_INFO("Initialized tool " << name()); + //everything worked out, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionTool::initializeCorrectionTools() +{ + //initialize all tools + ANA_CHECK(initializeTools("convertedPhotons", m_convertedPhotonConfFiles, m_convertedPhotonTools)); + ANA_CHECK(initializeTools("unconvertedPhotons", m_unconvertedPhotonConfFiles, m_unconvertedPhotonTools)); + ANA_CHECK(initializeTools("electrons", m_electronConfFiles, m_electronTools)); + // check if ApplyTo Flag matches with the tool holder + ANA_CHECK(applyToFlagMatchesToolHolder(m_convertedPhotonConfFiles, m_convertedPhotonTools, ElectronPhotonVariableCorrectionBase::EGammaObjects::convertedPhotons)); + ANA_CHECK(applyToFlagMatchesToolHolder(m_unconvertedPhotonConfFiles, m_unconvertedPhotonTools, ElectronPhotonVariableCorrectionBase::EGammaObjects::unconvertedPhotons)); + ANA_CHECK(applyToFlagMatchesToolHolder(m_electronConfFiles, m_electronTools, ElectronPhotonVariableCorrectionBase::EGammaObjects::allElectrons)); + + //everything worked out, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionTool::findAllConfigFiles( std::vector<std::string>& confFiles ) +{ + //loop over conf file vector, find conf file using path resolver + for( unsigned int confFile_itr = 0; confFile_itr < confFiles.size(); confFile_itr++ ) + { + std::string tmp_confFile = confFiles.at(confFile_itr); + confFiles.at(confFile_itr) = PathResolverFindCalibFile(confFiles.at(confFile_itr)); + if (confFiles.at(confFile_itr) == "") + { + ATH_MSG_ERROR("Could not locate configuration file " << tmp_confFile); + return StatusCode::FAILURE; + } + } + //everything worked out, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionTool::initializeTools( const std::string& name, const std::vector<std::string>& confFiles, std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>>& toolHolder ) +{ + // adapt size of toolHolder + toolHolder.resize(confFiles.size()); + // for each conf file, initialize one tool + for( unsigned int confFile_itr = 0; confFile_itr < confFiles.size(); confFile_itr++ ) + { + // name: supertool name + type name + variable name + std::string variable = ""; //get the name of the variable to be corrected + std::string confFileWithFullPath = PathResolverFindCalibFile(confFiles.at(confFile_itr)); + ANA_CHECK(getCorrectionVariableName(variable, confFileWithFullPath)); + TString toolname = TString::Format("%s_%s_%s", this->name().c_str(), name.c_str(), variable.c_str()); + ANA_MSG_DEBUG("Subtool name: " << toolname.Data()); + toolHolder.at(confFile_itr) = std::make_unique<ElectronPhotonVariableCorrectionBase>(toolname.Data()); + ANA_CHECK(toolHolder.at(confFile_itr)->setProperty("ConfigFile", confFiles.at(confFile_itr))); + ANA_CHECK(toolHolder.at(confFile_itr)->setProperty("OutputLevel", this->msg().level())); + ANA_CHECK(toolHolder.at(confFile_itr)->initialize()); + } + //everything worked out, so + return StatusCode::SUCCESS; +} + +const StatusCode ElectronPhotonVariableCorrectionTool::applyToFlagMatchesToolHolder( const std::vector<std::string>& confFiles, const std::vector<std::unique_ptr<ElectronPhotonVariableCorrectionBase>>& toolHolder, ElectronPhotonVariableCorrectionBase::EGammaObjects toolHolderType ) +{ + // loop over conf file holder + for (unsigned int tool_itr = 0; tool_itr < toolHolder.size(); tool_itr++) + { + // get ApplyTo flag + ElectronPhotonVariableCorrectionBase::EGammaObjects confFileType = toolHolder.at(tool_itr)->isAppliedTo(); + // skip all further tests if should be applied to all objects + if (confFileType == ElectronPhotonVariableCorrectionBase::EGammaObjects::allEGammaObjects) continue; + // continue if ApplyTo flag matches toolholder + if (toolHolderType == ElectronPhotonVariableCorrectionBase::EGammaObjects::convertedPhotons && toolHolder.at(tool_itr)->applyToConvertedPhotons()) continue; + if (toolHolderType == ElectronPhotonVariableCorrectionBase::EGammaObjects::unconvertedPhotons && toolHolder.at(tool_itr)->applyToUnconvertedPhotons()) continue; + if (toolHolderType == ElectronPhotonVariableCorrectionBase::EGammaObjects::allElectrons && toolHolder.at(tool_itr)->applyToElectrons()) continue; + // if this point is reached, something is wrong, so + ATH_MSG_ERROR("In conf " << confFiles.at(tool_itr) << ": The ApplyTo flag does not match with the container type from the conf file."); + return StatusCode::FAILURE; + } + //everything worked out, so + return StatusCode::SUCCESS; +} + +/* ======================================== + * Apply correction + * ======================================== */ + +const CP::CorrectionCode ElectronPhotonVariableCorrectionTool::applyCorrection( xAOD::Photon& photon ) const +{ + bool isConvertedPhoton = xAOD::EgammaHelpers::isConvertedPhoton(&photon); + + // check if need to run over converted or unconverted photons + if (isConvertedPhoton) + { + // correct variables on the converted photon + for (unsigned int convertedPhotonTool_itr = 0; convertedPhotonTool_itr < m_convertedPhotonTools.size(); convertedPhotonTool_itr++) + { + ATH_MSG_VERBOSE("Running tool " << m_convertedPhotonTools.at(convertedPhotonTool_itr)->name()); + if ((m_convertedPhotonTools.at(convertedPhotonTool_itr)->applyCorrection(photon)) != CP::CorrectionCode::Ok) + { + ATH_MSG_ERROR("Could not apply correction to converted photon object."); + return CP::CorrectionCode::Error; + } + } + } + else + { + // correct variables on the converted photon + for (unsigned int unconvertedPhotonTool_itr = 0; unconvertedPhotonTool_itr < m_unconvertedPhotonTools.size(); unconvertedPhotonTool_itr++) + { + ATH_MSG_VERBOSE("Running tool " << m_unconvertedPhotonTools.at(unconvertedPhotonTool_itr)->name()); + if ((m_unconvertedPhotonTools.at(unconvertedPhotonTool_itr)->applyCorrection(photon)) != CP::CorrectionCode::Ok) + { + ATH_MSG_ERROR("Could not apply correction to unconverted photon object."); + return CP::CorrectionCode::Error; + } + } + } + + // everything worked out, so + return CP::CorrectionCode::Ok; +} +const CP::CorrectionCode ElectronPhotonVariableCorrectionTool::applyCorrection( xAOD::Electron& electron ) const +{ + // correct variables on the electron + for (unsigned int electronTool_itr = 0; electronTool_itr < m_electronTools.size(); electronTool_itr++) + { + ATH_MSG_VERBOSE("Running tool " << m_electronTools.at(electronTool_itr)->name()); + if ((m_electronTools.at(electronTool_itr)->applyCorrection(electron)) != CP::CorrectionCode::Ok) + { + ATH_MSG_ERROR("Could not apply correction to electron object."); + return CP::CorrectionCode::Error; + } + } + // everything worked out, so + return CP::CorrectionCode::Ok; +} + +/* ======================================== + * Corrected Copy + * ======================================== */ + +const CP::CorrectionCode ElectronPhotonVariableCorrectionTool::correctedCopy( const xAOD::Photon& in_photon, xAOD::Photon*& out_photon ) const +{ + ATH_MSG_VERBOSE("Will correct photon " << &in_photon << " of pT, eta = " << in_photon.pt() << " " << in_photon.eta()); + out_photon = new xAOD::Photon(in_photon); + return applyCorrection(*out_photon); +} + +const CP::CorrectionCode ElectronPhotonVariableCorrectionTool::correctedCopy( const xAOD::Electron& in_electron, xAOD::Electron*& out_electron) const +{ + ATH_MSG_VERBOSE("Will correct electron " << &in_electron << " of pT, eta = " << in_electron.pt() << " " << in_electron.eta()); + out_electron = new xAOD::Electron(in_electron); + return applyCorrection(*out_electron); +} + +/* ======================================== + * Helper functions + * ======================================== */ + +const StatusCode ElectronPhotonVariableCorrectionTool::getCorrectionVariableName( std::string &variableName, const std::string& confFile ) const +{ + // Retreive properties from configuration file, using TEnv class + TEnv env(confFile.c_str()); + // Send warning if duplicates found in conf file + env.IgnoreDuplicates(false); + + // retrieve variable name + if (env.Lookup("Variable")) + { + variableName = env.GetValue("Variable", ""); + } + else + { + ATH_MSG_ERROR("In conf file " << confFile << ": Correction variable is empty or not in configuration file."); + return StatusCode::FAILURE; + } + //everything worked out, so + return StatusCode::SUCCESS; +} diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/EGammaVariableCorrectionTool_ExampleConf.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/EGammaVariableCorrectionTool_ExampleConf.conf new file mode 100644 index 0000000000000000000000000000000000000000..b99d311da8a205b76d7ef0936eb6858b909093e9 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/EGammaVariableCorrectionTool_ExampleConf.conf @@ -0,0 +1,10 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of a (un)converted photon or electron +# simultaneously in all to be corrected variables + +ElectronConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Eratio.conf; ++ElectronConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Rhad.conf +ConvertedPhotonConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Eratio.conf; ++ConvertedPhotonConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Rhad.conf +UnconvertedPhotonConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Eratio.conf; ++UnconvertedPhotonConfigs: EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Rhad.conf \ No newline at end of file diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Eratio.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Eratio.conf new file mode 100644 index 0000000000000000000000000000000000000000..f4f67885dfa9a4d01bab6c48012ece77e5372b5a --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Eratio.conf @@ -0,0 +1,20 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to photons + +Variable: Eratio +Function: [0] * x + [1] +nFunctionParameters: 2 + +ApplyTo: convertedPhotons + +UncorrectedDiscontinuities: 0.; 1. + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. + +EtaBins: 0.0; 1.5; 3. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Rhad.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Rhad.conf new file mode 100644 index 0000000000000000000000000000000000000000..36f75f5380fc6284b25f2f3539de45a803e19c63 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Rhad.conf @@ -0,0 +1,20 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to photons + +Variable: Rhad +Function: [0] * x + [1] +nFunctionParameters: 2 + +ApplyTo: convertedPhotons + +UncorrectedDiscontinuities: 0.; 1. + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. + +EtaBins: 0.0; 1.5; 3. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Eratio.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Eratio.conf new file mode 100644 index 0000000000000000000000000000000000000000..04e812b46ec864e0ef99a0dec9a704bbe63b599d --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Eratio.conf @@ -0,0 +1,31 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to electrons + +ApplyTo: allElectrons + +Variable: Eratio + +UncorrectedDiscontinuities: 0.; 1. + +Function: [0] + [1] * x + [2] * x**2 +nFunctionParameters: 3 + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. +Parameter0Interpolate: TRUE + +EtaBins: 0.0; 1.; 2. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. + +Parameter2Type: EtaTimesPtBinned +Parameter2Interpolate: TRUE +#pT \ eta >0.0 >1. >2. +Parameter2Values: 0.5; 1.; 1.5; #0 - 10 GeV# ++Parameter2Values: 2.; 2.5; 3; #10 - 50 GeV# ++Parameter2Values: 3.1; 3.2; 3.3; #50 - 100 GeV# ++Parameter2Values: 3.4; 3.5; 4. #> 100 GeV# + diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Rhad.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Rhad.conf new file mode 100644 index 0000000000000000000000000000000000000000..fe80f9e62c5d88bb320b85cc43559ad44b38aff5 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Rhad.conf @@ -0,0 +1,21 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to electrons + +ApplyTo: allElectrons + +Variable: Rhad + +UncorrectedDiscontinuities: 0.; 1. + +Function: [0] * x + [1] +nFunctionParameters: 2 + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. + +EtaBins: 0.0; 1.5; 3. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Eratio.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Eratio.conf new file mode 100644 index 0000000000000000000000000000000000000000..b0553a664c7813edb3db406d1d1d7fef4e36d174 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Eratio.conf @@ -0,0 +1,20 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to photons + +Variable: Eratio +Function: [0] * x + [1] +nFunctionParameters: 2 + +ApplyTo: unconvertedPhotons + +UncorrectedDiscontinuities: 0.; 1. + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. + +EtaBins: 0.0; 1.5; 3. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Rhad.conf b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Rhad.conf new file mode 100644 index 0000000000000000000000000000000000000000..f0b10f7fcaed5bba0b371853e7d159625dd0dd22 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/data/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Rhad.conf @@ -0,0 +1,20 @@ +# This is a ROOT configuration file to be read via TEnv +# Purpose is to steer the correction of an aux var via the apply tool + +# This is a configuration file to test the application of corrections to photons + +Variable: Rhad +Function: [0] * x + [1] +nFunctionParameters: 2 + +ApplyTo: unconvertedPhotons + +UncorrectedDiscontinuities: 0.; 1. + +PtBins: 0.0; 10000; 50000; 100000 +Parameter0Type: PtBinned +Parameter0Values: 1.; 2.; 3.; 4. + +EtaBins: 0.0; 1.5; 3. +Parameter1Type: EtaBinned +Parameter1Values: 1.; 2.; 3. diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/src/components/EGammaVariableCorrection_entries.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/src/components/EGammaVariableCorrection_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c1ffd7aa303934d6d66ef81ac8d09c0b84d0ddf3 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/src/components/EGammaVariableCorrection_entries.cxx @@ -0,0 +1,6 @@ +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool.h" + +DECLARE_COMPONENT( ElectronPhotonVariableCorrectionTool ) + + + diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionBase.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionBase.cxx new file mode 100644 index 0000000000000000000000000000000000000000..93a77a0055b34b3146f232df14b212ee5fd88be0 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionBase.cxx @@ -0,0 +1,243 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + * @brief Test if ElectronPhotonVariableCorrectionBase runs fine on electrons and (un-)converted photons + * @author Nils Gillwald (DESY) nils.gillwald@desy.de + * @date February 2020 + **/ + +// ROOT includes +#include "TFile.h" +#include "TString.h" + +//EDM includes +#include "xAODEgamma/ElectronContainer.h" +#include "xAODEgamma/PhotonContainer.h" +#include "xAODEgamma/PhotonAuxContainer.h" +#include "xAODEgamma/Electron.h" +#include "xAODEgamma/Photon.h" + +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h" +// +#include "AsgMessaging/AsgMessaging.h" + +////infrastructure includes +#ifdef ROOTCORE +#include "xAODRootAccess/TEvent.h" +//#include "xAODRootAccess/TStore.h" +#endif //ROOTCORE + +//main test code + +int main (int argc, char* argv[]) +{ + using namespace asg::msgUserCode; + // The application's name: + const char* APP_NAME = argv[0]; + + // Check if we received a file name: + if (argc < 2) + { + ANA_MSG_ERROR("No file name received!" ); + ANA_MSG_ERROR( " Usage: %s [xAOD file name] %d Num of events to process %d (0 photons , 1 electrons)"); + return EXIT_FAILURE; + } + + // Check if we want to process Electron/Photon shifts + bool isPhoton = false; + bool isElectron = false; + if (argc < 4) + { + ANA_MSG_INFO(APP_NAME << " By default looking at Photon" ); + isPhoton=true; + } + else + { + int argv1 = atoi(argv[3]); + if (argv1 == 0) + { + isPhoton = true; + ANA_MSG_INFO (APP_NAME << " We are looking at photon shifts" ); + } + else if (argv1 == 1) + { + isElectron = true; + ANA_MSG_INFO(APP_NAME << " We are looking at electron shifts"); + } + else + { + ANA_MSG_ERROR ("Usage: %s [xAOD file name] %d Num of events to process %d (0 photons , 1 electrons)"); + return EXIT_FAILURE; + } + } + + //open input file + const TString fileName = argv[1]; + ANA_MSG_INFO(APP_NAME << " Opening file: " << fileName.Data()); + std::unique_ptr<TFile> inputFile(TFile::Open(fileName, "READ")); //probably rather "READ"? + ANA_CHECK(inputFile.get()); + + // Create a TEvent object (persistent store) + xAOD::TEvent pers(xAOD::TEvent::kClassAccess); + // Create a TStore object (transient store) + xAOD::TStore trans; + ANA_CHECK(pers.readFrom(inputFile.get())); + // + ANA_MSG_INFO("Number of events in the file: " << pers.getEntries()); + + // Decide how many events to run over: + Long64_t entries = pers.getEntries(); + if (argc > 2) + { + const Long64_t userInputEntries = atoll(argv[2]); + if (userInputEntries < entries) + { + entries = userInputEntries; + ANA_MSG_INFO("Running over " << userInputEntries << " events."); + } + } + + // =============================================== + // Photon test + // =============================================== + if (isPhoton) + { + //initialise the tool + //converted photons + std::string configFilePathConverted = "ElectronPhotonShowerShapeFudgeTool/ElectronPhotonVariableCorrectionBase_ExampleConvertedPhotonConf_Eratio.conf"; + ElectronPhotonVariableCorrectionBase CorrectConvertedPhotonTool("CorrectConvertedPhotonTool"); + ANA_CHECK(CorrectConvertedPhotonTool.setProperty("ConfigFile",configFilePathConverted)); + ANA_CHECK(CorrectConvertedPhotonTool.initialize()); + //unconverted photons + std::string configFilePathUnconverted = "ElectronPhotonShowerShapeFudgeTool/ElectronPhotonVariableCorrectionBase_ExampleUnconvertedPhotonConf_Eratio.conf"; + ElectronPhotonVariableCorrectionBase CorrectUnconvertedPhotonTool("CorrectUnconvertedPhotonTool"); + ANA_CHECK(CorrectUnconvertedPhotonTool.setProperty("ConfigFile",configFilePathUnconverted)); + ANA_CHECK(CorrectUnconvertedPhotonTool.initialize()); + + std::string correctionVariable = CorrectConvertedPhotonTool.getCorrectionVariable(); + ANA_MSG_INFO("Correcting Variable: " << correctionVariable); + + //loop over the events + for (Long64_t entry = 0; entry < entries; entry++) + { + //get entry + pers.getEntry(entry); + ANA_MSG_INFO("============================"); + ANA_MSG_INFO("Event: " << entry); + + //get photon container + const xAOD::PhotonContainer* photons; + ANA_CHECK(pers.retrieve(photons, "Photons")); + + // Make a deep copy of the photon container + // Create the new container and its auxiliary store. + auto photons_copy = std::make_unique<xAOD::PhotonContainer>(); + auto photons_copy_aux = std::make_unique<xAOD::AuxContainerBase>(); + photons_copy->setStore (photons_copy_aux.get()); //< Connect the two + + //copy photons over + for (auto photon : *photons) { + // Copy this photon to the output container: + xAOD::Photon* photon_copy = new xAOD::Photon(); + photons_copy->push_back (photon_copy); // photon acquires the photon_copy auxstore + *photon_copy = *photon; // copies auxdata from one auxstore to the other + } + // end make deep copy + + //loop over deep copy of photon container + for (unsigned int photon_itr = 0; photon_itr < photons_copy->size(); photon_itr++) + { + ANA_MSG_INFO("---------------------------"); + ANA_MSG_INFO("Photon: " << photon_itr); + xAOD::Photon* photon = photons_copy->at(photon_itr); + + //apply correction + if (xAOD::EgammaHelpers::isConvertedPhoton(photon)) + { + ANA_MSG_INFO("Converted Photon."); + ANA_CHECK(CorrectConvertedPhotonTool.applyCorrection(*photon)); + } + else + { + ANA_MSG_INFO("Unconverted Photon."); + ANA_CHECK(CorrectUnconvertedPhotonTool.applyCorrection(*photon)); + } + + //get original and corrected value + SG::AuxElement::Accessor<float> VariableToCorrect(correctionVariable + "_original"); + SG::AuxElement::Accessor<float> CorrectedVariable(correctionVariable); + + //print results + ANA_MSG_INFO("Original value: " << VariableToCorrect(*photon)); + ANA_MSG_INFO("Corrected value: " << CorrectedVariable(*photon)); + + } // loop over deep copy of photon container + } // loop over events + } // if isPhoton + + // =============================================== + // Electron test + // =============================================== + if (isElectron) + { + std::string configFilePath = "ElectronPhotonShowerShapeFudgeTool/ElectronPhotonVariableCorrectionBase_ExampleElectronConf_Eratio.conf"; + ElectronPhotonVariableCorrectionBase CorrectElectronTool("CorrectElectronTool"); + ANA_CHECK(CorrectElectronTool.setProperty("ConfigFile",configFilePath)); + ANA_CHECK(CorrectElectronTool.initialize()); + + std::string correctionVariable = CorrectElectronTool.getCorrectionVariable(); + ANA_MSG_INFO("Correcting Variable: " << correctionVariable); + + //loop over the events + for (Long64_t entry = 0; entry < entries; entry++) + { + //get entry + pers.getEntry(entry); + ANA_MSG_INFO("============================"); + ANA_MSG_INFO("Event: " << entry); + + //get electron container + const xAOD::ElectronContainer* electrons; + ANA_CHECK(pers.retrieve(electrons, "Electrons")); + + // Make a deep copy of the electron container + // Create the new container and its auxiliary store. + auto electrons_copy = std::make_unique<xAOD::ElectronContainer>(); + auto electrons_copy_aux = std::make_unique<xAOD::AuxContainerBase>(); + electrons_copy->setStore (electrons_copy_aux.get()); //< Connect the two + + //copy electrons over + for (auto electron : *electrons) { + // Copy this electron to the output container: + xAOD::Electron* electron_copy = new xAOD::Electron(); + electrons_copy->push_back (electron_copy); // electron acquires the electron_copy auxstore + *electron_copy = *electron; // copies auxdata from one auxstore to the other + } + // end make deep copy + + //loop over deep copy of electron container + for (unsigned int electron_itr = 0; electron_itr < electrons_copy->size(); electron_itr++) + { + ANA_MSG_INFO("---------------------------"); + ANA_MSG_INFO("Electron: " << electron_itr); + xAOD::Electron* electron = electrons_copy->at(electron_itr); + + //apply correction + ANA_CHECK(CorrectElectronTool.applyCorrection(*electron)); + + //get original and corrected value + SG::AuxElement::Accessor<float> VariableToCorrect(correctionVariable + "_original"); + SG::AuxElement::Accessor<float> CorrectedVariable(correctionVariable); + + //print results + ANA_MSG_INFO("Original value: " << VariableToCorrect(*electron)); + ANA_MSG_INFO("Corrected value: " << CorrectedVariable(*electron)); + + } // loop over deep copy of electron container + } // loop over events + } // if isElectron + + return 0; +} diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..40338d003b75fb12cda2937067ade9eb42623ec2 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool.cxx @@ -0,0 +1,157 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + * @brief Test if ElectronPhotonVariableCorrectionTool runs fine on electrons and (un-)converted photons + * @author Nils Gillwald (DESY) nils.gillwald@desy.de + * @date February 2020 + **/ + +// ROOT includes +#include "TFile.h" +#include "TString.h" + +//EDM includes +#include "xAODEgamma/ElectronContainer.h" +#include "xAODEgamma/PhotonContainer.h" +#include "xAODEgamma/PhotonAuxContainer.h" +#include "xAODEgamma/Electron.h" +#include "xAODEgamma/Photon.h" + +#include "AsgTools/AnaToolHandle.h" +#include "EgammaAnalysisInterfaces/IElectronPhotonShowerShapeFudgeTool.h" + +#include "AsgMessaging/AsgMessaging.h" + +////infrastructure includes +#ifdef ROOTCORE +#include "xAODRootAccess/TEvent.h" +#include "xAODRootAccess/TStore.h" +#endif //ROOTCORE + +// main test code + +int main(int argc, char* argv[]) +{ + using namespace asg::msgUserCode; + + // The application's name: + const char* APP_NAME = argv[0]; + + // Check if we received a file name: + if (argc < 2) + { + ANA_MSG_ERROR("No file name received!" ); + ANA_MSG_ERROR( " Usage: %s [xAOD file name] %d Num of events to process %d (0 photons , 1 electrons)"); + return EXIT_FAILURE; + } + + //open input file + const TString fileName = argv[1]; + ANA_MSG_INFO(APP_NAME << " Opening file: " << fileName.Data()); + std::unique_ptr<TFile> inputFile(TFile::Open(fileName, "READ")); //probably rather "READ"? + ANA_CHECK(inputFile.get()); + + // Create a TEvent object (persistent store) + xAOD::TEvent pers(xAOD::TEvent::kClassAccess); + // Create a TStore object (transient store) + xAOD::TStore trans; + ANA_CHECK(pers.readFrom(inputFile.get())); + // + ANA_MSG_INFO("Number of events in the file: " << pers.getEntries()); + + // Decide how many events to run over: + Long64_t entries = pers.getEntries(); + if (argc > 2) + { + const Long64_t userInputEntries = atoll(argv[2]); + if (userInputEntries < entries) + { + entries = userInputEntries; + } + } + ANA_MSG_INFO("Running over " << entries << " events."); + + // initialize the tool + asg::AnaToolHandle<IElectronPhotonShowerShapeFudgeTool> myTool("ElectronPhotonVariableCorrectionTool/myTestTool"); + std::string configFile = "EGammaVariableCorrection/ElectronPhotonVariableCorrectionTool_ExampleConf.conf"; + ANA_CHECK(myTool.setProperty("ConfigFile", configFile)); + ANA_CHECK(myTool.initialize()); + + // loop over the events + for (Long64_t entry = 0; entry < entries; entry++) + { + //get entry + pers.getEntry(entry); + + // ==================================== + // Photons + // ==================================== + + //get photon container + const xAOD::PhotonContainer* photons; + ANA_CHECK(pers.retrieve(photons, "Photons")); + + // Make a deep copy of the photon container + // Create the new container and its auxiliary store. + auto photons_copy = std::make_unique<xAOD::PhotonContainer>(); + auto photons_copy_aux = std::make_unique<xAOD::AuxContainerBase>(); + photons_copy->setStore (photons_copy_aux.get()); //< Connect the two + + //copy photons over + for (auto photon : *photons) { + // Copy this photon to the output container: + xAOD::Photon* photon_copy = new xAOD::Photon(); + photons_copy->push_back (photon_copy); // photon acquires the photon_copy auxstore + *photon_copy = *photon; // copies auxdata from one auxstore to the other + } + // end make deep copy + + //loop over deep copy of photon container + for (unsigned int photon_itr = 0; photon_itr < photons_copy->size(); photon_itr++) + { + xAOD::Photon* photon = photons_copy->at(photon_itr); + + //apply correction + ANA_CHECK(myTool->applyCorrection(*photon)); + + } // loop over deep copy of photon container + + // ==================================== + // Electrons + // ==================================== + + //get electron container + const xAOD::ElectronContainer* electrons; + ANA_CHECK(pers.retrieve(electrons, "Electrons")); + + // Make a deep copy of the electron container + // Create the new container and its auxiliary store. + auto electrons_copy = std::make_unique<xAOD::ElectronContainer>(); + auto electrons_copy_aux = std::make_unique<xAOD::AuxContainerBase>(); + electrons_copy->setStore (electrons_copy_aux.get()); //< Connect the two + + //copy electrons over + for (auto electron : *electrons) { + // Copy this electron to the output container: + xAOD::Electron* electron_copy = new xAOD::Electron(); + electrons_copy->push_back (electron_copy); // electron acquires the electron_copy auxstore + *electron_copy = *electron; // copies auxdata from one auxstore to the other + } + // end make deep copy + + //loop over deep copy of electron container + for (unsigned int electron_itr = 0; electron_itr < electrons_copy->size(); electron_itr++) + { + xAOD::Electron* electron = electrons_copy->at(electron_itr); + + //apply correction + ANA_CHECK(myTool->applyCorrection(*electron)); + + } // loop over deep copy of electrone container + + } // loop over events + + return 0; +} diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool_DictionaryToolHandle.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool_DictionaryToolHandle.cxx new file mode 100644 index 0000000000000000000000000000000000000000..ec38bf5fe53f45443331975747571b0d47ea64fb --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testElectronPhotonVariableCorrectionTool_DictionaryToolHandle.cxx @@ -0,0 +1,38 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + @brief Test code to test ElectronPhotonVariableCorrectionTool Dictionaries + @author Nils Gillwald (DESY) nils.gillwald@desy.de + @date February 2020 +**/ + +// ROOT include(s): +#include <TError.h> + +// EDM include(s): +#include "AsgTools/AnaToolHandle.h" + +#include "EgammaAnalysisInterfaces/IElectronPhotonShowerShapeFudgeTool.h" + +#define MSGSOURCE "testElectronPhotonVariableCorrectionTool_DictionaryToolHandle" + +int main(/*int argc, char* argv[]*/) //unused variable warnings!! +{ + Info(MSGSOURCE, "Configuring the ElectronPhotonVariableCorrectionTool"); + asg::AnaToolHandle<IElectronPhotonShowerShapeFudgeTool> myTool("ElectronPhotonVariableCorrectionTool/myTool"); + myTool.setProperty("ConfigFile", "EGammaVariableCorrection/EGammaVariableCorrectionTool_ExampleConf.conf").ignore(); + if(myTool.initialize() != StatusCode::SUCCESS) + { + Error(MSGSOURCE, "Unable to initialize the ElectronPhotonVariableCorrectionTool!"); + return 1; + } + else + { + Info(MSGSOURCE, "Initialized the ElectronPhotonVariableCorrectionTool!"); + } + + + return 0; +} //end main diff --git a/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testIsoCorrection.cxx b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testIsoCorrection.cxx new file mode 100644 index 0000000000000000000000000000000000000000..edd58d251b633fa636e859a2c271a18ca3df2e1a --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EGammaVariableCorrection/util/testIsoCorrection.cxx @@ -0,0 +1,119 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +/** + @brief Test code to compare output of ElectronPhotonVariableCorrectionBase correction to the isolation correction tool + @author Nils Gillwald (DESY) nils.gillwald@desy.de + @date February 2020 +**/ + +// ROOT includes +#include "TFile.h" +#include "TString.h" + +//EDM includes +#include "xAODEgamma/PhotonContainer.h" +#include "xAODEgamma/PhotonAuxContainer.h" +#include "xAODEgamma/Photon.h" + +// must be changed to not be a relative path! +#include "EGammaVariableCorrection/ElectronPhotonVariableCorrectionBase.h" +// +#include "AsgMessaging/AsgMessaging.h" + +//infrastructure includes +#ifdef ROOTCORE +#include "xAODRootAccess/TEvent.h" +#include "xAODRootAccess/TStore.h" +#endif //ROOTCORE +#include "IsolationCorrections/IsolationCorrectionTool.h" +#include "xAODCore/ShallowCopy.h" + +//main test code + +int main (int argc, char* argv[]) +{ + using namespace asg::msgUserCode; + // The application's name: + const char* APP_NAME = argv[0]; + // Check if we received a file name: + if (argc < 2) + { + ANA_MSG_ERROR("No file name received!" ); + ANA_MSG_ERROR( " Usage: %s [xAOD file name] %d Num of events to process %d (0 photons , 1 electrons)"); + return EXIT_FAILURE; + } + + //open input file + const TString fileName = argv[1]; + ANA_MSG_INFO(APP_NAME << " Opening file: " << fileName.Data()); + std::unique_ptr<TFile> inputFile(TFile::Open(fileName, "READ")); + ANA_CHECK(inputFile.get()); + + // Create a TEvent object (persistent store) + xAOD::TEvent pers(xAOD::TEvent::kClassAccess); + // Create a TStore object (transient store) + xAOD::TStore trans; + ANA_CHECK(pers.readFrom(inputFile.get())); + + ANA_MSG_INFO("Number of events in the file: " << pers.getEntries()); + + // Decide how many events to run over: + Long64_t entries = pers.getEntries(); + if (argc > 2) + { + const Long64_t userInputEntries = atoll(argv[2]); + if (userInputEntries < entries) + { + entries = userInputEntries; + } + } + + //initialise the tool + ElectronPhotonVariableCorrectionBase fudgeTool("fudgeTool"); + ANA_CHECK(fudgeTool.setProperty("ConfigFile","ElectronPhotonShowerShapeFudgeTool/ElectronPhotonVariableCorrectionBase_ExampleIsoCorrectionConf.conf")); + ANA_CHECK(fudgeTool.initialize()); + + //crosscheck with isolationcorrectiontool + // check these StatusCodes once update to a newer AnalysisBase version (not 21.2.56, probably > 21.2.94) + CP::IsolationCorrectionTool isoCorrToolFull("isoCorrToolFull"); + ANA_CHECK(isoCorrToolFull.setProperty("Apply_etaEDPar_mc_correction", true)); + ANA_CHECK(isoCorrToolFull.setProperty("Apply_etaEDParPU_correction", true)); + ANA_CHECK(isoCorrToolFull.initialize()); + + CP::IsolationCorrectionTool isoCorrToolStep("isoCorrToolStep"); + ANA_CHECK(isoCorrToolStep.setProperty("Apply_etaEDPar_mc_correction", false)); + ANA_CHECK(isoCorrToolStep.setProperty("Apply_etaEDParPU_correction", true)); + ANA_CHECK(isoCorrToolStep.initialize()); + + //loop over the events + for (Long64_t entry = 0; entry < entries; entry++) + { + //get entry + pers.getEntry(entry); + ANA_MSG_INFO("============================"); + ANA_MSG_INFO("Event: " << entry); + + //get photon container + const xAOD::PhotonContainer* photons; + ANA_CHECK(pers.retrieve(photons, "Photons")); + std::pair< xAOD::PhotonContainer*, xAOD::ShallowAuxContainer* > photons_isocorr = xAOD::shallowCopyContainer( *photons ); + std::pair< xAOD::PhotonContainer*, xAOD::ShallowAuxContainer* > photons_fudge = xAOD::shallowCopyContainer( *photons ); + + //loop over photon container + for (unsigned int idx = 0; idx < photons->size(); idx++) + { + ANA_MSG_INFO("---------------------------"); + ANA_MSG_INFO("Photon: " << idx); + ANA_CHECK(isoCorrToolFull.applyCorrection(*(photons_isocorr.first->at(idx)))); + ANA_CHECK(isoCorrToolStep.applyCorrection(*(photons_fudge.first->at(idx)))); + ANA_CHECK(fudgeTool.applyCorrection(*(photons_fudge.first->at(idx)))); + ANA_MSG_INFO("topoetcone40 fudge before applyCorrection: " << photons_fudge.first->at(idx)->auxdata<float>("topoetcone40_original")); + ANA_MSG_INFO("topoetcone40 fudge after applyCorrection : " << photons_fudge.first->at(idx)->auxdata<float>("topoetcone40")); + ANA_MSG_INFO("topoetcone40 with IsolationCorrectionTool : " << photons_isocorr.first->at(idx)->auxdata<float>("topoetcone40")); + } // loop over photon container + } // loop over events + + return 0; +} diff --git a/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/EgammaAnalysisHelpers/AsgEGammaConfigHelper.h b/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/EgammaAnalysisHelpers/AsgEGammaConfigHelper.h new file mode 100644 index 0000000000000000000000000000000000000000..b21125e1d9696cb8fb0b864cdcd06177de4127a9 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/EgammaAnalysisHelpers/AsgEGammaConfigHelper.h @@ -0,0 +1,30 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + + +// Dear emacs, this is -*-c++-*- +#ifndef ASGELEGAMMACONFIGHELPER_H +#define ASGELEGAMMACONFIGHELPER_H + +/** + @brief Tool to simplify the configuration of some egamma tools based on TEnv input + @author Christos Anastopoulos (2017) Nils Gillwald (DESY, 2020) nils.gillwald@desy.de + @date February 2020 +*/ +#include <string> +#include <vector> +#include <map> +class TEnv; + +namespace AsgConfigHelper{ + std::string findConfigFile (std::string input, const std::map<std::string,std::string>& configmap); + unsigned int findMask (std::string input, const std::map<std::string,unsigned int>& maskmap); + std::vector<double> HelperDouble(const std::string& input, TEnv& env); + std::vector<float> HelperFloat(const std::string& input, TEnv& env); + std::vector<int> HelperInt(const std::string& input, TEnv& env); + std::vector<std::string> HelperString(const std::string& input, TEnv& env); +} + + +#endif //ASGEGAMMACONFIGHELPER_H diff --git a/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/Root/AsgEGammaConfigHelper.cxx b/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/Root/AsgEGammaConfigHelper.cxx new file mode 100644 index 0000000000000000000000000000000000000000..09211e9ad376cc05a6bf782a28715cf8d44f7910 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/EgammaAnalysisHelpers/Root/AsgEGammaConfigHelper.cxx @@ -0,0 +1,117 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#include "EgammaAnalysisHelpers/AsgEGammaConfigHelper.h" +#include "AsgMessaging/AsgMessaging.h" +#include "TEnv.h" +#include <iostream> +#include <sstream> + +namespace AsgConfigHelper{ + + std::string findConfigFile (std::string input, const std::map<std::string,std::string>& configmap){ + auto confFile_itr=configmap.find(input); + if(confFile_itr == configmap.end()){ + static const asg::AsgMessaging msg("Egamma::AsgConfigHelper"); + msg.msg(MSG::WARNING)<<"Key "<<input<<" not found in map, no config file returned"<<endmsg; + return ""; + } + return confFile_itr->second; + } + + unsigned int findMask (std::string input, const std::map<std::string,unsigned int>& maskmap){ + auto mask_itr=maskmap.find(input); + if(mask_itr==maskmap.end()){ + static const asg::AsgMessaging msg("Egamma::AsgConfigHelper"); + msg.msg(MSG::WARNING)<<"Key "<<input<<" not found in map, egammaPID::EgPidUndefined mask returned"<<endmsg; + // mask has the choice to default to all 1 or all 0 bits, choose the former + return std::numeric_limits<unsigned int>::max(); + } + return static_cast < unsigned int> (mask_itr->second); + } + + template <typename T> + bool strtof(const std::string& input, T& f){ + int diff = 0 ; + std::string tmp = input; + std::string::size_type first(0); + std::string::size_type last(0); + first = ( input.find("#") ) ; + + //if we do not find a comment character "#" we are fine + if (first == std::string::npos) { + std::istringstream buffer (tmp); + buffer>>f; + return true; + } + else { + //if we have found comment character check if it is inlined between two "#" + last = (input.find("#",first+1) ); + //if nor error + if (last == std::string::npos) { + static const asg::AsgMessaging msg("Egamma::AsgConfigHelper"); + msg.msg(MSG::WARNING)<<" Improper comment format , inline comment should be enclosed between two # "<<endmsg; + return false; + } + //else if between two "#" remove this part + diff = last - first ; + tmp= tmp.erase(first,diff+1); + std::istringstream buffer (tmp); + buffer>>f; + return true; + } + } + + template <typename T> + std::vector<T> Helper (const std::string& input, TEnv& env){ + std::vector<T> CutVector; + std::string env_input(env.GetValue(input.c_str(), "")); + if (env_input.size() > 0) { + std::string::size_type end; + do { + end = env_input.find(";"); + T myValue(0); + if(AsgConfigHelper::strtof(env_input.substr(0,end),myValue)){ + CutVector.push_back(myValue); + } + if (end != std::string::npos) { + env_input= env_input.substr(end+1); + } + } while (end != std::string::npos); + } + return CutVector; + } +} + +//use the specializations +std::vector<double> AsgConfigHelper::HelperDouble(const std::string& input, TEnv& env){ + return AsgConfigHelper::Helper<double> ( input,env); +} +std::vector<float> AsgConfigHelper::HelperFloat(const std::string& input, TEnv& env){ + return AsgConfigHelper::Helper<float> (input, env); +} +std::vector<int> AsgConfigHelper::HelperInt(const std::string& input, TEnv& env){ + return AsgConfigHelper::Helper<int> (input, env); +} +// template does not work for std::string because of the T myValue(0); declaration, so implement it again for std::string +std::vector<std::string> AsgConfigHelper::HelperString(const std::string& input, TEnv& env){ + std::vector<std::string> CutVector; + std::string env_input(env.GetValue(input.c_str(), "")); + if (env_input.size() > 0) { + std::string::size_type end; + do { + end = env_input.find(";"); + std::string myValue(""); + if(AsgConfigHelper::strtof(env_input.substr(0,end),myValue)){ + CutVector.push_back(myValue); + } + if (end != std::string::npos) { + env_input= env_input.substr(end+1); + } + } while (end != std::string::npos); + } + return CutVector; +} + + diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool.h b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool.h index 5454a2bb9337fa6c60f81e4e22206d16608a9acd..6c3bf32a613e8eb04d992d62adb96596d3b97a82 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool.h +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool.h @@ -16,11 +16,7 @@ */ // Framework include(s): #include "AsgTools/AsgTool.h" -#ifdef USE_NEW_TOOL #include "ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h" -#else -#include "ElectronPhotonShowerShapeFudgeTool/FudgeMCTool.h" -#endif #include "ElectronPhotonShowerShapeFudgeTool/TElectronMCShifterTool.h" #include "EgammaAnalysisInterfaces/IElectronPhotonShowerShapeFudgeTool.h" #include "TEnv.h" @@ -57,16 +53,13 @@ public: private: -#ifdef USE_NEW_TOOL TPhotonMCShifterTool* m_ph_rootTool; -#else - FudgeMCTool* m_ph_rootTool; -#endif TElectronMCShifterTool* m_el_rootTool; int m_preselection; std::string m_configFile; + std::string m_ffFile; /** Copied over from the configuration helper so that the selector tools do not need to be included */ std::vector<float> GetFloatVector(const std::string& input, TEnv& env); diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h index 68f49fa3d83f81dd414e46685a5531d1ebda4a31..e5980262b886c6e804ef7d4aaca67ede5c06162b 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h @@ -2,380 +2,454 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -#ifndef TPhotonMCShifterTool_h -#define TPhotonMCShifterTool_h - -/* @author Giovanni Marchiori <giovanni.marchiori@cern.ch> - * - * Reimplementation of FudgeMCTool with photon fudge factors - * stored in ROOT files instead of a C++ file - * - */ - -#include <iostream> -#include <stdio.h> -#include <cmath> -#include "TH2D.h" -#include "TGraphErrors.h" -#include "TSystem.h" - -namespace IDVAR -{ - const int RHAD1 = 0; - const int RHAD = 1; - const int E277 = 2; - const int RETA = 3; - const int RPHI = 4; - const int WETA2 = 5; - const int F1 = 6; - const int FSIDE = 7; - const int WTOT = 8; - const int W1 = 9; - const int DE = 10; - const int ERATIO = 11; -} - - -class TPhotonMCShifterTool -{ - - public: - - TPhotonMCShifterTool(); - ~TPhotonMCShifterTool(); - - //fudge shower shapes for preselection chosen. - void FudgeShowers( double pt , - double eta2 , - double& rhad1 , - double& rhad , - double& e277 , - double& reta , - double& rphi , - double& weta2 , - double& f1 , - double& fside , - double& wtot , - double& w1 , - double& deltae , - double& eratio , - int isConv , - int preselection=-999); - - //fudge shower shapes for preselection chosen. - void FudgeShowers( float pt , - float eta2 , - float& rhad1 , - float& rhad , - float& e277 , - float& reta , - float& rphi , - float& weta2 , - float& f1 , - float& fside , - float& wtot , - float& w1 , - float& deltae , - float& eratio , - int isConv , - int preselection=-999); - - // fudge showers using D3PD vectors (except for eratio) - void FudgeShowers( std::vector<float> clE, - std::vector<float> eta2 , - std::vector<float>& rhad1 , - std::vector<float>& rhad , - std::vector<float>& e277 , - std::vector<float>& reta , - std::vector<float>& rphi , - std::vector<float>& weta2 , - std::vector<float>& f1 , - std::vector<float>& fside , - std::vector<float>& wtot , - std::vector<float>& w1 , - std::vector<float>& deltae , - std::vector<float>& eratio , - std::vector<int> isConv , - int preselection=-999); - - /** Calculate Eratio (the only variable which is not stored in the egamma/PhotonD3PD) - * given emaxs1 and Emax2. Return a pointer to a vector<float> - **/ - static std::vector<float>* getEratio(std::vector<float> emaxs1, std::vector<float> Emax2) - { - std::vector<float> *eratio = new std::vector<float>(); - for (unsigned int i = 0; i < emaxs1.size(); ++i) - eratio->push_back(emaxs1[i] + Emax2[i] == 0 ? 0 : - (emaxs1[i] - Emax2[i])/(emaxs1[i] + Emax2[i])); - return eratio; - } - - // get fudge factor for predefined preselection - double GetFF_Rhad1 (double pt, double eta2, int conv){ - if (conv) return (h_c_rhad1->GetBinContent(h_u_rhad1->FindBin(pt, fabs(eta2)))); - else return (h_u_rhad1->GetBinContent(h_u_rhad1->FindBin(pt, fabs(eta2))));}; - - double GetFF_Rhad (double pt, double eta2, int conv){ - if (conv) return (h_c_rhad->GetBinContent(h_c_rhad->FindBin(pt, fabs(eta2)))); - else return (h_u_rhad->GetBinContent(h_u_rhad->FindBin(pt, fabs(eta2))));}; - - double GetFF_E277 (double pt, double eta2, int conv){ - if (conv) return (h_c_e277->GetBinContent(h_c_e277->FindBin(pt, fabs(eta2)))); - else return (h_u_e277->GetBinContent(h_u_e277->FindBin(pt, fabs(eta2))));}; - - double GetFF_Reta (double pt, double eta2, int conv){ - if (conv) return (h_c_reta->GetBinContent(h_c_reta->FindBin(pt, fabs(eta2)))); - else return (h_u_reta->GetBinContent(h_u_reta->FindBin(pt, fabs(eta2))));}; - - double GetFF_Rphi (double pt, double eta2, int conv){ - if (conv) return (h_c_rphi->GetBinContent(h_c_rphi->FindBin(pt, fabs(eta2)))); - else return (h_u_rphi->GetBinContent(h_u_rphi->FindBin(pt, fabs(eta2))));}; - - double GetFF_Weta2 (double pt, double eta2, int conv){ - if (conv) return (h_c_weta2->GetBinContent(h_c_weta2->FindBin(pt, fabs(eta2)))); - else return (h_u_weta2->GetBinContent(h_u_weta2->FindBin(pt, fabs(eta2))));}; - - double GetFF_F1 (double pt, double eta2, int conv){ - if (conv) return (h_c_f1->GetBinContent(h_c_f1->FindBin(pt, fabs(eta2)))); - else return (h_u_f1->GetBinContent(h_u_f1->FindBin(pt, fabs(eta2))));}; - - double GetFF_DE (double pt, double eta2, int conv){ - if (conv) return (h_c_de->GetBinContent(h_c_de->FindBin(pt, fabs(eta2)))); - else return (h_u_de->GetBinContent(h_u_de->FindBin(pt, fabs(eta2))));}; - - double GetFF_Eratio(double pt, double eta2, int conv){ - if (conv) return (h_c_eratio->GetBinContent(h_c_eratio->FindBin(pt, fabs(eta2)))); - else return (h_u_eratio->GetBinContent(h_u_eratio->FindBin(pt, fabs(eta2))));}; - - double GetFF_Fside (double pt, double eta2, int conv){ - if (conv) return (h_c_fside->GetBinContent(h_c_fside->FindBin(pt, fabs(eta2)))); - else return (h_u_fside->GetBinContent(h_u_fside->FindBin(pt, fabs(eta2))));}; - - double GetFF_Wtot (double pt, double eta2, int conv){ - if (conv) return (h_c_wtot->GetBinContent(h_c_wtot->FindBin(pt, fabs(eta2)))); - else return (h_u_wtot->GetBinContent(h_u_wtot->FindBin(pt, fabs(eta2))));}; - - double GetFF_W1 (double pt, double eta2, int conv){ - if (conv) return (h_c_w1->GetBinContent(h_c_w1->FindBin(pt, fabs(eta2)))); - else return (h_u_w1->GetBinContent(h_u_w1->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Rhad1 (double pt, double eta2, int conv){ - if (conv) return (h_c_rhad1->GetBinError(h_u_rhad1->FindBin(pt, fabs(eta2)))); - else return (h_u_rhad1->GetBinError(h_u_rhad1->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Rhad (double pt, double eta2, int conv){ - if (conv) return (h_c_rhad->GetBinError(h_c_rhad->FindBin(pt, fabs(eta2)))); - else return (h_u_rhad->GetBinError(h_u_rhad->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_E277 (double pt, double eta2, int conv){ - if (conv) return (h_c_e277->GetBinError(h_c_e277->FindBin(pt, fabs(eta2)))); - else return (h_u_e277->GetBinError(h_u_e277->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Reta (double pt, double eta2, int conv){ - if (conv) return (h_c_reta->GetBinError(h_c_reta->FindBin(pt, fabs(eta2)))); - else return (h_u_reta->GetBinError(h_u_reta->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Rphi (double pt, double eta2, int conv){ - if (conv) return (h_c_rphi->GetBinError(h_c_rphi->FindBin(pt, fabs(eta2)))); - else return (h_u_rphi->GetBinError(h_u_rphi->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Weta2 (double pt, double eta2, int conv){ - if (conv) return (h_c_weta2->GetBinError(h_c_weta2->FindBin(pt, fabs(eta2)))); - else return (h_u_weta2->GetBinError(h_u_weta2->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_F1 (double pt, double eta2, int conv){ - if (conv) return (h_c_f1->GetBinError(h_c_f1->FindBin(pt, fabs(eta2)))); - else return (h_u_f1->GetBinError(h_u_f1->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_DE (double pt, double eta2, int conv){ - if (conv) return (h_c_de->GetBinError(h_c_de->FindBin(pt, fabs(eta2)))); - else return (h_u_de->GetBinError(h_u_de->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Eratio(double pt, double eta2, int conv){ - if (conv) return (h_c_eratio->GetBinError(h_c_eratio->FindBin(pt, fabs(eta2)))); - else return (h_u_eratio->GetBinError(h_u_eratio->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Fside (double pt, double eta2, int conv){ - if (conv) return (h_c_fside->GetBinError(h_c_fside->FindBin(pt, fabs(eta2)))); - else return (h_u_fside->GetBinError(h_u_fside->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_Wtot (double pt, double eta2, int conv){ - if (conv) return (h_c_wtot->GetBinError(h_c_wtot->FindBin(pt, fabs(eta2)))); - else return (h_u_wtot->GetBinError(h_u_wtot->FindBin(pt, fabs(eta2))));}; - - double GetFFerr_W1 (double pt, double eta2, int conv){ - if (conv) return (h_c_w1->GetBinError(h_c_w1->FindBin(pt, fabs(eta2)))); - else return (h_u_w1->GetBinError(h_u_w1->FindBin(pt, fabs(eta2))));}; - - double GetFF (int var, double pt, double eta2, int conv){ - switch (var) { - case IDVAR::RHAD1: return GetFF_Rhad1( pt, eta2, conv ); - case IDVAR::RHAD: return GetFF_Rhad( pt, eta2, conv ); - case IDVAR::E277: return GetFF_E277( pt, eta2, conv ); - case IDVAR::RETA: return GetFF_Reta( pt, eta2, conv ); - case IDVAR::RPHI: return GetFF_Rphi( pt, eta2, conv ); - case IDVAR::WETA2: return GetFF_Weta2( pt, eta2, conv ); - case IDVAR::F1: return GetFF_F1( pt, eta2, conv ); - case IDVAR::FSIDE: return GetFF_Fside( pt, eta2, conv ); - case IDVAR::WTOT: return GetFF_Wtot( pt, eta2, conv ); - case IDVAR::W1: return GetFF_W1( pt, eta2, conv ); - case IDVAR::DE: return GetFF_DE( pt, eta2, conv ); - case IDVAR::ERATIO: return GetFF_Eratio( pt, eta2, conv ); - default: return 0.0; - } - } - - // fudge a specific variable - double Fudge_Rhad1 ( double rhad1, double pt, double eta2, int conv){ return ( rhad1 + GetFF_Rhad1 ( pt, eta2, conv ) ); } - double Fudge_Rhad ( double rhad, double pt, double eta2, int conv){ return ( rhad + GetFF_Rhad ( pt, eta2, conv ) ); } - double Fudge_E277 ( double e277, double pt, double eta2, int conv){ return ( e277 + GetFF_E277 ( pt, eta2, conv ) ); } - double Fudge_Reta ( double reta, double pt, double eta2, int conv){ return ( reta + GetFF_Reta ( pt, eta2, conv ) ); } - double Fudge_Rphi ( double rphi, double pt, double eta2, int conv){ return ( rphi + GetFF_Rphi ( pt, eta2, conv ) ); } - double Fudge_Weta2 ( double weta2, double pt, double eta2, int conv){ return ( weta2 + GetFF_Weta2 ( pt, eta2, conv ) ); } - double Fudge_F1 ( double f1, double pt, double eta2, int conv){ return ( f1 + GetFF_F1 ( pt, eta2, conv ) ); } - double Fudge_DE ( double deltae, double pt, double eta2, int conv){ return ( deltae + GetFF_DE ( pt, eta2, conv ) ); } - double Fudge_Eratio( double eratio, double pt, double eta2, int conv){ return ( eratio + GetFF_Eratio ( pt, eta2, conv ) ); } - double Fudge_Fside ( double fside, double pt, double eta2, int conv){ return ( fside + GetFF_Fside ( pt, eta2, conv ) ); } - double Fudge_Wtot ( double wtot, double pt, double eta2, int conv){ return ( wtot + GetFF_Wtot ( pt, eta2, conv ) ); } - double Fudge_W1 ( double w1, double pt, double eta2, int conv){ return ( w1 + GetFF_W1 ( pt, eta2, conv ) ); } - - // set shower preselection cuts - // *** 2010 - // 0 = tight isolated - // 1 = loose isolated - // 2 = tightx isolated - // 3 = tighty isolated - // 4 = tight isolated . Distorted material - // *** 2011 - // 5 = tight isolated - // 6 = tight non-isolated - // *** - // 10= tight isolated. Old menu rel15 (tune 3 in PhotonIDTool). - // .. - void SetVerbose( bool verbose=true ){ m_verbose=verbose; } - int GetPreselection(){ return m_preselection; }; - - //drawing methods - TGraphErrors* GetFFmap (int var, double eta, int isConv, int preselection); - - TGraphErrors* GetFFmap_Rhad1 (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Rhad (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_E277 (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Reta (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Rphi (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Weta2 (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_F1 (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Fside (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Wtot (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_W1 (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_DE (double eta, int isConv, int preselection); - TGraphErrors* GetFFmap_Eratio(double eta, int isConv, int preselection); - // TH2D* GetFFTH2D(int var, int isConv, int preselection); - void LoadFFs(int preselection); - - private: - bool m_verbose; - - //shower preselection to extract FFs - int m_preselection; - - // collections of fudge factors - TH2D* h_u_rhad1; - TH2D* h_u_rhad; - TH2D* h_u_e277; - TH2D* h_u_reta; - TH2D* h_u_rphi; - TH2D* h_u_weta2; - TH2D* h_u_f1; - TH2D* h_u_fside; - TH2D* h_u_wtot; - TH2D* h_u_w1; - TH2D* h_u_de; - TH2D* h_u_eratio; - TH2D* h_c_rhad1; - TH2D* h_c_rhad; - TH2D* h_c_e277; - TH2D* h_c_reta; - TH2D* h_c_rphi; - TH2D* h_c_weta2; - TH2D* h_c_f1; - TH2D* h_c_fside; - TH2D* h_c_wtot; - TH2D* h_c_w1; - TH2D* h_c_de; - TH2D* h_c_eratio; - - void SafeDeleteHistograms(); -}; - -#endif // #ifdef TPhotonMCShifterTool_h - -#ifdef TPhotonMCShifterTool_cxx - -TPhotonMCShifterTool::TPhotonMCShifterTool() -{ - m_preselection = -1; - h_u_rhad1 = 0; - h_u_rhad = 0; - h_u_e277 = 0; - h_u_reta = 0; - h_u_rphi = 0; - h_u_weta2 = 0; - h_u_f1 = 0; - h_u_fside = 0; - h_u_wtot = 0; - h_u_w1 = 0; - h_u_de = 0; - h_u_eratio = 0; - h_c_rhad1 = 0; - h_c_rhad = 0; - h_c_e277 = 0; - h_c_reta = 0; - h_c_rphi = 0; - h_c_weta2 = 0; - h_c_f1 = 0; - h_c_fside = 0; - h_c_wtot = 0; - h_c_w1 = 0; - h_c_de = 0; - h_c_eratio = 0; - - m_verbose = false; - // std::cout << "Initiliazing the tool ---->>> preselection " << preselection << std::endl; - -} - -TPhotonMCShifterTool::~TPhotonMCShifterTool() -{ - SafeDeleteHistograms(); -} - -void TPhotonMCShifterTool::SafeDeleteHistograms() { - if (h_u_rhad1) {delete h_u_rhad1; h_u_rhad1 = 0;} - if (h_u_rhad) {delete h_u_rhad; h_u_rhad = 0;} - if (h_u_e277) {delete h_u_e277; h_u_e277 = 0;} - if (h_u_reta) {delete h_u_reta; h_u_reta = 0;} - if (h_u_rphi) {delete h_u_rphi; h_u_rphi = 0;} - if (h_u_weta2) {delete h_u_weta2; h_u_weta2 = 0;} - if (h_u_f1) {delete h_u_f1; h_u_f1 = 0;} - if (h_u_fside) {delete h_u_fside; h_u_fside = 0;} - if (h_u_wtot) {delete h_u_wtot; h_u_wtot = 0;} - if (h_u_w1) {delete h_u_w1; h_u_w1 = 0;} - if (h_u_de) {delete h_u_de; h_u_de = 0;} - if (h_u_eratio) {delete h_u_eratio; h_u_eratio = 0;} - if (h_c_rhad1) {delete h_c_rhad1; h_c_rhad1 = 0;} - if (h_c_rhad) {delete h_c_rhad; h_c_rhad = 0;} - if (h_c_e277) {delete h_c_e277; h_c_e277 = 0;} - if (h_c_reta) {delete h_c_reta; h_c_reta = 0;} - if (h_c_rphi) {delete h_c_rphi; h_c_rphi = 0;} - if (h_c_weta2) {delete h_c_weta2; h_c_weta2 = 0;} - if (h_c_f1) {delete h_c_f1; h_c_f1 = 0;} - if (h_c_fside) {delete h_c_fside; h_c_fside = 0;} - if (h_c_wtot) {delete h_c_wtot; h_c_wtot = 0;} - if (h_c_w1) {delete h_c_w1; h_c_w1 = 0;} - if (h_c_de) {delete h_c_de; h_c_de = 0;} - if (h_c_eratio) {delete h_c_eratio; h_c_eratio = 0;} -} - -#endif// #ifdef TPhotonMCShifterTool_cxx - +#ifndef TPhotonMCShifterTool_h +#define TPhotonMCShifterTool_h + +/* @author Giovanni Marchiori <giovanni.marchiori@cern.ch> + * + * Reimplementation of FudgeMCTool with photon fudge factors + * stored in ROOT files instead of a C++ file + * + */ + +#include <iostream> +#include <stdio.h> +#include <cmath> +#include "TH2D.h" +#include "TGraphErrors.h" +#include "TSystem.h" + +namespace IDVAR +{ + const int RHAD1 = 0; + const int RHAD = 1; + const int E277 = 2; + const int RETA = 3; + const int RPHI = 4; + const int WETA2 = 5; + const int F1 = 6; + const int FSIDE = 7; + const int WTOT = 8; + const int W1 = 9; + const int DE = 10; + const int ERATIO = 11; +} + + +class TPhotonMCShifterTool +{ + + public: + + TPhotonMCShifterTool(); + ~TPhotonMCShifterTool(); + + //fudge shower shapes for preselection chosen. + void FudgeShowers( double pt , + double eta2 , + double& rhad1 , + double& rhad , + double& e277 , + double& reta , + double& rphi , + double& weta2 , + double& f1 , + double& fside , + double& wtot , + double& w1 , + double& deltae , + double& eratio , + int isConv , + int preselection=-999); + + //fudge shower shapes for preselection chosen. + void FudgeShowers( float pt , + float eta2 , + float& rhad1 , + float& rhad , + float& e277 , + float& reta , + float& rphi , + float& weta2 , + float& f1 , + float& fside , + float& wtot , + float& w1 , + float& deltae , + float& eratio , + int isConv , + int preselection=-999); + + // fudge showers using D3PD vectors (except for eratio) + void FudgeShowers( std::vector<float> clE, + std::vector<float> eta2 , + std::vector<float>& rhad1 , + std::vector<float>& rhad , + std::vector<float>& e277 , + std::vector<float>& reta , + std::vector<float>& rphi , + std::vector<float>& weta2 , + std::vector<float>& f1 , + std::vector<float>& fside , + std::vector<float>& wtot , + std::vector<float>& w1 , + std::vector<float>& deltae , + std::vector<float>& eratio , + std::vector<int> isConv , + int preselection=-999); + + /** Calculate Eratio (the only variable which is not stored in the egamma/PhotonD3PD) + * given emaxs1 and Emax2. Return a pointer to a vector<float> + **/ + static std::vector<float>* getEratio(std::vector<float> emaxs1, std::vector<float> Emax2) + { + std::vector<float> *eratio = new std::vector<float>(); + for (unsigned int i = 0; i < emaxs1.size(); ++i) + eratio->push_back(emaxs1[i] + Emax2[i] == 0 ? 0 : + (emaxs1[i] - Emax2[i])/(emaxs1[i] + Emax2[i])); + return eratio; + } + + // get fudge factor for predefined preselection + double GetFF_Rhad1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_rhad1->GetBinContent(h_u_rhad1->FindBin(pt, fabs(eta2))); } + else + { return h_u_rhad1->GetBinContent(h_u_rhad1->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Rhad (double pt, double eta2, int conv){ + if (conv) + { return h_c_rhad->GetBinContent(h_c_rhad->FindBin(pt, fabs(eta2))); } + else + { return h_u_rhad->GetBinContent(h_u_rhad->FindBin(pt, fabs(eta2))); } + } + + double GetFF_E277 (double pt, double eta2, int conv){ + if (conv) + { return h_c_e277->GetBinContent(h_c_e277->FindBin(pt, fabs(eta2))); } + else + { return h_u_e277->GetBinContent(h_u_e277->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Reta (double pt, double eta2, int conv){ + if (conv) + { return h_c_reta->GetBinContent(h_c_reta->FindBin(pt, fabs(eta2))); } + else + { return h_u_reta->GetBinContent(h_u_reta->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Rphi (double pt, double eta2, int conv){ + if (conv) + { return h_c_rphi->GetBinContent(h_c_rphi->FindBin(pt, fabs(eta2))); } + else + { return h_u_rphi->GetBinContent(h_u_rphi->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Weta2 (double pt, double eta2, int conv){ + if (conv) + { return h_c_weta2->GetBinContent(h_c_weta2->FindBin(pt, fabs(eta2))); } + else + { return h_u_weta2->GetBinContent(h_u_weta2->FindBin(pt, fabs(eta2))); } + } + + double GetFF_F1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_f1->GetBinContent(h_c_f1->FindBin(pt, fabs(eta2))); } + else + { return h_u_f1->GetBinContent(h_u_f1->FindBin(pt, fabs(eta2))); } + } + + double GetFF_DE (double pt, double eta2, int conv){ + if (conv) + { return h_c_de->GetBinContent(h_c_de->FindBin(pt, fabs(eta2))); } + else + { return h_u_de->GetBinContent(h_u_de->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Eratio(double pt, double eta2, int conv){ + if (conv) + { return h_c_eratio->GetBinContent(h_c_eratio->FindBin(pt, fabs(eta2))); } + else + { return h_u_eratio->GetBinContent(h_u_eratio->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Fside (double pt, double eta2, int conv){ + if (conv) + { return h_c_fside->GetBinContent(h_c_fside->FindBin(pt, fabs(eta2))); } + else + { return h_u_fside->GetBinContent(h_u_fside->FindBin(pt, fabs(eta2))); } + } + + double GetFF_Wtot (double pt, double eta2, int conv){ + if (conv) + { return h_c_wtot->GetBinContent(h_c_wtot->FindBin(pt, fabs(eta2))); } + else + { return h_u_wtot->GetBinContent(h_u_wtot->FindBin(pt, fabs(eta2))); } + } + + double GetFF_W1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_w1->GetBinContent(h_c_w1->FindBin(pt, fabs(eta2))); } + else + { return h_u_w1->GetBinContent(h_u_w1->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Rhad1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_rhad1->GetBinError(h_u_rhad1->FindBin(pt, fabs(eta2))); } + else + { return h_u_rhad1->GetBinError(h_u_rhad1->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Rhad (double pt, double eta2, int conv){ + if (conv) + { return h_c_rhad->GetBinError(h_c_rhad->FindBin(pt, fabs(eta2))); } + else + { return h_u_rhad->GetBinError(h_u_rhad->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_E277 (double pt, double eta2, int conv){ + if (conv) + { return h_c_e277->GetBinError(h_c_e277->FindBin(pt, fabs(eta2))); } + else + { return h_u_e277->GetBinError(h_u_e277->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Reta (double pt, double eta2, int conv){ + if (conv) + { return h_c_reta->GetBinError(h_c_reta->FindBin(pt, fabs(eta2))); } + else + { return h_u_reta->GetBinError(h_u_reta->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Rphi (double pt, double eta2, int conv){ + if (conv) + { return h_c_rphi->GetBinError(h_c_rphi->FindBin(pt, fabs(eta2))); } + else + { return h_u_rphi->GetBinError(h_u_rphi->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Weta2 (double pt, double eta2, int conv){ + if (conv) + { return h_c_weta2->GetBinError(h_c_weta2->FindBin(pt, fabs(eta2))); } + else + { return h_u_weta2->GetBinError(h_u_weta2->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_F1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_f1->GetBinError(h_c_f1->FindBin(pt, fabs(eta2))); } + else + { return h_u_f1->GetBinError(h_u_f1->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_DE (double pt, double eta2, int conv){ + if (conv) + { return h_c_de->GetBinError(h_c_de->FindBin(pt, fabs(eta2))); } + else + { return h_u_de->GetBinError(h_u_de->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Eratio(double pt, double eta2, int conv){ + if (conv) + { return h_c_eratio->GetBinError(h_c_eratio->FindBin(pt, fabs(eta2))); } + else + { return h_u_eratio->GetBinError(h_u_eratio->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Fside (double pt, double eta2, int conv){ + if (conv) + { return h_c_fside->GetBinError(h_c_fside->FindBin(pt, fabs(eta2))); } + else + { return h_u_fside->GetBinError(h_u_fside->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_Wtot (double pt, double eta2, int conv){ + if (conv) + { return h_c_wtot->GetBinError(h_c_wtot->FindBin(pt, fabs(eta2))); } + else + { return h_u_wtot->GetBinError(h_u_wtot->FindBin(pt, fabs(eta2))); } + } + + double GetFFerr_W1 (double pt, double eta2, int conv){ + if (conv) + { return h_c_w1->GetBinError(h_c_w1->FindBin(pt, fabs(eta2))); } + else + { return h_u_w1->GetBinError(h_u_w1->FindBin(pt, fabs(eta2))); } + } + + double GetFF (int var, double pt, double eta2, int conv){ + switch (var) { + case IDVAR::RHAD1: return GetFF_Rhad1( pt, eta2, conv ); + case IDVAR::RHAD: return GetFF_Rhad( pt, eta2, conv ); + case IDVAR::E277: return GetFF_E277( pt, eta2, conv ); + case IDVAR::RETA: return GetFF_Reta( pt, eta2, conv ); + case IDVAR::RPHI: return GetFF_Rphi( pt, eta2, conv ); + case IDVAR::WETA2: return GetFF_Weta2( pt, eta2, conv ); + case IDVAR::F1: return GetFF_F1( pt, eta2, conv ); + case IDVAR::FSIDE: return GetFF_Fside( pt, eta2, conv ); + case IDVAR::WTOT: return GetFF_Wtot( pt, eta2, conv ); + case IDVAR::W1: return GetFF_W1( pt, eta2, conv ); + case IDVAR::DE: return GetFF_DE( pt, eta2, conv ); + case IDVAR::ERATIO: return GetFF_Eratio( pt, eta2, conv ); + default: return 0.0; + } + } + + // fudge a specific variable + double Fudge_Rhad1 ( double rhad1, double pt, double eta2, int conv){ return ( rhad1 + GetFF_Rhad1 ( pt, eta2, conv ) ); } + double Fudge_Rhad ( double rhad, double pt, double eta2, int conv){ return ( rhad + GetFF_Rhad ( pt, eta2, conv ) ); } + double Fudge_E277 ( double e277, double pt, double eta2, int conv){ return ( e277 + GetFF_E277 ( pt, eta2, conv ) ); } + double Fudge_Reta ( double reta, double pt, double eta2, int conv){ return ( reta + GetFF_Reta ( pt, eta2, conv ) ); } + double Fudge_Rphi ( double rphi, double pt, double eta2, int conv){ return ( rphi + GetFF_Rphi ( pt, eta2, conv ) ); } + double Fudge_Weta2 ( double weta2, double pt, double eta2, int conv){ return ( weta2 + GetFF_Weta2 ( pt, eta2, conv ) ); } + double Fudge_F1 ( double f1, double pt, double eta2, int conv){ return ( f1 + GetFF_F1 ( pt, eta2, conv ) ); } + double Fudge_DE ( double deltae, double pt, double eta2, int conv){ return ( deltae + GetFF_DE ( pt, eta2, conv ) ); } + double Fudge_Eratio( double eratio, double pt, double eta2, int conv){ return ( eratio + GetFF_Eratio ( pt, eta2, conv ) ); } + double Fudge_Fside ( double fside, double pt, double eta2, int conv){ return ( fside + GetFF_Fside ( pt, eta2, conv ) ); } + double Fudge_Wtot ( double wtot, double pt, double eta2, int conv){ return ( wtot + GetFF_Wtot ( pt, eta2, conv ) ); } + double Fudge_W1 ( double w1, double pt, double eta2, int conv){ return ( w1 + GetFF_W1 ( pt, eta2, conv ) ); } + + // set shower preselection cuts + // *** 2010 + // 0 = tight isolated + // 1 = loose isolated + // 2 = tightx isolated + // 3 = tighty isolated + // 4 = tight isolated . Distorted material + // *** 2011 + // 5 = tight isolated + // 6 = tight non-isolated + // *** + // 10= tight isolated. Old menu rel15 (tune 3 in PhotonIDTool). + // .. + void SetVerbose( bool verbose=true ){ m_verbose=verbose; } + int GetPreselection(){ return m_preselection; }; + + //drawing methods + TGraphErrors* GetFFmap (int var, double eta, int isConv, int preselection); + + TGraphErrors* GetFFmap_Rhad1 (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Rhad (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_E277 (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Reta (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Rphi (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Weta2 (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_F1 (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Fside (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Wtot (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_W1 (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_DE (double eta, int isConv, int preselection); + TGraphErrors* GetFFmap_Eratio(double eta, int isConv, int preselection); + // TH2D* GetFFTH2D(int var, int isConv, int preselection); + void LoadFFs(int preselection, std::string file); + + private: + bool m_verbose; + + //shower preselection to extract FFs + int m_preselection; + + std::string m_corr_file; + + // collections of fudge factors + TH2D* h_u_rhad1; + TH2D* h_u_rhad; + TH2D* h_u_e277; + TH2D* h_u_reta; + TH2D* h_u_rphi; + TH2D* h_u_weta2; + TH2D* h_u_f1; + TH2D* h_u_fside; + TH2D* h_u_wtot; + TH2D* h_u_w1; + TH2D* h_u_de; + TH2D* h_u_eratio; + TH2D* h_c_rhad1; + TH2D* h_c_rhad; + TH2D* h_c_e277; + TH2D* h_c_reta; + TH2D* h_c_rphi; + TH2D* h_c_weta2; + TH2D* h_c_f1; + TH2D* h_c_fside; + TH2D* h_c_wtot; + TH2D* h_c_w1; + TH2D* h_c_de; + TH2D* h_c_eratio; + + void SafeDeleteHistograms(); +}; + +#endif // #ifdef TPhotonMCShifterTool_h + +#ifdef TPhotonMCShifterTool_cxx + +TPhotonMCShifterTool::TPhotonMCShifterTool() +{ + m_preselection = -1; + h_u_rhad1 = 0; + h_u_rhad = 0; + h_u_e277 = 0; + h_u_reta = 0; + h_u_rphi = 0; + h_u_weta2 = 0; + h_u_f1 = 0; + h_u_fside = 0; + h_u_wtot = 0; + h_u_w1 = 0; + h_u_de = 0; + h_u_eratio = 0; + h_c_rhad1 = 0; + h_c_rhad = 0; + h_c_e277 = 0; + h_c_reta = 0; + h_c_rphi = 0; + h_c_weta2 = 0; + h_c_f1 = 0; + h_c_fside = 0; + h_c_wtot = 0; + h_c_w1 = 0; + h_c_de = 0; + h_c_eratio = 0; + + m_verbose = false; + // std::cout << "Initiliazing the tool ---->>> preselection " << preselection << std::endl; + +} + +TPhotonMCShifterTool::~TPhotonMCShifterTool() +{ + SafeDeleteHistograms(); +} + +void TPhotonMCShifterTool::SafeDeleteHistograms() { + if (h_u_rhad1) {delete h_u_rhad1; h_u_rhad1 = 0;} + if (h_u_rhad) {delete h_u_rhad; h_u_rhad = 0;} + if (h_u_e277) {delete h_u_e277; h_u_e277 = 0;} + if (h_u_reta) {delete h_u_reta; h_u_reta = 0;} + if (h_u_rphi) {delete h_u_rphi; h_u_rphi = 0;} + if (h_u_weta2) {delete h_u_weta2; h_u_weta2 = 0;} + if (h_u_f1) {delete h_u_f1; h_u_f1 = 0;} + if (h_u_fside) {delete h_u_fside; h_u_fside = 0;} + if (h_u_wtot) {delete h_u_wtot; h_u_wtot = 0;} + if (h_u_w1) {delete h_u_w1; h_u_w1 = 0;} + if (h_u_de) {delete h_u_de; h_u_de = 0;} + if (h_u_eratio) {delete h_u_eratio; h_u_eratio = 0;} + if (h_c_rhad1) {delete h_c_rhad1; h_c_rhad1 = 0;} + if (h_c_rhad) {delete h_c_rhad; h_c_rhad = 0;} + if (h_c_e277) {delete h_c_e277; h_c_e277 = 0;} + if (h_c_reta) {delete h_c_reta; h_c_reta = 0;} + if (h_c_rphi) {delete h_c_rphi; h_c_rphi = 0;} + if (h_c_weta2) {delete h_c_weta2; h_c_weta2 = 0;} + if (h_c_f1) {delete h_c_f1; h_c_f1 = 0;} + if (h_c_fside) {delete h_c_fside; h_c_fside = 0;} + if (h_c_wtot) {delete h_c_wtot; h_c_wtot = 0;} + if (h_c_w1) {delete h_c_w1; h_c_w1 = 0;} + if (h_c_de) {delete h_c_de; h_c_de = 0;} + if (h_c_eratio) {delete h_c_eratio; h_c_eratio = 0;} +} + +#endif// #ifdef TPhotonMCShifterTool_cxx + diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/ElectronPhotonShowerShapeFudgeTool.cxx b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/ElectronPhotonShowerShapeFudgeTool.cxx index e9c964fb55acd84fb52359043c3fa4d8308b7262..64af52d65c7d26ceb593e60bc29cf171bc7b85f7 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/ElectronPhotonShowerShapeFudgeTool.cxx +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/ElectronPhotonShowerShapeFudgeTool.cxx @@ -28,15 +28,12 @@ ElectronPhotonShowerShapeFudgeTool::ElectronPhotonShowerShapeFudgeTool(std::stri m_configFile("") { + declareProperty("FFCalibFile", m_ffFile="ElectronPhotonShowerShapeFudgeTool/v2/PhotonFudgeFactors.root", "Calib path file for Photon MC corrections"); declareProperty("Preselection",m_preselection=-999); declareProperty("ConfigFile",m_configFile="","The config file to use for the Electron Shifter"); // Create an instance of the underlying ROOT tool -#ifdef USE_NEW_TOOL m_ph_rootTool = new TPhotonMCShifterTool(); -#else - m_ph_rootTool = new FudgeMCTool(); -#endif m_el_rootTool = new TElectronMCShifterTool(); } @@ -102,6 +99,8 @@ StatusCode ElectronPhotonShowerShapeFudgeTool::initialize() m_el_rootTool->Widths[ElePIDNames::Var::e277] = GetFloatVector("width_e277", env); m_el_rootTool->Widths[ElePIDNames::Var::DeltaE] = GetFloatVector("width_DeltaE", env); + m_ph_rootTool->LoadFFs(m_preselection, m_ffFile); + return StatusCode::SUCCESS; } diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/TPhotonMCShifterTool.cxx b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/TPhotonMCShifterTool.cxx index 2383927c9e6dfcb75bb35dba67ab20ef2d238adf..4d4fa5301f3374011854eb3ab36b2813f0a4d44f 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/TPhotonMCShifterTool.cxx +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonShowerShapeFudgeTool/Root/TPhotonMCShifterTool.cxx @@ -37,7 +37,7 @@ void TPhotonMCShifterTool::FudgeShowers( double pt , if (preselection<0) preselection = m_preselection; if (preselection>0 && preselection!=m_preselection) { m_preselection = preselection; - this->LoadFFs(preselection); + this->LoadFFs(preselection,m_corr_file); } //fudge showers @@ -77,7 +77,7 @@ void TPhotonMCShifterTool::FudgeShowers( float pt , if (preselection<0) preselection = m_preselection; if (preselection>0 && preselection!=m_preselection) { m_preselection = preselection; - this->LoadFFs(preselection); + this->LoadFFs(preselection,m_corr_file); } //fudge showers @@ -93,6 +93,7 @@ void TPhotonMCShifterTool::FudgeShowers( float pt , fside = Fudge_Fside(fside, pt, eta2, isConv); w1 = Fudge_W1(w1, pt, eta2, isConv); eratio = Fudge_Eratio(eratio, pt, eta2, isConv); + } void TPhotonMCShifterTool::FudgeShowers( std::vector<float> clE, @@ -133,17 +134,21 @@ void TPhotonMCShifterTool::FudgeShowers( std::vector<float> clE, } -void TPhotonMCShifterTool::LoadFFs(int preselection) +void TPhotonMCShifterTool::LoadFFs(int preselection, std::string file) { if (preselection == m_preselection) return; m_preselection = preselection; - static const std::string FILE_NAME = - "ElectronPhotonShowerShapeFudgeTool/PhotonFudgeFactors.root"; - const std::string filePath = PathResolverFindCalibFile( FILE_NAME ); - std::unique_ptr< TFile > f( TFile::Open( filePath.c_str(), "READ" ) ); + const std::string filename = PathResolverFindCalibFile( file ); + if (filename.empty()){ + throw std::runtime_error("FudgeMCTool ERROR Could NOT resolve file name "+ file); + } else{ + Info("FudgeMCTool", "Path found = %s", filename.c_str()); + } + m_corr_file = filename; + std::unique_ptr< TFile > f( TFile::Open( m_corr_file.c_str(), "READ" ) ); if( ( ! f.get() ) || f->IsZombie() ) { - throw std::runtime_error( "Couldn't open file: " + FILE_NAME ); + throw std::runtime_error( "Couldn't open file: " + m_corr_file ); } if (!f->FindKey(Form("TUNE%d",preselection))) { std::cout << "Directory TUNE" << preselection << " does not exist in fudge factor file. Aborting" << std::endl; @@ -225,123 +230,11 @@ void TPhotonMCShifterTool::LoadFFs(int preselection) f->Close(); } -/* -TH2D* TPhotonMCShifterTool::GetFFTH2D(int var, int isConv, int preselection){ - - // read the FFs - if (preselection>=0 && preselection != m_preselection) - this->LoadFFs(preselection); - - // create the histogram - const int ptbins = m_rhad1_ff.GetPtBins()+1; - double *ptarray = m_rhad1_ff.GetPtArray(); - const int etabins = m_rhad1_ff.GetEtaBins()+1; - double *etaarray = m_rhad1_ff.GetEtaArray(); - //std::cout << ptbins << " " << etabins << std::endl; - double x[100]; - double y[100]; - - for (int i=0; i<ptbins; i++){ - if(i==0) - x[0]=0.; - else - x[i] = ptarray[i-1]; - //std::cout << i << " " << x[i] << std::endl; - } - for (int i=0; i<etabins; i++){ - if(i==0) - y[0]=0.; - else - y[i] = etaarray[i-1]; - //std::cout << i << " " << y[i] << std::endl; - } - - - - TH2D* hff = new TH2D("hff","hff",ptbins-1,x,etabins-1,y); - hff->GetXaxis()->SetTitle("p_{T} [MeV]"); - hff->GetYaxis()->SetTitle("|#eta_{S2}|"); - - // fill the histogram - double pt, eta; - for (int i=1; i<ptbins; i++){ - pt = hff->GetXaxis()->GetBinCenter(i); - for (int j=1; j<etabins; j++){ - eta = hff->GetYaxis()->GetBinCenter(j); - //std::cout << "pT = " << pt << ", eta = " << eta << std::endl; - - double ff, fferr; - - switch (var) { - case IDVAR::RHAD1: - ff = GetFF_Rhad1( pt, eta ); - fferr = GetFFerr_Rhad1( pt, eta ); - break; - case IDVAR::RHAD: - ff = GetFF_Rhad( pt, eta ); - fferr = GetFFerr_Rhad( pt, eta ); - break; - case IDVAR::E277: - ff = GetFF_E277( pt, eta ); - fferr = GetFFerr_E277( pt, eta ); - break; - case IDVAR::RETA: - ff = GetFF_Reta( pt, eta ); - fferr = GetFFerr_Reta( pt, eta ); - break; - case IDVAR::RPHI: - ff = GetFF_Rphi( pt, eta ); - fferr = GetFFerr_Rphi( pt, eta ); - break; - case IDVAR::WETA2: - ff = GetFF_Weta2( pt, eta ); - fferr = GetFFerr_Weta2( pt, eta ); - break; - case IDVAR::F1: - ff = GetFF_F1( pt, eta ); - fferr = GetFFerr_F1( pt, eta ); - break; - case IDVAR::FSIDE: - ff = GetFF_Fside( pt, eta ); - fferr = GetFFerr_Fside( pt, eta ); - break; - case IDVAR::WTOT: - ff = GetFF_Wtot( pt, eta ); - fferr = GetFFerr_Wtot( pt, eta ); - break; - case IDVAR::W1: - ff = GetFF_W1( pt, eta ); - fferr = GetFFerr_W1( pt, eta ); - break; - case IDVAR::DE: - ff = GetFF_DE( pt, eta ); - fferr = GetFFerr_DE( pt, eta ); - break; - case IDVAR::ERATIO: - ff = GetFF_Eratio( pt, eta ); - fferr = GetFFerr_Eratio( pt, eta ); - break; - default: - ff = 0.; - fferr = 0.; - break; - } - //std::cout << "FF = " << ff << " +- " << fferr << std::endl; - hff->SetBinContent(i,j,ff); - hff->SetBinError(i,j,fferr); - } - } - - // return the histogram - return hff; -} -*/ - //Get graph of FFs TGraphErrors* TPhotonMCShifterTool::GetFFmap(int var, double eta, int isConv, int preselection ){ if (preselection>=0 && preselection != m_preselection) - this->LoadFFs(preselection); + this->LoadFFs(preselection,m_corr_file); if (h_u_rhad1==0) return 0;