SETSeedFinder Class Reference

#include <SETSeedFinder.h>

Inheritance diagram for SETSeedFinder:

Public Types
typedef MuonTransientTrackingRecHit::MuonRecHitContainer	MuonRecHitContainer

Public Member Functions
std::pair< int, int >	checkAngleDeviation (double dPhi_1, double dPhi_2) const
void	estimateMomentum (const MuonRecHitContainer &validSet, CLHEP::Hep3Vector &momentum, int &charge) const
std::vector< SeedCandidate >	fillSeedCandidates (std::vector< MuonRecHitContainer > &allValidSets)
std::vector< MuonRecHitContainer >	findAllValidSets (const std::vector< MuonRecHitContainer > &MuonRecHitContainer_perLayer)
void	limitCombinatorics (std::vector< MuonRecHitContainer > &MuonRecHitContainer_perLayer)
TrajectorySeed	makeSeed (const TrajectoryStateOnSurface &tsos, const TransientTrackingRecHit::ConstRecHitContainer &hits) const
void	pre_prune (MuonRecHitContainer &validSet) const
virtual void	seeds (const MuonRecHitContainer &cluster, std::vector< TrajectorySeed > &result)
virtual void	setBField (const MagneticField *field)
	ignore - uses MuonServiceProxy More...
	SETSeedFinder (const edm::ParameterSet &pset)
void	setServiceProxy (MuonServiceProxy *service)
std::vector< MuonRecHitContainer >	sortByLayer (MuonRecHitContainer &cluster) const
void	validSetsPrePruning (std::vector< MuonRecHitContainer > &allValidSets)
virtual	~SETSeedFinder ()
Public Member Functions inherited from MuonSeedVFinder
void	setBeamSpot (const GlobalVector &gv)
virtual	~MuonSeedVFinder ()

Private Attributes
bool	apply_prePruning
MuonServiceProxy *	theService
bool	useSegmentsInTrajectory

Additional Inherited Members
Protected Attributes inherited from MuonSeedVFinder
MuonSeedPtExtractor *	thePtExtractor

Detailed Description

I. Bloch, E. James, S. Stoynev

Definition at line 12 of file SETSeedFinder.h.

Member Typedef Documentation

typedef MuonTransientTrackingRecHit::MuonRecHitContainer SETSeedFinder::MuonRecHitContainer

Definition at line 15 of file SETSeedFinder.h.

Constructor & Destructor Documentation

SETSeedFinder::SETSeedFinder ( const edm::ParameterSet &  pset )

explicit

Definition at line 24 of file SETSeedFinder.cc.

References apply_prePruning, edm::ParameterSet::getParameter(), and MuonSeedVFinder::thePtExtractor.

: MuonSeedVFinder()
{
  // Parameter set for the Builder
  ParameterSet trajectoryBuilderParameters = parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters");
  apply_prePruning = trajectoryBuilderParameters.getParameter<bool>("Apply_prePruning");
  // load pT seed parameters
  thePtExtractor = new MuonSeedPtExtractor(trajectoryBuilderParameters);

} 

virtual SETSeedFinder::~SETSeedFinder ( )

inlinevirtual

Definition at line 18 of file SETSeedFinder.h.

References MuonSeedVFinder::thePtExtractor.

{delete thePtExtractor;}

Member Function Documentation

std::pair< int, int > SETSeedFinder::checkAngleDeviation	(	double	dPhi_1,
		double	dPhi_2
	)		const

Definition at line 317 of file SETSeedFinder.cc.

References VarParsing::mult.

Referenced by pre_prune().

{
  // Two conditions:
  // a) deviations should be only to one side (above some absolute value cut to avoid
  //    material effects; this should be refined)
  // b) deviatiation in preceding steps should be bigger due to higher magnetic field
  //    (again - a minimal value cut should be in place; this also should account for
  //     the small (Z) distances in overlaping CSC chambers)

  double mult = dPhi_1 * dPhi_2;
  int signVal = 1;
  if(fabs(dPhi_1)<fabs(dPhi_2)){
    signVal = -1;
  }
  int signMult = -1;
  if(mult>0) signMult = 1;
  std::pair <int, int> sign;
  sign = make_pair (signVal, signMult);

  return sign;
}

void SETSeedFinder::estimateMomentum	(	const MuonRecHitContainer &	validSet,
		CLHEP::Hep3Vector &	momentum,
		int &	charge
	)		const

Definition at line 484 of file SETSeedFinder.cc.

