Inheritance diagram for CTPPSProtonReconstructionEfficiencyEstimatorMC:

Classes
struct	PlotGroup

Public Member Functions
	CTPPSProtonReconstructionEfficiencyEstimatorMC (const edm::ParameterSet &)

Public Member Functions inherited from edm::one::EDAnalyzer<>
	EDAnalyzer ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDAnalyzerBase () override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &) override

void	endJob () override

Private Attributes
std::string	lhcInfoLabel_

std::string	outputFile_

std::map< unsigned int, std::map< unsigned int, PlotGroup > >	plots_

std::map< unsigned int, unsigned int >	rpDecId_far_

std::map< unsigned int, unsigned int >	rpDecId_near_

unsigned int	rpId_45_F_

unsigned int	rpId_45_N_

unsigned int	rpId_56_F_

unsigned int	rpId_56_N_

edm::EDGetTokenT< edm::HepMCProduct >	tokenHepMCAfterSmearing_

edm::EDGetTokenT< std::map< int, edm::DetSetVector< CTPPSPixelRecHit > > >	tokenPixelRecHitsPerParticle_

edm::EDGetTokenT< reco::ForwardProtonCollection >	tokenRecoProtonsMultiRP_

edm::EDGetTokenT< std::map< int, edm::DetSetVector< TotemRPRecHit > > >	tokenStripRecHitsPerParticle_

edm::EDGetTokenT< CTPPSLocalTrackLiteCollection >	tracksToken_

unsigned int	verbosity_

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 39 of file CTPPSProtonReconstructionEfficiencyEstimatorMC.cc.

Constructor & Destructor Documentation

CTPPSProtonReconstructionEfficiencyEstimatorMC::CTPPSProtonReconstructionEfficiencyEstimatorMC ( const edm::ParameterSet & iConfig )

explicit

Definition at line 95 of file CTPPSProtonReconstructionEfficiencyEstimatorMC.cc.

References protons_cff::arm, np, plots_, rpDecId_far_, rpDecId_near_, rpId_45_F_, rpId_45_N_, rpId_56_F_, and rpId_56_N_.

     : tokenHepMCAfterSmearing_(
           consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCAfterSmearing"))),
       tokenStripRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<TotemRPRecHit>>>(
           iConfig.getParameter<edm::InputTag>("tagStripRecHitsPerParticle"))),
       tokenPixelRecHitsPerParticle_(consumes<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>>(
           iConfig.getParameter<edm::InputTag>("tagPixelRecHitsPerParticle"))),
       tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
       tokenRecoProtonsMultiRP_(
           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
       lhcInfoLabel_(iConfig.getParameter<std::string>("lhcInfoLabel")),
 
       rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
       rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
 
       outputFile_(iConfig.getParameter<string>("outputFile")),
 
       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0))
 {
   rpDecId_near_[0] = rpId_45_N_;
   rpDecId_far_[0] = rpId_45_F_;
 
   rpDecId_near_[1] = rpId_56_N_;
   rpDecId_far_[1] = rpId_56_F_;
 
   // book plots
   for (unsigned int arm = 0; arm < 2; ++arm) {
     for (unsigned int np = 1; np <= 5; ++np)
       plots_[arm][np] = PlotGroup();
   }
 }

Member Function Documentation

void CTPPSProtonReconstructionEfficiencyEstimatorMC::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivate

Definition at line 132 of file CTPPSProtonReconstructionEfficiencyEstimatorMC.cc.

