36 return (hit_1->globalPosition().mag2()>hit_2->globalPosition().mag2());
39 return (fabs(hit_1->globalPosition().z())>fabs(hit_2->globalPosition().z()));
61 <<
"SETFilter destructor called"<<endl;
89 <<
"Unrecognized module type in incrementChamberCounters";
98 std::vector < SeedCandidate> & validSegmentsSet_out){
100 validSegmentsSet_out.clear();
107 bool validStep =
true;
108 std::vector <double> chi2AllCombinations(validSegmentsSet_in.size());
109 std::vector < TrajectoryStateOnSurface > lastUpdatedTSOS_Vect(validSegmentsSet_in.size());
111 for(
unsigned int iSet = 0; iSet<validSegmentsSet_in.size(); ++iSet){
114 CLHEP::Hep3Vector origin (0.,0.,0.);
117 chi2AllCombinations[iSet] =
findMinChi2(iSet, origin, validSegmentsSet_in[iSet], lastUpdatedTSOS_Vect,
118 trajectoryMeasurementsInTheSet_tmp);
121 std::vector < double >::iterator itMin = min_element(chi2AllCombinations.begin(),chi2AllCombinations.end());
123 int positionMin = itMin - chi2AllCombinations.begin();
126 validSegmentsSet_out.push_back(validSegmentsSet_in[positionMin]);
134 bool validTrajectory =
true;
137 finalCandidate.clear();
159 validTrajectory =
false;
162 return validTrajectory;
173 for(
int iMeas = measurements_segments.size() - 1; iMeas>-1;--iMeas){
176 for(
unsigned int jMeas = 0; jMeas<measurements_segments[iMeas].recHit()->transientHits().size();++jMeas){
177 if(measurements_segments[iMeas].recHit()->transientHits().at(jMeas)->transientHits().size()>1){
180 for(
unsigned int kMeas = 0;kMeas<measurements_segments[iMeas].recHit()->transientHits().at(jMeas)->transientHits().size();
182 sortedHits.push_back(
183 measurements_segments[iMeas].recHit()->transientHits().
at(jMeas)->transientHits().
at(kMeas));
187 sortedHits = measurements_segments[iMeas].recHit()->transientHits();
192 hitContainer.insert(hitContainer.end(),sortedHits.begin(),sortedHits.end());
196 FreeTrajectoryState ftsStart = *(measurements_segments.at(measurements_segments.size()-1).forwardPredictedState().freeState());
200 DetId detId_last = hitContainer[0]->hit()->geographicalId();
205 firstTSOS = tSOSDest;
223 for(
unsigned int iMeas = 0; iMeas<measurements_segments.size();++iMeas){
224 hitContainer.push_back(measurements_segments[iMeas].recHit());
228 for(
int iMeas = measurements_segments.size() - 1; iMeas>-1;--iMeas){
229 hitContainer.push_back(measurements_segments[iMeas].recHit());
234 firstTSOS = measurements_segments.at(0).forwardPredictedState();
240 const CLHEP::Hep3Vector& r3T,
244 bool detailedOutput){
250 trajectoryMeasurementsInTheSet.clear();
262 ftsStart.setCurvilinearError(cov);
266 double chi2_loc = 0.;
267 for(
unsigned int iMeas = 0; iMeas <muonCandidate.
theSet.size(); ++iMeas){
269 DetId detId = muonRecHit->hit()->geographicalId();
279 lastUpdatedTSOS = tSOSDest;
285 chi2_loc = 9999999999.;
290 LocalPoint locHitPos = muonRecHit->localPosition();
291 LocalError locHitErr = muonRecHit->localPositionError();
292 const GlobalPoint globPropPos = ftsStart.position();
297 CLHEP::HepMatrix dist(1,2);
298 double chi2_intermed = -9;
300 dist(1,1) = locPropPos.
x() - locHitPos.
x();
301 dist(1,2) = locPropPos.
y() - locHitPos.
y();
302 CLHEP::HepMatrix IC(2,2);
303 IC(1,1) = locHitErr.
xx();
304 IC(2,1) = locHitErr.
xy();
305 IC(2,2) = locHitErr.
yy();
309 if(4!=muonRecHit->hit()->dimension()){
310 for(
int iE = 1;iE<3;++iE){
312 if ( fabs(IC(iE,iE)) < 0.00000001){
313 IC(iE,iE) = 62500./12.;
322 chi2_intermed =
pow(dist(1,1),2)*IC(1,1) + 2.*dist(1,1)*dist(1,2)*IC(1,2) +
pow(dist(1,2),2)*IC(2,2);
324 chi2_intermed = 9999999999.;
326 chi2_loc += chi2_intermed;
330 DetId detId = muonRecHit->hit()->geographicalId();
354 std::vector < TrajectoryStateOnSurface > &lastUpdatedTSOS_Vect,
361 bool detailedOutput =
false;
363 double cosTheta = muonCandidate.
momentum.cosTheta();
364 double theta = acos(cosTheta);
366 double pMag = muonCandidate.
momentum.mag();
378 double invP = 1./pMag;
385 const double reflect = 1;
386 const double expand = -0.5;
387 const double contract = 0.25;
391 std::vector <CLHEP::Hep3Vector> feet(nDim+1);
392 std::vector <double> chi2Feet(nDim+1);
393 std::vector <TrajectoryStateOnSurface*> lastUpdatedTSOS_pointer(nDim+1);
397 CLHEP::Hep3Vector foot1(invP, theta, phi);
400 trajectoryMeasurementsInTheSet, detailedOutput);
403 std::vector <CLHEP::Hep3Vector> morePoints =
405 feet[1] = morePoints[0];
406 feet[2] = morePoints[1];
407 feet[3] = morePoints[2];
410 for(
unsigned int iFoot = 1; iFoot<feet.size();++iFoot){
412 chi2Feet[iFoot] =
chi2AtSpecificStep(feet[iFoot], r3T, muonCandidate, lastUpdatedTSOS,
413 trajectoryMeasurementsInTheSet, detailedOutput);
417 unsigned int high, second_high, low;
427 trajectoryMeasurementsInTheSet, detailedOutput);
430 if(chi2Feet[high] <chi2Feet[low]){
434 trajectoryMeasurementsInTheSet, detailedOutput);
438 else if( chi2Feet[high] > chi2Feet[second_high]){
439 double worstChi2 = chi2Feet[high];
443 trajectoryMeasurementsInTheSet, detailedOutput);
446 if(chi2Feet[high] <worstChi2){
448 for(
unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
450 chi2Feet[iFoot] =
chi2AtSpecificStep(feet[iFoot], r3T, muonCandidate, lastUpdatedTSOS,
451 trajectoryMeasurementsInTheSet, detailedOutput);
461 int bestFitElement = min_element(chi2Feet.begin(),chi2Feet.end()) - chi2Feet.begin();
464 detailedOutput =
true;
465 chi2Feet[bestFitElement] =
chi2AtSpecificStep(feet[bestFitElement], r3T, muonCandidate, lastUpdatedTSOS,
466 trajectoryMeasurementsInTheSet, detailedOutput);
467 minChi2 = chi2Feet[bestFitElement];
469 double pMag_updated = 1./feet[bestFitElement].x();
470 double sin_theta =
sin(feet[bestFitElement].
y());
471 double cos_theta =
cos(feet[bestFitElement].
y());
472 double sin_phi =
sin(feet[bestFitElement].
z());
473 double cos_phi =
cos(feet[bestFitElement].
z());
475 double best_pX = pMag_updated*sin_theta*cos_phi;
476 double best_pY = pMag_updated*sin_theta*sin_phi;
477 double best_pZ = pMag_updated*cos_theta;
478 CLHEP::Hep3Vector bestP(best_pX, best_pY, best_pZ);
486 lastUpdatedTSOS_Vect[iSet]= *(lastUpdatedTSOS_pointer[bestFitElement]);
495 const CLHEP::Hep3Vector& r3T,
499 bool detailedOutput){
501 double chi2 = 999999999999.;
503 double pMag_updated = 1./foot.x();
504 double sin_theta =
sin(foot.y());
505 double cos_theta =
cos(foot.y());
506 double sin_phi =
sin(foot.z());
507 double cos_phi =
cos(foot.z());
509 double pX = pMag_updated*sin_theta*cos_phi;
510 double pY = pMag_updated*sin_theta*sin_phi;
511 double pZ = pMag_updated*cos_theta;
513 muonCandidate, lastUpdatedTSOS,
514 trajectoryMeasurementsInTheSet, detailedOutput);
521 const CLHEP::Hep3Vector& r3T,
525 std::vector <CLHEP::Hep3Vector> morePoints;
526 double invP = key_foot.x();
527 double theta = key_foot.y();
528 double phi = key_foot.z();
531 double deltaPhi_init = 0.005;
532 double deltaTheta_init = 0.005;
534 double deltaInvP_init = 0.1 * invP;
538 bool optimized =
true;
545 bool detailedOutput =
false;
548 std::vector < double > pInv_init(3);
549 std::vector < double > theta_init(3);
550 std::vector < double > phi_init(3);
551 std::vector < double > chi2_init(3);
553 pInv_init[0] = invP - deltaInvP_init;
555 pInv_init[2] = invP + deltaInvP_init;
558 theta_init[0] = theta-deltaTheta_init;
559 theta_init[1] =
theta;
560 theta_init[2] = theta+deltaTheta_init;
562 phi_init[0] = phi-deltaPhi_init;
564 phi_init[2] = phi+deltaPhi_init;
566 double sin_theta_nom =
sin(theta_init[1]);
567 double cos_theta_nom =
cos(theta_init[1]);
568 double sin_phi_nom =
sin(phi_init[1]);
569 double cos_phi_nom =
cos(phi_init[1]);
570 double pMag_updated_nom = 1./pInv_init[1];
573 for(
int i=0;
i<3;++
i){
574 double pMag_updated = 1./pInv_init[
i];
575 double pX = pMag_updated*sin_theta_nom*cos_phi_nom;
576 double pY = pMag_updated*sin_theta_nom*sin_phi_nom;
577 double pZ = pMag_updated*cos_theta_nom;
579 muonCandidate, lastUpdatedTSOS,
580 trajectoryMeasurementsInTheSet, detailedOutput);
582 std::pair <double,double> result_pInv =
586 for(
int i=0;
i<3;++
i){
587 double sin_theta =
sin(theta_init[
i]);
588 double cos_theta =
cos(theta_init[i]);
589 double pX = pMag_updated_nom*sin_theta*cos_phi_nom;
590 double pY = pMag_updated_nom*sin_theta*sin_phi_nom;
591 double pZ = pMag_updated_nom*cos_theta;
593 muonCandidate, lastUpdatedTSOS,
594 trajectoryMeasurementsInTheSet, detailedOutput);
596 std::pair <double,double> result_theta =
600 for(
int i=0;
i<3;++
i){
601 double sin_phi =
sin(phi_init[
i]);
602 double cos_phi =
cos(phi_init[i]);
603 double pX = pMag_updated_nom*sin_theta_nom*cos_phi;
604 double pY = pMag_updated_nom*sin_theta_nom*sin_phi;
605 double pZ = pMag_updated_nom*cos_theta_nom;
607 muonCandidate, lastUpdatedTSOS,
608 trajectoryMeasurementsInTheSet, detailedOutput);
610 std::pair <double,double> result_phi =
613 double newPhi = result_phi.first;
614 if(fabs(newPhi - phi)<0.02){
617 CLHEP::Hep3Vector foot2(invP, theta, result_phi.first);
618 CLHEP::Hep3Vector foot3(invP, result_theta.first , phi);
619 double newInvP = result_pInv.first;
620 if(fabs(newInvP - invP)<0.001){
621 newInvP = invP + 0.001;
623 CLHEP::Hep3Vector foot4(result_pInv.first, theta, phi);
624 morePoints.push_back(foot2);
625 morePoints.push_back(foot3);
626 morePoints.push_back(foot4);
630 CLHEP::Hep3Vector foot2(invP, theta, phi-deltaPhi_init);
631 CLHEP::Hep3Vector foot3(invP, theta-deltaTheta_init, phi);
632 CLHEP::Hep3Vector foot4(invP-deltaInvP_init, theta, phi);
633 morePoints.push_back(foot2);
634 morePoints.push_back(foot3);
635 morePoints.push_back(foot4);
641 std::vector <double> &quadratic_chi2){
645 double paramAtMin = -99.;
646 std::vector <double> quadratic_param(3);
649 CLHEP::HepMatrix enumerator_1(3,3);
650 CLHEP::HepMatrix enumerator_2(3,3);
651 CLHEP::HepMatrix enumerator_3(3,3);
653 for(
int iCol=1;iCol<4;++iCol){
658 enumerator_1(1,iCol) = quadratic_chi2.at(iCol-1);
663 enumerator_2(2,iCol) = quadratic_chi2.at(iCol-1);
668 enumerator_3(3,iCol) = quadratic_chi2.at(iCol-1);
670 const double mult = 5.;
676 std::vector <CLHEP::HepMatrix> enumerator;
677 enumerator.push_back(enumerator_1);
678 enumerator.push_back(enumerator_2);
679 enumerator.push_back(enumerator_3);
681 double determinant = denominator.determinant();
682 if(fabs(determinant)>0.00000000001){
683 for(
int iPar=0;iPar<3;++iPar){
684 quadratic_param.at(iPar) = enumerator.at(iPar).determinant()/determinant;
690 if(quadratic_param.at(2)>0){
691 paramAtMin = - quadratic_param.at(1)/quadratic_param.at(2)/2;
695 paramAtMin = quadratic_var.at(1);
697 double chi2_min = quadratic_param.at(0) + quadratic_param.at(1)*paramAtMin + quadratic_param.at(2)*
pow(paramAtMin,2);
698 std::pair <double, double>
result;
699 result = std::make_pair ( paramAtMin, chi2_min);
704 unsigned int & high,
unsigned int & second_high,
unsigned int & low){
706 std::vector <double> chi2Feet_tmp = chi2Feet;
707 std::vector <double>::iterator minEl = min_element(chi2Feet.begin(), chi2Feet.end());
708 std::vector <double>::iterator maxEl = max_element(chi2Feet.begin(), chi2Feet.end());
709 high = maxEl - chi2Feet.begin();
710 low = minEl - chi2Feet.begin();
711 chi2Feet_tmp[high] = chi2Feet_tmp[low];
712 std::vector <double>::iterator second_maxEl = max_element(chi2Feet_tmp.begin(), chi2Feet_tmp.end());
713 second_high = second_maxEl - chi2Feet_tmp.begin();
719 unsigned int key_foot,
double scale ){
721 CLHEP::Hep3Vector newPosition(0.,0.,0.);
726 CLHEP::Hep3Vector centroid(0,0,0);
727 for(
unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
731 centroid += feet[iFoot];
733 centroid/=(feet.size()-1);
734 CLHEP::Hep3Vector displacement = 2.*(centroid - feet[key_foot]);
735 newPosition = feet[key_foot] + scale * displacement;
740 for(
unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
745 feet[iFoot] += feet[low];
void nDimContract(std::vector< CLHEP::Hep3Vector > &feet, unsigned int low)
std::string thePropagatorName
the propagator name
bool transform(Trajectory::DataContainer &measurements_segments, TransientTrackingRecHit::ConstRecHitContainer &hitContainer, TrajectoryStateOnSurface &firstTSOS)
transforms "segment trajectory" to "rechit container"
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
TrajectoryStateOnSurface theLastUpdatedTSOS
the trajectory state on the last available surface
double chi2AtSpecificStep(CLHEP::Hep3Vector &foot, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, TrajectoryStateOnSurface &lastUpdatedTSOS, Trajectory::DataContainer &trajectoryMeasurementsInTheSet, bool detailedOutput)
bool transformLight(Trajectory::DataContainer &measurements_segments, TransientTrackingRecHit::ConstRecHitContainer &hitContainer, TrajectoryStateOnSurface &firstTSOS)
transforms "segment trajectory" to "segment container"
const Propagator * propagator() const
access at the propagator
CLHEP::Hep3Vector momentum
bool operator()(TransientTrackingRecHit::ConstRecHitPointer hit_1, TransientTrackingRecHit::ConstRecHitPointer hit_2)
Sin< T >::type sin(const T &t)
TrajectoryStateOnSurface theLastButOneUpdatedTSOS
the trajectory state on the last but one available surface
Geom::Theta< T > theta() const
Trajectory::DataContainer trajectoryMeasurementsInTheSet
virtual SubDetector subDetector() const =0
The type of detector (PixelBarrel, PixelEndcap, TIB, TOB, TID, TEC, CSC, DT, RPCBarrel, RPCEndcap)
bool buildTrajectoryMeasurements(SeedCandidate *validSegmentsSet, Trajectory::DataContainer &finalCandidate)
from SeedCandidate to DataContainer only
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
double findMinChi2(unsigned int iSet, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, std::vector< TrajectoryStateOnSurface > &lastUpdatedTSOS_Vect, Trajectory::DataContainer &trajectoryMeasurementsInTheSet)
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
TrajectoryStateOnSurface lastUpdatedTSOS() const
the Trajectory state on the last surface of the fitting
std::vector< TrajectoryMeasurement > DataContainer
bool fwfit_SET(std::vector< SeedCandidate > &validSegmentsSet_in, std::vector< SeedCandidate > &validSegmentsSet_out)
Perform the SET inner-outward fitting.
FreeTrajectoryState * freeState(bool withErrors=true) const
virtual ~SETFilter()
Destructor.
Cos< T >::type cos(const T &t)
bool useSegmentsInTrajectory
SETFilter(const edm::ParameterSet &par, const MuonServiceProxy *service)
Constructor.
std::vector< CLHEP::Hep3Vector > find3MoreStartingPoints(CLHEP::Hep3Vector &key_foot, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
void pickElements(std::vector< double > &chi2Feet, unsigned int &high, unsigned int &second_high, unsigned int &low)
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
std::vector< ConstRecHitPointer > ConstRecHitContainer
const MuonServiceProxy * theService
double findChi2(double pX, double pY, double pZ, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, TrajectoryStateOnSurface &lastUpdatedTSOS, Trajectory::DataContainer &trajectoryMeasurementsInTheSet, bool detailedOutput)
void incrementChamberCounters(const DetLayer *layer)
Increment the DT,CSC,RPC counters.
const BoundPlane & surface() const
The nominal surface of the GeomDet.
MuonTransientTrackingRecHit::MuonRecHitContainer theSet
std::vector< const DetLayer * > theDetLayers
the det layer used in the reconstruction
Power< A, B >::type pow(const A &a, const B &b)
std::pair< double, double > findParabolaMinimum(std::vector< double > &quadratic_var, std::vector< double > &quadratic_chi2)
virtual void setEvent(const edm::Event &event)
Pass the Event to the algo at each event.
CLHEP::Hep3Vector reflectFoot(std::vector< CLHEP::Hep3Vector > &feet, unsigned int key_foot, double scale)