SETFilter Class Reference

#include <SETFilter.h>

Public Member Functions
bool	buildTrajectoryMeasurements (SeedCandidate *validSegmentsSet, Trajectory::DataContainer &finalCandidate)
	from SeedCandidate to DataContainer only More...
bool	fwfit_SET (std::vector< SeedCandidate > &validSegmentsSet_in, std::vector< SeedCandidate > &validSegmentsSet_out)
	Perform the SET inner-outward fitting. More...
int	getCSCChamberUsed () const
int	getDTChamberUsed () const
int	getRPCChamberUsed () const
int	getTotalChamberUsed () const
bool	goodState () const
const DetLayer *	lastDetLayer () const
	return the last det layer More...
FreeTrajectoryState	lastUpdatedFTS () const
	the last free trajectory state More...
TrajectoryStateOnSurface	lastUpdatedTSOS () const
	the Trajectory state on the last surface of the fitting More...
std::vector< const DetLayer * >	layers () const
	return the layer used for the refit More...
PropagationDirection	propagationDirection () const
	Return the propagation direction. More...
void	refit (const TrajectoryStateOnSurface &initialState, const DetLayer *, Trajectory &trajectory)
	Perform the inner-outward fitting. More...
void	reset ()
virtual void	setEvent (const edm::Event &event)
	Pass the Event to the algo at each event. More...
	SETFilter (const edm::ParameterSet &par, const MuonServiceProxy *service)
	Constructor. More...
bool	transform (Trajectory::DataContainer &measurements_segments, TransientTrackingRecHit::ConstRecHitContainer &hitContainer, TrajectoryStateOnSurface &firstTSOS)
	transforms "segment trajectory" to "rechit container" More...
bool	transformLight (Trajectory::DataContainer &measurements_segments, TransientTrackingRecHit::ConstRecHitContainer &hitContainer, TrajectoryStateOnSurface &firstTSOS)
	transforms "segment trajectory" to "segment container" More...
virtual	~SETFilter ()
	Destructor. More...

Private Member Functions
double	chi2AtSpecificStep (CLHEP::Hep3Vector &foot, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, TrajectoryStateOnSurface &lastUpdatedTSOS, Trajectory::DataContainer &trajectoryMeasurementsInTheSet, bool detailedOutput)
std::vector< CLHEP::Hep3Vector >	find3MoreStartingPoints (CLHEP::Hep3Vector &key_foot, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate)
double	findChi2 (double pX, double pY, double pZ, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, TrajectoryStateOnSurface &lastUpdatedTSOS, Trajectory::DataContainer &trajectoryMeasurementsInTheSet, bool detailedOutput)
double	findMinChi2 (unsigned int iSet, const CLHEP::Hep3Vector &r3T, SeedCandidate &muonCandidate, std::vector< TrajectoryStateOnSurface > &lastUpdatedTSOS_Vect, Trajectory::DataContainer &trajectoryMeasurementsInTheSet)
std::pair< double, double >	findParabolaMinimum (std::vector< double > &quadratic_var, std::vector< double > &quadratic_chi2)
void	incrementChamberCounters (const DetLayer *layer)
	Increment the DT,CSC,RPC counters. More...
void	nDimContract (std::vector< CLHEP::Hep3Vector > &feet, unsigned int low)
void	pickElements (std::vector< double > &chi2Feet, unsigned int &high, unsigned int &second_high, unsigned int &low)
const Propagator *	propagator () const
	access at the propagator More...
CLHEP::Hep3Vector	reflectFoot (std::vector< CLHEP::Hep3Vector > &feet, unsigned int key_foot, double scale)
void	setLastButOneUpdatedTSOS (TrajectoryStateOnSurface tsos)
	Set the last but one TSOS. More...
void	setLastUpdatedTSOS (TrajectoryStateOnSurface tsos)
	Set the last TSOS. More...

Private Attributes
int	cscChambers
int	dtChambers
int	rpcChambers
edm::ESHandle< TrajectoryFitter >	theBWLightFitter
	used in the SET BW fit More...
std::string	theBWLightFitterName
std::vector< const DetLayer * >	theDetLayers
	the det layer used in the reconstruction More...
NavigationDirection	theFitDirection
	the propagation direction More...
TrajectoryStateOnSurface	theLastButOneUpdatedTSOS
	the trajectory state on the last but one available surface More...
TrajectoryStateOnSurface	theLastUpdatedTSOS
	the trajectory state on the last available surface More...
std::string	thePropagatorName
	the propagator name More...
const MuonServiceProxy *	theService
int	totalChambers
bool	useSegmentsInTrajectory

Detailed Description

I. Bloch, E. James, S. Stoynev

Definition at line 49 of file SETFilter.h.

Constructor & Destructor Documentation

SETFilter::SETFilter	(	const edm::ParameterSet &	par,
		const MuonServiceProxy *	service
	)

Constructor.

Definition at line 52 of file SETFilter.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), thePropagatorName, and useSegmentsInTrajectory.

  :theService(service)//,
 //theOverlappingChambersFlag(true)
{
  thePropagatorName = par.getParameter<string>("Propagator");
  useSegmentsInTrajectory = par.getUntrackedParameter<bool>("UseSegmentsInTrajectory");
}

SETFilter::~SETFilter ( )

virtual

Destructor.

Definition at line 61 of file SETFilter.cc.

References LogTrace.

                     {

  LogTrace("Muon|RecoMuon|SETFilter")
    <<"SETFilter destructor called"<<endl;
  
}

Member Function Documentation

bool SETFilter::buildTrajectoryMeasurements	(	SeedCandidate *	validSegmentsSet,
		Trajectory::DataContainer &	finalCandidate
	)

from SeedCandidate to DataContainer only

Definition at line 135 of file SETFilter.cc.

References incrementChamberCounters(), reset(), theDetLayers, theLastUpdatedTSOS, and SeedCandidate::trajectoryMeasurementsInTheSet.

