169 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"beginning analyze from HitEff" << endl;
172 using namespace reco;
175 int run_nr = e.
id().
run();
184 if (lumiScalers->begin() != lumiScalers->end()) {
185 instLumi = lumiScalers->begin()->instantLumi();
186 PU = lumiScalers->begin()->pileup();
220 es.get<
TkStripCPERecord>().
get(
"StripCPEfromTrackAngle", parameterestimator);
289 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"number ckf tracks found = " << tracksCKF->size() << endl;
291 if (!tracksCKF->empty()) {
297 #ifdef ExtendedCALIBTree 300 if (e.
getByLabel(
"dedxMedianCTF", dEdxUncalibHandle)) {
304 dedx = dEdxTrackUncalib[itTrack].dEdx();
305 dedxNOM = dEdxTrackUncalib[itTrack].numberOfMeasurements();
307 dedx = -999.0; dedxNOM = -999;
314 if(e.
getByLabel(
"muonsWitht0Correction",muH)){
316 if(!muonsT0.empty()) {
320 timeDTDOF = mt0.
nDof;
322 bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
323 if (hasCaloEnergyInfo) timeECAL = muonsT0[0].calEnergy().ecal_time;
326 timeDT = -999.0; timeDTErr = -999.0; timeDTDOF = -999; timeECAL = -999.0;
334 it!=trajTrackAssociationHandle->
end();
341 nHits = itraj->foundHits();
342 #ifdef ExtendedCALIBTree 343 nLostHits = itraj->lostHits();
344 chi2 = (itraj->chiSquared()/itraj->ndof());
345 p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
347 pT =
sqrt( ( itraj->lastMeasurement().updatedState().globalMomentum().x() *
348 itraj->lastMeasurement().updatedState().globalMomentum().x()) +
349 ( itraj->lastMeasurement().updatedState().globalMomentum().y() *
350 itraj->lastMeasurement().updatedState().globalMomentum().y()) );
357 std::vector<TrajectoryMeasurement> TMeas=itraj->measurements();
358 vector<TrajectoryMeasurement>::iterator itm;
363 double angleX = -999.;
364 double angleY = -999.;
365 double xglob,yglob,zglob;
369 bool hasMissingHits=
false;
370 for (itm=TMeas.begin();itm!=TMeas.end();itm++){
371 auto theHit = (*itm).recHit();
372 if(theHit->getType()==TrackingRecHit::Type::missing) hasMissingHits=
true;
379 for (itm=TMeas.begin();itm!=TMeas.end();itm++){
380 auto theInHit = (*itm).recHit();
382 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"theInHit is valid = " << theInHit->isValid() << endl;
384 unsigned int iidd = theInHit->geographicalId().rawId();
386 unsigned int TKlayers =
checkLayer(iidd, tTopo);
387 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"TKlayer from trajectory: " << TKlayers <<
" from module = " << iidd <<
" matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
392 bool isFirstMeas = (itm==(TMeas.end()-1));
393 bool isLastMeas = (itm==(TMeas.begin()));
404 if ( TKlayers == 10 || TKlayers == 22 ) {
405 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"skipping original TM for TOB 6 or TEC 9" << endl;
410 std::vector<TrajectoryAtInvalidHit> TMs;
428 LogDebug(
"SiStripHitEfficiency:HitEff") <<
" found a hit with a missing partner" << endl;
440 bool isValid = theInHit->isValid();
441 bool isLast = (itm==(TMeas.end()-1));
442 bool isLastTOB5 =
true;
444 if (
checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9 ) isLastTOB5 =
false;
445 else isLastTOB5 =
true;
449 if ( TKlayers==9 && isValid && isLastTOB5 ) {
452 std::vector< BarrelDetLayer const*> barrelTOBLayers = measurementTrackerHandle->geometricSearchTracker()->tobLayers() ;
453 const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size()-1];
457 vector<TrajectoryMeasurement>
tmp = theLayerMeasurements->
measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
460 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"size of TM from propagation = " << tmp.size() << endl;
466 const auto& tob6Hit = tob6TM.
recHit();
468 if (tob6Hit->geographicalId().rawId()!=0) {
469 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"tob6 hit actually being added to TM vector" << endl;
475 bool isLastTEC8 =
true;
477 if (
checkLayer((++itm)->
recHit()->geographicalId().rawId(), tTopo) == 21 ) isLastTEC8 =
false;
478 else isLastTEC8 =
true;
482 if ( TKlayers==21 && isValid && isLastTEC8 ) {
484 std::vector< const ForwardDetLayer*> posTecLayers = measurementTrackerHandle->geometricSearchTracker()->posTecLayers() ;
485 const DetLayer* tec9pos = posTecLayers[posTecLayers.size()-1];
486 std::vector< const ForwardDetLayer*> negTecLayers = measurementTrackerHandle->geometricSearchTracker()->negTecLayers() ;
487 const DetLayer* tec9neg = negTecLayers[negTecLayers.size()-1];
497 vector<TrajectoryMeasurement>
tmp;
498 if ( tTopo->
tecSide(iidd) == 1 ) {
499 tmp = theLayerMeasurements->
measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
502 if ( tTopo->
tecSide(iidd) == 2 ) {
503 tmp = theLayerMeasurements->
measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
512 const auto& tec9Hit = tec9TM.
recHit();
514 unsigned int tec9id = tec9Hit->geographicalId().rawId();
515 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"tec9id = " << tec9id <<
" is Double sided = " <<
isDoubleSided(tec9id, tTopo) <<
" and 0x3 = " << (tec9id & 0x3) << endl;
517 if (tec9Hit->geographicalId().rawId()!=0) {
518 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"tec9 hit actually being added to TM vector" << endl;
534 for(std::vector<TrajectoryAtInvalidHit>::const_iterator TM=TMs.begin();TM!=TMs.end();++TM) {
537 iidd = TM->monodet_id();
538 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"setting iidd = " << iidd <<
" before checking efficiency and ";
543 angleX = atan( TM->localDxDz() );
544 angleY = atan( TM->localDyDz() );
548 xglob = TM->globalX();
549 yglob = TM->globalY();
550 zglob = TM->globalZ();
551 xErr = TM->localErrorX();
552 yErr = TM->localErrorY();
565 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Looking at layer under study" << endl;
574 if (!input.
empty() ) {
575 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Checking clusters with size = " << input.
size() << endl;
577 std::vector< std::vector<float> > VCluster_info;
581 unsigned int ClusterId = DSViter->id();
582 if (ClusterId == iidd) {
583 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"found (ClusterId == iidd) with ClusterId = " << ClusterId <<
" and iidd = " << iidd << endl;
584 DetId ClusterDetId(ClusterId);
595 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
596 std::array<const float, 4>
const & parameterTrap = (*trapezoidalBounds).
parameters();
597 hbedge = parameterTrap[0];
598 htedge = parameterTrap[1];
599 hapoth = parameterTrap[3];
600 uylfac = (htedge-hbedge)/(htedge+hbedge)/hapoth;
601 uxlden = 1 + yloc*uylfac;
605 if( TrajStrip==-1 ) {
608 TrajStrip = xloc/pitch + nstrips/2.0;
611 float TrajLocXMid = xloc / (1 + (htedge-hbedge)*yloc/(htedge+hbedge)/hapoth) ;
612 TrajStrip = TrajLocXMid/pitch + nstrips/2.0;
621 float res = (parameters.first.x() - xloc);
630 res = parameters.first.x() - xloc/uxlden ;
631 sigma =
abs(res) /
sqrt(parameters.second.xx() + xErr*xErr/uxlden/uxlden + yErr*yErr*xloc*xloc*uylfac*uylfac/uxlden/uxlden/uxlden/uxlden);
638 std::vector< float > cluster_info;
639 cluster_info.push_back(res);
640 cluster_info.push_back(sigma);
641 cluster_info.push_back(parameters.first.x());
642 cluster_info.push_back(
sqrt(parameters.second.xx()));
643 cluster_info.push_back(parameters.first.y());
644 cluster_info.push_back(
sqrt(parameters.second.yy()));
645 cluster_info.push_back( clusterInfo.signalOverNoise() );
647 VCluster_info.push_back(cluster_info);
649 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Have ID match. residual = " << VCluster_info.back()[0] <<
" res sigma = " << VCluster_info.back()[1] << endl;
650 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
651 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"hit position = " << parameters.first.x() <<
" hit error = " <<
sqrt(parameters.second.xx()) <<
" trajectory position = " << xloc <<
" traj error = " << xErr << endl;
655 float FinalResSig = 1000.0;
656 float FinalCluster[7]= {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
658 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"found clusters > 0" << endl;
661 vector< vector<float> >::iterator
ires;
662 for (ires=VCluster_info.begin(); ires!=VCluster_info.end(); ires++){
663 if (
abs((*ires)[1]) <
abs(FinalResSig)) {
664 FinalResSig = (*ires)[1];
665 for (
unsigned int i = 0;
i<ires->size();
i++) {
666 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"filling final cluster. i = " <<
i <<
" before fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
667 FinalCluster[
i] = (*ires)[
i];
668 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"filling final cluster. i = " <<
i <<
" after fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
671 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"iresidual = " << (*ires)[0] <<
" isigma = " << (*ires)[1] <<
" and FinalRes = " << FinalCluster[0] << endl;
675 FinalResSig = VCluster_info.at(0)[1];
676 for (
unsigned int i = 0;
i<VCluster_info.at(0).size();
i++) {
677 FinalCluster[
i] = VCluster_info.at(0)[
i];
681 VCluster_info.clear();
684 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Final residual in X = " << FinalCluster[0] <<
"+-" << (FinalCluster[0]/FinalResSig) << endl;
685 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Checking location of trajectory: abs(yloc) = " <<
abs(yloc) <<
" abs(xloc) = " <<
abs(xloc) << endl;
686 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Checking location of cluster hit: yloc = " << FinalCluster[4] <<
"+-" << FinalCluster[5] <<
" xloc = " << FinalCluster[2] <<
"+-" << FinalCluster[3] << endl;
687 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Final cluster signal to noise = " << FinalCluster[6] << endl;
689 float exclusionWidth = 0.4;
690 float TOBexclusion = 0.0;
691 float TECexRing5 = -0.89;
692 float TECexRing6 = -0.56;
693 float TECexRing7 = 0.60;
696 int ringnumber = ((iidd>>5) & 0x7);
699 if((TKlayers >= 5 && TKlayers < 11)||((subdetector == 6)&&( (ringnumber >= 5)&&(ringnumber <=7) ))) {
701 float highzone = 0.0;
703 float higherr = yloc + 5.0*yErr;
704 float lowerr = yloc - 5.0*yErr;
705 if(TKlayers >= 5 && TKlayers < 11) {
707 highzone = TOBexclusion + exclusionWidth;
708 lowzone = TOBexclusion - exclusionWidth;
709 }
else if (ringnumber == 5) {
711 highzone = TECexRing5 + exclusionWidth;
712 lowzone = TECexRing5 - exclusionWidth;
713 }
else if (ringnumber == 6) {
715 highzone = TECexRing6 + exclusionWidth;
716 lowzone = TECexRing6 - exclusionWidth;
717 }
else if (ringnumber == 7) {
719 highzone = TECexRing7 + exclusionWidth;
720 lowzone = TECexRing7 - exclusionWidth;
745 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"strip is bad from SiStripQuality" << endl;
748 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"strip is good from SiStripQuality" << endl;
752 for (
unsigned int ii=0;
ii< fedErrorIds->
size();
ii++) {
753 if (iidd == (*fedErrorIds)[
ii].rawId() )
761 ResX = FinalCluster[0];
763 if (FinalResSig != FinalCluster[1])
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"Problem with best cluster selection because FinalResSig = " << FinalResSig <<
" and FinalCluster[1] = " << FinalCluster[1] << endl;
772 if(commonModeDigis.
isValid() && TrajStrip>=0 && TrajStrip<=768) {
774 if(digiframe != commonModeDigis->end())
775 if( (
unsigned) TrajStrip/128 < digiframe->data.
size())
776 commonMode = digiframe->data.at(TrajStrip/128).adc();
779 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"before check good" << endl;
781 if ( FinalResSig < 999.0) {
783 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"hit being counted as good " << FinalCluster[0] <<
" FinalRecHit " <<
784 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
785 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
790 LogDebug(
"SiStripHitEfficiency:HitEff") <<
"hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0] <<
" FinalRecHit " <<
791 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
792 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
796 LogDebug(
"SiStripHitEfficiency:HitEff") <<
" RPhi Error " <<
sqrt(xErr*xErr + yErr*yErr) <<
" ErrorX " << xErr <<
" yErr " << yErr << endl;
797 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"after good location check" << endl;
798 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"after list of clusters" << endl;
799 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"After layers=TKLayers if" << endl;
800 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"After looping over TrajAtValidHit list" << endl;
801 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"end TMeasurement loop" << endl;
802 }
LogDebug(
"SiStripHitEfficiency:HitEff") <<
"end of trajectories loop" << endl;
EventNumber_t event() const
std::vector< TrajectoryMeasurement > measurements(const DetLayer &layer, const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
bool check2DPartner(unsigned int iidd, const std::vector< TrajectoryMeasurement > &traj)
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
ConstRecHitPointer const & recHit() const
virtual const std::array< const float, 4 > parameters() const
friend struct const_iterator
const_iterator end() const
last iterator over the map (read only)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
bool useAllHitsFromTracksWithMissingHits_
std::vector< Track > TrackCollection
collection of Tracks
int bunchCrossing() const
const Bounds & bounds() const
std::pair< LocalPoint, LocalError > LocalValues
data_type const * const_iterator
key_type key() const
Accessor for product key.
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
const Plane & surface() const
The nominal surface of the GeomDet.
std::vector< Muon > MuonCollection
collection of Muon objects
const edm::EDGetTokenT< MeasurementTrackerEvent > trackerEvent_token_
static std::string const input
virtual float width() const =0
const TrackerGeomDet * idToDetUnit(DetId) const override
Return the pointer to the GeomDetUnit corresponding to a given DetId.
bool isDoubleSided(unsigned int iidd, const TrackerTopology *tTopo) const
short getBadApvs(const uint32_t &detid) const
const edm::EDGetTokenT< DetIdCollection > digis_token_
int nDof
number of muon stations used
Abs< T >::type abs(const T &t)
unsigned int trajHitValid
const edm::EDGetTokenT< TrajTrackAssociationCollection > trajTrackAsso_token_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
size_type size() const
Return the number of contained DetSets.
double checkConsistency(const StripClusterParameterEstimator::LocalValues ¶meters, double xx, double xerr)
T const * product() const
virtual int nstrips() const =0
unsigned int trackMultiplicityCut_
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
unsigned int checkLayer(unsigned int iidd, const TrackerTopology *tTopo)
const edm::EDGetTokenT< reco::TrackCollection > combinatorialTracks_token_
unsigned int SiStripQualBad
std::vector< std::vector< double > > tmp
const_iterator begin() const
first iterator over the map (read only)
const edm::EDGetTokenT< edmNew::DetSetVector< SiStripCluster > > clusters_token_
collection_type::const_iterator const_iterator
const edm::EDGetTokenT< LumiScalersCollection > scalerToken_
const edm::EDGetTokenT< std::vector< Trajectory > > trajectories_token_
T const * product() const
const edm::EDGetTokenT< edm::DetSetVector< SiStripRawDigi > > commonModeToken_
const_iterator begin(bool update=false) const
unsigned int tecSide(const DetId &id) const