From a7b7dd1143c966161f81dda42635de985aa7ae53 Mon Sep 17 00:00:00 2001
From: Peter Loch <lochp@cern.ch>
Date: Fri, 6 Jul 2012 11:21:24 +0200
Subject: [PATCH] Tagged Checkreq fixes (EventShapeUtils-00-00-05)
---
.../EventShapeUtils/EventShapeCalculators.h | 394 ++++++++++++++++++
.../EventShapeUtils/EventShapeHelpers.h | 150 +++++++
.../EventShapeUtils/cmt/requirements | 28 ++
.../src/EventShapeCalculators.cxx | 305 ++++++++++++++
4 files changed, 877 insertions(+)
create mode 100644 Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeCalculators.h
create mode 100644 Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeHelpers.h
create mode 100755 Reconstruction/EventShapes/EventShapeUtils/cmt/requirements
create mode 100644 Reconstruction/EventShapes/EventShapeUtils/src/EventShapeCalculators.cxx
diff --git a/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeCalculators.h b/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeCalculators.h
new file mode 100644
index 000000000000..f6001e8c4b10
--- /dev/null
+++ b/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeCalculators.h
@@ -0,0 +1,394 @@
+// -*- c++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EVENTSHAPEUTILS_EVENTSHAPECALCULATORS_H
+#define EVENTSHAPEUTILS_EVENTSHAPECALCULATORS_H
+
+#include "TF1.h"
+
+#include <vector>
+
+namespace EventShapeCalculators
+{
+ class ProfileShape
+ {
+ public:
+
+ /*! @brief Default constructor */
+ ProfileShape();
+ /*! @brief Constructor with list of parameters
+ *
+ * @param[in] pars @c const reference to non-modifiable list of parameters
+ */
+ ProfileShape(const std::vector<double>& pars);
+ /*! @brief Constructor with list parameters and list of errors
+ *
+ * @param[in] pars @c const reference to non-modifiable list of parameters
+ * @param[in] errs @c const reference to non-modifiable list of errors
+ */
+ ProfileShape(const std::vector<double>& pars,
+ const std::vector<double>& errs);
+ /*! @brief Copy constructor */
+ ProfileShape(const ProfileShape& shape);
+ /*! @brief Base class destructor */
+ virtual ~ProfileShape();
+
+ /*! @brief Interface supporting @c ROOT
+ *
+ * For compatibility with @c ROOT . Forwards to
+ * @c ProfileShape::evaluate(xArr[0]) after setting the function
+ * parameters to the values provided in @c pArr . Useful for fitting
+ * within the @c ROOT environment.
+ *
+ * @param[in] xArr pointer to array of variables
+ * @param[in] pArr pointer to array of parameters
+ */
+ double operator()(double* xArr,double* pArr);
+
+ /*! @brief Access to function values
+ *
+ * @return Function value evaluated at @c x
+ *
+ * @param[in] x value at which function is evaluated
+ */
+ double evaluate(double x) const;
+ /*! @brief Access to function values
+ *
+ * Forwards to @c ProfileShape::evaluate(x)
+ *
+ * @return Function value evaluated at @x
+ *
+ * @param[in] x value at which function is evaluated
+ */
+ double operator()(double x) const { return this->evaluate(x); }
+
+ /*! @brief Returns number of function parameters */
+ size_t np();
+
+ /*! @brief Sets parameter value and error in function
+ *
+ * @param[in] index index of parameter in function
+ * @param[in] value parameter value
+ * @param[in] error parameter error
+ */
+ bool setParameter(size_t index,double value,double error=0.);
+
+ /*! @brief Retrieves value of parameter of function
+ *
+ * @return Value of parameter at index @c idx
+ *
+ * @param[in] idx index of parameter
+ */
+ double parameter(size_t idx) const;
+ /*! @brief Retrieves error associated with parameter of function
+ *
+ * @return Error of parameter at index @c idx
+ *
+ * @param[in] idx index of parameter
+ */
+ double error(size_t idx) const;
+ /*! @brief Set values and errors of specific parameter: plateau height
+ *
+ * Sets the (relative) transverse energy density in plateau region.
+ * Typical a reference value like 1.
+ *
+ * @param[in] value value to be set
+ * @param[in] error (optional) associated error
+ */
+ bool setPlateau(double value,double error=0.);
+ /*! @brief Set values and errors of specific parameter: plateau width
+ *
+ * Sets the size of the plateau region @f$ \eta_{\mathrm{plateau}} @f$.
+ * The plateau is given by a symmetric pseudorapidity range
+ * @f$ \left|\eta\right| < \eta_{\mathrm{plateau}} @f$
+ *
+ * @param[in] value value to be set
+ * @param[in] error (optional) associated error
+ */
+ bool setRange(double value,double error=0.);
+ /*! @brief Set values and errors of specific parameter: Gaussian slope
+ *
+ * Sets the slope of the Gaussian shoulders left and right of the
+ * plateau region. The slope is measured by @f$ \sigma_{\eta} @f$.
+ *
+ * @param[in] value value to be set
+ * @param[in] error (optional) associated error
+ */
+ bool setShape(double value,double error=0.);
+ /*! @brief Set values and errors of specific parameter: forward baseline
+ *
+ * Sets the density baseline in the forward direction
+ * @f$ \rho_{\mathrm{base}} @f$.
+ *
+ * @param[in] value value to be set
+ * @param[in] error (optional) associated error
+ */
+ bool setBase(double value,double error=0.);
+
+ /*! @brief Retrieve specific parameter value: plateau height */
+ double plateau() const;
+ /*! @brief Retrieve specific parameter value: symmetric plateau range */
+ double range() const;
+ /*! @brief Retrieve specific parameter value: shape of Gaussian shoulder */
+ double shape() const;
+ /*! @brief Retrieve specific parameter value: forward baseline */
+ double base() const;
+
+ /*! @brief Retrieve specific parameter error: plateau height */
+ double plateauError() const;
+ /*! @brief Retrieve specific parameter error: symmetric plateau range */
+ double rangeError() const;
+ /*! @brief Retrieve specific parameter error: shape of Gaussian shoulder */
+ double shapeError() const;
+ /*! @brief Retrieve specific parameter error: forward baseline */
+ double baseError() const;
+
+ private:
+
+ /*! @brief Container for parameter values list */
+ std::vector<double> m_parameters;
+ /*! @brief Container for parameter errors list */
+ std::vector<double> m_errors;
+
+ /*! @brief Index of plateau height parameter in list */
+ static size_t m_plateauIdx;
+ /*! @brief Index of plateau range parameter in list */
+ static size_t m_rangeIdx;
+ /*! @brief Index of Gaussian shoulder parameter in list */
+ static size_t m_shapeIdx;
+ /*! @brief Index of forward baseline in list */
+ static size_t m_baseIdx;
+ };
+ /*! @class ProfileShape
+ *
+ * @brief Eta profile shape descriptor and calculator.
+ *
+ * This class implements the actual profile calculation. The functional
+ * form best fitting the 2011 minimum bias data is
+ *
+ * @f[
+ * \rho(\eta) =
+ * \left\{\begin{array}{ll}
+ * \rho_{\mathrm{event}} &
+ * \left|\eta\right| < \eta_{\mathrm{plateau}} \\
+ * (\rho_{\mathrm{event}}-\rho_{\mathrm{base}})
+ * \exp((\eta-\mathrm{sgn}(\eta)
+ * \eta_{\mathrm{plateau}})^{2}/
+ * \sigma_{\eta}^{2}) + \rho_{\mathrm{base}} &
+ * \left|\eta\right| > \eta_{\mathrm{plateau}}
+ * \end{array}
+ * \right.
+ * @f]
+ *
+ */
+
+
+ class ProfileFunctions
+ {
+ public:
+
+ /*! @brief Default constructor */
+ ProfileFunctions();
+ /*! @brief Constructor with a given @f$ \eta @f$ range */
+ ProfileFunctions(double etaMin,double etaMax);
+ /*! @brief Constructor with @c ProfileShape */
+ ProfileFunctions(const ProfileShape* shape);
+ /*! @brief Copy constructor */
+ ProfileFunctions(const ProfileFunctions& profile);
+
+ /*! @brief Base class destructor */
+ virtual ~ProfileFunctions();
+
+ /*! @brief Event transverse energy calculator
+ *
+ * Implements the density calculation depending parameters.
+ *
+ * @return Event transverse energy density evaluated under certain
+ * pileup conditions and at a given direction. The returned value
+ * can be relative (@c rhoRef @c = @c 1 ) or absolute
+ * (@c rhoRef @c = @f$ \rho_{\mathrm{event}} @f$ ).
+ *
+ * @param[in] npv number of reconstructed primary vertices
+ * @param[in] mu average number of interactions
+ * @param[in] eta direction at which density is evaluated
+ * @param[in] rhoRef (optional) reference value for central density
+ * @f$ \rho_{\mathrm{event}} @f$. Defaults to 1.
+ * @param[in] window (optional) window size for profile calculation,
+ * defaults to 1.6.
+ *
+ */
+ double rho(double npv,double mu,double eta,
+ double rhoRef=1.,double window=1.6);
+ /*! @brief Event transverse energy calculator
+ *
+ * @return Event transverse energy density evaluated under certain
+ * previously set
+ * pileup conditions and at a given direction. The returned value
+ * can be relative (@c rhoRef @c = @c 1 ) or absolute
+ * (@c rhoRef @c = @f$ \rho_{\mathrm{event}} @f$ ).
+ *
+ */
+ double rho(double eta,double rhoRef);
+
+ /*! @brief Retrieve @c ROOT function of profile
+ *
+ * @return Pointer to @c ROOT function object of type @c TF1 . The profile
+ * shape depends on pileup conditions and the window size
+ * used in the calculation of the profile.
+ *
+ * @param[in] npv number of reconstructed primary vertices
+ * @param[in] mu average number of interactions
+ * @param[in] window (optional) window size for profile calculation,
+ * defaults to 1.6.
+ * @param[in] name (optional) name of function object, default is
+ * @c SlidingWindow .
+ *
+ * @note The @c TF1 function object pointed to is not owned by the
+ * @c ProfileFunctions instance. The client is responsible
+ * to delete the object whenever needed.
+ */
+ TF1* profile(double npv,double mu,double window=1.6,
+ const std::string& name="SlidingWindow");
+
+ /*! @brief Rerieve pointer to @c ProfileShape used
+ *
+ * @return @c const pointer to non-modifiable @c PofileShape used in the
+ * instance of @c ProfileFunctions .
+ *
+ * @param[in] npv number of reconstructed primary vertices
+ * @param[in] mu average number of interactions
+ * @param[in] window (optional) window size for profile calculation,
+ * defaults to 1.6.
+ */
+ const ProfileShape* shape(double npv,double mu,double window=1.6);
+
+ private:
+
+ /*! @brief @c ProfileShape object used to evaluate density */
+ mutable ProfileShape* f_shape;
+
+ /*! @brief Helper function (2nd order polynominal)*/
+ double f_poly2(double p0,double p1,double x);
+ /*! @brief Helper function (3rd order polynominal)*/
+ double f_poly3(double p0,double p1,double p2,double x);
+
+ ///////////////////////////////////////
+ // Range only depends on window size //
+ ///////////////////////////////////////
+
+ /*! @defgroup range_p Function parameters for range
+ *
+ * @brief Parametrization of range @f$ \eta_{\mathrm{plateau}} @f$
+ *
+ * The range @f$ \eta_{\mathrm{plateau}} @f$ is only a function
+ * of the window size @f$ w @f$ used to calculate the profile:
+ *
+ * @f[
+ * \eta_{\mathrm{plateau}}(w) = p_{0,\mathrm{plateau}}
+ * + p_{1,\mathrm{plateau}} \cdot w
+ * + p_{2,\mathrm{plateau}} \cdot w^{2}
+ * @f]
+ *
+ * @{
+ */
+ static double p_range_window_p0;
+ /*!< @brief Parameter @f$ p_{0,\mathrm{plateau}} @f$ */
+ static double p_range_window_p1;
+ /*!< @brief Parameter @f$ p_{1,\mathrm{plateau}} @f$ */
+ static double p_range_window_p2;
+ /*!< @brief Parameter @f$ p_{2,\mathrm{plateau}} @f$ */
+ /*! @}*/
+
+ ///////////////////////////////////////
+ // Shape only depends on window size //
+ ///////////////////////////////////////
+
+ /*! @defgroup shape_p Function parameters for shape
+ *
+ * @brief Parametrization of shape @f$ \sigma_{\eta} @f$
+ *
+ * The shape @f$ \sigma_{\eta} @f$ is only a function
+ * of the window size @f$ w @f$ used to calculate the profile:
+ *
+ * @f[
+ * \sigma_{\eta}(w) = p_{0,\mathrm{shape}}
+ * + p_{1,\mathrm{shape}} \cdot w
+ * + p_{2,\mathrm{shape}} \cdot w^{2}
+ * @f]
+ *
+ * @{
+ */
+ static double p_shape_window_p0;
+ /*!< @brief Parameter @f$ p_{0,\mathrm{shape}} @f$ */
+ static double p_shape_window_p1;
+ /*!< @brief Parameter @f$ p_{1,\mathrm{shape}} @f$ */
+ static double p_shape_window_p2;
+ /*!< @brief Parameter @f$ p_{2,\mathrm{shape}} @f$ */
+ /*! @}*/
+
+ ////////////////////////////
+ // Baseline depends on mu //
+ ////////////////////////////
+
+ /*! @defgroup base_p Function parameters for forward baseline
+ *
+ * @brief Parametrization of forward baseline @f$ \rho_{\mathrm{base}} @f$
+ *
+ * The shape @f$ \rho_{\mathrm{base}} @f$ is a function
+ * of the average number of interactions @f$ \mu @f$ and the number
+ * of reconstructed vertices @f$ N_{\mathrm{PV}} @f$:
+ *
+ * @f[
+ * \rho_{\mathrm{base}}(N_{\mathrm{PV}},\mu) =
+ * \underbrace{( p_{0,\mathrm{base}}
+ * + p_{1,\mathrm{base}} \cdot \mu )}_{\mathrm{offset}}
+ * +
+ * \underbrace{( q_{0,\mathrm{base}}
+ * + q_{1,\mathrm{base}} \cdot \mu )}_{\mathrm{slope}}
+ * \cdot N_{\mathrm{PV}}
+ * @f]
+ *
+ * @{
+ */
+ static double p_baseline_mu_offset_p0;
+ /*!< @brief Parameter @f$ p_{0,\mathrm{base}} @f$ */
+ static double p_baseline_mu_offset_p1;
+ /*!< @brief Parameter @f$ p_{1,\mathrm{base}} @f$ */
+ static double p_baseline_mu_slope_p0;
+ /*!< @brief Parameter @f$ q_{0,\mathrm{base}} @f$ */
+ static double p_baseline_mu_slope_p1;
+ /*!< @brief Parameter @f$ q_{1,\mathrm{base}} @f$ */
+ /*! @}*/
+
+ double m_window; /*!< @brief Actual window size */
+ double m_npv; /*!< @brief Actual number of reconstructed vertices */
+ double m_mu; /*!< @brief Actual average number of collisions */
+ double m_etaMin; /*!< @brief Actual @f$ \eta_{\mathrm{min}} @f$ */
+ double m_etaMax; /*!< @brief Actual @f$ \eta_{\mathrm{max}} @f$ */
+
+ };
+ /*! @class ProfileFunctions
+ *
+ * @brief Implements functions describing the @f$ \eta @f$ density profile.
+ */
+}
+
+inline double
+EventShapeCalculators::ProfileFunctions::f_poly2(double p0,double p1,double x)
+{ return p0+p1*x; }
+
+inline double
+EventShapeCalculators::ProfileFunctions::f_poly3(double p0,double p1,double p2,
+ double x)
+{ return this->f_poly2(p0,p1,x)+p2*x*x; }
+
+/*! @namespace EventShapeCalculators
+ *
+ * @brief Collection of function objects useful for shape calculations.
+ *
+ */
+
+#endif
diff --git a/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeHelpers.h b/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeHelpers.h
new file mode 100644
index 000000000000..704231ef10fe
--- /dev/null
+++ b/Reconstruction/EventShapes/EventShapeUtils/EventShapeUtils/EventShapeHelpers.h
@@ -0,0 +1,150 @@
+// -*- c++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EVENTSHAPEUTILS_EVENTSHAPEHELPERS_H
+#define EVENTSHAPEUTILS_EVENTSHAPEHELPERS_H
+
+#include "GaudiKernel/Bootstrap.h"
+
+#include "GaudiKernel/ISvcLocator.h"
+
+#include "StoreGate/StoreGateSvc.h"
+
+#include "NavFourMom/INavigable4MomentumCollection.h"
+
+#include <string>
+
+struct EventShapeHelpers
+{
+ /*! @brief Format string: reduced to fixed length
+ *
+ * This function shortens a string to a specified length @c len . If
+ * the string is of shorter or equal length, a copy is returned. If
+ * the string length exceeds @c len , the string is shortened to
+ * @c len and the last @c nf characters as replaced by a fill character
+ * specified in @c fill .
+ *
+ * @return Modified string formatted to match given length.
+ *
+ * @param[in] str @c const reference to input string (not modified)
+ * @param[in] len length of returned string
+ * @param[in] fill (optional) fill character, default is "."
+ * @param[in] nf (optional) number of characters to be replaced by
+ * @c fill at end of string. Default is 3.
+ */
+ static std::string stringFormatter(const std::string& str,size_t len,
+ const std::string& fill=".",
+ size_t nf=3)
+ {
+ if ( str.length() <= len) return str;
+ std::string outStr = str.substr(len);
+ for ( size_t i(len-nf); i<len; ++i ) { outStr.insert(i,fill); }
+ return outStr;
+ }
+ /*! @brief Retrieve object from @c StoreGate
+ *
+ * This function requires an external pointer to the @c StoreGateSvc .
+ *
+ * @param[in] pStoreGate pointer to @c StoreGate service
+ * @param[in] sgKey key of object in @c StoreGate
+ * @param[out] pObj reference to variable holding @c const pointer
+ * to retrieved object
+ *
+ * @return @c StatusCode::SUCCESS if retrieval successful, else
+ * @c StatusCode::FAILURE (can happen if object with specified key is not
+ * in transient store, or pointer to transient store is invalid)
+ */
+ template<class C>
+ static StatusCode retrieveFromSG(StoreGateSvc* pStoreGate,
+ const std::string& sgKey,
+ const C*& pObj)
+ {
+ pObj = 0;
+ if ( pStoreGate == 0 )
+ { return StatusCode::FAILURE; }
+ else
+ { return pStoreGate->retrieve(pObj,sgKey); }
+ }
+
+ /*! @brief Retrieve object from @c StoreGate
+ *
+ * This function requires an external pointer to the @c StoreGateSvc .
+ *
+ * @param[in] pStoreGate pointer to @c StoreGate service
+ * @param[in] sgKey key of object in @c StoreGate
+ * @param[out] pObj reference to variable holding pointer
+ * to retrieved object
+ *
+ * @return @c StatusCode::SUCCESS if retrieval successful, else
+ * @c StatusCode::FAILURE (can happen if object with specified key is not
+ * in transient store, or pointer to transient store is invalid)
+ */
+ template<class C>
+ static StatusCode retrieveFromSG(StoreGateSvc* pStoreGate,
+ const std::string& sgKey,
+ C*& pObj)
+ {
+ pObj = 0;
+ if ( pStoreGate == 0 )
+ { return StatusCode::FAILURE; }
+ else
+ { return pStoreGate->retrieve(pObj,sgKey); }
+ }
+
+ /*! @brief Retrieve object from @c StoreGate
+ *
+ * This function retrieves the pointer to @c StoreGate itself.
+ *
+ * @param[in] pStoreGate pointer to @c StoreGate service
+ * @param[in] sgKey key of object in @c StoreGate
+ * @param[out] pObj reference to variable holding @c const pointer
+ * to retrieved object
+ *
+ * @return @c StatusCode::SUCCESS if retrieval successful, else
+ * @c StatusCode::FAILURE (can happen if object with specified key is not
+ * in transient store, or pointer to transient store is invalid)
+ */
+ template<class C>
+ static StatusCode retrieveFromSG(const std::string& sgKey,const C*& pObj)
+ {
+ StoreGateSvc* pStoreGate = 0;
+ if ( Gaudi::svcLocator()->service("StoreGate",pStoreGate).isFailure() )
+ { return StatusCode::FAILURE; }
+ else
+ { return retrieveFromSG(pStoreGate,sgKey,pObj); }
+ }
+
+
+ /*! @brief Retrieve object from @c StoreGate
+ *
+ * This function retrieves the pointer to @c StoreGate itself.
+ *
+ * @param[in] pStoreGate pointer to @c StoreGate service
+ * @param[in] sgKey key of object in @c StoreGate
+ * @param[out] pObj reference to variable holding pointer
+ * to retrieved object
+ *
+ * @return @c StatusCode::SUCCESS if retrieval successful, else
+ * @c StatusCode::FAILURE (can happen if object with specified key is not
+ * in transient store, or pointer to transient store is invalid)
+ */
+ template<class C>
+ static StatusCode retrieveFromSG(const std::string& sgKey,C*& pObj)
+ {
+ StoreGateSvc* pStoreGate = 0;
+ if ( Gaudi::svcLocator()->service("StoreGate",pStoreGate).isFailure() )
+ { return StatusCode::FAILURE; }
+ else
+ { return retrieveFromSG(pStoreGate,sgKey,pObj); }
+ }
+};
+
+/*! @struct EventShapeHelpers
+ *
+ * @brief Collection of static functions supporting @c EventShape
+ * reconstruction
+ */
+#endif
diff --git a/Reconstruction/EventShapes/EventShapeUtils/cmt/requirements b/Reconstruction/EventShapes/EventShapeUtils/cmt/requirements
new file mode 100755
index 000000000000..0aa83e1a3cb6
--- /dev/null
+++ b/Reconstruction/EventShapes/EventShapeUtils/cmt/requirements
@@ -0,0 +1,28 @@
+package EventShapeUtils
+
+author Peter Loch <loch@physics.arizona.edu>
+
+use AtlasPolicy AtlasPolicy-*
+
+use AtlasROOT AtlasROOT-* External
+
+#use CLIDSvc CLIDSvc-00-* Control
+
+use StoreGate StoreGate-* Control
+
+use NavFourMom NavFourMom-* Event
+
+#use EventKernel EventKernel-* Event
+
+use GaudiInterface GaudiInterface-* External
+
+library EventShapeUtils *.cxx
+
+apply_pattern installed_library
+#
+#private
+#use AtlasReflex AtlasReflex-00-* External -no_auto_imports
+#
+#apply_pattern lcgdict dict=EventShapes selectionfile=selection.xml \
+# headerfiles="../EventShapeEvent/EventShapeEventDict.h"
+
diff --git a/Reconstruction/EventShapes/EventShapeUtils/src/EventShapeCalculators.cxx b/Reconstruction/EventShapes/EventShapeUtils/src/EventShapeCalculators.cxx
new file mode 100644
index 000000000000..55705cf1ed22
--- /dev/null
+++ b/Reconstruction/EventShapes/EventShapeUtils/src/EventShapeCalculators.cxx
@@ -0,0 +1,305 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "EventShapeUtils/EventShapeCalculators.h"
+
+#include "TMath.h"
+
+#include <cstdlib>
+#include <iostream>
+
+size_t EventShapeCalculators::ProfileShape::m_plateauIdx = 0;
+size_t EventShapeCalculators::ProfileShape::m_rangeIdx = 1;
+size_t EventShapeCalculators::ProfileShape::m_baseIdx = 2;
+size_t EventShapeCalculators::ProfileShape::m_shapeIdx = 3;
+
+EventShapeCalculators::ProfileShape::ProfileShape()
+ : m_parameters(4,0.)
+ , m_errors(4,0.)
+{ }
+
+EventShapeCalculators::ProfileShape::ProfileShape(const std::vector<double>&
+ pars)
+ : m_parameters(pars)
+ , m_errors(4,0.)
+{ }
+
+EventShapeCalculators::ProfileShape::ProfileShape(const std::vector<double>&
+ pars,
+ const std::vector<double>&
+ errs)
+ : m_parameters(pars)
+ , m_errors(errs)
+{ }
+
+EventShapeCalculators::ProfileShape::ProfileShape(const ProfileShape& shape)
+ : m_parameters(shape.m_parameters)
+ , m_errors(shape.m_errors)
+{ }
+
+EventShapeCalculators::ProfileShape::~ProfileShape()
+{ }
+
+// interface to ROOT
+double EventShapeCalculators::ProfileShape::operator()(double* xArr,
+ double* pArr)
+{
+ // save parameters
+ for ( size_t i(0); i<m_parameters.size(); ++i )
+ { m_parameters[i] = pArr[i]; }
+ return this->evaluate(xArr[0]);
+}
+
+double EventShapeCalculators::ProfileShape::evaluate(double x) const
+{
+ // central part is flat
+ if ( fabs(x) < m_parameters[m_rangeIdx] )
+ { return m_parameters[m_plateauIdx]; }
+
+ // near Gaussian part`
+ return
+ x > 0.
+ ? (m_parameters[m_plateauIdx]-m_parameters[m_baseIdx])
+ * TMath::Gaus(x,m_parameters[m_rangeIdx],m_parameters[m_shapeIdx]) +
+ m_parameters[m_baseIdx]
+ : (m_parameters[m_plateauIdx]-m_parameters[m_baseIdx])
+ * TMath::Gaus(x,-m_parameters[m_rangeIdx],m_parameters[m_shapeIdx]) +
+ m_parameters[m_baseIdx];
+}
+
+bool EventShapeCalculators::ProfileShape::setParameter(size_t idx,
+ double value,
+ double error)
+{
+ if ( idx >= m_parameters.size() ) return false;
+ m_parameters[idx] = value;
+ m_errors[idx] = error;
+ return true;
+}
+
+double EventShapeCalculators::ProfileShape::parameter(size_t idx) const
+{
+ if ( idx >= m_parameters.size() ) return 0.;
+ return m_parameters.at(idx);
+}
+
+double EventShapeCalculators::ProfileShape::error(size_t idx) const
+{
+ if ( idx >= m_errors.size() ) return 0.;
+ return m_errors.at(idx);
+}
+
+size_t EventShapeCalculators::ProfileShape::np()
+{ return m_parameters.size(); }
+
+bool EventShapeCalculators::ProfileShape::setPlateau(double value,double error)
+{ return this->setParameter(m_plateauIdx,value,error); }
+
+bool EventShapeCalculators::ProfileShape::setRange(double value,double error)
+{ return this->setParameter(m_rangeIdx,value,error); }
+
+bool EventShapeCalculators::ProfileShape::setShape(double value,double error)
+{ return this->setParameter(m_shapeIdx,value,error); }
+
+bool EventShapeCalculators::ProfileShape::setBase(double value,double error)
+{ return this->setParameter(m_baseIdx,value,error); }
+
+double EventShapeCalculators::ProfileShape::plateau() const
+{ return this->parameter(m_plateauIdx); }
+
+double EventShapeCalculators::ProfileShape::range() const
+{ return this->parameter(m_rangeIdx); }
+
+double EventShapeCalculators::ProfileShape::shape() const
+{ return this->parameter(m_shapeIdx); }
+
+double EventShapeCalculators::ProfileShape::base() const
+{ return this->parameter(m_baseIdx); }
+
+double EventShapeCalculators::ProfileShape::plateauError() const
+{ return this->error(m_plateauIdx); }
+
+double EventShapeCalculators::ProfileShape::rangeError() const
+{ return this->error(m_rangeIdx); }
+
+double EventShapeCalculators::ProfileShape::shapeError() const
+{ return this->error(m_shapeIdx); }
+
+double EventShapeCalculators::ProfileShape::baseError() const
+{ return this->error(m_baseIdx); }
+
+//////
+// The following function parameters are from 2011 MB data. They need to be
+// moved to the conditions DB at some point in time.
+//////
+
+#ifndef _P_RANGE_WINDOW_P0
+#define _P_RANGE_WINDOW_P0 2.21153e+00
+#endif
+
+#ifndef _P_RANGE_WINDOW_P1
+#define _P_RANGE_WINDOW_P1 -6.43196e-02
+#endif
+
+#ifndef _P_RANGE_WINDOW_P2
+#define _P_RANGE_WINDOW_P2 -1.07936e-01
+#endif
+
+#ifndef _P_SHAPE_WINDOW_P0
+#define _P_SHAPE_WINDOW_P0 3.43037e-01
+#endif
+
+#ifndef _P_SHAPE_WINDOW_P1
+#define _P_SHAPE_WINDOW_P1 3.81656e-02
+#endif
+
+#ifndef _P_SHAPE_WINDOW_P2
+#define _P_SHAPE_WINDOW_P2 1.43255e-01
+#endif
+
+#ifndef _P_BASELINE_MU_OFFSET_P0
+#define _P_BASELINE_MU_OFFSET_P0 2.5694e-02
+#endif
+
+#ifndef _P_BASELINE_MU_OFFSET_P1
+#define _P_BASELINE_MU_OFFSET_P1 6.92623e-03
+#endif
+
+#ifndef _P_BASELINE_MU_SLOPE_P0
+#define _P_BASELINE_MU_SLOPE_P0 3.63131e-02
+#endif
+
+#ifndef _P_BASELINE_MU_SLOPE_P1
+#define _P_BASELINE_MU_SLOPE_P1 1.80655e-03
+#endif
+
+double EventShapeCalculators::ProfileFunctions::p_range_window_p0
+= _P_RANGE_WINDOW_P0 ;
+double EventShapeCalculators::ProfileFunctions::p_range_window_p1
+= _P_RANGE_WINDOW_P1 ;
+double EventShapeCalculators::ProfileFunctions::p_range_window_p2
+= _P_RANGE_WINDOW_P2 ;
+
+double EventShapeCalculators::ProfileFunctions::p_shape_window_p0
+= _P_SHAPE_WINDOW_P0 ;
+double EventShapeCalculators::ProfileFunctions::p_shape_window_p1
+= _P_SHAPE_WINDOW_P1 ;
+double EventShapeCalculators::ProfileFunctions::p_shape_window_p2
+= _P_SHAPE_WINDOW_P2 ;
+
+double EventShapeCalculators::ProfileFunctions::p_baseline_mu_offset_p0
+= _P_BASELINE_MU_OFFSET_P0 ;
+double EventShapeCalculators::ProfileFunctions::p_baseline_mu_offset_p1
+= _P_BASELINE_MU_OFFSET_P1 ;
+
+double EventShapeCalculators::ProfileFunctions::p_baseline_mu_slope_p0
+= _P_BASELINE_MU_SLOPE_P0;
+double EventShapeCalculators::ProfileFunctions::p_baseline_mu_slope_p1
+= _P_BASELINE_MU_SLOPE_P1;
+
+EventShapeCalculators::ProfileFunctions::ProfileFunctions()
+ : f_shape((ProfileShape*)0)
+ , m_window(1.6)
+ , m_npv(6.)
+ , m_mu(8.)
+ , m_etaMin(-5.)
+ , m_etaMax(5.)
+{
+ std::vector<double> preLoad(4,0.);
+ f_shape = new EventShapeCalculators::ProfileShape(preLoad);
+ f_shape->setPlateau(1.);
+ this->shape(m_npv,m_mu,m_window);
+}
+
+EventShapeCalculators::ProfileFunctions::ProfileFunctions(double etaMin,
+ double etaMax)
+ : f_shape((ProfileShape*)0)
+ , m_window(1.6)
+ , m_npv(6.)
+ , m_mu(8.)
+ , m_etaMin(etaMin)
+ , m_etaMax(etaMax)
+{
+ std::vector<double> preLoad(4,0.);
+ f_shape = new EventShapeCalculators::ProfileShape(preLoad);
+ f_shape->setPlateau(1.);
+ this->shape(m_npv,m_mu,m_window);
+}
+
+EventShapeCalculators::ProfileFunctions::ProfileFunctions(const ProfileShape*
+ shape)
+ : f_shape(new ProfileShape(*shape))
+ , m_window(1.6)
+ , m_npv(6.)
+ , m_mu(8.)
+ , m_etaMin(-5.)
+ , m_etaMax(5.)
+{
+ this->shape(m_npv,m_mu,m_window);
+}
+
+EventShapeCalculators::ProfileFunctions::ProfileFunctions(const
+ ProfileFunctions&
+ profile)
+ : f_shape(new ProfileShape(*(profile.f_shape)))
+ , m_window(profile.m_window)
+ , m_npv(profile.m_npv)
+ , m_mu(profile.m_mu)
+ , m_etaMin(profile.m_etaMin)
+ , m_etaMax(profile.m_etaMax)
+{ }
+
+EventShapeCalculators::ProfileFunctions::~ProfileFunctions()
+{ if ( f_shape != 0 ) delete f_shape; }
+
+
+double EventShapeCalculators::ProfileFunctions::rho(double npv,double mu,
+ double eta,
+ double rhoRef,double window)
+{
+ if ( f_shape == 0 ) return 0.;
+ return rhoRef*(this->shape(npv,mu,window)->operator())(eta);
+}
+
+double EventShapeCalculators::ProfileFunctions::rho(double eta,double rhoRef)
+{
+ if ( f_shape == 0 ) return 0.;
+ return rhoRef*(f_shape->operator())(eta);
+}
+
+const EventShapeCalculators::ProfileShape*
+EventShapeCalculators::ProfileFunctions::shape(double npv,double mu,
+ double window)
+{
+ // calculate parameterized function parameters
+ double range(this->f_poly3(p_range_window_p0,p_range_window_p1,
+ p_range_window_p2,window));
+ double shape(this->f_poly3(p_shape_window_p0,p_shape_window_p1,
+ p_shape_window_p2,window));
+ double baseOS(this->f_poly2(p_baseline_mu_offset_p0,
+ p_baseline_mu_offset_p1,mu));
+ double baseSP(this->f_poly2(p_baseline_mu_slope_p0,
+ p_baseline_mu_slope_p1,mu));
+ double base(this->f_poly2(baseOS,baseSP,npv));
+
+ // set parameters
+ if ( !f_shape->setRange(range) )
+ printf("[ProfileShape::shape(...)] - WARNING problem setting new range value %6.3f\n",range);
+ if ( !f_shape->setShape(shape) )
+ printf("[ProfileShape::shape(...)] - WARNING problem setting new shape value %6.3f\n",shape);
+ if ( !f_shape->setBase(base) )
+ printf("[ProfileShape::shape(...)] - WARNING problem setting new base value %6.3f\n",base);
+ return f_shape;
+}
+
+TF1* EventShapeCalculators::ProfileFunctions::profile(double npv,double mu,
+ double window,
+ const std::string& name)
+{
+ this->shape(npv,mu,window);
+ TF1* f = new TF1(name.c_str(),f_shape,m_etaMin,m_etaMax,f_shape->np());
+ if ( f != 0 ) f->SetNpx(2000);
+ return f;
+}
--
GitLab