diff --git a/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc b/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc index 953f0caa26ba9..7f4f178bccd39 100644 --- a/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc +++ b/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc @@ -44,6 +44,8 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT tracksToken_; + bool pixelDiscardBXShiftedTracks_; + std::string lhcInfoLabel_; std::string opticsLabel_; @@ -59,25 +61,32 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> { struct AssociationCuts { bool x_cut_apply; - double x_cut_value; + double x_cut_mean, x_cut_value; bool y_cut_apply; - double y_cut_value; + double y_cut_mean, y_cut_value; bool xi_cut_apply; - double xi_cut_value; + double xi_cut_mean, xi_cut_value; bool th_y_cut_apply; - double th_y_cut_value; + double th_y_cut_mean, th_y_cut_value; double ti_tr_min; double ti_tr_max; void load(const edm::ParameterSet &ps) { x_cut_apply = ps.getParameter("x_cut_apply"); + x_cut_mean = ps.getParameter("x_cut_mean"); x_cut_value = ps.getParameter("x_cut_value"); + y_cut_apply = ps.getParameter("y_cut_apply"); + y_cut_mean = ps.getParameter("y_cut_mean"); y_cut_value = ps.getParameter("y_cut_value"); + xi_cut_apply = ps.getParameter("xi_cut_apply"); + xi_cut_mean = ps.getParameter("xi_cut_mean"); xi_cut_value = ps.getParameter("xi_cut_value"); + th_y_cut_apply = ps.getParameter("th_y_cut_apply"); + th_y_cut_mean = ps.getParameter("th_y_cut_mean"); th_y_cut_value = ps.getParameter("th_y_cut_value"); ti_tr_min = ps.getParameter("ti_tr_min"); @@ -88,12 +97,19 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> { edm::ParameterSetDescription desc; desc.add("x_cut_apply", false)->setComment("whether to apply track-association cut in x"); + desc.add("x_cut_mean", 0E-6)->setComment("mean of track-association cut in x, mm"); desc.add("x_cut_value", 800E-6)->setComment("threshold of track-association cut in x, mm"); + desc.add("y_cut_apply", false)->setComment("whether to apply track-association cut in y"); + desc.add("y_cut_mean", 0E-6)->setComment("mean of track-association cut in y, mm"); desc.add("y_cut_value", 600E-6)->setComment("threshold of track-association cut in y, mm"); + desc.add("xi_cut_apply", true)->setComment("whether to apply track-association cut in xi"); + desc.add("xi_cut_mean", 0.)->setComment("mean of track-association cut in xi"); desc.add("xi_cut_value", 0.013)->setComment("threshold of track-association cut in xi"); + desc.add("th_y_cut_apply", true)->setComment("whether to apply track-association cut in th_y"); + desc.add("th_y_cut_mean", 0E-6)->setComment("mean of track-association cut in th_y, rad"); desc.add("th_y_cut_value", 20E-6)->setComment("threshold of track-association cut in th_y, rad"); desc.add("ti_tr_min", -1.)->setComment("minimum value for timing-tracking association cut"); @@ -106,6 +122,7 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> { std::map association_cuts_; // map: arm -> AssociationCuts unsigned int max_n_timing_tracks_; + double default_time_; ProtonReconstructionAlgorithm algorithm_; @@ -117,6 +134,9 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> { CTPPSProtonProducer::CTPPSProtonProducer(const edm::ParameterSet &iConfig) : tracksToken_(consumes(iConfig.getParameter("tagLocalTrackLite"))), + + pixelDiscardBXShiftedTracks_(iConfig.getParameter("pixelDiscardBXShiftedTracks")), + lhcInfoLabel_(iConfig.getParameter("lhcInfoLabel")), opticsLabel_(iConfig.getParameter("opticsLabel")), verbosity_(iConfig.getUntrackedParameter("verbosity", 0)), @@ -131,6 +151,7 @@ CTPPSProtonProducer::CTPPSProtonProducer(const edm::ParameterSet &iConfig) localAngleYMax_(iConfig.getParameter("localAngleYMax")), max_n_timing_tracks_(iConfig.getParameter("max_n_timing_tracks")), + default_time_(iConfig.getParameter("default_time")), algorithm_(iConfig.getParameter("fitVtxY"), iConfig.getParameter("useImprovedInitialEstimate"), @@ -157,6 +178,9 @@ void CTPPSProtonProducer::fillDescriptions(edm::ConfigurationDescriptions &descr desc.add("tagLocalTrackLite", edm::InputTag("ctppsLocalTrackLiteProducer")) ->setComment("specification of the input lite-track collection"); + desc.add("pixelDiscardBXShiftedTracks", false) + ->setComment("whether to discard pixel tracks built from BX-shifted planes"); + desc.add("lhcInfoLabel", "")->setComment("label of the LHCInfo record"); desc.add("opticsLabel", "")->setComment("label of the optics record"); @@ -187,6 +211,8 @@ void CTPPSProtonProducer::fillDescriptions(edm::ConfigurationDescriptions &descr desc.add("max_n_timing_tracks", 5)->setComment("maximum number of timing tracks per RP"); + desc.add("default_time", 0.)->setComment("proton time to be used when no timing information available"); + desc.add("fitVtxY", true) ->setComment("for multi-RP reconstruction, flag whether y* should be free fit parameter"); @@ -253,6 +279,12 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe tr.ty() > localAngleYMax_) continue; + if (pixelDiscardBXShiftedTracks_) { + if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes || + tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes) + continue; + } + const CTPPSDetId rpId(tr.rpId()); if (verbosity_) @@ -298,32 +330,36 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe // do multi-RP reco if chosen if (doMultiRPReconstruction_ && rpIds.size() == 2) { - // find matching track pairs from different tracking RPs + // find matching track pairs from different tracking RPs, ordered: i=near, j=far RP std::vector> idx_pairs; std::map idx_pair_multiplicity; for (const auto &i : indices) { for (const auto &j : indices) { - if (j <= i) - continue; - const auto &tr_i = hTracks->at(i); const auto &tr_j = hTracks->at(j); + const double z_i = hGeometry->rpTranslation(tr_i.rpId()).z(); + const double z_j = hGeometry->rpTranslation(tr_j.rpId()).z(); + const auto &pr_i = singleRPResultsIndexed[i]; const auto &pr_j = singleRPResultsIndexed[j]; if (tr_i.rpId() == tr_j.rpId()) continue; + if (std::abs(z_i) >= std::abs(z_j)) + continue; + bool matching = true; - if (ac.x_cut_apply && std::abs(tr_i.x() - tr_j.x()) > ac.x_cut_value) + if (ac.x_cut_apply && std::abs(tr_i.x() - tr_j.x() - ac.x_cut_mean) > ac.x_cut_value) matching = false; - if (ac.y_cut_apply && std::abs(tr_i.y() - tr_j.y()) > ac.y_cut_value) + else if (ac.y_cut_apply && std::abs(tr_i.y() - tr_j.y() - ac.y_cut_mean) > ac.y_cut_value) matching = false; - if (ac.xi_cut_apply && std::abs(pr_i.xi() - pr_j.xi()) > ac.xi_cut_value) + else if (ac.xi_cut_apply && std::abs(pr_i.xi() - pr_j.xi() - ac.xi_cut_mean) > ac.xi_cut_value) matching = false; - if (ac.th_y_cut_apply && std::abs(pr_i.thetaY() - pr_j.thetaY()) > ac.th_y_cut_value) + else if (ac.th_y_cut_apply && + std::abs(pr_i.thetaY() - pr_j.thetaY() - ac.th_y_cut_mean) > ac.th_y_cut_value) matching = false; if (!matching) @@ -432,7 +468,7 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe swt += w * tr.time(); } - float time = 0., time_unc = 0.; + float time = default_time_, time_unc = 0.; if (sw > 0.) { time = swt / sw; time_unc = 1. / sqrt(sw); diff --git a/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py b/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py index 6dd251723c601..03b01f2299540 100644 --- a/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py +++ b/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py @@ -20,28 +20,73 @@ def applyDefaultSettings(ctppsProtons): ctppsProtons.association_cuts_45.xi_cut_apply = True ctppsProtons.association_cuts_45.xi_cut_value = 0.010 ctppsProtons.association_cuts_45.th_y_cut_apply = False + ctppsProtons.association_cuts_45.ti_tr_min = -1.5 + ctppsProtons.association_cuts_45.ti_tr_max = 2.0 ctppsProtons.association_cuts_56.x_cut_apply = False ctppsProtons.association_cuts_56.y_cut_apply = False ctppsProtons.association_cuts_56.xi_cut_apply = True ctppsProtons.association_cuts_56.xi_cut_value = 0.015 ctppsProtons.association_cuts_56.th_y_cut_apply = False + ctppsProtons.association_cuts_56.ti_tr_min = -1.5 + ctppsProtons.association_cuts_56.ti_tr_max = 2.0 + + ctppsProtons.pixelDiscardBXShiftedTracks = True + ctppsProtons.default_time = -999. ctpps_2016.toModify(ctppsProtons, applyDefaultSettings) # applied for all Run2 years (2016, 2017 and 2018) def apply2017Settings(ctppsProtons): - ctppsProtons.association_cuts_45.xi_cut_value = 0.010 - ctppsProtons.association_cuts_56.xi_cut_value = 0.015 + ctppsProtons.association_cuts_45.x_cut_apply = False + ctppsProtons.association_cuts_45.y_cut_apply = False + + ctppsProtons.association_cuts_45.xi_cut_apply = True + ctppsProtons.association_cuts_45.xi_cut_value = 5. * 0.00121 + ctppsProtons.association_cuts_45.xi_cut_mean = +6.0695e-5 + + ctppsProtons.association_cuts_45.th_y_cut_apply = False + + ctppsProtons.association_cuts_56.x_cut_apply = False + + ctppsProtons.association_cuts_56.y_cut_apply = True + ctppsProtons.association_cuts_56.y_cut_value = 5. * 0.14777 + ctppsProtons.association_cuts_56.y_cut_mean = -0.022612 + + ctppsProtons.association_cuts_56.xi_cut_apply = True + ctppsProtons.association_cuts_56.xi_cut_value = 5. * 0.0020627 + ctppsProtons.association_cuts_56.xi_cut_mean = +8.012857e-5 + + ctppsProtons.association_cuts_56.th_y_cut_apply = False ctpps_2017.toModify(ctppsProtons, apply2017Settings) def apply2018Settings(ctppsProtons): - ctppsProtons.association_cuts_45.xi_cut_value = 0.013 - ctppsProtons.association_cuts_45.th_y_cut_apply = True - ctppsProtons.association_cuts_45.th_y_cut_value = 30E-6 + ctppsProtons.association_cuts_45.x_cut_apply = True + ctppsProtons.association_cuts_45.x_cut_value = 4. * 0.16008188 + ctppsProtons.association_cuts_45.x_cut_mean = -0.065194856 + + ctppsProtons.association_cuts_45.y_cut_apply = True + ctppsProtons.association_cuts_45.y_cut_value = 4. * 0.1407986 + ctppsProtons.association_cuts_45.y_cut_mean = +0.10973631 - ctppsProtons.association_cuts_56.xi_cut_value = 0.013 - ctppsProtons.association_cuts_56.th_y_cut_apply = True - ctppsProtons.association_cuts_56.th_y_cut_value = 20E-6 + ctppsProtons.association_cuts_45.xi_cut_apply = True + ctppsProtons.association_cuts_45.xi_cut_value = 4. * 0.0012403586 + ctppsProtons.association_cuts_45.xi_cut_mean = +3.113062e-5 + + ctppsProtons.association_cuts_45.th_y_cut_apply = False + + ctppsProtons.association_cuts_56.x_cut_apply = True + ctppsProtons.association_cuts_56.x_cut_value = 5. * 0.18126434 + ctppsProtons.association_cuts_56.x_cut_mean = +0.073016431 + + ctppsProtons.association_cuts_56.y_cut_apply = True + ctppsProtons.association_cuts_56.y_cut_value = 5. * 0.14990802 + ctppsProtons.association_cuts_56.y_cut_mean = +0.064261029 + + ctppsProtons.association_cuts_56.xi_cut_apply = True + ctppsProtons.association_cuts_56.xi_cut_value = 5. * 0.002046409 + ctppsProtons.association_cuts_56.xi_cut_mean = -1.1852528e-5 + + ctppsProtons.association_cuts_56.th_y_cut_apply = False ctpps_2018.toModify(ctppsProtons, apply2018Settings) diff --git a/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc b/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc index 7af961a039085..f54408ecc6b24 100644 --- a/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc +++ b/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc @@ -382,7 +382,8 @@ reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromMultiRP(const // print results if (verbosity_) - os << " xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2 + os << " fit valid=" << valid << std::endl + << " xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2 << std::endl; // save reco candidate diff --git a/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc b/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc index 4f53821b2906d..0ae2082463077 100644 --- a/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc +++ b/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc @@ -192,7 +192,7 @@ void CTPPSDirectProtonSimulation::produce(edm::Event &iEvent, const edm::EventSe std::unique_ptr> pTracks(new std::vector()); // loop over event vertices - auto evt = new HepMC::GenEvent(*hepmc_prod->GetEvent()); + auto evt = hepmc_prod->GetEvent(); for (auto it_vtx = evt->vertices_begin(); it_vtx != evt->vertices_end(); ++it_vtx) { auto vtx = *(it_vtx); diff --git a/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc b/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc index 99379d6564e0e..977b47235c240 100644 --- a/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc +++ b/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc @@ -34,6 +34,9 @@ class CTPPSLHCInfoPlotter : public edm::one::EDAnalyzer<> { TH1D *h_beamEnergy_; TH1D *h_xangle_; TH1D *h_betaStar_; + + TH1D *h_fill_; + TH1D *h_run_; }; //---------------------------------------------------------------------------------------------------- @@ -49,7 +52,10 @@ CTPPSLHCInfoPlotter::CTPPSLHCInfoPlotter(const edm::ParameterSet &iConfig) h_beamEnergy_(new TH1D("h_beamEnergy", ";beam energy (GeV)", 81, -50., 8050.)), h_xangle_(new TH1D("h_xangle", ";(half) crossing angle (#murad)", 201, -0.5, 200.5)), - h_betaStar_(new TH1D("h_betaStar", ";#beta^{*} (m)", 101, -0.005, 1.005)) {} + h_betaStar_(new TH1D("h_betaStar", ";#beta^{*} (m)", 101, -0.005, 1.005)), + + h_fill_(new TH1D("h_fill", ";fill", 4001, 3999.5, 8000.5)), + h_run_(new TH1D("h_run", ";run", 6000, 270E3, 330E3)) {} //---------------------------------------------------------------------------------------------------- @@ -60,6 +66,9 @@ void CTPPSLHCInfoPlotter::analyze(const edm::Event &iEvent, const edm::EventSetu h_beamEnergy_->Fill(hLHCInfo->energy()); h_xangle_->Fill(hLHCInfo->crossingAngle()); h_betaStar_->Fill(hLHCInfo->betaStar()); + + h_fill_->Fill(hLHCInfo->fillNumber()); + h_run_->Fill(iEvent.id().run()); } //---------------------------------------------------------------------------------------------------- @@ -70,6 +79,9 @@ void CTPPSLHCInfoPlotter::endJob() { h_beamEnergy_->Write(); h_xangle_->Write(); h_betaStar_->Write(); + + h_fill_->Write(); + h_run_->Write(); } //---------------------------------------------------------------------------------------------------- diff --git a/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc b/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc index a4a37399e1573..158ae44ae16d6 100644 --- a/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc +++ b/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc @@ -34,6 +34,9 @@ class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> { std::string opticsLabel_; + unsigned int rpId_45_N_, rpId_45_F_; + unsigned int rpId_56_N_, rpId_56_F_; + std::string outputFile_; struct RPPlots { @@ -75,6 +78,24 @@ class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> { }; std::map rp_plots_; + + struct ArmPlots { + unsigned int id_N, id_F; + + std::unique_ptr g_de_x_vs_x_disp, g_de_y_vs_x_disp; + + ArmPlots() : g_de_x_vs_x_disp(new TGraph), g_de_y_vs_x_disp(new TGraph) {} + + void write() const { + g_de_x_vs_x_disp->SetTitle(";x_N (cm);x_F - x_N (cm)"); + g_de_x_vs_x_disp->Write("g_de_x_vs_x_disp"); + + g_de_y_vs_x_disp->SetTitle(";x_N (cm);y_F - y_N (cm)"); + g_de_y_vs_x_disp->Write("g_de_y_vs_x_disp"); + } + }; + + std::map arm_plots_; }; //---------------------------------------------------------------------------------------------------- @@ -86,7 +107,19 @@ using namespace edm; CTPPSOpticsPlotter::CTPPSOpticsPlotter(const edm::ParameterSet& iConfig) : opticsLabel_(iConfig.getParameter("opticsLabel")), - outputFile_(iConfig.getParameter("outputFile")) {} + + rpId_45_N_(iConfig.getParameter("rpId_45_N")), + rpId_45_F_(iConfig.getParameter("rpId_45_F")), + rpId_56_N_(iConfig.getParameter("rpId_56_N")), + rpId_56_F_(iConfig.getParameter("rpId_56_F")), + + outputFile_(iConfig.getParameter("outputFile")) { + arm_plots_[0].id_N = rpId_45_N_; + arm_plots_[0].id_F = rpId_45_F_; + + arm_plots_[1].id_N = rpId_56_N_; + arm_plots_[1].id_F = rpId_56_F_; +} //---------------------------------------------------------------------------------------------------- @@ -103,7 +136,7 @@ void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup if (hOpticalFunctions->empty()) return; - // make plots + // make per-RP plots for (const auto& it : *hOpticalFunctions) { CTPPSDetId rpId(it.first); unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); @@ -151,6 +184,48 @@ void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup pl.h_y_vs_x_disp->SetPoint(idx, k_out_xi.x - k_out_beam.x, k_out_xi.y - k_out_beam.y); } } + + // make per-arm plots + for (const auto& ap : arm_plots_) { + // find optics objects + const LHCInterpolatedOpticalFunctionsSet *opt_N = nullptr, *opt_F = nullptr; + + for (const auto& it : *hOpticalFunctions) { + CTPPSDetId rpId(it.first); + unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + if (rpDecId == ap.second.id_N) + opt_N = &it.second; + if (rpDecId == ap.second.id_F) + opt_F = &it.second; + } + + if (!opt_N || !opt_F) { + edm::LogError("CTPPSOpticsPlotter::analyze") << "Cannot find optics objects for arm " << ap.first; + continue; + } + + LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.}; + + LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam_N, k_out_beam_F; + opt_N->transport(k_in_beam, k_out_beam_N); + opt_F->transport(k_in_beam, k_out_beam_F); + + for (double xi = 0.; xi < 0.30001; xi += 0.001) { + LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi}; + + LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_N, k_out_xi_F; + opt_N->transport(k_in_xi, k_out_xi_N); + opt_F->transport(k_in_xi, k_out_xi_F); + + int idx = ap.second.g_de_x_vs_x_disp->GetN(); + + ap.second.g_de_x_vs_x_disp->SetPoint( + idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.x - k_out_beam_F.x) - (k_out_xi_N.x - k_out_beam_N.x)); + ap.second.g_de_y_vs_x_disp->SetPoint( + idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.y - k_out_beam_F.y) - (k_out_xi_N.y - k_out_beam_N.y)); + } + } } //---------------------------------------------------------------------------------------------------- @@ -162,6 +237,11 @@ void CTPPSOpticsPlotter::endJob() { gDirectory = f_out->mkdir(Form("%u", p.first)); p.second.write(); } + + for (const auto& p : arm_plots_) { + gDirectory = f_out->mkdir(Form("arm %u", p.first)); + p.second.write(); + } } //---------------------------------------------------------------------------------------------------- diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc new file mode 100644 index 0000000000000..754d5e7de1cc7 --- /dev/null +++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc @@ -0,0 +1,681 @@ +/**************************************************************************** + * Authors: + * Jan Kašpar + * Mauricio Thiel + ****************************************************************************/ + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ESWatcher.h" + +#include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h" +#include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h" + +#include "DataFormats/ProtonReco/interface/ForwardProton.h" +#include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h" +#include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h" +#include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h" + +#include "TFile.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TProfile.h" +#include "TF1.h" +#include "TGraph.h" + +#include +#include +#include + +//---------------------------------------------------------------------------------------------------- + +class CTPPSProtonReconstructionEfficiencyEstimatorData : public edm::one::EDAnalyzer<> { +public: + explicit CTPPSProtonReconstructionEfficiencyEstimatorData(const edm::ParameterSet &); + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void analyze(const edm::Event &, const edm::EventSetup &) override; + void endJob() override; + + edm::EDGetTokenT tokenTracks_; + + edm::EDGetTokenT tokenRecoProtonsMultiRP_; + + bool pixelDiscardBXShiftedTracks_; + + double localAngleXMin_, localAngleXMax_, localAngleYMin_, localAngleYMax_; + + std::string opticsLabel_; + + unsigned int n_prep_events_; + + unsigned int n_exp_prot_max_; + + std::vector n_sigmas_; + + std::string outputFile_; + + unsigned int verbosity_; + + edm::ESWatcher opticsWatcher_; + + struct ArmData { + unsigned int rpId_N, rpId_F; + + std::shared_ptr s_x_to_xi_N; + + unsigned int n_events_with_tracks; + + std::unique_ptr h_de_x, h_de_y; + std::unique_ptr h2_de_y_vs_de_x; + + double de_x_mean, de_x_sigma; + double de_y_mean, de_y_sigma; + + std::vector n_sigmas; + + unsigned int n_exp_prot_max; + + struct EffPlots { + std::unique_ptr p_eff1_vs_x_N; + std::unique_ptr p_eff1_vs_xi_N; + + std::unique_ptr p_eff2_vs_x_N; + std::unique_ptr p_eff2_vs_xi_N; + + EffPlots() + : p_eff1_vs_x_N(new TProfile("", ";x_{N} (mm);efficiency", 50, 0., 25.)), + p_eff1_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)), + + p_eff2_vs_x_N(new TProfile("", ";x_{N} (mm);efficiency", 50, 0., 25.)), + p_eff2_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)) {} + + void write() const { + p_eff1_vs_x_N->Write("p_eff1_vs_x_N"); + p_eff1_vs_xi_N->Write("p_eff1_vs_xi_N"); + + p_eff2_vs_x_N->Write("p_eff2_vs_x_N"); + p_eff2_vs_xi_N->Write("p_eff2_vs_xi_N"); + } + }; + + // (n exp protons, index in n_sigmas) --> plots + // n exp protons = 0 --> no condition (summary case) + std::map> effPlots; + + ArmData() + : n_events_with_tracks(0), + + h_de_x(new TH1D("", ";x_{F} - x_{N} (mm)", 200, -1., +1.)), + h_de_y(new TH1D("", ";y_{F} - y_{N} (mm)", 200, -1., +1.)), + h2_de_y_vs_de_x(new TH2D("", ";x_{F} - x_{N} (mm);y_{F} - y_{N} (mm)", 100, -1., +1., 100, -1., +1.)), + + de_x_mean(0.), + de_x_sigma(0.), + de_y_mean(0.), + de_y_sigma(0.) { + for (unsigned int nep = 0; nep <= n_exp_prot_max; ++nep) { + for (unsigned int nsi = 0; nsi < n_sigmas.size(); ++nsi) { + effPlots[nep][nsi] = EffPlots(); + } + } + } + + void write() const { + h_de_x->Write("h_de_x"); + h_de_y->Write("h_de_y"); + h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x"); + + char buf[100]; + + for (const auto &n_si : n_sigmas) { + sprintf(buf, "g_de_y_vs_de_x_n_si=%.1f", n_si); + TGraph *g = new TGraph(); + g->SetName(buf); + + g->SetPoint(0, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma); + g->SetPoint(1, de_x_mean + n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma); + g->SetPoint(2, de_x_mean + n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma); + g->SetPoint(3, de_x_mean - n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma); + g->SetPoint(4, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma); + + g->Write(); + } + + TDirectory *d_top = gDirectory; + + for (const auto &nep_p : effPlots) { + if (nep_p.first == 0) + sprintf(buf, "exp prot all"); + else + sprintf(buf, "exp prot %u", nep_p.first); + + TDirectory *d_nep = d_top->mkdir(buf); + + for (const auto &nsi_p : nep_p.second) { + sprintf(buf, "nsi = %.1f", n_sigmas[nsi_p.first]); + + TDirectory *d_nsi = d_nep->mkdir(buf); + + gDirectory = d_nsi; + + nsi_p.second.write(); + } + } + + gDirectory = d_top; + } + + void UpdateOptics(const LHCInterpolatedOpticalFunctionsSetCollection &ofc) { + const LHCInterpolatedOpticalFunctionsSet *of = nullptr; + + for (const auto &p : ofc) { + CTPPSDetId rpId(p.first); + unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + if (decRPId == rpId_N) { + of = &p.second; + break; + } + } + + if (!of) { + edm::LogError("ArmData::UpdateOptics") << "Cannot find optical functions for RP " << rpId_N; + return; + } + + std::vector xiValues = + of->getXiValues(); // local copy made since the TSpline constructor needs non-const parameters + std::vector xDValues = of->getFcnValues()[LHCOpticalFunctionsSet::exd]; + s_x_to_xi_N = std::make_shared("", xDValues.data(), xiValues.data(), xiValues.size()); + } + }; + + std::map data_; + + std::unique_ptr ff_; +}; + +//---------------------------------------------------------------------------------------------------- + +using namespace std; +using namespace edm; + +//---------------------------------------------------------------------------------------------------- + +CTPPSProtonReconstructionEfficiencyEstimatorData::CTPPSProtonReconstructionEfficiencyEstimatorData( + const edm::ParameterSet &iConfig) + : tokenTracks_(consumes(iConfig.getParameter("tagTracks"))), + + tokenRecoProtonsMultiRP_( + consumes(iConfig.getParameter("tagRecoProtonsMultiRP"))), + + pixelDiscardBXShiftedTracks_(iConfig.getParameter("pixelDiscardBXShiftedTracks")), + + localAngleXMin_(iConfig.getParameter("localAngleXMin")), + localAngleXMax_(iConfig.getParameter("localAngleXMax")), + localAngleYMin_(iConfig.getParameter("localAngleYMin")), + localAngleYMax_(iConfig.getParameter("localAngleYMax")), + + opticsLabel_(iConfig.getParameter("opticsLabel")), + n_prep_events_(iConfig.getParameter("n_prep_events")), + n_exp_prot_max_(iConfig.getParameter("n_exp_prot_max")), + n_sigmas_(iConfig.getParameter>("n_sigmas")), + + outputFile_(iConfig.getParameter("outputFile")), + + verbosity_(iConfig.getUntrackedParameter("verbosity")), + + ff_(new TF1("ff", "[0] * exp(- (x-[1])*(x-[1]) / 2 / [2]/[2]) + [4]")) { + data_[0].n_exp_prot_max = n_exp_prot_max_; + data_[0].n_sigmas = n_sigmas_; + data_[0].rpId_N = iConfig.getParameter("rpId_45_N"); + data_[0].rpId_F = iConfig.getParameter("rpId_45_F"); + + data_[1].n_exp_prot_max = n_exp_prot_max_; + data_[1].n_sigmas = n_sigmas_; + data_[1].rpId_N = iConfig.getParameter("rpId_56_N"); + data_[1].rpId_F = iConfig.getParameter("rpId_56_F"); +} + +//---------------------------------------------------------------------------------------------------- + +void CTPPSProtonReconstructionEfficiencyEstimatorData::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + + desc.add("tagTracks", edm::InputTag())->setComment("input tag for local lite tracks"); + desc.add("tagRecoProtonsMultiRP", edm::InputTag())->setComment("input tag for multi-RP reco protons"); + + desc.add("pixelDiscardBXShiftedTracks", false) + ->setComment("whether to discard pixel tracks built from BX-shifted planes"); + + desc.add("localAngleXMin", -0.03)->setComment("minimal accepted value of local horizontal angle (rad)"); + desc.add("localAngleXMax", +0.03)->setComment("maximal accepted value of local horizontal angle (rad)"); + desc.add("localAngleYMin", -0.04)->setComment("minimal accepted value of local vertical angle (rad)"); + desc.add("localAngleYMax", +0.04)->setComment("maximal accepted value of local vertical angle (rad)"); + + desc.add("opticsLabel", "")->setComment("label of the optics records"); + + desc.add("n_prep_events", 1000) + ->setComment("number of preparatory events (to determine de x and de y window)"); + + desc.add("n_exp_prot_max", 5)->setComment("maximum number of expected protons per event and per arm"); + + desc.add>("n_sigmas", {3., 5., 7.})->setComment("list of n_sigma values"); + + desc.add("rpId_45_N", 0)->setComment("decimal RP id for 45 near"); + desc.add("rpId_45_F", 0)->setComment("decimal RP id for 45 far"); + desc.add("rpId_56_N", 0)->setComment("decimal RP id for 56 near"); + desc.add("rpId_56_F", 0)->setComment("decimal RP id for 56 far"); + + desc.add("outputFile", "output.root")->setComment("output file name"); + + desc.addUntracked("verbosity", 0)->setComment("verbosity level"); + + descriptions.add("ctppsProtonReconstructionEfficiencyEstimatorData", desc); +} + +//---------------------------------------------------------------------------------------------------- + +void CTPPSProtonReconstructionEfficiencyEstimatorData::analyze(const edm::Event &iEvent, + const edm::EventSetup &iSetup) { + std::ostringstream os; + + // get conditions + edm::ESHandle hOpticalFunctions; + iSetup.get().get(opticsLabel_, hOpticalFunctions); + + // check optics change + if (opticsWatcher_.check(iSetup)) { + data_[0].UpdateOptics(*hOpticalFunctions); + data_[1].UpdateOptics(*hOpticalFunctions); + } + + // get input + edm::Handle hTracks; + iEvent.getByToken(tokenTracks_, hTracks); + + Handle hRecoProtonsMultiRP; + iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP); + + // pre-selection of tracks + std::vector sel_track_idc; + for (unsigned int idx = 0; idx < hTracks->size(); ++idx) { + const auto &tr = hTracks->at(idx); + + if (tr.tx() < localAngleXMin_ || tr.tx() > localAngleXMax_ || tr.ty() < localAngleYMin_ || + tr.ty() > localAngleYMax_) + continue; + + if (pixelDiscardBXShiftedTracks_) { + if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes || + tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes) + continue; + } + + sel_track_idc.push_back(idx); + } + + // debug print + if (verbosity_ > 1) { + os << "* tracks (pre-selected):" << std::endl; + + for (const auto idx : sel_track_idc) { + const auto &tr = hTracks->at(idx); + CTPPSDetId rpId(tr.rpId()); + unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + os << " [" << idx << "] RP=" << decRPId << ", x=" << tr.x() << ", y=" << tr.y() << std::endl; + } + + os << "* protons:" << std::endl; + for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) { + const auto &pr = (*hRecoProtonsMultiRP)[i]; + os << " [" << i << "] track indices: "; + for (const auto &trr : pr.contributingLocalTracks()) + os << trr.key() << ", "; + os << std::endl; + } + } + + // make de_x and de_y plots + map hasTracksInArm; + + for (const auto idx_i : sel_track_idc) { + const auto &tr_i = hTracks->at(idx_i); + CTPPSDetId rpId_i(tr_i.rpId()); + unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp(); + + for (const auto idx_j : sel_track_idc) { + const auto &tr_j = hTracks->at(idx_j); + CTPPSDetId rpId_j(tr_j.rpId()); + unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp(); + + // check whether desired RP combination + unsigned int arm = 123; + for (const auto &ap : data_) { + if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j) + arm = ap.first; + } + + if (arm > 1) + continue; + + // update flag + hasTracksInArm[arm] = true; + + // update plots + auto &ad = data_[arm]; + const double de_x = tr_j.x() - tr_i.x(); + const double de_y = tr_j.y() - tr_i.y(); + + if (ad.n_events_with_tracks < n_prep_events_) { + ad.h_de_x->Fill(de_x); + ad.h_de_y->Fill(de_y); + } + + ad.h2_de_y_vs_de_x->Fill(de_x, de_y); + } + } + + // update event counter + for (auto &ap : data_) { + if (hasTracksInArm[ap.first]) + ap.second.n_events_with_tracks++; + } + + // if threshold reached do fits + for (auto &ap : data_) { + auto &ad = ap.second; + + if (ad.n_events_with_tracks == n_prep_events_) { + if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.) + continue; + + std::vector> m; + m.emplace_back(0, ad.h_de_x.get()); + m.emplace_back(1, ad.h_de_y.get()); + + for (const auto &p : m) { + double max_pos = -1E100, max_val = -1.; + for (int bi = 1; bi < p.second->GetNbinsX(); ++bi) { + const double pos = p.second->GetBinCenter(bi); + const double val = p.second->GetBinContent(bi); + + if (val > max_val) { + max_val = val; + max_pos = pos; + } + } + + const double sig = 0.2; + + ff_->SetParameters(max_val, max_pos, sig, 0.); + p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig); + p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig); + + if (p.first == 0) { + ad.de_x_mean = ff_->GetParameter(1); + ad.de_x_sigma = fabs(ff_->GetParameter(2)); + } + if (p.first == 1) { + ad.de_y_mean = ff_->GetParameter(1); + ad.de_y_sigma = fabs(ff_->GetParameter(2)); + } + } + + if (verbosity_) { + os << "* fitting arm " << ap.first << std::endl + << " de_x: mean = " << ad.de_x_mean << ", sigma = " << ad.de_x_sigma << std::endl + << " de_y: mean = " << ad.de_y_mean << ", sigma = " << ad.de_y_sigma; + } + } + } + + // data structures for efficiency analysis + struct ArmEventData { + std::map> matched_track_idc; + + std::set reco_proton_idc; + + std::map> matched_track_with_prot_idc, matched_track_without_prot_idc; + }; + + std::map armEventData; + + // determine the number of expected protons + for (const auto idx_i : sel_track_idc) { + const auto &tr_i = hTracks->at(idx_i); + CTPPSDetId rpId_i(tr_i.rpId()); + unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp(); + + for (const auto idx_j : sel_track_idc) { + const auto &tr_j = hTracks->at(idx_j); + CTPPSDetId rpId_j(tr_j.rpId()); + unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp(); + + // check whether desired RP combination + unsigned int arm = 123; + for (const auto &ap : data_) { + if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j) + arm = ap.first; + } + + if (arm > 1) + continue; + + // check near-far matching + auto &ad = data_[arm]; + const double de_x = tr_j.x() - tr_i.x(); + const double de_y = tr_j.y() - tr_i.y(); + + for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) { + const double n_si = ad.n_sigmas[nsi]; + const bool match_x = fabs(de_x - ad.de_x_mean) < n_si * ad.de_x_sigma; + const bool match_y = fabs(de_y - ad.de_y_mean) < n_si * ad.de_y_sigma; + if (match_x && match_y) + armEventData[arm].matched_track_idc[nsi].insert(idx_i); + } + } + } + + // determine the number of reconstructed protons + for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) { + const auto &proton = (*hRecoProtonsMultiRP)[i]; + + CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId()); + unsigned int arm = rpId.arm(); + + if (proton.validFit()) + armEventData[arm].reco_proton_idc.insert(i); + } + + // compare matched tracks with reco protons + if (verbosity_ > 1) + os << "* cmp matched tracks vs. reco protons" << std::endl; + + for (auto &ap : armEventData) { + auto &ad = data_[ap.first]; + + if (verbosity_ > 1) + os << " arm " << ap.first << std::endl; + + for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) { + if (verbosity_ > 1) + os << " nsi = " << nsi << std::endl; + + for (const auto &tri : ap.second.matched_track_idc[nsi]) { + const auto &track = hTracks->at(tri); + + bool some_proton_matching = false; + + if (verbosity_ > 1) + os << " tri = " << tri << std::endl; + + for (const auto &pri : ap.second.reco_proton_idc) { + const auto &proton = (*hRecoProtonsMultiRP)[pri]; + + bool proton_matching = false; + + for (const auto &pr_tr : proton.contributingLocalTracks()) { + CTPPSDetId rpId(pr_tr->rpId()); + unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); + + if (decRPId != ad.rpId_N) + continue; + + const double x = pr_tr->x(); + const double y = pr_tr->y(); + const double th = 1E-3; // 1 um + + const bool match = (fabs(x - track.x()) < th) && (fabs(y - track.y()) < th); + + if (verbosity_ > 1) + os << " pri = " << pri << ": x_tr = " << track.x() << ", x_pr = " << x + << ", match = " << match << std::endl; + + if (match) { + proton_matching = true; + break; + } + } + + if (proton_matching) { + some_proton_matching = true; + break; + } + } + + if (verbosity_ > 1) + os << " --> some_proton_matching " << some_proton_matching << std::endl; + + if (some_proton_matching) + ap.second.matched_track_with_prot_idc[nsi].insert(tri); + else + ap.second.matched_track_without_prot_idc[nsi].insert(tri); + } + } + } + + // debug print + if (verbosity_ > 1) { + for (auto &ap : armEventData) { + auto &ad = data_[ap.first]; + + if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.) + continue; + + os << "* results for arm " << ap.first << std::endl; + + os << " reco_proton_idc: "; + for (const auto &idx : ap.second.reco_proton_idc) + os << idx << ", "; + os << std::endl; + + for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) { + os << " n_si = " << ad.n_sigmas[nsi] << std::endl; + + os << " matched_track_idc: "; + for (const auto &idx : ap.second.matched_track_idc[nsi]) + os << idx << ", "; + os << std::endl; + + os << " matched_track_with_prot_idc: "; + for (const auto &idx : ap.second.matched_track_with_prot_idc[nsi]) + os << idx << ", "; + os << std::endl; + + os << " matched_track_without_prot_idc: "; + for (const auto &idx : ap.second.matched_track_without_prot_idc[nsi]) + os << idx << ", "; + os << std::endl; + } + } + } + + // update efficiency plots + for (auto &ap : armEventData) { + auto &ad = data_[ap.first]; + + // stop if sigmas not yet determined + if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.) + continue; + + // loop over n_sigma choices + for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) { + const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size(); + const unsigned int n_rec_prot = ap.second.reco_proton_idc.size(); + + // stop if N(expected protons) out of range + if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max) + continue; + + // update method 1 plots + const double eff = double(n_rec_prot) / n_exp_prot; + + for (unsigned int tri : ap.second.matched_track_idc[nsi]) { + const double x_N = hTracks->at(tri).x(); + const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1); // conversion mm to cm + + ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff); + ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff); + + ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff); + ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff); + } + + // update method 2 plots + for (const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) { + const double x_N = hTracks->at(tri).x(); + const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1); // conversion mm to cm + + ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.); + ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.); + + ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.); + ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.); + } + + for (const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) { + const double x_N = hTracks->at(tri).x(); + const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1); // conversion mm to cm + + ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.); + ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.); + + ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.); + ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.); + } + } + } + + if (verbosity_) + edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str(); +} + +//---------------------------------------------------------------------------------------------------- + +void CTPPSProtonReconstructionEfficiencyEstimatorData::endJob() { + auto f_out = std::make_unique(outputFile_.c_str(), "recreate"); + + for (const auto &ait : data_) { + char buf[100]; + sprintf(buf, "arm %u", ait.first); + TDirectory *d_arm = f_out->mkdir(buf); + gDirectory = d_arm; + + ait.second.write(); + } +} + +//---------------------------------------------------------------------------------------------------- + +DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorData); diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc similarity index 86% rename from Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc rename to Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc index a96e9dc24307c..0b6fb03619c41 100644 --- a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc +++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc @@ -36,9 +36,9 @@ //---------------------------------------------------------------------------------------------------- -class CTPPSProtonReconstructionEfficiencyEstimator : public edm::one::EDAnalyzer<> { +class CTPPSProtonReconstructionEfficiencyEstimatorMC : public edm::one::EDAnalyzer<> { public: - explicit CTPPSProtonReconstructionEfficiencyEstimator(const edm::ParameterSet &); + explicit CTPPSProtonReconstructionEfficiencyEstimatorMC(const edm::ParameterSet &); private: void analyze(const edm::Event &, const edm::EventSetup &) override; @@ -62,6 +62,8 @@ class CTPPSProtonReconstructionEfficiencyEstimator : public edm::one::EDAnalyzer std::string outputFile_; + unsigned int verbosity_; + struct PlotGroup { std::unique_ptr p_eff_vs_xi; @@ -90,7 +92,7 @@ using namespace HepMC; //---------------------------------------------------------------------------------------------------- -CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficiencyEstimator( +CTPPSProtonReconstructionEfficiencyEstimatorMC::CTPPSProtonReconstructionEfficiencyEstimatorMC( const edm::ParameterSet &iConfig) : tokenHepMCAfterSmearing_( consumes(iConfig.getParameter("tagHepMCAfterSmearing"))), @@ -108,7 +110,9 @@ CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficienc rpId_56_N_(iConfig.getParameter("rpId_56_N")), rpId_56_F_(iConfig.getParameter("rpId_56_F")), - outputFile_(iConfig.getParameter("outputFile")) { + outputFile_(iConfig.getParameter("outputFile")), + + verbosity_(iConfig.getUntrackedParameter("verbosity", 0)) { rpDecId_near_[0] = rpId_45_N_; rpDecId_far_[0] = rpId_45_F_; @@ -124,17 +128,8 @@ CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficienc //---------------------------------------------------------------------------------------------------- -void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { - bool verbosity = false; - - const auto eid = iEvent.id().event(); - //if (eid == 46 || eid == 741 || eid == 1649 || eid == 4223 || eid == 4279) - // verbosity = true; - - if (verbosity) { - printf("--------------------------------------------------\n"); - printf("event %llu\n", eid); - } +void CTPPSProtonReconstructionEfficiencyEstimatorMC::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { + std::ostringstream os; // get conditions edm::ESHandle hLHCInfo; @@ -220,8 +215,8 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv std::map isStripRPNear, isStripRPFar; for (auto &pp : particleInfo) { - if (verbosity) - printf("* barcode=%i, arm=%u, xi=%.3f\n", pp.first, pp.second.arm, pp.second.xi); + 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); @@ -247,18 +242,15 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv pp.second.inAcceptanceFar = true; } - if (verbosity) - printf(" RP %u: %u hits\n", rpDecId, rpp.second); + if (verbosity_) + os << " RP " << rpDecId << ": " << rpp.second << " hits" << std::endl; } pp.second.inAcceptance = pp.second.inAcceptanceNear && pp.second.inAcceptanceFar; - if (verbosity) { - printf(" inAcceptance: near=%u, far=%u, global=%u\n", - pp.second.inAcceptanceNear, - pp.second.inAcceptanceFar, - pp.second.inAcceptance); - } + if (verbosity_) + os << " inAcceptance: near=" << pp.second.inAcceptanceNear << ", far=" << pp.second.inAcceptanceFar + << ", global=" << pp.second.inAcceptance << std::endl; } // count particles in acceptance @@ -308,16 +300,10 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv const auto &npa = nParticlesInAcceptance[arm]; const auto &nrt = nReconstructedTracks[arm]; - if (verbosity) { - printf("* arm %u: nRecoProtons=%u (tracks near=%u, far=%u), nAcc=%u (near=%u, far=%u)\n", - arm, - nReconstructedProtons[arm], - nReconstructedTracks[arm].near, - nReconstructedTracks[arm].far, - npa.global, - npa.near, - npa.far); - } + 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) @@ -334,8 +320,8 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv const double eff = double(nReconstructedProtons[arm]) / npa.global; - if (verbosity) - printf(" eff=%.3f\n", eff); + if (verbosity_) + os << " eff=" << eff << std::endl; for (auto &pp : particleInfo) { if (pp.second.arm != arm || !pp.second.inAcceptance) @@ -344,11 +330,14 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv p.p_eff_vs_xi->Fill(pp.second.xi, eff); } } + + if (verbosity_) + edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str(); } //---------------------------------------------------------------------------------------------------- -void CTPPSProtonReconstructionEfficiencyEstimator::endJob() { +void CTPPSProtonReconstructionEfficiencyEstimatorMC::endJob() { auto f_out = std::make_unique(outputFile_.c_str(), "recreate"); for (const auto &ait : plots_) { @@ -368,4 +357,4 @@ void CTPPSProtonReconstructionEfficiencyEstimator::endJob() { //---------------------------------------------------------------------------------------------------- -DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimator); +DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorMC); diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc index 4faae5d3c4ae2..564c692740680 100644 --- a/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc +++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc @@ -46,6 +46,18 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { unsigned int rpId_45_N_, rpId_45_F_; unsigned int rpId_56_N_, rpId_56_F_; + struct AssociationCuts { + double ti_tr_min; + double ti_tr_max; + + void load(const edm::ParameterSet &ps) { + ti_tr_min = ps.getParameter("ti_tr_min"); + ti_tr_max = ps.getParameter("ti_tr_max"); + } + }; + + std::map association_cuts_; + std::string outputFile_; signed int maxNonEmptyEvents_; @@ -82,17 +94,19 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { std::unique_ptr p_th_y_vs_xi; std::map m_h_xi_nTracks; + std::unique_ptr h_xi_n1f1; SingleRPPlots() : h_multiplicity(new TH1D("", ";reconstructed protons", 11, -0.5, 10.5)), h_xi(new TH1D("", ";#xi", 100, 0., 0.3)), h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)), - p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3)) { + p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3)), + h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)) { for (unsigned int n = 2; n <= 10; ++n) m_h_xi_nTracks[n] = new TH1D(*h_xi); } - void fill(const reco::ForwardProton &p, unsigned int nTracks) { + void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) { if (p.validFit()) { h_xi->Fill(p.xi()); @@ -103,6 +117,9 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { auto it = m_h_xi_nTracks.find(nTracks); if (it != m_h_xi_nTracks.end()) it->second->Fill(p.xi()); + + if (n1f1) + h_xi_n1f1->Fill(p.xi()); } } @@ -123,6 +140,8 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { } gDirectory = d_top; + + h_xi_n1f1->Write("h_xi_n1f1"); } }; @@ -140,9 +159,11 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { std::unique_ptr h2_timing_tracks_vs_prot_mult; std::map m_h_xi_nTracks; + std::unique_ptr h_xi_n1f1; std::unique_ptr h_de_x_timing_vs_tracking, h_de_x_rel_timing_vs_tracking, h_de_x_match_timing_vs_tracking; - std::unique_ptr h_de_x_rel_timing_vs_tracking_ClCo; + std::unique_ptr h_de_x_timing_vs_tracking_ClCo, h_de_x_rel_timing_vs_tracking_ClCo, + h_de_x_match_timing_vs_tracking_ClCo; std::unique_ptr h2_y_vs_x_tt0_ClCo, h2_y_vs_x_tt1_ClCo, h2_y_vs_x_ttm_ClCo; @@ -150,7 +171,7 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { : h_multiplicity(new TH1D("", ";reconstructed protons per event", 11, -0.5, 10.5)), h_xi(new TH1D("", ";#xi", 100, 0., 0.3)), h_th_x(new TH1D("", ";#theta_{x} (rad)", 250, -500E-6, +500E-6)), - h_th_y(new TH1D("", ";#theta_{y} (rad)", 250, -500E-6, +500E-6)), + h_th_y(new TH1D("", ";#theta_{y} (rad)", 500, -1000E-6, +1000E-6)), h_vtx_y(new TH1D("", ";vtx_{y} (cm)", 100, -100E-3, +100E-3)), h_chi_sq(new TH1D("", ";#chi^{2}", 100, 0., 10.)), h_log_chi_sq(new TH1D("", ";log_{10} #chi^{2}", 100, -20., 5.)), @@ -166,10 +187,17 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { p_vtx_y_vs_xi(new TProfile("", ";#xi;vtx_{y} (cm)", 100, 0., 0.3)), h2_timing_tracks_vs_prot_mult( new TH2D("", ";reco protons per event;timing tracks per event", 11, -0.5, 10.5, 11, -0.5, 10.5)), + h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)), + h_de_x_timing_vs_tracking(new TH1D("", ";#Delta x (mm)", 200, -1., +1.)), h_de_x_rel_timing_vs_tracking(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)), h_de_x_match_timing_vs_tracking(new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)), + + h_de_x_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x (mm)", 200, -1., +1.)), h_de_x_rel_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)), + h_de_x_match_timing_vs_tracking_ClCo( + new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)), + h2_y_vs_x_tt0_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)), h2_y_vs_x_tt1_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)), h2_y_vs_x_ttm_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)) { @@ -191,7 +219,7 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { m_h_xi_nTracks[n] = new TH1D(*h_xi); } - void fill(const reco::ForwardProton &p, unsigned int nTracks) { + void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) { if (!p.validFit()) return; @@ -246,6 +274,9 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { auto it = m_h_xi_nTracks.find(nTracks); if (it != m_h_xi_nTracks.end()) it->second->Fill(p.xi()); + + if (n1f1) + h_xi_n1f1->Fill(p.xi()); } void write() const { @@ -306,10 +337,15 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> { gDirectory = d_top; + h_xi_n1f1->Write("h_xi_n1f1"); + h_de_x_timing_vs_tracking->Write("h_de_x_timing_vs_tracking"); h_de_x_rel_timing_vs_tracking->Write("h_de_x_rel_timing_vs_tracking"); h_de_x_match_timing_vs_tracking->Write("h_de_x_match_timing_vs_tracking"); + + h_de_x_timing_vs_tracking_ClCo->Write("h_de_x_timing_vs_tracking_ClCo"); h_de_x_rel_timing_vs_tracking_ClCo->Write("h_de_x_rel_timing_vs_tracking_ClCo"); + h_de_x_match_timing_vs_tracking_ClCo->Write("h_de_x_match_timing_vs_tracking_ClCo"); h2_y_vs_x_tt0_ClCo->Write("h2_y_vs_x_tt0_ClCo"); h2_y_vs_x_tt1_ClCo->Write("h2_y_vs_x_tt1_ClCo"); @@ -444,7 +480,12 @@ CTPPSProtonReconstructionPlotter::CTPPSProtonReconstructionPlotter(const edm::Pa p_y_L_diffNF_vs_y_L_N_(new TProfile("p_y_L_diffNF_vs_y_L_N", ";y_{LN};y_{LF} - y_{LN}", 100, -20., +20.)), p_y_R_diffNF_vs_y_R_N_(new TProfile("p_y_R_diffNF_vs_y_R_N", ";y_{RN};y_{RF} - y_{RN}", 100, -20., +20.)), - n_non_empty_events_(0) {} + n_non_empty_events_(0) { + for (const std::string §or : {"45", "56"}) { + const unsigned int arm = (sector == "45") ? 0 : 1; + association_cuts_[arm].load(ps.getParameterSet("association_cuts_" + sector)); + } +} //---------------------------------------------------------------------------------------------------- @@ -500,19 +541,28 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed const CTPPSLocalTrackLite *tr_R_N = nullptr; const CTPPSLocalTrackLite *tr_R_F = nullptr; std::map armTrackCounter, armTimingTrackCounter; + std::map armTrackCounter_N, armTrackCounter_F; for (const auto &tr : *hTracks) { CTPPSDetId rpId(tr.rpId()); unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); - if (decRPId == rpId_45_N_) + if (decRPId == rpId_45_N_) { tr_L_N = &tr; - if (decRPId == rpId_45_F_) + armTrackCounter_N[0]++; + } + if (decRPId == rpId_45_F_) { tr_L_F = &tr; - if (decRPId == rpId_56_N_) + armTrackCounter_F[0]++; + } + if (decRPId == rpId_56_N_) { tr_R_N = &tr; - if (decRPId == rpId_56_F_) + armTrackCounter_N[1]++; + } + if (decRPId == rpId_56_F_) { tr_R_F = &tr; + armTrackCounter_F[1]++; + } armTrackCounter[rpId.arm()]++; @@ -542,7 +592,10 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed for (const auto &proton : *hRecoProtonsSingleRP) { CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId()); unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); - singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()]); + + const bool n1f1 = (armTrackCounter_N[rpId.arm()] == 1 && armTrackCounter_F[rpId.arm()] == 1); + + singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()], n1f1); if (proton.validFit()) singleRPMultiplicity[decRPId]++; @@ -555,7 +608,10 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed for (const auto &proton : *hRecoProtonsMultiRP) { CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId()); unsigned int armId = rpId.arm(); - multiRPPlots_[armId].fill(proton, armTrackCounter[armId]); + + const bool n1f1 = (armTrackCounter_N[armId] == 1 && armTrackCounter_F[armId] == 1); + + multiRPPlots_[armId].fill(proton, armTrackCounter[armId], n1f1); if (proton.validFit()) multiRPMultiplicity[armId]++; @@ -594,15 +650,18 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed double de_x = 0., de_x_unc = 0.; CalculateTimingTrackingDistance(proton, tr, *hGeometry, de_x, de_x_unc); - double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10; + const double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10; + const auto &ac = association_cuts_[armId]; + const bool match = (ac.ti_tr_min <= fabs(rd) && fabs(rd) <= ac.ti_tr_max); pl.h_de_x_timing_vs_tracking->Fill(de_x); pl.h_de_x_rel_timing_vs_tracking->Fill(rd); - pl.h_de_x_match_timing_vs_tracking->Fill(fabs(de_x / de_x_unc) <= 1. ? 1. : 0.); + pl.h_de_x_match_timing_vs_tracking->Fill(match ? 1. : 0.); - if (clCo[armId]) { - if (armTimingTrackCounter[armId] == 1) - pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd); + if (clCo[armId] && armTimingTrackCounter[armId] == 1) { + pl.h_de_x_timing_vs_tracking_ClCo->Fill(de_x); + pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd); + pl.h_de_x_match_timing_vs_tracking_ClCo->Fill(match ? 1. : 0.); } } } @@ -645,13 +704,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed const reco::ForwardProton *p_s_N = nullptr, *p_s_F = nullptr; - /* - printf("multi-RP candidate: "); - for (const auto &tr_m : proton_m.contributingLocalTracks()) - printf("%u, ", tr_m.key()); - printf("\n"); - */ - for (const auto &proton_s : *hRecoProtonsSingleRP) { // skip if source of single-RP reco not included in multi-RP reco const auto key_s = proton_s.contributingLocalTracks()[0].key(); @@ -663,8 +715,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed } } - //printf(" key_s = %u --> compatible = %u\n", key_s, compatible); - if (!compatible) continue; @@ -673,8 +723,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed const unsigned int idx = rpId_s.arm() * 1000 + rpId_s.station() * 100 + rpId_s.rp() * 10 + rpId_s.arm(); singleMultiCorrelationPlots_[idx].fill(proton_s, proton_m); - //printf(" arm_s = %u, arm_m = %u\n", rpId_s.arm(), arm); - // select protons for arm-correlation plots const unsigned int rpDecId_s = rpId_s.arm() * 100 + rpId_s.station() * 10 + rpId_s.rp(); if (rpDecId_s == rpId_45_N_ || rpDecId_s == rpId_56_N_) diff --git a/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc b/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc index 732bc56feb9bd..611274ebb965a 100644 --- a/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc +++ b/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc @@ -23,6 +23,7 @@ #include "TH2D.h" #include "TProfile.h" #include "TProfile2D.h" +#include "TGraph.h" #include @@ -44,6 +45,9 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { std::string outputFile_; + unsigned int events_total_; + std::map events_per_arm_; + struct RPPlots { bool initialized; @@ -51,13 +55,14 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { std::unique_ptr p_y_vs_x; std::unique_ptr h_x; std::unique_ptr h_y; + std::unique_ptr h_time; RPPlots() : initialized(false) {} void init(bool pixel, double pitch) { const double bin_size_x = (pixel) ? pitch * cos(18.4 / 180. * M_PI) : 100E-3; - h2_y_vs_x.reset(new TH2D("", "", 300, -10., +70., 300, -30., +30.)); + h2_y_vs_x.reset(new TH2D("", "", 300, -10., +70., 600, -30., +30.)); p_y_vs_x.reset(new TProfile("", "", 300, -10., +70.)); int n_mi = ceil(10. / bin_size_x); @@ -67,14 +72,17 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { h_y.reset(new TH1D("", "", 300, -15., +15.)); + h_time.reset(new TH1D("", ";time", 500, -50., +50.)); + initialized = true; } - void fill(double x, double y) { + void fill(double x, double y, double time) { h2_y_vs_x->Fill(x, y); p_y_vs_x->Fill(x, y); h_x->Fill(x); h_y->Fill(y); + h_time->Fill(time); } void write() const { @@ -82,15 +90,26 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { p_y_vs_x->Write("p_y_vs_x"); h_x->Write("h_x"); h_y->Write("h_y"); + h_time->Write("h_time"); } }; std::map rpPlots; struct ArmPlots { + unsigned int rpId_N, rpId_F; + std::unique_ptr h_de_x, h_de_y; std::unique_ptr p_de_x_vs_x, p_de_y_vs_x; std::unique_ptr p2_de_x_vs_x_y, p2_de_y_vs_x_y; + std::unique_ptr h2_de_x_vs_x, h2_de_y_vs_x; + std::unique_ptr h2_de_y_vs_de_x; + + struct Stat { + unsigned int sN = 0, s1 = 0; + }; + + std::map> m_stat; ArmPlots() : h_de_x(new TH1D("", ";x^{F} - x^{N}", 100, -1., +1.)), @@ -98,7 +117,10 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { p_de_x_vs_x(new TProfile("", ";x^{N};x^{F} - x^{N}", 40, 0., 40.)), p_de_y_vs_x(new TProfile("", ";x^{N};y^{F} - y^{N}", 40, 0., 40.)), p2_de_x_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)), - p2_de_y_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)) {} + p2_de_y_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)), + h2_de_x_vs_x(new TH2D("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., 100, -1., +1.)), + h2_de_y_vs_x(new TH2D("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., 100, -1., +1.)), + h2_de_y_vs_de_x(new TH2D("", ";x^{F} - x^{N};y^{F} - y^{N}", 100, -1., +1., 100, -1., +1.)) {} void fill(double x_N, double y_N, double x_F, double y_F) { h_de_x->Fill(x_F - x_N); @@ -109,6 +131,11 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { p2_de_x_vs_x_y->Fill(x_N, y_N, x_F - x_N); p2_de_y_vs_x_y->Fill(x_N, y_N, y_F - y_N); + + h2_de_x_vs_x->Fill(x_N, x_F - x_N); + h2_de_y_vs_x->Fill(x_N, y_F - y_N); + + h2_de_y_vs_de_x->Fill(x_F - x_N, y_F - y_N); } void write() const { @@ -120,6 +147,27 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { p2_de_x_vs_x_y->Write("p2_de_x_vs_x_y"); p2_de_y_vs_x_y->Write("p2_de_y_vs_x_y"); + + h2_de_x_vs_x->Write("h2_de_x_vs_x"); + h2_de_y_vs_x->Write("h2_de_y_vs_x"); + + h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x"); + + for (const auto& rp : m_stat) { + TGraph* g = new TGraph(); + + char buf[100]; + sprintf(buf, "g_mean_track_mult_run_%u", rp.first); + g->SetName(buf); + + for (const auto& lsp : rp.second) { + const int idx = g->GetN(); + const double m = (lsp.second.s1 > 0) ? double(lsp.second.sN) / lsp.second.s1 : 0.; + g->SetPoint(idx, lsp.first, m); + } + + g->Write(); + } } }; @@ -131,7 +179,14 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> { CTPPSTrackDistributionPlotter::CTPPSTrackDistributionPlotter(const edm::ParameterSet& iConfig) : tracksToken_(consumes(iConfig.getParameter("tagTracks"))), x_pitch_pixels_(iConfig.getUntrackedParameter("x_pitch_pixels", 150E-3)), - outputFile_(iConfig.getParameter("outputFile")) {} + outputFile_(iConfig.getParameter("outputFile")), + events_total_(0) { + armPlots[0].rpId_N = iConfig.getParameter("rpId_45_N"); + armPlots[0].rpId_F = iConfig.getParameter("rpId_45_F"); + + armPlots[1].rpId_N = iConfig.getParameter("rpId_56_N"); + armPlots[1].rpId_F = iConfig.getParameter("rpId_56_F"); +} //---------------------------------------------------------------------------------------------------- @@ -141,6 +196,8 @@ void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm: iEvent.getByToken(tracksToken_, tracks); // process tracks + std::map m_mult; + for (const auto& trk : *tracks) { CTPPSDetId rpId(trk.rpId()); unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp(); @@ -150,27 +207,53 @@ void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm: if (!pl.initialized) pl.init(rpPixel, x_pitch_pixels_); - pl.fill(trk.x(), trk.y()); + pl.fill(trk.x(), trk.y(), trk.time()); + + m_mult[rpId.arm()]++; + } + + for (unsigned int arm = 0; arm < 2; ++arm) { + auto& st = armPlots[arm].m_stat[iEvent.id().run()][iEvent.id().luminosityBlock()]; + st.s1++; + st.sN += m_mult[arm]; } for (const auto& t1 : *tracks) { - CTPPSDetId rpId1(t1.rpId()); + const CTPPSDetId rpId1(t1.rpId()); for (const auto& t2 : *tracks) { - CTPPSDetId rpId2(t2.rpId()); + const CTPPSDetId rpId2(t2.rpId()); if (rpId1.arm() != rpId2.arm()) continue; - if (rpId1.station() == 0 && rpId2.station() == 2) - armPlots[rpId1.arm()].fill(t1.x(), t1.y(), t2.x(), t2.y()); + auto& ap = armPlots[rpId1.arm()]; + + const unsigned int rpDecId1 = rpId1.arm() * 100 + rpId1.station() * 10 + rpId1.rp(); + const unsigned int rpDecId2 = rpId2.arm() * 100 + rpId2.station() * 10 + rpId2.rp(); + + if (rpDecId1 == ap.rpId_N && rpDecId2 == ap.rpId_F) + ap.fill(t1.x(), t1.y(), t2.x(), t2.y()); } } + + // update counters + events_total_++; + + if (m_mult[0] > 0) + events_per_arm_[0]++; + if (m_mult[1] > 0) + events_per_arm_[1]++; } //---------------------------------------------------------------------------------------------------- void CTPPSTrackDistributionPlotter::endJob() { + edm::LogInfo("CTPPSTrackDistributionPlotter") + << " events processed: " << events_total_ << " (" << std::scientific << double(events_total_) << ")\n" + << " events with tracks in sector 45: " << events_per_arm_[0] << " (" << double(events_per_arm_[0]) << ")\n" + << " events with tracks in sector 56: " << events_per_arm_[1] << " (" << double(events_per_arm_[1]) << ")"; + auto f_out = std::make_unique(outputFile_.c_str(), "recreate"); for (const auto& it : rpPlots) {