SiStripHitEfficiency/plugins/HitEff.cc

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 "Geometry/CommonDetUnit/interface/GeomDet.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0048 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0049 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0050 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0051 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0052 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0053 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0054 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0055 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0056 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0057 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0058 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0059 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0060
0061 // ROOT includes
0062 #include "TMath.h"
0063 #include "TH1F.h"
0064
0065 // custom made printout
0066 #define LOGPRINT edm::LogPrint("SiStripHitEfficiency:HitEff")
0067
0068 //
0069 // constructors and destructor
0070 //
0071
0072 using namespace std;
0073 HitEff::HitEff(const edm::ParameterSet& conf)
0074     : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))),
0075       metaDataToken_(consumes<OnlineLuminosityRecord>(conf.getParameter<edm::InputTag>("metadata"))),
0076       commonModeToken_(mayConsume<edm::DetSetVector<SiStripRawDigi> >(conf.getParameter<edm::InputTag>("commonMode"))),
0077       siStripClusterInfo_(consumesCollector()),
0078       combinatorialTracks_token_(
0079           consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))),
0080       trajectories_token_(consumes<std::vector<Trajectory> >(conf.getParameter<edm::InputTag>("trajectories"))),
0081       trajTrackAsso_token_(consumes<TrajTrackAssociationCollection>(conf.getParameter<edm::InputTag>("trajectories"))),
0082       clusters_token_(
0083           consumes<edmNew::DetSetVector<SiStripCluster> >(conf.getParameter<edm::InputTag>("siStripClusters"))),
0084       digis_token_(consumes<DetIdCollection>(conf.getParameter<edm::InputTag>("siStripDigis"))),
0085       trackerEvent_token_(consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("trackerEvent"))),
0086       topoToken_(esConsumes()),
0087       geomToken_(esConsumes()),
0088       cpeToken_(esConsumes(edm::ESInputTag("", "StripCPEfromTrackAngle"))),
0089       siStripQualityToken_(esConsumes()),
0090       magFieldToken_(esConsumes()),
0091       measurementTkToken_(esConsumes()),
0092       chi2MeasurementEstimatorToken_(esConsumes(edm::ESInputTag("", "Chi2"))),
0093       propagatorToken_(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))),
0094       conf_(conf) {
0095   compSettings = conf_.getUntrackedParameter<int>("CompressionSettings", -1);
0096   layers = conf_.getParameter<int>("Layer");
0097   DEBUG = conf_.getParameter<bool>("Debug");
0098   addLumi_ = conf_.getUntrackedParameter<bool>("addLumi", false);
0099   addCommonMode_ = conf_.getUntrackedParameter<bool>("addCommonMode", false);
0100   cutOnTracks_ = conf_.getUntrackedParameter<bool>("cutOnTracks", false);
0101   trackMultiplicityCut_ = conf.getUntrackedParameter<unsigned int>("trackMultiplicity", 100);
0102   useFirstMeas_ = conf_.getUntrackedParameter<bool>("useFirstMeas", false);
0103   useLastMeas_ = conf_.getUntrackedParameter<bool>("useLastMeas", false);
0104   useAllHitsFromTracksWithMissingHits_ =
0105       conf_.getUntrackedParameter<bool>("useAllHitsFromTracksWithMissingHits", false);
0106 }
0107
0108 void HitEff::beginJob() {
0109   edm::Service<TFileService> fs;
0110   if (compSettings > 0) {
0111     edm::LogInfo("SiStripHitEfficiency:HitEff") << "the compressions settings are:" << compSettings << std::endl;
0112     fs->file().SetCompressionSettings(compSettings);
0113   }
0114
0115   traj = fs->make<TTree>("traj", "tree of trajectory positions");
0116 #ifdef ExtendedCALIBTree
0117   traj->Branch("timeDT", &timeDT, "timeDT/F");
0118   traj->Branch("timeDTErr", &timeDTErr, "timeDTErr/F");
0119   traj->Branch("timeDTDOF", &timeDTDOF, "timeDTDOF/I");
0120   traj->Branch("timeECAL", &timeECAL, "timeECAL/F");
0121   traj->Branch("dedx", &dedx, "dedx/F");
0122   traj->Branch("dedxNOM", &dedxNOM, "dedxNOM/I");
0123   traj->Branch("nLostHits", &nLostHits, "nLostHits/I");
0124   traj->Branch("chi2", &chi2, "chi2/F");
0125   traj->Branch("p", &p, "p/F");
0126 #endif
0127   traj->Branch("TrajGlbX", &TrajGlbX, "TrajGlbX/F");
0128   traj->Branch("TrajGlbY", &TrajGlbY, "TrajGlbY/F");
0129   traj->Branch("TrajGlbZ", &TrajGlbZ, "TrajGlbZ/F");
0130   traj->Branch("TrajLocX", &TrajLocX, "TrajLocX/F");
0131   traj->Branch("TrajLocY", &TrajLocY, "TrajLocY/F");
0132   traj->Branch("TrajLocAngleX", &TrajLocAngleX, "TrajLocAngleX/F");
0133   traj->Branch("TrajLocAngleY", &TrajLocAngleY, "TrajLocAngleY/F");
0134   traj->Branch("TrajLocErrX", &TrajLocErrX, "TrajLocErrX/F");
0135   traj->Branch("TrajLocErrY", &TrajLocErrY, "TrajLocErrY/F");
0136   traj->Branch("ClusterLocX", &ClusterLocX, "ClusterLocX/F");
0137   traj->Branch("ClusterLocY", &ClusterLocY, "ClusterLocY/F");
0138   traj->Branch("ClusterLocErrX", &ClusterLocErrX, "ClusterLocErrX/F");
0139   traj->Branch("ClusterLocErrY", &ClusterLocErrY, "ClusterLocErrY/F");
0140   traj->Branch("ClusterStoN", &ClusterStoN, "ClusterStoN/F");
0141   traj->Branch("ResX", &ResX, "ResX/F");
0142   traj->Branch("ResXSig", &ResXSig, "ResXSig/F");
0143   traj->Branch("ModIsBad", &ModIsBad, "ModIsBad/i");
0144   traj->Branch("SiStripQualBad", &SiStripQualBad, "SiStripQualBad/i");
0145   traj->Branch("withinAcceptance", &withinAcceptance, "withinAcceptance/O");
0146   traj->Branch("nHits", &nHits, "nHits/I");
0147   traj->Branch("pT", &pT, "pT/F");
0148   traj->Branch("highPurity", &highPurity, "highPurity/O");
0149   traj->Branch("trajHitValid", &trajHitValid, "trajHitValid/i");
0150   traj->Branch("Id", &Id, "Id/i");
0151   traj->Branch("run", &run, "run/i");
0152   traj->Branch("event", &event, "event/i");
0153   traj->Branch("layer", &whatlayer, "layer/i");
0154   traj->Branch("tquality", &tquality, "tquality/I");
0155   traj->Branch("bunchx", &bunchx, "bunchx/I");
0156   if (addLumi_) {
0157     traj->Branch("instLumi", &instLumi, "instLumi/F");
0158     traj->Branch("PU", &PU, "PU/F");
0159   }
0160   if (addCommonMode_)
0161     traj->Branch("commonMode", &commonMode, "commonMode/F");
0162
0163   events = 0;
0164   EventTrackCKF = 0;
0165 }
0166
0167 void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
0168   //Retrieve tracker topology from geometry
0169   const TrackerTopology* tTopo = &es.getData(topoToken_);
0170   siStripClusterInfo_.initEvent(es);
0171
0172   //  bool DEBUG = false;
0173
0174   LogDebug("SiStripHitEfficiency:HitEff") << "beginning analyze from HitEff" << endl;
0175
0176   using namespace edm;
0177   using namespace reco;
0178   // Step A: Get Inputs
0179
0180   int run_nr = e.id().run();
0181   int ev_nr = e.id().event();
0182   int bunch_nr = e.bunchCrossing();
0183
0184   // Luminosity informations
0185   edm::Handle<LumiScalersCollection> lumiScalers = e.getHandle(scalerToken_);
0186   edm::Handle<OnlineLuminosityRecord> metaData = e.getHandle(metaDataToken_);
0187
0188   instLumi = 0;
0189   PU = 0;
0190   if (addLumi_) {
0191     if (lumiScalers.isValid() && !lumiScalers->empty()) {
0192       if (lumiScalers->begin() != lumiScalers->end()) {
0193         instLumi = lumiScalers->begin()->instantLumi();
0194         PU = lumiScalers->begin()->pileup();
0195       }
0196     } else if (metaData.isValid()) {
0197       instLumi = metaData->instLumi();
0198       PU = metaData->avgPileUp();
0199     } else {
0200       edm::LogWarning("SiStripHitEfficiencyWorker") << "could not find a source for the Luminosity and PU";
0201     }
0202   }
0203
0204   // CM
0205   edm::Handle<edm::DetSetVector<SiStripRawDigi> > commonModeDigis;
0206   if (addCommonMode_)
0207     e.getByToken(commonModeToken_, commonModeDigis);
0208
0209   //CombinatoriaTrack
0210   edm::Handle<reco::TrackCollection> trackCollectionCKF;
0211   //edm::InputTag TkTagCKF = conf_.getParameter<edm::InputTag>("combinatorialTracks");
0212   e.getByToken(combinatorialTracks_token_, trackCollectionCKF);
0213
0214   edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
0215   //edm::InputTag TkTrajCKF = conf_.getParameter<edm::InputTag>("trajectories");
0216   e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0217
0218   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0219   e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0220
0221   // Clusters
0222   // get the SiStripClusters from the event
0223   edm::Handle<edmNew::DetSetVector<SiStripCluster> > theClusters;
0224   //e.getByLabel("siStripClusters", theClusters);
0225   e.getByToken(clusters_token_, theClusters);
0226
0227   //get tracker geometry
0228   edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_);
0229   const TrackerGeometry* tkgeom = &(*tracker);
0230
0231   //get Cluster Parameter Estimator
0232   //std::string cpe = conf_.getParameter<std::string>("StripCPE");
0233   edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_);
0234   const StripClusterParameterEstimator& stripcpe(*parameterestimator);
0235
0236   // get the SiStripQuality records
0237   edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_);
0238
0239   const MagneticField* magField_ = &es.getData(magFieldToken_);
0240
0241   // get the list of module IDs with FED-detected errors
0242   edm::Handle<DetIdCollection> fedErrorIds;
0243   //e.getByLabel("siStripDigis", fedErrorIds );
0244   e.getByToken(digis_token_, fedErrorIds);
0245
0246   ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTkToken_);
0247
0248   edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
0249   //e.getByLabel("MeasurementTrackerEvent", measurementTrackerEvent);
0250   e.getByToken(trackerEvent_token_, measurementTrackerEvent);
0251
0252   const MeasurementEstimator* estimator = &es.getData(chi2MeasurementEstimatorToken_);
0253   const Propagator* thePropagator = &es.getData(propagatorToken_);
0254
0255   events++;
0256
0257   // *************** SiStripCluster Collection
0258   const edmNew::DetSetVector<SiStripCluster>& input = *theClusters;
0259
0260   //go through clusters to write out global position of good clusters for the layer understudy for comparison
0261   // Loop through clusters just to print out locations
0262   // Commented out to avoid discussion, should really be deleted.
0263   /*
0264   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end(); DSViter++) {
0265     // DSViter is a vector of SiStripClusters located on a single module
0266     unsigned int ClusterId = DSViter->id();
0267     DetId ClusterDetId(ClusterId);
0268     const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
0269
0270     edmNew::DetSet<SiStripCluster>::const_iterator begin=DSViter->begin();
0271     edmNew::DetSet<SiStripCluster>::const_iterator end  =DSViter->end();
0272     for(edmNew::DetSet<SiStripCluster>::const_iterator iter=begin;iter!=end;++iter) {
0273       //iter is a single SiStripCluster
0274       StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
0275
0276       const Surface* surface;
0277       surface = &(tracker->idToDet(ClusterDetId)->surface());
0278       LocalPoint lp = parameters.first;
0279       GlobalPoint gp = surface->toGlobal(lp);
0280       unsigned int layer = ::checkLayer(ClusterId, tTopo);
0281             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;
0282     }
0283   }
0284   */
0285
0286   // Tracking
0287   const reco::TrackCollection* tracksCKF = trackCollectionCKF.product();
0288   LogDebug("SiStripHitEfficiency:HitEff") << "number ckf tracks found = " << tracksCKF->size() << endl;
0289   //if (tracksCKF->size() == 1 ){
0290   if (!tracksCKF->empty()) {
0291     if (cutOnTracks_ && (tracksCKF->size() >= trackMultiplicityCut_))
0292       return;
0293     if (cutOnTracks_)
0294       LogDebug("SiStripHitEfficiency:HitEff")
0295           << "starting checking good event with < " << trackMultiplicityCut_ << " tracks" << endl;
0296
0297     EventTrackCKF++;
0298
0299 #ifdef ExtendedCALIBTree
0300     //get dEdx info if available
0301     Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
0302     if (e.getByLabel("dedxMedianCTF", dEdxUncalibHandle)) {
0303       const ValueMap<DeDxData> dEdxTrackUncalib = *dEdxUncalibHandle.product();
0304
0305       reco::TrackRef itTrack = reco::TrackRef(trackCollectionCKF, 0);
0306       dedx = dEdxTrackUncalib[itTrack].dEdx();
0307       dedxNOM = dEdxTrackUncalib[itTrack].numberOfMeasurements();
0308     } else {
0309       dedx = -999.0;
0310       dedxNOM = -999;
0311     }
0312
0313     //get muon and ecal timing info if available
0314     Handle<MuonCollection> muH;
0315     if (e.getByLabel("muonsWitht0Correction", muH)) {
0316       const MuonCollection& muonsT0 = *muH.product();
0317       if (!muonsT0.empty()) {
0318         MuonTime mt0 = muonsT0[0].time();
0319         timeDT = mt0.timeAtIpInOut;
0320         timeDTErr = mt0.timeAtIpInOutErr;
0321         timeDTDOF = mt0.nDof;
0322
0323         bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
0324         if (hasCaloEnergyInfo)
0325           timeECAL = muonsT0[0].calEnergy().ecal_time;
0326       }
0327     } else {
0328       timeDT = -999.0;
0329       timeDTErr = -999.0;
0330       timeDTDOF = -999;
0331       timeECAL = -999.0;
0332     }
0333
0334 #endif
0335     // actually should do a loop over all the tracks in the event here
0336
0337     // Looping over traj-track associations to be able to get traj & track informations
0338     for (TrajTrackAssociationCollection::const_iterator it = trajTrackAssociationHandle->begin();
0339          it != trajTrackAssociationHandle->end();
0340          it++) {
0341       edm::Ref<std::vector<Trajectory> > itraj = it->key;
0342       reco::TrackRef itrack = it->val;
0343
0344       // for each track, fill some variables such as number of hits and momentum
0345       nHits = itraj->foundHits();
0346 #ifdef ExtendedCALIBTree
0347       nLostHits = itraj->lostHits();
0348       chi2 = (itraj->chiSquared() / itraj->ndof());
0349       p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
0350 #endif
0351       pT = sqrt((itraj->lastMeasurement().updatedState().globalMomentum().x() *
0352                  itraj->lastMeasurement().updatedState().globalMomentum().x()) +
0353                 (itraj->lastMeasurement().updatedState().globalMomentum().y() *
0354                  itraj->lastMeasurement().updatedState().globalMomentum().y()));
0355
0356       // track quality
0357       highPurity = itrack->quality(reco::TrackBase::TrackQuality::highPurity);
0358
0359       std::vector<TrajectoryMeasurement> TMeas = itraj->measurements();
0360       vector<TrajectoryMeasurement>::iterator itm;
0361       double xloc = 0.;
0362       double yloc = 0.;
0363       double xErr = 0.;
0364       double yErr = 0.;
0365       double angleX = -999.;
0366       double angleY = -999.;
0367       double xglob, yglob, zglob;
0368
0369       // Check whether the trajectory has some missing hits
0370       bool hasMissingHits = false;
0371       for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
0372         auto theHit = (*itm).recHit();
0373         if (theHit->getType() == TrackingRecHit::Type::missing)
0374           hasMissingHits = true;
0375       }
0376
0377       // Loop on each measurement and take it into consideration
0378       //--------------------------------------------------------
0379
0380       for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
0381         auto theInHit = (*itm).recHit();
0382
0383         LogDebug("SiStripHitEfficiency:HitEff") << "theInHit is valid = " << theInHit->isValid() << endl;
0384
0385         unsigned int iidd = theInHit->geographicalId().rawId();
0386
0387         unsigned int TKlayers = ::checkLayer(iidd, tTopo);
0388         LogDebug("SiStripHitEfficiency:HitEff") << "TKlayer from trajectory: " << TKlayers << "  from module = " << iidd
0389                                                 << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0390                                                 << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
0391
0392         // Test first and last points of the trajectory
0393         // the list of measurements starts from outer layers  !!! This could change -> should add a check
0394         bool isFirstMeas = (itm == (TMeas.end() - 1));
0395         bool isLastMeas = (itm == (TMeas.begin()));
0396
0397         if (!useFirstMeas_ && isFirstMeas)
0398           continue;
0399         if (!useLastMeas_ && isLastMeas)
0400           continue;
0401
0402         // In case of missing hit in the track, check whether to use the other hits or not.
0403         if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
0404             !useAllHitsFromTracksWithMissingHits_)
0405           continue;
0406
0407         // If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
0408         if (TKlayers == 10 || TKlayers == 22) {
0409           LogDebug("SiStripHitEfficiency:HitEff") << "skipping original TM for TOB 6 or TEC 9" << endl;
0410           continue;
0411         }
0412
0413         // Make vector of TrajectoryAtInvalidHits to hold the trajectories
0414         std::vector<TrajectoryAtInvalidHit> TMs;
0415
0416         // Make AnalyticalPropagator to use in TAVH constructor
0417         AnalyticalPropagator propagator(magField_, anyDirection);
0418
0419         // for double sided layers check both sensors--if no hit was found on either sensor surface,
0420         // the trajectory measurements only have one invalid hit entry on the matched surface
0421         // so get the TrajectoryAtInvalidHit for both surfaces and include them in the study
0422         if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0423           // do hit eff check twice--once for each sensor
0424           //add a TM for each surface
0425           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0426           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0427         } else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
0428           // if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
0429           // the matched layer should be added to the study as well
0430           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0431           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0432           LogDebug("SiStripHitEfficiency:HitEff") << " found a hit with a missing partner" << endl;
0433         } else {
0434           //only add one TM for the single surface and the other will be added in the next iteration
0435           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator));
0436         }
0437
0438         //////////////////////////////////////////////
0439         //Now check for tracks at TOB6 and TEC9
0440
0441         // to make sure we only propagate on the last TOB5 hit check the next entry isn't also in TOB5
0442         // to avoid bias, make sure the TOB5 hit is valid (an invalid hit on TOB5 could only exist with a valid hit on TOB6)
0443
0444         bool isValid = theInHit->isValid();
0445         bool isLast = (itm == (TMeas.end() - 1));
0446         bool isLastTOB5 = true;
0447         if (!isLast) {
0448           if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9)
0449             isLastTOB5 = false;
0450           else
0451             isLastTOB5 = true;
0452           --itm;
0453         }
0454
0455         if (TKlayers == 9 && isValid && isLastTOB5) {
0456           //      if ( TKlayers==9 && itm==TMeas.rbegin()) {
0457           //      if ( TKlayers==9 && (itm==TMeas.back()) ) {     // to check for only the last entry in the trajectory for propagation
0458           std::vector<BarrelDetLayer const*> barrelTOBLayers =
0459               measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0460           const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size() - 1];
0461           const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0462           const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
0463           auto tmp = layerMeasurements.measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
0464
0465           if (!tmp.empty()) {
0466             LogDebug("SiStripHitEfficiency:HitEff") << "size of TM from propagation = " << tmp.size() << endl;
0467
0468             // take the last of the TMs, which is always an invalid hit
0469             // if no detId is available, ie detId==0, then no compatible layer was crossed
0470             // otherwise, use that TM for the efficiency measurement
0471             TrajectoryMeasurement tob6TM(tmp.back());
0472             const auto& tob6Hit = tob6TM.recHit();
0473
0474             if (tob6Hit->geographicalId().rawId() != 0) {
0475               LogDebug("SiStripHitEfficiency:HitEff") << "tob6 hit actually being added to TM vector" << endl;
0476               TMs.push_back(TrajectoryAtInvalidHit(tob6TM, tTopo, tkgeom, propagator));
0477             }
0478           }
0479         }
0480
0481         bool isLastTEC8 = true;
0482         if (!isLast) {
0483           if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 21)
0484             isLastTEC8 = false;
0485           else
0486             isLastTEC8 = true;
0487           --itm;
0488         }
0489
0490         if (TKlayers == 21 && isValid && isLastTEC8) {
0491           std::vector<const ForwardDetLayer*> posTecLayers =
0492               measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0493           const DetLayer* tec9pos = posTecLayers[posTecLayers.size() - 1];
0494           std::vector<const ForwardDetLayer*> negTecLayers =
0495               measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0496           const DetLayer* tec9neg = negTecLayers[negTecLayers.size() - 1];
0497           const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0498           const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
0499
0500           // check if track on positive or negative z
0501           if (!(iidd == StripSubdetector::TEC))
0502             LogDebug("SiStripHitEfficiency:HitEff") << "there is a problem with TEC 9 extrapolation" << endl;
0503
0504           //LOGPRINT << " tec9 id = " << iidd << " and side = " << tTopo->tecSide(iidd) << endl;
0505           vector<TrajectoryMeasurement> tmp;
0506           if (tTopo->tecSide(iidd) == 1) {
0507             tmp = layerMeasurements.measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
0508             //LOGPRINT << "on negative side" << endl;
0509           }
0510           if (tTopo->tecSide(iidd) == 2) {
0511             tmp = layerMeasurements.measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
0512             //LOGPRINT << "on positive side" << endl;
0513           }
0514
0515           if (!tmp.empty()) {
0516             // take the last of the TMs, which is always an invalid hit
0517             // if no detId is available, ie detId==0, then no compatible layer was crossed
0518             // otherwise, use that TM for the efficiency measurement
0519             TrajectoryMeasurement tec9TM(tmp.back());
0520             const auto& tec9Hit = tec9TM.recHit();
0521
0522             unsigned int tec9id = tec9Hit->geographicalId().rawId();
0523             LogDebug("SiStripHitEfficiency:HitEff")
0524                 << "tec9id = " << tec9id << " is Double sided = " << ::isDoubleSided(tec9id, tTopo)
0525                 << "  and 0x3 = " << (tec9id & 0x3) << endl;
0526
0527             if (tec9Hit->geographicalId().rawId() != 0) {
0528               LogDebug("SiStripHitEfficiency:HitEff") << "tec9 hit actually being added to TM vector" << endl;
0529               // in tec the hit can be single or doubled sided. whenever the invalid hit at the end of vector of TMs is
0530               // double sided it is always on the matched surface, so we need to split it into the true sensor surfaces
0531               if (::isDoubleSided(tec9id, tTopo)) {
0532                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 1));
0533                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 2));
0534               } else
0535                 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator));
0536             }
0537           }  //else LOGPRINT << "tec9 tmp empty" << endl;
0538         }
0539
0540         ////////////////////////////////////////////////////////
0541
0542         // Modules Constraints
0543
0544         for (std::vector<TrajectoryAtInvalidHit>::const_iterator TM = TMs.begin(); TM != TMs.end(); ++TM) {
0545           // --> Get trajectory from combinatedState
0546           iidd = TM->monodet_id();
0547           LogDebug("SiStripHitEfficiency:HitEff") << "setting iidd = " << iidd << " before checking efficiency and ";
0548
0549           xloc = TM->localX();
0550           yloc = TM->localY();
0551
0552           angleX = atan(TM->localDxDz());
0553           angleY = atan(TM->localDyDz());
0554
0555           TrajLocErrX = 0.0;
0556           TrajLocErrY = 0.0;
0557
0558           xglob = TM->globalX();
0559           yglob = TM->globalY();
0560           zglob = TM->globalZ();
0561           xErr = TM->localErrorX();
0562           yErr = TM->localErrorY();
0563
0564           TrajGlbX = 0.0;
0565           TrajGlbY = 0.0;
0566           TrajGlbZ = 0.0;
0567           withinAcceptance = TM->withinAcceptance();
0568
0569           trajHitValid = TM->validHit();
0570           int TrajStrip = -1;
0571
0572           // reget layer from iidd here, to account for TOB 6 and TEC 9 TKlayers being off
0573           TKlayers = ::checkLayer(iidd, tTopo);
0574
0575           if ((layers == TKlayers) || (layers == 0)) {  // Look at the layer not used to reconstruct the track
0576             whatlayer = TKlayers;
0577             LogDebug("SiStripHitEfficiency:HitEff") << "Looking at layer under study" << endl;
0578             ModIsBad = 2;
0579             Id = 0;
0580             SiStripQualBad = 0;
0581             run = 0;
0582             event = 0;
0583             TrajLocX = 0.0;
0584             TrajLocY = 0.0;
0585             TrajLocAngleX = -999.0;
0586             TrajLocAngleY = -999.0;
0587             ResX = 0.0;
0588             ResXSig = 0.0;
0589             ClusterLocX = 0.0;
0590             ClusterLocY = 0.0;
0591             ClusterLocErrX = 0.0;
0592             ClusterLocErrY = 0.0;
0593             ClusterStoN = 0.0;
0594             bunchx = 0;
0595             commonMode = -100;
0596
0597             // RPhi RecHit Efficiency
0598
0599             if (!input.empty()) {
0600               LogDebug("SiStripHitEfficiency:HitEff") << "Checking clusters with size = " << input.size() << endl;
0601               int nClusters = 0;
0602               std::vector<std::vector<float> >
0603                   VCluster_info;  //fill with X residual, X residual pull, local X, sig(X), local Y, sig(Y), StoN
0604               for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end();
0605                    DSViter++) {
0606                 // DSViter is a vector of SiStripClusters located on a single module
0607                 //if (DEBUG)      LOGPRINT << "the ID from the DSViter = " << DSViter->id() << endl;
0608                 unsigned int ClusterId = DSViter->id();
0609                 if (ClusterId == iidd) {
0610                   LogDebug("SiStripHitEfficiency:HitEff")
0611                       << "found  (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
0612                   DetId ClusterDetId(ClusterId);
0613                   const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
0614                   const StripTopology& Topo = stripdet->specificTopology();
0615
0616                   float hbedge = 0.0;
0617                   float htedge = 0.0;
0618                   float hapoth = 0.0;
0619                   float uylfac = 0.0;
0620                   float uxlden = 0.0;
0621                   if (TKlayers >= 11) {
0622                     const BoundPlane& plane = stripdet->surface();
0623                     const TrapezoidalPlaneBounds* trapezoidalBounds(
0624                         dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0625                     std::array<const float, 4> const& parameterTrap =
0626                         (*trapezoidalBounds).parameters();  // el bueno aqui
0627                     hbedge = parameterTrap[0];
0628                     htedge = parameterTrap[1];
0629                     hapoth = parameterTrap[3];
0630                     uylfac = (htedge - hbedge) / (htedge + hbedge) / hapoth;
0631                     uxlden = 1 + yloc * uylfac;
0632                   }
0633
0634                   // Need to know position of trajectory in strip number for selecting the right APV later
0635                   if (TrajStrip == -1) {
0636                     int nstrips = Topo.nstrips();
0637                     float pitch = stripdet->surface().bounds().width() / nstrips;
0638                     TrajStrip = xloc / pitch + nstrips / 2.0;
0639                     // Need additionnal corrections for endcap
0640                     if (TKlayers >= 11) {
0641                       float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0642                                                           hapoth);  // radialy extrapolated x loc position at middle
0643                       TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
0644                     }
0645                     //LOGPRINT<<" Layer "<<TKlayers<<" TrajStrip: "<<nstrips<<" "<<pitch<<" "<<TrajStrip<<endl;
0646                   }
0647
0648                   for (edmNew::DetSet<SiStripCluster>::const_iterator iter = DSViter->begin(); iter != DSViter->end();
0649                        ++iter) {
0650                     //iter is a single SiStripCluster
0651                     StripClusterParameterEstimator::LocalValues parameters = stripcpe.localParameters(*iter, *stripdet);
0652                     float res = (parameters.first.x() - xloc);
0653                     float sigma = ::checkConsistency(parameters, xloc, xErr);
0654                     // The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
0655                     // you need a TransientTrackingRecHit instead of the cluster
0656                     //theEstimator=       new Chi2MeasurementEstimator(30);
0657                     //const Chi2MeasurementEstimator *theEstimator(100);
0658                     //theEstimator->estimate(TM->tsos(), TransientTrackingRecHit);
0659
0660                     if (TKlayers >= 11) {
0661                       res = parameters.first.x() - xloc / uxlden;  // radialy extrapolated x loc position at middle
0662                       sigma = abs(res) /
0663                               sqrt(parameters.second.xx() + xErr * xErr / uxlden / uxlden +
0664                                    yErr * yErr * xloc * xloc * uylfac * uylfac / uxlden / uxlden / uxlden / uxlden);
0665                     }
0666
0667                     siStripClusterInfo_.setCluster(*iter, ClusterId);
0668                     // signal to noise from SiStripClusterInfo not working in 225. I'll fix this after the interface
0669                     // redesign in 300 -ku
0670                     //float cluster_info[7] = {res, sigma, parameters.first.x(), sqrt(parameters.second.xx()), parameters.first.y(), sqrt(parameters.second.yy()), signal_to_noise};
0671                     std::vector<float> cluster_info;
0672                     cluster_info.push_back(res);
0673                     cluster_info.push_back(sigma);
0674                     cluster_info.push_back(parameters.first.x());
0675                     cluster_info.push_back(sqrt(parameters.second.xx()));
0676                     cluster_info.push_back(parameters.first.y());
0677                     cluster_info.push_back(sqrt(parameters.second.yy()));
0678                     cluster_info.push_back(siStripClusterInfo_.signalOverNoise());
0679                     VCluster_info.push_back(cluster_info);
0680                     nClusters++;
0681                     LogDebug("SiStripHitEfficiency:HitEff") << "Have ID match. residual = " << VCluster_info.back()[0]
0682                                                             << "  res sigma = " << VCluster_info.back()[1] << endl;
0683                     LogDebug("SiStripHitEfficiency:HitEff")
0684                         << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
0685                     LogDebug("SiStripHitEfficiency:HitEff")
0686                         << "hit position = " << parameters.first.x() << "  hit error = " << sqrt(parameters.second.xx())
0687                         << "  trajectory position = " << xloc << "  traj error = " << xErr << endl;
0688                   }
0689                 }
0690               }
0691               float FinalResSig = 1000.0;
0692               float FinalCluster[7] = {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
0693               if (nClusters > 0) {
0694                 LogDebug("SiStripHitEfficiency:HitEff") << "found clusters > 0" << endl;
0695                 if (nClusters > 1) {
0696                   //get the smallest one
0697                   vector<vector<float> >::iterator ires;
0698                   for (ires = VCluster_info.begin(); ires != VCluster_info.end(); ires++) {
0699                     if (abs((*ires)[1]) < abs(FinalResSig)) {
0700                       FinalResSig = (*ires)[1];
0701                       for (unsigned int i = 0; i < ires->size(); i++) {
0702                         LogDebug("SiStripHitEfficiency:HitEff")
0703                             << "filling final cluster. i = " << i << " before fill FinalCluster[i]=" << FinalCluster[i]
0704                             << " and (*ires)[i] =" << (*ires)[i] << endl;
0705                         FinalCluster[i] = (*ires)[i];
0706                         LogDebug("SiStripHitEfficiency:HitEff")
0707                             << "filling final cluster. i = " << i << " after fill FinalCluster[i]=" << FinalCluster[i]
0708                             << " and (*ires)[i] =" << (*ires)[i] << endl;
0709                       }
0710                     }
0711                     LogDebug("SiStripHitEfficiency:HitEff")
0712                         << "iresidual = " << (*ires)[0] << "  isigma = " << (*ires)[1]
0713                         << "  and FinalRes = " << FinalCluster[0] << endl;
0714                   }
0715                 } else {
0716                   FinalResSig = VCluster_info.at(0)[1];
0717                   for (unsigned int i = 0; i < VCluster_info.at(0).size(); i++) {
0718                     FinalCluster[i] = VCluster_info.at(0)[i];
0719                   }
0720                 }
0721                 VCluster_info.clear();
0722               }
0723
0724               LogDebug("SiStripHitEfficiency:HitEff")
0725                   << "Final residual in X = " << FinalCluster[0] << "+-" << (FinalCluster[0] / FinalResSig) << endl;
0726               LogDebug("SiStripHitEfficiency:HitEff") << "Checking location of trajectory: abs(yloc) = " << abs(yloc)
0727                                                       << "  abs(xloc) = " << abs(xloc) << endl;
0728               LogDebug("SiStripHitEfficiency:HitEff")
0729                   << "Checking location of cluster hit: yloc = " << FinalCluster[4] << "+-" << FinalCluster[5]
0730                   << "  xloc = " << FinalCluster[2] << "+-" << FinalCluster[3] << endl;
0731               LogDebug("SiStripHitEfficiency:HitEff") << "Final cluster signal to noise = " << FinalCluster[6] << endl;
0732
0733               float exclusionWidth = 0.4;
0734               float TOBexclusion = 0.0;
0735               float TECexRing5 = -0.89;
0736               float TECexRing6 = -0.56;
0737               float TECexRing7 = 0.60;
0738               //Added by Chris Edelmaier to do TEC bonding exclusion
0739               int subdetector = ((iidd >> 25) & 0x7);
0740               int ringnumber = ((iidd >> 5) & 0x7);
0741
0742               //New TOB and TEC bonding region exclusion zone
0743               if ((TKlayers >= 5 && TKlayers < 11) ||
0744                   ((subdetector == 6) && ((ringnumber >= 5) && (ringnumber <= 7)))) {
0745                 //There are only 2 cases that we need to exclude for
0746                 float highzone = 0.0;
0747                 float lowzone = 0.0;
0748                 float higherr = yloc + 5.0 * yErr;
0749                 float lowerr = yloc - 5.0 * yErr;
0750                 if (TKlayers >= 5 && TKlayers < 11) {
0751                   //TOB zone
0752                   highzone = TOBexclusion + exclusionWidth;
0753                   lowzone = TOBexclusion - exclusionWidth;
0754                 } else if (ringnumber == 5) {
0755                   //TEC ring 5
0756                   highzone = TECexRing5 + exclusionWidth;
0757                   lowzone = TECexRing5 - exclusionWidth;
0758                 } else if (ringnumber == 6) {
0759                   //TEC ring 6
0760                   highzone = TECexRing6 + exclusionWidth;
0761                   lowzone = TECexRing6 - exclusionWidth;
0762                 } else if (ringnumber == 7) {
0763                   //TEC ring 7
0764                   highzone = TECexRing7 + exclusionWidth;
0765                   lowzone = TECexRing7 - exclusionWidth;
0766                 }
0767                 //Now that we have our exclusion region, we just have to properly identify it
0768                 if ((highzone <= higherr) && (highzone >= lowerr))
0769                   withinAcceptance = false;
0770                 if ((lowzone >= lowerr) && (lowzone <= higherr))
0771                   withinAcceptance = false;
0772                 if ((higherr <= highzone) && (higherr >= lowzone))
0773                   withinAcceptance = false;
0774                 if ((lowerr >= lowzone) && (lowerr <= highzone))
0775                   withinAcceptance = false;
0776               }
0777
0778               // fill ntuple varibles
0779               //get global position from module id number iidd
0780               TrajGlbX = xglob;
0781               TrajGlbY = yglob;
0782               TrajGlbZ = zglob;
0783
0784               TrajLocErrX = xErr;
0785               TrajLocErrY = yErr;
0786
0787               Id = iidd;
0788               run = run_nr;
0789               event = ev_nr;
0790               bunchx = bunch_nr;
0791               //if ( SiStripQuality_->IsModuleBad(iidd) ) {
0792               if (SiStripQuality_->getBadApvs(iidd) != 0) {
0793                 SiStripQualBad = 1;
0794                 LogDebug("SiStripHitEfficiency:HitEff") << "strip is bad from SiStripQuality" << endl;
0795               } else {
0796                 SiStripQualBad = 0;
0797                 LogDebug("SiStripHitEfficiency:HitEff") << "strip is good from SiStripQuality" << endl;
0798               }
0799
0800               //check for FED-detected errors and include those in SiStripQualBad
0801               for (unsigned int ii = 0; ii < fedErrorIds->size(); ii++) {
0802                 if (iidd == (*fedErrorIds)[ii].rawId())
0803                   SiStripQualBad = 1;
0804               }
0805
0806               TrajLocX = xloc;
0807               TrajLocY = yloc;
0808               TrajLocAngleX = angleX;
0809               TrajLocAngleY = angleY;
0810               ResX = FinalCluster[0];
0811               ResXSig = FinalResSig;
0812               if (FinalResSig != FinalCluster[1])
0813                 LogDebug("SiStripHitEfficiency:HitEff")
0814                     << "Problem with best cluster selection because FinalResSig = " << FinalResSig
0815                     << " and FinalCluster[1] = " << FinalCluster[1] << endl;
0816               ClusterLocX = FinalCluster[2];
0817               ClusterLocY = FinalCluster[4];
0818               ClusterLocErrX = FinalCluster[3];
0819               ClusterLocErrY = FinalCluster[5];
0820               ClusterStoN = FinalCluster[6];
0821
0822               // CM of APV crossed by traj
0823               if (addCommonMode_)
0824                 if (commonModeDigis.isValid() && TrajStrip >= 0 && TrajStrip <= 768) {
0825                   edm::DetSetVector<SiStripRawDigi>::const_iterator digiframe = commonModeDigis->find(iidd);
0826                   if (digiframe != commonModeDigis->end())
0827                     if ((unsigned)TrajStrip / 128 < digiframe->data.size())
0828                       commonMode = digiframe->data.at(TrajStrip / 128).adc();
0829                 }
0830
0831               LogDebug("SiStripHitEfficiency:HitEff") << "before check good" << endl;
0832
0833               if (FinalResSig < 999.0) {  //could make requirement on track/hit consistency, but for
0834                 //now take anything with a hit on the module
0835                 LogDebug("SiStripHitEfficiency:HitEff")
0836                     << "hit being counted as good " << FinalCluster[0] << " FinalRecHit " << iidd << "   TKlayers  "
0837                     << TKlayers << " xloc " << xloc << " yloc  " << yloc << " module " << iidd
0838                     << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1) << "/"
0839                     << ((iidd & 0x3) == 2) << endl;
0840                 ModIsBad = 0;
0841                 traj->Fill();
0842               } else {
0843                 LogDebug("SiStripHitEfficiency:HitEff")
0844                     << "hit being counted as bad   ######### Invalid RPhi FinalResX " << FinalCluster[0]
0845                     << " FinalRecHit " << iidd << "   TKlayers  " << TKlayers << " xloc " << xloc << " yloc  " << yloc
0846                     << " module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0847                     << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
0848                 ModIsBad = 1;
0849                 traj->Fill();
0850
0851                 LogDebug("SiStripHitEfficiency:HitEff") << " RPhi Error " << sqrt(xErr * xErr + yErr * yErr)
0852                                                         << " ErrorX " << xErr << " yErr " << yErr << endl;
0853               }
0854               LogDebug("SiStripHitEfficiency:HitEff") << "after good location check" << endl;
0855             }
0856             LogDebug("SiStripHitEfficiency:HitEff") << "after list of clusters" << endl;
0857           }
0858           LogDebug("SiStripHitEfficiency:HitEff") << "After layers=TKLayers if" << endl;
0859         }
0860         LogDebug("SiStripHitEfficiency:HitEff") << "After looping over TrajAtValidHit list" << endl;
0861       }
0862       LogDebug("SiStripHitEfficiency:HitEff") << "end TMeasurement loop" << endl;
0863     }
0864     LogDebug("SiStripHitEfficiency:HitEff") << "end of trajectories loop" << endl;
0865   }
0866 }
0867
0868 void HitEff::endJob() {
0869   traj->GetDirectory()->cd();
0870   traj->Write();
0871
0872   LogDebug("SiStripHitEfficiency:HitEff") << " Events Analysed             " << events << endl;
0873   LogDebug("SiStripHitEfficiency:HitEff") << " Number Of Tracked events    " << EventTrackCKF << endl;
0874 }
0875
0876 //define this as a plug-in
0877 DEFINE_FWK_MODULE(HitEff);