132 std::ostringstream os;
155 struct ParticleInfo {
156 unsigned int arm = 2;
158 std::map<unsigned int, unsigned int> recHitsPerRP;
159 bool inAcceptanceNear =
false, inAcceptanceFar =
false, inAcceptance =
false;
162 std::map<int, ParticleInfo> particleInfo;
165 for (
auto it = hepMCEventAfterSmearing->particles_begin(); it != hepMCEventAfterSmearing->particles_end(); ++it) {
166 const auto &
part = *it;
169 if (
part->pdg_id() != 2212)
172 if (
part->status() != 1)
178 const auto &mom =
part->momentum();
185 info.arm = (mom.z() > 0.) ? 0 : 1;
187 const double p_nom = lhcInfo.energy();
188 info.xi = (p_nom - mom.rho()) / p_nom;
194 for (
const auto &
pp : *hStripRecHitsPerParticle) {
195 const auto barcode =
pp.first;
197 for (
const auto &ds :
pp.second) {
200 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
204 for (
const auto &
pp : *hPixelRecHitsPerParticle) {
205 const auto barcode =
pp.first;
207 for (
const auto &ds :
pp.second) {
210 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
214 std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
216 for (
auto &
pp : particleInfo) {
218 os <<
"* barcode=" <<
pp.first <<
", arm=" <<
pp.second.arm <<
", xi=" <<
pp.second.xi << std::endl;
220 for (
const auto &rpp :
pp.second.recHitsPerRP) {
222 unsigned int needed_rec_hits = 1000;
228 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
232 isStripRPNear[rpId.arm()] =
true;
234 isStripRPFar[rpId.arm()] =
true;
237 if (rpp.second >= needed_rec_hits) {
239 pp.second.inAcceptanceNear =
true;
241 pp.second.inAcceptanceFar =
true;
245 os <<
" RP " << rpDecId <<
": " << rpp.second <<
" hits" << std::endl;
248 pp.second.inAcceptance =
pp.second.inAcceptanceNear &&
pp.second.inAcceptanceFar;
251 os <<
" inAcceptance: near=" <<
pp.second.inAcceptanceNear <<
", far=" <<
pp.second.inAcceptanceFar
252 <<
", global=" <<
pp.second.inAcceptance << std::endl;
257 unsigned int near = 0, far = 0, global = 0;
259 std::map<unsigned int, ArmCounter> nParticlesInAcceptance;
260 for (
auto &
pp : particleInfo) {
261 auto &
counter = nParticlesInAcceptance[
pp.second.arm];
262 if (
pp.second.inAcceptanceNear)
264 if (
pp.second.inAcceptanceFar)
266 if (
pp.second.inAcceptance)
271 std::map<unsigned int, ArmCounter> nReconstructedTracks;
272 for (
const auto &tr : *
tracks) {
274 unsigned int rpDecId = rpId.
arm() * 100 + rpId.station() * 10 + rpId.rp();
277 nReconstructedTracks[rpId.arm()].near++;
279 nReconstructedTracks[rpId.arm()].far++;
283 std::map<unsigned int, unsigned int> nReconstructedProtons;
284 for (
const auto &
pr : *hRecoProtonsMultiRP) {
288 unsigned int arm = 2;
294 nReconstructedProtons[
arm]++;
298 for (
unsigned int arm = 0;
arm < 2;
arm++) {
299 const auto &npa = nParticlesInAcceptance[
arm];
300 const auto &nrt = nReconstructedTracks[
arm];
303 os <<
"* arm " <<
arm <<
": nRecoProtons=" << nReconstructedProtons[
arm]
304 <<
" (tracks near=" << nReconstructedTracks[
arm].near <<
", far=" << nReconstructedTracks[
arm].far
305 <<
"), nAcc=" << npa.global <<
" (near=" << npa.near <<
", far=" << npa.far <<
")" << std::endl;
313 p.h_n_part_acc_nr->Fill(npa.near);
314 p.h_n_part_acc_fr->Fill(npa.far);
317 if (nrt.near != npa.near || nrt.far != npa.far)
320 const double eff = double(nReconstructedProtons[
arm]) / npa.global;
323 os <<
" eff=" << eff << std::endl;
325 for (
auto &
pp : particleInfo) {
326 if (
pp.second.arm !=
arm || !
pp.second.inAcceptance)
329 p.p_eff_vs_xi->Fill(
pp.second.xi, eff);
334 edm::LogInfo(
"CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tracksToken_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
std::map< unsigned int, std::map< unsigned int, PlotGroup > > plots_
std::map< unsigned int, unsigned int > rpDecId_near_
edm::EDGetTokenT< reco::ForwardProtonCollection > tokenRecoProtonsMultiRP_
edm::EDGetTokenT< std::map< int, edm::DetSetVector< CTPPSPixelRecHit > > > tokenPixelRecHitsPerParticle_
edm::EDGetTokenT< edm::HepMCProduct > tokenHepMCAfterSmearing_
const HepMC::GenEvent * GetEvent() const
Log< level::Info, false > LogInfo
Base class for CTPPS detector IDs.
std::map< unsigned int, unsigned int > rpDecId_far_
edm::EDGetTokenT< std::map< int, edm::DetSetVector< TotemRPRecHit > > > tokenStripRecHitsPerParticle_
edm::ESGetToken< LHCInfo, LHCInfoRcd > lhcInfoESToken_