165 std::ostringstream os;
189 struct ParticleInfo {
190 unsigned int arm = 2;
192 std::map<unsigned int, unsigned int> recHitsPerRP;
193 bool inAcceptanceNear =
false, inAcceptanceFar =
false, inAcceptance =
false;
196 std::map<int, ParticleInfo> particleInfo;
199 for (
auto it = hepMCEventAfterSmearing->particles_begin();
it != hepMCEventAfterSmearing->particles_end(); ++
it) {
203 if (
part->pdg_id() != 2212)
206 if (
part->status() != 1)
212 const auto &mom =
part->momentum();
219 info.arm = (mom.z() > 0.) ? 0 : 1;
221 const double p_nom = lhcInfoCombined.energy;
222 info.xi = (p_nom - mom.rho()) / p_nom;
228 for (
const auto &
pp : *hStripRecHitsPerParticle) {
229 const auto barcode =
pp.first;
231 for (
const auto &ds :
pp.second) {
234 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
238 for (
const auto &
pp : *hPixelRecHitsPerParticle) {
239 const auto barcode =
pp.first;
241 for (
const auto &ds :
pp.second) {
244 particleInfo[barcode].recHitsPerRP[rpId] += ds.size();
248 std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
250 for (
auto &
pp : particleInfo) {
252 os <<
"* barcode=" <<
pp.first <<
", arm=" <<
pp.second.arm <<
", xi=" <<
pp.second.xi << std::endl;
254 for (
const auto &rpp :
pp.second.recHitsPerRP) {
256 unsigned int needed_rec_hits = 1000;
262 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
266 isStripRPNear[rpId.arm()] =
true;
268 isStripRPFar[rpId.arm()] =
true;
271 if (rpp.second >= needed_rec_hits) {
273 pp.second.inAcceptanceNear =
true;
275 pp.second.inAcceptanceFar =
true;
279 os <<
" RP " << rpDecId <<
": " << rpp.second <<
" hits" << std::endl;
282 pp.second.inAcceptance =
pp.second.inAcceptanceNear &&
pp.second.inAcceptanceFar;
285 os <<
" inAcceptance: near=" <<
pp.second.inAcceptanceNear <<
", far=" <<
pp.second.inAcceptanceFar
286 <<
", global=" <<
pp.second.inAcceptance << std::endl;
291 unsigned int near = 0, far = 0, global = 0;
293 std::map<unsigned int, ArmCounter> nParticlesInAcceptance;
294 for (
auto &
pp : particleInfo) {
295 auto &
counter = nParticlesInAcceptance[
pp.second.arm];
296 if (
pp.second.inAcceptanceNear)
298 if (
pp.second.inAcceptanceFar)
300 if (
pp.second.inAcceptance)
305 std::map<unsigned int, ArmCounter> nReconstructedTracks;
306 for (
const auto &tr : *
tracks) {
308 unsigned int rpDecId = rpId.
arm() * 100 + rpId.station() * 10 + rpId.rp();
311 nReconstructedTracks[rpId.arm()].near++;
313 nReconstructedTracks[rpId.arm()].far++;
317 std::map<unsigned int, unsigned int> nReconstructedProtons;
318 for (
const auto &
pr : *hRecoProtonsMultiRP) {
322 unsigned int arm = 2;
328 nReconstructedProtons[
arm]++;
332 for (
unsigned int arm = 0;
arm < 2;
arm++) {
333 const auto &npa = nParticlesInAcceptance[
arm];
334 const auto &nrt = nReconstructedTracks[
arm];
337 os <<
"* arm " <<
arm <<
": nRecoProtons=" << nReconstructedProtons[
arm]
338 <<
" (tracks near=" << nReconstructedTracks[
arm].near <<
", far=" << nReconstructedTracks[
arm].far
339 <<
"), nAcc=" << npa.global <<
" (near=" << npa.near <<
", far=" << npa.far <<
")" << std::endl;
347 p.h_n_part_acc_nr->Fill(npa.near);
348 p.h_n_part_acc_fr->Fill(npa.far);
351 if (nrt.near != npa.near || nrt.far != npa.far)
354 const double eff = double(nReconstructedProtons[
arm]) / npa.global;
357 os <<
" eff=" << eff << std::endl;
359 for (
auto &
pp : particleInfo) {
360 if (
pp.second.arm !=
arm || !
pp.second.inAcceptance)
363 p.p_eff_vs_xi->Fill(
pp.second.xi, eff);
368 edm::LogInfo(
"CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tracksToken_
const bool useNewLHCInfo_
std::map< unsigned int, std::map< unsigned int, PlotGroup > > plots_
std::map< unsigned int, unsigned int > rpDecId_near_
const edm::ESGetToken< LHCInfoPerLS, LHCInfoPerLSRcd > lhcInfoPerLSToken_
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.
const edm::ESGetToken< LHCInfo, LHCInfoRcd > lhcInfoToken_
std::map< unsigned int, unsigned int > rpDecId_far_
const edm::ESGetToken< LHCInfoPerFill, LHCInfoPerFillRcd > lhcInfoPerFillToken_
edm::EDGetTokenT< std::map< int, edm::DetSetVector< TotemRPRecHit > > > tokenStripRecHitsPerParticle_
if(threadIdxLocalY==0 &&threadIdxLocalX==0)