AlignmentMonitorTracksFromTrajectories.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     CommonAlignmentProducer
// Class  :     AlignmentMonitorTracksFromTrajectories
//
// Implementation:
//     <Notes on implementation>
//
// Original Author:  Jim Pivarski
//         Created:  Thu Jun 28 01:38:33 CDT 2007
//

// system include files
#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "TH1.h"
#include "TObject.h"
#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"

#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
#include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

#include <fstream>

// user include files

//
// class definition
//

class AlignmentMonitorTracksFromTrajectories : public AlignmentMonitorBase {
public:
  AlignmentMonitorTracksFromTrajectories(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
  ~AlignmentMonitorTracksFromTrajectories(){};

  void book();
  void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection &iTrajTracks);
  void afterAlignment();

private:
  MuonServiceProxy *theMuonServiceProxy;
  MuonUpdatorAtVertex *theMuonUpdatorAtVertex;
  const bool m_vertexConstraint;
  const edm::InputTag m_beamSpot;
  const edm::EDGetTokenT<reco::BeamSpot> bsToken_;

  TH1F *m_diMuon_Z;
  TH1F *m_diMuon_Zforward;
  TH1F *m_diMuon_Zbarrel;
  TH1F *m_diMuon_Zbackward;
  TH1F *m_diMuon_Ups;
  TH1F *m_diMuon_Jpsi;
  TH1F *m_diMuon_log;
  TH1F *m_chi2_100;
  TH1F *m_chi2_10000;
  TH1F *m_chi2_1000000;
  TH1F *m_chi2_log;
  TH1F *m_chi2DOF_5;
  TH1F *m_chi2DOF_100;
  TH1F *m_chi2DOF_1000;
  TH1F *m_chi2DOF_log;
  TH1F *m_chi2_improvement;
  TH1F *m_chi2DOF_improvement;
  TH1F *m_pt[36];
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//

// AlignmentMonitorTracksFromTrajectories::AlignmentMonitorTracksFromTrajectories(const AlignmentMonitorTracksFromTrajectories& rhs)
// {
//    // do actual copying here;
// }

//
// assignment operators
//
// const AlignmentMonitorTracksFromTrajectories& AlignmentMonitorTracksFromTrajectories::operator=(const AlignmentMonitorTracksFromTrajectories& rhs)
// {
//   //An exception safe implementation is
//   AlignmentMonitorTracksFromTrajectories temp(rhs);
//   swap(rhs);
//
//   return *this;
// }

AlignmentMonitorTracksFromTrajectories::AlignmentMonitorTracksFromTrajectories(const edm::ParameterSet &cfg,
                                                                               edm::ConsumesCollector iC)
    : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorTracksFromTrajectories"),
      m_vertexConstraint(cfg.getParameter<bool>("vertexConstraint")),
      m_beamSpot(cfg.getParameter<edm::InputTag>("beamSpot")),
      bsToken_(iC.consumes<reco::BeamSpot>(m_beamSpot)) {
  theMuonServiceProxy = new MuonServiceProxy(cfg.getParameter<edm::ParameterSet>("ServiceParameters"));
  theMuonUpdatorAtVertex = new MuonUpdatorAtVertex(cfg.getParameter<edm::ParameterSet>("MuonUpdatorAtVertexParameters"),
                                                   theMuonServiceProxy);
}

//
// member functions
//

// book()

void AlignmentMonitorTracksFromTrajectories::book() {
  m_diMuon_Z = book1D("/iterN/", "diMuon_Z", "Di-muon mass (GeV)", 200, 90. - 50., 90. + 50.);
  m_diMuon_Zforward = book1D("/iterN/", "diMuon_Zforward", "Di-muon mass (GeV) eta > 1.4", 200, 90. - 50., 90. + 50.);
  m_diMuon_Zbarrel =
      book1D("/iterN/", "diMuon_Zbarrel", "Di-muon mass (GeV) -1.4 < eta < 1.4", 200, 90. - 50., 90. + 50.);
  m_diMuon_Zbackward =
      book1D("/iterN/", "diMuon_Zbackward", "Di-muon mass (GeV) eta < -1.4", 200, 90. - 50., 90. + 50.);
  m_diMuon_Ups = book1D("/iterN/", "diMuon_Ups", "Di-muon mass (GeV)", 200, 9. - 3., 9. + 3.);
  m_diMuon_Jpsi = book1D("/iterN/", "diMuon_Jpsi", "Di-muon mass (GeV)", 200, 3. - 1., 3. + 1.);
  m_diMuon_log = book1D("/iterN/", "diMuon_log", "Di-muon mass (log GeV)", 300, 0, 3);
  m_chi2_100 = book1D("/iterN/", "m_chi2_100", "Track chi^2", 100, 0, 100);
  m_chi2_10000 = book1D("/iterN/", "m_chi2_10000", "Track chi^2", 100, 0, 10000);
  m_chi2_1000000 = book1D("/iterN/", "m_chi2_1000000", "Track chi^2", 100, 1, 1000000);
  m_chi2_log = book1D("/iterN/", "m_chi2_log", "Log track chi^2", 100, -3, 7);
  m_chi2DOF_5 = book1D("/iterN/", "m_chi2DOF_5", "Track chi^2/nDOF", 100, 0, 5);
  m_chi2DOF_100 = book1D("/iterN/", "m_chi2DOF_100", "Track chi^2/nDOF", 100, 0, 100);
  m_chi2DOF_1000 = book1D("/iterN/", "m_chi2DOF_1000", "Track chi^2/nDOF", 100, 0, 1000);
  m_chi2DOF_log = book1D("/iterN/", "m_chi2DOF_log", "Log track chi^2/nDOF", 100, -3, 7);
  m_chi2_improvement = book1D("/iterN/", "m_chi2_improvement", "Track-by-track chi^2/original chi^2", 100, 0., 10.);
  m_chi2DOF_improvement = book1D(
      "/iterN/", "m_chi2DOF_improvement", "Track-by-track (chi^2/DOF)/(original chi^2/original DOF)", 100, 0., 10.);
  for (int i = 0; i < 36; i++) {
    char name[100], title[100];
    snprintf(name, sizeof(name), "m_pt_phi%d", i);
    snprintf(title, sizeof(title), "Track pt (GeV) in phi bin %d/36", i);
    m_pt[i] = book1D("/iterN/", name, title, 100, 0, 100);
  }
}

// event()

void AlignmentMonitorTracksFromTrajectories::event(const edm::Event &iEvent,
                                                   const edm::EventSetup &iSetup,
                                                   const ConstTrajTrackPairCollection &tracks) {
  theMuonServiceProxy->update(iSetup);

  const edm::Handle<reco::BeamSpot> &beamSpot = iEvent.getHandle(bsToken_);

  GlobalVector p1, p2;
  double e1 = 0.;
  double e2 = 0.;

  for (ConstTrajTrackPairCollection::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
    const Trajectory *traj = it->first;
    const reco::Track *track = it->second;

    int nDOF = 0;
    TrajectoryStateOnSurface closestTSOS;
    double closest = 10000.;

    std::vector<TrajectoryMeasurement> measurements = traj->measurements();
    for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end(); ++im) {
      const TrajectoryMeasurement meas = *im;
      const TransientTrackingRecHit *hit = &(*meas.recHit());
      //     const DetId id = hit->geographicalId();

      nDOF += hit->dimension();

      TrajectoryStateOnSurface TSOS = meas.backwardPredictedState();
      GlobalPoint where = TSOS.surface().toGlobal(LocalPoint(0, 0, 0));

      if (where.mag() < closest) {
        closest = where.mag();
        closestTSOS = TSOS;
      }

    }  // end loop over measurements
    nDOF -= 5;

    if (closest != 10000.) {
      std::pair<bool, FreeTrajectoryState> state;

      if (m_vertexConstraint) {
        state = theMuonUpdatorAtVertex->propagateWithUpdate(closestTSOS, *beamSpot);
        // add in chi^2 contribution from vertex contratint?
      } else {
        state = theMuonUpdatorAtVertex->propagate(closestTSOS, *beamSpot);
      }

      if (state.first) {
        double chi2 = traj->chiSquared();
        double chi2DOF = chi2 / double(nDOF);

        m_chi2_100->Fill(chi2);
        m_chi2_10000->Fill(chi2);
        m_chi2_1000000->Fill(chi2);
        m_chi2_log->Fill(log10(chi2));

        m_chi2DOF_5->Fill(chi2DOF);
        m_chi2DOF_100->Fill(chi2DOF);
        m_chi2DOF_1000->Fill(chi2DOF);
        m_chi2DOF_log->Fill(log10(chi2DOF));
        m_chi2_improvement->Fill(chi2 / track->chi2());
        m_chi2DOF_improvement->Fill(chi2DOF / track->normalizedChi2());

        GlobalVector momentum = state.second.momentum();
        double energy = momentum.mag();

        if (energy > e1) {
          e2 = e1;
          e1 = energy;
          p2 = p1;
          p1 = momentum;
        } else if (energy > e2) {
          e2 = energy;
          p2 = momentum;
        }

        double pt = momentum.perp();
        double phi = momentum.phi();
        while (phi < -M_PI)
          phi += 2. * M_PI;
        while (phi > M_PI)
          phi -= 2. * M_PI;

        double phibin = (phi + M_PI) / (2. * M_PI) * 36.;

        for (int i = 0; i < 36; i++) {
          if (phibin < i + 1) {
            m_pt[i]->Fill(pt);
            break;
          }
        }

      }  // end if propagate was successful
    }    // end if we have a closest TSOS
  }      // end loop over tracks

  if (e1 > 0. && e2 > 0.) {
    double energy_tot = e1 + e2;
    GlobalVector momentum_tot = p1 + p2;
    double mass = sqrt(energy_tot * energy_tot - momentum_tot.mag2());
    double eta = momentum_tot.eta();

    m_diMuon_Z->Fill(mass);
    if (eta > 1.4)
      m_diMuon_Zforward->Fill(mass);
    else if (eta > -1.4)
      m_diMuon_Zbarrel->Fill(mass);
    else
      m_diMuon_Zbackward->Fill(mass);

    m_diMuon_Ups->Fill(mass);
    m_diMuon_Jpsi->Fill(mass);
    m_diMuon_log->Fill(log10(mass));
  }  // end if we have two tracks
}

// afterAlignment()

void AlignmentMonitorTracksFromTrajectories::afterAlignment() {}

//
// const member functions
//

//
// static member functions
//

//
// SEAL definitions
//

DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory,
                  AlignmentMonitorTracksFromTrajectories,
                  "AlignmentMonitorTracksFromTrajectories");