ME0SegAlgoMM Class Reference

#include <ME0SegAlgoMM.h>

Inheritance diagram for ME0SegAlgoMM:

Public Types
typedef std::vector< const ME0RecHit * >	EnsembleHitContainer
	Typedefs. More...
typedef std::vector < EnsembleHitContainer >	ProtoSegments
Public Types inherited from ME0SegmentAlgorithm
typedef std::pair< const ME0EtaPartition *, std::map < uint32_t, const ME0EtaPartition * > >	ME0Ensemble

Public Member Functions
	ME0SegAlgoMM (const edm::ParameterSet &ps)
	Constructor. More...
std::vector< ME0Segment >	run (const ME0Ensemble &ensemble, const EnsembleHitContainer &rechits)
virtual	~ME0SegAlgoMM ()
	Destructor. More...
Public Member Functions inherited from ME0SegmentAlgorithm
	ME0SegmentAlgorithm (const edm::ParameterSet &)
	Constructor. More...
virtual std::vector< ME0Segment >	run (const ME0Ensemble &ensemble, const std::vector< const ME0RecHit * > &rechits)=0
virtual	~ME0SegmentAlgorithm ()
	Destructor. More...

Private Member Functions
std::vector< ME0Segment >	buildSegments (const EnsembleHitContainer &rechits)
ProtoSegments	chainHits (const EnsembleHitContainer &rechits)
ProtoSegments	clusterHits (const EnsembleHitContainer &rechits)
	Utility functions. More...
bool	isGoodToMerge (const EnsembleHitContainer &newChain, const EnsembleHitContainer &oldChain)

Private Attributes
bool	debug
double	dEtaChainBoxMax
double	dPhiChainBoxMax
double	dTimeChainBoxMax
double	dXclusBoxMax
double	dYclusBoxMax
int	maxRecHitsInCluster
unsigned int	minHitsPerSegment
const std::string	myName
bool	preClustering
bool	preClustering_useChaining
EnsembleHitContainer	proto_segment
std::unique_ptr< ME0SegFit >	sfit_
ME0Ensemble	theEnsemble

Static Private Attributes
static float	running_max =999999.

Detailed Description

This algorithm is very basic no attemp to deal with ambiguities , noise etc. The ME0 track segments is built out of the rechit's in a the 6 ME0 Layer denoted as the ME0 Ensabmle .

Authors

Marcello Maggi

Definition at line 24 of file ME0SegAlgoMM.h.

Member Typedef Documentation

typedef std::vector<const ME0RecHit*> ME0SegAlgoMM::EnsembleHitContainer

Typedefs.

Definition at line 30 of file ME0SegAlgoMM.h.

typedef std::vector<EnsembleHitContainer> ME0SegAlgoMM::ProtoSegments

Definition at line 31 of file ME0SegAlgoMM.h.

Constructor & Destructor Documentation

ME0SegAlgoMM::ME0SegAlgoMM ( const edm::ParameterSet &  ps )

explicit

Constructor.

Definition at line 25 of file ME0SegAlgoMM.cc.

