GEMSimTrackMatch.cc

Go to the documentation of this file.

#include "Validation/MuonGEMHits/interface/GEMSimTrackMatch.h"
#include "Validation/MuonGEMHits/interface/GEMDetLabel.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/GEMDigi/interface/GEMDigiCollection.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include <TMath.h>
#include <TH1F.h>

using namespace std;
using namespace GEMDetLabel;
GEMSimTrackMatch::GEMSimTrackMatch(const edm::ParameterSet& ps) : GEMTrackMatch(ps) {
  std::string simInputLabel_ = ps.getUntrackedParameter<std::string>("simInputLabel");

  simHitsToken_ = consumes<edm::PSimHitContainer>(edm::InputTag(simInputLabel_, "MuonGEMHits"));
  simTracksToken_ = consumes<edm::SimTrackContainer>(ps.getParameter<edm::InputTag>("simTrackCollection"));
  simVerticesToken_ = consumes<edm::SimVertexContainer>(ps.getParameter<edm::InputTag>("simVertexCollection"));
  detailPlot_ = ps.getParameter<bool>("detailPlot");

  cfg_ = ps;
}

void GEMSimTrackMatch::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& run, edm::EventSetup const& iSetup) {
  // Mandantory
  edm::ESHandle<GEMGeometry> hGeom;
  iSetup.get<MuonGeometryRecord>().get(hGeom);
  const GEMGeometry& geom = *hGeom;
  setGeometry(geom);

  ibooker.setCurrentFolder("MuonGEMHitsV/GEMHitsTask");
  LogDebug("GEMSimTrackMatch") << "ibooker set current folder\n";

  const float PI = TMath::Pi();

  nstation = geom.regions()[0]->stations().size();
  if (detailPlot_) {
    for (unsigned int j = 0; j < nstation; j++) {
      string track_eta_name = string("track_eta") + s_suffix[j];
      string track_eta_title = string("track_eta") + ";SimTrack |#eta|;# of tracks";
      track_eta[j] = ibooker.book1D(track_eta_name.c_str(), track_eta_title.c_str(), 140, minEta_, maxEta_);

      for (unsigned int k = 0; k < 3; k++) {
        string suffix = string(s_suffix[j]) + string(c_suffix[k]);
        string track_phi_name = string("track_phi") + suffix;
        string track_phi_title = string("track_phi") + suffix + ";SimTrack #phi;# of tracks";
        track_phi[j][k] = ibooker.book1D(track_phi_name.c_str(), track_phi_title.c_str(), 200, -PI, PI);
      }

      for (unsigned int i = 0; i < 4; i++) {
        string suffix = string(s_suffix[j]) + l_suffix[i];

        string sh_eta_name = string("sh_eta") + suffix;
        string sh_eta_title = sh_eta_name + "; tracks |#eta|; # of tracks";
        sh_eta[i][j] = ibooker.book1D(sh_eta_name.c_str(), sh_eta_title.c_str(), 140, minEta_, maxEta_);

        for (unsigned int k = 0; k < 3; k++) {
          suffix = string(s_suffix[j]) + string(l_suffix[i]) + string(c_suffix[k]);
          string sh_phi_name = string("sh_phi") + suffix;
          string sh_phi_title = sh_phi_name + "; tracks #phi; # of tracks";
          sh_phi[i][j][k] = ibooker.book1D((sh_phi_name).c_str(), sh_phi_title.c_str(), 200, -PI, PI);
        }
      }
    }
  }
}

GEMSimTrackMatch::~GEMSimTrackMatch() {}

void GEMSimTrackMatch::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::ESHandle<GEMGeometry> hGeom;
  iSetup.get<MuonGeometryRecord>().get(hGeom);
  const GEMGeometry& geom = *hGeom;

  MySimTrack track_;
  edm::Handle<edm::PSimHitContainer> simhits;
  edm::Handle<edm::SimTrackContainer> sim_tracks;
  edm::Handle<edm::SimVertexContainer> sim_vertices;

  iEvent.getByToken(simHitsToken_, simhits);
  iEvent.getByToken(simTracksToken_, sim_tracks);
  iEvent.getByToken(simVerticesToken_, sim_vertices);
  if (!simhits.isValid() || !sim_tracks.isValid() || !sim_vertices.isValid())
    return;

  //const edm::SimVertexContainer & sim_vert = *sim_vertices.product();
  const edm::SimTrackContainer& sim_trks = *sim_tracks.product();

  if (detailPlot_) {
    for (auto& t : sim_trks) {
      if (!isSimTrackGood(t)) {
        continue;
      }

      track_.pt = t.momentum().pt();
      track_.phi = t.momentum().phi();
      track_.eta = t.momentum().eta();
      std::fill(std::begin(track_.hitOdd), std::end(track_.hitOdd), false);
      std::fill(std::begin(track_.hitEven), std::end(track_.hitEven), false);
      for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 2; j++) {
          track_.gem_sh[i][j] = false;
        }
      }

      // match hits to this SimTrack
      const SimHitMatcher match_sh =
          SimHitMatcher(t, iEvent, geom, cfg_, simHitsToken_, simTracksToken_, simVerticesToken_);

      // check for hit chambers
      const auto gem_sh_ids_ch = match_sh.chamberIdsGEM();

      for (auto d : gem_sh_ids_ch) {
        const GEMDetId id(d);
        if (id.chamber() % 2 == 0)
          track_.hitEven[id.station() - 1] = true;
        else if (id.chamber() % 2 == 1)
          track_.hitOdd[id.station() - 1] = true;
        else {
          std::cout << "Error to get chamber id" << std::endl;
        }
        track_.gem_sh[id.station() - 1][(id.layer() - 1)] = true;
      }
      FillWithTrigger(track_eta, fabs(track_.eta));
      FillWithTrigger(track_phi, fabs(track_.eta), track_.phi, track_.hitOdd, track_.hitEven);
      FillWithTrigger(sh_eta, track_.gem_sh, fabs(track_.eta));
      FillWithTrigger(sh_phi, track_.gem_sh, fabs(track_.eta), track_.phi, track_.hitOdd, track_.hitEven);
    }
  }
}