289 std::ostringstream os;
308 std::vector<unsigned int> sel_track_idc;
309 for (
unsigned int idx = 0;
idx < hTracks->size(); ++
idx) {
310 const auto &tr = hTracks->at(
idx);
322 sel_track_idc.push_back(
idx);
327 os <<
"* tracks (pre-selected):" << std::endl;
329 for (
const auto idx : sel_track_idc) {
330 const auto &tr = hTracks->at(
idx);
334 os <<
" [" <<
idx <<
"] RP=" <<
decRPId <<
", x=" << tr.x() <<
", y=" << tr.y() << std::endl;
337 os <<
"* protons:" << std::endl;
338 for (
unsigned int i = 0;
i < hRecoProtonsMultiRP->size(); ++
i) {
339 const auto &
pr = (*hRecoProtonsMultiRP)[
i];
340 os <<
" [" <<
i <<
"] track indices: ";
341 for (
const auto &trr :
pr.contributingLocalTracks())
342 os << trr.key() <<
", ";
348 map<unsigned int, bool> hasTracksInArm;
350 for (
const auto idx_i : sel_track_idc) {
351 const auto &tr_i = hTracks->at(idx_i);
353 unsigned int decRPId_i = rpId_i.
arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
355 for (
const auto idx_j : sel_track_idc) {
356 const auto &tr_j = hTracks->at(idx_j);
358 unsigned int decRPId_j = rpId_j.
arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
361 unsigned int arm = 123;
362 for (
const auto &ap :
data_) {
363 if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
371 hasTracksInArm[
arm] =
true;
375 const double de_x = tr_j.x() - tr_i.x();
376 const double de_y = tr_j.y() - tr_i.y();
379 ad.h_de_x->Fill(
de_x);
380 ad.h_de_y->Fill(
de_y);
383 ad.h2_de_y_vs_de_x->Fill(
de_x,
de_y);
388 for (
auto &ap :
data_) {
389 if (hasTracksInArm[ap.first])
390 ap.second.n_events_with_tracks++;
394 for (
auto &ap :
data_) {
395 auto &ad = ap.second;
398 if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.)
401 std::vector<std::pair<unsigned int, TH1D *>>
m;
402 m.emplace_back(0, ad.h_de_x.get());
403 m.emplace_back(1, ad.h_de_y.get());
405 for (
const auto &
p :
m) {
406 double max_pos = -1E100, max_val = -1.;
407 for (
int bi = 1; bi <
p.second->GetNbinsX(); ++bi) {
408 const double pos =
p.second->GetBinCenter(bi);
409 const double val =
p.second->GetBinContent(bi);
417 const double sig = 0.2;
419 ff_->SetParameters(max_val, max_pos, sig, 0.);
420 p.second->Fit(
ff_.get(),
"Q",
"", max_pos - 3. * sig, max_pos + 3. * sig);
421 p.second->Fit(
ff_.get(),
"Q",
"", max_pos - 3. * sig, max_pos + 3. * sig);
424 ad.de_x_mean =
ff_->GetParameter(1);
425 ad.de_x_sigma = fabs(
ff_->GetParameter(2));
428 ad.de_y_mean =
ff_->GetParameter(1);
429 ad.de_y_sigma = fabs(
ff_->GetParameter(2));
434 os <<
"* fitting arm " << ap.first << std::endl
435 <<
" de_x: mean = " << ad.de_x_mean <<
", sigma = " << ad.de_x_sigma << std::endl
436 <<
" de_y: mean = " << ad.de_y_mean <<
", sigma = " << ad.de_y_sigma;
442 struct ArmEventData {
443 std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
445 std::set<unsigned int> reco_proton_idc;
447 std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
450 std::map<unsigned int, ArmEventData> armEventData;
453 for (
const auto idx_i : sel_track_idc) {
454 const auto &tr_i = hTracks->at(idx_i);
456 unsigned int decRPId_i = rpId_i.
arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
458 for (
const auto idx_j : sel_track_idc) {
459 const auto &tr_j = hTracks->at(idx_j);
461 unsigned int decRPId_j = rpId_j.
arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
464 unsigned int arm = 123;
465 for (
const auto &ap :
data_) {
466 if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
475 const double de_x = tr_j.x() - tr_i.x();
476 const double de_y = tr_j.y() - tr_i.y();
478 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
479 const double n_si = ad.n_sigmas[nsi];
480 const bool match_x = fabs(
de_x - ad.de_x_mean) < n_si * ad.de_x_sigma;
481 const bool match_y = fabs(
de_y - ad.de_y_mean) < n_si * ad.de_y_sigma;
482 if (match_x && match_y)
483 armEventData[
arm].matched_track_idc[nsi].insert(idx_i);
489 for (
unsigned int i = 0;
i < hRecoProtonsMultiRP->size(); ++
i) {
490 const auto &proton = (*hRecoProtonsMultiRP)[
i];
495 if (proton.validFit())
496 armEventData[
arm].reco_proton_idc.insert(
i);
501 os <<
"* cmp matched tracks vs. reco protons" << std::endl;
503 for (
auto &ap : armEventData) {
504 auto &ad =
data_[ap.first];
507 os <<
" arm " << ap.first << std::endl;
509 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
511 os <<
" nsi = " << nsi << std::endl;
513 for (
const auto &tri : ap.second.matched_track_idc[nsi]) {
514 const auto &
track = hTracks->at(tri);
516 bool some_proton_matching =
false;
519 os <<
" tri = " << tri << std::endl;
521 for (
const auto &pri : ap.second.reco_proton_idc) {
522 const auto &proton = (*hRecoProtonsMultiRP)[pri];
524 bool proton_matching =
false;
526 for (
const auto &pr_tr : proton.contributingLocalTracks()) {
533 const double x = pr_tr->x();
534 const double y = pr_tr->y();
535 const double th = 1E-3;
540 os <<
" pri = " << pri <<
": x_tr = " <<
track.x() <<
", x_pr = " <<
x
541 <<
", match = " <<
match << std::endl;
544 proton_matching =
true;
549 if (proton_matching) {
550 some_proton_matching =
true;
556 os <<
" --> some_proton_matching " << some_proton_matching << std::endl;
558 if (some_proton_matching)
559 ap.second.matched_track_with_prot_idc[nsi].insert(tri);
561 ap.second.matched_track_without_prot_idc[nsi].insert(tri);
568 for (
auto &ap : armEventData) {
569 auto &ad =
data_[ap.first];
571 if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
574 os <<
"* results for arm " << ap.first << std::endl;
576 os <<
" reco_proton_idc: ";
577 for (
const auto &
idx : ap.second.reco_proton_idc)
581 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
582 os <<
" n_si = " << ad.n_sigmas[nsi] << std::endl;
584 os <<
" matched_track_idc: ";
585 for (
const auto &
idx : ap.second.matched_track_idc[nsi])
589 os <<
" matched_track_with_prot_idc: ";
590 for (
const auto &
idx : ap.second.matched_track_with_prot_idc[nsi])
594 os <<
" matched_track_without_prot_idc: ";
595 for (
const auto &
idx : ap.second.matched_track_without_prot_idc[nsi])
603 for (
auto &ap : armEventData) {
604 auto &ad =
data_[ap.first];
607 if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
611 for (
unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
612 const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size();
613 const unsigned int n_rec_prot = ap.second.reco_proton_idc.size();
616 if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
620 const double eff = double(n_rec_prot) / n_exp_prot;
622 for (
unsigned int tri : ap.second.matched_track_idc[nsi]) {
623 const double x_N = hTracks->at(tri).x();
624 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
626 ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
627 ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
629 ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
630 ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
634 for (
const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) {
635 const double x_N = hTracks->at(tri).x();
636 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
638 ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
639 ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
641 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
642 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
645 for (
const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) {
646 const double x_N = hTracks->at(tri).x();
647 const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);
649 ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
650 ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
652 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
653 ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
659 edm::LogInfo(
"CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str();