References debug, dEtaChainBoxMax, dPhiChainBoxMax, dTimeChainBoxMax, dXclusBoxMax, dYclusBoxMax, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maxRecHitsInCluster, minHitsPerSegment, preClustering, and preClustering_useChaining.

                                                    : ME0SegmentAlgorithm(ps), myName("ME0SegAlgoMM")
{
  debug                     = ps.getUntrackedParameter<bool>("ME0Debug");
  minHitsPerSegment         = ps.getParameter<unsigned int>("minHitsPerSegment");
  preClustering             = ps.getParameter<bool>("preClustering");
  dXclusBoxMax              = ps.getParameter<double>("dXclusBoxMax");
  dYclusBoxMax              = ps.getParameter<double>("dYclusBoxMax");
  preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
  dPhiChainBoxMax           = ps.getParameter<double>("dPhiChainBoxMax");
  dEtaChainBoxMax           = ps.getParameter<double>("dEtaChainBoxMax");
  dTimeChainBoxMax          = ps.getParameter<double>("dTimeChainBoxMax");
  maxRecHitsInCluster       = ps.getParameter<int>("maxRecHitsInCluster");

  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::ctor] Parameters to build segments :: "
                   << "preClustering = "<<preClustering<<" preClustering_useChaining = "<<preClustering_useChaining
                   <<" dPhiChainBoxMax = "<<dPhiChainBoxMax<<" dEtaChainBoxMax = "<<dEtaChainBoxMax<<" dTimeChainBoxMax = "<<dTimeChainBoxMax
                   <<" minHitsPerSegment = "<<minHitsPerSegment<<" maxRecHitsInCluster = "<<maxRecHitsInCluster;
}

ME0SegAlgoMM::~ME0SegAlgoMM ( )

virtual

Destructor.

Definition at line 47 of file ME0SegAlgoMM.cc.

                            {
}

Member Function Documentation

std::vector< ME0Segment > ME0SegAlgoMM::buildSegments ( const EnsembleHitContainer &  rechits )

private

Definition at line 323 of file ME0SegAlgoMM.cc.

References minHitsPerSegment, funct::pow(), proto_segment, sfit_, mathSSE::sqrt(), theEnsemble, and tmp.

Referenced by run().

                                                                                     {
  std::vector<ME0Segment> me0segs;

  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] will now try to fit a ME0Segment from collection of "<<rechits.size()<<" ME0 RecHits";
  #ifdef EDM_ML_DEBUG // have lines below only compiled when in debug mode 
    for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
      auto me0id = (*rh)->me0Id();
      auto rhLP = (*rh)->localPosition();
      edm::LogVerbatim("ME0SegAlgoMM") << "[RecHit :: Loc x = "<<std::showpos<<std::setw(9)<<rhLP.x()<<" Loc y = "<<std::showpos<<std::setw(9)<<rhLP.y()
                       <<" Time = "<<std::showpos<<(*rh)->tof()<<" -- "<<me0id.rawId()<<" = "<<me0id<<" ]";
      }
  #endif

  proto_segment.clear();
  // select hits from the ensemble and sort it 
  for (auto rh=rechits.begin(); rh!=rechits.end();rh++){
    proto_segment.push_back(*rh);
  }
  if (proto_segment.size() < minHitsPerSegment){
    return me0segs;
  }

  // The actual fit on all hits of the vector of the selected Tracking RecHits:
  sfit_ = std::unique_ptr<ME0SegFit>(new ME0SegFit(theEnsemble.second, rechits));
  sfit_->fit();
  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] ME0Segment fit done";

  // obtain all information necessary to make the segment:
  LocalPoint protoIntercept      = sfit_->intercept();
  LocalVector protoDirection     = sfit_->localdir();
  AlgebraicSymMatrix protoErrors = sfit_->covarianceMatrix(); 
  double protoChi2               = sfit_->chi2();

  // Calculate the central value and uncertainty of the segment time
  float averageTime=0.;
  for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
    averageTime += (*rh)->tof();                                          
  }
  if(rechits.size() != 0) averageTime=averageTime/(rechits.size());
  float timeUncrt=0.;
  for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
    timeUncrt += pow((*rh)->tof()-averageTime,2);
  }
  if(rechits.size() > 1) timeUncrt=timeUncrt/(rechits.size()-1);
  timeUncrt = sqrt(timeUncrt);

  // save all information inside GEMCSCSegment
  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] will now try to make ME0Segment from collection of "<<rechits.size()<<" ME0 RecHits";
  ME0Segment tmp(proto_segment, protoIntercept, protoDirection, protoErrors, protoChi2, averageTime, timeUncrt);

  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] ME0Segment made";
  edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegAlgoMM::buildSegments] "<<tmp;

  me0segs.push_back(tmp);
  return me0segs;
}