Referenced by SETMuonSeedProducer::produce().

                                                                                                              {
  // this is the SET algorithm fit
  bool validTrajectory = true;
  // reset the fitter 
  reset(); // the layer counters
  finalCandidate.clear();

  //---- Check if (only last?) TSOS is valid and build a trajectory (for the backward filter) 

  if(finalMuon->trajectoryMeasurementsInTheSet.size() &&
     finalMuon->trajectoryMeasurementsInTheSet.back().forwardPredictedState().isValid()){
    // loop over all measurements in the set
    for(unsigned int iMeas =0; iMeas<finalMuon->trajectoryMeasurementsInTheSet.size();++iMeas){
      // strore the measurements 
      finalCandidate.push_back(finalMuon->trajectoryMeasurementsInTheSet[iMeas]);
      const DetLayer *layer = finalMuon->trajectoryMeasurementsInTheSet[iMeas].layer();

      incrementChamberCounters(layer);

      theDetLayers.push_back(layer);

    }
    theLastUpdatedTSOS = finalMuon->trajectoryMeasurementsInTheSet.at(finalMuon->trajectoryMeasurementsInTheSet.size()-1).forwardPredictedState();
    //std::cout<<"  THE OUTPUT FROM FW FILTER: |P| = "<<finalMuon->momentum.mag()<<
    //" theta = "<<finalMuon->momentum.theta()<<" phi = "<<finalMuon->momentum.phi()<<std::endl;
  }
  else{
    validTrajectory = false;
    //std::cout<<" TSOS not valid; no trajectory build"<<std::endl;
  }
  return validTrajectory;
}

double SETFilter::chi2AtSpecificStep ( CLHEP::Hep3Vector &  foot, const CLHEP::Hep3Vector &  r3T, SeedCandidate &  muonCandidate, TrajectoryStateOnSurface &  lastUpdatedTSOS, Trajectory::DataContainer &  trajectoryMeasurementsInTheSet, bool  detailedOutput  )

private

Definition at line 497 of file SETFilter.cc.

References funct::cos(), findChi2(), and funct::sin().

Referenced by findMinChi2().

                                       {
  // specific input parameters - find chi2
  double chi2 = 999999999999.;
  if(foot.x()>0){ // this is |P|; maybe return a flag too?
    double pMag_updated = 1./foot.x();
    double sin_theta = sin(foot.y());
    double cos_theta = cos(foot.y());
    double sin_phi = sin(foot.z());
    double cos_phi = cos(foot.z());
      
    double pX = pMag_updated*sin_theta*cos_phi;
    double pY = pMag_updated*sin_theta*sin_phi;
    double pZ = pMag_updated*cos_theta;
    chi2 = findChi2(pX, pY, pZ, r3T, 
                    muonCandidate, lastUpdatedTSOS,
                    trajectoryMeasurementsInTheSet, detailedOutput);
  }
  return chi2;
}   

std::vector< CLHEP::Hep3Vector > SETFilter::find3MoreStartingPoints ( CLHEP::Hep3Vector &  key_foot, const CLHEP::Hep3Vector &  r3T, SeedCandidate &  muonCandidate  )

private

Definition at line 523 of file SETFilter.cc.

References funct::cos(), findChi2(), findParabolaMinimum(), i, lastUpdatedTSOS(), phi, funct::sin(), and theta().

