CTPPS/plugins/CTPPSProtonReconstructionValidator.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/ESWatcher.h"
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014
0015 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0016
0017 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0018 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0019
0020 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0021 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0022 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0023
0024 #include "TFile.h"
0025 #include "TH1D.h"
0026
0027 #include <map>
0028 #include <string>
0029
0030 //----------------------------------------------------------------------------------------------------
0031
0032 class CTPPSProtonReconstructionValidator : public edm::one::EDAnalyzer<> {
0033 public:
0034   explicit CTPPSProtonReconstructionValidator(const edm::ParameterSet&);
0035
0036 private:
0037   void analyze(const edm::Event&, const edm::EventSetup&) override;
0038   void endJob() override;
0039
0040   edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtons_;
0041   edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsESToken_;
0042   double chiSqCut_;
0043   std::string outputFile_;
0044
0045   struct RPPlots {
0046     std::unique_ptr<TH1D> h_de_x, h_de_y;
0047     RPPlots()
0048         : h_de_x(new TH1D("", ";#Deltax   (mm)", 1000, -1., +1.)),
0049           h_de_y(new TH1D("", ";#Deltay   (mm)", 1000, -1., +1.)) {}
0050
0051     void fill(double de_x, double de_y) {
0052       h_de_x->Fill(de_x);
0053       h_de_y->Fill(de_y);
0054     }
0055
0056     void write() const {
0057       h_de_x->Write("h_de_x");
0058       h_de_y->Write("h_de_y");
0059     }
0060   };
0061
0062   std::map<unsigned int, RPPlots> rp_plots_;
0063 };
0064
0065 //----------------------------------------------------------------------------------------------------
0066
0067 using namespace std;
0068 using namespace edm;
0069
0070 //----------------------------------------------------------------------------------------------------
0071
0072 CTPPSProtonReconstructionValidator::CTPPSProtonReconstructionValidator(const edm::ParameterSet& iConfig)
0073     : tokenRecoProtons_(consumes<reco::ForwardProtonCollection>(iConfig.getParameter<edm::InputTag>("tagRecoProtons"))),
0074       opticsESToken_(esConsumes()),
0075       chiSqCut_(iConfig.getParameter<double>("chiSqCut")),
0076       outputFile_(iConfig.getParameter<string>("outputFile")) {}
0077
0078 //----------------------------------------------------------------------------------------------------
0079
0080 void CTPPSProtonReconstructionValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0081   // get conditions
0082   const auto& opticalFunctions = iSetup.getData(opticsESToken_);
0083
0084   // stop if conditions invalid
0085   if (opticalFunctions.empty())
0086     return;
0087
0088   // get input
0089   Handle<reco::ForwardProtonCollection> hRecoProtons;
0090   iEvent.getByToken(tokenRecoProtons_, hRecoProtons);
0091
0092   // process tracks
0093   for (const auto& pr : *hRecoProtons) {
0094     if (!pr.validFit())
0095       continue;
0096
0097     if (pr.chi2() > chiSqCut_)
0098       continue;
0099
0100     for (const auto& tr : pr.contributingLocalTracks()) {
0101       CTPPSDetId rpId(tr->rpId());
0102       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0103
0104       // skip other than tracking RPs
0105       if (rpId.subdetId() != CTPPSDetId::sdTrackingStrip && rpId.subdetId() != CTPPSDetId::sdTrackingPixel)
0106         continue;
0107
0108       // try to get optics for the RP
0109       if (opticalFunctions.count(rpId) == 0)
0110         continue;
0111       const auto& func = opticalFunctions.at(rpId);
0112
0113       // do propagation
0114       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
0115       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam;
0116       func.transport(k_in_beam, k_out_beam);
0117
0118       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
0119           -pr.vx(), -pr.thetaX(), pr.vy(), pr.thetaY(), pr.xi()};  // conversions: CMS --> LHC convention
0120       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
0121       func.transport(k_in, k_out);
0122
0123       // fill plots
0124       const double de_x = (k_out.x - k_out_beam.x) * 10. - tr->x();  // conversions: cm --> mm
0125       const double de_y = (k_out.y - k_out_beam.y) * 10. - tr->y();  // conversions: cm --> mm
0126
0127       rp_plots_[rpDecId].fill(de_x, de_y);
0128     }
0129   }
0130 }
0131
0132 //----------------------------------------------------------------------------------------------------
0133
0134 void CTPPSProtonReconstructionValidator::endJob() {
0135   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0136
0137   for (const auto& p : rp_plots_) {
0138     gDirectory = f_out->mkdir(Form("%u", p.first));
0139     p.second.write();
0140   }
0141 }
0142
0143 //----------------------------------------------------------------------------------------------------
0144
0145 DEFINE_FWK_MODULE(CTPPSProtonReconstructionValidator);