Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-11 03:33:48 

0001 ////////////////////////////////////////////////////////////////////////////////
0002 // Package:          CalibTracker/SiStripHitEfficiency
0003 // Class:            HitEff
0004 // Original Author:  Keith Ulmer--University of Colorado
0005 //                   keith.ulmer@colorado.edu
0006 //
0007 ///////////////////////////////////////////////////////////////////////////////
0008 
0009 // system include files
0010 #include <memory>
0011 #include <string>
0012 #include <iostream>
0013 
0014 // user includes
0015 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0016 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0017 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
0018 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0019 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h"
0020 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
0021 #include "CalibTracker/SiStripHitEfficiency/plugins/HitEff.h"
0022 #include "DataFormats/Common/interface/DetSetVector.h"
0023 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0024 #include "DataFormats/Common/interface/Handle.h"
0025 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0027 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0028 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0029 #include "DataFormats/MuonReco/interface/Muon.h"
0030 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0031 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0033 #include "DataFormats/TrackReco/interface/DeDxData.h"
0034 #include "DataFormats/TrackReco/interface/Track.h"
0035 #include "DataFormats/TrackReco/interface/TrackBase.h"
0036 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0037 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0038 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0039 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0040 #include "FWCore/Framework/interface/ESHandle.h"
0041 #include "FWCore/Framework/interface/Event.h"
0042 #include "FWCore/Framework/interface/EventSetup.h"
0043 #include "FWCore/Framework/interface/Frameworkfwd.h"
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0046 #include "FWCore/Utilities/interface/Exception.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0049 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0051 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0052 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0053 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0054 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0055 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0056 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0057 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0058 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0059 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0060 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0061 
0062 // ROOT includes
0063 #include "TMath.h"
0064 #include "TH1F.h"
0065 
0066 // custom made printout
0067 #define LOGPRINT edm::LogPrint("SiStripHitEfficiency:HitEff")
0068 
0069 //
0070 // constructors and destructor
0071 //
0072 
0073 using namespace std;
0074 HitEff::HitEff(const edm::ParameterSet& conf)
0075     : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))),
0076       metaDataToken_(consumes<OnlineLuminosityRecord>(conf.getParameter<edm::InputTag>("metadata"))),
0077       commonModeToken_(mayConsume<edm::DetSetVector<SiStripRawDigi> >(conf.getParameter<edm::InputTag>("commonMode"))),
0078       siStripClusterInfo_(consumesCollector()),
0079       combinatorialTracks_token_(
0080           consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))),
0081       trajectories_token_(consumes<std::vector<Trajectory> >(conf.getParameter<edm::InputTag>("trajectories"))),
0082       trajTrackAsso_token_(consumes<TrajTrackAssociationCollection>(conf.getParameter<edm::InputTag>("trajectories"))),
0083       clusters_token_(
0084           consumes<edmNew::DetSetVector<SiStripCluster> >(conf.getParameter<edm::InputTag>("siStripClusters"))),
0085       digisCol_token_(consumes(conf.getParameter<edm::InputTag>("siStripDigis"))),
0086       digisVec_token_(consumes(conf.getParameter<edm::InputTag>("siStripDigis"))),
0087       trackerEvent_token_(consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("trackerEvent"))),
0088       topoToken_(esConsumes()),
0089       geomToken_(esConsumes()),
0090       cpeToken_(esConsumes(edm::ESInputTag("", "StripCPEfromTrackAngle"))),
0091       siStripQualityToken_(esConsumes()),
0092       magFieldToken_(esConsumes()),
0093       measurementTkToken_(esConsumes()),
0094       chi2MeasurementEstimatorToken_(esConsumes(edm::ESInputTag("", "Chi2"))),
0095       propagatorToken_(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))),
0096       conf_(conf) {
0097   usesResource(TFileService::kSharedResource);
0098   compSettings = conf_.getUntrackedParameter<int>("CompressionSettings", -1);
0099   layers = conf_.getParameter<int>("Layer");
0100   DEBUG = conf_.getParameter<bool>("Debug");
0101   addLumi_ = conf_.getUntrackedParameter<bool>("addLumi", false);
0102   addCommonMode_ = conf_.getUntrackedParameter<bool>("addCommonMode", false);
0103   cutOnTracks_ = conf_.getUntrackedParameter<bool>("cutOnTracks", false);
0104   trackMultiplicityCut_ = conf.getUntrackedParameter<unsigned int>("trackMultiplicity", 100);
0105   useFirstMeas_ = conf_.getUntrackedParameter<bool>("useFirstMeas", false);
0106   useLastMeas_ = conf_.getUntrackedParameter<bool>("useLastMeas", false);
0107   useAllHitsFromTracksWithMissingHits_ =
0108       conf_.getUntrackedParameter<bool>("useAllHitsFromTracksWithMissingHits", false);
0109   doMissingHitsRecovery_ = conf_.getUntrackedParameter<bool>("doMissingHitsRecovery", false);
0110 
0111   hitRecoveryCounters.resize(k_END_OF_LAYERS, 0);
0112   hitTotalCounters.resize(k_END_OF_LAYERS, 0);
0113 }
0114 
0115 void HitEff::beginJob() {
0116   edm::Service<TFileService> fs;
0117   if (compSettings > 0) {
0118     edm::LogInfo("SiStripHitEfficiency:HitEff") << "the compressions settings are:" << compSettings << std::endl;
0119     fs->file().SetCompressionSettings(compSettings);
0120   }
0121 
0122   traj = fs->make<TTree>("traj", "tree of trajectory positions");
0123 #ifdef ExtendedCALIBTree
0124   traj->Branch("timeDT", &timeDT, "timeDT/F");
0125   traj->Branch("timeDTErr", &timeDTErr, "timeDTErr/F");
0126   traj->Branch("timeDTDOF", &timeDTDOF, "timeDTDOF/I");
0127   traj->Branch("timeECAL", &timeECAL, "timeECAL/F");
0128   traj->Branch("dedx", &dedx, "dedx/F");
0129   traj->Branch("dedxNOM", &dedxNOM, "dedxNOM/I");
0130   traj->Branch("nLostHits", &nLostHits, "nLostHits/I");
0131   traj->Branch("chi2", &chi2, "chi2/F");
0132   traj->Branch("p", &p, "p/F");
0133 #endif
0134   traj->Branch("TrajGlbX", &TrajGlbX, "TrajGlbX/F");
0135   traj->Branch("TrajGlbY", &TrajGlbY, "TrajGlbY/F");
0136   traj->Branch("TrajGlbZ", &TrajGlbZ, "TrajGlbZ/F");
0137   traj->Branch("TrajLocX", &TrajLocX, "TrajLocX/F");
0138   traj->Branch("TrajLocY", &TrajLocY, "TrajLocY/F");
0139   traj->Branch("TrajLocAngleX", &TrajLocAngleX, "TrajLocAngleX/F");
0140   traj->Branch("TrajLocAngleY", &TrajLocAngleY, "TrajLocAngleY/F");
0141   traj->Branch("TrajLocErrX", &TrajLocErrX, "TrajLocErrX/F");
0142   traj->Branch("TrajLocErrY", &TrajLocErrY, "TrajLocErrY/F");
0143   traj->Branch("ClusterLocX", &ClusterLocX, "ClusterLocX/F");
0144   traj->Branch("ClusterLocY", &ClusterLocY, "ClusterLocY/F");
0145   traj->Branch("ClusterLocErrX", &ClusterLocErrX, "ClusterLocErrX/F");
0146   traj->Branch("ClusterLocErrY", &ClusterLocErrY, "ClusterLocErrY/F");
0147   traj->Branch("ClusterStoN", &ClusterStoN, "ClusterStoN/F");
0148   traj->Branch("ResX", &ResX, "ResX/F");
0149   traj->Branch("ResXSig", &ResXSig, "ResXSig/F");
0150   traj->Branch("ModIsBad", &ModIsBad, "ModIsBad/i");
0151   traj->Branch("SiStripQualBad", &SiStripQualBad, "SiStripQualBad/i");
0152   traj->Branch("withinAcceptance", &withinAcceptance, "withinAcceptance/O");
0153   traj->Branch("nHits", &nHits, "nHits/I");
0154   traj->Branch("pT", &pT, "pT/F");
0155   traj->Branch("highPurity", &highPurity, "highPurity/O");
0156   traj->Branch("trajHitValid", &trajHitValid, "trajHitValid/i");
0157   traj->Branch("Id", &Id, "Id/i");
0158   traj->Branch("run", &run, "run/i");
0159   traj->Branch("event", &event, "event/i");
0160   traj->Branch("layer", &whatlayer, "layer/i");
0161   traj->Branch("tquality", &tquality, "tquality/I");
0162   traj->Branch("bunchx", &bunchx, "bunchx/I");
0163   if (addLumi_) {
0164     traj->Branch("instLumi", &instLumi, "instLumi/F");
0165     traj->Branch("PU", &PU, "PU/F");
0166   }
0167   if (addCommonMode_)
0168     traj->Branch("commonMode", &commonMode, "commonMode/F");
0169 
0170   events = 0;
0171   EventTrackCKF = 0;
0172 
0173   totalNbHits = 0;
0174   missHitPerLayer.resize(k_END_OF_LAYERS, 0);
0175 }
0176 
0177 void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
0178   //Retrieve tracker topology from geometry
0179   const TrackerTopology* tTopo = &es.getData(topoToken_);
0180   siStripClusterInfo_.initEvent(es);
0181 
0182   //  bool DEBUG = false;
0183 
0184   LogDebug("SiStripHitEfficiency:HitEff") << "beginning analyze from HitEff" << endl;
0185 
0186   using namespace edm;
0187   using namespace reco;
0188   // Step A: Get Inputs
0189 
0190   int run_nr = e.id().run();
0191   int ev_nr = e.id().event();
0192   int bunch_nr = e.bunchCrossing();
0193 
0194   // Luminosity informations
0195   edm::Handle<LumiScalersCollection> lumiScalers = e.getHandle(scalerToken_);
0196   edm::Handle<OnlineLuminosityRecord> metaData = e.getHandle(metaDataToken_);
0197 
0198   instLumi = 0;
0199   PU = 0;
0200   if (addLumi_) {
0201     if (lumiScalers.isValid() && !lumiScalers->empty()) {
0202       if (lumiScalers->begin() != lumiScalers->end()) {
0203         instLumi = lumiScalers->begin()->instantLumi();
0204         PU = lumiScalers->begin()->pileup();
0205       }
0206     } else if (metaData.isValid()) {
0207       instLumi = metaData->instLumi();
0208       PU = metaData->avgPileUp();
0209     } else {
0210       edm::LogWarning("SiStripHitEfficiencyWorker") << "could not find a source for the Luminosity and PU";
0211     }
0212   }
0213 
0214   // CM
0215   edm::Handle<edm::DetSetVector<SiStripRawDigi> > commonModeDigis;
0216   if (addCommonMode_)
0217     e.getByToken(commonModeToken_, commonModeDigis);
0218 
0219   //CombinatoriaTrack
0220   edm::Handle<reco::TrackCollection> trackCollectionCKF;
0221   //edm::InputTag TkTagCKF = conf_.getParameter<edm::InputTag>("combinatorialTracks");
0222   e.getByToken(combinatorialTracks_token_, trackCollectionCKF);
0223 
0224   edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
0225   //edm::InputTag TkTrajCKF = conf_.getParameter<edm::InputTag>("trajectories");
0226   e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0227 
0228   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0229   e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0230 
0231   // Clusters
0232   // get the SiStripClusters from the event
0233   edm::Handle<edmNew::DetSetVector<SiStripCluster> > theClusters;
0234   //e.getByLabel("siStripClusters", theClusters);
0235   e.getByToken(clusters_token_, theClusters);
0236 
0237   //get tracker geometry
0238   edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_);
0239   const TrackerGeometry* tkgeom = &(*tracker);
0240 
0241   //get Cluster Parameter Estimator
0242   //std::string cpe = conf_.getParameter<std::string>("StripCPE");
0243   edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_);
0244   const StripClusterParameterEstimator& stripcpe(*parameterestimator);
0245 
0246   // get the SiStripQuality records
0247   edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_);
0248 
0249   const MagneticField* magField_ = &es.getData(magFieldToken_);
0250 
0251   // get the list of module IDs with FED-detected errors
0252   //  - In Aug-2023, the data format was changed from DetIdCollection to DetIdVector.
0253   //  - To provide some level of backward-compatibility,
0254   //    the plugin checks for both types giving preference to the new format.
0255   //  - If only the old format is available, the collection is
0256   //    converted to the new format, then used downstream.
0257   auto const& fedErrorIdsCol_h = e.getHandle(digisCol_token_);
0258   auto const& fedErrorIdsVec_h = e.getHandle(digisVec_token_);
0259   if (not fedErrorIdsCol_h.isValid() and not fedErrorIdsVec_h.isValid()) {
0260     throw cms::Exception("InvalidProductSiStripDetIdsWithFEDErrors")
0261         << "no valid product for SiStrip DetIds with FED errors (see parameter \"siStripDigis\"), "
0262            "neither for new format (DetIdVector) nor old format (DetIdCollection)";
0263   }
0264   auto const& fedErrorIds = fedErrorIdsVec_h.isValid() ? *fedErrorIdsVec_h : fedErrorIdsCol_h->as_vector();
0265 
0266   edm::ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTkToken_);
0267 
0268   edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
0269   //e.getByLabel("MeasurementTrackerEvent", measurementTrackerEvent);
0270   e.getByToken(trackerEvent_token_, measurementTrackerEvent);
0271 
0272   const MeasurementEstimator* estimator = &es.getData(chi2MeasurementEstimatorToken_);
0273   const Propagator* thePropagator = &es.getData(propagatorToken_);
0274 
0275   events++;
0276 
0277   // *************** SiStripCluster Collection
0278   const edmNew::DetSetVector<SiStripCluster>& input = *theClusters;
0279 
0280   //go through clusters to write out global position of good clusters for the layer understudy for comparison
0281   // Loop through clusters just to print out locations
0282   // Commented out to avoid discussion, should really be deleted.
0283   /*
0284   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end(); DSViter++) {
0285     // DSViter is a vector of SiStripClusters located on a single module
0286     unsigned int ClusterId = DSViter->id();
0287     DetId ClusterDetId(ClusterId);
0288     const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
0289     
0290     edmNew::DetSet<SiStripCluster>::const_iterator begin=DSViter->begin();
0291     edmNew::DetSet<SiStripCluster>::const_iterator end  =DSViter->end();
0292     for(edmNew::DetSet<SiStripCluster>::const_iterator iter=begin;iter!=end;++iter) {
0293       //iter is a single SiStripCluster
0294       StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
0295       
0296       const Surface* surface;
0297       surface = &(tracker->idToDet(ClusterDetId)->surface());
0298       LocalPoint lp = parameters.first;
0299       GlobalPoint gp = surface->toGlobal(lp);
0300       unsigned int layer = ::checkLayer(ClusterId, tTopo);
0301             if(DEBUG) LOGPRINT << "Found hit in cluster collection layer = " << layer << " with id = " << ClusterId << "   local X position = " << lp.x() << " +- " << sqrt(parameters.second.xx()) << "   matched/stereo/rphi = " << ((ClusterId & 0x3)==0) << "/" << ((ClusterId & 0x3)==1) << "/" << ((ClusterId & 0x3)==2) << endl;
0302     }
0303   }
0304   */
0305 
0306   // Tracking
0307   const reco::TrackCollection* tracksCKF = trackCollectionCKF.product();
0308   LogDebug("SiStripHitEfficiency:HitEff") << "number ckf tracks found = " << tracksCKF->size() << endl;
0309   //if (tracksCKF->size() == 1 ){
0310   if (!tracksCKF->empty()) {
0311     if (cutOnTracks_ && (tracksCKF->size() >= trackMultiplicityCut_))
0312       return;
0313     if (cutOnTracks_)
0314       LogDebug("SiStripHitEfficiency:HitEff")
0315           << "starting checking good event with < " << trackMultiplicityCut_ << " tracks" << endl;
0316 
0317     EventTrackCKF++;
0318 
0319 #ifdef ExtendedCALIBTree
0320     //get dEdx info if available
0321     edm::Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
0322     if (e.getByLabel("dedxMedianCTF", dEdxUncalibHandle)) {
0323       const ValueMap<DeDxData> dEdxTrackUncalib = *dEdxUncalibHandle.product();
0324 
0325       reco::TrackRef itTrack = reco::TrackRef(trackCollectionCKF, 0);
0326       dedx = dEdxTrackUncalib[itTrack].dEdx();
0327       dedxNOM = dEdxTrackUncalib[itTrack].numberOfMeasurements();
0328     } else {
0329       dedx = -999.0;
0330       dedxNOM = -999;
0331     }
0332 
0333     //get muon and ecal timing info if available
0334     edm::Handle<MuonCollection> muH;
0335     if (e.getByLabel("muonsWitht0Correction", muH)) {
0336       const MuonCollection& muonsT0 = *muH.product();
0337       if (!muonsT0.empty()) {
0338         MuonTime mt0 = muonsT0[0].time();
0339         timeDT = mt0.timeAtIpInOut;
0340         timeDTErr = mt0.timeAtIpInOutErr;
0341         timeDTDOF = mt0.nDof;
0342 
0343         bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
0344         if (hasCaloEnergyInfo)
0345           timeECAL = muonsT0[0].calEnergy().ecal_time;
0346       }
0347     } else {
0348       timeDT = -999.0;
0349       timeDTErr = -999.0;
0350       timeDTDOF = -999;
0351       timeECAL = -999.0;
0352     }
0353 
0354 #endif
0355     // actually should do a loop over all the tracks in the event here
0356 
0357     // Looping over traj-track associations to be able to get traj & track informations
0358     for (TrajTrackAssociationCollection::const_iterator it = trajTrackAssociationHandle->begin();
0359          it != trajTrackAssociationHandle->end();
0360          it++) {
0361       edm::Ref<std::vector<Trajectory> > itraj = it->key;
0362       reco::TrackRef itrack = it->val;
0363 
0364       // for each track, fill some variables such as number of hits and momentum
0365       nHits = itraj->foundHits();
0366 #ifdef ExtendedCALIBTree
0367       nLostHits = itraj->lostHits();
0368       chi2 = (itraj->chiSquared() / itraj->ndof());
0369       p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
0370 #endif
0371       pT = sqrt((itraj->lastMeasurement().updatedState().globalMomentum().x() *
0372                  itraj->lastMeasurement().updatedState().globalMomentum().x()) +
0373                 (itraj->lastMeasurement().updatedState().globalMomentum().y() *
0374                  itraj->lastMeasurement().updatedState().globalMomentum().y()));
0375 
0376       // track quality
0377       highPurity = itrack->quality(reco::TrackBase::TrackQuality::highPurity);
0378 
0379       std::vector<TrajectoryMeasurement> TMeas = itraj->measurements();
0380       totalNbHits += int(TMeas.size());
0381       vector<TrajectoryMeasurement>::iterator itm;
0382       double xloc = 0.;
0383       double yloc = 0.;
0384       double xErr = 0.;
0385       double yErr = 0.;
0386       double angleX = -999.;
0387       double angleY = -999.;
0388       double xglob, yglob, zglob;
0389 
0390       // Check whether the trajectory has some missing hits
0391       bool hasMissingHits = false;
0392       int previous_layer = 999;
0393       vector<unsigned int> missedLayers;
0394       for (const auto& itm : TMeas) {
0395         auto theHit = itm.recHit();
0396         unsigned int iidd = theHit->geographicalId().rawId();
0397         int layer = ::checkLayer(iidd, tTopo);
0398         int missedLayer = (layer + 1);
0399         int previousMissedLayer = (layer + 2);
0400         int diffPreviousLayer = (layer - previous_layer);
0401         if (doMissingHitsRecovery_) {
0402           //Layers from TIB + TOB
0403           if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
0404             missHitPerLayer[missedLayer] += 1;
0405             hasMissingHits = true;
0406           }
0407           //Layers from TID
0408           else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
0409             missHitPerLayer[missedLayer] += 1;
0410             hasMissingHits = true;
0411           }
0412           //Layers from TEC
0413           else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
0414             missHitPerLayer[missedLayer] += 1;
0415             hasMissingHits = true;
0416           }
0417 
0418           //##### TID Layer 11 in transition TID -> TIB : layer is in TIB, previous layer  = 12
0419           if ((layer > k_LayersStart && layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
0420             missHitPerLayer[11] += 1;
0421             hasMissingHits = true;
0422           }
0423 
0424           //##### TEC Layer 14 in transition TEC -> TOB : layer is in TOB, previous layer =  15
0425           if ((layer > k_LayersAtTIBEnd && layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
0426             missHitPerLayer[14] += 1;
0427             hasMissingHits = true;
0428           }
0429 
0430           //####### Consecutive missing hits case #######
0431 
0432           //##### Layers from TIB + TOB
0433           if (diffPreviousLayer == -3 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd &&
0434               previousMissedLayer > k_LayersStart && previousMissedLayer < k_LayersAtTOBEnd) {
0435             missHitPerLayer[missedLayer] += 1;
0436             missHitPerLayer[previousMissedLayer] += 1;
0437             hasMissingHits = true;
0438           }
0439 
0440           //##### Layers from TEC
0441           else if (diffPreviousLayer == -3 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd &&
0442                    previousMissedLayer > k_LayersAtTIDEnd && previousMissedLayer <= k_LayersAtTECEnd) {
0443             missHitPerLayer[missedLayer] += 1;
0444             missHitPerLayer[previousMissedLayer] += 1;
0445             hasMissingHits = true;
0446           }
0447         }
0448 
0449         if (theHit->getType() == TrackingRecHit::Type::missing)
0450           hasMissingHits = true;
0451 
0452         if (hasMissingHits)
0453           missedLayers.push_back(layer);
0454         previous_layer = layer;
0455       }
0456 
0457       // Loop on each measurement and take it into consideration
0458       //--------------------------------------------------------
0459       unsigned int prev_TKlayers = 0;
0460       for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
0461         auto theInHit = (*itm).recHit();
0462 
0463         LogDebug("SiStripHitEfficiency:HitEff") << "theInHit is valid = " << theInHit->isValid() << endl;
0464 
0465         unsigned int iidd = theInHit->geographicalId().rawId();
0466         bool foundConsMissingHits = false;
0467         unsigned int TKlayers = ::checkLayer(iidd, tTopo);
0468         LogDebug("SiStripHitEfficiency:HitEff") << "TKlayer from trajectory: " << TKlayers << "  from module = " << iidd
0469                                                 << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0470                                                 << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
0471 
0472         // Test first and last points of the trajectory
0473         // the list of measurements starts from outer layers  !!! This could change -> should add a check
0474         bool isFirstMeas = (itm == (TMeas.end() - 1));
0475         bool isLastMeas = (itm == (TMeas.begin()));
0476 
0477         if (!useFirstMeas_ && isFirstMeas)
0478           continue;
0479         if (!useLastMeas_ && isLastMeas)
0480           continue;
0481 
0482         // In case of missing hit in the track, check whether to use the other hits or not.
0483         if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
0484             !useAllHitsFromTracksWithMissingHits_)
0485           continue;
0486 
0487         // If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
0488         if (TKlayers == 10 || TKlayers == 22) {
0489           LogDebug("SiStripHitEfficiency:HitEff") << "skipping original TM for TOB 6 or TEC 9" << endl;
0490           continue;
0491         }
0492 
0493         // Make vector of TrajectoryAtInvalidHits to hold the trajectories
0494         std::vector<TrajectoryAtInvalidHit> TMs;
0495 
0496         // Make AnalyticalPropagator to use in TAVH constructor
0497         AnalyticalPropagator propagator(magField_, anyDirection);
0498 
0499         // for double sided layers check both sensors--if no hit was found on either sensor surface,
0500         // the trajectory measurements only have one invalid hit entry on the matched surface
0501         // so get the TrajectoryAtInvalidHit for both surfaces and include them in the study
0502         if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0503           // do hit eff check twice--once for each sensor
0504           //add a TM for each surface
0505           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0506           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0507         } else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
0508           // if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
0509           // the matched layer should be added to the study as well
0510           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0511           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0512           LogDebug("SiStripHitEfficiency:HitEff") << " found a hit with a missing partner";
0513         } else {
0514           //only add one TM for the single surface and the other will be added in the next iteration
0515           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator));
0516         }
0517         bool missingHitAdded = false;
0518 
0519         vector<TrajectoryMeasurement> tmpTmeas, prev_tmpTmeas;
0520         unsigned int misLayer = TKlayers + 1;
0521         unsigned int previousMisLayer = TKlayers + 2;
0522         //Use bool doMissingHitsRecovery to add possible missing hits based on actual/previous hit
0523         if (doMissingHitsRecovery_) {
0524           if (int(TKlayers) - int(prev_TKlayers) == -2) {
0525             const DetLayer* detlayer = itm->layer();
0526             const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0527             const TrajectoryStateOnSurface tsos = itm->updatedState();
0528             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0529 
0530             if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) {  //TEC
0531               std::vector<ForwardDetLayer const*> negTECLayers =
0532                   measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0533               std::vector<ForwardDetLayer const*> posTECLayers =
0534                   measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0535               const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0536               const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0537               if (tTopo->tecSide(iidd) == 1) {
0538                 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
0539               } else if (tTopo->tecSide(iidd) == 2) {
0540                 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
0541               }
0542             }
0543 
0544             else if (misLayer == (k_LayersAtTIDEnd - 1) ||
0545                      misLayer == k_LayersAtTIDEnd) {  // This is for TID layers 12 and 13
0546 
0547               std::vector<ForwardDetLayer const*> negTIDLayers =
0548                   measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
0549               std::vector<ForwardDetLayer const*> posTIDLayers =
0550                   measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
0551               const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0552               const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0553 
0554               if (tTopo->tidSide(iidd) == 1) {
0555                 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
0556               } else if (tTopo->tidSide(iidd) == 2) {
0557                 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
0558               }
0559             }
0560 
0561             if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) {  // Barrel
0562 
0563               std::vector<BarrelDetLayer const*> barrelTIBLayers =
0564                   measurementTrackerHandle->geometricSearchTracker()->tibLayers();
0565               std::vector<BarrelDetLayer const*> barrelTOBLayers =
0566                   measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0567 
0568               if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
0569                 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
0570                 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, *thePropagator, *estimator);
0571               } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
0572                 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
0573                 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, *thePropagator, *estimator);
0574               }
0575             }
0576           }
0577           if ((int(TKlayers) > k_LayersStart && int(TKlayers) <= k_LayersAtTIBEnd) && int(prev_TKlayers) == 12) {
0578             const DetLayer* detlayer = itm->layer();
0579             const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0580             const TrajectoryStateOnSurface tsos = itm->updatedState();
0581             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0582             std::vector<ForwardDetLayer const*> negTIDLayers =
0583                 measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
0584             std::vector<ForwardDetLayer const*> posTIDLayers =
0585                 measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
0586 
0587             const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
0588             const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
0589             if (tTopo->tidSide(iidd) == 1) {
0590               tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
0591             } else if (tTopo->tidSide(iidd) == 2) {
0592               tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
0593             }
0594           }
0595 
0596           if ((int(TKlayers) > k_LayersAtTIBEnd && int(TKlayers) <= k_LayersAtTOBEnd) && int(prev_TKlayers) == 15) {
0597             const DetLayer* detlayer = itm->layer();
0598             const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0599             const TrajectoryStateOnSurface tsos = itm->updatedState();
0600             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0601 
0602             std::vector<ForwardDetLayer const*> negTECLayers =
0603                 measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0604             std::vector<ForwardDetLayer const*> posTECLayers =
0605                 measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0606 
0607             const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
0608             const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
0609             if (tTopo->tecSide(iidd) == 1) {
0610               tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
0611             } else if (tTopo->tecSide(iidd) == 2) {
0612               tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
0613             }
0614           }
0615 
0616           //Test for two consecutive missing hits
0617           if (int(TKlayers) - int(prev_TKlayers) == -3) {
0618             foundConsMissingHits = true;
0619             const DetLayer* detlayer = itm->layer();
0620             const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0621             const TrajectoryStateOnSurface tsos = itm->updatedState();
0622             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0623 
0624             if (misLayer > k_LayersStart && misLayer <= k_LayersAtTOBEnd && previousMisLayer > k_LayersStart &&
0625                 previousMisLayer <= k_LayersAtTOBEnd) {  //Barrel case
0626               std::vector<BarrelDetLayer const*> barrelTIBLayers =
0627                   measurementTrackerHandle->geometricSearchTracker()->tibLayers();
0628               std::vector<BarrelDetLayer const*> barrelTOBLayers =
0629                   measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0630               if (misLayer > k_LayersStart && misLayer < k_LayersAtTIBEnd) {
0631                 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
0632                 const DetLayer* prevTibLayer = barrelTIBLayers[previousMisLayer - k_LayersStart - 1];
0633 
0634                 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, *thePropagator, *estimator);
0635                 prev_tmpTmeas = layerMeasurements.measurements(*prevTibLayer, tsos, *thePropagator, *estimator);
0636               } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
0637                 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
0638                 const DetLayer* prevTobLayer = barrelTOBLayers[previousMisLayer - k_LayersAtTIBEnd - 1];
0639                 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, *thePropagator, *estimator);
0640                 prev_tmpTmeas = layerMeasurements.measurements(*prevTobLayer, tsos, *thePropagator, *estimator);
0641               }
0642             } else if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd &&
0643                        previousMisLayer > k_LayersAtTIDEnd && previousMisLayer < k_LayersAtTECEnd) {  //TEC
0644               std::vector<ForwardDetLayer const*> negTECLayers =
0645                   measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0646               std::vector<ForwardDetLayer const*> posTECLayers =
0647                   measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0648 
0649               const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0650               const DetLayer* prevTecLayerneg = negTECLayers[previousMisLayer - k_LayersAtTIDEnd - 1];
0651 
0652               const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0653               const DetLayer* prevTecLayerpos = posTECLayers[previousMisLayer - k_LayersAtTIDEnd - 1];
0654 
0655               if (tTopo->tecSide(iidd) == 1) {
0656                 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
0657                 prev_tmpTmeas = layerMeasurements.measurements(*prevTecLayerneg, tsos, *thePropagator, *estimator);
0658               } else if (tTopo->tecSide(iidd) == 2) {
0659                 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
0660                 prev_tmpTmeas = layerMeasurements.measurements(*prevTecLayerpos, tsos, *thePropagator, *estimator);
0661               }
0662             }
0663           }
0664           if (!tmpTmeas.empty() && !foundConsMissingHits) {
0665             TrajectoryMeasurement TM_tmp(tmpTmeas.back());
0666             unsigned int iidd_tmp = TM_tmp.recHit()->geographicalId().rawId();
0667             if (iidd_tmp != 0) {
0668               LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0669               if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0670                 TMs.clear();
0671               if (::isDoubleSided(iidd_tmp, tTopo)) {
0672                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 1));
0673                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 2));
0674               } else
0675                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator));
0676               missingHitAdded = true;
0677               hitRecoveryCounters[misLayer] += 1;
0678             }
0679           }
0680 
0681           if (!tmpTmeas.empty() && !prev_tmpTmeas.empty() &&
0682               foundConsMissingHits) {  //Found two consecutive missing hits
0683             TrajectoryMeasurement TM_tmp1(tmpTmeas.back());
0684             TrajectoryMeasurement TM_tmp2(prev_tmpTmeas.back());
0685             //Inner and outer hits module IDs
0686             unsigned int modIdInner = TM_tmp1.recHit()->geographicalId().rawId();
0687             unsigned int modIdOuter = TM_tmp2.recHit()->geographicalId().rawId();
0688             bool innerModInactive = false, outerModInactive = false;
0689             for (const auto& tm : tmpTmeas) {  //Check if inner module is inactive
0690               unsigned int tmModId = tm.recHit()->geographicalId().rawId();
0691               if (tmModId == modIdInner && tm.recHit()->getType() == 2) {
0692                 innerModInactive = true;
0693                 break;
0694               }
0695             }
0696             for (const auto& tm : prev_tmpTmeas) {  //Check if outer module is inactive
0697               unsigned int tmModId = tm.recHit()->geographicalId().rawId();
0698               if (tmModId == modIdOuter && tm.recHit()->getType() == 2) {
0699                 outerModInactive = true;
0700                 break;  //Found the inactive module
0701               }
0702             }
0703 
0704             if (outerModInactive) {  //If outer missing hit is in inactive module, recover the inner one
0705               if (modIdInner != 0) {
0706                 LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0707                 if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0708                   TMs.clear();
0709                 if (::isDoubleSided(modIdInner, tTopo)) {
0710                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator, 1));
0711                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator, 2));
0712                 } else
0713                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator));
0714                 missingHitAdded = true;
0715                 hitRecoveryCounters[misLayer] += 1;
0716               }
0717             }
0718             if (innerModInactive) {  //If inner missing hit is in inactive module, recover the outer one
0719               if (modIdOuter != 0) {
0720                 LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0721                 if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0722                   TMs.clear();
0723                 if (::isDoubleSided(modIdOuter, tTopo)) {
0724                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator, 1));
0725                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator, 2));
0726                 } else
0727                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator));
0728                 missingHitAdded = true;
0729                 hitRecoveryCounters[previousMisLayer] += 1;
0730               }
0731             }
0732           }
0733         }
0734 
0735         prev_TKlayers = TKlayers;
0736         if (!useFirstMeas_ && isFirstMeas && !missingHitAdded)
0737           continue;
0738         if (!useLastMeas_ && isLastMeas)
0739           continue;
0740         bool hitsWithBias = false;
0741         for (auto ilayer : missedLayers) {
0742           if (ilayer < TKlayers)
0743             hitsWithBias = true;
0744         }
0745         if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing && !missingHitAdded &&
0746             hitsWithBias && !useAllHitsFromTracksWithMissingHits_) {
0747           continue;
0748         }
0749         //////////////////////////////////////////////
0750         //Now check for tracks at TOB6 and TEC9
0751 
0752         // to make sure we only propagate on the last TOB5 hit check the next entry isn't also in TOB5
0753         // to avoid bias, make sure the TOB5 hit is valid (an invalid hit on TOB5 could only exist with a valid hit on TOB6)
0754 
0755         bool isValid = theInHit->isValid();
0756         bool isLast = (itm == (TMeas.end() - 1));
0757         bool isLastTOB5 = true;
0758         if (!isLast) {
0759           if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9)
0760             isLastTOB5 = false;
0761           else
0762             isLastTOB5 = true;
0763           --itm;
0764         }
0765 
0766         if (TKlayers == 9 && isValid && isLastTOB5) {
0767           //      if ( TKlayers==9 && itm==TMeas.rbegin()) {
0768           //      if ( TKlayers==9 && (itm==TMeas.back()) ) {     // to check for only the last entry in the trajectory for propagation
0769           std::vector<BarrelDetLayer const*> barrelTOBLayers =
0770               measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0771           const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size() - 1];
0772           const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0773           const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
0774           auto tmp = layerMeasurements.measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
0775 
0776           if (!tmp.empty()) {
0777             LogDebug("SiStripHitEfficiency:HitEff") << "size of TM from propagation = " << tmp.size() << endl;
0778 
0779             // take the last of the TMs, which is always an invalid hit
0780             // if no detId is available, ie detId==0, then no compatible layer was crossed
0781             // otherwise, use that TM for the efficiency measurement
0782             TrajectoryMeasurement tob6TM(tmp.back());
0783             const auto& tob6Hit = tob6TM.recHit();
0784 
0785             if (tob6Hit->geographicalId().rawId() != 0) {
0786               LogDebug("SiStripHitEfficiency:HitEff") << "tob6 hit actually being added to TM vector" << endl;
0787               TMs.push_back(TrajectoryAtInvalidHit(tob6TM, tTopo, tkgeom, propagator));
0788             }
0789           }
0790         }
0791 
0792         bool isLastTEC8 = true;
0793         if (!isLast) {
0794           if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 21)
0795             isLastTEC8 = false;
0796           else
0797             isLastTEC8 = true;
0798           --itm;
0799         }
0800 
0801         if (TKlayers == 21 && isValid && isLastTEC8) {
0802           std::vector<const ForwardDetLayer*> posTecLayers =
0803               measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0804           const DetLayer* tec9pos = posTecLayers[posTecLayers.size() - 1];
0805           std::vector<const ForwardDetLayer*> negTecLayers =
0806               measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0807           const DetLayer* tec9neg = negTecLayers[negTecLayers.size() - 1];
0808           const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0809           const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
0810 
0811           // check if track on positive or negative z
0812           if (!(iidd == StripSubdetector::TEC))
0813             LogDebug("SiStripHitEfficiency:HitEff") << "there is a problem with TEC 9 extrapolation" << endl;
0814 
0815           //LOGPRINT << " tec9 id = " << iidd << " and side = " << tTopo->tecSide(iidd) << endl;
0816           vector<TrajectoryMeasurement> tmp;
0817           if (tTopo->tecSide(iidd) == 1) {
0818             tmp = layerMeasurements.measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
0819             //LOGPRINT << "on negative side" << endl;
0820           }
0821           if (tTopo->tecSide(iidd) == 2) {
0822             tmp = layerMeasurements.measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
0823             //LOGPRINT << "on positive side" << endl;
0824           }
0825 
0826           if (!tmp.empty()) {
0827             // take the last of the TMs, which is always an invalid hit
0828             // if no detId is available, ie detId==0, then no compatible layer was crossed
0829             // otherwise, use that TM for the efficiency measurement
0830             TrajectoryMeasurement tec9TM(tmp.back());
0831             const auto& tec9Hit = tec9TM.recHit();
0832 
0833             unsigned int tec9id = tec9Hit->geographicalId().rawId();
0834             LogDebug("SiStripHitEfficiency:HitEff")
0835                 << "tec9id = " << tec9id << " is Double sided = " << ::isDoubleSided(tec9id, tTopo)
0836                 << "  and 0x3 = " << (tec9id & 0x3) << endl;
0837 
0838             if (tec9Hit->geographicalId().rawId() != 0) {
0839               LogDebug("SiStripHitEfficiency:HitEff") << "tec9 hit actually being added to TM vector" << endl;
0840               // in tec the hit can be single or doubled sided. whenever the invalid hit at the end of vector of TMs is
0841               // double sided it is always on the matched surface, so we need to split it into the true sensor surfaces
0842               if (::isDoubleSided(tec9id, tTopo)) {
0843                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 1));
0844                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 2));
0845               } else
0846                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator));
0847             }
0848           }  //else LOGPRINT << "tec9 tmp empty" << endl;
0849         }
0850         hitTotalCounters[TKlayers] += 1;
0851 
0852         ////////////////////////////////////////////////////////
0853 
0854         // Modules Constraints
0855 
0856         for (std::vector<TrajectoryAtInvalidHit>::const_iterator TM = TMs.begin(); TM != TMs.end(); ++TM) {
0857           // --> Get trajectory from combinatedState
0858           iidd = TM->monodet_id();
0859           LogDebug("SiStripHitEfficiency:HitEff") << "setting iidd = " << iidd << " before checking efficiency and ";
0860 
0861           xloc = TM->localX();
0862           yloc = TM->localY();
0863 
0864           angleX = atan(TM->localDxDz());
0865           angleY = atan(TM->localDyDz());
0866 
0867           TrajLocErrX = 0.0;
0868           TrajLocErrY = 0.0;
0869 
0870           xglob = TM->globalX();
0871           yglob = TM->globalY();
0872           zglob = TM->globalZ();
0873           xErr = TM->localErrorX();
0874           yErr = TM->localErrorY();
0875 
0876           TrajGlbX = 0.0;
0877           TrajGlbY = 0.0;
0878           TrajGlbZ = 0.0;
0879           withinAcceptance = TM->withinAcceptance();
0880 
0881           trajHitValid = TM->validHit();
0882           int TrajStrip = -1;
0883 
0884           // reget layer from iidd here, to account for TOB 6 and TEC 9 TKlayers being off
0885           TKlayers = ::checkLayer(iidd, tTopo);
0886 
0887           if ((layers == TKlayers) || (layers == 0)) {  // Look at the layer not used to reconstruct the track
0888             whatlayer = TKlayers;
0889             LogDebug("SiStripHitEfficiency:HitEff") << "Looking at layer under study" << endl;
0890             ModIsBad = 2;
0891             Id = 0;
0892             SiStripQualBad = 0;
0893             run = 0;
0894             event = 0;
0895             TrajLocX = 0.0;
0896             TrajLocY = 0.0;
0897             TrajLocAngleX = -999.0;
0898             TrajLocAngleY = -999.0;
0899             ResX = 0.0;
0900             ResXSig = 0.0;
0901             ClusterLocX = 0.0;
0902             ClusterLocY = 0.0;
0903             ClusterLocErrX = 0.0;
0904             ClusterLocErrY = 0.0;
0905             ClusterStoN = 0.0;
0906             bunchx = 0;
0907             commonMode = -100;
0908 
0909             // RPhi RecHit Efficiency
0910 
0911             if (!input.empty()) {
0912               LogDebug("SiStripHitEfficiency:HitEff") << "Checking clusters with size = " << input.size() << endl;
0913               int nClusters = 0;
0914               std::vector<std::vector<float> >
0915                   VCluster_info;  //fill with X residual, X residual pull, local X, sig(X), local Y, sig(Y), StoN
0916               for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end();
0917                    DSViter++) {
0918                 // DSViter is a vector of SiStripClusters located on a single module
0919                 //if (DEBUG)      LOGPRINT << "the ID from the DSViter = " << DSViter->id() << endl;
0920                 unsigned int ClusterId = DSViter->id();
0921                 if (ClusterId == iidd) {
0922                   LogDebug("SiStripHitEfficiency:HitEff")
0923                       << "found  (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
0924                   DetId ClusterDetId(ClusterId);
0925                   const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
0926                   const StripTopology& Topo = stripdet->specificTopology();
0927 
0928                   float hbedge = 0.0;
0929                   float htedge = 0.0;
0930                   float hapoth = 0.0;
0931                   float uylfac = 0.0;
0932                   float uxlden = 0.0;
0933                   if (TKlayers >= 11) {
0934                     const BoundPlane& plane = stripdet->surface();
0935                     const TrapezoidalPlaneBounds* trapezoidalBounds(
0936                         dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0937                     std::array<const float, 4> const& parameterTrap =
0938                         (*trapezoidalBounds).parameters();  // el bueno aqui
0939                     hbedge = parameterTrap[0];
0940                     htedge = parameterTrap[1];
0941                     hapoth = parameterTrap[3];
0942                     uylfac = (htedge - hbedge) / (htedge + hbedge) / hapoth;
0943                     uxlden = 1 + yloc * uylfac;
0944                   }
0945 
0946                   // Need to know position of trajectory in strip number for selecting the right APV later
0947                   if (TrajStrip == -1) {
0948                     int nstrips = Topo.nstrips();
0949                     float pitch = stripdet->surface().bounds().width() / nstrips;
0950                     TrajStrip = xloc / pitch + nstrips / 2.0;
0951                     // Need additionnal corrections for endcap
0952                     if (TKlayers >= 11) {
0953                       float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0954                                                           hapoth);  // radialy extrapolated x loc position at middle
0955                       TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
0956                     }
0957                     //LOGPRINT<<" Layer "<<TKlayers<<" TrajStrip: "<<nstrips<<" "<<pitch<<" "<<TrajStrip<<endl;
0958                   }
0959 
0960                   for (edmNew::DetSet<SiStripCluster>::const_iterator iter = DSViter->begin(); iter != DSViter->end();
0961                        ++iter) {
0962                     //iter is a single SiStripCluster
0963                     StripClusterParameterEstimator::LocalValues parameters = stripcpe.localParameters(*iter, *stripdet);
0964                     float res = (parameters.first.x() - xloc);
0965                     float sigma = ::checkConsistency(parameters, xloc, xErr);
0966                     // The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
0967                     // you need a TransientTrackingRecHit instead of the cluster
0968                     //theEstimator=       new Chi2MeasurementEstimator(30);
0969                     //const Chi2MeasurementEstimator *theEstimator(100);
0970                     //theEstimator->estimate(TM->tsos(), TransientTrackingRecHit);
0971 
0972                     if (TKlayers >= 11) {
0973                       res = parameters.first.x() - xloc / uxlden;  // radialy extrapolated x loc position at middle
0974                       sigma = abs(res) /
0975                               sqrt(parameters.second.xx() + xErr * xErr / uxlden / uxlden +
0976                                    yErr * yErr * xloc * xloc * uylfac * uylfac / uxlden / uxlden / uxlden / uxlden);
0977                     }
0978 
0979                     siStripClusterInfo_.setCluster(*iter, ClusterId);
0980                     // signal to noise from SiStripClusterInfo not working in 225. I'll fix this after the interface
0981                     // redesign in 300 -ku
0982                     //float cluster_info[7] = {res, sigma, parameters.first.x(), sqrt(parameters.second.xx()), parameters.first.y(), sqrt(parameters.second.yy()), signal_to_noise};
0983                     std::vector<float> cluster_info;
0984                     cluster_info.push_back(res);
0985                     cluster_info.push_back(sigma);
0986                     cluster_info.push_back(parameters.first.x());
0987                     cluster_info.push_back(sqrt(parameters.second.xx()));
0988                     cluster_info.push_back(parameters.first.y());
0989                     cluster_info.push_back(sqrt(parameters.second.yy()));
0990                     cluster_info.push_back(siStripClusterInfo_.signalOverNoise());
0991                     VCluster_info.push_back(cluster_info);
0992                     nClusters++;
0993                     LogDebug("SiStripHitEfficiency:HitEff") << "Have ID match. residual = " << VCluster_info.back()[0]
0994                                                             << "  res sigma = " << VCluster_info.back()[1] << endl;
0995                     LogDebug("SiStripHitEfficiency:HitEff")
0996                         << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
0997                     LogDebug("SiStripHitEfficiency:HitEff")
0998                         << "hit position = " << parameters.first.x() << "  hit error = " << sqrt(parameters.second.xx())
0999                         << "  trajectory position = " << xloc << "  traj error = " << xErr << endl;
1000                   }
1001                 }
1002               }
1003               float FinalResSig = 1000.0;
1004               float FinalCluster[7] = {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1005               if (nClusters > 0) {
1006                 LogDebug("SiStripHitEfficiency:HitEff") << "found clusters > 0" << endl;
1007                 if (nClusters > 1) {
1008                   //get the smallest one
1009                   vector<vector<float> >::iterator ires;
1010                   for (ires = VCluster_info.begin(); ires != VCluster_info.end(); ires++) {
1011                     if (abs((*ires)[1]) < abs(FinalResSig)) {
1012                       FinalResSig = (*ires)[1];
1013                       for (unsigned int i = 0; i < ires->size(); i++) {
1014                         LogDebug("SiStripHitEfficiency:HitEff")
1015                             << "filling final cluster. i = " << i << " before fill FinalCluster[i]=" << FinalCluster[i]
1016                             << " and (*ires)[i] =" << (*ires)[i] << endl;
1017                         FinalCluster[i] = (*ires)[i];
1018                         LogDebug("SiStripHitEfficiency:HitEff")
1019                             << "filling final cluster. i = " << i << " after fill FinalCluster[i]=" << FinalCluster[i]
1020                             << " and (*ires)[i] =" << (*ires)[i] << endl;
1021                       }
1022                     }
1023                     LogDebug("SiStripHitEfficiency:HitEff")
1024                         << "iresidual = " << (*ires)[0] << "  isigma = " << (*ires)[1]
1025                         << "  and FinalRes = " << FinalCluster[0] << endl;
1026                   }
1027                 } else {
1028                   FinalResSig = VCluster_info.at(0)[1];
1029                   for (unsigned int i = 0; i < VCluster_info.at(0).size(); i++) {
1030                     FinalCluster[i] = VCluster_info.at(0)[i];
1031                   }
1032                 }
1033                 VCluster_info.clear();
1034               }
1035 
1036               LogDebug("SiStripHitEfficiency:HitEff")
1037                   << "Final residual in X = " << FinalCluster[0] << "+-" << (FinalCluster[0] / FinalResSig) << endl;
1038               LogDebug("SiStripHitEfficiency:HitEff") << "Checking location of trajectory: abs(yloc) = " << abs(yloc)
1039                                                       << "  abs(xloc) = " << abs(xloc) << endl;
1040               LogDebug("SiStripHitEfficiency:HitEff")
1041                   << "Checking location of cluster hit: yloc = " << FinalCluster[4] << "+-" << FinalCluster[5]
1042                   << "  xloc = " << FinalCluster[2] << "+-" << FinalCluster[3] << endl;
1043               LogDebug("SiStripHitEfficiency:HitEff") << "Final cluster signal to noise = " << FinalCluster[6] << endl;
1044 
1045               float exclusionWidth = 0.4;
1046               float TOBexclusion = 0.0;
1047               float TECexRing5 = -0.89;
1048               float TECexRing6 = -0.56;
1049               float TECexRing7 = 0.60;
1050               //Added by Chris Edelmaier to do TEC bonding exclusion
1051               int subdetector = ((iidd >> 25) & 0x7);
1052               int ringnumber = ((iidd >> 5) & 0x7);
1053 
1054               //New TOB and TEC bonding region exclusion zone
1055               if ((TKlayers >= 5 && TKlayers < 11) ||
1056                   ((subdetector == 6) && ((ringnumber >= 5) && (ringnumber <= 7)))) {
1057                 //There are only 2 cases that we need to exclude for
1058                 float highzone = 0.0;
1059                 float lowzone = 0.0;
1060                 float higherr = yloc + 5.0 * yErr;
1061                 float lowerr = yloc - 5.0 * yErr;
1062                 if (TKlayers >= 5 && TKlayers < 11) {
1063                   //TOB zone
1064                   highzone = TOBexclusion + exclusionWidth;
1065                   lowzone = TOBexclusion - exclusionWidth;
1066                 } else if (ringnumber == 5) {
1067                   //TEC ring 5
1068                   highzone = TECexRing5 + exclusionWidth;
1069                   lowzone = TECexRing5 - exclusionWidth;
1070                 } else if (ringnumber == 6) {
1071                   //TEC ring 6
1072                   highzone = TECexRing6 + exclusionWidth;
1073                   lowzone = TECexRing6 - exclusionWidth;
1074                 } else if (ringnumber == 7) {
1075                   //TEC ring 7
1076                   highzone = TECexRing7 + exclusionWidth;
1077                   lowzone = TECexRing7 - exclusionWidth;
1078                 }
1079                 //Now that we have our exclusion region, we just have to properly identify it
1080                 if ((highzone <= higherr) && (highzone >= lowerr))
1081                   withinAcceptance = false;
1082                 if ((lowzone >= lowerr) && (lowzone <= higherr))
1083                   withinAcceptance = false;
1084                 if ((higherr <= highzone) && (higherr >= lowzone))
1085                   withinAcceptance = false;
1086                 if ((lowerr >= lowzone) && (lowerr <= highzone))
1087                   withinAcceptance = false;
1088               }
1089 
1090               // fill ntuple varibles
1091               //get global position from module id number iidd
1092               TrajGlbX = xglob;
1093               TrajGlbY = yglob;
1094               TrajGlbZ = zglob;
1095 
1096               TrajLocErrX = xErr;
1097               TrajLocErrY = yErr;
1098 
1099               Id = iidd;
1100               run = run_nr;
1101               event = ev_nr;
1102               bunchx = bunch_nr;
1103               //if ( SiStripQuality_->IsModuleBad(iidd) ) {
1104               if (SiStripQuality_->getBadApvs(iidd) != 0) {
1105                 SiStripQualBad = 1;
1106                 LogDebug("SiStripHitEfficiency:HitEff") << "strip is bad from SiStripQuality" << endl;
1107               } else {
1108                 SiStripQualBad = 0;
1109                 LogDebug("SiStripHitEfficiency:HitEff") << "strip is good from SiStripQuality" << endl;
1110               }
1111 
1112               //check for FED-detected errors and include those in SiStripQualBad
1113               for (unsigned int ii = 0; ii < fedErrorIds.size(); ii++) {
1114                 if (iidd == fedErrorIds[ii].rawId())
1115                   SiStripQualBad = 1;
1116               }
1117 
1118               TrajLocX = xloc;
1119               TrajLocY = yloc;
1120               TrajLocAngleX = angleX;
1121               TrajLocAngleY = angleY;
1122               ResX = FinalCluster[0];
1123               ResXSig = FinalResSig;
1124               if (FinalResSig != FinalCluster[1])
1125                 LogDebug("SiStripHitEfficiency:HitEff")
1126                     << "Problem with best cluster selection because FinalResSig = " << FinalResSig
1127                     << " and FinalCluster[1] = " << FinalCluster[1] << endl;
1128               ClusterLocX = FinalCluster[2];
1129               ClusterLocY = FinalCluster[4];
1130               ClusterLocErrX = FinalCluster[3];
1131               ClusterLocErrY = FinalCluster[5];
1132               ClusterStoN = FinalCluster[6];
1133 
1134               // CM of APV crossed by traj
1135               if (addCommonMode_)
1136                 if (commonModeDigis.isValid() && TrajStrip >= 0 && TrajStrip <= 768) {
1137                   edm::DetSetVector<SiStripRawDigi>::const_iterator digiframe = commonModeDigis->find(iidd);
1138                   if (digiframe != commonModeDigis->end())
1139                     if ((unsigned)TrajStrip / 128 < digiframe->data.size())
1140                       commonMode = digiframe->data.at(TrajStrip / 128).adc();
1141                 }
1142 
1143               LogDebug("SiStripHitEfficiency:HitEff") << "before check good" << endl;
1144 
1145               if (FinalResSig < 999.0) {  //could make requirement on track/hit consistency, but for
1146                 //now take anything with a hit on the module
1147                 LogDebug("SiStripHitEfficiency:HitEff")
1148                     << "hit being counted as good " << FinalCluster[0] << " FinalRecHit " << iidd << "   TKlayers  "
1149                     << TKlayers << " xloc " << xloc << " yloc  " << yloc << " module " << iidd
1150                     << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1) << "/"
1151                     << ((iidd & 0x3) == 2) << endl;
1152                 ModIsBad = 0;
1153                 traj->Fill();
1154               } else {
1155                 LogDebug("SiStripHitEfficiency:HitEff")
1156                     << "hit being counted as bad   ######### Invalid RPhi FinalResX " << FinalCluster[0]
1157                     << " FinalRecHit " << iidd << "   TKlayers  " << TKlayers << " xloc " << xloc << " yloc  " << yloc
1158                     << " module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
1159                     << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
1160                 ModIsBad = 1;
1161                 traj->Fill();
1162 
1163                 LogDebug("SiStripHitEfficiency:HitEff") << " RPhi Error " << sqrt(xErr * xErr + yErr * yErr)
1164                                                         << " ErrorX " << xErr << " yErr " << yErr << endl;
1165               }
1166               LogDebug("SiStripHitEfficiency:HitEff") << "after good location check" << endl;
1167             }
1168             LogDebug("SiStripHitEfficiency:HitEff") << "after list of clusters" << endl;
1169           }
1170           LogDebug("SiStripHitEfficiency:HitEff") << "After layers=TKLayers if" << endl;
1171         }
1172         LogDebug("SiStripHitEfficiency:HitEff") << "After looping over TrajAtValidHit list" << endl;
1173       }
1174       LogDebug("SiStripHitEfficiency:HitEff") << "end TMeasurement loop" << endl;
1175     }
1176     LogDebug("SiStripHitEfficiency:HitEff") << "end of trajectories loop" << endl;
1177   }
1178 }
1179 
1180 void HitEff::endJob() {
1181   traj->GetDirectory()->cd();
1182   traj->Write();
1183 
1184   LogDebug("SiStripHitEfficiency:HitEff") << " Events Analysed             " << events << endl;
1185   LogDebug("SiStripHitEfficiency:HitEff") << " Number Of Tracked events    " << EventTrackCKF << endl;
1186 
1187   if (doMissingHitsRecovery_) {
1188     float totTIB = 0.0;
1189     float totTOB = 0.0;
1190     float totTID = 0.0;
1191     float totTEC = 0.0;
1192 
1193     float totTIBrepro = 0.0;
1194     float totTOBrepro = 0.0;
1195     float totTIDrepro = 0.0;
1196     float totTECrepro = 0.0;
1197 
1198     edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TIB :";
1199     for (int i = 0; i <= k_LayersAtTIBEnd; i++) {
1200       edm::LogInfo("SiStripHitEfficiency:HitEff")
1201           << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1202           << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1203       totTIB += missHitPerLayer[i];
1204       edm::LogInfo("SiStripHitEfficiency:HitEff")
1205           << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1206           << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1207           << " % of missing hit";
1208       totTIBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1209     }
1210     edm::LogInfo("SiStripHitEfficiency:HitEff")
1211         << "TOTAL % of missing hits within TIB :" << (totTIB * 1.0 / totalNbHits) * 100 << "%";
1212     edm::LogInfo("SiStripHitEfficiency:HitEff")
1213         << "AFTER repropagation :" << (totTIBrepro * 1.0 / totalNbHits) * 100 << "%";
1214 
1215     edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TOB :";
1216     for (int i = k_LayersAtTIBEnd + 1; i <= k_LayersAtTOBEnd; i++) {
1217       edm::LogInfo("SiStripHitEfficiency:HitEff")
1218           << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1219           << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1220       totTOB += missHitPerLayer[i];
1221       edm::LogInfo("SiStripHitEfficiency:HitEff")
1222           << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1223           << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1224           << " % of missing hit";
1225       totTOBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1226     }
1227     edm::LogInfo("SiStripHitEfficiency:HitEff")
1228         << "TOTAL % of missing hits within TOB :" << (totTOB * 1.0 / totalNbHits) * 100 << "%";
1229     edm::LogInfo("SiStripHitEfficiency:HitEff")
1230         << "AFTER repropagation :" << (totTOBrepro * 1.0 / totalNbHits) * 100 << "%";
1231 
1232     edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TID :";
1233     for (int i = k_LayersAtTOBEnd + 1; i <= k_LayersAtTIDEnd; i++) {
1234       edm::LogInfo("SiStripHitEfficiency:HitEff")
1235           << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1236           << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1237       totTID += missHitPerLayer[i];
1238       edm::LogInfo("SiStripHitEfficiency:HitEff")
1239           << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1240           << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1241           << " % of missing hit";
1242       totTIDrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1243     }
1244     edm::LogInfo("SiStripHitEfficiency:HitEff")
1245         << "TOTAL % of missing hits within TID :" << (totTID * 1.0 / totalNbHits) * 100 << "%";
1246     edm::LogInfo("SiStripHitEfficiency:HitEff")
1247         << "AFTER repropagation :" << (totTIDrepro * 1.0 / totalNbHits) * 100 << "%";
1248 
1249     edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TEC :";
1250     for (int i = k_LayersAtTIDEnd + 1; i < k_END_OF_LAYERS; i++) {
1251       edm::LogInfo("SiStripHitEfficiency:HitEff")
1252           << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1253           << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1254       totTEC += missHitPerLayer[i];
1255       edm::LogInfo("SiStripHitEfficiency:HitEff")
1256           << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1257           << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1258           << " % of missing hit";
1259       totTECrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1260     }
1261     edm::LogInfo("SiStripHitEfficiency:HitEff")
1262         << "TOTAL % of missing hits within TEC :" << (totTEC * 1.0 / totalNbHits) * 100 << "%";
1263     edm::LogInfo("SiStripHitEfficiency:HitEff")
1264         << "AFTER repropagation :" << (totTECrepro * 1.0 / totalNbHits) * 100 << "%";
1265 
1266     edm::LogInfo("SiStripHitEfficiency:HitEff") << " Hit recovery summary:";
1267 
1268     for (int ilayer = 0; ilayer < k_END_OF_LAYERS; ilayer++) {
1269       edm::LogInfo("SiStripHitEfficiency:HitEff")
1270           << " layer " << ilayer << ": " << hitRecoveryCounters[ilayer] << " / " << hitTotalCounters[ilayer];
1271     }
1272   }
1273 }
1274 
1275 //define this as a plug-in
1276 DEFINE_FWK_MODULE(HitEff);
[ 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!