References DeDxDiscriminatorTools::charge(), pat::helper::ParametrizationHelper::dimension(), python.cmstools::loop(), Gflash::par, funct::pow(), MuonSeedPtExtractor::pT_extract(), and MuonSeedVFinder::thePtExtractor.

Referenced by fillSeedCandidates().

{
  int firstMeasurement = -1;
  int lastMeasurement = -1;

  // don't use 2D measurements for momentum estimation 

  //if( 4==allValidSets[iSet].front()->dimension() &&
  //(allValidSets[iSet].front()->isCSC() || allValidSets[iSet].front()->isDT())){
  //firstMeasurement = 0;
  //}
  //else{
  // which is the "first" hit (4D)?
  for(unsigned int iMeas = 0;iMeas<validSet.size();++iMeas){
    if(4==validSet[iMeas]->dimension() &&
       (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT())){
      firstMeasurement = iMeas;
      break;
    }
  }
    //}

  std::vector<double> momentum_estimate;
  double pT = 0.;
  MuonTransientTrackingRecHit::ConstMuonRecHitPointer  firstHit;
  MuonTransientTrackingRecHit::ConstMuonRecHitPointer  secondHit;
  // which is the second hit?
  for(int loop = 0; loop<2; ++loop){// it is actually not used; to be removed
    // this is the last measurement
    if(!loop){// this is what is used currently
      // 23.04.09 : it becomes a problem with introduction of ME42 chambers -
      // the initial pT parametrization is incorrect for them
      for(int iMeas = validSet.size()-1;iMeas>-1;--iMeas){
        if(4==validSet[iMeas]->dimension() &&
           (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT()) &&
    // below is a fix saying "don't use ME4 chambers for initial pT estimation";
    // not using ME41 should not be a big loss too (and is more "symmetric" solution)
      fabs(validSet[iMeas]->globalPosition().z())<1000.){
          lastMeasurement = iMeas;
          break;
        }
      }
    }
    else{
      // this is the second measurement
      for(unsigned int iMeas = 1;iMeas<validSet.size();++iMeas){
        if(4==validSet[iMeas]->dimension() &&
           (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT())){
          lastMeasurement = iMeas;
          break;
        }
      }
    }
    // only 2D measurements (it should have been already abandoned)
    if(-1==lastMeasurement && -1==firstMeasurement){
       firstMeasurement = 0;
       lastMeasurement = validSet.size()-1;
    }
    // because of the ME42 above lastMeasurement could be -1
    else if(-1==lastMeasurement){
      lastMeasurement = firstMeasurement;
    }
   else if(-1==firstMeasurement){
     firstMeasurement = lastMeasurement;    
   }

    firstHit = validSet[firstMeasurement];
    secondHit = validSet[lastMeasurement];
      if(firstHit->isRPC() && secondHit->isRPC()){ // remove all RPCs from here?
        momentum_estimate.push_back(300.);
     momentum_estimate.push_back(300.);
    }
    else{
      if(firstHit->isRPC()){
        firstHit = secondHit;
      }
      else if(secondHit->isRPC()){
        secondHit  = firstHit;
      }
      //---- estimate pT given two hits
      //std::cout<<"   hits for initial pT estimate: first -> dim = "<<firstHit->dimension()<<" pos = "<<firstHit->globalPosition()<<
      //" , second -> "<<" dim = "<<secondHit->dimension()<<" pos = "<<secondHit->globalPosition()<<std::endl;
     //---- pT throws exception if hits are MB4 
     // (no coding for them - 2D hits in the outer station) 
     if(2==firstHit->dimension() && 2==secondHit->dimension()){
        momentum_estimate.push_back(999999999.);
        momentum_estimate.push_back(999999999.);
      }
      else{
        momentum_estimate = thePtExtractor->pT_extract(firstHit, secondHit);
      } 
    }
    pT = fabs(momentum_estimate[0]);
    if(1 || pT>40.){ //it is skipped; we have to look at least into number of hits in the chamber actually...
      // and then decide which segment to use
      // use the last measurement, otherwise use the second; this is to be investigated
      break;
    }
  }

  const float pT_min = 1.99;// many hardcoded - remove them!
  if(pT>3000.){
    pT=3000.;
  }
  else if(pT<pT_min ){
    if(pT>0){
      pT=pT_min ;
    }
    else if(pT>(-1)*pT_min ){
      pT=(-1)*pT_min ;
    }
    else if (pT<-3000.){
      pT= -3000;
    }
  }
  //std::cout<<"  THE pT from the parametrization: "<<momentum_estimate[0]<<std::endl;
  // estimate the charge of the track candidate from the delta phi of two segments:
  //int charge      = dPhi > 0 ? 1 : -1; // what we want is: dphi < 0 => charge = -1
  charge =  momentum_estimate[0]> 0 ? 1 : -1;

  // we have the pT - get the 3D momentum estimate as well

  // this is already final info:
  double xHit     = validSet[firstMeasurement]->globalPosition().x();
  double yHit     = validSet[firstMeasurement]->globalPosition().y();
  double rHit     = TMath::Sqrt(pow(xHit,2) + pow(yHit,2));

  double thetaInner = validSet[firstMeasurement]->globalPosition().theta();
  // if some of the segments is missing r-phi measurement then we should
  // use only the 4D phi estimate (also use 4D eta estimate only)
  // the direction is not so important (it will be corrected)

  double rTrack   = (pT /(0.3*3.8))*100.; //times 100 for conversion to cm!

  double par      = -1.*(2./charge)*(TMath::ASin(rHit/(2*rTrack)));
  double sinPar     = TMath::Sin( par );
  double cosPar     = TMath::Cos( par );

  // calculate phi at coordinate origin (0,0,0).
  double sinPhiH  = 1./(2.*charge*rTrack)*(xHit + ((sinPar)/(cosPar-1.))*yHit);
  double cosPhiH  = -1./(2.*charge*rTrack)*(((sinPar)/(1.-cosPar))*xHit + yHit);

  // finally set the return vector

  // try out the reco info:
  momEstimate = CLHEP::Hep3Vector(pT*cosPhiH, pT*sinPhiH, pT/TMath::Tan(thetaInner)); // should used into to theta directly here (rather than tan(atan2(...)))
  //Hep3Vector momEstimate(6.97961,      5.89732,     -50.0855);
  const float minMomenum = 5.; //hardcoded - remove it! same in SETFilter
  if (momEstimate.mag()<minMomenum){
    int sign = (pT<0.) ? -1: 1;
    pT = sign * (fabs(pT)+1);
    CLHEP::Hep3Vector momEstimate2(pT*cosPhiH, pT*sinPhiH, pT/TMath::Tan(thetaInner));
    momEstimate = momEstimate2;
    if (momEstimate.mag()<minMomenum){
      pT = sign * (fabs(pT)+1);
      CLHEP::Hep3Vector momEstimate3(pT*cosPhiH, pT*sinPhiH, pT/TMath::Tan(thetaInner));
      momEstimate = momEstimate3;
      if (momEstimate.mag()<minMomenum){
        pT = sign * (fabs(pT)+1);
        CLHEP::Hep3Vector momEstimate4(pT*cosPhiH, pT*sinPhiH, pT/TMath::Tan(thetaInner));
        momEstimate = momEstimate4;
      }
    }
  }
}

