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;
  bool m_vertexConstraint;
  edm::InputTag m_beamSpot;
 
  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")) {
  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);
 
  edm::Handle<reco::BeamSpot> beamSpot;
  iEvent.getByLabel(m_beamSpot, beamSpot);
 
  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");