File indexing completed on 2024-04-06 12:33:10
0001
0002
0003
0004
0005
0006
0007 #include "Validation/RecoMuon/src/MuonTrackAnalyzer.h"
0008
0009
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017
0018 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0019 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0020 #include "DataFormats/Math/interface/deltaR.h"
0021
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0023
0024 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0025
0026 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028 #include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
0029
0030 #include "Validation/RecoMuon/src/Histograms.h"
0031 #include "Validation/RecoMuon/src/HTrack.h"
0032
0033 #include "TFile.h"
0034 #include "TH1F.h"
0035 #include "TH2F.h"
0036
0037 using namespace std;
0038 using namespace edm;
0039
0040
0041 MuonTrackAnalyzer::MuonTrackAnalyzer(const ParameterSet &ps) {
0042
0043 pset = ps;
0044 ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
0045
0046 theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0047
0048 theSimTracksLabel = edm::InputTag("g4SimHits");
0049 theSimTracksToken = consumes<edm::SimTrackContainer>(theSimTracksLabel);
0050
0051 theTracksLabel = pset.getParameter<InputTag>("Tracks");
0052 theTracksToken = consumes<reco::TrackCollection>(theTracksLabel);
0053 doTracksAnalysis = pset.getUntrackedParameter<bool>("DoTracksAnalysis", true);
0054
0055 doSeedsAnalysis = pset.getUntrackedParameter<bool>("DoSeedsAnalysis", false);
0056 if (doSeedsAnalysis) {
0057 theSeedsLabel = pset.getParameter<InputTag>("MuonSeed");
0058 theSeedsToken = consumes<TrajectorySeedCollection>(theSeedsLabel);
0059 ParameterSet updatorPar = pset.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
0060 theSeedPropagatorName = updatorPar.getParameter<string>("Propagator");
0061
0062 theUpdator = new MuonUpdatorAtVertex(updatorPar, theService);
0063 }
0064
0065 theCSCSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
0066 theDTSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
0067 theRPCSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
0068 theCSCSimHitToken = consumes<std::vector<PSimHit> >(theCSCSimHitLabel);
0069 theDTSimHitToken = consumes<std::vector<PSimHit> >(theDTSimHitLabel);
0070 theRPCSimHitToken = consumes<std::vector<PSimHit> >(theRPCSimHitLabel);
0071
0072 theEtaRange = (EtaRange)pset.getParameter<int>("EtaRange");
0073
0074
0075 numberOfSimTracks = 0;
0076
0077 numberOfRecTracks = 0;
0078
0079 dbe_ = edm::Service<DQMStore>().operator->();
0080 out = pset.getUntrackedParameter<string>("rootFileName");
0081 dirName_ = pset.getUntrackedParameter<std::string>("dirName");
0082 subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0083 }
0084
0085
0086 MuonTrackAnalyzer::~MuonTrackAnalyzer() {
0087 if (theService)
0088 delete theService;
0089 }
0090
0091 void MuonTrackAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0092 edm::Run const &iRun,
0093 edm::EventSetup const & ) {
0094 ibooker.cd();
0095
0096 InputTag algo = theTracksLabel;
0097 string dirName = dirName_;
0098 if (!algo.process().empty())
0099 dirName += algo.process() + "_";
0100 if (!algo.label().empty())
0101 dirName += algo.label() + "_";
0102 if (!algo.instance().empty())
0103 dirName += algo.instance() + "";
0104 if (dirName.find("Tracks") < dirName.length()) {
0105 dirName.replace(dirName.find("Tracks"), 6, "");
0106 }
0107 std::replace(dirName.begin(), dirName.end(), ':', '_');
0108 ibooker.setCurrentFolder(dirName);
0109
0110
0111 std::string simName = dirName;
0112 simName += "/SimTracks";
0113 hSimTracks = new HTrackVariables(ibooker, simName, "SimTracks");
0114
0115 ibooker.cd();
0116 ibooker.setCurrentFolder(dirName);
0117
0118
0119
0120
0121 if (doSeedsAnalysis) {
0122 ibooker.cd();
0123 ibooker.setCurrentFolder(dirName);
0124 hRecoSeedInner = new HTrack(ibooker, dirName, "RecoSeed", "Inner");
0125 hRecoSeedPCA = new HTrack(ibooker, dirName, "RecoSeed", "PCA");
0126 }
0127
0128 if (doTracksAnalysis) {
0129 ibooker.cd();
0130 ibooker.setCurrentFolder(dirName);
0131 hRecoTracksPCA = new HTrack(ibooker, dirName, "RecoTracks", "PCA");
0132 hRecoTracksInner = new HTrack(ibooker, dirName, "RecoTracks", "Inner");
0133 hRecoTracksOuter = new HTrack(ibooker, dirName, "RecoTracks", "Outer");
0134
0135 ibooker.cd();
0136 ibooker.setCurrentFolder(dirName);
0137
0138
0139
0140 hChi2 = ibooker.book1D("chi2", "#chi^2", 200, 0, 200);
0141 hChi2VsEta = ibooker.book2D("chi2VsEta", "#chi^2 VS #eta", 120, -3., 3., 200, 0, 200);
0142
0143 hChi2Norm = ibooker.book1D("chi2Norm", "Normalized #chi^2", 400, 0, 100);
0144 hChi2NormVsEta = ibooker.book2D("chi2NormVsEta", "Normalized #chi^2 VS #eta", 120, -3., 3., 400, 0, 100);
0145
0146 hHitsPerTrack = ibooker.book1D("HitsPerTrack", "Number of hits per track", 55, 0, 55);
0147 hHitsPerTrackVsEta =
0148 ibooker.book2D("HitsPerTrackVsEta", "Number of hits per track VS #eta", 120, -3., 3., 55, 0, 55);
0149
0150 hDof = ibooker.book1D("dof", "Number of Degree of Freedom", 55, 0, 55);
0151 hDofVsEta = ibooker.book2D("dofVsEta", "Number of Degree of Freedom VS #eta", 120, -3., 3., 55, 0, 55);
0152
0153 hChi2Prob = ibooker.book1D("chi2Prob", "#chi^2 probability", 200, 0, 1);
0154 hChi2ProbVsEta = ibooker.book2D("chi2ProbVsEta", "#chi^2 probability VS #eta", 120, -3., 3., 200, 0, 1);
0155
0156 hNumberOfTracks = ibooker.book1D("NumberOfTracks", "Number of reconstructed tracks per event", 200, 0, 200);
0157 hNumberOfTracksVsEta = ibooker.book2D(
0158 "NumberOfTracksVsEta", "Number of reconstructed tracks per event VS #eta", 120, -3., 3., 10, 0, 10);
0159
0160 hChargeVsEta = ibooker.book2D("ChargeVsEta", "Charge vs #eta gen", 120, -3., 3., 4, -2., 2.);
0161 hChargeVsPt = ibooker.book2D("ChargeVsPt", "Charge vs P_{T} gen", 250, 0, 200, 4, -2., 2.);
0162 hPtRecVsPtGen = ibooker.book2D("PtRecVsPtGen", "P_{T} rec vs P_{T} gen", 250, 0, 200, 250, 0, 200);
0163
0164 hDeltaPtVsEta = ibooker.book2D("DeltaPtVsEta", "#Delta P_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
0165 hDeltaPt_In_Out_VsEta =
0166 ibooker.book2D("DeltaPt_In_Out_VsEta_", "P^{in}_{t} - P^{out}_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
0167 }
0168 }
0169
0170 void MuonTrackAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0171 LogDebug("MuonTrackAnalyzer") << "Run: " << event.id().run() << " Event: " << event.id().event();
0172
0173
0174 theService->update(eventSetup);
0175
0176 Handle<SimTrackContainer> simTracks;
0177 event.getByToken(theSimTracksToken, simTracks);
0178 fillPlots(event, simTracks);
0179
0180 if (doTracksAnalysis)
0181 tracksAnalysis(event, eventSetup, simTracks);
0182
0183 if (doSeedsAnalysis)
0184 seedsAnalysis(event, eventSetup, simTracks);
0185 }
0186
0187 void MuonTrackAnalyzer::seedsAnalysis(const Event &event,
0188 const EventSetup &eventSetup,
0189 Handle<SimTrackContainer> simTracks) {
0190 MuonPatternRecoDumper debug;
0191
0192
0193 Handle<TrajectorySeedCollection> seeds;
0194 event.getByToken(theSeedsToken, seeds);
0195
0196 LogTrace("MuonTrackAnalyzer") << "Number of reconstructed seeds: " << seeds->size() << endl;
0197
0198 for (TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed) {
0199 TrajectoryStateOnSurface seedTSOS = getSeedTSOS(*seed);
0200 pair<SimTrack, double> sim = getSimTrack(seedTSOS, simTracks);
0201 fillPlots(seedTSOS, sim.first, hRecoSeedInner, debug);
0202
0203 std::pair<bool, FreeTrajectoryState> propSeed = theUpdator->propagateToNominalLine(seedTSOS);
0204 if (propSeed.first)
0205 fillPlots(propSeed.second, sim.first, hRecoSeedPCA, debug);
0206 else
0207 LogTrace("MuonTrackAnalyzer") << "Error in seed propagation" << endl;
0208 }
0209 }
0210
0211 void MuonTrackAnalyzer::tracksAnalysis(const Event &event,
0212 const EventSetup &eventSetup,
0213 Handle<SimTrackContainer> simTracks) {
0214 MuonPatternRecoDumper debug;
0215
0216
0217 Handle<reco::TrackCollection> tracks;
0218 event.getByToken(theTracksToken, tracks);
0219
0220 LogTrace("MuonTrackAnalyzer") << "Reconstructed tracks: " << tracks->size() << endl;
0221 hNumberOfTracks->Fill(tracks->size());
0222
0223 if (!tracks->empty())
0224 numberOfRecTracks++;
0225
0226
0227 for (reco::TrackCollection::const_iterator t = tracks->begin(); t != tracks->end(); ++t) {
0228 reco::TransientTrack track(*t, &*theService->magneticField(), theService->trackingGeometry());
0229
0230 TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
0231 TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0232 TrajectoryStateOnSurface pcaTSOS = track.impactPointState();
0233
0234 pair<SimTrack, double> sim = getSimTrack(pcaTSOS, simTracks);
0235 SimTrack simTrack = sim.first;
0236 hNumberOfTracksVsEta->Fill(simTrack.momentum().eta(), tracks->size());
0237 fillPlots(track, simTrack);
0238
0239 LogTrace("MuonTrackAnalyzer") << "State at the outer surface: " << endl;
0240 fillPlots(outerTSOS, simTrack, hRecoTracksOuter, debug);
0241
0242 LogTrace("MuonTrackAnalyzer") << "State at the inner surface: " << endl;
0243 fillPlots(innerTSOS, simTrack, hRecoTracksInner, debug);
0244
0245 LogTrace("MuonTrackAnalyzer") << "State at PCA: " << endl;
0246 fillPlots(pcaTSOS, simTrack, hRecoTracksPCA, debug);
0247
0248 double deltaPt_in_out = innerTSOS.globalMomentum().perp() - outerTSOS.globalMomentum().perp();
0249 hDeltaPt_In_Out_VsEta->Fill(simTrack.momentum().eta(), deltaPt_in_out);
0250
0251 double deltaPt_pca_sim = pcaTSOS.globalMomentum().perp() - sqrt(simTrack.momentum().Perp2());
0252 hDeltaPtVsEta->Fill(simTrack.momentum().eta(), deltaPt_pca_sim);
0253
0254 hChargeVsEta->Fill(simTrack.momentum().eta(), pcaTSOS.charge());
0255
0256 hChargeVsPt->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.charge());
0257
0258 hPtRecVsPtGen->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.globalMomentum().perp());
0259 }
0260 LogTrace("MuonTrackAnalyzer") << "--------------------------------------------" << endl;
0261 }
0262
0263 void MuonTrackAnalyzer::fillPlots(const Event &event, edm::Handle<edm::SimTrackContainer> &simTracks) {
0264 if (!checkMuonSimHitPresence(event, simTracks))
0265 return;
0266
0267
0268 SimTrackContainer::const_iterator simTrack;
0269 LogTrace("MuonTrackAnalyzer") << "Simulated tracks: " << simTracks->size() << endl;
0270
0271 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
0272 if (abs((*simTrack).type()) == 13) {
0273 if (!isInTheAcceptance((*simTrack).momentum().eta()))
0274 continue;
0275
0276 numberOfSimTracks++;
0277
0278 LogTrace("MuonTrackAnalyzer") << "Simualted muon:" << endl;
0279 LogTrace("MuonTrackAnalyzer") << "Sim pT: " << sqrt((*simTrack).momentum().perp2()) << endl;
0280 LogTrace("MuonTrackAnalyzer") << "Sim Eta: " << (*simTrack).momentum().eta() << endl;
0281
0282 hSimTracks->Fill((*simTrack).momentum().mag(),
0283 sqrt((*simTrack).momentum().perp2()),
0284 (*simTrack).momentum().eta(),
0285 (*simTrack).momentum().phi(),
0286 -(*simTrack).type() / abs((*simTrack).type()));
0287 LogTrace("MuonTrackAnalyzer") << "hSimTracks filled" << endl;
0288 }
0289
0290 LogTrace("MuonTrackAnalyzer") << endl;
0291 }
0292
0293 void MuonTrackAnalyzer::fillPlots(reco::TransientTrack &track, SimTrack &simTrack) {
0294 LogTrace("MuonTrackAnalyzer") << "Analizer: New track, chi2: " << track.chi2() << " dof: " << track.ndof() << endl;
0295 hChi2->Fill(track.chi2());
0296 hDof->Fill(track.ndof());
0297 hChi2Norm->Fill(track.normalizedChi2());
0298 hHitsPerTrack->Fill(track.recHitsSize());
0299
0300 hChi2Prob->Fill(ChiSquaredProbability(track.chi2(), track.ndof()));
0301
0302 hChi2VsEta->Fill(simTrack.momentum().eta(), track.chi2());
0303 hChi2NormVsEta->Fill(simTrack.momentum().eta(), track.normalizedChi2());
0304 hChi2ProbVsEta->Fill(simTrack.momentum().eta(), ChiSquaredProbability(track.chi2(), track.ndof()));
0305 hHitsPerTrackVsEta->Fill(simTrack.momentum().eta(), track.recHitsSize());
0306 hDofVsEta->Fill(simTrack.momentum().eta(), track.ndof());
0307 }
0308
0309 void MuonTrackAnalyzer::fillPlots(TrajectoryStateOnSurface &recoTSOS,
0310 SimTrack &simTrack,
0311 HTrack *histo,
0312 MuonPatternRecoDumper &debug) {
0313 LogTrace("MuonTrackAnalyzer") << debug.dumpTSOS(recoTSOS) << endl;
0314 histo->Fill(recoTSOS);
0315
0316 GlobalVector tsosVect = recoTSOS.globalMomentum();
0317 math::XYZVectorD reco(tsosVect.x(), tsosVect.y(), tsosVect.z());
0318 double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
0319 histo->FillDeltaR(deltaRVal);
0320
0321 histo->computeResolutionAndPull(recoTSOS, simTrack);
0322 }
0323
0324 void MuonTrackAnalyzer::fillPlots(FreeTrajectoryState &recoFTS,
0325 SimTrack &simTrack,
0326 HTrack *histo,
0327 MuonPatternRecoDumper &debug) {
0328 LogTrace("MuonTrackAnalyzer") << debug.dumpFTS(recoFTS) << endl;
0329 histo->Fill(recoFTS);
0330
0331 GlobalVector ftsVect = recoFTS.momentum();
0332 math::XYZVectorD reco(ftsVect.x(), ftsVect.y(), ftsVect.z());
0333 double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
0334 histo->FillDeltaR(deltaRVal);
0335
0336 histo->computeResolutionAndPull(recoFTS, simTrack);
0337 }
0338
0339 pair<SimTrack, double> MuonTrackAnalyzer::getSimTrack(TrajectoryStateOnSurface &tsos,
0340 Handle<SimTrackContainer> simTracks) {
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356 SimTrackContainer::const_iterator simTrack;
0357
0358 SimTrack result;
0359
0360 double bestDeltaR = 10e5;
0361 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0362 if (abs((*simTrack).type()) != 13)
0363 continue;
0364
0365
0366 GlobalVector tsosVect = tsos.globalMomentum();
0367 math::XYZVectorD vect(tsosVect.x(), tsosVect.y(), tsosVect.z());
0368 double newDeltaR = deltaR<double>(vect.eta(), vect.phi(), simTrack->momentum().eta(), simTrack->momentum().phi());
0369
0370 if (newDeltaR < bestDeltaR) {
0371 LogTrace("MuonTrackAnalyzer") << "Matching Track with DeltaR = " << newDeltaR << endl;
0372 bestDeltaR = newDeltaR;
0373 result = *simTrack;
0374 }
0375 }
0376 return pair<SimTrack, double>(result, bestDeltaR);
0377 }
0378
0379 bool MuonTrackAnalyzer::isInTheAcceptance(double eta) {
0380 switch (theEtaRange) {
0381 case all:
0382 return (abs(eta) <= 2.4) ? true : false;
0383 case barrel:
0384 return (abs(eta) < 1.1) ? true : false;
0385 case endcap:
0386 return (abs(eta) >= 1.1 && abs(eta) <= 2.4) ? true : false;
0387 default: {
0388 LogTrace("MuonTrackAnalyzer") << "No correct Eta range selected!! " << endl;
0389 return false;
0390 }
0391 }
0392 }
0393
0394 bool MuonTrackAnalyzer::checkMuonSimHitPresence(const Event &event, edm::Handle<edm::SimTrackContainer> simTracks) {
0395
0396 Handle<PSimHitContainer> dtSimHits;
0397 event.getByToken(theDTSimHitToken, dtSimHits);
0398
0399 Handle<PSimHitContainer> cscSimHits;
0400 event.getByToken(theCSCSimHitToken, cscSimHits);
0401
0402 Handle<PSimHitContainer> rpcSimHits;
0403 event.getByToken(theRPCSimHitToken, rpcSimHits);
0404
0405 map<unsigned int, vector<const PSimHit *> > mapOfMuonSimHits;
0406
0407 for (PSimHitContainer::const_iterator simhit = dtSimHits->begin(); simhit != dtSimHits->end(); ++simhit) {
0408 if (abs(simhit->particleType()) != 13)
0409 continue;
0410 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0411 }
0412
0413 for (PSimHitContainer::const_iterator simhit = cscSimHits->begin(); simhit != cscSimHits->end(); ++simhit) {
0414 if (abs(simhit->particleType()) != 13)
0415 continue;
0416 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0417 }
0418
0419 for (PSimHitContainer::const_iterator simhit = rpcSimHits->begin(); simhit != rpcSimHits->end(); ++simhit) {
0420 if (abs(simhit->particleType()) != 13)
0421 continue;
0422 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0423 }
0424
0425 bool presence = false;
0426
0427 for (SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0428 if (abs(simTrack->type()) != 13)
0429 continue;
0430
0431 map<unsigned int, vector<const PSimHit *> >::const_iterator mapIterator =
0432 mapOfMuonSimHits.find(simTrack->trackId());
0433
0434 if (mapIterator != mapOfMuonSimHits.end())
0435 presence = true;
0436 }
0437
0438 return presence;
0439 }
0440
0441 TrajectoryStateOnSurface MuonTrackAnalyzer::getSeedTSOS(const TrajectorySeed &seed) {
0442
0443 PTrajectoryStateOnDet pTSOD = seed.startingState();
0444
0445
0446
0447 DetId seedDetId(pTSOD.detId());
0448
0449 const GeomDet *gdet = theService->trackingGeometry()->idToDet(seedDetId);
0450
0451 TrajectoryStateOnSurface initialState =
0452 trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*theService->magneticField());
0453
0454
0455 const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer(seedDetId);
0456
0457 PropagationDirection detLayerOrder = oppositeToMomentum;
0458
0459
0460 vector<const DetLayer *> detLayers;
0461 detLayers =
0462 theService->muonNavigationSchool()->compatibleLayers(*initialLayer, *initialState.freeState(), detLayerOrder);
0463
0464 TrajectoryStateOnSurface result = initialState;
0465 if (!detLayers.empty()) {
0466 const DetLayer *finalLayer = detLayers.back();
0467 const TrajectoryStateOnSurface propagatedState =
0468 theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->surface());
0469 if (propagatedState.isValid())
0470 result = propagatedState;
0471 }
0472
0473 return result;
0474 }