ParticleFlowReco/src/PFBlock.cc

0001 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0004
0005 #include <iomanip>
0006 #include <sstream>
0007
0008 using namespace std;
0009 using namespace reco;
0010
0011 void PFBlock::addElement(PFBlockElement* element) {
0012   element->setIndex(elements_.size());
0013   element->lock();
0014   elements_.push_back(element->clone());
0015 }
0016
0017 void PFBlock::bookLinkData() {}
0018
0019 void PFBlock::setLink(unsigned i1, unsigned i2, double Dist, LinkData& linkData, LinkTest test) const {
0020   assert(test < LINKTEST_ALL);
0021
0022   unsigned index = 0;
0023   bool ok = matrix2vector(i1, i2, index);
0024
0025   if (ok) {
0026     //ignore the  -1, -1 pair
0027     if (Dist > -0.5) {
0028       Link& l = linkData[index];
0029       l.distance = Dist;
0030       l.test |= (1 << test);
0031     } else  //delete if existing
0032     {
0033       LinkData::iterator it = linkData.find(index);
0034       if (it != linkData.end())
0035         linkData.erase(it);
0036     }
0037
0038   } else {
0039     assert(0);
0040   }
0041 }
0042
0043 // void PFBlock::lock(unsigned i, LinkData& linkData ) const {
0044
0045 //   assert( linkData.size() == linkDataSize() );
0046
0047 //   for(unsigned j=0; j<elements_.size(); j++) {
0048
0049 //     if(i==j) continue;
0050
0051 //     unsigned index = 0;
0052 //     bool ok =  matrix2vector(i,j, index);
0053 //     if(ok)
0054 //       linkData[index] = -1;
0055 //     else
0056 //       assert(0);
0057 //   }
0058 // }
0059
0060 void PFBlock::associatedElements(unsigned i,
0061                                  const LinkData& linkData,
0062                                  multimap<double, unsigned>& sortedAssociates,
0063                                  PFBlockElement::Type type,
0064                                  LinkTest test) const {
0065   sortedAssociates.clear();
0066
0067   // i is too large
0068   if (i > elements_.size())
0069     return;
0070   // assert(i>=0); i >= 0, since i is unsigned
0071
0072   for (unsigned ie = 0; ie < elements_.size(); ie++) {
0073     // considered element itself
0074     if (ie == i) {
0075       continue;
0076     }
0077     // not the right type
0078     if (type != PFBlockElement::NONE && elements_[ie].type() != type) {
0079       continue;
0080     }
0081
0082     // Order the elements by increasing distance !
0083
0084     unsigned index = 0;
0085     if (!matrix2vector(i, ie, index))
0086       continue;
0087
0088     double c2 = -1;
0089     LinkData::const_iterator it = linkData.find(index);
0090     if (it != linkData.end() && (((1 << test) & it->second.test) != 0 || (test == LINKTEST_ALL)))
0091       c2 = it->second.distance;
0092
0093     // not associated
0094     if (c2 < 0) {
0095       continue;
0096     }
0097
0098     sortedAssociates.insert(pair<double, unsigned>(c2, ie));
0099   }
0100 }
0101
0102 bool PFBlock::matrix2vector(unsigned iindex, unsigned jindex, unsigned& index) const {
0103   unsigned size = elements_.size();
0104   if (iindex == jindex || iindex >= size || jindex >= size) {
0105     return false;
0106   }
0107
0108   if (iindex > jindex)
0109     std::swap(iindex, jindex);
0110
0111   index = jindex - iindex - 1;
0112
0113   if (iindex > 0) {
0114     index += iindex * size;
0115     unsigned missing = iindex * (iindex + 1) / 2;
0116     index -= missing;
0117   }
0118
0119   return true;
0120 }
0121
0122 double PFBlock::dist(unsigned ie1, unsigned ie2, const LinkData& linkData) const {
0123   double Dist = -1;
0124
0125   unsigned index = 0;
0126   if (!matrix2vector(ie1, ie2, index))
0127     return Dist;
0128   LinkData::const_iterator it = linkData.find(index);
0129   if (it != linkData.end())
0130     Dist = it->second.distance;
0131
0132   return Dist;
0133 }
0134
0135 ostream& reco::operator<<(ostream& out, const reco::PFBlock& block) {
0136   if (!out)
0137     return out;
0138   const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
0139   out << "\t--- PFBlock ---  " << endl;
0140   out << "\tnumber of elements: " << elements.size() << endl;
0141
0142   // Build element label (string) : elid from type, layer and occurence number
0143   // use stringstream instead of sprintf to concatenate string and integer into string
0144
0145   vector<string> elid;
0146   string s;
0147   stringstream ss;
0148   int iel = 0;
0149   int iTK = 0;
0150   int iGSF = 0;
0151   int iBREM = 0;
0152   int iPS1 = 0;
0153   int iPS2 = 0;
0154   int iEE = 0;
0155   int iEB = 0;
0156   int iHE = 0;
0157   int iHB = 0;
0158   int iHFEM = 0;
0159   int iHFHAD = 0;
0160   int iSC = 0;
0161   int iHO = 0;
0162
0163   // for each element in turn
0164   std::vector<bool> toPrint(elements.size(), static_cast<bool>(true));
0165   for (unsigned ie = 0; ie < elements.size(); ie++) {
0166     PFBlockElement::Type type = elements[ie].type();
0167     std::multimap<double, unsigned> ecalElems;
0168     switch (type) {
0169       case PFBlockElement::TRACK:
0170         iTK++;
0171         ss << "TK" << iTK;
0172         break;
0173       case PFBlockElement::GSF:
0174         iGSF++;
0175         ss << "GSF" << iGSF;
0176         break;
0177       case PFBlockElement::BREM:
0178         block.associatedElements(
0179             ie, block.linkData(), ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
0180         iBREM++;
0181         if (!ecalElems.empty()) {
0182           ss << "BR" << iBREM;
0183         } else {
0184           toPrint[ie] = false;
0185         }
0186         break;
0187       case PFBlockElement::SC:
0188         iSC++;
0189         ss << "SC" << iSC;
0190         break;
0191       default: {
0192         PFClusterRef clusterref = elements[ie].clusterRef();
0193         int layer = clusterref->layer();
0194         switch (layer) {
0195           case PFLayer::PS1:
0196             iPS1++;
0197             ss << "PV" << iPS1;
0198             break;
0199           case PFLayer::PS2:
0200             iPS2++;
0201             ss << "PH" << iPS2;
0202             break;
0203           case PFLayer::ECAL_ENDCAP:
0204             iEE++;
0205             ss << "EE" << iEE;
0206             break;
0207           case PFLayer::ECAL_BARREL:
0208             iEB++;
0209             ss << "EB" << iEB;
0210             break;
0211           case PFLayer::HCAL_ENDCAP:
0212             iHE++;
0213             ss << "HE" << iHE;
0214             break;
0215           case PFLayer::HCAL_BARREL1:
0216             iHB++;
0217             ss << "HB" << iHB;
0218             break;
0219           case PFLayer::HCAL_BARREL2:
0220             iHO++;
0221             ss << "HO" << iHO;
0222             break;
0223           case PFLayer::HF_EM:
0224             iHFEM++;
0225             ss << "FE" << iHFEM;
0226             break;
0227           case PFLayer::HF_HAD:
0228             iHFHAD++;
0229             ss << "FH" << iHFHAD;
0230             break;
0231           default:
0232             iel++;
0233             ss << "??" << iel;
0234             break;
0235         }
0236         break;
0237       }
0238     }
0239     s = ss.str();
0240     elid.push_back(s);
0241     // clear stringstream
0242     ss.str("");
0243
0244     if (toPrint[ie])
0245       out << "\t" << s << " " << elements[ie] << endl;
0246   }
0247
0248   out << endl;
0249
0250   int width = 6;
0251   if (!block.linkData().empty()) {
0252     out << endl << "\tlink data (distance x 1000): " << endl;
0253     out << setiosflags(ios::right);
0254     out << "\t" << setw(width) << " ";
0255     for (unsigned ie = 0; ie < elid.size(); ie++)
0256       if (toPrint[ie])
0257         out << setw(width) << elid[ie];
0258     out << endl;
0259     out << setiosflags(ios::fixed);
0260     out << setprecision(1);
0261
0262     for (unsigned i = 0; i < block.elements().size(); i++) {
0263       if (!toPrint[i])
0264         continue;
0265       out << "\t";
0266       out << setw(width) << elid[i];
0267       for (unsigned j = 0; j < block.elements().size(); j++) {
0268         if (!toPrint[j])
0269           continue;
0270         double Dist = block.dist(i, j, block.linkData());  //,PFBlock::LINKTEST_ALL);
0271
0272         // out<<setw(width)<< Dist*1000.;
0273         if (Dist > -0.5)
0274           out << setw(width) << Dist * 1000.;
0275         else
0276           out << setw(width) << " ";
0277       }
0278       out << endl;
0279     }
0280
0281     out << endl << "\tlink data (distance x 1000) for tracking links : " << endl;
0282     out << setiosflags(ios::right);
0283     out << "\t" << setw(width) << " ";
0284     for (unsigned ie = 0; ie < elid.size(); ie++)
0285       if (toPrint[ie] &&
0286           (block.elements()[ie].type() == PFBlockElement::TRACK || block.elements()[ie].type() == PFBlockElement::GSF))
0287         out << setw(width) << elid[ie];
0288     out << endl;
0289     out << setiosflags(ios::fixed);
0290     out << setprecision(1);
0291
0292     for (unsigned i = 0; i < block.elements().size(); i++) {
0293       if (!toPrint[i] ||
0294           (block.elements()[i].type() != PFBlockElement::TRACK && block.elements()[i].type() != PFBlockElement::GSF))
0295         continue;
0296       out << "\t";
0297       out << setw(width) << elid[i];
0298       for (unsigned j = 0; j < block.elements().size(); j++) {
0299         if (!toPrint[j] ||
0300             (block.elements()[j].type() != PFBlockElement::TRACK && block.elements()[j].type() != PFBlockElement::GSF))
0301           continue;
0302         double Dist = block.dist(i, j, block.linkData());  //,PFBlock::LINKTEST_ALL);
0303
0304         // out<<setw(width)<< Dist*1000.;
0305         if (Dist > -0.5)
0306           out << setw(width) << Dist * 1000.;
0307         else
0308           out << setw(width) << " ";
0309       }
0310       out << endl;
0311     }
0312
0313     out << setprecision(3);
0314     out << resetiosflags(ios::right | ios::fixed);
0315
0316   } else {
0317     out << "\tno links." << endl;
0318   }
0319
0320   return out;
0321 }
0322
0323 unsigned PFBlock::linkDataSize() const {
0324   unsigned n = elements_.size();
0325
0326   // number of possible undirected links between n elements.
0327   // reflective links impossible.
0328
0329   return n * (n - 1) / 2;
0330 }