ME0SegAlgoMM::ProtoSegments ME0SegAlgoMM::chainHits ( const EnsembleHitContainer &  rechits )

private

Definition at line 204 of file ME0SegAlgoMM.cc.

References begin, end, i, isGoodToMerge(), and groupFilesInBlocks::temp.

Referenced by run().

                                                           {

  ProtoSegments rechits_chains; 
  EnsembleHitContainer temp;
  ProtoSegments seeds;

  std::vector <bool> usedCluster;

  // split rechits into subvectors and return vector of vectors:
  // Loop over rechits
  // Create one seed per hit
  for(unsigned int i = 0; i < rechits.size(); ++i) {
    temp.clear();
    temp.push_back(rechits[i]);
    seeds.push_back(temp);
    usedCluster.push_back(false);
  }

  // merge chains that are too close ("touch" each other)
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    for(size_t MMM = NNN+1; MMM < seeds.size(); ++MMM) {
      if(usedCluster[MMM] || usedCluster[NNN]){
        continue;
      }
      // all is in the way we define "good";
      // try not to "cluster" the hits but to "chain" them;
      // it does the clustering but also does a better job
      // for inclined tracks (not clustering them together;
      // crossed tracks would be still clustered together) 
      // 22.12.09: In fact it is not much more different 
      // than the "clustering", we just introduce another
      // variable in the game - Z. And it makes sense 
      // to re-introduce Y (or actually wire group mumber)
      // in a similar way as for the strip number - see
      // the code below.
      bool goodToMerge  = isGoodToMerge(seeds[NNN], seeds[MMM]);
      if(goodToMerge){
        // merge chains!
        // merge by adding seed NNN to seed MMM and erasing seed NNN

        // add seed NNN to MMM (lower to larger number)
        seeds[MMM].insert(seeds[MMM].end(),seeds[NNN].begin(),seeds[NNN].end());

        // mark seed NNN as used
        usedCluster[NNN] = true;
        // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
        // next seed (NNN+1)
        break;
      }

    }
  }

  // hand over the final seeds to the output
  // would be more elegant if we could do the above step with
  // erasing the merged ones, rather than the

  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    if(usedCluster[NNN]) continue; //skip seeds that have been marked as used up in merging
    rechits_chains.push_back(seeds[NNN]);
  }

  //***************************************************************

  return rechits_chains;
}

ME0SegAlgoMM::ProtoSegments ME0SegAlgoMM::clusterHits ( const EnsembleHitContainer &  rechits )

private

Utility functions.

Definition at line 104 of file ME0SegAlgoMM.cc.

References begin, dXclusBoxMax, dYclusBoxMax, end, i, LogDebug, running_max, findQualityFiles::size, groupFilesInBlocks::temp, x, and y.

