Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-03 01:57:38

0001 // -*- C++ -*-
0002 //
0003 // Package:     CommonAlignmentProducer
0004 // Class  :     AlignmentMonitorTracksFromTrajectories
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:  Jim Pivarski
0010 //         Created:  Thu Jun 28 01:38:33 CDT 2007
0011 //
0012 
0013 // system include files
0014 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0018 #include "DataFormats/GeometrySurface/interface/Surface.h"
0019 #include "TH1.h"
0020 #include "TObject.h"
0021 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
0022 
0023 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0024 #include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
0025 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0026 
0027 #include <fstream>
0028 
0029 // user include files
0030 
0031 //
0032 // class definition
0033 //
0034 
0035 class AlignmentMonitorTracksFromTrajectories : public AlignmentMonitorBase {
0036 public:
0037   AlignmentMonitorTracksFromTrajectories(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
0038   ~AlignmentMonitorTracksFromTrajectories(){};
0039 
0040   void book();
0041   void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection &iTrajTracks);
0042   void afterAlignment();
0043 
0044 private:
0045   MuonServiceProxy *theMuonServiceProxy;
0046   MuonUpdatorAtVertex *theMuonUpdatorAtVertex;
0047   bool m_vertexConstraint;
0048   edm::InputTag m_beamSpot;
0049 
0050   TH1F *m_diMuon_Z;
0051   TH1F *m_diMuon_Zforward;
0052   TH1F *m_diMuon_Zbarrel;
0053   TH1F *m_diMuon_Zbackward;
0054   TH1F *m_diMuon_Ups;
0055   TH1F *m_diMuon_Jpsi;
0056   TH1F *m_diMuon_log;
0057   TH1F *m_chi2_100;
0058   TH1F *m_chi2_10000;
0059   TH1F *m_chi2_1000000;
0060   TH1F *m_chi2_log;
0061   TH1F *m_chi2DOF_5;
0062   TH1F *m_chi2DOF_100;
0063   TH1F *m_chi2DOF_1000;
0064   TH1F *m_chi2DOF_log;
0065   TH1F *m_chi2_improvement;
0066   TH1F *m_chi2DOF_improvement;
0067   TH1F *m_pt[36];
0068 };
0069 
0070 //
0071 // constants, enums and typedefs
0072 //
0073 
0074 //
0075 // static data member definitions
0076 //
0077 
0078 //
0079 // constructors and destructor
0080 //
0081 
0082 // AlignmentMonitorTracksFromTrajectories::AlignmentMonitorTracksFromTrajectories(const AlignmentMonitorTracksFromTrajectories& rhs)
0083 // {
0084 //    // do actual copying here;
0085 // }
0086 
0087 //
0088 // assignment operators
0089 //
0090 // const AlignmentMonitorTracksFromTrajectories& AlignmentMonitorTracksFromTrajectories::operator=(const AlignmentMonitorTracksFromTrajectories& rhs)
0091 // {
0092 //   //An exception safe implementation is
0093 //   AlignmentMonitorTracksFromTrajectories temp(rhs);
0094 //   swap(rhs);
0095 //
0096 //   return *this;
0097 // }
0098 
0099 AlignmentMonitorTracksFromTrajectories::AlignmentMonitorTracksFromTrajectories(const edm::ParameterSet &cfg,
0100                                                                                edm::ConsumesCollector iC)
0101     : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorTracksFromTrajectories"),
0102       m_vertexConstraint(cfg.getParameter<bool>("vertexConstraint")),
0103       m_beamSpot(cfg.getParameter<edm::InputTag>("beamSpot")) {
0104   theMuonServiceProxy = new MuonServiceProxy(cfg.getParameter<edm::ParameterSet>("ServiceParameters"));
0105   theMuonUpdatorAtVertex = new MuonUpdatorAtVertex(cfg.getParameter<edm::ParameterSet>("MuonUpdatorAtVertexParameters"),
0106                                                    theMuonServiceProxy);
0107 }
0108 
0109 //
0110 // member functions
0111 //
0112 
0113 //////////////////////////////////////////////////////////////////////
0114 // book()
0115 //////////////////////////////////////////////////////////////////////
0116 
0117 void AlignmentMonitorTracksFromTrajectories::book() {
0118   m_diMuon_Z = book1D("/iterN/", "diMuon_Z", "Di-muon mass (GeV)", 200, 90. - 50., 90. + 50.);
0119   m_diMuon_Zforward = book1D("/iterN/", "diMuon_Zforward", "Di-muon mass (GeV) eta > 1.4", 200, 90. - 50., 90. + 50.);
0120   m_diMuon_Zbarrel =
0121       book1D("/iterN/", "diMuon_Zbarrel", "Di-muon mass (GeV) -1.4 < eta < 1.4", 200, 90. - 50., 90. + 50.);
0122   m_diMuon_Zbackward =
0123       book1D("/iterN/", "diMuon_Zbackward", "Di-muon mass (GeV) eta < -1.4", 200, 90. - 50., 90. + 50.);
0124   m_diMuon_Ups = book1D("/iterN/", "diMuon_Ups", "Di-muon mass (GeV)", 200, 9. - 3., 9. + 3.);
0125   m_diMuon_Jpsi = book1D("/iterN/", "diMuon_Jpsi", "Di-muon mass (GeV)", 200, 3. - 1., 3. + 1.);
0126   m_diMuon_log = book1D("/iterN/", "diMuon_log", "Di-muon mass (log GeV)", 300, 0, 3);
0127   m_chi2_100 = book1D("/iterN/", "m_chi2_100", "Track chi^2", 100, 0, 100);
0128   m_chi2_10000 = book1D("/iterN/", "m_chi2_10000", "Track chi^2", 100, 0, 10000);
0129   m_chi2_1000000 = book1D("/iterN/", "m_chi2_1000000", "Track chi^2", 100, 1, 1000000);
0130   m_chi2_log = book1D("/iterN/", "m_chi2_log", "Log track chi^2", 100, -3, 7);
0131   m_chi2DOF_5 = book1D("/iterN/", "m_chi2DOF_5", "Track chi^2/nDOF", 100, 0, 5);
0132   m_chi2DOF_100 = book1D("/iterN/", "m_chi2DOF_100", "Track chi^2/nDOF", 100, 0, 100);
0133   m_chi2DOF_1000 = book1D("/iterN/", "m_chi2DOF_1000", "Track chi^2/nDOF", 100, 0, 1000);
0134   m_chi2DOF_log = book1D("/iterN/", "m_chi2DOF_log", "Log track chi^2/nDOF", 100, -3, 7);
0135   m_chi2_improvement = book1D("/iterN/", "m_chi2_improvement", "Track-by-track chi^2/original chi^2", 100, 0., 10.);
0136   m_chi2DOF_improvement = book1D(
0137       "/iterN/", "m_chi2DOF_improvement", "Track-by-track (chi^2/DOF)/(original chi^2/original DOF)", 100, 0., 10.);
0138   for (int i = 0; i < 36; i++) {
0139     char name[100], title[100];
0140     snprintf(name, sizeof(name), "m_pt_phi%d", i);
0141     snprintf(title, sizeof(title), "Track pt (GeV) in phi bin %d/36", i);
0142     m_pt[i] = book1D("/iterN/", name, title, 100, 0, 100);
0143   }
0144 }
0145 
0146 //////////////////////////////////////////////////////////////////////
0147 // event()
0148 //////////////////////////////////////////////////////////////////////
0149 
0150 void AlignmentMonitorTracksFromTrajectories::event(const edm::Event &iEvent,
0151                                                    const edm::EventSetup &iSetup,
0152                                                    const ConstTrajTrackPairCollection &tracks) {
0153   theMuonServiceProxy->update(iSetup);
0154 
0155   edm::Handle<reco::BeamSpot> beamSpot;
0156   iEvent.getByLabel(m_beamSpot, beamSpot);
0157 
0158   GlobalVector p1, p2;
0159   double e1 = 0.;
0160   double e2 = 0.;
0161 
0162   for (ConstTrajTrackPairCollection::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0163     const Trajectory *traj = it->first;
0164     const reco::Track *track = it->second;
0165 
0166     int nDOF = 0;
0167     TrajectoryStateOnSurface closestTSOS;
0168     double closest = 10000.;
0169 
0170     std::vector<TrajectoryMeasurement> measurements = traj->measurements();
0171     for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end(); ++im) {
0172       const TrajectoryMeasurement meas = *im;
0173       const TransientTrackingRecHit *hit = &(*meas.recHit());
0174       //     const DetId id = hit->geographicalId();
0175 
0176       nDOF += hit->dimension();
0177 
0178       TrajectoryStateOnSurface TSOS = meas.backwardPredictedState();
0179       GlobalPoint where = TSOS.surface().toGlobal(LocalPoint(0, 0, 0));
0180 
0181       if (where.mag() < closest) {
0182         closest = where.mag();
0183         closestTSOS = TSOS;
0184       }
0185 
0186     }  // end loop over measurements
0187     nDOF -= 5;
0188 
0189     if (closest != 10000.) {
0190       std::pair<bool, FreeTrajectoryState> state;
0191 
0192       if (m_vertexConstraint) {
0193         state = theMuonUpdatorAtVertex->propagateWithUpdate(closestTSOS, *beamSpot);
0194         // add in chi^2 contribution from vertex contratint?
0195       } else {
0196         state = theMuonUpdatorAtVertex->propagate(closestTSOS, *beamSpot);
0197       }
0198 
0199       if (state.first) {
0200         double chi2 = traj->chiSquared();
0201         double chi2DOF = chi2 / double(nDOF);
0202 
0203         m_chi2_100->Fill(chi2);
0204         m_chi2_10000->Fill(chi2);
0205         m_chi2_1000000->Fill(chi2);
0206         m_chi2_log->Fill(log10(chi2));
0207 
0208         m_chi2DOF_5->Fill(chi2DOF);
0209         m_chi2DOF_100->Fill(chi2DOF);
0210         m_chi2DOF_1000->Fill(chi2DOF);
0211         m_chi2DOF_log->Fill(log10(chi2DOF));
0212         m_chi2_improvement->Fill(chi2 / track->chi2());
0213         m_chi2DOF_improvement->Fill(chi2DOF / track->normalizedChi2());
0214 
0215         GlobalVector momentum = state.second.momentum();
0216         double energy = momentum.mag();
0217 
0218         if (energy > e1) {
0219           e2 = e1;
0220           e1 = energy;
0221           p2 = p1;
0222           p1 = momentum;
0223         } else if (energy > e2) {
0224           e2 = energy;
0225           p2 = momentum;
0226         }
0227 
0228         double pt = momentum.perp();
0229         double phi = momentum.phi();
0230         while (phi < -M_PI)
0231           phi += 2. * M_PI;
0232         while (phi > M_PI)
0233           phi -= 2. * M_PI;
0234 
0235         double phibin = (phi + M_PI) / (2. * M_PI) * 36.;
0236 
0237         for (int i = 0; i < 36; i++) {
0238           if (phibin < i + 1) {
0239             m_pt[i]->Fill(pt);
0240             break;
0241           }
0242         }
0243 
0244       }  // end if propagate was successful
0245     }    // end if we have a closest TSOS
0246   }      // end loop over tracks
0247 
0248   if (e1 > 0. && e2 > 0.) {
0249     double energy_tot = e1 + e2;
0250     GlobalVector momentum_tot = p1 + p2;
0251     double mass = sqrt(energy_tot * energy_tot - momentum_tot.mag2());
0252     double eta = momentum_tot.eta();
0253 
0254     m_diMuon_Z->Fill(mass);
0255     if (eta > 1.4)
0256       m_diMuon_Zforward->Fill(mass);
0257     else if (eta > -1.4)
0258       m_diMuon_Zbarrel->Fill(mass);
0259     else
0260       m_diMuon_Zbackward->Fill(mass);
0261 
0262     m_diMuon_Ups->Fill(mass);
0263     m_diMuon_Jpsi->Fill(mass);
0264     m_diMuon_log->Fill(log10(mass));
0265   }  // end if we have two tracks
0266 }
0267 
0268 //////////////////////////////////////////////////////////////////////
0269 // afterAlignment()
0270 //////////////////////////////////////////////////////////////////////
0271 
0272 void AlignmentMonitorTracksFromTrajectories::afterAlignment() {}
0273 
0274 //
0275 // const member functions
0276 //
0277 
0278 //
0279 // static member functions
0280 //
0281 
0282 //
0283 // SEAL definitions
0284 //
0285 
0286 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory,
0287                   AlignmentMonitorTracksFromTrajectories,
0288                   "AlignmentMonitorTracksFromTrajectories");