Referenced by findMinChi2().

                                                 {
  // SIMPLEX uses nDim + 1 starting points; 
  // so here we need 3 more (one we already have)
  std::vector <CLHEP::Hep3Vector> morePoints;// again - CLHEP::Hep3Vector is not a good choice here
  double invP = key_foot.x();
  double theta = key_foot.y();
  double phi = key_foot.z();


  double deltaPhi_init = 0.005;
  double deltaTheta_init = 0.005;
  //double deltaInvP_init = 1.1 * invP;
  double deltaInvP_init = 0.1 * invP;
  //deltaInvP_init = 0.5 * invP;

  // try to chose better point to start with
  bool optimized = true;
  if(optimized){
    //---- Find a minimum chi2 for every variable (others are fixed) by supposing chi2 is a parabola
    //---- Then these points ("minima") are probably better starting points for the real minimization

    TrajectoryStateOnSurface  lastUpdatedTSOS;// fake here
    Trajectory::DataContainer trajectoryMeasurementsInTheSet;// fake here
    bool detailedOutput = false;// fake here


    std::vector < double > pInv_init(3);
    std::vector < double > theta_init(3);
    std::vector < double > phi_init(3);
    std::vector < double > chi2_init(3);

    pInv_init[0] = invP - deltaInvP_init;
    pInv_init[1] = invP;
    pInv_init[2] = invP + deltaInvP_init;
    //pInv_init[2] = invP +  0.1 * invP;

    theta_init[0] = theta-deltaTheta_init;
    theta_init[1] = theta;
    theta_init[2] = theta+deltaTheta_init;

    phi_init[0] = phi-deltaPhi_init;
    phi_init[1] = phi;
    phi_init[2] = phi+deltaPhi_init;

    double sin_theta_nom = sin(theta_init[1]);
    double cos_theta_nom = cos(theta_init[1]);
    double sin_phi_nom = sin(phi_init[1]);
    double cos_phi_nom = cos(phi_init[1]);
    double pMag_updated_nom =  1./pInv_init[1];

    //--- invP
    for(int i=0;i<3;++i){
      double pMag_updated = 1./pInv_init[i];
      double pX = pMag_updated*sin_theta_nom*cos_phi_nom;
      double pY = pMag_updated*sin_theta_nom*sin_phi_nom;
      double pZ = pMag_updated*cos_theta_nom;
      chi2_init[i] = findChi2(pX, pY, pZ, r3T,
                         muonCandidate, lastUpdatedTSOS,
                         trajectoryMeasurementsInTheSet, detailedOutput);
    }
    std::pair <double,double> result_pInv =
      findParabolaMinimum(pInv_init, chi2_init);

    //---- theta
    for(int i=0;i<3;++i){
      double sin_theta = sin(theta_init[i]);
      double cos_theta = cos(theta_init[i]);
      double pX = pMag_updated_nom*sin_theta*cos_phi_nom;
      double pY = pMag_updated_nom*sin_theta*sin_phi_nom;
      double pZ = pMag_updated_nom*cos_theta;
        chi2_init[i] =  findChi2(pX, pY, pZ, r3T,
                         muonCandidate, lastUpdatedTSOS,
                         trajectoryMeasurementsInTheSet, detailedOutput);
    }
    std::pair <double,double> result_theta =
      findParabolaMinimum(theta_init, chi2_init);

    //---- phi
    for(int i=0;i<3;++i){
      double sin_phi = sin(phi_init[i]);
      double cos_phi = cos(phi_init[i]);
      double pX = pMag_updated_nom*sin_theta_nom*cos_phi;
      double pY = pMag_updated_nom*sin_theta_nom*sin_phi;
      double pZ = pMag_updated_nom*cos_theta_nom;
      chi2_init[i] =  findChi2(pX, pY, pZ, r3T,
                         muonCandidate, lastUpdatedTSOS,
                         trajectoryMeasurementsInTheSet, detailedOutput);
    }
    std::pair <double,double> result_phi =
      findParabolaMinimum(phi_init, chi2_init);
    // should we use that?
    double newPhi = result_phi.first;
    if(fabs(newPhi - phi)<0.02){// too close?
      newPhi = phi + 0.02;
    }
    CLHEP::Hep3Vector foot2(invP, theta, result_phi.first);
    CLHEP::Hep3Vector foot3(invP, result_theta.first , phi);
    double newInvP = result_pInv.first;
    if(fabs(newInvP - invP)<0.001){//too close
      newInvP = invP + 0.001;
    }
    CLHEP::Hep3Vector foot4(result_pInv.first, theta, phi);
    morePoints.push_back(foot2);
    morePoints.push_back(foot3);
    morePoints.push_back(foot4);
  }
  else{
    // the points
    CLHEP::Hep3Vector foot2(invP, theta, phi-deltaPhi_init);
    CLHEP::Hep3Vector foot3(invP, theta-deltaTheta_init, phi);
    CLHEP::Hep3Vector foot4(invP-deltaInvP_init, theta, phi);
    morePoints.push_back(foot2);
    morePoints.push_back(foot3);
    morePoints.push_back(foot4);
  }
  return morePoints;
}

double SETFilter::findChi2 ( double  pX, double  pY, double  pZ, const CLHEP::Hep3Vector &  r3T, SeedCandidate &  muonCandidate, TrajectoryStateOnSurface &  lastUpdatedTSOS, Trajectory::DataContainer &  trajectoryMeasurementsInTheSet, bool  detailedOutput  )

private

Definition at line 242 of file SETFilter.cc.

References RecoTauCleanerPlugins::charge, SeedCandidate::charge, TrackingRecHit::dimension(), alignCSCRings::e, TrajectoryStateOnSurface::freeState(), TrackingRecHit::geographicalId(), GenericTransientTrackingRecHit::hit(), TrajectoryStateOnSurface::isValid(), GenericTransientTrackingRecHit::localPosition(), GenericTransientTrackingRecHit::localPositionError(), funct::pow(), Propagator::propagate(), propagator(), GeomDet::surface(), theService, SeedCandidate::theSet, GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), and LocalError::yy().

