PFBlock.cc

Go to the documentation of this file.

#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFLayer.h"

#include <iomanip>
#include <sstream>

using namespace std;
using namespace reco;


void PFBlock::addElement(  PFBlockElement* element ) {
  element->setIndex( elements_.size() );
  element->lock();
  elements_.push_back( element->clone() ); 

}


void PFBlock::bookLinkData() {

}



void PFBlock::setLink(unsigned i1, 
              unsigned i2, 
              double Dist,
                      LinkData& linkData, 
              LinkTest test) const {
  
  assert( test<LINKTEST_ALL );
  
  unsigned index = 0;
  bool ok =  matrix2vector(i1,i2, index);

  if(ok) {
    //ignore the  -1, -1 pair
    if ( Dist > -0.5 ) {
      Link & l = linkData[index];
      l.distance=Dist;
      l.test |= (1 << test);
    }     else  //delete if existing
      {
    LinkData::iterator it = linkData.find(index);
    if(it!=linkData.end()) linkData.erase(it);
      }

  } else {
    assert(0);
  }
  
}


// void PFBlock::lock(unsigned i, LinkData& linkData ) const {
  
//   assert( linkData.size() == linkDataSize() );
  
//   for(unsigned j=0; j<elements_.size(); j++) {
    
//     if(i==j) continue;
    
//     unsigned index = 0;
//     bool ok =  matrix2vector(i,j, index);
//     if(ok)
//       linkData[index] = -1;
//     else 
//       assert(0);
//   }
// }



void PFBlock::associatedElements( unsigned i, 
                                  const LinkData& linkData, 
                                  multimap<double, unsigned>& sortedAssociates,
                                  PFBlockElement::Type type,
                  LinkTest 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;
    }
    // not the right type
    if(type !=  PFBlockElement::NONE && 
       elements_[ie].type() != type ) {
      continue;
    }

    // Order the elements by increasing distance !

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

    double c2=-1;
    LinkData::const_iterator it =  linkData.find(index);
    if ( it!=linkData.end() && 
     ( ( (1 << test ) & it->second.test) !=0 || (test == LINKTEST_ALL) ) ) 
      c2= it->second.distance;

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

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

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

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

  
  index = jindex-iindex-1;

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


double PFBlock::dist( unsigned ie1, unsigned ie2,
                      const LinkData& linkData) const {

  double Dist = -1;

  unsigned index = 0;
  if( !matrix2vector(ie1, ie2, index) ) return Dist;
  LinkData::const_iterator it =  linkData.find(index);
  if( it!=linkData.end() ) Dist= it->second.distance;

  return Dist;

}