std::vector< SeedCandidate > SETSeedFinder::fillSeedCandidates

(

std::vector< MuonRecHitContainer > &

allValidSets

)

Definition at line 451 of file SETSeedFinder.cc.

References SeedCandidate::charge, estimateMomentum(), SeedCandidate::momentum, and SeedCandidate::theSet.

Referenced by SETMuonSeedProducer::produce().

                                                                  {
  //---- we have the valid sets constructed; transform the information in an
  //---- apropriate form; meanwhile - estimate the momentum for a given set

  // RPCs should not be used (no parametrization)
  std::vector <SeedCandidate> seedCandidates_inCluster;
  // calculate and fill the inputs needed
  // loop over all valid sets
  for(unsigned int iSet = 0;iSet<allValidSets.size();++iSet){
    //
    //std::cout<<"  This is SET number : "<<iSet<<std::endl; 
    //for(unsigned int iHit = 0;iHit<allValidSets[iSet].size();++iHit){
    //std::cout<<"   measurements in the SET:  iHit = "<<iHit<<" pos = "<<allValidSets[iSet][iHit]->globalPosition()<<
    //" dim = "<<allValidSets[iSet][iHit]->dimension()<<std::endl;
    //}


    CLHEP::Hep3Vector momEstimate;
    int chargeEstimate;
    estimateMomentum(allValidSets[iSet], momEstimate, chargeEstimate);
    MuonRecHitContainer MuonRecHitContainer_theSet_prep;
    // currently hardcoded - will be in proper loop of course:

    SeedCandidate seedCandidates_inCluster_prep;
    seedCandidates_inCluster_prep.theSet   = allValidSets[iSet];
    seedCandidates_inCluster_prep.momentum = momEstimate;
    seedCandidates_inCluster_prep.charge   = chargeEstimate;
    seedCandidates_inCluster.push_back(seedCandidates_inCluster_prep);
    // END estimateMomentum
  }
  return seedCandidates_inCluster;
}