Referenced by chi2AtSpecificStep(), and find3MoreStartingPoints().

                                   {
  //---- actual chi2 calculations; only the measurement error is taken into accout!
  //---- chi2 is to compare an extrapolated point to various measurements so
  //---- the extrapolation error is not an issue (is it?)

  if(detailedOutput){
    trajectoryMeasurementsInTheSet.clear();
  }

  int charge =  muonCandidate.charge;
  GlobalVector p3GV(pX,pY,pZ);
  GlobalPoint r3GP(r3T.x(), r3T.y(), r3T.z());
  //---- how to disable error propagation?
  // VI: just not set it!
  FreeTrajectoryState ftsStart(r3GP, p3GV, charge, &*(theService->magneticField()));
  // VI let's be backward compatible...
  if(detailedOutput) {
    AlgebraicSymMatrix55 cov; cov*=1e-20;
    ftsStart.setCurvilinearError(cov);
  }
  TrajectoryStateOnSurface tSOSDest;
    
  double chi2_loc = 0.;
  for(unsigned int iMeas = 0; iMeas <muonCandidate.theSet.size(); ++iMeas){
    MuonTransientTrackingRecHit::MuonRecHitPointer muonRecHit =  muonCandidate.theSet[iMeas];
    DetId detId = muonRecHit->hit()->geographicalId();
    const GeomDet* layer = theService->trackingGeometry()->idToDet(detId);

    //---- propagate the muon starting from position r3T and momentum p3T 

    //    bool radX0CorrectionMode_ = false; // not needed here?
    //if (radX0CorrectionMode_ ){
    //} else {

    tSOSDest = propagator()->propagate(ftsStart, layer->surface());
    lastUpdatedTSOS = tSOSDest;
    if (tSOSDest.isValid()){
      //---- start next step ("adding" measurement) from the last TSOS
      ftsStart = *tSOSDest.freeState();
    } else{
      //std::cout<<"... not valid TSOS"<<std::endl;
      chi2_loc = 9999999999.;
      break;
    }
    //}

    LocalPoint locHitPos = muonRecHit->localPosition();
    LocalError locHitErr = muonRecHit->localPositionError();
    const GlobalPoint globPropPos = ftsStart.position();
    LocalPoint locPropPos = layer->toLocal(globPropPos);

    //
    //---- chi2 calculated in local system; correlation taken into accont
    CLHEP::HepMatrix dist(1,2);//, distT(2,1);
    double  chi2_intermed = -9;
    int ierr = 0;
    dist(1,1) = locPropPos.x() - locHitPos.x();
    dist(1,2) = locPropPos.y() - locHitPos.y();
    CLHEP::HepMatrix IC(2,2);
    IC(1,1) = locHitErr.xx();
    IC(2,1) = locHitErr.xy();
    IC(2,2) = locHitErr.yy();
    IC(1,2) = IC(2,1);

    //---- Special care is needed for "1-dim measurements" (RPCs, outer DT(?))
    if(4!=muonRecHit->hit()->dimension()){
      for(int iE = 1;iE<3;++iE){
    // this is bellow is a DT fix; hopefully it will not be needed
        if ( fabs(IC(iE,iE)) < 0.00000001){ 
          IC(iE,iE) = 62500./12.;// error squared - ( 250 cm /sqrt(12) )^2; large 
        }
      }
    }
    //---- Invert covariance matrix
    IC.invert(ierr);
    //if (ierr != 0) {
    //std::cout << "failed to invert covariance matrix (2x2) =\n" << IC << std::endl;;
    //}
    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);
    if(chi2_intermed<0){// should we check?
       chi2_intermed = 9999999999.;
    }
    chi2_loc += chi2_intermed;

    // that's for the last call; we don't need to construct a TrajectoryMeasurement at every chi2 step
    if(detailedOutput){
      DetId detId = muonRecHit->hit()->geographicalId();
      const DetLayer *layer = theService->detLayerGeometry()->idToLayer( detId);
      //std::cout<<"    seg pos in traj : "<<lastUpdatedTSOS.globalPosition()<<std::endl;
      // put the measurement into the set
      // VI set the error as the fit needs it... (it is nonsense anyhow...)   
      // (do it on the tsos)
      /*
      if (!lastUpdatedTSOS.hasError()){
    AlgebraicSymMatrix55 cov; cov*=1e6;
    lastUpdatedTSOS.freeTrajectoryState().setCurvilinearError(cov);
      }
      */
      trajectoryMeasurementsInTheSet.push_back( TrajectoryMeasurement
                        ( lastUpdatedTSOS,
                          muonRecHit,
                          chi2_intermed,
                          layer ) );
    }
  }
  return chi2_loc;
}

double SETFilter::findMinChi2 ( unsigned int  iSet, const CLHEP::Hep3Vector &  r3T, SeedCandidate &  muonCandidate, std::vector< TrajectoryStateOnSurface > &  lastUpdatedTSOS_Vect, Trajectory::DataContainer &  trajectoryMeasurementsInTheSet  )

private

Definition at line 355 of file SETFilter.cc.

References chi2AtSpecificStep(), funct::cos(), find3MoreStartingPoints(), lastUpdatedTSOS(), SeedCandidate::momentum, nDimContract(), phi, pickElements(), reflectFoot(), funct::sin(), theta(), SeedCandidate::trajectoryMeasurementsInTheSet, y, and z.

