141 if (
DEBUG)
cout <<
"beginning analyze from HitEff" << endl;
144 using namespace reco;
147 int run_nr = e.
id().
run();
158 e.
getByLabel(TkTrajCKF,TrajectoryCollectionCKF);
173 es.get<
TkStripCPERecord>().
get(
"StripCPEfromTrackAngle", parameterestimator);
192 e.
getByLabel(
"MeasurementTrackerEvent", measurementTrackerEvent);
199 const Propagator* thePropagator = prop.product();
234 if (
DEBUG)
cout <<
"number ckf tracks found = " << tracksCKF->size() << endl;
236 if (tracksCKF->size() > 0 && tracksCKF->size()<100) {
237 if (
DEBUG)
cout <<
"starting checking good event with < 100 tracks" << endl;
259 if(e.
getByLabel(
"muonsWitht0Correction",muH)){
261 if(muonsT0.size()!=0) {
267 bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
268 if (hasCaloEnergyInfo)
timeECAL = muonsT0[0].calEnergy().ecal_time;
276 for (vector<Trajectory>::const_iterator itraj = TrajectoryCollectionCKF.product()->begin();
277 itraj != TrajectoryCollectionCKF.product()->end();
281 nHits = itraj->foundHits();
283 chi2 = (itraj->chiSquared()/itraj->ndof());
284 pT =
sqrt( ( itraj->lastMeasurement().updatedState().globalMomentum().x() *
285 itraj->lastMeasurement().updatedState().globalMomentum().x()) +
286 ( itraj->lastMeasurement().updatedState().globalMomentum().y() *
287 itraj->lastMeasurement().updatedState().globalMomentum().y()) );
288 p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
293 std::vector<TrajectoryMeasurement> TMeas=itraj->measurements();
294 vector<TrajectoryMeasurement>::iterator itm;
299 double angleX = -999.;
300 double angleY = -999.;
301 double xglob,yglob,zglob;
303 for (itm=TMeas.begin();itm!=TMeas.end();itm++){
304 auto theInHit = (*itm).recHit();
306 if(
DEBUG)
cout <<
"theInHit is valid = " << theInHit->isValid() << endl;
308 unsigned int iidd = theInHit->geographicalId().rawId();
310 unsigned int TKlayers =
checkLayer(iidd, tTopo);
311 if (
DEBUG)
cout <<
"TKlayer from trajectory: " << TKlayers <<
" from module = " << iidd <<
" matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
314 if ( TKlayers == 10 || TKlayers == 22 ) {
315 if (
DEBUG)
cout <<
"skipping original TM for TOB 6 or TEC 9" << endl;
320 std::vector<TrajectoryAtInvalidHit> TMs;
338 if (
DEBUG)
cout <<
" found a hit with a missing partner" << endl;
350 bool isValid = theInHit->isValid();
351 bool isLast = (itm==(TMeas.end()-1));
352 bool isLastTOB5 =
true;
354 if (
checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9 ) isLastTOB5 =
false;
355 else isLastTOB5 =
true;
359 if ( TKlayers==9 && isValid && isLastTOB5 ) {
362 std::vector< BarrelDetLayer const*> barrelTOBLayers = measurementTrackerHandle->geometricSearchTracker()->tobLayers() ;
363 const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size()-1];
367 vector<TrajectoryMeasurement>
tmp = theLayerMeasurements->
measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
370 if (
DEBUG)
cout <<
"size of TM from propagation = " << tmp.size() << endl;
376 auto tob6Hit = tob6TM.
recHit();
378 if (tob6Hit->geographicalId().rawId()!=0) {
379 if (
DEBUG)
cout <<
"tob6 hit actually being added to TM vector" << endl;
385 bool isLastTEC8 =
true;
387 if (
checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 21 ) isLastTEC8 =
false;
388 else isLastTEC8 =
true;
392 if ( TKlayers==21 && isValid && isLastTEC8 ) {
394 std::vector< const ForwardDetLayer*> posTecLayers = measurementTrackerHandle->geometricSearchTracker()->posTecLayers() ;
395 const DetLayer* tec9pos = posTecLayers[posTecLayers.size()-1];
396 std::vector< const ForwardDetLayer*> negTecLayers = measurementTrackerHandle->geometricSearchTracker()->negTecLayers() ;
397 const DetLayer* tec9neg = negTecLayers[negTecLayers.size()-1];
407 vector<TrajectoryMeasurement>
tmp;
408 if ( tTopo->
tecSide(iidd) == 1 ) {
409 tmp = theLayerMeasurements->
measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
412 if ( tTopo->
tecSide(iidd) == 2 ) {
413 tmp = theLayerMeasurements->
measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
422 auto tec9Hit = tec9TM.
recHit();
424 unsigned int tec9id = tec9Hit->geographicalId().rawId();
425 if (
DEBUG)
cout <<
"tec9id = " << tec9id <<
" is Double sided = " <<
isDoubleSided(tec9id, tTopo) <<
" and 0x3 = " << (tec9id & 0x3) << endl;
427 if (tec9Hit->geographicalId().rawId()!=0) {
428 if (
DEBUG)
cout <<
"tec9 hit actually being added to TM vector" << endl;
444 for(std::vector<TrajectoryAtInvalidHit>::const_iterator TM=TMs.begin();TM!=TMs.end();++TM) {
447 iidd = TM->monodet_id();
448 if (
DEBUG)
cout <<
"setting iidd = " << iidd <<
" before checking efficiency and ";
453 xErr = TM->localErrorX();
454 yErr = TM->localErrorY();
456 angleX = atan( TM->localDxDz() );
457 angleY = atan( TM->localDyDz() );
459 xglob = TM->globalX();
460 yglob = TM->globalY();
461 zglob = TM->globalZ();
471 if (
DEBUG)
cout <<
"Looking at layer under study" << endl;
480 if (
input.size() > 0 ) {
481 if (
DEBUG)
cout <<
"Checking clusters with size = " <<
input.size() << endl;
483 std::vector< std::vector<float> > VCluster_info;
487 unsigned int ClusterId = DSViter->id();
488 if (ClusterId == iidd) {
489 if (
DEBUG)
cout <<
"found (ClusterId == iidd) with ClusterId = " << ClusterId <<
" and iidd = " << iidd << endl;
490 DetId ClusterDetId(ClusterId);
496 float res = (parameters.first.x() - xloc);
508 std::vector< float > cluster_info;
509 cluster_info.push_back(res);
510 cluster_info.push_back(sigma);
511 cluster_info.push_back(parameters.first.x());
512 cluster_info.push_back(
sqrt(parameters.second.xx()));
513 cluster_info.push_back(parameters.first.y());
514 cluster_info.push_back(
sqrt(parameters.second.yy()));
517 VCluster_info.push_back(cluster_info);
519 if (
DEBUG)
cout <<
"Have ID match. residual = " << VCluster_info.back()[0] <<
" res sigma = " << VCluster_info.back()[1] << endl;
520 if (
DEBUG)
cout <<
"trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
521 if (
DEBUG)
cout <<
"hit position = " << parameters.first.x() <<
" hit error = " <<
sqrt(parameters.second.xx()) <<
" trajectory position = " << xloc <<
" traj error = " << xErr << endl;
525 float FinalResSig = 1000.0;
526 float FinalCluster[7]= {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
528 if (
DEBUG)
cout <<
"found clusters > 0" << endl;
531 vector< vector<float> >::iterator
ires;
532 for (ires=VCluster_info.begin(); ires!=VCluster_info.end(); ires++){
533 if (
abs((*ires)[1]) <
abs(FinalResSig)) {
534 FinalResSig = (*ires)[1];
535 for (
unsigned int i = 0;
i<ires->size();
i++) {
536 if (
DEBUG)
cout <<
"filling final cluster. i = " <<
i <<
" before fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
537 FinalCluster[
i] = (*ires)[
i];
538 if (
DEBUG)
cout <<
"filling final cluster. i = " <<
i <<
" after fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
541 if (
DEBUG)
cout <<
"iresidual = " << (*ires)[0] <<
" isigma = " << (*ires)[1] <<
" and FinalRes = " << FinalCluster[0] << endl;
545 FinalResSig = VCluster_info.at(0)[1];
546 for (
unsigned int i = 0;
i<VCluster_info.at(0).size();
i++) {
547 FinalCluster[
i] = VCluster_info.at(0)[
i];
551 VCluster_info.clear();
554 if (
DEBUG)
cout <<
"Final residual in X = " << FinalCluster[0] <<
"+-" << (FinalCluster[0]/FinalResSig) << endl;
555 if (
DEBUG)
cout <<
"Checking location of trajectory: abs(yloc) = " <<
abs(yloc) <<
" abs(xloc) = " <<
abs(xloc) << endl;
556 if (
DEBUG)
cout <<
"Checking location of cluster hit: yloc = " << FinalCluster[4] <<
"+-" << FinalCluster[5] <<
" xloc = " << FinalCluster[2] <<
"+-" << FinalCluster[3] << endl;
557 if (
DEBUG)
cout <<
"Final cluster signal to noise = " << FinalCluster[6] << endl;
559 float exclusionWidth = 0.4;
560 float TOBexclusion = 0.0;
561 float TECexRing5 = -0.89;
562 float TECexRing6 = -0.56;
563 float TECexRing7 = 0.60;
565 int subdetector = ((iidd>>25) & 0x7);
566 int ringnumber = ((iidd>>5) & 0x7);
569 if((TKlayers >= 5 && TKlayers < 11)||((subdetector == 6)&&( (ringnumber >= 5)&&(ringnumber <=7) ))) {
571 float highzone = 0.0;
573 float higherr = yloc + 5.0*yErr;
574 float lowerr = yloc - 5.0*yErr;
575 if(TKlayers >= 5 && TKlayers < 11) {
577 highzone = TOBexclusion + exclusionWidth;
578 lowzone = TOBexclusion - exclusionWidth;
579 }
else if (ringnumber == 5) {
581 highzone = TECexRing5 + exclusionWidth;
582 lowzone = TECexRing5 - exclusionWidth;
583 }
else if (ringnumber == 6) {
585 highzone = TECexRing6 + exclusionWidth;
586 lowzone = TECexRing6 - exclusionWidth;
587 }
else if (ringnumber == 7) {
589 highzone = TECexRing7 + exclusionWidth;
590 lowzone = TECexRing7 - exclusionWidth;
609 if ( SiStripQuality_->getBadApvs(iidd)!=0 ) {
611 if(
DEBUG)
cout <<
"strip is bad from SiStripQuality" << endl;
614 if(
DEBUG)
cout <<
"strip is good from SiStripQuality" << endl;
618 for (
unsigned int ii=0;
ii< fedErrorIds->size();
ii++) {
619 if (iidd == (*fedErrorIds)[
ii].rawId() )
629 ResX = FinalCluster[0];
631 if (FinalResSig != FinalCluster[1])
if (
DEBUG)
cout <<
"Problem with best cluster selection because FinalResSig = " << FinalResSig <<
" and FinalCluster[1] = " << FinalCluster[1] << endl;
638 if (
DEBUG)
cout <<
"before check good" << endl;
640 if ( FinalResSig < 999.0) {
642 if (
DEBUG)
cout <<
"hit being counted as good " << FinalCluster[0] <<
" FinalRecHit " <<
643 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
644 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
649 if (
DEBUG)
cout <<
"hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0] <<
" FinalRecHit " <<
650 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
651 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
655 if (
DEBUG)
cout <<
" RPhi Error " <<
sqrt(xErr*xErr + yErr*yErr) <<
" ErrorX " << xErr <<
" yErr " << yErr << endl;
656 }
if (
DEBUG)
cout <<
"after good location check" << endl;
657 }
if (
DEBUG)
cout <<
"after list of clusters" << endl;
658 }
if (
DEBUG)
cout <<
"After layers=TKLayers if" << endl;
659 }
if (
DEBUG)
cout <<
"After looping over TrajAtValidHit list" << endl;
660 }
if (
DEBUG)
cout <<
"end TMeasurement loop" << endl;
661 }
if (
DEBUG)
cout <<
"end of trajectories loop" << endl;
T getParameter(std::string const &) const
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
ConstRecHitPointer const & recHit() const
std::vector< Track > TrackCollection
collection of Tracks
int bunchCrossing() const
data_type const * const_iterator
std::vector< Muon > MuonCollection
collection of Muon objects
static std::string const input
float signalOverNoise() const
bool isDoubleSided(unsigned int iidd, const TrackerTopology *tTopo) const
int nDof
number of muon stations used
Abs< T >::type abs(const T &t)
unsigned int trajHitValid
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
double checkConsistency(const StripClusterParameterEstimator::LocalValues ¶meters, double xx, double xerr)
T const * product() const
unsigned int checkLayer(unsigned int iidd, const TrackerTopology *tTopo)
T const * product() const
std::pair< LocalPoint, LocalError > LocalValues
unsigned int SiStripQualBad
std::vector< std::vector< double > > tmp
tuple AnalyticalPropagator
unsigned int tecSide(const DetId &id) const