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;
}