86 traj = fs->
make<TTree>(
"traj",
"tree of trajectory positions");
101 traj->Branch(
"ResX",&
ResX,
"ResX/F");
108 traj->Branch(
"chi2",&
chi2,
"chi2/F");
109 traj->Branch(
"p",&
p,
"p/F");
110 traj->Branch(
"pT",&
pT,
"pT/F");
112 traj->Branch(
"Id",&
Id,
"Id/i");
113 traj->Branch(
"run",&
run,
"run/i");
120 traj->Branch(
"dedx",&
dedx,
"dedx/F");
140 if (
DEBUG)
cout <<
"beginning analyze from HitEff" << endl;
143 using namespace reco;
146 int run_nr = e.
id().
run();
157 e.
getByLabel(TkTrajCKF,TrajectoryCollectionCKF);
200 const Propagator* thePropagator = prop.product();
235 if (
DEBUG)
cout <<
"number ckf tracks found = " << tracksCKF->size() << endl;
237 if (tracksCKF->size() > 0 && tracksCKF->size()<100) {
238 if (
DEBUG)
cout <<
"starting checking good event with < 100 tracks" << endl;
260 if(e.
getByLabel(
"muonsWitht0Correction",muH)){
262 if(muonsT0.size()!=0) {
268 bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
269 if (hasCaloEnergyInfo)
timeECAL = muonsT0[0].calEnergy().ecal_time;
277 for (vector<Trajectory>::const_iterator itraj = TrajectoryCollectionCKF.product()->begin();
278 itraj != TrajectoryCollectionCKF.product()->end();
282 nHits = itraj->foundHits();
284 chi2 = (itraj->chiSquared()/itraj->ndof());
285 pT =
sqrt( ( itraj->lastMeasurement().updatedState().globalMomentum().x() *
286 itraj->lastMeasurement().updatedState().globalMomentum().x()) +
287 ( itraj->lastMeasurement().updatedState().globalMomentum().y() *
288 itraj->lastMeasurement().updatedState().globalMomentum().y()) );
289 p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
294 std::vector<TrajectoryMeasurement> TMeas=itraj->measurements();
295 vector<TrajectoryMeasurement>::iterator itm;
300 double angleX = -999.;
301 double angleY = -999.;
302 double xglob,yglob,zglob;
304 for (itm=TMeas.begin();itm!=TMeas.end();itm++){
306 theInHit = (*itm).recHit();
308 if(
DEBUG)
cout <<
"theInHit is valid = " << theInHit->isValid() << endl;
310 unsigned int iidd = theInHit->geographicalId().rawId();
312 unsigned int TKlayers =
checkLayer(iidd, tTopo);
313 if (
DEBUG)
cout <<
"TKlayer from trajectory: " << TKlayers <<
" from module = " << iidd <<
" matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
316 if ( TKlayers == 10 || TKlayers == 22 ) {
317 if (
DEBUG)
cout <<
"skipping original TM for TOB 6 or TEC 9" << endl;
322 std::vector<TrajectoryAtInvalidHit> TMs;
340 if (
DEBUG)
cout <<
" found a hit with a missing partner" << endl;
352 bool isValid = theInHit->isValid();
353 bool isLast = (itm==(TMeas.end()-1));
354 bool isLastTOB5 =
true;
356 if (
checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 9 ) isLastTOB5 =
false;
357 else isLastTOB5 =
true;
361 if ( TKlayers==9 && isValid && isLastTOB5 ) {
364 std::vector< BarrelDetLayer*> barrelTOBLayers = measurementTrackerHandle->geometricSearchTracker()->tobLayers() ;
365 const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size()-1];
369 vector<TrajectoryMeasurement>
tmp = theLayerMeasurements->
measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
372 if (
DEBUG)
cout <<
"size of TM from propagation = " << tmp.size() << endl;
379 tob6Hit = tob6TM.recHit();
381 if (tob6Hit->geographicalId().rawId()!=0) {
382 if (
DEBUG)
cout <<
"tob6 hit actually being added to TM vector" << endl;
388 bool isLastTEC8 =
true;
390 if (
checkLayer((++itm)->recHit()->geographicalId().rawId(), tTopo) == 21 ) isLastTEC8 =
false;
391 else isLastTEC8 =
true;
395 if ( TKlayers==21 && isValid && isLastTEC8 ) {
397 std::vector< ForwardDetLayer*> posTecLayers = measurementTrackerHandle->geometricSearchTracker()->posTecLayers() ;
398 const DetLayer* tec9pos = posTecLayers[posTecLayers.size()-1];
399 std::vector< ForwardDetLayer*> negTecLayers = measurementTrackerHandle->geometricSearchTracker()->negTecLayers() ;
400 const DetLayer* tec9neg = negTecLayers[negTecLayers.size()-1];
410 vector<TrajectoryMeasurement>
tmp;
411 if ( tTopo->
tecSide(iidd) == 1 ) {
412 tmp = theLayerMeasurements->
measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
415 if ( tTopo->
tecSide(iidd) == 2 ) {
416 tmp = theLayerMeasurements->
measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
426 tec9Hit = tec9TM.recHit();
428 unsigned int tec9id = tec9Hit->geographicalId().rawId();
429 if (
DEBUG)
cout <<
"tec9id = " << tec9id <<
" is Double sided = " <<
isDoubleSided(tec9id, tTopo) <<
" and 0x3 = " << (tec9id & 0x3) << endl;
431 if (tec9Hit->geographicalId().rawId()!=0) {
432 if (
DEBUG)
cout <<
"tec9 hit actually being added to TM vector" << endl;
448 for(std::vector<TrajectoryAtInvalidHit>::const_iterator TM=TMs.begin();TM!=TMs.end();++TM) {
451 iidd = TM->monodet_id();
452 if (
DEBUG)
cout <<
"setting iidd = " << iidd <<
" before checking efficiency and ";
457 xErr = TM->localErrorX();
458 yErr = TM->localErrorY();
460 angleX = atan( TM->localDxDz() );
461 angleY = atan( TM->localDyDz() );
463 xglob = TM->globalX();
464 yglob = TM->globalY();
465 zglob = TM->globalZ();
475 if (
DEBUG)
cout <<
"Looking at layer under study" << endl;
484 if (
input.size() > 0 ) {
485 if (
DEBUG)
cout <<
"Checking clusters with size = " <<
input.size() << endl;
487 std::vector< std::vector<float> > VCluster_info;
491 unsigned int ClusterId = DSViter->id();
492 if (ClusterId == iidd) {
493 if (
DEBUG)
cout <<
"found (ClusterId == iidd) with ClusterId = " << ClusterId <<
" and iidd = " << iidd << endl;
494 DetId ClusterDetId(ClusterId);
500 float res = (parameters.first.x() - xloc);
512 std::vector< float > cluster_info;
513 cluster_info.push_back(res);
514 cluster_info.push_back(sigma);
515 cluster_info.push_back(parameters.first.x());
516 cluster_info.push_back(
sqrt(parameters.second.xx()));
517 cluster_info.push_back(parameters.first.y());
518 cluster_info.push_back(
sqrt(parameters.second.yy()));
521 VCluster_info.push_back(cluster_info);
523 if (
DEBUG)
cout <<
"Have ID match. residual = " << VCluster_info.back()[0] <<
" res sigma = " << VCluster_info.back()[1] << endl;
524 if (
DEBUG)
cout <<
"trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
525 if (
DEBUG)
cout <<
"hit position = " << parameters.first.x() <<
" hit error = " <<
sqrt(parameters.second.xx()) <<
" trajectory position = " << xloc <<
" traj error = " << xErr << endl;
529 float FinalResSig = 1000.0;
530 float FinalCluster[7]= {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
532 if (
DEBUG)
cout <<
"found clusters > 0" << endl;
535 vector< vector<float> >::iterator
ires;
536 for (ires=VCluster_info.begin(); ires!=VCluster_info.end(); ires++){
537 if (
abs((*ires)[1]) <
abs(FinalResSig)) {
538 FinalResSig = (*ires)[1];
539 for (
unsigned int i = 0;
i<ires->size();
i++) {
540 if (
DEBUG)
cout <<
"filling final cluster. i = " <<
i <<
" before fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
541 FinalCluster[
i] = (*ires)[
i];
542 if (
DEBUG)
cout <<
"filling final cluster. i = " <<
i <<
" after fill FinalCluster[i]=" << FinalCluster[
i] <<
" and (*ires)[i] =" << (*ires)[
i] << endl;
545 if (
DEBUG)
cout <<
"iresidual = " << (*ires)[0] <<
" isigma = " << (*ires)[1] <<
" and FinalRes = " << FinalCluster[0] << endl;
549 FinalResSig = VCluster_info.at(0)[1];
550 for (
unsigned int i = 0;
i<VCluster_info.at(0).size();
i++) {
551 FinalCluster[
i] = VCluster_info.at(0)[
i];
555 VCluster_info.clear();
558 if (
DEBUG)
cout <<
"Final residual in X = " << FinalCluster[0] <<
"+-" << (FinalCluster[0]/FinalResSig) << endl;
559 if (
DEBUG)
cout <<
"Checking location of trajectory: abs(yloc) = " <<
abs(yloc) <<
" abs(xloc) = " <<
abs(xloc) << endl;
560 if (
DEBUG)
cout <<
"Checking location of cluster hit: yloc = " << FinalCluster[4] <<
"+-" << FinalCluster[5] <<
" xloc = " << FinalCluster[2] <<
"+-" << FinalCluster[3] << endl;
561 if (
DEBUG)
cout <<
"Final cluster signal to noise = " << FinalCluster[6] << endl;
563 float exclusionWidth = 0.4;
564 float TOBexclusion = 0.0;
565 float TECexRing5 = -0.89;
566 float TECexRing6 = -0.56;
567 float TECexRing7 = 0.60;
569 int subdetector = ((iidd>>25) & 0x7);
570 int ringnumber = ((iidd>>5) & 0x7);
573 if((TKlayers >= 5 && TKlayers < 11)||((subdetector == 6)&&( (ringnumber >= 5)&&(ringnumber <=7) ))) {
575 float highzone = 0.0;
577 float higherr = yloc + 5.0*yErr;
578 float lowerr = yloc - 5.0*yErr;
579 if(TKlayers >= 5 && TKlayers < 11) {
581 highzone = TOBexclusion + exclusionWidth;
582 lowzone = TOBexclusion - exclusionWidth;
583 }
else if (ringnumber == 5) {
585 highzone = TECexRing5 + exclusionWidth;
586 lowzone = TECexRing5 - exclusionWidth;
587 }
else if (ringnumber == 6) {
589 highzone = TECexRing6 + exclusionWidth;
590 lowzone = TECexRing6 - exclusionWidth;
591 }
else if (ringnumber == 7) {
593 highzone = TECexRing7 + exclusionWidth;
594 lowzone = TECexRing7 - exclusionWidth;
613 if ( SiStripQuality_->getBadApvs(iidd)!=0 ) {
615 if(
DEBUG)
cout <<
"strip is bad from SiStripQuality" << endl;
618 if(
DEBUG)
cout <<
"strip is good from SiStripQuality" << endl;
622 for (
unsigned int ii=0;
ii< fedErrorIds->size();
ii++) {
623 if (iidd == (*fedErrorIds)[
ii].rawId() )
633 ResX = FinalCluster[0];
635 if (FinalResSig != FinalCluster[1])
if (
DEBUG)
cout <<
"Problem with best cluster selection because FinalResSig = " << FinalResSig <<
" and FinalCluster[1] = " << FinalCluster[1] << endl;
642 if (
DEBUG)
cout <<
"before check good" << endl;
644 if ( FinalResSig < 999.0) {
646 if (
DEBUG)
cout <<
"hit being counted as good " << FinalCluster[0] <<
" FinalRecHit " <<
647 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
648 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
653 if (
DEBUG)
cout <<
"hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0] <<
" FinalRecHit " <<
654 iidd <<
" TKlayers " << TKlayers <<
" xloc " << xloc <<
" yloc " << yloc <<
" module " << iidd <<
655 " matched/stereo/rphi = " << ((iidd & 0x3)==0) <<
"/" << ((iidd & 0x3)==1) <<
"/" << ((iidd & 0x3)==2) << endl;
659 if (
DEBUG)
cout <<
" RPhi Error " <<
sqrt(xErr*xErr + yErr*yErr) <<
" ErrorX " << xErr <<
" yErr " << yErr << endl;
660 }
if (
DEBUG)
cout <<
"after good location check" << endl;
661 }
if (
DEBUG)
cout <<
"after list of clusters" << endl;
662 }
if (
DEBUG)
cout <<
"After layers=TKLayers if" << endl;
663 }
if (
DEBUG)
cout <<
"After looping over TrajAtValidHit list" << endl;
664 }
if (
DEBUG)
cout <<
"end TMeasurement loop" << endl;
665 }
if (
DEBUG)
cout <<
"end of trajectories loop" << endl;
670 traj->GetDirectory()->cd();
678 double error =
sqrt(parameters.second.xx() + xerr*xerr);
679 double separation =
abs(parameters.first.x() - xx);
680 double consistency = separation/
error;
686 unsigned int subid=strip.
subdetId();
687 unsigned int layer = 0;
691 if (layer == 1 || layer == 2)
return true;
697 if (layer == 5 || layer == 6)
return true;
702 layer = tTopo->
tidRing(iidd) + 10;
703 if (layer == 11 || layer == 12)
return true;
708 layer = tTopo->
tecRing(iidd) + 13 ;
709 if (layer == 14 || layer == 15 || layer == 18)
return true;
717 unsigned int partner_iidd = 0;
718 bool found2DPartner =
false;
720 if ((iidd & 0x3)==1) partner_iidd = iidd+1;
721 if ((iidd & 0x3)==2) partner_iidd = iidd-1;
724 for (std::vector<TrajectoryMeasurement>::const_iterator iTM=traj.begin(); iTM!=traj.end(); ++iTM) {
725 if (iTM->recHit()->geographicalId().rawId()==partner_iidd) {
726 found2DPartner =
true;
729 return found2DPartner;
734 unsigned int subid=strip.
subdetId();
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
unsigned int tibLayer(const DetId &id) const
unsigned int tidRing(const DetId &id) const
virtual void analyze(const edm::Event &e, const edm::EventSetup &c)
#define DEFINE_FWK_MODULE(type)
unsigned int tecRing(const DetId &id) const
ring id
std::vector< Track > TrackCollection
collection of Tracks
int bunchCrossing() const
unsigned int tidWheel(const DetId &id) const
HitEff(const edm::ParameterSet &conf)
data_type const * const_iterator
std::pair< LocalPoint, LocalError > LocalValues
std::vector< Muon > MuonCollection
collection of Muon objects
virtual LocalValues localParameters(const T &, const GeomDetUnit &) const =0
float signalOverNoise() const
bool isDoubleSided(unsigned int iidd, const TrackerTopology *tTopo) const
int nDof
number of muon stations used
unsigned int trajHitValid
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
double checkConsistency(const StripClusterParameterEstimator::LocalValues ¶meters, double xx, double xerr)
unsigned int checkLayer(unsigned int iidd, const TrackerTopology *tTopo)
T const * product() const
T const * product() const
unsigned int SiStripQualBad
std::vector< std::vector< double > > tmp
T * make() const
make new ROOT object
unsigned int tecWheel(const DetId &id) const
tuple AnalyticalPropagator
unsigned int tobLayer(const DetId &id) const
unsigned int tecSide(const DetId &id) const