Referenced by run().

                                                             {

  ProtoSegments rechits_clusters; // this is a collection of groups of rechits

  float dXclus_box = 0.0;
  float dYclus_box = 0.0;

  EnsembleHitContainer temp;
  ProtoSegments seeds;

  std::vector<float> running_meanX;
  std::vector<float> running_meanY;

  std::vector<float> seed_minX;
  std::vector<float> seed_maxX;
  std::vector<float> seed_minY;
  std::vector<float> seed_maxY;

  // split rechits into subvectors and return vector of vectors:
  // Loop over rechits 
  // Create one seed per hit
  for(unsigned int i = 0; i < rechits.size(); ++i) {
    temp.clear();
    temp.push_back(rechits[i]);
    seeds.push_back(temp);

    // First added hit in seed defines the mean to which the next hit is compared
    // for this seed.

    running_meanX.push_back( rechits[i]->localPosition().x() );
    running_meanY.push_back( rechits[i]->localPosition().y() );
    
    // set min/max X and Y for box containing the hits in the precluster:
    seed_minX.push_back( rechits[i]->localPosition().x() );
    seed_maxX.push_back( rechits[i]->localPosition().x() );
    seed_minY.push_back( rechits[i]->localPosition().y() );
    seed_maxY.push_back( rechits[i]->localPosition().y() );
  }
    
  // merge clusters that are too close
  // measure distance between final "running mean"
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    for(size_t MMM = NNN+1; MMM < seeds.size(); ++MMM) {
      if(running_meanX[MMM] == running_max || running_meanX[NNN] == running_max ) {
    LogDebug("ME0SegAlgoMM") << "[ME0SegAlgoMM::clusterHits]: ALARM! Skipping used seeds, this should not happen - inform developers!";
    continue; //skip seeds that have been used 
      }
      
      // calculate cut criteria for simple running mean distance cut:
      //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
      //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);
      // calculate minmal distance between precluster boxes containing the hits:
      if ( running_meanX[NNN] > running_meanX[MMM] ) dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
      else                                           dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
      if ( running_meanY[NNN] > running_meanY[MMM] ) dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
      else                                           dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
      
      
      if( dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax ) {
    // merge clusters!
    // merge by adding seed NNN to seed MMM and erasing seed NNN
        
    // calculate running mean for the merged seed:
    if(seeds[NNN].size()+seeds[MMM].size() != 0) {
      running_meanX[MMM] = (running_meanX[NNN]*seeds[NNN].size() + running_meanX[MMM]*seeds[MMM].size()) / (seeds[NNN].size()+seeds[MMM].size());
      running_meanY[MMM] = (running_meanY[NNN]*seeds[NNN].size() + running_meanY[MMM]*seeds[MMM].size()) / (seeds[NNN].size()+seeds[MMM].size());
    }

    // update min/max X and Y for box containing the hits in the merged cluster:
    if ( seed_minX[NNN] <  seed_minX[MMM] ) seed_minX[MMM] = seed_minX[NNN];
    if ( seed_maxX[NNN] >  seed_maxX[MMM] ) seed_maxX[MMM] = seed_maxX[NNN];
    if ( seed_minY[NNN] <  seed_minY[MMM] ) seed_minY[MMM] = seed_minY[NNN];
    if ( seed_maxY[NNN] >  seed_maxY[MMM] ) seed_maxY[MMM] = seed_maxY[NNN];
        
    // add seed NNN to MMM (lower to larger number)
    seeds[MMM].insert(seeds[MMM].end(),seeds[NNN].begin(),seeds[NNN].end());
        
    // mark seed NNN as used (at the moment just set running mean to 999999.)
    running_meanX[NNN] = running_max;
    running_meanY[NNN] = running_max;
    // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to 
    // next seed (NNN+1)
    break;
      }
    }
  }

  // hand over the final seeds to the output
  // would be more elegant if we could do the above step with 
  // erasing the merged ones, rather than the 
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    if(running_meanX[NNN] == running_max) continue; //skip seeds that have been marked as used up in merging
    rechits_clusters.push_back(seeds[NNN]);
  }

  return rechits_clusters; 
}

bool ME0SegAlgoMM::isGoodToMerge ( const EnsembleHitContainer &  newChain, const EnsembleHitContainer &  oldChain  )

private

Definition at line 271 of file ME0SegAlgoMM.cc.

References funct::abs(), reco::deltaPhi(), dEtaChainBoxMax, dPhiChainBoxMax, dTimeChainBoxMax, PV3DBase< T, PVType, FrameType >::eta(), PV3DBase< T, PVType, FrameType >::phi(), and theEnsemble.