Referenced by fwfit_SET().

                                                                                                  {
  // a chi2 minimization procedure 

  //---- Which three variables to use? 
  //---- (1/|P|, theta, phi) ? in that case many sin() and cos() operations :-/
  //---- maybe vary directly sin() and cos()?
  bool detailedOutput = false;

  double cosTheta = muonCandidate.momentum.cosTheta();
  double theta = acos(cosTheta);
  double phi = muonCandidate.momentum.phi();
  double pMag = muonCandidate.momentum.mag();

  double minChi2 = -999.;
  TrajectoryStateOnSurface  lastUpdatedTSOS;

  //---- Fit Parameters

  if(pMag<5.){// hardcoded - remove it! same in SETSeedFinder
    pMag = 5.;// GeV
  }
  //---- This offset helps the minimization to go faster (in the specific case)
  pMag *=1.2;
  double invP = 1./pMag;
  //std::cout<<"    INIT pMag = "<<pMag<<" invP = "<<invP<<" theta = "<<theta<<" phi = "<<phi<<std::endl;

  //---- apply downhill SIMPLEX minimization (also "amoeba" method; thus the "feet" below are  amoeba's feet)

  //std::cout<<"    SIMPLEX minimization"<<std::endl;
  //---- parameters ; the should be hardcoded   
  const double reflect = 1;
  const double expand = -0.5;
  const double contract = 0.25;

  const int nDim = 3; // invP, theta, phi
  //---- Now choose nDim + 1 points
  std::vector <CLHEP::Hep3Vector> feet(nDim+1);
  std::vector <double> chi2Feet(nDim+1);
  std::vector <TrajectoryStateOnSurface*> lastUpdatedTSOS_pointer(nDim+1);// probably not needed; to be deleted
   
  //---- The minimization procedure strats from these nDim+1 points (feet)

  CLHEP::Hep3Vector foot1(invP, theta, phi);// well obviuosly it is not a real Hep3Vector; better change it to a simple vector
  feet[0] = foot1;
  chi2Feet[0] = chi2AtSpecificStep(feet[0], r3T, muonCandidate, lastUpdatedTSOS,
                                       trajectoryMeasurementsInTheSet, detailedOutput);
  lastUpdatedTSOS_pointer[0] = &lastUpdatedTSOS;

  std::vector <CLHEP::Hep3Vector> morePoints =
    find3MoreStartingPoints(feet[0], r3T, muonCandidate);
  feet[1] = morePoints[0];
  feet[2] = morePoints[1];
  feet[3] = morePoints[2];

  //--- SIMPLEX initial step(s)
  for(unsigned int iFoot = 1; iFoot<feet.size();++iFoot){

    chi2Feet[iFoot] = chi2AtSpecificStep(feet[iFoot], r3T, muonCandidate, lastUpdatedTSOS,
                                         trajectoryMeasurementsInTheSet, detailedOutput);
    lastUpdatedTSOS_pointer[iFoot] = &lastUpdatedTSOS;
  }

  unsigned int high, second_high, low;  
  //const unsigned int iterations = 1;//---- to be replaced by something better
  int iCalls = 0;
  //for(unsigned int iIt = 0;iIt<iterations;++iIt){
  while(iCalls<3.){
    //---- SIMPLEX step 1
    pickElements(chi2Feet, high, second_high, low);
    ++iCalls;
    feet[high] = reflectFoot(feet, high, reflect );
    chi2Feet[high] = chi2AtSpecificStep(feet[high], r3T, muonCandidate, lastUpdatedTSOS,
                                        trajectoryMeasurementsInTheSet, detailedOutput);
    lastUpdatedTSOS_pointer[high] = &lastUpdatedTSOS;
    //---- SIMPLEX step 2.1
    if(chi2Feet[high] <chi2Feet[low]){
      ++iCalls; 
      feet[high] = reflectFoot(feet, high, expand);
      chi2Feet[high] = chi2AtSpecificStep(feet[high], r3T, muonCandidate, lastUpdatedTSOS,
                                          trajectoryMeasurementsInTheSet, detailedOutput);
      lastUpdatedTSOS_pointer[high] = &lastUpdatedTSOS;
    }
    //---- SIMPLEX step 2.2
    else if( chi2Feet[high] > chi2Feet[second_high]){
      double worstChi2 = chi2Feet[high];
      ++iCalls;
      feet[high] =  reflectFoot(feet, high, contract);
      chi2Feet[high] = chi2AtSpecificStep(feet[high], r3T, muonCandidate, lastUpdatedTSOS,
                                          trajectoryMeasurementsInTheSet, detailedOutput);
      lastUpdatedTSOS_pointer[high] = &lastUpdatedTSOS;
      //---- SIMPLEX step 2.2.1
      if(chi2Feet[high] <worstChi2){
        nDimContract(feet, low);
        for(unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
      ++iCalls; 
          chi2Feet[iFoot] = chi2AtSpecificStep(feet[iFoot], r3T, muonCandidate, lastUpdatedTSOS,
                                               trajectoryMeasurementsInTheSet, detailedOutput);
          lastUpdatedTSOS_pointer[iFoot] = &lastUpdatedTSOS;
        }
    --iCalls;// one of the above is not changed
      }
    }
  }
  //---- Here the SIMPLEX minimization ends

  // this is the minimum found
  int bestFitElement = min_element(chi2Feet.begin(),chi2Feet.end()) - chi2Feet.begin();

  //---- repeat to get the trajectoryMeasurementsInTheSet (optimize?)
  detailedOutput = true;
  chi2Feet[bestFitElement] = chi2AtSpecificStep(feet[bestFitElement], r3T, muonCandidate, lastUpdatedTSOS,
                                                trajectoryMeasurementsInTheSet, detailedOutput);
  minChi2 = chi2Feet[bestFitElement];

  double pMag_updated = 1./feet[bestFitElement].x();
  double sin_theta = sin(feet[bestFitElement].y());
  double cos_theta = cos(feet[bestFitElement].y());
  double sin_phi = sin(feet[bestFitElement].z());
  double cos_phi = cos(feet[bestFitElement].z());

  double best_pX = pMag_updated*sin_theta*cos_phi;
  double best_pY = pMag_updated*sin_theta*sin_phi;
  double best_pZ = pMag_updated*cos_theta;
  CLHEP::Hep3Vector bestP(best_pX, best_pY, best_pZ);
  //---- update the best momentum estimate  
  //if(minChi2<999999. && pMag_updated>0.){//fit failed?! check
  muonCandidate.momentum = bestP;
  //}
  //---- update the trajectory
  muonCandidate.trajectoryMeasurementsInTheSet = trajectoryMeasurementsInTheSet;
  // do we need that?
  lastUpdatedTSOS_Vect[iSet]= *(lastUpdatedTSOS_pointer[bestFitElement]);

  //std::cout<<"   FINAL:  P = "<<muonCandidate.momentum.mag()<<" theta = "<<muonCandidate.momentum.theta()<<
  //" phi = "<<muonCandidate.momentum.phi()<<"   chi = "<<chi2Feet[bestFitElement]<<std::endl;
  return minChi2;
}

std::pair< double, double > SETFilter::findParabolaMinimum ( std::vector< double > &  quadratic_var, std::vector< double > &  quadratic_chi2  )

private

Definition at line 643 of file SETFilter.cc.

References ztail::d, cuy::denominator, VarParsing::mult, funct::pow(), and query::result.