References protons_cff::arm, CTPPSDetId::arm(), LHCInfo::energy(), edm::EventSetup::get(), edm::Event::getByToken(), edm::HepMCProduct::GetEvent(), info(), lhcInfoLabel_, eostools::move(), AlCaHLTBitMon_ParallelJobs::p, plots_, createTree::pp, CTPPSDetId::rp(), rpDecId_far_, rpDecId_near_, year_2016_postTS2_cff::rpId, CTPPSDetId::sdTrackingPixel, CTPPSDetId::sdTrackingStrip, reco::ForwardProton::sector45, reco::ForwardProton::sector56, CTPPSDetId::station(), DetId::subdetId(), tokenHepMCAfterSmearing_, tokenPixelRecHitsPerParticle_, tokenRecoProtonsMultiRP_, tokenStripRecHitsPerParticle_, l1t::tracks, tracksToken_, verbosity_, and protons_cff::xi.

                                                                                                                 {
   std::ostringstream os;
 
   // get conditions
   edm::ESHandle<LHCInfo> hLHCInfo;
   iSetup.get<LHCInfoRcd>().get(lhcInfoLabel_, hLHCInfo);
 
   // get input
   edm::Handle<edm::HepMCProduct> hHepMCAfterSmearing;
   iEvent.getByToken(tokenHepMCAfterSmearing_, hHepMCAfterSmearing);
   HepMC::GenEvent *hepMCEventAfterSmearing = (HepMC::GenEvent *)hHepMCAfterSmearing->GetEvent();
 
   edm::Handle<std::map<int, edm::DetSetVector<TotemRPRecHit>>> hStripRecHitsPerParticle;
   iEvent.getByToken(tokenStripRecHitsPerParticle_, hStripRecHitsPerParticle);
 
   edm::Handle<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>> hPixelRecHitsPerParticle;
   iEvent.getByToken(tokenPixelRecHitsPerParticle_, hPixelRecHitsPerParticle);
 
   edm::Handle<CTPPSLocalTrackLiteCollection> tracks;
   iEvent.getByToken(tracksToken_, tracks);
 
   Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
   iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
 
   // buffer for particle information
   struct ParticleInfo {
     unsigned int arm = 2;
     double xi = 0.;
     std::map<unsigned int, unsigned int> recHitsPerRP;
     bool inAcceptanceNear = false, inAcceptanceFar = false, inAcceptance = false;
   };
 
   std::map<int, ParticleInfo> particleInfo;  // barcode --> info
 
   // process HepMC
   for (auto it = hepMCEventAfterSmearing->particles_begin(); it != hepMCEventAfterSmearing->particles_end(); ++it) {
     const auto &part = *it;
 
     // accept only stable non-beam protons
     if (part->pdg_id() != 2212)
       continue;
 
     if (part->status() != 1)
       continue;
 
     if (part->is_beam())
       continue;
 
     const auto &mom = part->momentum();
 
     if (mom.e() < 4500.)
       continue;
 
     ParticleInfo info;
 
     info.arm = (mom.z() > 0.) ? 0 : 1;
 
     const double p_nom = hLHCInfo->energy();
     info.xi = (p_nom - mom.rho()) / p_nom;
 
     particleInfo[part->barcode()] = std::move(info);
   }
 
   // check acceptance
   for (const auto &pp : *hStripRecHitsPerParticle) {
     const auto barcode = pp.first;
 
     for (const auto &ds : pp.second) {
       CTPPSDetId detId(ds.detId());
       CTPPSDetId rpId = detId.getRPId();
       particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
     }
   }
 
   for (const auto &pp : *hPixelRecHitsPerParticle) {
     const auto barcode = pp.first;
 
     for (const auto &ds : pp.second) {
       CTPPSDetId detId(ds.detId());
       CTPPSDetId rpId = detId.getRPId();
       particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
     }
   }
 
   std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
 
   for (auto &pp : particleInfo) {
     if (verbosity_)
       os << "* barcode=" << pp.first << ", arm=" << pp.second.arm << ", xi=" << pp.second.xi << std::endl;
 
     for (const auto &rpp : pp.second.recHitsPerRP) {
       CTPPSDetId rpId(rpp.first);
       unsigned int needed_rec_hits = 1000;
       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip)
         needed_rec_hits = 6;
       if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel)
         needed_rec_hits = 3;
 
       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) {
         if (rpDecId == rpDecId_near_[rpId.arm()])
           isStripRPNear[rpId.arm()] = true;
         if (rpDecId == rpDecId_far_[rpId.arm()])
           isStripRPFar[rpId.arm()] = true;
       }
 
       if (rpp.second >= needed_rec_hits) {
         if (rpDecId == rpDecId_near_[rpId.arm()])
           pp.second.inAcceptanceNear = true;
         if (rpDecId == rpDecId_far_[rpId.arm()])
           pp.second.inAcceptanceFar = true;
       }
 
       if (verbosity_)
         os << "    RP " << rpDecId << ": " << rpp.second << " hits" << std::endl;
     }
 
     pp.second.inAcceptance = pp.second.inAcceptanceNear && pp.second.inAcceptanceFar;
 
     if (verbosity_)
       os << "    inAcceptance: near=" << pp.second.inAcceptanceNear << ", far=" << pp.second.inAcceptanceFar << ", global=" << pp.second.inAcceptance << std::endl;
   }
 
   // count particles in acceptance
   struct ArmCounter {
     unsigned int near = 0, far = 0, global = 0;
   };
   std::map<unsigned int, ArmCounter> nParticlesInAcceptance;
   for (auto &pp : particleInfo) {
     auto &counter = nParticlesInAcceptance[pp.second.arm];
     if (pp.second.inAcceptanceNear)
       counter.near++;
     if (pp.second.inAcceptanceFar)
       counter.far++;
     if (pp.second.inAcceptance)
       counter.global++;
   }
 
   // count reconstructed tracks
   std::map<unsigned int, ArmCounter> nReconstructedTracks;
   for (const auto &tr : *tracks) {
     CTPPSDetId rpId(tr.getRPId());
     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
     if (rpDecId == rpDecId_near_[rpId.arm()])
       nReconstructedTracks[rpId.arm()].near++;
     if (rpDecId == rpDecId_far_[rpId.arm()])
       nReconstructedTracks[rpId.arm()].far++;
   }
 
   // count reconstructed protons
   std::map<unsigned int, unsigned int> nReconstructedProtons;
   for (const auto &pr : *hRecoProtonsMultiRP) {
     if (!pr.validFit())
       continue;
 
     unsigned int arm = 2;
     if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector45)
       arm = 0;
     if (pr.lhcSector() == reco::ForwardProton::LHCSector::sector56)
       arm = 1;
 
     nReconstructedProtons[arm]++;
   }
 
   // fill plots
   for (unsigned int arm = 0; arm < 2; arm++) {
     const auto &npa = nParticlesInAcceptance[arm];
     const auto &nrt = nReconstructedTracks[arm];
 
     if (verbosity_)
       os << "* arm " << arm << ": nRecoProtons="
         << nReconstructedProtons[arm] << " (tracks near=" << nReconstructedTracks[arm].near << ", far=" << nReconstructedTracks[arm].far
         << "), nAcc=" << npa.global << " (near=" << npa.near << ", far=" << npa.far << ")" << std::endl;
 
     // skip event if no track in global acceptance
     if (npa.global < 1)
       continue;
 
     const auto &p = plots_[arm][npa.global];
 
     p.h_n_part_acc_nr->Fill(npa.near);
     p.h_n_part_acc_fr->Fill(npa.far);
 
     // skip events with some local reconstruction inefficiency
     if (nrt.near != npa.near || nrt.far != npa.far)
       continue;
 
     const double eff = double(nReconstructedProtons[arm]) / npa.global;
 
     if (verbosity_)
       os << "    eff=" << eff << std::endl;
 
     for (auto &pp : particleInfo) {
       if (pp.second.arm != arm || !pp.second.inAcceptance)
         continue;
 
       p.p_eff_vs_xi->Fill(pp.second.xi, eff);
     }
   }
 
   if (verbosity_)
     edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
 }

void CTPPSProtonReconstructionEfficiencyEstimatorMC::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 340 of file CTPPSProtonReconstructionEfficiencyEstimatorMC.cc.

References DEFINE_FWK_MODULE, outputFile_, and plots_.

                                                             {
   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
 
   for (const auto &ait : plots_) {
     char buf[100];
     sprintf(buf, "arm%u", ait.first);
     TDirectory *d_arm = f_out->mkdir(buf);
 
     for (const auto &npit : ait.second) {
       sprintf(buf, "%u", npit.first);
       TDirectory *d_np = d_arm->mkdir(buf);
       gDirectory = d_np;
 
       npit.second.write();
     }
   }
 }