Referenced by chainHits().

                                                                                                           {

  std::vector<float> phi_new, eta_new, time_new, phi_old, eta_old, time_old;
  std::vector<int> layer_new, layer_old;

  for(size_t iRH_new = 0;iRH_new<newChain.size();++iRH_new){
    GlobalPoint pos_new = theEnsemble.first->toGlobal(newChain[iRH_new]->localPosition());
    layer_new.push_back(newChain[iRH_new]->me0Id().layer());
    phi_new.push_back(pos_new.phi());
    eta_new.push_back(pos_new.eta());
    time_new.push_back(newChain[iRH_new]->tof());
  }  
  for(size_t iRH_old = 0;iRH_old<oldChain.size();++iRH_old){
    GlobalPoint pos_old = theEnsemble.first->toGlobal(oldChain[iRH_old]->localPosition());
    layer_old.push_back(oldChain[iRH_old]->me0Id().layer());
    phi_old.push_back(pos_old.phi());
    eta_old.push_back(pos_old.eta());
    time_old.push_back(oldChain[iRH_old]->tof());
  }

  for(size_t jRH_new = 0; jRH_new<phi_new.size(); ++jRH_new){
    for(size_t jRH_old = 0; jRH_old<phi_old.size(); ++jRH_old){

      // to be chained, two hits need to be in neighbouring layers...
      // or better allow few missing layers (upto 3 to avoid inefficiencies);
      // however we'll not make an angle correction because it
      // worsen the situation in some of the "regular" cases 
      // (not making the correction means that the conditions for
      // forming a cluster are different if we have missing layers -
      // this could affect events at the boundaries ) 

      // to be chained, two hits need also to be "close" in phi and eta
      bool phiRequirementOK = std::abs(reco::deltaPhi(phi_new[jRH_new],phi_old[jRH_old])) < dPhiChainBoxMax;
      bool etaRequirementOK = fabs(eta_new[jRH_new]-eta_old[jRH_old]) < dEtaChainBoxMax;
      // and the difference in layer index should be < (nlayers-1)
      bool layerRequirementOK = abs(layer_new[jRH_new]-layer_old[jRH_old]) < (theEnsemble.first->id().nlayers()-1);
      // and they should have a time difference compatible with the hypothesis 
      // that the rechits originate from the same particle, but were detected in different layers
      bool timeRequirementOK = fabs(time_new[jRH_new] - time_old[jRH_old]) < dTimeChainBoxMax;

      if(layerRequirementOK && phiRequirementOK && etaRequirementOK && timeRequirementOK){
        return true;
      }
    }
  } 
  return false;
}

std::vector< ME0Segment > ME0SegAlgoMM::run	(	const ME0Ensemble &	ensemble,
		const EnsembleHitContainer &	rechits
	)

Build segments for all desired groups of hits

Definition at line 51 of file ME0SegAlgoMM.cc.

References buildSegments(), chainHits(), clusterHits(), preClustering, preClustering_useChaining, and theEnsemble.

                                                                                                        {

  theEnsemble = ensemble;

  #ifdef EDM_ML_DEBUG // have lines below only compiled when in debug mode
  ME0DetId chId((theEnsamble.first)->id());  
  edm::LogVerbatim("GEMSegAlgoMM") << "[ME0SegAlgoMM::run] build segments in chamber " << chId << " which contains "<<rechits.size()<<" rechits";
  for (auto rh=rechits.begin(); rh!=rechits.end(); ++rh){
    auto me0id = (*rh)->me0Id();
    auto rhLP = (*rh)->localPosition();
    edm::LogVerbatim("ME0SegAlgoMM") << "[RecHit :: Loc x = "<<std::showpos<<std::setw(9)<<rhLP.x()<<" Loc y = "<<std::showpos<<std::setw(9)<<rhLP.y()
                                     <<" Time = "<<std::showpos<<(*rh)->tof()<<" -- "<<me0id.rawId()<<" = "<<me0id<<" ]";
  }
  #endif

  // pre-cluster rechits and loop over all sub clusters separately
  std::vector<ME0Segment>          segments_temp;
  std::vector<ME0Segment>          segments;
  ProtoSegments rechits_clusters; // this is a collection of groups of rechits
  
  if(preClustering) {
    // run a pre-clusterer on the given rechits to split obviously separated segment seeds:
    if(preClustering_useChaining){
      // it uses X,Y,Z information; there are no configurable parameters used;
      // the X, Y, Z "cuts" are just (much) wider than reasonable high pt segments
      rechits_clusters = this->chainHits( rechits );
    }
    else{
      // it uses X,Y information + configurable parameters
      rechits_clusters = this->clusterHits(rechits );
    }
    // loop over the found clusters:
    for(auto sub_rechits = rechits_clusters.begin(); sub_rechits !=  rechits_clusters.end(); ++sub_rechits ) {
      // clear the buffer for the subset of segments:
      segments_temp.clear();
      // build the subset of segments:
      segments_temp = this->buildSegments( (*sub_rechits) );
      // add the found subset of segments to the collection of all segments in this chamber:
      segments.insert( segments.end(), segments_temp.begin(), segments_temp.end() );
    }
  

    return segments;
  }
  else {
    segments = this->buildSegments(rechits);
    return segments;
  }
}