Referenced by find3MoreStartingPoints().

                                                                                                       {

  // quadratic equation minimization

  double paramAtMin = -99.;
  std::vector <double> quadratic_param(3);

  math::Matrix<3,3>::type denominator;
  math::Matrix<3,3>::type enumerator_1;
  math::Matrix<3,3>::type enumerator_2;
  math::Matrix<3,3>::type enumerator_3;

  for(int iCol=0;iCol<3;++iCol){
    denominator(0,iCol) = 1;
    denominator(1,iCol) = quadratic_var.at(iCol);
    denominator(2,iCol) = pow(quadratic_var.at(iCol),2);

    enumerator_1(0,iCol) = quadratic_chi2.at(iCol);
    enumerator_1(1,iCol) = denominator(1,iCol);
    enumerator_1(2,iCol) = denominator(2,iCol);

    enumerator_2(0,iCol) = denominator(0,iCol);
    enumerator_2(1,iCol) = quadratic_chi2.at(iCol);
    enumerator_2(2,iCol) = denominator(2,iCol);

    enumerator_3(0,iCol) = denominator(0,iCol);
    enumerator_3(1,iCol) = denominator(1,iCol);
    enumerator_3(2,iCol) = quadratic_chi2.at(iCol);
  }
  const double mult = 5.;// "save" accuracy"; result is independent on "mult"
  denominator *=mult;
  enumerator_1*=mult;
  enumerator_2*=mult;
  enumerator_3*=mult;

  std::vector <math::Matrix<3,3>::type> enumerator;
  enumerator.push_back(enumerator_1);
  enumerator.push_back(enumerator_2);
  enumerator.push_back(enumerator_3);

  double determinant=0;
  denominator.Det2(determinant);
  if(fabs(determinant)>0.00000000001){
    for(int iPar=0;iPar<3;++iPar){
      double d = 0.;
      enumerator.at(iPar).Det2(d);
      quadratic_param.at(iPar) = d/determinant;
    }
  }
  else{
    //std::cout<<" determinant 0? Check the code..."<<std::endl;
  }
  if(quadratic_param.at(2)>0){
    paramAtMin = - quadratic_param.at(1)/quadratic_param.at(2)/2;
  }
  else {
    //std::cout<<" parabola has a maximum or division by zero... Using initial value. "<<std::endl;
    paramAtMin = quadratic_var.at(1);
  }
  double chi2_min = quadratic_param.at(0) + quadratic_param.at(1)*paramAtMin + quadratic_param.at(2)*pow(paramAtMin,2);
  std::pair <double, double> result;
  result =  std::make_pair ( paramAtMin, chi2_min);
  return result;
}

bool SETFilter::fwfit_SET	(	std::vector< SeedCandidate > &	validSegmentsSet_in,
		std::vector< SeedCandidate > &	validSegmentsSet_out
	)

Perform the SET inner-outward fitting.

Definition at line 100 of file SETFilter.cc.

References findMinChi2().

Referenced by SETMuonSeedProducer::produce().

                                                                     {
  // this is the SET algorithm fit
  validSegmentsSet_out.clear();

  //---- It is supposed to be called within a loop over "segment clusters"; 
  //---- so std::vector < SeedCandidate> consists of "valid" combinations (sets) within a "cluster"
  //---- "seed" above has nothing to do with the "Seed" in the STA code
 
  // a trajectory is not really build but a TSOS is build and is checked (below)
  bool validStep = true;
  std::vector <double> chi2AllCombinations(validSegmentsSet_in.size());
  std::vector < TrajectoryStateOnSurface > lastUpdatedTSOS_Vect(validSegmentsSet_in.size());
  // loop over all valid sets
  for(unsigned int iSet = 0; iSet<validSegmentsSet_in.size(); ++iSet){
    //std::cout<<"   iSET = "<<iSet<<std::endl;
    //---- start fit from the origin
    CLHEP::Hep3Vector origin (0.,0.,0.);
    Trajectory::DataContainer trajectoryMeasurementsInTheSet_tmp;
    //---- Find minimum chi2 (corresponding to a specific 3D-momentum)    
    chi2AllCombinations[iSet]  = findMinChi2(iSet, origin, validSegmentsSet_in[iSet], lastUpdatedTSOS_Vect,
                                             trajectoryMeasurementsInTheSet_tmp);
  }
  //---- Find the best muon candidate (min chi2) in the cluster; find more candidates?
  std::vector < double >::iterator itMin = min_element(chi2AllCombinations.begin(),chi2AllCombinations.end());

  int positionMin = itMin - chi2AllCombinations.begin();

  // "the best" set; have to find reasonable conditions to include more than one set   
  validSegmentsSet_out.push_back(validSegmentsSet_in[positionMin]);

  return validStep;
}

int SETFilter::getCSCChamberUsed ( ) const

inline

Definition at line 97 of file SETFilter.h.

References cscChambers.

{return cscChambers;}

int SETFilter::getDTChamberUsed ( ) const

inline

Definition at line 96 of file SETFilter.h.

References dtChambers.

{return dtChambers;}

int SETFilter::getRPCChamberUsed ( ) const

inline

Definition at line 98 of file SETFilter.h.

References rpcChambers.

{return rpcChambers;}

int SETFilter::getTotalChamberUsed ( ) const

inline

Definition at line 95 of file SETFilter.h.

References totalChambers.

{return totalChambers;}

bool SETFilter::goodState ( ) const

inline

Definition at line 100 of file SETFilter.h.

References cscChambers, dtChambers, and totalChambers.

                                {return totalChambers >= 2 && 
                   ((dtChambers + cscChambers) >0 );}

void SETFilter::incrementChamberCounters ( const DetLayer *  layer )

private

Increment the DT,CSC,RPC counters.

Definition at line 85 of file SETFilter.cc.

References GeomDetEnumerators::CSC, cscChambers, GeomDetEnumerators::DT, dtChambers, GeomDetEnumerators::RPCBarrel, rpcChambers, GeomDetEnumerators::RPCEndcap, DetLayer::subDetector(), and totalChambers.

Referenced by buildTrajectoryMeasurements().

                                                             {

  if(layer->subDetector()==GeomDetEnumerators::DT) dtChambers++; 
  else if(layer->subDetector()==GeomDetEnumerators::CSC) cscChambers++; 
  else if(layer->subDetector()==GeomDetEnumerators::RPCBarrel || layer->subDetector()==GeomDetEnumerators::RPCEndcap) rpcChambers++; 
  else 
    LogError("Muon|RecoMuon|SETFilter")
      << "Unrecognized module type in incrementChamberCounters";
  // FIXME:
  //   << layer->module() << " " <<layer->Part() << endl;
  
  totalChambers++;
}

const DetLayer* SETFilter::lastDetLayer ( ) const

inline

return the last det layer

Definition at line 107 of file SETFilter.h.

References theDetLayers.

{return theDetLayers.back();}

