CTPPSProtonReconstructionValidator.cc

Go to the documentation of this file.

/****************************************************************************
 * Authors:
 *   Jan Kašpar
 ****************************************************************************/

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"

#include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
#include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"

#include "DataFormats/ProtonReco/interface/ForwardProton.h"
#include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
#include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"

#include "TFile.h"
#include "TH1D.h"

#include <map>
#include <string>

//----------------------------------------------------------------------------------------------------

class CTPPSProtonReconstructionValidator : public edm::one::EDAnalyzer<> {
public:
  explicit CTPPSProtonReconstructionValidator(const edm::ParameterSet&);

private:
  void analyze(const edm::Event&, const edm::EventSetup&) override;
  void endJob() override;

  edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtons_;
  edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsESToken_;
  double chiSqCut_;
  std::string outputFile_;

  struct RPPlots {
    std::unique_ptr<TH1D> h_de_x, h_de_y;
    RPPlots()
        : h_de_x(new TH1D("", ";#Deltax   (mm)", 1000, -1., +1.)),
          h_de_y(new TH1D("", ";#Deltay   (mm)", 1000, -1., +1.)) {}

    void fill(double de_x, double de_y) {
      h_de_x->Fill(de_x);
      h_de_y->Fill(de_y);
    }

    void write() const {
      h_de_x->Write("h_de_x");
      h_de_y->Write("h_de_y");
    }
  };

  std::map<unsigned int, RPPlots> rp_plots_;
};

//----------------------------------------------------------------------------------------------------

using namespace std;
using namespace edm;

//----------------------------------------------------------------------------------------------------

CTPPSProtonReconstructionValidator::CTPPSProtonReconstructionValidator(const edm::ParameterSet& iConfig)
    : tokenRecoProtons_(consumes<reco::ForwardProtonCollection>(iConfig.getParameter<edm::InputTag>("tagRecoProtons"))),
      opticsESToken_(esConsumes()),
      chiSqCut_(iConfig.getParameter<double>("chiSqCut")),
      outputFile_(iConfig.getParameter<string>("outputFile")) {}

//----------------------------------------------------------------------------------------------------

void CTPPSProtonReconstructionValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // get conditions
  const auto& opticalFunctions = iSetup.getData(opticsESToken_);

  // stop if conditions invalid
  if (opticalFunctions.empty())
    return;

  // get input
  Handle<reco::ForwardProtonCollection> hRecoProtons;
  iEvent.getByToken(tokenRecoProtons_, hRecoProtons);

  // process tracks
  for (const auto& pr : *hRecoProtons) {
    if (!pr.validFit())
      continue;

    if (pr.chi2() > chiSqCut_)
      continue;

    for (const auto& tr : pr.contributingLocalTracks()) {
      CTPPSDetId rpId(tr->rpId());
      unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();

      // skip other than tracking RPs
      if (rpId.subdetId() != CTPPSDetId::sdTrackingStrip && rpId.subdetId() != CTPPSDetId::sdTrackingPixel)
        continue;

      // try to get optics for the RP
      if (opticalFunctions.count(rpId) == 0)
        continue;
      const auto& func = opticalFunctions.at(rpId);

      // do propagation
      LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
      LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam;
      func.transport(k_in_beam, k_out_beam);

      LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
          -pr.vx(), -pr.thetaX(), pr.vy(), pr.thetaY(), pr.xi()};  // conversions: CMS --> LHC convention
      LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
      func.transport(k_in, k_out);

      // fill plots
      const double de_x = (k_out.x - k_out_beam.x) * 10. - tr->x();  // conversions: cm --> mm
      const double de_y = (k_out.y - k_out_beam.y) * 10. - tr->y();  // conversions: cm --> mm

      rp_plots_[rpDecId].fill(de_x, de_y);
    }
  }
}

//----------------------------------------------------------------------------------------------------

void CTPPSProtonReconstructionValidator::endJob() {
  auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");

  for (const auto& p : rp_plots_) {
    gDirectory = f_out->mkdir(Form("%u", p.first));
    p.second.write();
  }
}

//----------------------------------------------------------------------------------------------------

DEFINE_FWK_MODULE(CTPPSProtonReconstructionValidator);