std::vector< SETSeedFinder::MuonRecHitContainer > SETSeedFinder::findAllValidSets

(

const std::vector< MuonRecHitContainer > &

MuonRecHitContainer_perLayer

)

Definition at line 166 of file SETSeedFinder.cc.

References findQualityFiles::size.

Referenced by SETMuonSeedProducer::produce().

{
  std::vector <MuonRecHitContainer> allValidSets;
  // build all possible combinations (i.e valid sets; the algorithm name is after this feature -
  // SET algorithm)
  //
  // ugly... use recursive function?!
  // or implement Ingo's suggestion (a la ST)
  unsigned nLayers = MuonRecHitContainer_perLayer.size();
  if(1==nLayers){
    return allValidSets;
  }
  MuonRecHitContainer validSet;
  unsigned int iPos0 = 0;
  std::vector <unsigned int> iLayer(12);// could there be more than 11 layers?
  std::vector <unsigned int> size(12);
  if(iPos0<nLayers){
    size.at(iPos0) =  MuonRecHitContainer_perLayer.at(iPos0).size();
    for(iLayer[iPos0] = 0; iLayer[iPos0]<size[iPos0];++iLayer[iPos0]){
      validSet.clear();
      validSet.push_back(MuonRecHitContainer_perLayer[iPos0][iLayer[iPos0]]);
      unsigned int iPos1 = 1;
      if(iPos1<nLayers){
        size.at(iPos1) =  MuonRecHitContainer_perLayer.at(iPos1).size();
        for(iLayer[iPos1] = 0; iLayer[iPos1]<size[iPos1];++iLayer[iPos1]){
          validSet.resize(iPos1);
          validSet.push_back(MuonRecHitContainer_perLayer[iPos1][iLayer[iPos1]]);
          unsigned int iPos2 = 2;
          if(iPos2<nLayers){
            size.at(iPos2) =  MuonRecHitContainer_perLayer.at(iPos2).size();
            for(iLayer[iPos2] = 0; iLayer[iPos2]<size[iPos2];++iLayer[iPos2]){
              validSet.resize(iPos2);
              validSet.push_back(MuonRecHitContainer_perLayer[iPos2][iLayer[iPos2]]);
              unsigned int iPos3 = 3;
              if(iPos3<nLayers){
                size.at(iPos3) =  MuonRecHitContainer_perLayer.at(iPos3).size();
                for(iLayer[iPos3] = 0; iLayer[iPos3]<size[iPos3];++iLayer[iPos3]){
                  validSet.resize(iPos3);
                  validSet.push_back(MuonRecHitContainer_perLayer[iPos3][iLayer[iPos3]]);
                  unsigned int iPos4 = 4;
                  if(iPos4<nLayers){
                    size.at(iPos4) =  MuonRecHitContainer_perLayer.at(iPos4).size();
                    for(iLayer[iPos4] = 0; iLayer[iPos4]<size[iPos4];++iLayer[iPos4]){
                      validSet.resize(iPos4);
                      validSet.push_back(MuonRecHitContainer_perLayer[iPos4][iLayer[iPos4]]);
                      unsigned int iPos5 = 5;
                      if(iPos5<nLayers){
                        size.at(iPos5) =  MuonRecHitContainer_perLayer.at(iPos5).size();
                        for(iLayer[iPos5] = 0; iLayer[iPos5]<size[iPos5];++iLayer[iPos5]){
                          validSet.resize(iPos5);
                          validSet.push_back(MuonRecHitContainer_perLayer[iPos5][iLayer[iPos5]]);
                          unsigned int iPos6 = 6;
                          if(iPos6<nLayers){
                            size.at(iPos6) =  MuonRecHitContainer_perLayer.at(iPos6).size();
                            for(iLayer[iPos6] = 0; iLayer[iPos6]<size[iPos6];++iLayer[iPos6]){
                              validSet.resize(iPos6);
                              validSet.push_back(MuonRecHitContainer_perLayer[iPos6][iLayer[iPos6]]);
                              unsigned int iPos7 = 7;
                              if(iPos7<nLayers){
                                size.at(iPos7) =  MuonRecHitContainer_perLayer.at(iPos7).size();
                                for(iLayer[iPos7] = 0; iLayer[iPos7]<size[iPos7];++iLayer[iPos7]){
                                  validSet.resize(iPos7);
                                  validSet.push_back(MuonRecHitContainer_perLayer[iPos7][iLayer[iPos7]]);
                                  unsigned int iPos8 = 8;
                                  if(iPos8<nLayers){
                                    size.at(iPos8) =  MuonRecHitContainer_perLayer.at(iPos8).size();
                                    for(iLayer[iPos8] = 0; iLayer[iPos8]<size[iPos8];++iLayer[iPos8]){
                                      validSet.resize(iPos8);
                                      validSet.push_back(MuonRecHitContainer_perLayer[iPos8][iLayer[iPos8]]);
                                      unsigned int iPos9 = 9;
                                      if(iPos9<nLayers){
                                        size.at(iPos9) =  MuonRecHitContainer_perLayer.at(iPos9).size();
                                        for(iLayer[iPos9] = 0; iLayer[iPos9]<size[iPos9];++iLayer[iPos9]){
                                          validSet.resize(iPos9);
                                          validSet.push_back(MuonRecHitContainer_perLayer[iPos9][iLayer[iPos9]]);
                                          unsigned int iPos10 = 10;
                                          if(iPos10<nLayers){
                                            size.at(iPos10) =  MuonRecHitContainer_perLayer.at(iPos10).size();
                                            for(iLayer[iPos10] = 0; iLayer[iPos10]<size[iPos10];++iLayer[iPos10]){
                                              validSet.resize(iPos10);
                                              validSet.push_back(MuonRecHitContainer_perLayer[iPos10][iLayer[iPos10]]);
                                              unsigned int iPos11 = 11;// more?
                                              if(iPos11<nLayers){
                                                size.at(iPos11) =  MuonRecHitContainer_perLayer.at(iPos11).size();
                                                for(iLayer[iPos11] = 0; iLayer[iPos11]<size[iPos11];++iLayer[iPos11]){
                                               }
                                              }
                                              else{
                                                allValidSets.push_back(validSet);

                                              }
                                            }
                                          }
                                          else{
                                            allValidSets.push_back(validSet);
                                          }
                                        }
                                      }
                                      else{
                                        allValidSets.push_back(validSet);
                                      }
                                    }
                                  }
                                  else{
                                    allValidSets.push_back(validSet);
                                  }
                                }
                              }
                              else{
                                allValidSets.push_back(validSet);
                              }
                            }
                          }
                          else{
                            allValidSets.push_back(validSet);
                          }
                        }
                      }
                      else{
                        allValidSets.push_back(validSet);
                      }
                    }
                  }
                  else{
                    allValidSets.push_back(validSet);
                  }
                }
              }
              else{
                allValidSets.push_back(validSet);
              }
            }
          }
          else{
            allValidSets.push_back(validSet);
          }
        }
      }
      else{
        allValidSets.push_back(validSet);
      }
    }
  }
  else{
    allValidSets.push_back(validSet);
  }
  return allValidSets;
}