FreeTrajectoryState SETFilter::lastUpdatedFTS ( ) const

inline

the last free trajectory state

Definition at line 64 of file SETFilter.h.

References TrajectoryStateOnSurface::freeTrajectoryState(), and theLastUpdatedTSOS.

{return *theLastUpdatedTSOS.freeTrajectoryState();}

TrajectoryStateOnSurface SETFilter::lastUpdatedTSOS ( ) const

inline

the Trajectory state on the last surface of the fitting

Definition at line 67 of file SETFilter.h.

References theLastUpdatedTSOS.

Referenced by find3MoreStartingPoints(), and findMinChi2().

{return theLastUpdatedTSOS;}

std::vector<const DetLayer*> SETFilter::layers ( ) const

inline

return the layer used for the refit

Definition at line 104 of file SETFilter.h.

References theDetLayers.

Referenced by SETMuonSeedProducer::produce().

{return theDetLayers;}

void SETFilter::nDimContract ( std::vector< CLHEP::Hep3Vector > &  feet, unsigned int  low  )

private

Definition at line 745 of file SETFilter.cc.

Referenced by findMinChi2().

                                                                                {
  for(unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
    // a SIMPLEX function
    if(iFoot==low){
      continue;
    }
    feet[iFoot] +=  feet[low];
    feet[iFoot] /=2.;
  }
  return;
}

void SETFilter::pickElements ( std::vector< double > &  chi2Feet, unsigned int &  high, unsigned int &  second_high, unsigned int &  low  )

private

Definition at line 709 of file SETFilter.cc.

Referenced by findMinChi2().

                                                                                                           {
  // a SIMPLEX function
  std::vector <double> chi2Feet_tmp = chi2Feet;
  std::vector <double>::iterator minEl = min_element(chi2Feet.begin(), chi2Feet.end());
  std::vector <double>::iterator maxEl = max_element(chi2Feet.begin(), chi2Feet.end());
  high = maxEl - chi2Feet.begin();
  low = minEl - chi2Feet.begin();
  chi2Feet_tmp[high] = chi2Feet_tmp[low]; 
  std::vector <double>::iterator second_maxEl = max_element(chi2Feet_tmp.begin(), chi2Feet_tmp.end());
  second_high = second_maxEl - chi2Feet_tmp.begin();

  return;
}

PropagationDirection SETFilter::propagationDirection

(

)

const

Return the propagation direction.

const Propagator * SETFilter::propagator ( ) const

private

access at the propagator

Definition at line 80 of file SETFilter.cc.

References thePropagatorName, and theService.

Referenced by findChi2(), and transform().

                                              {
  return &*theService->propagator(thePropagatorName);
}

void SETFilter::refit	(	const TrajectoryStateOnSurface &	initialState,
		const DetLayer *	,
		Trajectory &	trajectory
	)

Perform the inner-outward fitting.

CLHEP::Hep3Vector SETFilter::reflectFoot ( std::vector< CLHEP::Hep3Vector > &  feet, unsigned int  key_foot, double  scale  )

private

Definition at line 724 of file SETFilter.cc.

Referenced by findMinChi2().

                                                                                      {
  // a SIMPLEX function
  CLHEP::Hep3Vector newPosition(0.,0.,0.);
  if(scale==0.5){
    //std::cout<<" STA muon: scale parameter for simplex method incorrect : "<<scale<<std::endl;
    return newPosition;
  }
  CLHEP::Hep3Vector centroid(0,0,0);
  for(unsigned int iFoot = 0; iFoot<feet.size();++iFoot){
    if(iFoot==key_foot){
      continue;
    }
    centroid += feet[iFoot];
  }
  centroid/=(feet.size()-1);
  CLHEP::Hep3Vector displacement = 2.*(centroid - feet[key_foot]);
  newPosition = feet[key_foot] + scale * displacement;
  return newPosition;
}

void SETFilter::reset

(

void

)

Definition at line 71 of file SETFilter.cc.

References cscChambers, dtChambers, rpcChambers, theDetLayers, theLastButOneUpdatedTSOS, theLastUpdatedTSOS, and totalChambers.

Referenced by buildTrajectoryMeasurements(), and SETMuonSeedProducer::produce().

                     {
  totalChambers = dtChambers = cscChambers = rpcChambers = 0;
  
  theLastUpdatedTSOS =  theLastButOneUpdatedTSOS = TrajectoryStateOnSurface();

  theDetLayers.clear();
}

void SETFilter::setEvent ( const edm::Event &  event )

virtual

Pass the Event to the algo at each event.

Definition at line 68 of file SETFilter.cc.

Referenced by SETMuonSeedProducer::setEvent().

                                          {
}

void SETFilter::setLastButOneUpdatedTSOS ( TrajectoryStateOnSurface  tsos )

inlineprivate

Set the last but one TSOS.

Definition at line 120 of file SETFilter.h.

References theLastButOneUpdatedTSOS.

{ theLastButOneUpdatedTSOS = tsos;}

void SETFilter::setLastUpdatedTSOS ( TrajectoryStateOnSurface  tsos )

inlineprivate

Set the last TSOS.

Definition at line 117 of file SETFilter.h.

References theLastUpdatedTSOS.

{ theLastUpdatedTSOS = tsos;}

bool SETFilter::transform	(	Trajectory::DataContainer &	measurements_segments,
		TransientTrackingRecHit::ConstRecHitContainer &	hitContainer,
		TrajectoryStateOnSurface &	firstTSOS
	)

transforms "segment trajectory" to "rechit container"

Definition at line 169 of file SETFilter.cc.

References TrajectoryStateOnSurface::isValid(), Propagator::propagate(), propagator(), python.multivaluedict::sort(), sortRadius, summarizeEdmComparisonLogfiles::success, GeomDet::surface(), and theService.

