Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:33

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);