void SETSeedFinder::limitCombinatorics

(

std::vector< MuonRecHitContainer > &

MuonRecHitContainer_perLayer

)

Definition at line 129 of file SETSeedFinder.cc.

References i, and mathSSE::return().

Referenced by SETMuonSeedProducer::produce().

                                                                                                     {
  const int maximumNumberOfCombinations = 1000000;
  unsigned nLayers = MuonRecHitContainer_perLayer.size();
  if(1==nLayers){
    return ;
  }
  // maximal number of (segment) layers would be upto ~12; see next function
  // below is just a quick fix for a rare "overflow"
  if(MuonRecHitContainer_perLayer.size()>15){
    MuonRecHitContainer_perLayer.resize(1);
    return;
  }
        
  std::vector <double> sizeOfLayer(nLayers);
  //std::cout<<" nLayers = "<<nLayers<<std::endl;
  double nAllCombinations = 1.;
  for(unsigned int i = 0;i<nLayers;++i ){
    //std::cout<<" i = "<<i<<" size = "<<MuonRecHitContainer_perLayer.at(i).size()<<std::endl;
    sizeOfLayer.at(i) = MuonRecHitContainer_perLayer.at(i).size();
    nAllCombinations*=MuonRecHitContainer_perLayer.at(i).size();
  }
  //std::cout<<"nAllCombinations = "<<nAllCombinations<<std::endl;
  //---- Erase most busy detector layers until we get less than maximumNumberOfCombinations combinations
  int iCycle = 0;
  while(nAllCombinations > float(maximumNumberOfCombinations)){
    ++iCycle;
    std::vector <double>::iterator maxEl_it = max_element(sizeOfLayer.begin(),sizeOfLayer.end());
    int maxEl = maxEl_it - sizeOfLayer.begin();
    nAllCombinations/=MuonRecHitContainer_perLayer.at(maxEl).size();
    //std::cout<<" iCycle = "<<iCycle<<" nAllCombinations = "<<nAllCombinations<<std::endl;
    MuonRecHitContainer_perLayer.erase(MuonRecHitContainer_perLayer.begin()+maxEl);
    sizeOfLayer.erase(sizeOfLayer.begin()+maxEl);
  }
  return;
}

