CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
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_;  // map: arm, n(particles in acceptance) --> 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   // book plots
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   // get conditions
0135   const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0136
0137   // get input
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   // buffer for particle information
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;  // barcode --> info
0163
0164   // process HepMC
0165   for (auto it = hepMCEventAfterSmearing->particles_begin(); it != hepMCEventAfterSmearing->particles_end(); ++it) {
0166     const auto &part = *it;
0167
0168     // accept only stable non-beam protons
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   // check acceptance
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   // count particles in acceptance
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   // count reconstructed tracks
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   // count reconstructed protons
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   // fill plots
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     // skip event if no track in global acceptance
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     // skip events with some local reconstruction inefficiency
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);