ostream& reco::operator<<(  ostream& out, 
                            const reco::PFBlock& block ) {

  if(! out) return out;
  const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
  out<<"\t--- PFBlock ---  "<<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
 
  vector <string> elid;
  string s;
  stringstream ss;
  int iel = 0;
  int iTK =0;
  int iGSF =0;
  int iBREM=0;
  int iPS1 = 0;
  int iPS2 = 0;
  int iEE = 0;
  int iEB = 0;
  int iHE = 0;
  int iHB = 0;
  int iHFEM = 0;
  int iHFHAD = 0;
  int iSC = 0;
  int iHO = 0;

  // for each element in turn
  std::vector<bool> toPrint(elements.size(),static_cast<bool>(true));
  for(unsigned ie=0; ie<elements.size(); ie++) {
    
    PFBlockElement::Type type = elements[ie].type();
    std::multimap<double, unsigned> ecalElems;
    switch(type){
    case PFBlockElement::TRACK:
      iTK++;
      ss << "TK" << iTK;
      break;
    case PFBlockElement::GSF:
      iGSF++;
      ss << "GSF" << iGSF;
      break;
    case PFBlockElement::BREM:
      block.associatedElements( ie,  block.linkData(),
                ecalElems ,
                reco::PFBlockElement::ECAL,
                reco::PFBlock::LINKTEST_ALL );
      iBREM++;
      if ( !ecalElems.empty() ) { 
    ss << "BR" << iBREM;
      } else {
    toPrint[ie] = false;
      }
      break;    
    case PFBlockElement::SC:
      iSC++;
      ss << "SC" << iSC;
      break;
    default:{
      PFClusterRef clusterref = elements[ie].clusterRef();
      int layer = clusterref->layer();
      switch (layer){
      case PFLayer::PS1:
        iPS1++;
        ss << "PV" << iPS1;
        break;
      case PFLayer::PS2:
        iPS2++;
        ss << "PH" << iPS2;
        break;
      case PFLayer::ECAL_ENDCAP:
        iEE++;
        ss << "EE" << iEE;
        break;
      case PFLayer::ECAL_BARREL:
        iEB++;
        ss << "EB" << iEB;
        break;
      case PFLayer::HCAL_ENDCAP:
        iHE++;
        ss << "HE" << iHE;
        break;
      case PFLayer::HCAL_BARREL1:
        iHB++;
        ss << "HB" << iHB;
        break;
      case PFLayer::HCAL_BARREL2:
        iHO++;
        ss << "HO" << iHO;
        break;
      case PFLayer::HF_EM:
        iHFEM++;
        ss << "FE" << iHFEM;
        break;
      case PFLayer::HF_HAD:
        iHFHAD++;
        ss << "FH" << iHFHAD;
    break;
      default:
        iel++;   
        ss << "??" << iel;
        break;
      }
      break;
    }
    }
    s = ss.str();
    elid.push_back( s );
    // clear stringstream
    ss.str("");

    if ( toPrint[ie] ) out<<"\t"<< s <<" "<<elements[ie] <<endl;
  }
  
  out<<endl;

  int width = 6;
  if( !block.linkData().empty() ) {
    out<<endl<<"\tlink data (distance x 1000): "<<endl;
    out<<setiosflags(ios::right);
    out<<"\t" << setw(width) << " ";
    for(unsigned ie=0; ie<elid.size(); ie++) 
      if ( toPrint[ie] ) out <<setw(width)<< elid[ie];
    out<<endl;  
    out<<setiosflags(ios::fixed);
    out<<setprecision(1);      
  
    for(unsigned i=0; i<block.elements().size(); i++) {
      if ( !toPrint[i] ) continue;
      out<<"\t";
      out <<setw(width) << elid[i];
      for(unsigned j=0; j<block.elements().size(); j++) {
    if ( !toPrint[j] ) continue;
        double Dist = block.dist(i,j, block.linkData());//,PFBlock::LINKTEST_ALL);

    // out<<setw(width)<< Dist*1000.;
    if (Dist > -0.5) out<<setw(width)<< Dist*1000.; 
        else  out <<setw(width)<< " ";
      }
      out<<endl;
    }

        
    out<<endl<<"\tlink data (distance x 1000) for tracking links : "<<endl;
    out<<setiosflags(ios::right);
    out<<"\t" << setw(width) << " ";
    for(unsigned ie=0; ie<elid.size(); ie++) 
      if ( toPrint[ie] && 
       ( block.elements()[ie].type() == PFBlockElement::TRACK ||
         block.elements()[ie].type() == PFBlockElement::GSF )) 
    out <<setw(width)<< elid[ie];
    out<<endl;  
    out<<setiosflags(ios::fixed);
    out<<setprecision(1);      
  
    for(unsigned i=0; i<block.elements().size(); i++) {
      if ( !toPrint[i] || 
       (block.elements()[i].type() != PFBlockElement::TRACK &&
        block.elements()[i].type() != PFBlockElement::GSF )) continue;
      out<<"\t";
      out <<setw(width) << elid[i];
      for(unsigned j=0; j<block.elements().size(); j++) {
    if ( !toPrint[j] || 
         (block.elements()[j].type() != PFBlockElement::TRACK &&
          block.elements()[j].type() != PFBlockElement::GSF )) continue;
    double Dist = block.dist(i,j, block.linkData());//,PFBlock::LINKTEST_ALL);

    // out<<setw(width)<< Dist*1000.;
    if (Dist > -0.5) out<<setw(width)<< Dist*1000.; 
    else  out <<setw(width)<< " ";
      }
      out<<endl;
    }
    
    out<<setprecision(3);  
    out<<resetiosflags(ios::right|ios::fixed);

  }
  else {
    out<<"\tno links."<<endl;
  }
  
  return out;
}
 
 
unsigned PFBlock::linkDataSize() const {
  unsigned n = elements_.size();
  
  // number of possible undirected links between n elements.
  // reflective links impossible.
 
  return n*(n-1)/2; 
}