TrajectorySeed SETSeedFinder::makeSeed	(	const TrajectoryStateOnSurface &	tsos,
		const TransientTrackingRecHit::ConstRecHitContainer &	hits
	)		const

Definition at line 653 of file SETSeedFinder.cc.

References alongMomentum, dir, oppositeToMomentum, TrajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), TrajectorySeed::recHits(), and useSegmentsInTrajectory.

{
  edm::OwnVector<TrackingRecHit> recHitsContainer;
  for(unsigned int iHit = 0;iHit < hits.size();++iHit){
    recHitsContainer.push_back(hits.at(iHit)->hit()->clone());
  }
  PropagationDirection dir = oppositeToMomentum;
  if(useSegmentsInTrajectory){
    dir = alongMomentum;// why forward (for rechits) later?
  }
  TrajectoryStateTransform tsTransform;
  PTrajectoryStateOnDet *seedTSOS =
  tsTransform.persistentState( firstTSOS, hits.at(0)->geographicalId().rawId());
  TrajectorySeed seed(*seedTSOS,recHitsContainer,dir);
  TrajectorySeed::range range = seed.recHits();

  //MuonPatternRecoDumper debug;
  //std::cout<<" firstTSOS = "<<debug.dumpTSOS(firstTSOS)<<std::endl;
  //std::cout<<" iTraj = ???"<<" hits = "<<range.second-range.first<<std::endl;
  //std::cout<<" nhits = "<<hits.size()<<std::endl;
  //for(unsigned int iRH=0;iRH<hits.size();++iRH){
  //std::cout<<" RH = "<<iRH+1<<" globPos = "<<hits.at(iRH)->globalPosition()<<std::endl;
  //}
  return seed;
}

void SETSeedFinder::pre_prune

(

SETSeedFinder::MuonRecHitContainer &

validSet

)

const

Definition at line 352 of file SETSeedFinder.cc.

References checkAngleDeviation(), prof2calltree::count, pat::helper::ParametrizationHelper::dimension(), dPhi(), and phi.

Referenced by validSetsPrePruning().

