PPSObjects/src/PPSAssociationCuts.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  *   Grzegorz Sroka
0005  ****************************************************************************/
0006
0007 #include "CondFormats/PPSObjects/interface/PPSAssociationCuts.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009
0010 #include "TF1.h"
0011
0012 #include <iostream>
0013
0014 //----------------------------------------------------------------------------------------------------
0015
0016 PPSAssociationCuts::CutsPerArm::CutsPerArm(const edm::ParameterSet &iConfig, int sector) {
0017   const auto &association_cuts = iConfig.getParameterSet("association_cuts_" + std::to_string(sector));
0018
0019   const std::vector<std::string> names{"x", "y", "xi", "th_y"};
0020   for (std::size_t i = 0; i < names.size(); ++i) {
0021     std::string mean = association_cuts.getParameter<std::string>(names[i] + "_cut_mean");
0022     s_means_.push_back(mean);
0023
0024     std::string threshold = association_cuts.getParameter<std::string>(names[i] + "_cut_threshold");
0025     s_thresholds_.push_back(threshold);
0026   }
0027
0028   ti_tr_min_ = association_cuts.getParameter<double>("ti_tr_min");
0029   ti_tr_max_ = association_cuts.getParameter<double>("ti_tr_max");
0030 }
0031
0032 //----------------------------------------------------------------------------------------------------
0033
0034 void PPSAssociationCuts::CutsPerArm::buildFunctions() {
0035   f_means_.clear();
0036   for (const auto &s : s_means_)
0037     f_means_.push_back(std::make_shared<TF1>("f", s.c_str()));
0038
0039   f_thresholds_.clear();
0040   for (const auto &s : s_thresholds_)
0041     f_thresholds_.push_back(std::make_shared<TF1>("f", s.c_str()));
0042 }
0043
0044 //----------------------------------------------------------------------------------------------------
0045
0046 bool PPSAssociationCuts::CutsPerArm::isApplied(Quantities quantity) const {
0047   return (!s_thresholds_[quantity].empty()) && (!s_means_[quantity].empty());
0048 }
0049
0050 //----------------------------------------------------------------------------------------------------
0051
0052 bool PPSAssociationCuts::CutsPerArm::isSatisfied(
0053     Quantities quantity, double x_near, double y_near, double xangle, double q_NF_diff) const {
0054   // if cut not applied, then condition considered as satisfied
0055   if (!isApplied(quantity))
0056     return true;
0057
0058   // evaluate mean and threshold
0059   const double mean = evaluateExpression(f_means_[quantity], x_near, y_near, xangle);
0060   const double threshold = evaluateExpression(f_thresholds_[quantity], x_near, y_near, xangle);
0061
0062   // make comparison
0063   return fabs(q_NF_diff - mean) < threshold;
0064 }
0065
0066 //----------------------------------------------------------------------------------------------------
0067
0068 double PPSAssociationCuts::CutsPerArm::evaluateExpression(std::shared_ptr<TF1> expression,
0069                                                           double x_near,
0070                                                           double y_near,
0071                                                           double xangle) {
0072   expression->SetParameter("x_near", x_near);
0073   expression->SetParameter("y_near", y_near);
0074   expression->SetParameter("xangle", xangle);
0075   return expression->EvalPar(nullptr);
0076 }
0077
0078 //----------------------------------------------------------------------------------------------------
0079
0080 PPSAssociationCuts::PPSAssociationCuts(const edm::ParameterSet &iConfig) {
0081   unsigned int i = 0;
0082   for (const int &sector : {45, 56})
0083     association_cuts_[i++] = CutsPerArm(iConfig, sector);
0084
0085   initialize();
0086 }
0087
0088 //----------------------------------------------------------------------------------------------------
0089
0090 edm::ParameterSetDescription PPSAssociationCuts::getDefaultParameters() {
0091   edm::ParameterSetDescription desc;
0092
0093   desc.add<std::string>("x_cut_mean", "")->setComment("mean of track-association cut in x, mm");
0094   desc.add<std::string>("x_cut_threshold", "")->setComment("threshold of track-association cut in x, mm");
0095
0096   desc.add<std::string>("y_cut_mean", "")->setComment("mean of track-association cut in y, mm");
0097   desc.add<std::string>("y_cut_threshold", "")->setComment("threshold of track-association cut in y, mm");
0098
0099   desc.add<std::string>("xi_cut_mean", "")->setComment("mean of track-association cut in xi");
0100   desc.add<std::string>("xi_cut_threshold", "")->setComment("threshold of track-association cut in xi");
0101
0102   desc.add<std::string>("th_y_cut_mean", "")->setComment("mean of track-association cut in th_y, rad");
0103   desc.add<std::string>("th_y_cut_threshold", "")->setComment("threshold of track-association cut in th_y, rad");
0104
0105   desc.add<double>("ti_tr_min", -1.)->setComment("minimum value for timing-tracking association cut");
0106   desc.add<double>("ti_tr_max", +1.)->setComment("maximum value for timing-tracking association cut");
0107
0108   return desc;
0109 }
0110
0111 //----------------------------------------------------------------------------------------------------
0112
0113 bool PPSAssociationCuts::isValid() const {
0114   // valid if association_cuts_ have two entries, with keys 0 and 1
0115   if (association_cuts_.size() != 2)
0116     return false;
0117
0118   if (association_cuts_.find(0) == association_cuts_.end())
0119     return false;
0120   if (association_cuts_.find(1) == association_cuts_.end())
0121     return false;
0122
0123   return true;
0124 }
0125
0126 //----------------------------------------------------------------------------------------------------
0127
0128 void PPSAssociationCuts::initialize() {
0129   if (!isValid())
0130     throw cms::Exception("PPS") << "Invalid structure of PPSAssociationCuts.";
0131
0132   for (auto &p : association_cuts_)
0133     p.second.buildFunctions();
0134 }
0135
0136 //----------------------------------------------------------------------------------------------------
0137
0138 std::ostream &operator<<(std::ostream &os, const PPSAssociationCuts::CutsPerArm &cutsPerArm) {
0139   os << "CutsPerArm {" << std::endl;
0140
0141   os << "\tmeans {";
0142   for (auto const &value : cutsPerArm.getMeans()) {
0143     os << "\"" << value << "\", ";
0144   }
0145   os << "}" << std::endl << std::endl;
0146
0147   os << "\tthresholds {";
0148   for (auto const &value : cutsPerArm.getThresholds()) {
0149     os << "\"" << value << "\", ";
0150   }
0151   os << "}" << std::endl << std::endl;
0152
0153   os << "\tti_tr_min " << cutsPerArm.getTiTrMin() << std::endl;
0154   os << "\tti_tr_max " << cutsPerArm.getTiTrMax() << std::endl;
0155   os << "}" << std::endl;
0156
0157   return os;
0158 }
0159
0160 //----------------------------------------------------------------------------------------------------
0161
0162 std::ostream &operator<<(std::ostream &os, const PPSAssociationCuts &ppsAssociationCuts) {
0163   os << "PPSAssociationCuts {" << std::endl;
0164   os << "45" << std::endl;
0165   os << ppsAssociationCuts.getAssociationCuts(0);
0166   os << "56" << std::endl;
0167   os << ppsAssociationCuts.getAssociationCuts(1);
0168   os << "}" << std::endl;
0169
0170   return os;
0171 }