ProtonReconstruction/plugins/PPSFilteredProtonProducer.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kaspar
0004 ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0012 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0013 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0014 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0015
0016 #include <ostream>
0017
0018 //----------------------------------------------------------------------------------------------------
0019
0020 /// Module to apply Proton POG quality criteria.
0021 class PPSFilteredProtonProducer : public edm::stream::EDProducer<> {
0022 public:
0023   explicit PPSFilteredProtonProducer(const edm::ParameterSet &);
0024   ~PPSFilteredProtonProducer() override = default;
0025
0026   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0027
0028 private:
0029   void produce(edm::Event &, const edm::EventSetup &) override;
0030   void endStream() override;
0031
0032   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0033
0034   bool verbosity_;
0035
0036   double tracks_all_local_angle_x_max_, tracks_all_local_angle_y_max_;
0037
0038   std::vector<unsigned int> tracks_pixel_forbidden_RecoInfo_values_;
0039   unsigned int tracks_pixel_number_of_hits_min_;
0040   double tracks_pixel_normalised_chi_sq_max_;
0041
0042   bool protons_single_rp_include_;
0043   edm::EDGetTokenT<reco::ForwardProtonCollection> protons_single_rp_input_token_;
0044   std::string protons_single_rp_output_label_;
0045
0046   bool protons_multi_rp_include_;
0047   edm::EDGetTokenT<reco::ForwardProtonCollection> protons_multi_rp_input_token_;
0048   std::string protons_multi_rp_output_label_;
0049
0050   bool protons_multi_rp_check_valid_fit_;
0051   double protons_multi_rp_chi_sq_max_;
0052   double protons_multi_rp_normalised_chi_sq_max_;
0053
0054   /// counters
0055   unsigned int n_protons_single_rp_all, n_protons_single_rp_kept;
0056   unsigned int n_protons_multi_rp_all, n_protons_multi_rp_kept;
0057
0058   /// check one track
0059   bool IsTrackOK(const CTPPSLocalTrackLite &tr, unsigned int idx, std::ostringstream &log);
0060 };
0061
0062 //----------------------------------------------------------------------------------------------------
0063
0064 PPSFilteredProtonProducer::PPSFilteredProtonProducer(const edm::ParameterSet &iConfig)
0065     : verbosity_(iConfig.getUntrackedParameter<bool>("verbosity")),
0066       n_protons_single_rp_all(0),
0067       n_protons_single_rp_kept(0),
0068       n_protons_multi_rp_all(0),
0069       n_protons_multi_rp_kept(0) {
0070   const auto &tracks_all = iConfig.getParameterSet("tracks_all");
0071   tracks_all_local_angle_x_max_ = tracks_all.getParameter<double>("local_angle_x_max");
0072   tracks_all_local_angle_y_max_ = tracks_all.getParameter<double>("local_angle_y_max");
0073
0074   const auto &tracks_pixel = iConfig.getParameterSet("tracks_pixel");
0075   tracks_pixel_forbidden_RecoInfo_values_ =
0076       tracks_pixel.getParameter<std::vector<unsigned int>>("forbidden_RecoInfo_values");
0077   tracks_pixel_number_of_hits_min_ = tracks_pixel.getParameter<unsigned int>("number_of_hits_min");
0078   tracks_pixel_normalised_chi_sq_max_ = tracks_pixel.getParameter<double>("normalised_chi_sq_max");
0079
0080   const auto &protons_single_rp = iConfig.getParameterSet("protons_single_rp");
0081   protons_single_rp_include_ = protons_single_rp.getParameter<bool>("include");
0082   protons_single_rp_input_token_ =
0083       consumes<reco::ForwardProtonCollection>(protons_single_rp.getParameter<edm::InputTag>("input_tag"));
0084   protons_single_rp_output_label_ = protons_single_rp.getParameter<std::string>("output_label");
0085
0086   const auto &protons_multi_rp = iConfig.getParameterSet("protons_multi_rp");
0087   protons_multi_rp_include_ = protons_multi_rp.getParameter<bool>("include");
0088   protons_multi_rp_input_token_ =
0089       consumes<reco::ForwardProtonCollection>(protons_multi_rp.getParameter<edm::InputTag>("input_tag"));
0090   protons_multi_rp_output_label_ = protons_multi_rp.getParameter<std::string>("output_label");
0091   protons_multi_rp_check_valid_fit_ = protons_multi_rp.getParameter<bool>("check_valid_fit");
0092   protons_multi_rp_chi_sq_max_ = protons_multi_rp.getParameter<double>("chi_sq_max");
0093   protons_multi_rp_normalised_chi_sq_max_ = protons_multi_rp.getParameter<double>("normalised_chi_sq_max");
0094
0095   if (protons_single_rp_include_)
0096     produces<reco::ForwardProtonCollection>(protons_single_rp_output_label_);
0097
0098   if (protons_multi_rp_include_)
0099     produces<reco::ForwardProtonCollection>(protons_multi_rp_output_label_);
0100 }
0101
0102 //----------------------------------------------------------------------------------------------------
0103
0104 void PPSFilteredProtonProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0105   edm::ParameterSetDescription desc;
0106
0107   desc.addUntracked<bool>("verbosity", false)->setComment("verbosity");
0108
0109   edm::ParameterSetDescription tracks_all;
0110   tracks_all.add<double>("local_angle_x_max", 0.020)
0111       ->setComment("maximum absolute value of local horizontal angle, in rad");
0112   tracks_all.add<double>("local_angle_y_max", 0.020)
0113       ->setComment("maximum absolute value of local horizontal angle, in rad");
0114   desc.add<edm::ParameterSetDescription>("tracks_all", tracks_all)->setComment("settings for all tracks");
0115
0116   edm::ParameterSetDescription tracks_pixel;
0117   const std::vector<unsigned int> def_for_RecoInfo_vals = {
0118       (unsigned int)CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes,
0119       (unsigned int)CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes};
0120   tracks_pixel.add<std::vector<unsigned int>>("forbidden_RecoInfo_values", def_for_RecoInfo_vals)
0121       ->setComment("list of forbidden RecoInfo values");
0122   tracks_pixel.add<unsigned int>("number_of_hits_min", 0)->setComment("minimum required number of hits");
0123   tracks_pixel.add<double>("normalised_chi_sq_max", 1E100)->setComment("maximum tolerated chi square / ndof");
0124   desc.add<edm::ParameterSetDescription>("tracks_pixel", tracks_pixel)
0125       ->setComment("specific settings for pixel-RP tracks");
0126
0127   edm::ParameterSetDescription protons_single_rp;
0128   protons_single_rp.add<bool>("include", true)->setComment("flag whether single-RP protons should be processed");
0129   protons_single_rp.add<edm::InputTag>("input_tag", edm::InputTag("ctppsProtons", "singleRP"))->setComment("input tag");
0130   protons_single_rp.add<std::string>("output_label", "singleRP")->setComment("output label");
0131   desc.add<edm::ParameterSetDescription>("protons_single_rp", protons_single_rp)
0132       ->setComment("settings for single-RP protons");
0133
0134   edm::ParameterSetDescription protons_multi_rp;
0135   protons_multi_rp.add<bool>("include", true)->setComment("flag whether multi-RP protons should be processed");
0136   protons_multi_rp.add<edm::InputTag>("input_tag", edm::InputTag("ctppsProtons", "multiRP"))->setComment("input tag");
0137   protons_multi_rp.add<std::string>("output_label", "multiRP")->setComment("output label");
0138   protons_multi_rp.add<bool>("check_valid_fit", true)->setComment("flag whether validFit should be checked");
0139   protons_multi_rp.add<double>("chi_sq_max", 1E-4)->setComment("maximum tolerated value of chi square");
0140   protons_multi_rp.add<double>("normalised_chi_sq_max", 1E100)
0141       ->setComment("maximum tolerated value of chi square / ndof, applied only if ndof > 0");
0142   desc.add<edm::ParameterSetDescription>("protons_multi_rp", protons_multi_rp)
0143       ->setComment("settings for multi-RP protons");
0144
0145   descriptions.add("ppsFilteredProtonProducer", desc);
0146 }
0147
0148 //----------------------------------------------------------------------------------------------------
0149
0150 bool PPSFilteredProtonProducer::IsTrackOK(const CTPPSLocalTrackLite &tr, unsigned int idx, std::ostringstream &log) {
0151   bool ok = true;
0152
0153   // checks for all tracks
0154   ok &= (std::abs(tr.tx()) < tracks_all_local_angle_x_max_);
0155   ok &= (std::abs(tr.ty()) < tracks_all_local_angle_y_max_);
0156
0157   // pixel checks
0158   const CTPPSDetId rpId(tr.rpId());
0159   if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0160     ok &= (find(tracks_pixel_forbidden_RecoInfo_values_.begin(),
0161                 tracks_pixel_forbidden_RecoInfo_values_.end(),
0162                 (unsigned int)tr.pixelTrackRecoInfo()) == tracks_pixel_forbidden_RecoInfo_values_.end());
0163     ok &= (tr.numberOfPointsUsedForFit() >= tracks_pixel_number_of_hits_min_);
0164     ok &= (tr.chiSquaredOverNDF() <= tracks_pixel_normalised_chi_sq_max_);
0165   }
0166
0167   if (!ok && verbosity_)
0168     log << "track idx=" << idx << " does not fulfil criteria." << std::endl;
0169
0170   return ok;
0171 }
0172
0173 //----------------------------------------------------------------------------------------------------
0174
0175 void PPSFilteredProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0176   std::ostringstream ssLog;
0177
0178   // process single-RP protons
0179   if (protons_single_rp_include_) {
0180     reco::ForwardProtonCollection const &hInputProtons = iEvent.get(protons_single_rp_input_token_);
0181     std::unique_ptr<reco::ForwardProtonCollection> pOutputProtons(new reco::ForwardProtonCollection);
0182
0183     for (const auto &proton : hInputProtons) {
0184       bool keep = true;
0185
0186       // no specific checks for single-RP protons
0187
0188       // test contributing tracks
0189       for (const auto &tr_ref : proton.contributingLocalTracks()) {
0190         if (!keep)
0191           break;
0192
0193         keep &= IsTrackOK(*tr_ref, tr_ref.key(), ssLog);
0194       }
0195
0196       n_protons_single_rp_all++;
0197
0198       if (keep) {
0199         n_protons_single_rp_kept++;
0200         pOutputProtons->push_back(proton);
0201       } else {
0202         if (verbosity_)
0203           ssLog << "single-RP proton idx=" << n_protons_single_rp_all - 1 << " excluded." << std::endl;
0204       }
0205     }
0206
0207     iEvent.put(std::move(pOutputProtons), protons_single_rp_output_label_);
0208   }
0209
0210   // process multi-RP protons
0211   if (protons_multi_rp_include_) {
0212     reco::ForwardProtonCollection const &hInputProtons = iEvent.get(protons_multi_rp_input_token_);
0213     std::unique_ptr<reco::ForwardProtonCollection> pOutputProtons(new reco::ForwardProtonCollection);
0214
0215     for (const auto &proton : hInputProtons) {
0216       bool keep = true;
0217
0218       // multi-RP proton checks
0219       if (protons_multi_rp_check_valid_fit_)
0220         keep &= proton.validFit();
0221
0222       keep &= (proton.chi2() <= protons_multi_rp_chi_sq_max_);
0223
0224       if (proton.ndof() > 0)
0225         keep &= (proton.normalizedChi2() <= protons_multi_rp_normalised_chi_sq_max_);
0226
0227       // test contributing tracks
0228       for (const auto &tr_ref : proton.contributingLocalTracks()) {
0229         if (!keep)
0230           break;
0231
0232         keep &= IsTrackOK(*tr_ref, tr_ref.key(), ssLog);
0233       }
0234
0235       n_protons_multi_rp_all++;
0236
0237       if (keep) {
0238         n_protons_multi_rp_kept++;
0239         pOutputProtons->push_back(proton);
0240       } else {
0241         if (verbosity_)
0242           ssLog << "multi-RP proton idx=" << n_protons_multi_rp_all - 1 << " excluded." << std::endl;
0243       }
0244     }
0245
0246     iEvent.put(std::move(pOutputProtons), protons_multi_rp_output_label_);
0247   }
0248
0249   if (verbosity_ && !ssLog.str().empty())
0250     edm::LogInfo("PPS") << ssLog.str();
0251 }
0252
0253 //----------------------------------------------------------------------------------------------------
0254
0255 void PPSFilteredProtonProducer::endStream() {
0256   edm::LogInfo("PPS")
0257       << "single-RP protons: total=" << n_protons_single_rp_all << ", kept=" << n_protons_single_rp_kept
0258       << " --> keep rate="
0259       << ((n_protons_single_rp_all > 0) ? double(n_protons_single_rp_kept) / n_protons_single_rp_all * 100. : 0.)
0260       << "%\n"
0261       << "multi-RP protons: total=" << n_protons_multi_rp_all << ", kept=" << n_protons_multi_rp_kept
0262       << " --> keep rate="
0263       << ((n_protons_multi_rp_all > 0) ? double(n_protons_multi_rp_kept) / n_protons_multi_rp_all * 100. : 0.) << "%";
0264 }
0265
0266 //----------------------------------------------------------------------------------------------------
0267
0268 DEFINE_FWK_MODULE(PPSFilteredProtonProducer);