Member Data Documentation

bool ME0SegAlgoMM::debug

private

Definition at line 63 of file ME0SegAlgoMM.h.

Referenced by ME0SegAlgoMM().

double ME0SegAlgoMM::dEtaChainBoxMax

private

Definition at line 70 of file ME0SegAlgoMM.h.

Referenced by isGoodToMerge(), and ME0SegAlgoMM().

double ME0SegAlgoMM::dPhiChainBoxMax

private

Definition at line 69 of file ME0SegAlgoMM.h.

Referenced by isGoodToMerge(), and ME0SegAlgoMM().

double ME0SegAlgoMM::dTimeChainBoxMax

private

Definition at line 71 of file ME0SegAlgoMM.h.

Referenced by isGoodToMerge(), and ME0SegAlgoMM().

double ME0SegAlgoMM::dXclusBoxMax

private

Definition at line 66 of file ME0SegAlgoMM.h.

Referenced by clusterHits(), and ME0SegAlgoMM().

double ME0SegAlgoMM::dYclusBoxMax

private

Definition at line 67 of file ME0SegAlgoMM.h.

Referenced by clusterHits(), and ME0SegAlgoMM().

int ME0SegAlgoMM::maxRecHitsInCluster

private

Definition at line 72 of file ME0SegAlgoMM.h.

Referenced by ME0SegAlgoMM().

unsigned int ME0SegAlgoMM::minHitsPerSegment

private

Definition at line 64 of file ME0SegAlgoMM.h.

Referenced by buildSegments(), and ME0SegAlgoMM().

const std::string ME0SegAlgoMM::myName

private

Definition at line 59 of file ME0SegAlgoMM.h.

bool ME0SegAlgoMM::preClustering

private

Definition at line 65 of file ME0SegAlgoMM.h.

Referenced by ME0SegAlgoMM(), and run().

bool ME0SegAlgoMM::preClustering_useChaining

private

Definition at line 68 of file ME0SegAlgoMM.h.

Referenced by ME0SegAlgoMM(), and run().

EnsembleHitContainer ME0SegAlgoMM::proto_segment

private

Definition at line 75 of file ME0SegAlgoMM.h.

Referenced by buildSegments().

float ME0SegAlgoMM::running_max =999999.

staticprivate

Definition at line 78 of file ME0SegAlgoMM.h.

Referenced by clusterHits().

std::unique_ptr<ME0SegFit> ME0SegAlgoMM::sfit_

private

Definition at line 79 of file ME0SegAlgoMM.h.

Referenced by buildSegments().

ME0Ensemble ME0SegAlgoMM::theEnsemble

private

Definition at line 76 of file ME0SegAlgoMM.h.

Referenced by buildSegments(), isGoodToMerge(), and run().