{
  unsigned nHits = validSet.size();
  if(nHits>3){ // to decide we need at least 4 measurements
    // any information could be used to make a decision for pruning
    // maybe dPhi (delta Phi) is enough
    std::vector <double> dPhi;
    double dPhi_tmp;
    bool wildCandidate;
    int pruneHit_tmp;

    for(unsigned int iHit = 1;iHit<nHits;++iHit){
      dPhi_tmp = validSet[iHit]->globalPosition().phi() -
                 validSet[iHit-1]->globalPosition().phi();
      dPhi.push_back(dPhi_tmp);
    }
    std::vector <int> pruneHit;
    //---- loop over all the hits in a set

    for(unsigned int iHit = 0;iHit<nHits;++iHit){
      double dPHI_MIN = 0.02;//?? hardcoded - remove it
      if(iHit){
        // if we have to remove the very first hit (iHit == 0) then
        // we'll probably be in trouble already
        wildCandidate = false;
        // actually 2D is bad only if not r-phi... Should I refine it?
        // a hit is a candidate for pruning only if dPhi > dPHI_MIN;
        // pruning decision is based on combination of hits characteristics
        if(4==validSet[iHit-1]->dimension() && 4 == validSet[iHit]->dimension() &&
           fabs(validSet[iHit]->globalPosition().phi() -
                validSet[iHit-1]->globalPosition().phi())>dPHI_MIN ){
          wildCandidate = true;
        }
        pruneHit_tmp = -1;
        if(wildCandidate){
          // OK - this couple doesn't look good (and is from 4D segments); proceed...
          if(1==iHit){// the first  and the last hits are special case
            if(4==validSet[iHit+1]->dimension() && 4 == validSet[iHit+2]->dimension()){//4D?
              // is the picture better if we remove the second hit?
              dPhi_tmp = validSet[iHit+1]->globalPosition().phi() -
                validSet[iHit-1]->globalPosition().phi();
              // is the deviation what we expect (sign, not magnitude)?
              std::pair <int, int> sign = checkAngleDeviation(dPhi_tmp, dPhi[2]);
              if( 1==sign.first && 1==sign.second){
                pruneHit_tmp = iHit;// mark the hit 1 for removing
              }
            }
          }
          else if(iHit>1 && iHit<validSet.size()-1){
            if(4 == validSet[0]->dimension() && // we rely on the first (most important) couple
               4 == validSet[1]->dimension() &&
               pruneHit.back()!=int(iHit-1) && pruneHit.back()!=1){// check if hits are already marked
              // decide which of the two hits should be removed (if any; preferably the outer one i.e.
              // iHit rather than iHit-1); here - check what we get by removing iHit
              dPhi_tmp = validSet[iHit+1]->globalPosition().phi() -
                validSet[iHit-1]->globalPosition().phi();
              // first couple is most important anyway so again compare to it
              std::pair <int, int> sign = checkAngleDeviation(dPhi[0],dPhi_tmp);
              if( 1==sign.first && 1==sign.second){
                pruneHit_tmp = iHit; // mark the hit iHit for removing
              }
              else{ // iHit is not to be removed; proceed...
                // what if we remove (iHit - 1) instead of iHit?
                dPhi_tmp = validSet[iHit+1]->globalPosition().phi() -
                  validSet[iHit]->globalPosition().phi();
                std::pair <int, int> sign = checkAngleDeviation(dPhi[0],dPhi_tmp);
                if( 1==sign.first && 1==sign.second){
                  pruneHit_tmp = iHit-1;// mark the hit (iHit -1) for removing
                }
              }
            }
          }
          else{
            // the last hit: if picture is not good - remove it
            if(pruneHit.size()>1 && pruneHit[pruneHit.size()-1]<0 && pruneHit[pruneHit.size()-2]<0){
              std::pair <int, int> sign = checkAngleDeviation(dPhi[dPhi.size()-2], dPhi[dPhi.size()-1]);
              if( -1==sign.first && -1==sign.second){// here logic is a bit twisted
                pruneHit_tmp = iHit; // mark the last hit for removing
              }
            }
          }
        }
        pruneHit.push_back(pruneHit_tmp);
      }
    }
      // }
     // actual pruning
     for(unsigned int iHit = 1;iHit<nHits;++iHit){
        int count = 0;
      if(pruneHit[iHit-1]>0){
        validSet.erase(validSet.begin()+pruneHit[iHit-1]-count);
        ++count;
      }
    }
  }
}

void SETSeedFinder::seeds ( const MuonRecHitContainer &  cluster, std::vector< TrajectorySeed > &  result  )

virtual

The container sent in is expected to be a cluster, which isn't the same as a pattern. A cluster can have more than one hit on a layer. Internally, this method splits the cluster into patterns, and chooses the best one via a chi2. But it calculates the trajectoryMeasurements at the same time, so we can't really separate the steps.

Implements MuonSeedVFinder.

Definition at line 36 of file SETSeedFinder.cc.

{
}

virtual void SETSeedFinder::setBField ( const MagneticField *  field )

inlinevirtual

ignore - uses MuonServiceProxy

Implements MuonSeedVFinder.

Definition at line 20 of file SETSeedFinder.h.

{}

void SETSeedFinder::setServiceProxy ( MuonServiceProxy *  service )

inline

Definition at line 32 of file SETSeedFinder.h.

References theService.

Referenced by SETMuonSeedProducer::SETMuonSeedProducer().

{theService = service;}

std::vector< SETSeedFinder::MuonRecHitContainer > SETSeedFinder::sortByLayer

(

MuonRecHitContainer &

cluster

)

const

Definition at line 55 of file SETSeedFinder.cc.

References MuonSubdetId::CSC, MuonSubdetId::DT, Plane::localZ(), mag(), CSCDetId::ring(), sortSegRadius, DTChamberId::station(), CSCDetId::station(), DetId::subdetId(), GeomDet::surface(), and theService.

Referenced by SETMuonSeedProducer::produce().

