File indexing completed on 2023-03-17 11:27:03
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 "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
0136
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();
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
0324 const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0325 const auto &opticalFunctions = iSetup.getData(opticsESToken_);
0326 const auto &ppsAssociationCuts = iSetup.getData(ppsAssociationCutsToken_);
0327
0328
0329 if (opticsWatcher_.check(iSetup)) {
0330 data_[0].UpdateOptics(opticalFunctions);
0331 data_[1].UpdateOptics(opticalFunctions);
0332 }
0333
0334
0335 edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
0336 iEvent.getByToken(tokenTracks_, hTracks);
0337
0338 Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0339 iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0340
0341
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
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
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
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
0416 hasTracksInArm[arm] = true;
0417
0418
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
0433 for (auto &ap : data_) {
0434 if (hasTracksInArm[ap.first])
0435 ap.second.n_events_with_tracks++;
0436 }
0437
0438
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
0487
0488
0489
0490
0491
0492
0493
0494 struct ArmEventData {
0495
0496 std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
0497
0498
0499 std::set<unsigned int> reco_proton_idc;
0500
0501
0502 std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
0503
0504
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
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
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
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
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
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;
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
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
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
0724 for (auto &ap : armEventData) {
0725 const auto &arm = ap.first;
0726 auto &ad = data_[arm];
0727
0728
0729 if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
0730 continue;
0731
0732
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
0738 if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
0739 continue;
0740
0741
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);
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
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);
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);
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
0779 for (const auto &cand : ap.second.matched_track_idc[nsi])
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)
0788 {
0789
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);
