File indexing completed on 2021-02-14 23:31:33
0001
0002
0003
0004
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012
0013 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0014
0015 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0016 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0017
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0019
0020 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0021 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0022 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0023 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0024 #include "DataFormats/Common/interface/DetSetVector.h"
0025 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0026 #include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h"
0027
0028 #include "TFile.h"
0029 #include "TH1D.h"
0030 #include "TH2D.h"
0031 #include "TProfile.h"
0032 #include "TGraphErrors.h"
0033
0034 #include <map>
0035 #include <string>
0036
0037
0038
0039 class CTPPSProtonReconstructionEfficiencyEstimatorMC : public edm::one::EDAnalyzer<> {
0040 public:
0041 explicit CTPPSProtonReconstructionEfficiencyEstimatorMC(const edm::ParameterSet &);
0042
0043 private:
0044 void analyze(const edm::Event &, const edm::EventSetup &) override;
0045 void endJob() override;
0046
0047 edm::EDGetTokenT<edm::HepMCProduct> tokenHepMCAfterSmearing_;
0048
0049 edm::EDGetTokenT<std::map<int, edm::DetSetVector<TotemRPRecHit>>> tokenStripRecHitsPerParticle_;
0050 edm::EDGetTokenT<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>> tokenPixelRecHitsPerParticle_;
0051
0052 edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0053
0054 edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
0055
0056 edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoESToken_;
0057
0058 unsigned int rpId_45_N_, rpId_45_F_;
0059 unsigned int rpId_56_N_, rpId_56_F_;
0060
0061 std::map<unsigned int, unsigned int> rpDecId_near_, rpDecId_far_;
0062
0063 std::string outputFile_;
0064
0065 unsigned int verbosity_;
0066
0067 struct PlotGroup {
0068 std::unique_ptr<TProfile> p_eff_vs_xi;
0069
0070 std::unique_ptr<TH1D> h_n_part_acc_nr, h_n_part_acc_fr;
0071
0072 PlotGroup()
0073 : p_eff_vs_xi(new TProfile("", ";#xi_{simu};efficiency", 19, 0.015, 0.205)),
0074 h_n_part_acc_nr(new TH1D("", ";n particles in acceptance", 6, -0.5, +5.5)),
0075 h_n_part_acc_fr(new TH1D("", ";n particles in acceptance", 6, -0.5, +5.5)) {}
0076
0077 void write() const {
0078 p_eff_vs_xi->Write("p_eff_vs_xi");
0079 h_n_part_acc_nr->Write("h_n_part_acc_nr");
0080 h_n_part_acc_fr->Write("h_n_part_acc_fr");
0081 }
0082 };
0083
0084 std::map<unsigned int, std::map<unsigned int, PlotGroup>> plots_;
0085 };
0086
0087
0088
0089 using namespace std;
0090 using namespace edm;
0091 using namespace HepMC;
0092
0093
0094
0095 CTPPSProtonReconstructionEfficiencyEstimatorMC::CTPPSProtonReconstructionEfficiencyEstimatorMC(
0096 const edm::ParameterSet &iConfig)
0097 : tokenHepMCAfterSmearing_(
0098 consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCAfterSmearing"))),
0099 tokenStripRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<TotemRPRecHit>>>(
0100 iConfig.getParameter<edm::InputTag>("tagStripRecHitsPerParticle"))),
0101 tokenPixelRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>>(
0102 iConfig.getParameter<edm::InputTag>("tagPixelRecHitsPerParticle"))),
0103 tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0104 tokenRecoProtonsMultiRP_(
0105 consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
0106 lhcInfoESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0107
0108 rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
0109 rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
0110 rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
0111 rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
0112
0113 outputFile_(iConfig.getParameter<string>("outputFile")),
0114
0115 verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)) {
0116 rpDecId_near_[0] = rpId_45_N_;
0117 rpDecId_far_[0] = rpId_45_F_;
0118
0119 rpDecId_near_[1] = rpId_56_N_;
0120 rpDecId_far_[1] = rpId_56_F_;
0121
0122
0123 for (unsigned int arm = 0; arm < 2; ++arm) {
0124 for (unsigned int np = 1; np <= 5; ++np)
0125 plots_[arm][np] = PlotGroup();
0126 }
0127 }
0128
0129
0130
0131 void CTPPSProtonReconstructionEfficiencyEstimatorMC::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0132 std::ostringstream os;
0133
0134
0135 const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0136
0137
0138 edm::Handle<edm::HepMCProduct> hHepMCAfterSmearing;
0139 iEvent.getByToken(tokenHepMCAfterSmearing_, hHepMCAfterSmearing);
0140 HepMC::GenEvent *hepMCEventAfterSmearing = (HepMC::GenEvent *)hHepMCAfterSmearing->GetEvent();
0141
0142 edm::Handle<std::map<int, edm::DetSetVector<TotemRPRecHit>>> hStripRecHitsPerParticle;
0143 iEvent.getByToken(tokenStripRecHitsPerParticle_, hStripRecHitsPerParticle);
0144
0145 edm::Handle<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>> hPixelRecHitsPerParticle;
0146 iEvent.getByToken(tokenPixelRecHitsPerParticle_, hPixelRecHitsPerParticle);
0147
0148 edm::Handle<CTPPSLocalTrackLiteCollection> tracks;
0149 iEvent.getByToken(tracksToken_, tracks);
0150
0151 Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0152 iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0153
0154
0155 struct ParticleInfo {
0156 unsigned int arm = 2;
0157 double xi = 0.;
0158 std::map<unsigned int, unsigned int> recHitsPerRP;
0159 bool inAcceptanceNear = false, inAcceptanceFar = false, inAcceptance = false;
0160 };
0161
0162 std::map<int, ParticleInfo> particleInfo;
0163
0164
0165 for (auto it = hepMCEventAfterSmearing->particles_begin(); it != hepMCEventAfterSmearing->particles_end(); ++it) {
0166 const auto &part = *it;
0167
0168
0169 if (part->pdg_id() != 2212)
0170 continue;
0171
0172 if (part->status() != 1)
0173 continue;
0174
0175 if (part->is_beam())
0176 continue;
0177
0178 const auto &mom = part->momentum();
0179
0180 if (mom.e() < 4500.)
0181 continue;
0182
0183 ParticleInfo info;
0184
0185 info.arm = (mom.z() > 0.) ? 0 : 1;
0186
0187 const double p_nom = lhcInfo.energy();
0188 info.xi = (p_nom - mom.rho()) / p_nom;
0189
0190 particleInfo[part->barcode()] = std::move(info);
0191 }
0192
0193
0194 for (const auto &pp : *hStripRecHitsPerParticle) {
0195 const auto barcode = pp.first;
0196
0197 for (const auto &ds : pp.second) {
0198 CTPPSDetId detId(ds.detId());
0199 CTPPSDetId rpId = detId.rpId();
0200 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
0201 }
0202 }
0203
0204 for (const auto &pp : *hPixelRecHitsPerParticle) {
0205 const auto barcode = pp.first;
0206
0207 for (const auto &ds : pp.second) {
0208 CTPPSDetId detId(ds.detId());
0209 CTPPSDetId rpId = detId.rpId();
0210 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
0211 }
0212 }
0213
0214 std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
0215
0216 for (auto &pp : particleInfo) {
0217 if (verbosity_)
0218 os << "* barcode=" << pp.first << ", arm=" << pp.second.arm << ", xi=" << pp.second.xi << std::endl;
0219
0220 for (const auto &rpp : pp.second.recHitsPerRP) {
0221 CTPPSDetId rpId(rpp.first);
0222 unsigned int needed_rec_hits = 1000;
0223 if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip)
0224 needed_rec_hits = 6;
0225 if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel)
0226 needed_rec_hits = 3;
0227
0228 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0229
0230 if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0231 if (rpDecId == rpDecId_near_[rpId.arm()])
0232 isStripRPNear[rpId.arm()] = true;
0233 if (rpDecId == rpDecId_far_[rpId.arm()])
0234 isStripRPFar[rpId.arm()] = true;
0235 }
0236
0237 if (rpp.second >= needed_rec_hits) {
0238 if (rpDecId == rpDecId_near_[rpId.arm()])
0239 pp.second.inAcceptanceNear = true;
0240 if (rpDecId == rpDecId_far_[rpId.arm()])
0241 pp.second.inAcceptanceFar = true;
0242 }
0243
0244 if (verbosity_)
0245 os << " RP " << rpDecId << ": " << rpp.second << " hits" << std::endl;
0246 }
0247
0248 pp.second.inAcceptance = pp.second.inAcceptanceNear && pp.second.inAcceptanceFar;
0249
0250 if (verbosity_)
0251 os << " inAcceptance: near=" << pp.second.inAcceptanceNear << ", far=" << pp.second.inAcceptanceFar
0252 << ", global=" << pp.second.inAcceptance << std::endl;
0253 }
0254
0255
0256 struct ArmCounter {
0257 unsigned int near = 0, far = 0, global = 0;
0258 };
0259 std::map<unsigned int, ArmCounter> nParticlesInAcceptance;
0260 for (auto &pp : particleInfo) {
0261 auto &counter = nParticlesInAcceptance[pp.second.arm];
0262 if (pp.second.inAcceptanceNear)
0263 counter.near++;
0264 if (pp.second.inAcceptanceFar)
0265 counter.far++;
0266 if (pp.second.inAcceptance)
0267 counter.global++;
0268 }
0269
0270
0271 std::map<unsigned int, ArmCounter> nReconstructedTracks;
0272 for (const auto &tr : *tracks) {
0273 CTPPSDetId rpId(tr.rpId());
0274 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0275
0276 if (rpDecId == rpDecId_near_[rpId.arm()])
0277 nReconstructedTracks[rpId.arm()].near++;
0278 if (rpDecId == rpDecId_far_[rpId.arm()])
0279 nReconstructedTracks[rpId.arm()].far++;
0280 }
0281
0282
0283 std::map<unsigned int, unsigned int> nReconstructedProtons;
0284 for (const auto &pr : *hRecoProtonsMultiRP) {
0285 if (!pr.validFit())
0286 continue;
0287
0288 unsigned int arm = 2;
0289 if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector45)
0290 arm = 0;
0291 if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector56)
0292 arm = 1;
0293
0294 nReconstructedProtons[arm]++;
0295 }
0296
0297
0298 for (unsigned int arm = 0; arm < 2; arm++) {
0299 const auto &npa = nParticlesInAcceptance[arm];
0300 const auto &nrt = nReconstructedTracks[arm];
0301
0302 if (verbosity_)
0303 os << "* arm " << arm << ": nRecoProtons=" << nReconstructedProtons[arm]
0304 << " (tracks near=" << nReconstructedTracks[arm].near << ", far=" << nReconstructedTracks[arm].far
0305 << "), nAcc=" << npa.global << " (near=" << npa.near << ", far=" << npa.far << ")" << std::endl;
0306
0307
0308 if (npa.global < 1)
0309 continue;
0310
0311 const auto &p = plots_[arm][npa.global];
0312
0313 p.h_n_part_acc_nr->Fill(npa.near);
0314 p.h_n_part_acc_fr->Fill(npa.far);
0315
0316
0317 if (nrt.near != npa.near || nrt.far != npa.far)
0318 continue;
0319
0320 const double eff = double(nReconstructedProtons[arm]) / npa.global;
0321
0322 if (verbosity_)
0323 os << " eff=" << eff << std::endl;
0324
0325 for (auto &pp : particleInfo) {
0326 if (pp.second.arm != arm || !pp.second.inAcceptance)
0327 continue;
0328
0329 p.p_eff_vs_xi->Fill(pp.second.xi, eff);
0330 }
0331 }
0332
0333 if (verbosity_)
0334 edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
0335 }
0336
0337
0338
0339 void CTPPSProtonReconstructionEfficiencyEstimatorMC::endJob() {
0340 auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0341
0342 for (const auto &ait : plots_) {
0343 char buf[100];
0344 sprintf(buf, "arm%u", ait.first);
0345 TDirectory *d_arm = f_out->mkdir(buf);
0346
0347 for (const auto &npit : ait.second) {
0348 sprintf(buf, "%u", npit.first);
0349 TDirectory *d_np = d_arm->mkdir(buf);
0350 gDirectory = d_np;
0351
0352 npit.second.write();
0353 }
0354 }
0355 }
0356
0357
0358
0359 DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorMC);