CMSSW_4_4_3_patch1/src/DataFormats/ParticleFlowReco/src/PFDisplacedVertexCandidate.cc

Go to the documentation of this file.

#include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertexCandidate.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/Math/interface/Point3D.h"

using namespace std;
using namespace reco;


PFDisplacedVertexCandidate::PFDisplacedVertexCandidate(){}

void PFDisplacedVertexCandidate::addElement(const TrackBaseRef element) {
  elements_.push_back( element ); 
}


void PFDisplacedVertexCandidate::setLink(const unsigned i1, 
                                         const unsigned i2, 
                                         const float dist,
                                         const GlobalPoint& dcaPoint,
                                         const VertexLinkTest test){

  
  assert( test<LINKTEST_ALL );
  
  unsigned index = 0;
  bool ok =  matrix2vector(i1,i2, index);

  if(ok) {
    //ignore the  -1, -1 pair
    if ( dist > -0.5 ) {
      VertexLink & l = vertexLinkData_[index];
      l.distance_ = dist;
      l.dcaPoint_ = dcaPoint;
      l.test_ |= (1 << test);
    }     else  //delete if existing
      {
        VertexLinkData::iterator it = vertexLinkData_.find(index);
        if(it!=vertexLinkData_.end()) vertexLinkData_.erase(it);
      }

  } else {
    assert(0);
  }
  
}


void PFDisplacedVertexCandidate::associatedElements( const unsigned i, 
                                                     const VertexLinkData& vertexLinkData, 
                                                     multimap<float, unsigned>& sortedAssociates,
                                                     const VertexLinkTest test ) const {

  sortedAssociates.clear();
  
  // i is too large
  if( i > elements_.size() ) return;
  // assert(i>=0); // i >= 0, since i is unsigned
  for(unsigned ie=0; ie<elements_.size(); ie++) {
    
    // considered element itself
    if( ie == i ) {
      continue;
    }

    // Order the elements by increasing distance !

    unsigned index = 0;
    if( !matrix2vector(i, ie, index) ) continue;

    float c2=-1;
    VertexLinkData::const_iterator it =  vertexLinkData.find(index);
    if ( it!=vertexLinkData.end() && 
         ( ( (1 << test ) & it->second.test_) !=0 || (test == LINKTEST_ALL) ) ) 
      c2= it->second.distance_;

    // not associated
    if( c2 < 0 ) { 
      continue;
    }

    sortedAssociates.insert( pair<float,unsigned>(c2, ie) );
  }
} 







// -------- Provide useful information -------- //


PFDisplacedVertexCandidate::DistMap PFDisplacedVertexCandidate::r2Map() const {

  DistMap r2Map;

  for (unsigned ie1 = 0; ie1<elements_.size(); ie1++)
    for (unsigned ie2 = ie1+1; ie2<elements_.size(); ie2++){

      GlobalPoint P = dcaPoint(ie1, ie2);
      if (P.x() > 1e9) continue;

      float r2 = P.x()*P.x()+P.y()*P.y()+P.z()*P.z();

      r2Map.insert(pair<float, pair<int,int> >(r2, pair <int, int>(ie1, ie2)));
    }

  return r2Map;

}



PFDisplacedVertexCandidate::DistVector PFDisplacedVertexCandidate::r2Vector() const {

  DistVector r2Vector;

  for (unsigned ie1 = 0; ie1<elements_.size(); ie1++)
    for (unsigned ie2 = ie1+1; ie2<elements_.size(); ie2++){

      GlobalPoint P = dcaPoint(ie1, ie2);
      if (P.x() > 1e9) continue;

      float r2 = P.x()*P.x()+P.y()*P.y()+P.z()*P.z();

      r2Vector.push_back(r2);
    }

  return r2Vector;

}


PFDisplacedVertexCandidate::DistVector PFDisplacedVertexCandidate::distVector() const {

  DistVector distVector;


  for (unsigned ie1 = 0; ie1<elements_.size(); ie1++)
    for (unsigned ie2 = ie1+1; ie2<elements_.size(); ie2++){

      float d = dist(ie1, ie2);
      if (d < -0.5) continue;

      distVector.push_back(d);

    }

  return distVector;

}

const GlobalPoint PFDisplacedVertexCandidate::dcaPoint( unsigned ie1, unsigned ie2) const {

  GlobalPoint dcaPoint(1e10,1e10,1e10);

  unsigned index = 0;
  if( !matrix2vector(ie1, ie2, index) ) return dcaPoint;
  VertexLinkData::const_iterator it =  vertexLinkData_.find(index);
  if( it!=vertexLinkData_.end() ) dcaPoint = it->second.dcaPoint_;

  return dcaPoint;

}


// -------- Internal tools -------- //

bool PFDisplacedVertexCandidate::testLink(unsigned ie1, unsigned ie2) const {
  float d = dist( ie1, ie2);
  if (d < -0.5) return false;
  return true;
}


const float PFDisplacedVertexCandidate::dist( unsigned ie1, unsigned ie2) const {

  float dist = -1;

  unsigned index = 0;
  if( !matrix2vector(ie1, ie2, index) ) return dist;
  VertexLinkData::const_iterator it =  vertexLinkData_.find(index);
  if( it!=vertexLinkData_.end() ) dist= it->second.distance_;

  return dist;

}









// -------- Storage of the information -------- //


unsigned PFDisplacedVertexCandidate::vertexLinkDataSize() const {
  unsigned n = elements_.size();
  
  // number of possible undirected links between n elements.
  // reflective links impossible.
 
  return n*(n-1)/2; 
}


bool PFDisplacedVertexCandidate::matrix2vector( unsigned iindex, 
                                                unsigned jindex, 
                                                unsigned& index ) const {

  unsigned size = elements_.size();
  if( iindex == jindex || 
      iindex >=  size ||
      jindex >=  size ) {
    return false;
  }
  
  if( iindex > jindex ) 
    swap( iindex, jindex);

  
  index = jindex-iindex-1;

  if(iindex>0) {
    index += iindex*size;
    unsigned missing = iindex*(iindex+1)/2;
    index -= missing;
  }
  
  return true;
}

void PFDisplacedVertexCandidate::Dump( ostream& out ) const {
  if(! out ) return;

  const vector < TrackBaseRef >& elements = elements_;
  out<<"\t--- DisplacedVertexCandidate ---  "<<endl;
  out<<"\tnumber of elements: "<<elements.size()<<endl;
  
  // Build element label (string) : elid from type, layer and occurence number
  // use stringstream instead of sprintf to concatenate string and integer into string
  for(unsigned ie=0; ie<elements.size(); ie++) {

    math::XYZPoint Pi(elements[ie].get()->innerPosition());
    math::XYZPoint Po(elements[ie].get()->outerPosition());

    float innermost_radius = sqrt(Pi.x()*Pi.x() + Pi.y()*Pi.y() + Pi.z()*Pi.z());
    float outermost_radius = sqrt(Po.x()*Po.x() + Po.y()*Po.y() + Po.z()*Po.z());
    float innermost_rho = sqrt(Pi.x()*Pi.x() + Pi.y()*Pi.y());
    float outermost_rho = sqrt(Po.x()*Po.x() + Po.y()*Po.y());
    
    double pt = elements[ie]->pt();


    out<<"ie = " << elements[ie].key() << " pt = " << pt
       <<" innermost hit radius = " << innermost_radius << " rho = " << innermost_rho
       <<" outermost hit radius = " << outermost_radius << " rho = " << outermost_rho
       <<endl;
  }
   
  out<<endl;


}