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 "CondTools/RunInfo/interface/LHCInfoCombined.h"
0016
0017 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0018
0019 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0020 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0021 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0022 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0023 #include "DataFormats/Common/interface/DetSetVector.h"
0024 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0025 #include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h"
0026
0027 #include "TFile.h"
0028 #include "TH1D.h"
0029 #include "TH2D.h"
0030 #include "TProfile.h"
0031 #include "TGraphErrors.h"
0032
0033 #include <map>
0034 #include <string>
0035
0036 //----------------------------------------------------------------------------------------------------
0037
0038 class CTPPSProtonReconstructionEfficiencyEstimatorMC : public edm::one::EDAnalyzer<> {
0039 public:
0040   explicit CTPPSProtonReconstructionEfficiencyEstimatorMC(const edm::ParameterSet &);
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0043
0044 private:
0045   void analyze(const edm::Event &, const edm::EventSetup &) override;
0046   void endJob() override;
0047
0048   edm::EDGetTokenT<edm::HepMCProduct> tokenHepMCAfterSmearing_;
0049
0050   edm::EDGetTokenT<std::map<int, edm::DetSetVector<TotemRPRecHit>>> tokenStripRecHitsPerParticle_;
0051   edm::EDGetTokenT<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>> tokenPixelRecHitsPerParticle_;
0052
0053   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0054
0055   edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
0056
0057   const edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoToken_;
0058   const edm::ESGetToken<LHCInfoPerLS, LHCInfoPerLSRcd> lhcInfoPerLSToken_;
0059   const edm::ESGetToken<LHCInfoPerFill, LHCInfoPerFillRcd> lhcInfoPerFillToken_;
0060   const bool useNewLHCInfo_;
0061
0062   unsigned int rpId_45_N_, rpId_45_F_;
0063   unsigned int rpId_56_N_, rpId_56_F_;
0064
0065   std::map<unsigned int, unsigned int> rpDecId_near_, rpDecId_far_;
0066
0067   std::string outputFile_;
0068
0069   unsigned int verbosity_;
0070
0071   struct PlotGroup {
0072     std::unique_ptr<TProfile> p_eff_vs_xi;
0073
0074     std::unique_ptr<TH1D> h_n_part_acc_nr, h_n_part_acc_fr;
0075
0076     PlotGroup()
0077         : p_eff_vs_xi(new TProfile("", ";#xi_{simu};efficiency", 19, 0.015, 0.205)),
0078           h_n_part_acc_nr(new TH1D("", ";n particles in acceptance", 6, -0.5, +5.5)),
0079           h_n_part_acc_fr(new TH1D("", ";n particles in acceptance", 6, -0.5, +5.5)) {}
0080
0081     void write() const {
0082       p_eff_vs_xi->Write("p_eff_vs_xi");
0083       h_n_part_acc_nr->Write("h_n_part_acc_nr");
0084       h_n_part_acc_fr->Write("h_n_part_acc_fr");
0085     }
0086   };
0087
0088   std::map<unsigned int, std::map<unsigned int, PlotGroup>> plots_;  // map: arm, n(particles in acceptance) --> plots
0089 };
0090
0091 //----------------------------------------------------------------------------------------------------
0092
0093 using namespace std;
0094 using namespace edm;
0095 using namespace HepMC;
0096
0097 //----------------------------------------------------------------------------------------------------
0098
0099 CTPPSProtonReconstructionEfficiencyEstimatorMC::CTPPSProtonReconstructionEfficiencyEstimatorMC(
0100     const edm::ParameterSet &iConfig)
0101     : tokenHepMCAfterSmearing_(
0102           consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCAfterSmearing"))),
0103       tokenStripRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<TotemRPRecHit>>>(
0104           iConfig.getParameter<edm::InputTag>("tagStripRecHitsPerParticle"))),
0105       tokenPixelRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>>(
0106           iConfig.getParameter<edm::InputTag>("tagPixelRecHitsPerParticle"))),
0107       tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0108       tokenRecoProtonsMultiRP_(
0109           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
0110
0111       lhcInfoToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0112       lhcInfoPerLSToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerLSLabel")))),
0113       lhcInfoPerFillToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerFillLabel")))),
0114       useNewLHCInfo_(iConfig.getParameter<bool>("useNewLHCInfo")),
0115
0116       rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
0117       rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
0118       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
0119       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
0120
0121       outputFile_(iConfig.getParameter<string>("outputFile")),
0122
0123       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)) {
0124   rpDecId_near_[0] = rpId_45_N_;
0125   rpDecId_far_[0] = rpId_45_F_;
0126
0127   rpDecId_near_[1] = rpId_56_N_;
0128   rpDecId_far_[1] = rpId_56_F_;
0129
0130   // book plots
0131   for (unsigned int arm = 0; arm < 2; ++arm) {
0132     for (unsigned int np = 1; np <= 5; ++np)
0133       plots_[arm][np] = PlotGroup();
0134   }
0135 }
0136
0137 //----------------------------------------------------------------------------------------------------
0138
0139 void CTPPSProtonReconstructionEfficiencyEstimatorMC::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0140   edm::ParameterSetDescription desc;
0141
0142   desc.add<edm::InputTag>("tagTracks", edm::InputTag())->setComment("input tag for local lite tracks");
0143   desc.add<edm::InputTag>("tagRecoProtonsMultiRP", edm::InputTag())->setComment("input tag for multi-RP reco protons");
0144
0145   desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
0146   desc.add<std::string>("lhcInfoPerLSLabel", "")->setComment("label of the LHCInfoPerLS record");
0147   desc.add<std::string>("lhcInfoPerFillLabel", "")->setComment("label of the LHCInfoPerFill record");
0148   desc.add<bool>("useNewLHCInfo", false)->setComment("flag whether to use new LHCInfoPer* records or old LHCInfo");
0149
0150   desc.add<unsigned int>("rpId_45_N", 0)->setComment("decimal RP id for 45 near");
0151   desc.add<unsigned int>("rpId_45_F", 0)->setComment("decimal RP id for 45 far");
0152   desc.add<unsigned int>("rpId_56_N", 0)->setComment("decimal RP id for 56 near");
0153   desc.add<unsigned int>("rpId_56_F", 0)->setComment("decimal RP id for 56 far");
0154
0155   desc.add<std::string>("outputFile", "output.root")->setComment("output file name");
0156
0157   desc.addUntracked<unsigned int>("verbosity", 0)->setComment("verbosity level");
0158
0159   descriptions.add("ctppsProtonReconstructionEfficiencyEstimatorMCDefault", desc);
0160 }
0161
0162 //----------------------------------------------------------------------------------------------------
0163
0164 void CTPPSProtonReconstructionEfficiencyEstimatorMC::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0165   std::ostringstream os;
0166
0167   // get conditions
0168   const LHCInfoCombined lhcInfoCombined(
0169       iSetup, lhcInfoPerLSToken_, lhcInfoPerFillToken_, lhcInfoToken_, useNewLHCInfo_);
0170
0171   // get input
0172   edm::Handle<edm::HepMCProduct> hHepMCAfterSmearing;
0173   iEvent.getByToken(tokenHepMCAfterSmearing_, hHepMCAfterSmearing);
0174   HepMC::GenEvent *hepMCEventAfterSmearing = (HepMC::GenEvent *)hHepMCAfterSmearing->GetEvent();
0175
0176   edm::Handle<std::map<int, edm::DetSetVector<TotemRPRecHit>>> hStripRecHitsPerParticle;
0177   iEvent.getByToken(tokenStripRecHitsPerParticle_, hStripRecHitsPerParticle);
0178
0179   edm::Handle<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>> hPixelRecHitsPerParticle;
0180   iEvent.getByToken(tokenPixelRecHitsPerParticle_, hPixelRecHitsPerParticle);
0181
0182   edm::Handle<CTPPSLocalTrackLiteCollection> tracks;
0183   iEvent.getByToken(tracksToken_, tracks);
0184
0185   Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0186   iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0187
0188   // buffer for particle information
0189   struct ParticleInfo {
0190     unsigned int arm = 2;
0191     double xi = 0.;
0192     std::map<unsigned int, unsigned int> recHitsPerRP;
0193     bool inAcceptanceNear = false, inAcceptanceFar = false, inAcceptance = false;
0194   };
0195
0196   std::map<int, ParticleInfo> particleInfo;  // barcode --> info
0197
0198   // process HepMC
0199   for (auto it = hepMCEventAfterSmearing->particles_begin(); it != hepMCEventAfterSmearing->particles_end(); ++it) {
0200     const auto &part = *it;
0201
0202     // accept only stable non-beam protons
0203     if (part->pdg_id() != 2212)
0204       continue;
0205
0206     if (part->status() != 1)
0207       continue;
0208
0209     if (part->is_beam())
0210       continue;
0211
0212     const auto &mom = part->momentum();
0213
0214     if (mom.e() < 4500.)
0215       continue;
0216
0217     ParticleInfo info;
0218
0219     info.arm = (mom.z() > 0.) ? 0 : 1;
0220
0221     const double p_nom = lhcInfoCombined.energy;
0222     info.xi = (p_nom - mom.rho()) / p_nom;
0223
0224     particleInfo[part->barcode()] = std::move(info);
0225   }
0226
0227   // check acceptance
0228   for (const auto &pp : *hStripRecHitsPerParticle) {
0229     const auto barcode = pp.first;
0230
0231     for (const auto &ds : pp.second) {
0232       CTPPSDetId detId(ds.detId());
0233       CTPPSDetId rpId = detId.rpId();
0234       particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
0235     }
0236   }
0237
0238   for (const auto &pp : *hPixelRecHitsPerParticle) {
0239     const auto barcode = pp.first;
0240
0241     for (const auto &ds : pp.second) {
0242       CTPPSDetId detId(ds.detId());
0243       CTPPSDetId rpId = detId.rpId();
0244       particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
0245     }
0246   }
0247
0248   std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
0249
0250   for (auto &pp : particleInfo) {
0251     if (verbosity_)
0252       os << "* barcode=" << pp.first << ", arm=" << pp.second.arm << ", xi=" << pp.second.xi << std::endl;
0253
0254     for (const auto &rpp : pp.second.recHitsPerRP) {
0255       CTPPSDetId rpId(rpp.first);
0256       unsigned int needed_rec_hits = 1000;
0257       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip)
0258         needed_rec_hits = 6;
0259       if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel)
0260         needed_rec_hits = 3;
0261
0262       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0263
0264       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0265         if (rpDecId == rpDecId_near_[rpId.arm()])
0266           isStripRPNear[rpId.arm()] = true;
0267         if (rpDecId == rpDecId_far_[rpId.arm()])
0268           isStripRPFar[rpId.arm()] = true;
0269       }
0270
0271       if (rpp.second >= needed_rec_hits) {
0272         if (rpDecId == rpDecId_near_[rpId.arm()])
0273           pp.second.inAcceptanceNear = true;
0274         if (rpDecId == rpDecId_far_[rpId.arm()])
0275           pp.second.inAcceptanceFar = true;
0276       }
0277
0278       if (verbosity_)
0279         os << "    RP " << rpDecId << ": " << rpp.second << " hits" << std::endl;
0280     }
0281
0282     pp.second.inAcceptance = pp.second.inAcceptanceNear && pp.second.inAcceptanceFar;
0283
0284     if (verbosity_)
0285       os << "    inAcceptance: near=" << pp.second.inAcceptanceNear << ", far=" << pp.second.inAcceptanceFar
0286          << ", global=" << pp.second.inAcceptance << std::endl;
0287   }
0288
0289   // count particles in acceptance
0290   struct ArmCounter {
0291     unsigned int near = 0, far = 0, global = 0;
0292   };
0293   std::map<unsigned int, ArmCounter> nParticlesInAcceptance;
0294   for (auto &pp : particleInfo) {
0295     auto &counter = nParticlesInAcceptance[pp.second.arm];
0296     if (pp.second.inAcceptanceNear)
0297       counter.near++;
0298     if (pp.second.inAcceptanceFar)
0299       counter.far++;
0300     if (pp.second.inAcceptance)
0301       counter.global++;
0302   }
0303
0304   // count reconstructed tracks
0305   std::map<unsigned int, ArmCounter> nReconstructedTracks;
0306   for (const auto &tr : *tracks) {
0307     CTPPSDetId rpId(tr.rpId());
0308     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0309
0310     if (rpDecId == rpDecId_near_[rpId.arm()])
0311       nReconstructedTracks[rpId.arm()].near++;
0312     if (rpDecId == rpDecId_far_[rpId.arm()])
0313       nReconstructedTracks[rpId.arm()].far++;
0314   }
0315
0316   // count reconstructed protons
0317   std::map<unsigned int, unsigned int> nReconstructedProtons;
0318   for (const auto &pr : *hRecoProtonsMultiRP) {
0319     if (!pr.validFit())
0320       continue;
0321
0322     unsigned int arm = 2;
0323     if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector45)
0324       arm = 0;
0325     if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector56)
0326       arm = 1;
0327
0328     nReconstructedProtons[arm]++;
0329   }
0330
0331   // fill plots
0332   for (unsigned int arm = 0; arm < 2; arm++) {
0333     const auto &npa = nParticlesInAcceptance[arm];
0334     const auto &nrt = nReconstructedTracks[arm];
0335
0336     if (verbosity_)
0337       os << "* arm " << arm << ": nRecoProtons=" << nReconstructedProtons[arm]
0338          << " (tracks near=" << nReconstructedTracks[arm].near << ", far=" << nReconstructedTracks[arm].far
0339          << "), nAcc=" << npa.global << " (near=" << npa.near << ", far=" << npa.far << ")" << std::endl;
0340
0341     // skip event if no track in global acceptance
0342     if (npa.global < 1)
0343       continue;
0344
0345     const auto &p = plots_[arm][npa.global];
0346
0347     p.h_n_part_acc_nr->Fill(npa.near);
0348     p.h_n_part_acc_fr->Fill(npa.far);
0349
0350     // skip events with some local reconstruction inefficiency
0351     if (nrt.near != npa.near || nrt.far != npa.far)
0352       continue;
0353
0354     const double eff = double(nReconstructedProtons[arm]) / npa.global;
0355
0356     if (verbosity_)
0357       os << "    eff=" << eff << std::endl;
0358
0359     for (auto &pp : particleInfo) {
0360       if (pp.second.arm != arm || !pp.second.inAcceptance)
0361         continue;
0362
0363       p.p_eff_vs_xi->Fill(pp.second.xi, eff);
0364     }
0365   }
0366
0367   if (verbosity_)
0368     edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
0369 }
0370
0371 //----------------------------------------------------------------------------------------------------
0372
0373 void CTPPSProtonReconstructionEfficiencyEstimatorMC::endJob() {
0374   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0375
0376   for (const auto &ait : plots_) {
0377     char buf[100];
0378     sprintf(buf, "arm%u", ait.first);
0379     TDirectory *d_arm = f_out->mkdir(buf);
0380
0381     for (const auto &npit : ait.second) {
0382       sprintf(buf, "%u", npit.first);
0383       TDirectory *d_np = d_arm->mkdir(buf);
0384       gDirectory = d_np;
0385
0386       npit.second.write();
0387     }
0388   }
0389 }
0390
0391 //----------------------------------------------------------------------------------------------------
0392
0393 DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorMC);