SiPixelPhase1RecHitsV.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     SiPixelPhase1RecHitsV
// Class:       SiPixelPhase1RecHitsV
//

// Original Author: Marcel Schneider
// Additional Authors: Alexander Morton - modifying code for validation use

#include "Validation/SiPixelPhase1RecHitsV/interface/SiPixelPhase1RecHitsV.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"

#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

SiPixelPhase1RecHitsV::SiPixelPhase1RecHitsV(const edm::ParameterSet& iConfig)
    : SiPixelPhase1Base(iConfig),
      trackerHitAssociatorConfig_(iConfig, consumesCollector()),
      srcToken_(consumes<SiPixelRecHitCollection>(iConfig.getParameter<edm::InputTag>("src"))) {}

void SiPixelPhase1RecHitsV::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::Handle<SiPixelRecHitCollection> input;
  iEvent.getByToken(srcToken_, input);
  if (!input.isValid())
    return;

  TrackerHitAssociator associate(iEvent, trackerHitAssociatorConfig_);

  SiPixelRecHitCollection::const_iterator it;
  for (it = input->begin(); it != input->end(); ++it) {
    auto id = DetId(it->detId());

    for (SiPixelRecHit const& rechit : *it) {
      std::vector<PSimHit> associateSimHit;
      associateSimHit = associate.associateHit(rechit);
      std::vector<PSimHit>::const_iterator closestIt = associateSimHit.begin();

      LocalPoint lp = rechit.localPosition();
      float rechit_x = lp.x();
      float rechit_y = lp.y();

      LocalError lerr = rechit.localPositionError();
      float lerr_x = sqrt(lerr.xx());
      float lerr_y = sqrt(lerr.yy());

      // loop over associated sim hits and find the closest
      if (!associateSimHit.empty()) {
        float closestSimHit = 9999.9;

        for (std::vector<PSimHit>::const_iterator m = associateSimHit.begin(); m < associateSimHit.end(); m++) {
          float sim_x1((*m).entryPoint().x()), sim_x2((*m).exitPoint().x()), sim_xpos(0.5 * (sim_x1 + sim_x2));
          float sim_y1((*m).entryPoint().y()), sim_y2((*m).exitPoint().y()), sim_ypos(0.5 * (sim_y1 + sim_y2));

          float xres(std::abs(sim_xpos - rechit_x)), yres(std::abs(sim_ypos - rechit_y));
          float dist = std::sqrt(xres * xres + yres * yres);

          if (dist < closestSimHit) {
            closestSimHit = dist;
            closestIt = m;
          }
        }
      }

      // Sim Hit stuff

      if (!associateSimHit.empty()) {
        const PSimHit& simHit = *closestIt;

        int bunch = simHit.eventId().bunchCrossing();
        int event = simHit.eventId().event();

        float sim_x1(simHit.entryPoint().x()), sim_x2(simHit.exitPoint().x()), sim_xpos(0.5 * (sim_x1 + sim_x2));
        float sim_y1(simHit.entryPoint().y()), sim_y2(simHit.exitPoint().y()), sim_ypos(0.5 * (sim_y1 + sim_y2));

        float res_x = (rechit_x - sim_xpos) * 10000.0;
        float res_y = (rechit_y - sim_ypos) * 10000.0;

        float pull_x = (rechit_x - sim_xpos) / lerr_x;
        float pull_y = (rechit_y - sim_ypos) / lerr_y;

        // Now Plotting stuff
        if (bunch == 0)
          histo[IN_TIME_BUNCH].fill(bunch, id, &iEvent);
        if (bunch != 0)
          histo[OUT_TIME_BUNCH].fill(bunch, id, &iEvent);

        histo[NSIMHITS].fill(event, id, &iEvent);

        histo[RECHIT_X].fill(rechit_x, id, &iEvent);
        histo[RECHIT_Y].fill(rechit_y, id, &iEvent);

        histo[RES_X].fill(res_x, id, &iEvent);
        histo[RES_Y].fill(res_y, id, &iEvent);

        histo[ERROR_X].fill(lerr_x, id, &iEvent);
        histo[ERROR_Y].fill(lerr_y, id, &iEvent);

        histo[PULL_X].fill(pull_x, id, &iEvent);
        histo[PULL_Y].fill(pull_y, id, &iEvent);
      }
    }
  }
}

DEFINE_FWK_MODULE(SiPixelPhase1RecHitsV);