File indexing completed on 2022-06-16 22:18: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 "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
0062 #include "TMath.h"
0063 #include "TH1F.h"
0064
0065
0066 #define LOGPRINT edm::LogPrint("SiStripHitEfficiency:HitEff")
0067
0068
0069
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
0169 const TrackerTopology* tTopo = &es.getData(topoToken_);
0170 siStripClusterInfo_.initEvent(es);
0171
0172
0173
0174 LogDebug("SiStripHitEfficiency:HitEff") << "beginning analyze from HitEff" << endl;
0175
0176 using namespace edm;
0177 using namespace reco;
0178
0179
0180 int run_nr = e.id().run();
0181 int ev_nr = e.id().event();
0182 int bunch_nr = e.bunchCrossing();
0183
0184
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
0205 edm::Handle<edm::DetSetVector<SiStripRawDigi> > commonModeDigis;
0206 if (addCommonMode_)
0207 e.getByToken(commonModeToken_, commonModeDigis);
0208
0209
0210 edm::Handle<reco::TrackCollection> trackCollectionCKF;
0211
0212 e.getByToken(combinatorialTracks_token_, trackCollectionCKF);
0213
0214 edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
0215
0216 e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0217
0218 edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0219 e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0220
0221
0222
0223 edm::Handle<edmNew::DetSetVector<SiStripCluster> > theClusters;
0224
0225 e.getByToken(clusters_token_, theClusters);
0226
0227
0228 edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_);
0229 const TrackerGeometry* tkgeom = &(*tracker);
0230
0231
0232
0233 edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_);
0234 const StripClusterParameterEstimator& stripcpe(*parameterestimator);
0235
0236
0237 edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_);
0238
0239 const MagneticField* magField_ = &es.getData(magFieldToken_);
0240
0241
0242 edm::Handle<DetIdCollection> fedErrorIds;
0243
0244 e.getByToken(digis_token_, fedErrorIds);
0245
0246 ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTkToken_);
0247
0248 edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
0249
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
0258 const edmNew::DetSetVector<SiStripCluster>& input = *theClusters;
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287 const reco::TrackCollection* tracksCKF = trackCollectionCKF.product();
0288 LogDebug("SiStripHitEfficiency:HitEff") << "number ckf tracks found = " << tracksCKF->size() << endl;
0289
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
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
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
0336
0337
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
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
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
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
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
0393
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
0403 if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
0404 !useAllHitsFromTracksWithMissingHits_)
0405 continue;
0406
0407
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
0414 std::vector<TrajectoryAtInvalidHit> TMs;
0415
0416
0417 AnalyticalPropagator propagator(magField_, anyDirection);
0418
0419
0420
0421
0422 if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0423
0424
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
0429
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
0435 TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator));
0436 }
0437
0438
0439
0440
0441
0442
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
0457
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
0469
0470
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
0501 if (!(iidd == StripSubdetector::TEC))
0502 LogDebug("SiStripHitEfficiency:HitEff") << "there is a problem with TEC 9 extrapolation" << endl;
0503
0504
0505 vector<TrajectoryMeasurement> tmp;
0506 if (tTopo->tecSide(iidd) == 1) {
0507 tmp = layerMeasurements.measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
0508
0509 }
0510 if (tTopo->tecSide(iidd) == 2) {
0511 tmp = layerMeasurements.measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
0512
0513 }
0514
0515 if (!tmp.empty()) {
0516
0517
0518
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
0530
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 }
0538 }
0539
0540
0541
0542
0543
0544 for (std::vector<TrajectoryAtInvalidHit>::const_iterator TM = TMs.begin(); TM != TMs.end(); ++TM) {
0545
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
0573 TKlayers = ::checkLayer(iidd, tTopo);
0574
0575 if ((layers == TKlayers) || (layers == 0)) {
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
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;
0604 for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end();
0605 DSViter++) {
0606
0607
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();
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
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
0640 if (TKlayers >= 11) {
0641 float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0642 hapoth);
0643 TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
0644 }
0645
0646 }
0647
0648 for (edmNew::DetSet<SiStripCluster>::const_iterator iter = DSViter->begin(); iter != DSViter->end();
0649 ++iter) {
0650
0651 StripClusterParameterEstimator::LocalValues parameters = stripcpe.localParameters(*iter, *stripdet);
0652 float res = (parameters.first.x() - xloc);
0653 float sigma = ::checkConsistency(parameters, xloc, xErr);
0654
0655
0656
0657
0658
0659
0660 if (TKlayers >= 11) {
0661 res = parameters.first.x() - xloc / uxlden;
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
0669
0670
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
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
0739 int subdetector = ((iidd >> 25) & 0x7);
0740 int ringnumber = ((iidd >> 5) & 0x7);
0741
0742
0743 if ((TKlayers >= 5 && TKlayers < 11) ||
0744 ((subdetector == 6) && ((ringnumber >= 5) && (ringnumber <= 7)))) {
0745
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
0752 highzone = TOBexclusion + exclusionWidth;
0753 lowzone = TOBexclusion - exclusionWidth;
0754 } else if (ringnumber == 5) {
0755
0756 highzone = TECexRing5 + exclusionWidth;
0757 lowzone = TECexRing5 - exclusionWidth;
0758 } else if (ringnumber == 6) {
0759
0760 highzone = TECexRing6 + exclusionWidth;
0761 lowzone = TECexRing6 - exclusionWidth;
0762 } else if (ringnumber == 7) {
0763
0764 highzone = TECexRing7 + exclusionWidth;
0765 lowzone = TECexRing7 - exclusionWidth;
0766 }
0767
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
0779
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
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
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
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) {
0834
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
0877 DEFINE_FWK_MODULE(HitEff);