Referenced by SETMuonSeedProducer::produce().

                                                   {
  // transforms "segment trajectory" to "rechit container"
  //sort(measurements_segments.begin(),measurements_segments.end(),sortRadius);
  bool success = true;
  // loop over all segments in the trajectory
  for(int iMeas = measurements_segments.size() - 1; iMeas>-1;--iMeas){
    TransientTrackingRecHit ::ConstRecHitContainer sortedHits;
    // loop over the rechits contained in the segments
    for(unsigned int jMeas = 0; jMeas<measurements_segments[iMeas].recHit()->transientHits().size();++jMeas){
      if(measurements_segments[iMeas].recHit()->transientHits().at(jMeas)->transientHits().size()>1){
    // loop over the rechits contained in the rechits contained in the segments (OK, OK - this is for DT only;
    // the convention there is a bit different from the CSCs)
    for(unsigned int kMeas = 0;kMeas<measurements_segments[iMeas].recHit()->transientHits().at(jMeas)->transientHits().size();
        ++kMeas){
      sortedHits.push_back( 
                   measurements_segments[iMeas].recHit()->transientHits().at(jMeas)->transientHits().at(kMeas));
    }
      }
      else{
        sortedHits = measurements_segments[iMeas].recHit()->transientHits();
      }
    }
    // sort the rechits by radius (or z) and put them in a container
    sort(sortedHits.begin(),sortedHits.end(),sortRadius);
    hitContainer.insert(hitContainer.end(),sortedHits.begin(),sortedHits.end());    
  }

  // this is the last segment state
  FreeTrajectoryState ftsStart = *(measurements_segments.at(measurements_segments.size()-1).forwardPredictedState().freeState());

  // this is the last (from the IP) rechit
  TransientTrackingRecHit::ConstRecHitPointer muonRecHit =  hitContainer[0];
  DetId detId_last = hitContainer[0]->hit()->geographicalId();
  const GeomDet* layer_last = theService->trackingGeometry()->idToDet(detId_last);

  // get the last rechit TSOS
  TrajectoryStateOnSurface tSOSDest = propagator()->propagate(ftsStart, layer_last->surface());
  firstTSOS = tSOSDest;
  // ftsStart should be at the last rechit surface
  if (!tSOSDest.isValid()){
    success = false;
    //     ftsStart = *tSOSDest.freeState();
  }
  return success; 
}

bool SETFilter::transformLight	(	Trajectory::DataContainer &	measurements_segments,
		TransientTrackingRecHit::ConstRecHitContainer &	hitContainer,
		TrajectoryStateOnSurface &	firstTSOS
	)

transforms "segment trajectory" to "segment container"

Definition at line 217 of file SETFilter.cc.

References summarizeEdmComparisonLogfiles::success, and useSegmentsInTrajectory.

Referenced by SETMuonSeedProducer::produce().

                                                        {
  // transforms "segment trajectory" to "segment container"

  bool success = true;
  // loop over all segments in the trajectory
  if(useSegmentsInTrajectory){// if segments the "backword fit" (rechits)
                              // performed later is actually a forward one (?!) 
    for(unsigned int iMeas = 0; iMeas<measurements_segments.size();++iMeas){
      hitContainer.push_back(measurements_segments[iMeas].recHit());
    }
  }
  else{
    for(int iMeas = measurements_segments.size() - 1; iMeas>-1;--iMeas){
      hitContainer.push_back(measurements_segments[iMeas].recHit());
    }

  }
  // this is the last segment state
  firstTSOS = measurements_segments.at(0).forwardPredictedState();
  return success;
}

Member Data Documentation

int SETFilter::cscChambers

private

Definition at line 183 of file SETFilter.h.

Referenced by getCSCChamberUsed(), goodState(), incrementChamberCounters(), and reset().

int SETFilter::dtChambers

private

Definition at line 182 of file SETFilter.h.

Referenced by getDTChamberUsed(), goodState(), incrementChamberCounters(), and reset().

int SETFilter::rpcChambers

private

Definition at line 184 of file SETFilter.h.

Referenced by getRPCChamberUsed(), incrementChamberCounters(), and reset().

edm::ESHandle<TrajectoryFitter> SETFilter::theBWLightFitter

private

used in the SET BW fit

Definition at line 189 of file SETFilter.h.

std::string SETFilter::theBWLightFitterName

private

Definition at line 190 of file SETFilter.h.

std::vector<const DetLayer*> SETFilter::theDetLayers

private

the det layer used in the reconstruction

Definition at line 179 of file SETFilter.h.

Referenced by buildTrajectoryMeasurements(), lastDetLayer(), layers(), and reset().

NavigationDirection SETFilter::theFitDirection

private

the propagation direction

Definition at line 171 of file SETFilter.h.

TrajectoryStateOnSurface SETFilter::theLastButOneUpdatedTSOS

private

the trajectory state on the last but one available surface

Definition at line 176 of file SETFilter.h.

Referenced by reset(), and setLastButOneUpdatedTSOS().

TrajectoryStateOnSurface SETFilter::theLastUpdatedTSOS

private

the trajectory state on the last available surface

Definition at line 174 of file SETFilter.h.

Referenced by buildTrajectoryMeasurements(), lastUpdatedFTS(), lastUpdatedTSOS(), reset(), and setLastUpdatedTSOS().

std::string SETFilter::thePropagatorName

private

the propagator name

Definition at line 168 of file SETFilter.h.

Referenced by propagator(), and SETFilter().

const MuonServiceProxy* SETFilter::theService

private

Definition at line 192 of file SETFilter.h.

Referenced by findChi2(), propagator(), and transform().

int SETFilter::totalChambers

private

Definition at line 181 of file SETFilter.h.

Referenced by getTotalChamberUsed(), goodState(), incrementChamberCounters(), and reset().

bool SETFilter::useSegmentsInTrajectory

private

Definition at line 186 of file SETFilter.h.

Referenced by SETFilter(), and transformLight().