Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:10 

0001 #include "Validation/RecoMuon/src/MuonTrackResidualAnalyzer.h"
0002 
0003 // Collaborating Class Header
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011 
0012 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0013 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0014 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0015 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
0016 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
0017 
0018 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0019 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0020 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0021 
0022 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0023 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0024 #include "Validation/RecoMuon/src/Histograms.h"
0025 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0026 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0027 
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031 #include "TDirectory.h"
0032 
0033 using namespace std;
0034 using namespace edm;
0035 
0036 /// Constructor
0037 
0038 MuonTrackResidualAnalyzer::MuonTrackResidualAnalyzer(const edm::ParameterSet &ps) {
0039   pset = ps;
0040   // service parameters
0041   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
0042   // the services
0043   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0044 
0045   theMuonTrackLabel = pset.getParameter<InputTag>("MuonTrack");
0046   theMuonTrackToken = consumes<reco::TrackCollection>(theMuonTrackLabel);
0047 
0048   cscSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
0049   dtSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
0050   rpcSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
0051   theCSCSimHitToken = consumes<std::vector<PSimHit> >(cscSimHitLabel);
0052   theDTSimHitToken = consumes<std::vector<PSimHit> >(dtSimHitLabel);
0053   theRPCSimHitToken = consumes<std::vector<PSimHit> >(rpcSimHitLabel);
0054 
0055   dbe_ = edm::Service<DQMStore>().operator->();
0056   out = pset.getUntrackedParameter<string>("rootFileName");
0057   dirName_ = pset.getUntrackedParameter<std::string>("dirName");
0058   subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0059 
0060   // Sim or Real
0061   theDataType = pset.getParameter<InputTag>("DataType");
0062   if (theDataType.label() != "RealData" && theDataType.label() != "SimData")
0063     LogDebug("MuonTrackResidualAnalyzer") << "Error in Data Type!!";
0064   theDataTypeToken = consumes<edm::SimTrackContainer>(theDataType);
0065 
0066   theEtaRange = (EtaRange)pset.getParameter<int>("EtaRange");
0067 
0068   theUpdator = new KFUpdator();
0069   theEstimator = new Chi2MeasurementEstimator(100000.);
0070 
0071   theMuonSimHitNumberPerEvent = 0;
0072 }
0073 
0074 /// Destructor
0075 MuonTrackResidualAnalyzer::~MuonTrackResidualAnalyzer() {
0076   delete theUpdator;
0077   delete theEstimator;
0078   delete theService;
0079 }
0080 
0081 void MuonTrackResidualAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0082                                                edm::Run const &iRun,
0083                                                edm::EventSetup const & /* iSetup */) {
0084   LogDebug("MuonTrackResidualAnalyzer") << "Begin Run";
0085 
0086   ibooker.cd();
0087   InputTag algo = theMuonTrackLabel;
0088   string dirName = dirName_;
0089   if (!algo.process().empty())
0090     dirName += algo.process() + "_";
0091   if (!algo.label().empty())
0092     dirName += algo.label() + "_";
0093   if (!algo.instance().empty())
0094     dirName += algo.instance() + "";
0095   if (dirName.find("Tracks") < dirName.length()) {
0096     dirName.replace(dirName.find("Tracks"), 6, "");
0097   }
0098   std::replace(dirName.begin(), dirName.end(), ':', '_');
0099   ibooker.setCurrentFolder(dirName);
0100 
0101   hDPtRef = ibooker.book1D("DeltaPtRef", "P^{in}_{t}-P^{in ref}", 10000, -20, 20);
0102 
0103   // Resolution wrt the 1D Rec Hits
0104   //  h1DRecHitRes = new HResolution1DRecHit("TotalRec");
0105 
0106   // Resolution wrt the 1d Sim Hits
0107   //  h1DSimHitRes = new HResolution1DRecHit("TotalSim");
0108 
0109   hSimHitsPerTrack = ibooker.book1D("SimHitsPerTrack", "Number of sim hits per track", 55, 0, 55);
0110   hSimHitsPerTrackVsEta =
0111       ibooker.book2D("SimHitsPerTrackVsEta", "Number of sim hits per track VS #eta", 120, -3., 3., 55, 0, 55);
0112   hDeltaPtVsEtaSim =
0113       ibooker.book2D("DeltaPtVsEtaSim", "#Delta P_{t} vs #eta gen, sim quantity", 120, -3., 3., 500, -250., 250.);
0114   hDeltaPtVsEtaSim2 =
0115       ibooker.book2D("DeltaPtVsEtaSim2", "#Delta P_{t} vs #eta gen, sim quantity", 120, -3., 3., 500, -250., 250.);
0116 }
0117 
0118 void MuonTrackResidualAnalyzer::dqmEndRun(edm::Run const &, edm::EventSetup const &) {
0119   if (!out.empty() && dbe_)
0120     dbe_->save(out);
0121 }
0122 void MuonTrackResidualAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0123   LogDebug("MuonTrackResidualAnalyzer") << "Analyze";
0124 
0125   // Update the services
0126   theService->update(eventSetup);
0127   MuonPatternRecoDumper debug;
0128   theMuonSimHitNumberPerEvent = 0;
0129 
0130   // Get the SimHit collection from the event
0131   Handle<PSimHitContainer> dtSimHits;
0132   event.getByToken(theDTSimHitToken, dtSimHits);
0133 
0134   Handle<PSimHitContainer> cscSimHits;
0135   event.getByToken(theCSCSimHitToken, cscSimHits);
0136 
0137   Handle<PSimHitContainer> rpcSimHits;
0138   event.getByToken(theRPCSimHitToken, rpcSimHits);
0139 
0140   Handle<SimTrackContainer> simTracks;
0141 
0142   // FIXME Add the tracker one??
0143 
0144   // Map simhits per DetId
0145   map<DetId, const PSimHit *> muonSimHitsPerId = mapMuSimHitsPerId(dtSimHits, cscSimHits, rpcSimHits);
0146 
0147   hSimHitsPerTrack->Fill(theMuonSimHitNumberPerEvent);
0148 
0149   double etaSim = 0;
0150 
0151   if (theDataType.label() == "SimData") {
0152     // Get the SimTrack collection from the event
0153     event.getByToken(theDataTypeToken, simTracks);
0154 
0155     // Loop over the Sim tracks
0156     SimTrackContainer::const_iterator simTrack;
0157 
0158     for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
0159       if (abs((*simTrack).type()) == 13) {
0160         hSimHitsPerTrackVsEta->Fill((*simTrack).momentum().eta(), theMuonSimHitNumberPerEvent);
0161         etaSim = (*simTrack).momentum().eta();
0162         theSimTkId = (*simTrack).trackId();
0163       }
0164   }
0165 
0166   // Get the RecTrack collection from the event
0167   Handle<reco::TrackCollection> muonTracks;
0168   event.getByToken(theMuonTrackToken, muonTracks);
0169 
0170   reco::TrackCollection::const_iterator muonTrack;
0171 
0172   // Loop over the Rec tracks
0173   for (muonTrack = muonTracks->begin(); muonTrack != muonTracks->end(); ++muonTrack) {
0174     reco::TransientTrack track(*muonTrack, &*theService->magneticField(), theService->trackingGeometry());
0175 
0176     TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
0177     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0178 
0179     TransientTrackingRecHit::ConstRecHitContainer result;
0180 
0181     // SimHit Energy loss analysis
0182     double momAtEntry = -150., momAtExit = -150.;
0183 
0184     if (theSimHitContainer.size() > 1) {
0185       const GeomDet *geomDetA = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.front()->detUnitId()));
0186       double distA = geomDetA->toGlobal(theSimHitContainer.front()->localPosition()).mag();
0187 
0188       const GeomDet *geomDetB = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.back()->detUnitId()));
0189       double distB = geomDetB->toGlobal(theSimHitContainer.back()->localPosition()).mag();
0190 
0191       LogTrace("MuonTrackResidualAnalyzer")
0192           << "Inner SimHit: " << theSimHitContainer.front()->particleType()
0193           << " Pt: " << theSimHitContainer.front()->momentumAtEntry().perp()
0194           << " E: " << theSimHitContainer.front()->momentumAtEntry().perp() << " R: " << distA << endl;
0195       LogTrace("MuonTrackResidualAnalyzer")
0196           << "Outer SimHit: " << theSimHitContainer.back()->particleType()
0197           << " Pt: " << theSimHitContainer.back()->momentumAtEntry().perp()
0198           << " E: " << theSimHitContainer.front()->momentumAtEntry().perp() << " R: " << distB << endl;
0199 
0200       momAtEntry = theSimHitContainer.front()->momentumAtEntry().perp();
0201       momAtExit = theSimHitContainer.back()->momentumAtEntry().perp();
0202     }
0203 
0204     trackingRecHit_iterator rhFirst = track.recHitsBegin();
0205     trackingRecHit_iterator rhLast = track.recHitsEnd() - 1;
0206     map<DetId, const PSimHit *>::const_iterator itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
0207     map<DetId, const PSimHit *>::const_iterator itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
0208 
0209     double momAtEntry2 = -150, momAtExit2 = -150.;
0210     if (itFirst != muonSimHitsPerId.end())
0211       momAtEntry2 = itFirst->second->momentumAtEntry().perp();
0212     else {
0213       LogDebug("MuonTrackResidualAnalyzer") << "No first sim hit found";
0214       ++rhFirst;
0215       itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
0216       if (itFirst != muonSimHitsPerId.end())
0217         momAtEntry2 = itFirst->second->momentumAtEntry().perp();
0218       else {
0219         LogDebug("MuonTrackResidualAnalyzer") << "No second sim hit found";
0220         // continue;
0221       }
0222     }
0223 
0224     if (itLast != muonSimHitsPerId.end())
0225       momAtExit2 = itLast->second->momentumAtEntry().perp();
0226     else {
0227       LogDebug("MuonTrackResidualAnalyzer") << "No last sim hit found";
0228       --rhLast;
0229       itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
0230       if (itLast != muonSimHitsPerId.end())
0231         momAtExit2 = itLast->second->momentumAtEntry().perp();
0232       else {
0233         LogDebug("MuonTrackResidualAnalyzer") << "No last but one sim hit found";
0234         // continue;
0235       }
0236     }
0237 
0238     if (etaSim) {
0239       if (momAtEntry >= 0 && momAtExit >= 0)
0240         hDeltaPtVsEtaSim->Fill(etaSim, momAtEntry - momAtExit);
0241       if (momAtEntry2 >= 0 && momAtExit2 >= 0)
0242         hDeltaPtVsEtaSim2->Fill(etaSim, momAtEntry2 - momAtExit2);
0243     } else
0244       LogDebug("MuonTrackResidualAnalyzer") << "NO SimTrack'eta";
0245     //
0246 
0247     // computeResolution(trajectoryBW,muonSimHitsPerId,h1DSimHitRes);
0248     // computeResolution(smoothed,muonSimHitsPerId,h1DSimHitRes);
0249   }
0250 }
0251 
0252 bool MuonTrackResidualAnalyzer::isInTheAcceptance(double eta) {
0253   switch (theEtaRange) {
0254     case all:
0255       return (abs(eta) <= 2.4) ? true : false;
0256     case barrel:
0257       return (abs(eta) < 1.1) ? true : false;
0258     case endcap:
0259       return (abs(eta) >= 1.1 && abs(eta) <= 2.4) ? true : false;
0260     default: {
0261       LogDebug("MuonTrackResidualAnalyzer") << "No correct Eta range selected!! ";
0262       return false;
0263     }
0264   }
0265 }
0266 
0267 // map the muon simhits by id
0268 map<DetId, const PSimHit *> MuonTrackResidualAnalyzer::mapMuSimHitsPerId(Handle<PSimHitContainer> dtSimhits,
0269                                                                          Handle<PSimHitContainer> cscSimhits,
0270                                                                          Handle<PSimHitContainer> rpcSimhits) {
0271   MuonPatternRecoDumper debug;
0272 
0273   map<DetId, const PSimHit *> hitIdMap;
0274   theSimHitContainer.clear();
0275 
0276   mapMuSimHitsPerId(dtSimhits, hitIdMap);
0277   mapMuSimHitsPerId(cscSimhits, hitIdMap);
0278   mapMuSimHitsPerId(rpcSimhits, hitIdMap);
0279 
0280   if (theSimHitContainer.size() > 1)
0281     stable_sort(
0282         theSimHitContainer.begin(), theSimHitContainer.end(), RadiusComparatorInOut(theService->trackingGeometry()));
0283 
0284   LogDebug("MuonTrackResidualAnalyzer") << "Sim Hit list";
0285   int count = 1;
0286   for (vector<const PSimHit *>::const_iterator it = theSimHitContainer.begin(); it != theSimHitContainer.end(); ++it) {
0287     LogTrace("MuonTrackResidualAnalyzer")
0288         << count << " "
0289         << " Process Type: " << (*it)->processType() << " " << debug.dumpMuonId(DetId((*it)->detUnitId())) << endl;
0290   }
0291 
0292   return hitIdMap;
0293 }
0294 
0295 void MuonTrackResidualAnalyzer::mapMuSimHitsPerId(Handle<PSimHitContainer> simhits,
0296                                                   map<DetId, const PSimHit *> &hitIdMap) {
0297   for (PSimHitContainer::const_iterator simhit = simhits->begin(); simhit != simhits->end(); ++simhit) {
0298     if (abs(simhit->particleType()) != 13 && theSimTkId != simhit->trackId())
0299       continue;
0300 
0301     theSimHitContainer.push_back(&*simhit);
0302     DetId id = DetId(simhit->detUnitId());
0303 
0304     if (id.subdetId() == MuonSubdetId::DT) {
0305       DTLayerId lId(id.rawId());
0306       id = DetId(lId.rawId());
0307     }
0308 
0309     map<DetId, const PSimHit *>::const_iterator it = hitIdMap.find(id);
0310 
0311     if (it == hitIdMap.end())
0312       hitIdMap[id] = &*simhit;
0313     else
0314       LogDebug("MuonTrackResidualAnalyzer") << "TWO muons in the same sensible volume!!";
0315 
0316     ++theMuonSimHitNumberPerEvent;
0317   }
0318 }
0319 
0320 void MuonTrackResidualAnalyzer::computeResolution(Trajectory &trajectory,
0321                                                   map<DetId, const PSimHit *> &hitIdMap,
0322                                                   HResolution1DRecHit *histos) {
0323   Trajectory::DataContainer data = trajectory.measurements();
0324 
0325   for (Trajectory::DataContainer::const_iterator datum = data.begin(); datum != data.end(); ++datum) {
0326     GlobalPoint fitPoint = datum->updatedState().globalPosition();
0327 
0328     // FIXME!
0329     //     double errX = datum->updatedState().cartesianError().matrix()[0][0];
0330     //     double errY = datum->updatedState().cartesianError().matrix()[1][1];
0331     //     double errZ = datum->updatedState().cartesianError().matrix()[2][2];
0332     //
0333     double errX = datum->updatedState().localError().matrix()(3, 3);
0334     double errY = datum->updatedState().localError().matrix()(4, 4);
0335     double errZ = 1.;
0336 
0337     map<DetId, const PSimHit *>::const_iterator it = hitIdMap.find(datum->recHit()->geographicalId());
0338 
0339     if (it == hitIdMap.end())
0340       continue;  // FIXME! Put a counter
0341 
0342     const PSimHit *simhit = it->second;
0343 
0344     LocalPoint simHitPoint = simhit->localPosition();
0345 
0346     const GeomDet *geomDet = theService->trackingGeometry()->idToDet(DetId(simhit->detUnitId()));
0347 
0348     LocalPoint fitLocalPoint = geomDet->toLocal(fitPoint);
0349 
0350     LocalVector diff = fitLocalPoint - simHitPoint;
0351 
0352     cout << "SimHit position " << simHitPoint << endl;
0353     cout << "Fit position " << fitLocalPoint << endl;
0354     cout << "Fit position2 " << datum->updatedState().localPosition() << endl;
0355     cout << "Errors on the fit position: (" << errX << "," << errY << "," << errZ << ")" << endl;
0356     cout << "Resolution on x: " << diff.x() / abs(simHitPoint.x()) << endl;
0357     cout << "Resolution on y: " << diff.y() / abs(simHitPoint.y()) << endl;
0358     cout << "Resolution on z: " << diff.z() / abs(simHitPoint.z()) << endl;
0359 
0360     cout << "Eta direction: " << simhit->momentumAtEntry().eta() << " eta position: " << simHitPoint.eta() << endl;
0361     cout << "Phi direction: " << simhit->momentumAtEntry().phi() << " phi position: " << simHitPoint.phi() << endl;
0362 
0363     histos->Fill(simHitPoint.x(),
0364                  simHitPoint.y(),
0365                  simHitPoint.z(),
0366                  diff.x(),
0367                  diff.y(),
0368                  diff.z(),
0369                  errX,
0370                  errY,
0371                  errZ,
0372                  simhit->momentumAtEntry().eta(),
0373                  simhit->momentumAtEntry().phi());
0374     // simHitPoint.eta(), simHitPoint.phi() ); // FIXME!
0375   }
0376 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!