{
  stable_sort(cluster.begin(), cluster.end(),sortSegRadius);
    //---- group hits in detector layers (if in same layer); the idea is that
    //---- some hits could not belong to a track simultaneously - these will be in a
    //---- group; two hits from one and the same group will not go to the same track
  std::vector< MuonRecHitContainer > MuonRecHitContainer_perLayer;
  if(cluster.size()){
    int iHit =0;
    MuonRecHitContainer hitsInThisLayer;
    hitsInThisLayer.push_back(cluster[iHit]);
    DetId  detId = cluster[iHit]->hit()->geographicalId();
    const GeomDet* geomDet = theService->trackingGeometry()->idToDet( detId );
    while(iHit<int(cluster.size())-1){
      DetId  detId_2 = cluster[iHit+1]->hit()->geographicalId();
      const GlobalPoint gp_nextHit = cluster[iHit+1]->globalPosition();

      // this is the distance of the "second" hit to the "first" detector (containing the "first hit")
      float distanceToDetector = fabs(geomDet->surface().localZ(gp_nextHit));

      //---- hits from DT and CSC  could be very close in angle but incosistent with
      //---- belonging to a common track (and these are different surfaces);
      //---- also DT (and CSC now - 090822) hits from a station (in a pre-cluster) should be always in a group together;
      //---- take this into account and put such hits in a group together

      bool specialCase = false;
      if( detId.subdetId()   == MuonSubdetId::DT &&
          detId_2.subdetId() == MuonSubdetId::DT    ){
        DTChamberId dtCh(detId);
        DTChamberId dtCh_2(detId_2);
        specialCase =  (dtCh.station() == dtCh_2.station());
      }
      else if(detId.subdetId()   == MuonSubdetId::CSC &&
              detId_2.subdetId()   == MuonSubdetId::CSC){
        CSCDetId cscCh(detId);
        CSCDetId cscCh_2(detId_2);
        specialCase = (cscCh.station() == cscCh_2.station() && cscCh.ring() == cscCh_2.ring());
      } 

      if(distanceToDetector<0.001 || true==specialCase){ // hardcoded value - remove!
        hitsInThisLayer.push_back(cluster[iHit+1]);
      }
      else{
        specialCase = false;
        if(( (cluster[iHit]->isDT() &&
              cluster[iHit+1]->isCSC()) ||
             (cluster[iHit]->isCSC() &&
              cluster[iHit+1]->isDT())) &&
         //---- what is the minimal distance between a DT and a CSC hit belonging
         //---- to a common track? (well, with "reasonable" errors; put 10 cm for now)
           fabs(cluster[iHit+1]->globalPosition().mag() -
                cluster[iHit]->globalPosition().mag())<10.){
          hitsInThisLayer.push_back(cluster[iHit+1]);
          // change to Stoyan - now we also update the detID here... give it a try. IBL 080905
          detId = cluster[iHit+1]->hit()->geographicalId();
          geomDet = theService->trackingGeometry()->idToDet( detId );
        }
        else if(!specialCase){
          //---- put the group of hits in the vector (containing the groups of hits)
          //---- and continue with next layer (group)
          MuonRecHitContainer_perLayer.push_back(hitsInThisLayer);
          hitsInThisLayer.clear();
          hitsInThisLayer.push_back(cluster[iHit+1]);
          detId = cluster[iHit+1]->hit()->geographicalId();
          geomDet = theService->trackingGeometry()->idToDet( detId );
        }
      }
      ++iHit;
    }
    MuonRecHitContainer_perLayer.push_back(hitsInThisLayer);
  }
  return MuonRecHitContainer_perLayer;
}

void SETSeedFinder::validSetsPrePruning

(

std::vector< MuonRecHitContainer > &

allValidSets

)

Definition at line 340 of file SETSeedFinder.cc.

References pre_prune().

Referenced by SETMuonSeedProducer::produce().

{
  //---- this actually is a pre-pruning; it does not include any fit information;
  //---- it is intended to remove only very "wild" segments from a set;
  //---- no "good" segment is to be lost (otherwise - widen the parameters)

  for(unsigned int iSet = 0;iSet<allValidSets.size();++iSet){
    pre_prune(allValidSets[iSet]);
  }
}

Member Data Documentation

bool SETSeedFinder::apply_prePruning

private

Definition at line 61 of file SETSeedFinder.h.

Referenced by SETSeedFinder().

MuonServiceProxy* SETSeedFinder::theService

private

Definition at line 59 of file SETSeedFinder.h.

Referenced by setServiceProxy(), and sortByLayer().

bool SETSeedFinder::useSegmentsInTrajectory

private

Definition at line 62 of file SETSeedFinder.h.

Referenced by makeSeed().