290 std::ostringstream os;
298 data_[0].UpdateOptics(*hOpticalFunctions);
299 data_[1].UpdateOptics(*hOpticalFunctions);
310 std::vector<unsigned int> sel_track_idc;
311 for (
unsigned int idx = 0;
idx < hTracks->size(); ++
idx) {
312 const auto &tr = hTracks->at(
idx);
324 sel_track_idc.push_back(
idx);
329 os <<
"* tracks (pre-selected):" << std::endl;
331 for (
const auto idx : sel_track_idc) {
332 const auto &tr = hTracks->at(
idx);
336 os <<
" [" <<
idx <<
"] RP=" << decRPId <<
", x=" << tr.getX() <<
", y=" << tr.getY() << std::endl;
339 os <<
"* protons:" << std::endl;
340 for (
unsigned int i = 0;
i < hRecoProtonsMultiRP->size(); ++
i) {
341 const auto &
pr = (*hRecoProtonsMultiRP)[
i];
342 os <<
" [" <<
i <<
"] track indices: ";
343 for (
const auto &trr :
pr.contributingLocalTracks())
344 os << trr.key() <<
", ";
350 map<unsigned int, bool> hasTracksInArm;
352 for (
const auto idx_i : sel_track_idc) {
353 const auto &tr_i = hTracks->at(idx_i);
355 unsigned int decRPId_i = rpId_i.
arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
357 for (
const auto idx_j : sel_track_idc) {
358 const auto &tr_j = hTracks->at(idx_j);
360 unsigned int decRPId_j = rpId_j.
arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
363 unsigned int arm = 123;
364 for (
const auto &ap :
data_) {
365 if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
373 hasTracksInArm[
arm] =
true;
376 auto &ad = data_[
arm];
377 const double de_x = tr_j.getX() - tr_i.getX();
378 const double de_y = tr_j.getY() - tr_i.getY();
381 ad.h_de_x->Fill(de_x);
382 ad.h_de_y->Fill(de_y);
385 ad.h2_de_y_vs_de_x->Fill(de_x, de_y);
390 for (
auto &ap : data_) {
391 if (hasTracksInArm[ap.first])
392 ap.second.n_events_with_tracks++;
396 for (
auto &ap : data_) {
397 auto &ad = ap.second;
400 if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.)
403 std::vector<std::pair<unsigned int, TH1D *>>
m;
404 m.emplace_back(0, ad.h_de_x.get());
405 m.emplace_back(1, ad.h_de_y.get());
407 for (
const auto &
p : m) {
408 double max_pos = -1E100, max_val = -1.;
409 for (
int bi = 1; bi <
p.second->GetNbinsX(); ++bi) {
410 const double pos =
p.second->GetBinCenter(bi);
411 const double val =
p.second->GetBinContent(bi);
419 const double sig = 0.2;
421 ff_->SetParameters(max_val, max_pos, sig, 0.);
422 p.second->Fit(
ff_.get(),
"Q",
"", max_pos - 3. * sig, max_pos + 3. * sig);
423 p.second->Fit(
ff_.get(),
"Q",
"", max_pos - 3. * sig, max_pos + 3. * sig);
426 ad.de_x_mean =
ff_->GetParameter(1);
427 ad.de_x_sigma = fabs(
ff_->GetParameter(2));
430 ad.de_y_mean =
ff_->GetParameter(1);
431 ad.de_y_sigma = fabs(
ff_->GetParameter(2));
436 os <<
"* fitting arm " << ap.first << std::endl
437 <<
" de_x: mean = " << ad.de_x_mean <<
", sigma = " << ad.de_x_sigma << std::endl
438 <<
" de_y: mean = " << ad.de_y_mean <<
", sigma = " << ad.de_y_sigma;
444 struct ArmEventData {
445 std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
447 std::set<unsigned int> reco_proton_idc;
449 std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
452 std::map<unsigned int, ArmEventData> armEventData;
455 for (
const auto idx_i : sel_track_idc) {
456 const auto &tr_i = hTracks->at(idx_i);
458 unsigned int decRPId_i = rpId_i.
arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
460 for (
const auto idx_j : sel_track_idc) {
461 const auto &tr_j = hTracks->at(idx_j);
463 unsigned int decRPId_j = rpId_j.
arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
466 unsigned int arm = 123;
467 for (
const auto &ap : data_) {
468 if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
476 auto &ad = data_[
arm];
477 const double de_x = tr_j.getX() - tr_i.getX();
478 const double de_y = tr_j.getY() - tr_i.getY();
480 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
481 const double n_si = ad.n_sigmas[nsi];
482 const bool match_x = fabs(de_x - ad.de_x_mean) < n_si * ad.de_x_sigma;
483 const bool match_y = fabs(de_y - ad.de_y_mean) < n_si * ad.de_y_sigma;
484 if (match_x && match_y)
485 armEventData[
arm].matched_track_idc[nsi].insert(idx_i);
491 for (
unsigned int i = 0;
i < hRecoProtonsMultiRP->size(); ++
i) {
492 const auto &proton = (*hRecoProtonsMultiRP)[
i];
494 CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->getRPId());
495 unsigned int arm =
rpId.arm();
497 if (proton.validFit())
498 armEventData[arm].reco_proton_idc.insert(
i);
503 os <<
"* cmp matched tracks vs. reco protons" << std::endl;
505 for (
auto &ap : armEventData) {
506 auto &ad = data_[ap.first];
509 os <<
" arm " << ap.first << std::endl;
511 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
513 os <<
" nsi = " << nsi << std::endl;
515 for (
const auto &tri : ap.second.matched_track_idc[nsi]) {
516 const auto &
track = hTracks->at(tri);
518 bool some_proton_matching =
false;
521 os <<
" tri = " << tri << std::endl;
523 for (
const auto &pri : ap.second.reco_proton_idc) {
524 const auto &proton = (*hRecoProtonsMultiRP)[pri];
526 bool proton_matching =
false;
528 for (
const auto &pr_tr : proton.contributingLocalTracks()) {
530 unsigned int decRPId =
rpId.arm() * 100 +
rpId.station() * 10 +
rpId.rp();
532 if (decRPId != ad.rpId_N)
535 const double x = pr_tr->getX();
536 const double y = pr_tr->getY();
537 const double th = 1E-3;
539 const bool match = (fabs(x -
track.getX()) < th) && (fabs(y -
track.getY()) < th);
542 os <<
" pri = " << pri <<
": x_tr = " <<
track.getX() <<
", x_pr = " << x
543 <<
", match = " << match << std::endl;
546 proton_matching =
true;
551 if (proton_matching) {
552 some_proton_matching =
true;
558 os <<
" --> some_proton_matching " << some_proton_matching << std::endl;
560 if (some_proton_matching)
561 ap.second.matched_track_with_prot_idc[nsi].insert(tri);
563 ap.second.matched_track_without_prot_idc[nsi].insert(tri);
570 for (
auto &ap : armEventData) {
571 auto &ad = data_[ap.first];
573 if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
576 os <<
"* results for arm " << ap.first << std::endl;
578 os <<
" reco_proton_idc: ";
579 for (
const auto &
idx : ap.second.reco_proton_idc)
583 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
584 os <<
" n_si = " << ad.n_sigmas[nsi] << std::endl;
586 os <<
" matched_track_idc: ";
587 for (
const auto &idx : ap.second.matched_track_idc[nsi])
591 os <<
" matched_track_with_prot_idc: ";
592 for (
const auto &idx : ap.second.matched_track_with_prot_idc[nsi])
596 os <<
" matched_track_without_prot_idc: ";
597 for (
const auto &idx : ap.second.matched_track_without_prot_idc[nsi])
605 for (
auto &ap : armEventData) {
606 auto &ad = data_[ap.first];
609 if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
613 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
614 const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size();
615 const unsigned int n_rec_prot = ap.second.reco_proton_idc.size();
618 if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
622 const double eff = double(n_rec_prot) / n_exp_prot;
624 for (
unsigned int tri : ap.second.matched_track_idc[nsi]) {
625 const double x_N = hTracks->at(tri).getX();
626 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
628 ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
629 ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
631 ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
632 ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
636 for (
const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) {
637 const double x_N = hTracks->at(tri).getX();
638 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
640 ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
641 ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
643 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
644 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
647 for (
const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) {
648 const double x_N = hTracks->at(tri).getX();
649 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
651 ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
652 ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
654 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
655 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
661 edm::LogInfo(
"CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str();
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::map< unsigned int, ArmData > data_
edm::EDGetTokenT< reco::ForwardProtonCollection > tokenRecoProtonsMultiRP_
std::unique_ptr< TF1 > ff_
bool pixelDiscardBXShiftedTracks_
bool check(const edm::EventSetup &iSetup)
Base class for CTPPS detector IDs.
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tokenTracks_
unsigned int n_prep_events_
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
edm::ESWatcher< CTPPSInterpolatedOpticsRcd > opticsWatcher_