Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:08 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Jeremy Andrea
0005  */
0006 
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0009 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0010 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0011 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013 
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "DQMServices/Core/interface/DQMStore.h"
0017 #include "DQM/TrackingMonitor/interface/TrackEfficiencyMonitor.h"
0018 #include <string>
0019 
0020 // needed to compute the efficiency
0021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0022 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0023 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0024 #include "DataFormats/GeometrySurface/interface/Plane.h"
0025 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0026 #include "DataFormats/Math/interface/Vector3D.h"
0027 #include "DataFormats/MuonReco/interface/Muon.h"
0028 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0029 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0030 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0031 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0032 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0033 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
0034 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0035 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0036 
0037 //-----------------------------------------------------------------------------------
0038 TrackEfficiencyMonitor::TrackEfficiencyMonitor(const edm::ParameterSet& iConfig)
0039     : ttbToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0040       propToken_(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0041       navToken_(esConsumes(edm::ESInputTag("", "CosmicNavigationSchool"))),
0042       gstToken_(esConsumes()),
0043       mfToken_(esConsumes()),
0044       mtToken_(esConsumes())
0045 //-----------------------------------------------------------------------------------
0046 {
0047   dqmStore_ = edm::Service<DQMStore>().operator->();
0048 
0049   theRadius_ = iConfig.getParameter<double>("theRadius");
0050   theMaxZ_ = iConfig.getParameter<double>("theMaxZ");
0051   isBFieldOff_ = iConfig.getParameter<bool>("isBFieldOff");
0052   trackEfficiency_ = iConfig.getParameter<bool>("trackEfficiency");
0053   theTKTracksLabel_ = iConfig.getParameter<edm::InputTag>("TKTrackCollection");
0054   theSTATracksLabel_ = iConfig.getParameter<edm::InputTag>("STATrackCollection");
0055   muonToken_ = consumes<edm::View<reco::Muon> >(iConfig.getParameter<edm::InputTag>("muoncoll"));
0056 
0057   theTKTracksToken_ = consumes<reco::TrackCollection>(theTKTracksLabel_);
0058   theSTATracksToken_ = consumes<reco::TrackCollection>(theSTATracksLabel_);
0059 
0060   conf_ = iConfig;
0061 }
0062 
0063 //-----------------------------------------------------------------------------------
0064 TrackEfficiencyMonitor::~TrackEfficiencyMonitor()
0065 //-----------------------------------------------------------------------------------
0066 {}
0067 
0068 //-----------------------------------------------------------------------------------
0069 void TrackEfficiencyMonitor::bookHistograms(DQMStore::IBooker& ibooker,
0070                                             edm::Run const& /* iRun */,
0071                                             edm::EventSetup const& /* iSetup */)
0072 //-----------------------------------------------------------------------------------
0073 {
0074   std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0075   std::string AlgoName = conf_.getParameter<std::string>("AlgoName");
0076 
0077   ibooker.setCurrentFolder(MEFolderName);
0078 
0079   //
0080   int muonXBin = conf_.getParameter<int>("muonXBin");
0081   double muonXMin = conf_.getParameter<double>("muonXMin");
0082   double muonXMax = conf_.getParameter<double>("muonXMax");
0083 
0084   histname = "muonX_";
0085   muonX = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonXBin, muonXMin, muonXMax);
0086   muonX->setAxisTitle("");
0087 
0088   //
0089   int muonYBin = conf_.getParameter<int>("muonYBin");
0090   double muonYMin = conf_.getParameter<double>("muonYMin");
0091   double muonYMax = conf_.getParameter<double>("muonYMax");
0092 
0093   histname = "muonY_";
0094   muonY = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonYBin, muonYMin, muonYMax);
0095   muonY->setAxisTitle("");
0096 
0097   //
0098   int muonZBin = conf_.getParameter<int>("muonZBin");
0099   double muonZMin = conf_.getParameter<double>("muonZMin");
0100   double muonZMax = conf_.getParameter<double>("muonZMax");
0101 
0102   histname = "muonZ_";
0103   muonZ = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonZBin, muonZMin, muonZMax);
0104   muonZ->setAxisTitle("");
0105 
0106   //
0107   int muonEtaBin = conf_.getParameter<int>("muonEtaBin");
0108   double muonEtaMin = conf_.getParameter<double>("muonEtaMin");
0109   double muonEtaMax = conf_.getParameter<double>("muonEtaMax");
0110 
0111   histname = "muonEta_";
0112   muonEta = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonEtaBin, muonEtaMin, muonEtaMax);
0113   muonEta->setAxisTitle("");
0114 
0115   //
0116   int muonPhiBin = conf_.getParameter<int>("muonPhiBin");
0117   double muonPhiMin = conf_.getParameter<double>("muonPhiMin");
0118   double muonPhiMax = conf_.getParameter<double>("muonPhiMax");
0119 
0120   histname = "muonPhi_";
0121   muonPhi = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonPhiBin, muonPhiMin, muonPhiMax);
0122   muonPhi->setAxisTitle("");
0123 
0124   //
0125   int muonD0Bin = conf_.getParameter<int>("muonD0Bin");
0126   double muonD0Min = conf_.getParameter<double>("muonD0Min");
0127   double muonD0Max = conf_.getParameter<double>("muonD0Max");
0128 
0129   histname = "muonD0_";
0130   muonD0 = ibooker.book1D(histname + AlgoName, histname + AlgoName, muonD0Bin, muonD0Min, muonD0Max);
0131   muonD0->setAxisTitle("");
0132 
0133   //
0134   int muonCompatibleLayersBin = conf_.getParameter<int>("muonCompatibleLayersBin");
0135   double muonCompatibleLayersMin = conf_.getParameter<double>("muonCompatibleLayersMin");
0136   double muonCompatibleLayersMax = conf_.getParameter<double>("muonCompatibleLayersMax");
0137 
0138   histname = "muonCompatibleLayers_";
0139   muonCompatibleLayers = ibooker.book1D(histname + AlgoName,
0140                                         histname + AlgoName,
0141                                         muonCompatibleLayersBin,
0142                                         muonCompatibleLayersMin,
0143                                         muonCompatibleLayersMax);
0144   muonCompatibleLayers->setAxisTitle("");
0145 
0146   //------------------------------------------------------------------------------------
0147 
0148   //
0149   int trackXBin = conf_.getParameter<int>("trackXBin");
0150   double trackXMin = conf_.getParameter<double>("trackXMin");
0151   double trackXMax = conf_.getParameter<double>("trackXMax");
0152 
0153   histname = "trackX_";
0154   trackX = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackXBin, trackXMin, trackXMax);
0155   trackX->setAxisTitle("");
0156 
0157   //
0158   int trackYBin = conf_.getParameter<int>("trackYBin");
0159   double trackYMin = conf_.getParameter<double>("trackYMin");
0160   double trackYMax = conf_.getParameter<double>("trackYMax");
0161 
0162   histname = "trackY_";
0163   trackY = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackYBin, trackYMin, trackYMax);
0164   trackY->setAxisTitle("");
0165 
0166   //
0167   int trackZBin = conf_.getParameter<int>("trackZBin");
0168   double trackZMin = conf_.getParameter<double>("trackZMin");
0169   double trackZMax = conf_.getParameter<double>("trackZMax");
0170 
0171   histname = "trackZ_";
0172   trackZ = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackZBin, trackZMin, trackZMax);
0173   trackZ->setAxisTitle("");
0174 
0175   //
0176   int trackEtaBin = conf_.getParameter<int>("trackEtaBin");
0177   double trackEtaMin = conf_.getParameter<double>("trackEtaMin");
0178   double trackEtaMax = conf_.getParameter<double>("trackEtaMax");
0179 
0180   histname = "trackEta_";
0181   trackEta = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackEtaBin, trackEtaMin, trackEtaMax);
0182   trackEta->setAxisTitle("");
0183 
0184   //
0185   int trackPhiBin = conf_.getParameter<int>("trackPhiBin");
0186   double trackPhiMin = conf_.getParameter<double>("trackPhiMin");
0187   double trackPhiMax = conf_.getParameter<double>("trackPhiMax");
0188 
0189   histname = "trackPhi_";
0190   trackPhi = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackPhiBin, trackPhiMin, trackPhiMax);
0191   trackPhi->setAxisTitle("");
0192 
0193   //
0194   int trackD0Bin = conf_.getParameter<int>("trackD0Bin");
0195   double trackD0Min = conf_.getParameter<double>("trackD0Min");
0196   double trackD0Max = conf_.getParameter<double>("trackD0Max");
0197 
0198   histname = "trackD0_";
0199   trackD0 = ibooker.book1D(histname + AlgoName, histname + AlgoName, trackD0Bin, trackD0Min, trackD0Max);
0200   trackD0->setAxisTitle("");
0201 
0202   //
0203   int trackCompatibleLayersBin = conf_.getParameter<int>("trackCompatibleLayersBin");
0204   double trackCompatibleLayersMin = conf_.getParameter<double>("trackCompatibleLayersMin");
0205   double trackCompatibleLayersMax = conf_.getParameter<double>("trackCompatibleLayersMax");
0206 
0207   histname = "trackCompatibleLayers_";
0208   trackCompatibleLayers = ibooker.book1D(histname + AlgoName,
0209                                          histname + AlgoName,
0210                                          trackCompatibleLayersBin,
0211                                          trackCompatibleLayersMin,
0212                                          trackCompatibleLayersMax);
0213   trackCompatibleLayers->setAxisTitle("");
0214 
0215   //------------------------------------------------------------------------------------
0216 
0217   //
0218   int deltaXBin = conf_.getParameter<int>("deltaXBin");
0219   double deltaXMin = conf_.getParameter<double>("deltaXMin");
0220   double deltaXMax = conf_.getParameter<double>("deltaXMax");
0221 
0222   histname = "deltaX_";
0223   deltaX = ibooker.book1D(histname + AlgoName, histname + AlgoName, deltaXBin, deltaXMin, deltaXMax);
0224   deltaX->setAxisTitle("");
0225 
0226   //
0227   int deltaYBin = conf_.getParameter<int>("deltaYBin");
0228   double deltaYMin = conf_.getParameter<double>("deltaYMin");
0229   double deltaYMax = conf_.getParameter<double>("deltaYMax");
0230 
0231   histname = "deltaY_";
0232   deltaY = ibooker.book1D(histname + AlgoName, histname + AlgoName, deltaYBin, deltaYMin, deltaYMax);
0233   deltaY->setAxisTitle("");
0234 
0235   //
0236   int signDeltaXBin = conf_.getParameter<int>("signDeltaXBin");
0237   double signDeltaXMin = conf_.getParameter<double>("signDeltaXMin");
0238   double signDeltaXMax = conf_.getParameter<double>("signDeltaXMax");
0239 
0240   histname = "signDeltaX_";
0241   signDeltaX = ibooker.book1D(histname + AlgoName, histname + AlgoName, signDeltaXBin, signDeltaXMin, signDeltaXMax);
0242   signDeltaX->setAxisTitle("");
0243 
0244   //
0245   int signDeltaYBin = conf_.getParameter<int>("signDeltaYBin");
0246   double signDeltaYMin = conf_.getParameter<double>("signDeltaYMin");
0247   double signDeltaYMax = conf_.getParameter<double>("signDeltaYMax");
0248 
0249   histname = "signDeltaY_";
0250   signDeltaY = ibooker.book1D(histname + AlgoName, histname + AlgoName, signDeltaYBin, signDeltaYMin, signDeltaYMax);
0251   signDeltaY->setAxisTitle("");
0252 
0253   histname = "GlobalMuonPtEtaPhi_LowPt_";
0254   GlobalMuonPtEtaPhiLowPt = ibooker.book2D(histname + AlgoName, histname + AlgoName, 20, -2.4, 2.4, 20, -3.25, 3.25);
0255   GlobalMuonPtEtaPhiLowPt->setAxisTitle("");
0256 
0257   histname = "StandaloneMuonPtEtaPhi_LowPt_";
0258   StandaloneMuonPtEtaPhiLowPt =
0259       ibooker.book2D(histname + AlgoName, histname + AlgoName, 20, -2.4, 2.4, 20, -3.25, 3.25);
0260   StandaloneMuonPtEtaPhiLowPt->setAxisTitle("");
0261 
0262   histname = "GlobalMuonPtEtaPhi_HighPt_";
0263   GlobalMuonPtEtaPhiHighPt = ibooker.book2D(histname + AlgoName, histname + AlgoName, 20, -2.4, 2.4, 20, -3.25, 3.25);
0264   GlobalMuonPtEtaPhiHighPt->setAxisTitle("");
0265 
0266   histname = "StandaloneMuonPtEtaPhi_HighPt_";
0267   StandaloneMuonPtEtaPhiHighPt =
0268       ibooker.book2D(histname + AlgoName, histname + AlgoName, 20, -2.4, 2.4, 20, -3.25, 3.25);
0269   StandaloneMuonPtEtaPhiHighPt->setAxisTitle("");
0270 }
0271 
0272 //-----------------------------------------------------------------------------------
0273 void TrackEfficiencyMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
0274 //-----------------------------------------------------------------------------------
0275 {
0276   edm::Handle<reco::TrackCollection> tkTracks = iEvent.getHandle(theTKTracksToken_);
0277   edm::Handle<reco::TrackCollection> staTracks = iEvent.getHandle(theSTATracksToken_);
0278 
0279   //initialize values
0280   failedToPropagate = 0;
0281   nCompatibleLayers = 0;
0282   findDetLayer = false;
0283 
0284   const NavigationSchool& nav = iSetup.getData(navToken_);
0285   measurementTracker = &iSetup.getData(mtToken_);
0286   theTTrackBuilder = &iSetup.getData(ttbToken_);
0287   thePropagator = &iSetup.getData(propToken_);
0288   bField = &iSetup.getData(mfToken_);
0289   theTracker = &iSetup.getData(gstToken_);
0290   theNavigation = new DirectTrackerNavigation(theTracker);
0291 
0292   edm::Handle<edm::View<reco::Muon> > muons = iEvent.getHandle(muonToken_);
0293   if (!muons.isValid())
0294     return;
0295   for (edm::View<reco::Muon>::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0296     if ((*muon).pt() < 5)
0297       continue;
0298     if (fabs((*muon).eta()) > 2.4)
0299       continue;
0300     if ((*muon).vertexNormalizedChi2() > 10)
0301       continue;
0302     if ((*muon).isStandAloneMuon() and (*muon).isGlobalMuon()) {
0303       if ((*muon).pt() < 20)
0304         GlobalMuonPtEtaPhiLowPt->Fill((*muon).eta(), (*muon).phi());
0305       else
0306         GlobalMuonPtEtaPhiHighPt->Fill((*muon).eta(), (*muon).phi());
0307     }
0308     if ((*muon).isStandAloneMuon()) {
0309       if ((*muon).pt() < 20)
0310         StandaloneMuonPtEtaPhiLowPt->Fill((*muon).eta(), (*muon).phi());
0311       else
0312         StandaloneMuonPtEtaPhiHighPt->Fill((*muon).eta(), (*muon).phi());
0313     }
0314   }
0315   if (trackEfficiency_) {
0316     //---------------------------------------------------
0317     // Select muons with good quality
0318     // If B field is on, no up-down matching between the muons
0319     //---------------------------------------------------
0320     bool isGoodMuon = false;
0321     double mudd0 = 0., mudphi = 0., muddsz = 0., mudeta = 0.;
0322     if (isBFieldOff_) {
0323       if (staTracks->size() == 2) {
0324         for (unsigned int bindex = 0; bindex < staTracks->size(); ++bindex) {
0325           if (0 == bindex) {
0326             mudd0 += (*staTracks)[bindex].d0();
0327             mudphi += (*staTracks)[bindex].phi();
0328             muddsz += (*staTracks)[bindex].dsz();
0329             mudeta += (*staTracks)[bindex].eta();
0330           }
0331           if (1 == bindex) {
0332             mudd0 -= (*staTracks)[bindex].d0();
0333             mudphi -= (*staTracks)[bindex].phi();
0334             muddsz -= (*staTracks)[bindex].dsz();
0335             mudeta -= (*staTracks)[bindex].eta();
0336           }
0337         }
0338         if ((fabs(mudd0) < 15.0) && (fabs(mudphi) < 0.045) && (fabs(muddsz) < 20.0) && (fabs(mudeta) < 0.060))
0339           isGoodMuon = true;
0340       }
0341 
0342       if (isGoodMuon)
0343         testTrackerTracks(tkTracks, staTracks, nav);
0344 
0345     } else if (staTracks->size() == 1 || staTracks->size() == 2)
0346       testTrackerTracks(tkTracks, staTracks, nav);
0347   }
0348 
0349   if (!trackEfficiency_ && tkTracks->size() == 1) {
0350     if ((tkTracks->front()).normalizedChi2() < 5 && (tkTracks->front()).hitPattern().numberOfValidHits() > 8)
0351       testSTATracks(tkTracks, staTracks);
0352   }
0353 
0354   delete theNavigation;
0355 }
0356 
0357 //-----------------------------------------------------------------------------------
0358 void TrackEfficiencyMonitor::testTrackerTracks(edm::Handle<reco::TrackCollection> tkTracks,
0359                                                edm::Handle<reco::TrackCollection> staTracks,
0360                                                const NavigationSchool& navigationSchool)
0361 //-----------------------------------------------------------------------------------
0362 {
0363   const std::string metname = "testStandAloneMuonTracks";
0364 
0365   //---------------------------------------------------
0366   // get the index of the "up" muon
0367   // histograms will only be computed for the "up" muon
0368   //---------------------------------------------------
0369 
0370   int nUpMuon = 0;
0371   int idxUpMuon = -1;
0372   for (unsigned int i = 0; i < staTracks->size(); i++) {
0373     if (checkSemiCylinder((*staTracks)[i]) == TrackEfficiencyMonitor::Up) {
0374       nUpMuon++;
0375       idxUpMuon = i;
0376     }
0377   }
0378 
0379   if (nUpMuon == 1) {
0380     //---------------------------------------------------
0381     //get the muon informations
0382     //---------------------------------------------------
0383 
0384     reco::TransientTrack staTT = theTTrackBuilder->build((*staTracks)[idxUpMuon]);
0385     double ipX = staTT.impactPointState().globalPosition().x();
0386     double ipY = staTT.impactPointState().globalPosition().y();
0387     double ipZ = staTT.impactPointState().globalPosition().z();
0388     double eta = staTT.impactPointState().globalDirection().eta();
0389     double phi = staTT.impactPointState().globalDirection().phi();
0390     double d0 = (*staTracks)[idxUpMuon].d0();
0391 
0392     TrajectoryStateOnSurface theTSOS = staTT.outermostMeasurementState();
0393     TrajectoryStateOnSurface theTSOSCompLayers = staTT.outermostMeasurementState();
0394 
0395     //---------------------------------------------------
0396     //check if the muon is in the tracker acceptance
0397     //---------------------------------------------------
0398     bool isInTrackerAcceptance = false;
0399     isInTrackerAcceptance = trackerAcceptance(theTSOS, theRadius_, theMaxZ_);
0400 
0401     //---------------------------------------------------st
0402     //count the number of compatible layers
0403     //---------------------------------------------------
0404     nCompatibleLayers = compatibleLayers(navigationSchool, theTSOSCompLayers);
0405 
0406     if (isInTrackerAcceptance && (*staTracks)[idxUpMuon].hitPattern().numberOfValidHits() > 28) {
0407       //---------------------------------------------------
0408       //count the number of good muon candidates
0409       //---------------------------------------------------
0410 
0411       TrajectoryStateOnSurface staState;
0412       LocalVector diffLocal;
0413 
0414       bool isTrack = false;
0415       if (!tkTracks->empty()) {
0416         //---------------------------------------------------
0417         //look for associated tracks
0418         //---------------------------------------------------
0419         float DR2min = 1000;
0420         reco::TrackCollection::const_iterator closestTrk = tkTracks->end();
0421 
0422         for (reco::TrackCollection::const_iterator tkTrack = tkTracks->begin(); tkTrack != tkTracks->end(); ++tkTrack) {
0423           reco::TransientTrack tkTT = theTTrackBuilder->build(*tkTrack);
0424           TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
0425           staState = thePropagator->propagate(staTT.outermostMeasurementState(), tkInner.surface());
0426           failedToPropagate = 1;
0427 
0428           if (staState.isValid()) {
0429             failedToPropagate = 0;
0430             diffLocal = tkInner.localPosition() - staState.localPosition();
0431             double DR2 = diffLocal.x() * diffLocal.x() + diffLocal.y() * diffLocal.y();
0432             if (DR2 < DR2min) {
0433               DR2min = DR2;
0434               closestTrk = tkTrack;
0435             }
0436             if (pow(DR2, 0.5) < 100.)
0437               isTrack = true;
0438           }
0439         }
0440 
0441         if (DR2min != 1000) {
0442           reco::TransientTrack tkTT = theTTrackBuilder->build(*closestTrk);
0443           TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
0444           staState = thePropagator->propagate(staTT.outermostMeasurementState(), tkInner.surface());
0445           deltaX->Fill(diffLocal.x());
0446           deltaY->Fill(diffLocal.y());
0447           signDeltaX->Fill(diffLocal.x() /
0448                            (tkInner.localError().positionError().xx() + staState.localError().positionError().xx()));
0449           signDeltaY->Fill(diffLocal.y() /
0450                            (tkInner.localError().positionError().yy() + staState.localError().positionError().yy()));
0451         }
0452       }
0453 
0454       if (failedToPropagate == 0) {
0455         muonX->Fill(ipX);
0456         muonY->Fill(ipY);
0457         muonZ->Fill(ipZ);
0458         muonEta->Fill(eta);
0459         muonPhi->Fill(phi);
0460         muonD0->Fill(d0);
0461         muonCompatibleLayers->Fill(nCompatibleLayers);
0462 
0463         if (isTrack) {
0464           trackX->Fill(ipX);
0465           trackY->Fill(ipY);
0466           trackZ->Fill(ipZ);
0467           trackEta->Fill(eta);
0468           trackPhi->Fill(phi);
0469           trackD0->Fill(d0);
0470           trackCompatibleLayers->Fill(nCompatibleLayers);
0471         }
0472       }
0473     }
0474   }
0475 }
0476 
0477 //-----------------------------------------------------------------------------------
0478 void TrackEfficiencyMonitor::testSTATracks(edm::Handle<reco::TrackCollection> tkTracks,
0479                                            edm::Handle<reco::TrackCollection> staTracks)
0480 //-----------------------------------------------------------------------------------
0481 {
0482   reco::TransientTrack tkTT = theTTrackBuilder->build(tkTracks->front());
0483   double ipX = tkTT.impactPointState().globalPosition().x();
0484   double ipY = tkTT.impactPointState().globalPosition().y();
0485   double ipZ = tkTT.impactPointState().globalPosition().z();
0486   double eta = tkTT.impactPointState().globalDirection().eta();
0487   double phi = tkTT.impactPointState().globalDirection().phi();
0488   double d0 = (*tkTracks)[0].d0();
0489 
0490   TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
0491   LocalVector diffLocal;
0492   TrajectoryStateOnSurface staState;
0493   bool isTrack = false;
0494 
0495   if (!staTracks->empty()) {
0496     //---------------------------------------------------
0497     //look for associated muons
0498     //---------------------------------------------------
0499 
0500     float DR2min = 1000;
0501     reco::TrackCollection::const_iterator closestTrk = staTracks->end();
0502     //----------------------loop on tracker tracks:
0503     for (reco::TrackCollection::const_iterator staTrack = staTracks->begin(); staTrack != staTracks->end();
0504          ++staTrack) {
0505       if (checkSemiCylinder(*staTrack) == TrackEfficiencyMonitor::Up) {
0506         reco::TransientTrack staTT = theTTrackBuilder->build(*staTrack);
0507         failedToPropagate = 1;
0508         staState = thePropagator->propagate(staTT.outermostMeasurementState(), tkInner.surface());
0509 
0510         if (staState.isValid()) {
0511           failedToPropagate = 0;
0512           diffLocal = tkInner.localPosition() - staState.localPosition();
0513 
0514           double DR2 = diffLocal.x() * diffLocal.x() + diffLocal.y() * diffLocal.y();
0515           if (DR2 < DR2min) {
0516             DR2min = DR2;
0517             closestTrk = staTrack;
0518           }
0519           if (pow(DR2, 0.5) < 100.)
0520             isTrack = true;
0521         }
0522       }
0523     }
0524   }
0525 
0526   if (failedToPropagate == 0) {
0527     trackX->Fill(ipX);
0528     trackY->Fill(ipY);
0529     trackZ->Fill(ipZ);
0530     trackEta->Fill(eta);
0531     trackPhi->Fill(phi);
0532     trackD0->Fill(d0);
0533 
0534     if (isTrack) {
0535       muonX->Fill(ipX);
0536       muonY->Fill(ipY);
0537       muonZ->Fill(ipZ);
0538       muonEta->Fill(eta);
0539       muonPhi->Fill(phi);
0540       muonD0->Fill(d0);
0541     }
0542   }
0543 }
0544 
0545 //-----------------------------------------------------------------------------------
0546 TrackEfficiencyMonitor::SemiCylinder TrackEfficiencyMonitor::checkSemiCylinder(const reco::Track& tk)
0547 //-----------------------------------------------------------------------------------
0548 {
0549   return tk.innerPosition().phi() > 0 ? TrackEfficiencyMonitor::Up : TrackEfficiencyMonitor::Down;
0550 }
0551 
0552 //-----------------------------------------------------------------------------------
0553 bool TrackEfficiencyMonitor::trackerAcceptance(TrajectoryStateOnSurface theTSOS, double theRadius, double theMaxZ)
0554 //-----------------------------------------------------------------------------------
0555 {
0556   //---------------------------------------------------
0557   // check if the muon is in the tracker acceptance
0558   // check the compatibility with a cylinder of radius "theRadius"
0559   //---------------------------------------------------
0560 
0561   //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
0562 
0563   //Propagator*  theTmpPropagator = &*thePropagator;
0564   Propagator* theTmpPropagator = thePropagator->clone();
0565 
0566   if (theTSOS.globalPosition().y() < 0)
0567     theTmpPropagator->setPropagationDirection(oppositeToMomentum);
0568   else
0569     theTmpPropagator->setPropagationDirection(alongMomentum);
0570 
0571   Cylinder::PositionType pos0;
0572   Cylinder::RotationType rot0;
0573   const Cylinder::CylinderPointer cyl = Cylinder::build(pos0, rot0, theRadius);
0574   TrajectoryStateOnSurface tsosAtCyl = theTmpPropagator->propagate(*theTSOS.freeState(), *cyl);
0575   double accept = false;
0576   if (tsosAtCyl.isValid()) {
0577     if (fabs(tsosAtCyl.globalPosition().z()) < theMaxZ) {
0578       Cylinder::PositionType pos02;
0579       Cylinder::RotationType rot02;
0580       const Cylinder::CylinderPointer cyl2 = Cylinder::build(pos02, rot02, theRadius - 10);
0581       TrajectoryStateOnSurface tsosAtCyl2 = theTmpPropagator->propagate(*tsosAtCyl.freeState(), *cyl2);
0582       if (tsosAtCyl2.isValid()) {
0583         Cylinder::PositionType pos03;
0584         Cylinder::RotationType rot03;
0585         const Cylinder::CylinderPointer cyl3 = Cylinder::build(pos03, rot03, theRadius);
0586         TrajectoryStateOnSurface tsosAtCyl3 = theTmpPropagator->propagate(*tsosAtCyl2.freeState(), *cyl3);
0587         if (tsosAtCyl3.isValid()) {
0588           accept = true;
0589         }
0590       }
0591     }
0592   }
0593   delete theTmpPropagator;
0594   //muon propagated to the barrel cylinder
0595   return accept;
0596 }
0597 
0598 //-----------------------------------------------------------------------------------
0599 int TrackEfficiencyMonitor::compatibleLayers(const NavigationSchool& navigationSchool, TrajectoryStateOnSurface theTSOS)
0600 //-----------------------------------------------------------------------------------
0601 {
0602   //---------------------------------------------------
0603   // check the number of compatible layers
0604   //---------------------------------------------------
0605 
0606   std::vector<const BarrelDetLayer*> barrelTOBLayers = measurementTracker->geometricSearchTracker()->tobLayers();
0607 
0608   unsigned int layers = 0;
0609   for (unsigned int k = 0; k < barrelTOBLayers.size(); k++) {
0610     const DetLayer* firstLay = barrelTOBLayers[barrelTOBLayers.size() - 1 - k];
0611 
0612     //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
0613 
0614     Propagator* theTmpPropagator = thePropagator->clone();
0615     theTmpPropagator->setPropagationDirection(alongMomentum);
0616 
0617     TrajectoryStateOnSurface startTSOS = theTmpPropagator->propagate(*theTSOS.freeState(), firstLay->surface());
0618 
0619     std::vector<const DetLayer*> trackCompatibleLayers;
0620 
0621     findDetLayer = true;
0622     bool isUpMuon = false;
0623     bool firstdtep = true;
0624 
0625     if (startTSOS.isValid()) {
0626       if (firstdtep)
0627         layers++;
0628 
0629       int nwhile = 0;
0630 
0631       //for other compatible layers
0632       while (startTSOS.isValid() && firstLay && findDetLayer) {
0633         if (firstdtep && startTSOS.globalPosition().y() > 0)
0634           isUpMuon = true;
0635 
0636         if (firstdtep) {
0637           std::vector<const DetLayer*> firstCompatibleLayers;
0638           firstCompatibleLayers.push_back(firstLay);
0639           std::pair<TrajectoryStateOnSurface, const DetLayer*> nextLayer =
0640               findNextLayer(theTSOS, firstCompatibleLayers, isUpMuon);
0641           firstdtep = false;
0642         } else {
0643           trackCompatibleLayers = navigationSchool.nextLayers(*firstLay, *(startTSOS.freeState()), alongMomentum);
0644           if (!trackCompatibleLayers.empty()) {
0645             std::pair<TrajectoryStateOnSurface, const DetLayer*> nextLayer =
0646                 findNextLayer(startTSOS, trackCompatibleLayers, isUpMuon);
0647             if (firstLay != nextLayer.second) {
0648               firstLay = nextLayer.second;
0649               startTSOS = nextLayer.first;
0650               layers++;
0651             } else
0652               firstLay = nullptr;
0653           }
0654         }
0655         nwhile++;
0656         if (nwhile > 100)
0657           break;
0658       }
0659       delete theTmpPropagator;
0660       break;
0661     }
0662     delete theTmpPropagator;
0663   }
0664   return layers;
0665 }
0666 
0667 //-----------------------------------------------------------------------------------
0668 std::pair<TrajectoryStateOnSurface, const DetLayer*> TrackEfficiencyMonitor::findNextLayer(
0669     TrajectoryStateOnSurface sTSOS, const std::vector<const DetLayer*>& trackCompatibleLayers, bool isUpMuon)
0670 //-----------------------------------------------------------------------------------
0671 {
0672   //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
0673 
0674   Propagator* theTmpPropagator = thePropagator->clone();
0675   theTmpPropagator->setPropagationDirection(alongMomentum);
0676 
0677   std::vector<const DetLayer*>::const_iterator itl;
0678   findDetLayer = false;
0679   for (itl = trackCompatibleLayers.begin(); itl != trackCompatibleLayers.end(); ++itl) {
0680     TrajectoryStateOnSurface tsos = theTmpPropagator->propagate(*(sTSOS.freeState()), (**itl).surface());
0681     if (tsos.isValid()) {
0682       sTSOS = tsos;
0683       findDetLayer = true;
0684 
0685       break;
0686     }
0687   }
0688   std::pair<TrajectoryStateOnSurface, const DetLayer*> blabla;
0689   blabla.first = sTSOS;
0690   blabla.second = &**itl;
0691   delete theTmpPropagator;
0692   return blabla;
0693 }
0694 
0695 DEFINE_FWK_MODULE(TrackEfficiencyMonitor);
[ 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!