Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-27 02:53:21

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