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