Muon/src/MuonSeedsAnalyzer.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  $Date: 2008/03/25
0005  18:37:05 $
0006  *  \author G. Mila - INFN Torino
0007  */
0008
0009 #include "DQMOffline/Muon/interface/MuonSeedsAnalyzer.h"
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0016 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0017 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0018 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0019
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021
0022 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024
0025 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0026 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0027 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0028 #include "FWCore/Framework/interface/ConsumesCollector.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030
0031 #include <string>
0032
0033 using namespace std;
0034 using namespace edm;
0035
0036 MuonSeedsAnalyzer::MuonSeedsAnalyzer(const edm::ParameterSet& pSet) {
0037   parameters = pSet;
0038
0039   theService = new MuonServiceProxy(parameters.getParameter<ParameterSet>("ServiceParameters"), consumesCollector());
0040
0041   theSeedsCollectionLabel_ = consumes<TrajectorySeedCollection>(parameters.getParameter<InputTag>("SeedCollection"));
0042
0043   seedHitBin = parameters.getParameter<int>("RecHitBin");
0044   seedHitMin = parameters.getParameter<double>("RecHitMin");
0045   seedHitMax = parameters.getParameter<double>("RecHitMax");
0046   PhiBin = parameters.getParameter<int>("PhiBin");
0047   PhiMin = parameters.getParameter<double>("PhiMin");
0048   PhiMax = parameters.getParameter<double>("PhiMax");
0049   EtaBin = parameters.getParameter<int>("EtaBin");
0050   EtaMin = parameters.getParameter<double>("EtaMin");
0051   EtaMax = parameters.getParameter<double>("EtaMax");
0052   ThetaBin = parameters.getParameter<int>("ThetaBin");
0053   ThetaMin = parameters.getParameter<double>("ThetaMin");
0054   ThetaMax = parameters.getParameter<double>("ThetaMax");
0055   seedPtBin = parameters.getParameter<int>("seedPtBin");
0056   seedPtMin = parameters.getParameter<double>("seedPtMin");
0057   seedPtMax = parameters.getParameter<double>("seedPtMax");
0058   seedPxyzBin = parameters.getParameter<int>("seedPxyzBin");
0059   seedPxyzMin = parameters.getParameter<double>("seedPxyzMin");
0060   seedPxyzMax = parameters.getParameter<double>("seedPxyzMax");
0061   pErrBin = parameters.getParameter<int>("pErrBin");
0062   pErrMin = parameters.getParameter<double>("pErrMin");
0063   pErrMax = parameters.getParameter<double>("pErrMax");
0064   pxyzErrBin = parameters.getParameter<int>("pxyzErrBin");
0065   pxyzErrMin = parameters.getParameter<double>("pxyzErrMin");
0066   pxyzErrMax = parameters.getParameter<double>("pxyzErrMax");
0067   phiErrBin = parameters.getParameter<int>("phiErrBin");
0068   phiErrMin = parameters.getParameter<double>("phiErrMin");
0069   phiErrMax = parameters.getParameter<double>("phiErrMax");
0070   etaErrBin = parameters.getParameter<int>("etaErrBin");
0071   etaErrMin = parameters.getParameter<double>("etaErrMin");
0072   etaErrMax = parameters.getParameter<double>("etaErrMax");
0073 }
0074 MuonSeedsAnalyzer::~MuonSeedsAnalyzer() { delete theService; }
0075 void MuonSeedsAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0076   ibooker.cd();
0077   ibooker.setCurrentFolder("Muons/MuonSeedsAnalyzer");
0078
0079   string histname = "NumberOfRecHitsPerSeed_";
0080   NumberOfRecHitsPerSeed = ibooker.book1D(histname, "Number of seed recHits", seedHitBin, seedHitMin, seedHitMax);
0081
0082   histname = "seedPhi_";
0083   seedPhi = ibooker.book1D(histname, "Seed #phi", PhiBin, PhiMin, PhiMax);
0084   seedPhi->setAxisTitle("rad");
0085
0086   histname = "seedEta_";
0087   seedEta = ibooker.book1D(histname, "Seed #eta", EtaBin, EtaMin, EtaMax);
0088
0089   histname = "seedTheta_";
0090   seedTheta = ibooker.book1D(histname, "Seed #theta", ThetaBin, ThetaMin, ThetaMax);
0091   seedTheta->setAxisTitle("rad");
0092
0093   histname = "seedPt_";
0094   seedPt = ibooker.book1D(histname, "Seed p_{t}", seedPtBin, seedPtMin, seedPtMax);
0095   seedPt->setAxisTitle("GeV");
0096
0097   histname = "seedPx_";
0098   seedPx = ibooker.book1D(histname, "Seed p_{x}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
0099   seedPx->setAxisTitle("GeV");
0100   histname = "seedPy_";
0101   seedPy = ibooker.book1D(histname, "Seed p_{y}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
0102   seedPy->setAxisTitle("GeV");
0103   histname = "seedPz_";
0104   seedPz = ibooker.book1D(histname, "Seed p_{z}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
0105   seedPz->setAxisTitle("GeV");
0106
0107   histname = "seedPtErrOverPt_";
0108   seedPtErr = ibooker.book1D(histname, "Seed p_{t}Err/p_{t}", pErrBin, pErrMin, pErrMax);
0109   histname = "seedPtErrOverPtVsPhi_";
0110   seedPtErrVsPhi =
0111       ibooker.book2D(histname, "Seed p_{t}Err/p_{t} vs #phi", PhiBin, PhiMin, PhiMax, pErrBin, pErrMin, pErrMax);
0112   seedPtErrVsPhi->setAxisTitle("rad", 2);
0113   histname = "seedPtErrOverPtVsEta_";
0114   seedPtErrVsEta =
0115       ibooker.book2D(histname, "Seed p_{t}Err/p_{t} vs #eta", EtaBin, EtaMin, EtaMax, pErrBin, pErrMin, pErrMax);
0116   histname = "seedPtErrOverPtVsPt_";
0117   seedPtErrVsPt = ibooker.book2D(
0118       histname, "Seed p_{t}Err/p_{t} vs p_{t}", seedPtBin / 5, seedPtMin, seedPtMax, pErrBin, pErrMin, pErrMax);
0119   seedPtErrVsPt->setAxisTitle("GeV", 2);
0120   histname = "seedPErrOverP_";
0121   seedPErr = ibooker.book1D(histname, "Seed pErr/p", pErrBin, pErrMin, pErrMax);
0122   histname = "seedPErrOverPVsPhi_";
0123   seedPErrVsPhi = ibooker.book2D(histname, "Seed pErr/p vs #phi", PhiBin, PhiMin, PhiMax, pErrBin, pErrMin, pErrMax);
0124   seedPErrVsPhi->setAxisTitle("rad", 2);
0125   histname = "seedPErrOverPVsEta_";
0126   seedPErrVsEta = ibooker.book2D(histname, "Seed pErr/p vs #eta", EtaBin, EtaMin, EtaMax, pErrBin, pErrMin, pErrMax);
0127   histname = "seedPErrOverPVsPt_";
0128   seedPErrVsPt =
0129       ibooker.book2D(histname, "Seed pErr/p vs p_{t}", seedPtBin / 5, seedPtMin, seedPtMax, pErrBin, pErrMin, pErrMax);
0130   seedPErrVsPt->setAxisTitle("GeV", 2);
0131
0132   histname = "seedPxErrOverPx_";
0133   seedPxErr = ibooker.book1D(histname, "Seed p_{x}Err/p_{x}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
0134   histname = "seedPyErrOverPy_";
0135   seedPyErr = ibooker.book1D(histname, "Seed p_{y}Err/p_{y}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
0136   histname = "seedPzErrOverPz_";
0137   seedPzErr = ibooker.book1D(histname, "Seed p_{z}Err/p_{z}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
0138
0139   histname = "seedPhiErr_";
0140   seedPhiErr = ibooker.book1D(histname, "Seed #phi error", phiErrBin, phiErrMin, phiErrMax);
0141
0142   histname = "seedEtaErr_";
0143   seedEtaErr = ibooker.book1D(histname, "Seed #eta error", etaErrBin, etaErrMin, etaErrMax);
0144 }
0145
0146 void MuonSeedsAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0147   theService->update(iSetup);
0148
0149   // Take the seeds container
0150   edm::Handle<TrajectorySeedCollection> seeds;
0151   iEvent.getByToken(theSeedsCollectionLabel_, seeds);
0152
0153   // if not valid, skip
0154   if (!seeds.isValid())
0155     return;
0156
0157   for (TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed) {
0158     //    const TrajectorySeed sd = *seed;
0159
0160     // Get the Trajectory State on Det (persistent version of a TSOS) from the seed
0161     PTrajectoryStateOnDet pTSOD = seed->startingState();
0162
0163     // Transform it in a TrajectoryStateOnSurface
0164     DetId seedDetId(pTSOD.detId());
0165     const GeomDet* gdet = theService->trackingGeometry()->idToDet(seedDetId);
0166     TrajectoryStateOnSurface seedTSOS =
0167         trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*(theService)->magneticField());
0168     AlgebraicSymMatrix66 errors = seedTSOS.cartesianError().matrix();
0169     double partialPterror =
0170         errors(3, 3) * pow(seedTSOS.globalMomentum().x(), 2) + errors(4, 4) * pow(seedTSOS.globalMomentum().y(), 2);
0171
0172     LogTrace(metname) << "[MuonSeedAnalyzer] Filling the histos";
0173
0174     // nhits
0175     LogTrace(metname) << "Number od recHits per seed: " << seed->nHits();
0176     NumberOfRecHitsPerSeed->Fill(seed->nHits());
0177
0178     // pt
0179     LogTrace(metname) << "seed momentum: " << seedTSOS.globalMomentum().perp();
0180     seedPt->Fill(seedTSOS.globalMomentum().perp());
0181
0182     // px
0183     LogTrace(metname) << "seed px: " << seedTSOS.globalMomentum().x();
0184     seedPx->Fill(seedTSOS.globalMomentum().x());
0185
0186     // py
0187     LogTrace(metname) << "seed py: " << seedTSOS.globalMomentum().y();
0188     seedPy->Fill(seedTSOS.globalMomentum().y());
0189
0190     // pz
0191     LogTrace(metname) << "seed pz: " << seedTSOS.globalMomentum().z();
0192     seedPz->Fill(seedTSOS.globalMomentum().z());
0193
0194     // phi
0195     LogTrace(metname) << "seed phi: " << seedTSOS.globalMomentum().phi();
0196     seedPhi->Fill(seedTSOS.globalMomentum().phi());
0197
0198     // theta
0199     LogTrace(metname) << "seed theta: " << seedTSOS.globalMomentum().theta();
0200     seedTheta->Fill(seedTSOS.globalMomentum().theta());
0201
0202     // eta
0203     LogTrace(metname) << "seed eta: " << seedTSOS.globalMomentum().eta();
0204     seedEta->Fill(seedTSOS.globalMomentum().eta());
0205
0206     // pt err
0207     LogTrace(metname) << "seed pt error: " << sqrt(partialPterror) / seedTSOS.globalMomentum().perp();
0208     seedPtErr->Fill(sqrt(partialPterror) / seedTSOS.globalMomentum().perp());
0209
0210     // ptErr/pt Vs phi
0211     seedPtErrVsPhi->Fill(seedTSOS.globalMomentum().phi(), sqrt(partialPterror) / seedTSOS.globalMomentum().perp());
0212     // ptErr/pt Vs eta
0213     seedPtErrVsEta->Fill(seedTSOS.globalMomentum().eta(), sqrt(partialPterror) / seedTSOS.globalMomentum().perp());
0214     // ptErr/pt Vs pt
0215     seedPtErrVsPt->Fill(seedTSOS.globalMomentum().perp(), sqrt(partialPterror) / seedTSOS.globalMomentum().perp());
0216
0217     // px err
0218     LogTrace(metname) << "seed px error: " << sqrt(errors(3, 3)) / seedTSOS.globalMomentum().x();
0219     seedPxErr->Fill(sqrt(errors(3, 3)) / seedTSOS.globalMomentum().x());
0220
0221     // py err
0222     LogTrace(metname) << "seed py error: " << sqrt(errors(4, 4)) / seedTSOS.globalMomentum().y();
0223     seedPyErr->Fill(sqrt(errors(4, 4)) / seedTSOS.globalMomentum().y());
0224
0225     // pz err
0226     LogTrace(metname) << "seed pz error: " << sqrt(errors(5, 5)) / seedTSOS.globalMomentum().z();
0227     seedPzErr->Fill(sqrt(errors(5, 5)) / seedTSOS.globalMomentum().z());
0228
0229     // p err
0230     LogTrace(metname) << "seed p error: "
0231                       << sqrt(partialPterror + errors(5, 5) * pow(seedTSOS.globalMomentum().z(), 2)) /
0232                              seedTSOS.globalMomentum().mag();
0233     seedPErr->Fill(sqrt(partialPterror + errors(5, 5) * pow(seedTSOS.globalMomentum().z(), 2)) /
0234                    seedTSOS.globalMomentum().mag());
0235
0236     // pErr/p Vs phi
0237     seedPErrVsPhi->Fill(
0238         seedTSOS.globalMomentum().phi(),
0239         sqrt(partialPterror + errors(5, 5) * pow(seedTSOS.globalMomentum().z(), 2)) / seedTSOS.globalMomentum().mag());
0240     // pErr/p Vs eta
0241     seedPErrVsEta->Fill(
0242         seedTSOS.globalMomentum().eta(),
0243         sqrt(partialPterror + errors(5, 5) * pow(seedTSOS.globalMomentum().z(), 2)) / seedTSOS.globalMomentum().mag());
0244     // pErr/p Vs pt
0245     seedPErrVsPt->Fill(
0246         seedTSOS.globalMomentum().perp(),
0247         sqrt(partialPterror + errors(5, 5) * pow(seedTSOS.globalMomentum().z(), 2)) / seedTSOS.globalMomentum().mag());
0248
0249     // phi err
0250     LogTrace(metname) << "seed phi error: " << sqrt(seedTSOS.curvilinearError().matrix()(2, 2));
0251     seedPhiErr->Fill(sqrt(seedTSOS.curvilinearError().matrix()(2, 2)));
0252
0253     // eta err
0254     LogTrace(metname) << "seed eta error: "
0255                       << sqrt(seedTSOS.curvilinearError().matrix()(1, 1)) * abs(sin(seedTSOS.globalMomentum().theta()));
0256     seedEtaErr->Fill(sqrt(seedTSOS.curvilinearError().matrix()(1, 1)) * abs(sin(seedTSOS.globalMomentum().theta())));
0257   }
0258 }