File indexing completed on 2023-10-25 09:35:31
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include <memory>
0011 #include <string>
0012 #include <iostream>
0013
0014
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
0063 #include "TMath.h"
0064 #include "TH1F.h"
0065
0066
0067 #define LOGPRINT edm::LogPrint("SiStripHitEfficiency:HitEff")
0068
0069
0070
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
0179 const TrackerTopology* tTopo = &es.getData(topoToken_);
0180 siStripClusterInfo_.initEvent(es);
0181
0182
0183
0184 LogDebug("SiStripHitEfficiency:HitEff") << "beginning analyze from HitEff" << endl;
0185
0186 using namespace edm;
0187 using namespace reco;
0188
0189
0190 int run_nr = e.id().run();
0191 int ev_nr = e.id().event();
0192 int bunch_nr = e.bunchCrossing();
0193
0194
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
0215 edm::Handle<edm::DetSetVector<SiStripRawDigi> > commonModeDigis;
0216 if (addCommonMode_)
0217 e.getByToken(commonModeToken_, commonModeDigis);
0218
0219
0220 edm::Handle<reco::TrackCollection> trackCollectionCKF;
0221
0222 e.getByToken(combinatorialTracks_token_, trackCollectionCKF);
0223
0224 edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
0225
0226 e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0227
0228 edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0229 e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0230
0231
0232
0233 edm::Handle<edmNew::DetSetVector<SiStripCluster> > theClusters;
0234
0235 e.getByToken(clusters_token_, theClusters);
0236
0237
0238 edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_);
0239 const TrackerGeometry* tkgeom = &(*tracker);
0240
0241
0242
0243 edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_);
0244 const StripClusterParameterEstimator& stripcpe(*parameterestimator);
0245
0246
0247 edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_);
0248
0249 const MagneticField* magField_ = &es.getData(magFieldToken_);
0250
0251
0252
0253
0254
0255
0256
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
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
0278 const edmNew::DetSetVector<SiStripCluster>& input = *theClusters;
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307 const reco::TrackCollection* tracksCKF = trackCollectionCKF.product();
0308 LogDebug("SiStripHitEfficiency:HitEff") << "number ckf tracks found = " << tracksCKF->size() << endl;
0309
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
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
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
0356
0357
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
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
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
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 diffPreviousLayer = (layer - previous_layer);
0400 if (doMissingHitsRecovery_) {
0401
0402 if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
0403 missHitPerLayer[missedLayer] += 1;
0404 hasMissingHits = true;
0405 }
0406
0407 else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
0408 missHitPerLayer[missedLayer] += 1;
0409 hasMissingHits = true;
0410 }
0411
0412 else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
0413 missHitPerLayer[missedLayer] += 1;
0414 hasMissingHits = true;
0415 }
0416
0417
0418 if ((layer > k_LayersStart && layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
0419 missHitPerLayer[11] += 1;
0420 hasMissingHits = true;
0421 }
0422
0423
0424 if ((layer > k_LayersAtTIBEnd && layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
0425 missHitPerLayer[14] += 1;
0426 hasMissingHits = true;
0427 }
0428 }
0429 if (theHit->getType() == TrackingRecHit::Type::missing)
0430 hasMissingHits = true;
0431
0432 if (hasMissingHits)
0433 missedLayers.push_back(layer);
0434 previous_layer = layer;
0435 }
0436
0437
0438
0439 unsigned int prev_TKlayers = 0;
0440 for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
0441 auto theInHit = (*itm).recHit();
0442
0443 LogDebug("SiStripHitEfficiency:HitEff") << "theInHit is valid = " << theInHit->isValid() << endl;
0444
0445 unsigned int iidd = theInHit->geographicalId().rawId();
0446
0447 unsigned int TKlayers = ::checkLayer(iidd, tTopo);
0448 LogDebug("SiStripHitEfficiency:HitEff") << "TKlayer from trajectory: " << TKlayers << " from module = " << iidd
0449 << " matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0450 << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
0451
0452
0453
0454 bool isFirstMeas = (itm == (TMeas.end() - 1));
0455 bool isLastMeas = (itm == (TMeas.begin()));
0456
0457 if (!useFirstMeas_ && isFirstMeas)
0458 continue;
0459 if (!useLastMeas_ && isLastMeas)
0460 continue;
0461
0462
0463 if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
0464 !useAllHitsFromTracksWithMissingHits_)
0465 continue;
0466
0467
0468 if (TKlayers == 10 || TKlayers == 22) {
0469 LogDebug("SiStripHitEfficiency:HitEff") << "skipping original TM for TOB 6 or TEC 9" << endl;
0470 continue;
0471 }
0472
0473
0474 std::vector<TrajectoryAtInvalidHit> TMs;
0475
0476
0477 AnalyticalPropagator propagator(magField_, anyDirection);
0478
0479
0480
0481
0482 if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0483
0484
0485 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0486 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0487 } else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
0488
0489
0490 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
0491 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
0492 LogDebug("SiStripHitEfficiency:HitEff") << " found a hit with a missing partner";
0493 } else {
0494
0495 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator));
0496 }
0497 bool missingHitAdded = false;
0498
0499 vector<TrajectoryMeasurement> tmpTmeas;
0500 unsigned int misLayer = TKlayers + 1;
0501
0502 if (doMissingHitsRecovery_) {
0503 if (int(TKlayers) - int(prev_TKlayers) == -2) {
0504 const DetLayer* detlayer = itm->layer();
0505 const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0506 const TrajectoryStateOnSurface tsos = itm->updatedState();
0507 std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0508
0509 if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) {
0510 std::vector<ForwardDetLayer const*> negTECLayers =
0511 measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0512 std::vector<ForwardDetLayer const*> posTECLayers =
0513 measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0514 const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0515 const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0516 if (tTopo->tecSide(iidd) == 1) {
0517 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
0518 } else if (tTopo->tecSide(iidd) == 2) {
0519 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
0520 }
0521 }
0522
0523 else if (misLayer == (k_LayersAtTIDEnd - 1) ||
0524 misLayer == k_LayersAtTIDEnd) {
0525
0526 std::vector<ForwardDetLayer const*> negTIDLayers =
0527 measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
0528 std::vector<ForwardDetLayer const*> posTIDLayers =
0529 measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
0530 const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0531 const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0532
0533 if (tTopo->tidSide(iidd) == 1) {
0534 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
0535 } else if (tTopo->tidSide(iidd) == 2) {
0536 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
0537 }
0538 }
0539
0540 if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) {
0541
0542 std::vector<BarrelDetLayer const*> barrelTIBLayers =
0543 measurementTrackerHandle->geometricSearchTracker()->tibLayers();
0544 std::vector<BarrelDetLayer const*> barrelTOBLayers =
0545 measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0546
0547 if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
0548 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
0549 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, *thePropagator, *estimator);
0550 } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
0551 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
0552 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, *thePropagator, *estimator);
0553 }
0554 }
0555 }
0556 if ((int(TKlayers) > k_LayersStart && int(TKlayers) <= k_LayersAtTIBEnd) && int(prev_TKlayers) == 12) {
0557 const DetLayer* detlayer = itm->layer();
0558 const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0559 const TrajectoryStateOnSurface tsos = itm->updatedState();
0560 std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0561 std::vector<ForwardDetLayer const*> negTIDLayers =
0562 measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
0563 std::vector<ForwardDetLayer const*> posTIDLayers =
0564 measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
0565
0566 const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
0567 const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
0568 if (tTopo->tidSide(iidd) == 1) {
0569 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
0570 } else if (tTopo->tidSide(iidd) == 2) {
0571 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
0572 }
0573 }
0574
0575 if ((int(TKlayers) > k_LayersAtTIBEnd && int(TKlayers) <= k_LayersAtTOBEnd) && int(prev_TKlayers) == 15) {
0576 const DetLayer* detlayer = itm->layer();
0577 const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0578 const TrajectoryStateOnSurface tsos = itm->updatedState();
0579 std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
0580
0581 std::vector<ForwardDetLayer const*> negTECLayers =
0582 measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0583 std::vector<ForwardDetLayer const*> posTECLayers =
0584 measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0585
0586 const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
0587 const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
0588 if (tTopo->tecSide(iidd) == 1) {
0589 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
0590 } else if (tTopo->tecSide(iidd) == 2) {
0591 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
0592 }
0593 }
0594
0595 if (!tmpTmeas.empty()) {
0596 TrajectoryMeasurement TM_tmp(tmpTmeas.back());
0597 unsigned int iidd_tmp = TM_tmp.recHit()->geographicalId().rawId();
0598 if (iidd_tmp != 0) {
0599 LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0600 if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0601 TMs.clear();
0602 if (::isDoubleSided(iidd_tmp, tTopo)) {
0603 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 1));
0604 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 2));
0605 } else
0606 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator));
0607 missingHitAdded = true;
0608 hitRecoveryCounters[misLayer] += 1;
0609 }
0610 }
0611 }
0612
0613 prev_TKlayers = TKlayers;
0614 if (!useFirstMeas_ && isFirstMeas && !missingHitAdded)
0615 continue;
0616 if (!useLastMeas_ && isLastMeas)
0617 continue;
0618 bool hitsWithBias = false;
0619 for (auto ilayer : missedLayers) {
0620 if (ilayer < TKlayers)
0621 hitsWithBias = true;
0622 }
0623 if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing && !missingHitAdded &&
0624 hitsWithBias && !useAllHitsFromTracksWithMissingHits_) {
0625 continue;
0626 }
0627
0628
0629
0630
0631
0632
0633 bool isValid = theInHit->isValid();
0634 bool isLast = (itm == (TMeas.end() - 1));
0635 bool isLastTOB5 = true;
0636 if (!isLast) {
0637 if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9)
0638 isLastTOB5 = false;
0639 else
0640 isLastTOB5 = true;
0641 --itm;
0642 }
0643
0644 if (TKlayers == 9 && isValid && isLastTOB5) {
0645
0646
0647 std::vector<BarrelDetLayer const*> barrelTOBLayers =
0648 measurementTrackerHandle->geometricSearchTracker()->tobLayers();
0649 const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size() - 1];
0650 const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0651 const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
0652 auto tmp = layerMeasurements.measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
0653
0654 if (!tmp.empty()) {
0655 LogDebug("SiStripHitEfficiency:HitEff") << "size of TM from propagation = " << tmp.size() << endl;
0656
0657
0658
0659
0660 TrajectoryMeasurement tob6TM(tmp.back());
0661 const auto& tob6Hit = tob6TM.recHit();
0662
0663 if (tob6Hit->geographicalId().rawId() != 0) {
0664 LogDebug("SiStripHitEfficiency:HitEff") << "tob6 hit actually being added to TM vector" << endl;
0665 TMs.push_back(TrajectoryAtInvalidHit(tob6TM, tTopo, tkgeom, propagator));
0666 }
0667 }
0668 }
0669
0670 bool isLastTEC8 = true;
0671 if (!isLast) {
0672 if (::checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 21)
0673 isLastTEC8 = false;
0674 else
0675 isLastTEC8 = true;
0676 --itm;
0677 }
0678
0679 if (TKlayers == 21 && isValid && isLastTEC8) {
0680 std::vector<const ForwardDetLayer*> posTecLayers =
0681 measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
0682 const DetLayer* tec9pos = posTecLayers[posTecLayers.size() - 1];
0683 std::vector<const ForwardDetLayer*> negTecLayers =
0684 measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
0685 const DetLayer* tec9neg = negTecLayers[negTecLayers.size() - 1];
0686 const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
0687 const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
0688
0689
0690 if (!(iidd == StripSubdetector::TEC))
0691 LogDebug("SiStripHitEfficiency:HitEff") << "there is a problem with TEC 9 extrapolation" << endl;
0692
0693
0694 vector<TrajectoryMeasurement> tmp;
0695 if (tTopo->tecSide(iidd) == 1) {
0696 tmp = layerMeasurements.measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
0697
0698 }
0699 if (tTopo->tecSide(iidd) == 2) {
0700 tmp = layerMeasurements.measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
0701
0702 }
0703
0704 if (!tmp.empty()) {
0705
0706
0707
0708 TrajectoryMeasurement tec9TM(tmp.back());
0709 const auto& tec9Hit = tec9TM.recHit();
0710
0711 unsigned int tec9id = tec9Hit->geographicalId().rawId();
0712 LogDebug("SiStripHitEfficiency:HitEff")
0713 << "tec9id = " << tec9id << " is Double sided = " << ::isDoubleSided(tec9id, tTopo)
0714 << " and 0x3 = " << (tec9id & 0x3) << endl;
0715
0716 if (tec9Hit->geographicalId().rawId() != 0) {
0717 LogDebug("SiStripHitEfficiency:HitEff") << "tec9 hit actually being added to TM vector" << endl;
0718
0719
0720 if (::isDoubleSided(tec9id, tTopo)) {
0721 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 1));
0722 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator, 2));
0723 } else
0724 TMs.push_back(TrajectoryAtInvalidHit(tec9TM, tTopo, tkgeom, propagator));
0725 }
0726 }
0727 }
0728 hitTotalCounters[TKlayers] += 1;
0729
0730
0731
0732
0733
0734 for (std::vector<TrajectoryAtInvalidHit>::const_iterator TM = TMs.begin(); TM != TMs.end(); ++TM) {
0735
0736 iidd = TM->monodet_id();
0737 LogDebug("SiStripHitEfficiency:HitEff") << "setting iidd = " << iidd << " before checking efficiency and ";
0738
0739 xloc = TM->localX();
0740 yloc = TM->localY();
0741
0742 angleX = atan(TM->localDxDz());
0743 angleY = atan(TM->localDyDz());
0744
0745 TrajLocErrX = 0.0;
0746 TrajLocErrY = 0.0;
0747
0748 xglob = TM->globalX();
0749 yglob = TM->globalY();
0750 zglob = TM->globalZ();
0751 xErr = TM->localErrorX();
0752 yErr = TM->localErrorY();
0753
0754 TrajGlbX = 0.0;
0755 TrajGlbY = 0.0;
0756 TrajGlbZ = 0.0;
0757 withinAcceptance = TM->withinAcceptance();
0758
0759 trajHitValid = TM->validHit();
0760 int TrajStrip = -1;
0761
0762
0763 TKlayers = ::checkLayer(iidd, tTopo);
0764
0765 if ((layers == TKlayers) || (layers == 0)) {
0766 whatlayer = TKlayers;
0767 LogDebug("SiStripHitEfficiency:HitEff") << "Looking at layer under study" << endl;
0768 ModIsBad = 2;
0769 Id = 0;
0770 SiStripQualBad = 0;
0771 run = 0;
0772 event = 0;
0773 TrajLocX = 0.0;
0774 TrajLocY = 0.0;
0775 TrajLocAngleX = -999.0;
0776 TrajLocAngleY = -999.0;
0777 ResX = 0.0;
0778 ResXSig = 0.0;
0779 ClusterLocX = 0.0;
0780 ClusterLocY = 0.0;
0781 ClusterLocErrX = 0.0;
0782 ClusterLocErrY = 0.0;
0783 ClusterStoN = 0.0;
0784 bunchx = 0;
0785 commonMode = -100;
0786
0787
0788
0789 if (!input.empty()) {
0790 LogDebug("SiStripHitEfficiency:HitEff") << "Checking clusters with size = " << input.size() << endl;
0791 int nClusters = 0;
0792 std::vector<std::vector<float> >
0793 VCluster_info;
0794 for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end();
0795 DSViter++) {
0796
0797
0798 unsigned int ClusterId = DSViter->id();
0799 if (ClusterId == iidd) {
0800 LogDebug("SiStripHitEfficiency:HitEff")
0801 << "found (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
0802 DetId ClusterDetId(ClusterId);
0803 const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
0804 const StripTopology& Topo = stripdet->specificTopology();
0805
0806 float hbedge = 0.0;
0807 float htedge = 0.0;
0808 float hapoth = 0.0;
0809 float uylfac = 0.0;
0810 float uxlden = 0.0;
0811 if (TKlayers >= 11) {
0812 const BoundPlane& plane = stripdet->surface();
0813 const TrapezoidalPlaneBounds* trapezoidalBounds(
0814 dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0815 std::array<const float, 4> const& parameterTrap =
0816 (*trapezoidalBounds).parameters();
0817 hbedge = parameterTrap[0];
0818 htedge = parameterTrap[1];
0819 hapoth = parameterTrap[3];
0820 uylfac = (htedge - hbedge) / (htedge + hbedge) / hapoth;
0821 uxlden = 1 + yloc * uylfac;
0822 }
0823
0824
0825 if (TrajStrip == -1) {
0826 int nstrips = Topo.nstrips();
0827 float pitch = stripdet->surface().bounds().width() / nstrips;
0828 TrajStrip = xloc / pitch + nstrips / 2.0;
0829
0830 if (TKlayers >= 11) {
0831 float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0832 hapoth);
0833 TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
0834 }
0835
0836 }
0837
0838 for (edmNew::DetSet<SiStripCluster>::const_iterator iter = DSViter->begin(); iter != DSViter->end();
0839 ++iter) {
0840
0841 StripClusterParameterEstimator::LocalValues parameters = stripcpe.localParameters(*iter, *stripdet);
0842 float res = (parameters.first.x() - xloc);
0843 float sigma = ::checkConsistency(parameters, xloc, xErr);
0844
0845
0846
0847
0848
0849
0850 if (TKlayers >= 11) {
0851 res = parameters.first.x() - xloc / uxlden;
0852 sigma = abs(res) /
0853 sqrt(parameters.second.xx() + xErr * xErr / uxlden / uxlden +
0854 yErr * yErr * xloc * xloc * uylfac * uylfac / uxlden / uxlden / uxlden / uxlden);
0855 }
0856
0857 siStripClusterInfo_.setCluster(*iter, ClusterId);
0858
0859
0860
0861 std::vector<float> cluster_info;
0862 cluster_info.push_back(res);
0863 cluster_info.push_back(sigma);
0864 cluster_info.push_back(parameters.first.x());
0865 cluster_info.push_back(sqrt(parameters.second.xx()));
0866 cluster_info.push_back(parameters.first.y());
0867 cluster_info.push_back(sqrt(parameters.second.yy()));
0868 cluster_info.push_back(siStripClusterInfo_.signalOverNoise());
0869 VCluster_info.push_back(cluster_info);
0870 nClusters++;
0871 LogDebug("SiStripHitEfficiency:HitEff") << "Have ID match. residual = " << VCluster_info.back()[0]
0872 << " res sigma = " << VCluster_info.back()[1] << endl;
0873 LogDebug("SiStripHitEfficiency:HitEff")
0874 << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
0875 LogDebug("SiStripHitEfficiency:HitEff")
0876 << "hit position = " << parameters.first.x() << " hit error = " << sqrt(parameters.second.xx())
0877 << " trajectory position = " << xloc << " traj error = " << xErr << endl;
0878 }
0879 }
0880 }
0881 float FinalResSig = 1000.0;
0882 float FinalCluster[7] = {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
0883 if (nClusters > 0) {
0884 LogDebug("SiStripHitEfficiency:HitEff") << "found clusters > 0" << endl;
0885 if (nClusters > 1) {
0886
0887 vector<vector<float> >::iterator ires;
0888 for (ires = VCluster_info.begin(); ires != VCluster_info.end(); ires++) {
0889 if (abs((*ires)[1]) < abs(FinalResSig)) {
0890 FinalResSig = (*ires)[1];
0891 for (unsigned int i = 0; i < ires->size(); i++) {
0892 LogDebug("SiStripHitEfficiency:HitEff")
0893 << "filling final cluster. i = " << i << " before fill FinalCluster[i]=" << FinalCluster[i]
0894 << " and (*ires)[i] =" << (*ires)[i] << endl;
0895 FinalCluster[i] = (*ires)[i];
0896 LogDebug("SiStripHitEfficiency:HitEff")
0897 << "filling final cluster. i = " << i << " after fill FinalCluster[i]=" << FinalCluster[i]
0898 << " and (*ires)[i] =" << (*ires)[i] << endl;
0899 }
0900 }
0901 LogDebug("SiStripHitEfficiency:HitEff")
0902 << "iresidual = " << (*ires)[0] << " isigma = " << (*ires)[1]
0903 << " and FinalRes = " << FinalCluster[0] << endl;
0904 }
0905 } else {
0906 FinalResSig = VCluster_info.at(0)[1];
0907 for (unsigned int i = 0; i < VCluster_info.at(0).size(); i++) {
0908 FinalCluster[i] = VCluster_info.at(0)[i];
0909 }
0910 }
0911 VCluster_info.clear();
0912 }
0913
0914 LogDebug("SiStripHitEfficiency:HitEff")
0915 << "Final residual in X = " << FinalCluster[0] << "+-" << (FinalCluster[0] / FinalResSig) << endl;
0916 LogDebug("SiStripHitEfficiency:HitEff") << "Checking location of trajectory: abs(yloc) = " << abs(yloc)
0917 << " abs(xloc) = " << abs(xloc) << endl;
0918 LogDebug("SiStripHitEfficiency:HitEff")
0919 << "Checking location of cluster hit: yloc = " << FinalCluster[4] << "+-" << FinalCluster[5]
0920 << " xloc = " << FinalCluster[2] << "+-" << FinalCluster[3] << endl;
0921 LogDebug("SiStripHitEfficiency:HitEff") << "Final cluster signal to noise = " << FinalCluster[6] << endl;
0922
0923 float exclusionWidth = 0.4;
0924 float TOBexclusion = 0.0;
0925 float TECexRing5 = -0.89;
0926 float TECexRing6 = -0.56;
0927 float TECexRing7 = 0.60;
0928
0929 int subdetector = ((iidd >> 25) & 0x7);
0930 int ringnumber = ((iidd >> 5) & 0x7);
0931
0932
0933 if ((TKlayers >= 5 && TKlayers < 11) ||
0934 ((subdetector == 6) && ((ringnumber >= 5) && (ringnumber <= 7)))) {
0935
0936 float highzone = 0.0;
0937 float lowzone = 0.0;
0938 float higherr = yloc + 5.0 * yErr;
0939 float lowerr = yloc - 5.0 * yErr;
0940 if (TKlayers >= 5 && TKlayers < 11) {
0941
0942 highzone = TOBexclusion + exclusionWidth;
0943 lowzone = TOBexclusion - exclusionWidth;
0944 } else if (ringnumber == 5) {
0945
0946 highzone = TECexRing5 + exclusionWidth;
0947 lowzone = TECexRing5 - exclusionWidth;
0948 } else if (ringnumber == 6) {
0949
0950 highzone = TECexRing6 + exclusionWidth;
0951 lowzone = TECexRing6 - exclusionWidth;
0952 } else if (ringnumber == 7) {
0953
0954 highzone = TECexRing7 + exclusionWidth;
0955 lowzone = TECexRing7 - exclusionWidth;
0956 }
0957
0958 if ((highzone <= higherr) && (highzone >= lowerr))
0959 withinAcceptance = false;
0960 if ((lowzone >= lowerr) && (lowzone <= higherr))
0961 withinAcceptance = false;
0962 if ((higherr <= highzone) && (higherr >= lowzone))
0963 withinAcceptance = false;
0964 if ((lowerr >= lowzone) && (lowerr <= highzone))
0965 withinAcceptance = false;
0966 }
0967
0968
0969
0970 TrajGlbX = xglob;
0971 TrajGlbY = yglob;
0972 TrajGlbZ = zglob;
0973
0974 TrajLocErrX = xErr;
0975 TrajLocErrY = yErr;
0976
0977 Id = iidd;
0978 run = run_nr;
0979 event = ev_nr;
0980 bunchx = bunch_nr;
0981
0982 if (SiStripQuality_->getBadApvs(iidd) != 0) {
0983 SiStripQualBad = 1;
0984 LogDebug("SiStripHitEfficiency:HitEff") << "strip is bad from SiStripQuality" << endl;
0985 } else {
0986 SiStripQualBad = 0;
0987 LogDebug("SiStripHitEfficiency:HitEff") << "strip is good from SiStripQuality" << endl;
0988 }
0989
0990
0991 for (unsigned int ii = 0; ii < fedErrorIds.size(); ii++) {
0992 if (iidd == fedErrorIds[ii].rawId())
0993 SiStripQualBad = 1;
0994 }
0995
0996 TrajLocX = xloc;
0997 TrajLocY = yloc;
0998 TrajLocAngleX = angleX;
0999 TrajLocAngleY = angleY;
1000 ResX = FinalCluster[0];
1001 ResXSig = FinalResSig;
1002 if (FinalResSig != FinalCluster[1])
1003 LogDebug("SiStripHitEfficiency:HitEff")
1004 << "Problem with best cluster selection because FinalResSig = " << FinalResSig
1005 << " and FinalCluster[1] = " << FinalCluster[1] << endl;
1006 ClusterLocX = FinalCluster[2];
1007 ClusterLocY = FinalCluster[4];
1008 ClusterLocErrX = FinalCluster[3];
1009 ClusterLocErrY = FinalCluster[5];
1010 ClusterStoN = FinalCluster[6];
1011
1012
1013 if (addCommonMode_)
1014 if (commonModeDigis.isValid() && TrajStrip >= 0 && TrajStrip <= 768) {
1015 edm::DetSetVector<SiStripRawDigi>::const_iterator digiframe = commonModeDigis->find(iidd);
1016 if (digiframe != commonModeDigis->end())
1017 if ((unsigned)TrajStrip / 128 < digiframe->data.size())
1018 commonMode = digiframe->data.at(TrajStrip / 128).adc();
1019 }
1020
1021 LogDebug("SiStripHitEfficiency:HitEff") << "before check good" << endl;
1022
1023 if (FinalResSig < 999.0) {
1024
1025 LogDebug("SiStripHitEfficiency:HitEff")
1026 << "hit being counted as good " << FinalCluster[0] << " FinalRecHit " << iidd << " TKlayers "
1027 << TKlayers << " xloc " << xloc << " yloc " << yloc << " module " << iidd
1028 << " matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1) << "/"
1029 << ((iidd & 0x3) == 2) << endl;
1030 ModIsBad = 0;
1031 traj->Fill();
1032 } else {
1033 LogDebug("SiStripHitEfficiency:HitEff")
1034 << "hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0]
1035 << " FinalRecHit " << iidd << " TKlayers " << TKlayers << " xloc " << xloc << " yloc " << yloc
1036 << " module " << iidd << " matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
1037 << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2) << endl;
1038 ModIsBad = 1;
1039 traj->Fill();
1040
1041 LogDebug("SiStripHitEfficiency:HitEff") << " RPhi Error " << sqrt(xErr * xErr + yErr * yErr)
1042 << " ErrorX " << xErr << " yErr " << yErr << endl;
1043 }
1044 LogDebug("SiStripHitEfficiency:HitEff") << "after good location check" << endl;
1045 }
1046 LogDebug("SiStripHitEfficiency:HitEff") << "after list of clusters" << endl;
1047 }
1048 LogDebug("SiStripHitEfficiency:HitEff") << "After layers=TKLayers if" << endl;
1049 }
1050 LogDebug("SiStripHitEfficiency:HitEff") << "After looping over TrajAtValidHit list" << endl;
1051 }
1052 LogDebug("SiStripHitEfficiency:HitEff") << "end TMeasurement loop" << endl;
1053 }
1054 LogDebug("SiStripHitEfficiency:HitEff") << "end of trajectories loop" << endl;
1055 }
1056 }
1057
1058 void HitEff::endJob() {
1059 traj->GetDirectory()->cd();
1060 traj->Write();
1061
1062 LogDebug("SiStripHitEfficiency:HitEff") << " Events Analysed " << events << endl;
1063 LogDebug("SiStripHitEfficiency:HitEff") << " Number Of Tracked events " << EventTrackCKF << endl;
1064
1065 if (doMissingHitsRecovery_) {
1066 float totTIB = 0.0;
1067 float totTOB = 0.0;
1068 float totTID = 0.0;
1069 float totTEC = 0.0;
1070
1071 float totTIBrepro = 0.0;
1072 float totTOBrepro = 0.0;
1073 float totTIDrepro = 0.0;
1074 float totTECrepro = 0.0;
1075
1076 edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TIB :";
1077 for (int i = 0; i <= k_LayersAtTIBEnd; i++) {
1078 edm::LogInfo("SiStripHitEfficiency:HitEff")
1079 << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1080 << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1081 totTIB += missHitPerLayer[i];
1082 edm::LogInfo("SiStripHitEfficiency:HitEff")
1083 << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1084 << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1085 << " % of missing hit";
1086 totTIBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1087 }
1088 edm::LogInfo("SiStripHitEfficiency:HitEff")
1089 << "TOTAL % of missing hits within TIB :" << (totTIB * 1.0 / totalNbHits) * 100 << "%";
1090 edm::LogInfo("SiStripHitEfficiency:HitEff")
1091 << "AFTER repropagation :" << (totTIBrepro * 1.0 / totalNbHits) * 100 << "%";
1092
1093 edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TOB :";
1094 for (int i = k_LayersAtTIBEnd + 1; i <= k_LayersAtTOBEnd; i++) {
1095 edm::LogInfo("SiStripHitEfficiency:HitEff")
1096 << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1097 << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1098 totTOB += missHitPerLayer[i];
1099 edm::LogInfo("SiStripHitEfficiency:HitEff")
1100 << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1101 << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1102 << " % of missing hit";
1103 totTOBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1104 }
1105 edm::LogInfo("SiStripHitEfficiency:HitEff")
1106 << "TOTAL % of missing hits within TOB :" << (totTOB * 1.0 / totalNbHits) * 100 << "%";
1107 edm::LogInfo("SiStripHitEfficiency:HitEff")
1108 << "AFTER repropagation :" << (totTOBrepro * 1.0 / totalNbHits) * 100 << "%";
1109
1110 edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TID :";
1111 for (int i = k_LayersAtTOBEnd + 1; i <= k_LayersAtTIDEnd; i++) {
1112 edm::LogInfo("SiStripHitEfficiency:HitEff")
1113 << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1114 << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1115 totTID += missHitPerLayer[i];
1116 edm::LogInfo("SiStripHitEfficiency:HitEff")
1117 << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1118 << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1119 << " % of missing hit";
1120 totTIDrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1121 }
1122 edm::LogInfo("SiStripHitEfficiency:HitEff")
1123 << "TOTAL % of missing hits within TID :" << (totTID * 1.0 / totalNbHits) * 100 << "%";
1124 edm::LogInfo("SiStripHitEfficiency:HitEff")
1125 << "AFTER repropagation :" << (totTIDrepro * 1.0 / totalNbHits) * 100 << "%";
1126
1127 edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TEC :";
1128 for (int i = k_LayersAtTIDEnd + 1; i < k_END_OF_LAYERS; i++) {
1129 edm::LogInfo("SiStripHitEfficiency:HitEff")
1130 << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
1131 << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
1132 totTEC += missHitPerLayer[i];
1133 edm::LogInfo("SiStripHitEfficiency:HitEff")
1134 << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
1135 << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
1136 << " % of missing hit";
1137 totTECrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
1138 }
1139 edm::LogInfo("SiStripHitEfficiency:HitEff")
1140 << "TOTAL % of missing hits within TEC :" << (totTEC * 1.0 / totalNbHits) * 100 << "%";
1141 edm::LogInfo("SiStripHitEfficiency:HitEff")
1142 << "AFTER repropagation :" << (totTECrepro * 1.0 / totalNbHits) * 100 << "%";
1143
1144 edm::LogInfo("SiStripHitEfficiency:HitEff") << " Hit recovery summary:";
1145
1146 for (int ilayer = 0; ilayer < k_END_OF_LAYERS; ilayer++) {
1147 edm::LogInfo("SiStripHitEfficiency:HitEff")
1148 << " layer " << ilayer << ": " << hitRecoveryCounters[ilayer] << " / " << hitTotalCounters[ilayer];
1149 }
1150 }
1151 }
1152
1153
1154 DEFINE_FWK_MODULE(HitEff);