Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:27:03

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  *   Mauricio Thiel
0005  ****************************************************************************/
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/ESWatcher.h"
0014 
0015 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0016 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0017 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0018 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0019 #include "CondFormats/DataRecord/interface/PPSAssociationCutsRcd.h"
0020 #include "CondFormats/PPSObjects/interface/PPSAssociationCuts.h"
0021 
0022 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0023 
0024 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0025 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0026 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0027 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0028 
0029 #include "TFile.h"
0030 #include "TH1D.h"
0031 #include "TH2D.h"
0032 #include "TProfile.h"
0033 #include "TF1.h"
0034 #include "TGraph.h"
0035 
0036 #include <map>
0037 #include <string>
0038 #include <sstream>
0039 #include <utility>
0040 
0041 //----------------------------------------------------------------------------------------------------
0042 
0043 class CTPPSProtonReconstructionEfficiencyEstimatorData : public edm::one::EDAnalyzer<> {
0044 public:
0045   explicit CTPPSProtonReconstructionEfficiencyEstimatorData(const edm::ParameterSet &);
0046 
0047   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0048 
0049 private:
0050   void analyze(const edm::Event &, const edm::EventSetup &) override;
0051   void endJob() override;
0052 
0053   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tokenTracks_;
0054   edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
0055 
0056   edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoESToken_;
0057   edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsESToken_;
0058   edm::ESGetToken<PPSAssociationCuts, PPSAssociationCutsRcd> ppsAssociationCutsToken_;
0059 
0060   bool pixelDiscardBXShiftedTracks_;
0061 
0062   double localAngleXMin_, localAngleXMax_, localAngleYMin_, localAngleYMax_;
0063 
0064   unsigned int n_prep_events_;
0065 
0066   unsigned int n_exp_prot_max_;
0067 
0068   std::vector<double> n_sigmas_;
0069 
0070   std::string outputFile_;
0071 
0072   unsigned int verbosity_;
0073 
0074   edm::ESWatcher<CTPPSInterpolatedOpticsRcd> opticsWatcher_;
0075 
0076   struct ArmData {
0077     unsigned int rpId_N, rpId_F;
0078 
0079     std::shared_ptr<const TSpline3> s_x_to_xi_N;
0080 
0081     unsigned int n_events_with_tracks;
0082 
0083     std::unique_ptr<TH1D> h_de_x, h_de_y;
0084     std::unique_ptr<TH2D> h2_de_y_vs_de_x;
0085 
0086     double de_x_mean, de_x_sigma;
0087     double de_y_mean, de_y_sigma;
0088 
0089     std::vector<double> n_sigmas;
0090 
0091     unsigned int n_exp_prot_max;
0092 
0093     struct EffPlots {
0094       std::unique_ptr<TProfile> p_eff1_vs_x_N;
0095       std::unique_ptr<TProfile> p_eff1_vs_xi_N;
0096 
0097       std::unique_ptr<TProfile> p_eff2_vs_x_N;
0098       std::unique_ptr<TProfile> p_eff2_vs_xi_N;
0099 
0100       std::unique_ptr<TProfile> p_fr_match_unique, p_fr_match_non_unique;
0101       std::unique_ptr<TProfile> p_fr_signature_0, p_fr_signature_1, p_fr_signature_2, p_fr_signature_12;
0102 
0103       EffPlots()
0104           : p_eff1_vs_x_N(new TProfile("", ";x_{N}   (mm);efficiency", 50, 0., 25.)),
0105             p_eff1_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)),
0106 
0107             p_eff2_vs_x_N(new TProfile("", ";x_{N}   (mm);efficiency", 50, 0., 25.)),
0108             p_eff2_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)),
0109 
0110             p_fr_match_unique(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)),
0111             p_fr_match_non_unique(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)),
0112 
0113             p_fr_signature_0(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)),
0114             p_fr_signature_1(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)),
0115             p_fr_signature_2(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)),
0116             p_fr_signature_12(new TProfile("", ";x_{N}   (mm);fraction", 50, 0., 25.)) {}
0117 
0118       void write() const {
0119         p_eff1_vs_x_N->Write("p_eff1_vs_x_N");
0120         p_eff1_vs_xi_N->Write("p_eff1_vs_xi_N");
0121 
0122         p_eff2_vs_x_N->Write("p_eff2_vs_x_N");
0123         p_eff2_vs_xi_N->Write("p_eff2_vs_xi_N");
0124 
0125         p_fr_match_unique->Write("p_fr_match_unique");
0126         p_fr_match_non_unique->Write("p_fr_match_non_unique");
0127 
0128         p_fr_signature_0->Write("p_fr_signature_0");
0129         p_fr_signature_1->Write("p_fr_signature_1");
0130         p_fr_signature_2->Write("p_fr_signature_2");
0131         p_fr_signature_12->Write("p_fr_signature_12");
0132       }
0133     };
0134 
0135     // (n exp protons, index in n_sigmas) --> plots
0136     // n exp protons = 0 --> no condition (summary case)
0137     std::map<unsigned int, std::map<unsigned int, EffPlots>> effPlots;
0138 
0139     ArmData()
0140         : n_events_with_tracks(0),
0141 
0142           h_de_x(new TH1D("", ";x_{F} - x_{N}   (mm)", 200, -1., +1.)),
0143           h_de_y(new TH1D("", ";y_{F} - y_{N}   (mm)", 200, -1., +1.)),
0144           h2_de_y_vs_de_x(new TH2D("", ";x_{F} - x_{N}   (mm);y_{F} - y_{N}   (mm)", 100, -1., +1., 100, -1., +1.)),
0145 
0146           de_x_mean(0.),
0147           de_x_sigma(0.),
0148           de_y_mean(0.),
0149           de_y_sigma(0.) {
0150       for (unsigned int nep = 0; nep <= n_exp_prot_max; ++nep) {
0151         for (unsigned int nsi = 0; nsi < n_sigmas.size(); ++nsi) {
0152           effPlots[nep][nsi] = EffPlots();
0153         }
0154       }
0155     }
0156 
0157     void write() const {
0158       h_de_x->Write("h_de_x");
0159       h_de_y->Write("h_de_y");
0160       h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
0161 
0162       char buf[100];
0163 
0164       for (const auto &n_si : n_sigmas) {
0165         sprintf(buf, "g_de_y_vs_de_x_n_si=%.1f", n_si);
0166         TGraph *g = new TGraph();
0167         g->SetName(buf);
0168 
0169         g->SetPoint(0, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
0170         g->SetPoint(1, de_x_mean + n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
0171         g->SetPoint(2, de_x_mean + n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma);
0172         g->SetPoint(3, de_x_mean - n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma);
0173         g->SetPoint(4, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
0174 
0175         g->Write();
0176       }
0177 
0178       TDirectory *d_top = gDirectory;
0179 
0180       for (const auto &nep_p : effPlots) {
0181         if (nep_p.first == 0)
0182           sprintf(buf, "exp prot all");
0183         else
0184           sprintf(buf, "exp prot %u", nep_p.first);
0185 
0186         TDirectory *d_nep = d_top->mkdir(buf);
0187 
0188         for (const auto &nsi_p : nep_p.second) {
0189           sprintf(buf, "nsi = %.1f", n_sigmas[nsi_p.first]);
0190 
0191           TDirectory *d_nsi = d_nep->mkdir(buf);
0192 
0193           gDirectory = d_nsi;
0194 
0195           nsi_p.second.write();
0196         }
0197       }
0198 
0199       gDirectory = d_top;
0200     }
0201 
0202     void UpdateOptics(const LHCInterpolatedOpticalFunctionsSetCollection &ofc) {
0203       const LHCInterpolatedOpticalFunctionsSet *of = nullptr;
0204 
0205       for (const auto &p : ofc) {
0206         CTPPSDetId rpId(p.first);
0207         unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0208 
0209         if (decRPId == rpId_N) {
0210           of = &p.second;
0211           break;
0212         }
0213       }
0214 
0215       if (!of) {
0216         edm::LogError("ArmData::UpdateOptics") << "Cannot find optical functions for RP " << rpId_N;
0217         return;
0218       }
0219 
0220       std::vector<double> xiValues =
0221           of->getXiValues();  // local copy made since the TSpline constructor needs non-const parameters
0222       std::vector<double> xDValues = of->getFcnValues()[LHCOpticalFunctionsSet::exd];
0223       s_x_to_xi_N = std::make_shared<TSpline3>("", xDValues.data(), xiValues.data(), xiValues.size());
0224     }
0225   };
0226 
0227   std::map<unsigned int, ArmData> data_;
0228 
0229   std::unique_ptr<TF1> ff_;
0230 };
0231 
0232 //----------------------------------------------------------------------------------------------------
0233 
0234 using namespace std;
0235 using namespace edm;
0236 
0237 //----------------------------------------------------------------------------------------------------
0238 
0239 CTPPSProtonReconstructionEfficiencyEstimatorData::CTPPSProtonReconstructionEfficiencyEstimatorData(
0240     const edm::ParameterSet &iConfig)
0241     : tokenTracks_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0242 
0243       tokenRecoProtonsMultiRP_(
0244           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
0245 
0246       lhcInfoESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0247       opticsESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
0248       ppsAssociationCutsToken_(
0249           esConsumes(ESInputTag("", iConfig.getParameter<std::string>("ppsAssociationCutsLabel")))),
0250 
0251       pixelDiscardBXShiftedTracks_(iConfig.getParameter<bool>("pixelDiscardBXShiftedTracks")),
0252 
0253       localAngleXMin_(iConfig.getParameter<double>("localAngleXMin")),
0254       localAngleXMax_(iConfig.getParameter<double>("localAngleXMax")),
0255       localAngleYMin_(iConfig.getParameter<double>("localAngleYMin")),
0256       localAngleYMax_(iConfig.getParameter<double>("localAngleYMax")),
0257 
0258       n_prep_events_(iConfig.getParameter<unsigned int>("n_prep_events")),
0259       n_exp_prot_max_(iConfig.getParameter<unsigned int>("n_exp_prot_max")),
0260       n_sigmas_(iConfig.getParameter<std::vector<double>>("n_sigmas")),
0261 
0262       outputFile_(iConfig.getParameter<string>("outputFile")),
0263 
0264       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity")),
0265 
0266       ff_(new TF1("ff", "[0] * exp(- (x-[1])*(x-[1]) / 2 / [2]/[2]) + [4]")) {
0267   data_[0].n_exp_prot_max = n_exp_prot_max_;
0268   data_[0].n_sigmas = n_sigmas_;
0269   data_[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
0270   data_[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
0271 
0272   data_[1].n_exp_prot_max = n_exp_prot_max_;
0273   data_[1].n_sigmas = n_sigmas_;
0274   data_[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
0275   data_[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
0276 }
0277 
0278 //----------------------------------------------------------------------------------------------------
0279 
0280 void CTPPSProtonReconstructionEfficiencyEstimatorData::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0281   edm::ParameterSetDescription desc;
0282 
0283   desc.add<edm::InputTag>("tagTracks", edm::InputTag())->setComment("input tag for local lite tracks");
0284   desc.add<edm::InputTag>("tagRecoProtonsMultiRP", edm::InputTag())->setComment("input tag for multi-RP reco protons");
0285 
0286   desc.add<string>("lhcInfoLabel", "")->setComment("label of LHCInfo data");
0287   desc.add<string>("opticsLabel", "")->setComment("label of optics data");
0288   desc.add<string>("ppsAssociationCutsLabel", "")->setComment("label of PPSAssociationCuts data");
0289 
0290   desc.add<bool>("pixelDiscardBXShiftedTracks", false)
0291       ->setComment("whether to discard pixel tracks built from BX-shifted planes");
0292 
0293   desc.add<double>("localAngleXMin", -0.03)->setComment("minimal accepted value of local horizontal angle (rad)");
0294   desc.add<double>("localAngleXMax", +0.03)->setComment("maximal accepted value of local horizontal angle (rad)");
0295   desc.add<double>("localAngleYMin", -0.04)->setComment("minimal accepted value of local vertical angle (rad)");
0296   desc.add<double>("localAngleYMax", +0.04)->setComment("maximal accepted value of local vertical angle (rad)");
0297 
0298   desc.add<unsigned int>("n_prep_events", 1000)
0299       ->setComment("number of preparatory events (to determine de x and de y window)");
0300 
0301   desc.add<unsigned int>("n_exp_prot_max", 5)->setComment("maximum number of expected protons per event and per arm");
0302 
0303   desc.add<std::vector<double>>("n_sigmas", {3., 5., 7.})->setComment("list of n_sigma values");
0304 
0305   desc.add<unsigned int>("rpId_45_N", 0)->setComment("decimal RP id for 45 near");
0306   desc.add<unsigned int>("rpId_45_F", 0)->setComment("decimal RP id for 45 far");
0307   desc.add<unsigned int>("rpId_56_N", 0)->setComment("decimal RP id for 56 near");
0308   desc.add<unsigned int>("rpId_56_F", 0)->setComment("decimal RP id for 56 far");
0309 
0310   desc.add<std::string>("outputFile", "output.root")->setComment("output file name");
0311 
0312   desc.addUntracked<unsigned int>("verbosity", 0)->setComment("verbosity level");
0313 
0314   descriptions.add("ctppsProtonReconstructionEfficiencyEstimatorData", desc);
0315 }
0316 
0317 //----------------------------------------------------------------------------------------------------
0318 
0319 void CTPPSProtonReconstructionEfficiencyEstimatorData::analyze(const edm::Event &iEvent,
0320                                                                const edm::EventSetup &iSetup) {
0321   std::ostringstream os;
0322 
0323   // get conditions
0324   const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0325   const auto &opticalFunctions = iSetup.getData(opticsESToken_);
0326   const auto &ppsAssociationCuts = iSetup.getData(ppsAssociationCutsToken_);
0327 
0328   // check optics change
0329   if (opticsWatcher_.check(iSetup)) {
0330     data_[0].UpdateOptics(opticalFunctions);
0331     data_[1].UpdateOptics(opticalFunctions);
0332   }
0333 
0334   // get input
0335   edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
0336   iEvent.getByToken(tokenTracks_, hTracks);
0337 
0338   Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0339   iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0340 
0341   // pre-selection of tracks
0342   std::vector<unsigned int> sel_track_idc;
0343 
0344   std::map<unsigned int, std::vector<unsigned int>> trackingSelection;
0345 
0346   for (unsigned int idx = 0; idx < hTracks->size(); ++idx) {
0347     const auto &tr = hTracks->at(idx);
0348 
0349     if (tr.tx() < localAngleXMin_ || tr.tx() > localAngleXMax_ || tr.ty() < localAngleYMin_ ||
0350         tr.ty() > localAngleYMax_)
0351       continue;
0352 
0353     if (pixelDiscardBXShiftedTracks_) {
0354       if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
0355           tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
0356         continue;
0357     }
0358 
0359     sel_track_idc.push_back(idx);
0360 
0361     const CTPPSDetId rpId(tr.rpId());
0362 
0363     const bool trackerRP =
0364         (rpId.subdetId() == CTPPSDetId::sdTrackingStrip || rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0365 
0366     if (trackerRP)
0367       trackingSelection[rpId.arm()].push_back(idx);
0368   }
0369 
0370   // debug print
0371   if (verbosity_ > 1) {
0372     os << "* tracks (pre-selected):" << std::endl;
0373 
0374     for (const auto idx : sel_track_idc) {
0375       const auto &tr = hTracks->at(idx);
0376       CTPPSDetId rpId(tr.rpId());
0377       unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0378 
0379       os << "    [" << idx << "] RP=" << decRPId << ", x=" << tr.x() << ", y=" << tr.y() << std::endl;
0380     }
0381 
0382     os << "* protons:" << std::endl;
0383     for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
0384       const auto &pr = (*hRecoProtonsMultiRP)[i];
0385       os << "    [" << i << "] track indices: ";
0386       for (const auto &trr : pr.contributingLocalTracks())
0387         os << trr.key() << ", ";
0388       os << std::endl;
0389     }
0390   }
0391 
0392   // make de_x and de_y plots
0393   map<unsigned int, bool> hasTracksInArm;
0394 
0395   for (const auto idx_i : sel_track_idc) {
0396     const auto &tr_i = hTracks->at(idx_i);
0397     CTPPSDetId rpId_i(tr_i.rpId());
0398     unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
0399 
0400     for (const auto idx_j : sel_track_idc) {
0401       const auto &tr_j = hTracks->at(idx_j);
0402       CTPPSDetId rpId_j(tr_j.rpId());
0403       unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
0404 
0405       // check whether desired RP combination
0406       unsigned int arm = 123;
0407       for (const auto &ap : data_) {
0408         if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
0409           arm = ap.first;
0410       }
0411 
0412       if (arm > 1)
0413         continue;
0414 
0415       // update flag
0416       hasTracksInArm[arm] = true;
0417 
0418       // update plots
0419       auto &ad = data_[arm];
0420       const double de_x = tr_j.x() - tr_i.x();
0421       const double de_y = tr_j.y() - tr_i.y();
0422 
0423       if (ad.n_events_with_tracks < n_prep_events_) {
0424         ad.h_de_x->Fill(de_x);
0425         ad.h_de_y->Fill(de_y);
0426       }
0427 
0428       ad.h2_de_y_vs_de_x->Fill(de_x, de_y);
0429     }
0430   }
0431 
0432   // update event counter
0433   for (auto &ap : data_) {
0434     if (hasTracksInArm[ap.first])
0435       ap.second.n_events_with_tracks++;
0436   }
0437 
0438   // if threshold reached do fits
0439   for (auto &ap : data_) {
0440     auto &ad = ap.second;
0441 
0442     if (ad.n_events_with_tracks == n_prep_events_) {
0443       if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.)
0444         continue;
0445 
0446       std::vector<std::pair<unsigned int, TH1D *>> m;
0447       m.emplace_back(0, ad.h_de_x.get());
0448       m.emplace_back(1, ad.h_de_y.get());
0449 
0450       for (const auto &p : m) {
0451         double max_pos = -1E100, max_val = -1.;
0452         for (int bi = 1; bi < p.second->GetNbinsX(); ++bi) {
0453           const double pos = p.second->GetBinCenter(bi);
0454           const double val = p.second->GetBinContent(bi);
0455 
0456           if (val > max_val) {
0457             max_val = val;
0458             max_pos = pos;
0459           }
0460         }
0461 
0462         const double sig = 0.2;
0463 
0464         ff_->SetParameters(max_val, max_pos, sig, 0.);
0465         p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
0466         p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
0467 
0468         if (p.first == 0) {
0469           ad.de_x_mean = ff_->GetParameter(1);
0470           ad.de_x_sigma = fabs(ff_->GetParameter(2));
0471         }
0472         if (p.first == 1) {
0473           ad.de_y_mean = ff_->GetParameter(1);
0474           ad.de_y_sigma = fabs(ff_->GetParameter(2));
0475         }
0476       }
0477 
0478       if (verbosity_) {
0479         os << "* fitting arm " << ap.first << std::endl
0480            << "    de_x: mean = " << ad.de_x_mean << ", sigma = " << ad.de_x_sigma << std::endl
0481            << "    de_y: mean = " << ad.de_y_mean << ", sigma = " << ad.de_y_sigma;
0482       }
0483     }
0484   }
0485 
0486   // logic below:
0487   //    1) evaluate "proton candiates" or "expected protons" = local tracks in the near RP which have a compatible
0488   //       partner in the far RP; the compatibility is evaluated at multiple n_sigma choices
0489   //    2) check how often each "proton candidate" can be matched to a reconstructed proton
0490   //    3) re-evaluate association cuts --> check details for each "proton candidate"
0491 
0492   // data structures for efficiency analysis
0493 
0494   struct ArmEventData {
0495     // n_sigma --> list of "proton candidates" (local track indices in the near RP)
0496     std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
0497 
0498     // list of valid reco-proton indices (in the chosen arm)
0499     std::set<unsigned int> reco_proton_idc;
0500 
0501     // n_sigma --> list of "proton candicates" which have/don't have matching reco proton
0502     std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
0503 
0504     // for re-evaluation of association cuts
0505     struct EvaluatedPair {
0506       unsigned idx_N, idx_F;
0507       bool x_cut, y_cut;
0508       bool match = false;
0509       bool unique = false;
0510     };
0511 
0512     vector<EvaluatedPair> evaluatedPairs;
0513   };
0514 
0515   std::map<unsigned int, ArmEventData> armEventData;
0516 
0517   // determine the number of proton candiates or expected protons
0518   for (const auto idx_i : sel_track_idc) {
0519     const auto &tr_i = hTracks->at(idx_i);
0520     CTPPSDetId rpId_i(tr_i.rpId());
0521     unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
0522 
0523     for (const auto idx_j : sel_track_idc) {
0524       const auto &tr_j = hTracks->at(idx_j);
0525       CTPPSDetId rpId_j(tr_j.rpId());
0526       unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
0527 
0528       // check whether desired RP combination
0529       unsigned int arm = 123;
0530       for (const auto &ap : data_) {
0531         if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
0532           arm = ap.first;
0533       }
0534 
0535       if (arm > 1)
0536         continue;
0537 
0538       // check near-far matching
0539       auto &ad = data_[arm];
0540       const double de_x = tr_j.x() - tr_i.x();
0541       const double de_y = tr_j.y() - tr_i.y();
0542 
0543       for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
0544         const double n_si = ad.n_sigmas[nsi];
0545         const bool match_x = fabs(de_x - ad.de_x_mean) < n_si * ad.de_x_sigma;
0546         const bool match_y = fabs(de_y - ad.de_y_mean) < n_si * ad.de_y_sigma;
0547         if (match_x && match_y)
0548           armEventData[arm].matched_track_idc[nsi].insert(idx_i);
0549       }
0550     }
0551   }
0552 
0553   // determine the number of reconstructed protons
0554   for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
0555     const auto &proton = (*hRecoProtonsMultiRP)[i];
0556 
0557     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
0558     unsigned int arm = rpId.arm();
0559 
0560     if (proton.validFit())
0561       armEventData[arm].reco_proton_idc.insert(i);
0562   }
0563 
0564   // compare matched tracks with reco protons
0565   if (verbosity_ > 1)
0566     os << "* cmp matched tracks vs. reco protons" << std::endl;
0567 
0568   for (auto &ap : armEventData) {
0569     auto &ad = data_[ap.first];
0570 
0571     if (verbosity_ > 1)
0572       os << "    arm " << ap.first << std::endl;
0573 
0574     for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
0575       if (verbosity_ > 1)
0576         os << "        nsi = " << nsi << std::endl;
0577 
0578       for (const auto &tri : ap.second.matched_track_idc[nsi]) {
0579         const auto &track = hTracks->at(tri);
0580 
0581         bool some_proton_matching = false;
0582 
0583         if (verbosity_ > 1)
0584           os << "            tri = " << tri << std::endl;
0585 
0586         for (const auto &pri : ap.second.reco_proton_idc) {
0587           const auto &proton = (*hRecoProtonsMultiRP)[pri];
0588 
0589           bool proton_matching = false;
0590 
0591           for (const auto &pr_tr : proton.contributingLocalTracks()) {
0592             CTPPSDetId rpId(pr_tr->rpId());
0593             unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0594 
0595             if (decRPId != ad.rpId_N)
0596               continue;
0597 
0598             const double x = pr_tr->x();
0599             const double y = pr_tr->y();
0600             const double th = 1E-3;  // 1 um
0601 
0602             const bool match = (fabs(x - track.x()) < th) && (fabs(y - track.y()) < th);
0603 
0604             if (verbosity_ > 1)
0605               os << "                pri = " << pri << ": x_tr = " << track.x() << ", x_pr = " << x
0606                  << ", match = " << match << std::endl;
0607 
0608             if (match) {
0609               proton_matching = true;
0610               break;
0611             }
0612           }
0613 
0614           if (proton_matching) {
0615             some_proton_matching = true;
0616             break;
0617           }
0618         }
0619 
0620         if (verbosity_ > 1)
0621           os << "            --> some_proton_matching " << some_proton_matching << std::endl;
0622 
0623         if (some_proton_matching)
0624           ap.second.matched_track_with_prot_idc[nsi].insert(tri);
0625         else
0626           ap.second.matched_track_without_prot_idc[nsi].insert(tri);
0627       }
0628     }
0629   }
0630 
0631   // redo association cuts
0632   for (auto &ap : armEventData) {
0633     const auto &arm = ap.first;
0634 
0635     auto &ad = data_[arm];
0636 
0637     auto &aed = ap.second;
0638 
0639     auto &ass_cut = ppsAssociationCuts.getAssociationCuts(arm);
0640 
0641     const auto &indices = trackingSelection[arm];
0642 
0643     map<unsigned int, unsigned> matching_multiplicity;
0644 
0645     for (const auto &i : indices) {
0646       for (const auto &j : indices) {
0647         const auto &tr_i = hTracks->at(i);
0648         const auto &tr_j = hTracks->at(j);
0649 
0650         const CTPPSDetId id_i(tr_i.rpId());
0651         const CTPPSDetId id_j(tr_j.rpId());
0652 
0653         const unsigned rpDecId_i = id_i.arm() * 100 + id_i.station() * 10 + id_i.rp();
0654         const unsigned rpDecId_j = id_j.arm() * 100 + id_j.station() * 10 + id_j.rp();
0655 
0656         if (rpDecId_i != ad.rpId_N || rpDecId_j != ad.rpId_F)
0657           continue;
0658 
0659         ArmEventData::EvaluatedPair evp;
0660         evp.idx_N = i;
0661         evp.idx_F = j;
0662 
0663         evp.x_cut = ass_cut.isSatisfied(ass_cut.qX, tr_i.x(), tr_i.y(), lhcInfo.crossingAngle(), tr_i.x() - tr_j.x());
0664         evp.y_cut = ass_cut.isSatisfied(ass_cut.qY, tr_i.x(), tr_i.y(), lhcInfo.crossingAngle(), tr_i.y() - tr_j.y());
0665 
0666         evp.match = evp.x_cut && evp.y_cut;
0667 
0668         aed.evaluatedPairs.push_back(evp);
0669 
0670         if (evp.match) {
0671           matching_multiplicity[i]++;
0672           matching_multiplicity[j]++;
0673         }
0674       }
0675     }
0676 
0677     for (auto &evp : aed.evaluatedPairs)
0678       evp.unique = (matching_multiplicity[evp.idx_N] == 1) && (matching_multiplicity[evp.idx_F] == 1);
0679   }
0680 
0681   // debug print
0682   if (verbosity_ > 1) {
0683     for (auto &ap : armEventData) {
0684       auto &ad = data_[ap.first];
0685 
0686       if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
0687         continue;
0688 
0689       os << "* results for arm " << ap.first << std::endl;
0690 
0691       os << "    reco_proton_idc: ";
0692       for (const auto &idx : ap.second.reco_proton_idc)
0693         os << idx << ", ";
0694       os << std::endl;
0695 
0696       for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
0697         os << "    n_si = " << ad.n_sigmas[nsi] << std::endl;
0698 
0699         os << "        matched_track_idc: ";
0700         for (const auto &idx : ap.second.matched_track_idc[nsi])
0701           os << idx << ", ";
0702         os << std::endl;
0703 
0704         os << "        matched_track_with_prot_idc: ";
0705         for (const auto &idx : ap.second.matched_track_with_prot_idc[nsi])
0706           os << idx << ", ";
0707         os << std::endl;
0708 
0709         os << "        matched_track_without_prot_idc: ";
0710         for (const auto &idx : ap.second.matched_track_without_prot_idc[nsi])
0711           os << idx << ", ";
0712         os << std::endl;
0713       }
0714 
0715       os << "    evaluated pairs: ";
0716       for (const auto &p : ap.second.evaluatedPairs)
0717         os << "(" << p.idx_N << "-" << p.idx_F << ",M" << p.match << ",U" << p.unique << ")"
0718            << ", ";
0719       os << std::endl;
0720     }
0721   }
0722 
0723   // update efficiency plots
0724   for (auto &ap : armEventData) {
0725     const auto &arm = ap.first;
0726     auto &ad = data_[arm];
0727 
0728     // stop if sigmas not yet determined
0729     if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
0730       continue;
0731 
0732     // loop over n_sigma choices
0733     for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
0734       const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size();
0735       const unsigned int n_rec_prot = ap.second.reco_proton_idc.size();
0736 
0737       // stop if N(expected protons) out of range
0738       if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
0739         continue;
0740 
0741       // update method 1 plots
0742       const double eff = double(n_rec_prot) / n_exp_prot;
0743 
0744       for (unsigned int tri : ap.second.matched_track_idc[nsi]) {
0745         const double x_N = hTracks->at(tri).x();
0746         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
0747 
0748         ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
0749         ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
0750 
0751         ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
0752         ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
0753       }
0754 
0755       // update method 2 plots
0756       for (const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) {
0757         const double x_N = hTracks->at(tri).x();
0758         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
0759 
0760         ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
0761         ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
0762 
0763         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
0764         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
0765       }
0766 
0767       for (const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) {
0768         const double x_N = hTracks->at(tri).x();
0769         const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
0770 
0771         ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
0772         ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
0773 
0774         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
0775         ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
0776       }
0777 
0778       // update association cut plots
0779       for (const auto &cand : ap.second.matched_track_idc[nsi])  // loop over candidates
0780       {
0781         const double x_N = hTracks->at(cand).x();
0782 
0783         auto &plots = ad.effPlots[n_exp_prot][nsi];
0784 
0785         bool hasMatchingAndUniquePair = false;
0786         bool hasMatchingAndNonUniquePair = false;
0787         for (const auto &pair : ap.second.evaluatedPairs)  // loop over pairs
0788         {
0789           // stop if cand not compatible with pair
0790           if (cand != pair.idx_N)
0791             continue;
0792 
0793           if (pair.match && pair.unique)
0794             hasMatchingAndUniquePair = true;
0795 
0796           if (pair.match && !pair.unique)
0797             hasMatchingAndNonUniquePair = true;
0798 
0799           unsigned int signature = 999999;
0800           if (!pair.x_cut && !pair.y_cut)
0801             signature = 0;
0802           if (pair.x_cut && !pair.y_cut)
0803             signature = 1;
0804           if (!pair.x_cut && pair.y_cut)
0805             signature = 2;
0806           if (pair.x_cut && pair.y_cut)
0807             signature = 12;
0808 
0809           plots.p_fr_signature_0->Fill(x_N, (signature == 0) ? 1 : 0);
0810           plots.p_fr_signature_1->Fill(x_N, (signature == 1) ? 1 : 0);
0811           plots.p_fr_signature_2->Fill(x_N, (signature == 2) ? 1 : 0);
0812           plots.p_fr_signature_12->Fill(x_N, (signature == 12) ? 1 : 0);
0813         }
0814 
0815         plots.p_fr_match_unique->Fill(x_N, (hasMatchingAndUniquePair) ? 1 : 0);
0816         plots.p_fr_match_non_unique->Fill(x_N, (hasMatchingAndNonUniquePair) ? 1 : 0);
0817       }
0818     }
0819   }
0820 
0821   if (verbosity_)
0822     edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str();
0823 }
0824 
0825 //----------------------------------------------------------------------------------------------------
0826 
0827 void CTPPSProtonReconstructionEfficiencyEstimatorData::endJob() {
0828   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0829 
0830   for (const auto &ait : data_) {
0831     char buf[100];
0832     sprintf(buf, "arm %u", ait.first);
0833     TDirectory *d_arm = f_out->mkdir(buf);
0834     gDirectory = d_arm;
0835 
0836     ait.second.write();
0837   }
0838 }
0839 
0840 //----------------------------------------------------------------------------------------------------
0841 
0842 DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorData);