Index: trunk/fftjet/RecombinedJet.icc =================================================================== --- trunk/fftjet/RecombinedJet.icc (revision 44) +++ trunk/fftjet/RecombinedJet.icc (revision 45) @@ -1,58 +1,60 @@ #include namespace fftjet { template inline double RecombinedJet::magnitude() const { const double x(vec_.px()); const double y(vec_.py()); return sqrt(x*x + y*y); } template inline bool RecombinedJet::operator==( const RecombinedJet& r) const { return peak_ == r.peak_ && vec_ == r.vec_ && ncells_ == r.ncells_ && etSum_ == r.etSum_ && centroidEta_ == r.centroidEta_ && centroidPhi_ == r.centroidPhi_ && etaWidth_ == r.etaWidth_ && phiWidth_ == r.phiWidth_ && etaPhiCorr_ == r.etaPhiCorr_ && fuzziness_ == r.fuzziness_ && + isolation_ == r.isolation_ && + entropy_ == r.entropy_ && pileup_ == r.pileup_ && uncertainty_ == r.uncertainty_ && convergenceDistance_ == r.convergenceDistance_; } template inline bool RecombinedJet::operator!=( const RecombinedJet& r) const { return !(*this == r); } template inline bool RecombinedJet::operator<( const RecombinedJet& r) const { return magnitude() < r.magnitude(); } template inline bool RecombinedJet::operator>( const RecombinedJet& r) const { return magnitude() > r.magnitude(); } template inline RecombinedJet RecombinedJet::dummy() { return RecombinedJet(Peak::dummy(), VectorLike(), - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0); } } Index: trunk/fftjet/RecombinedJet.hh =================================================================== --- trunk/fftjet/RecombinedJet.hh (revision 44) +++ trunk/fftjet/RecombinedJet.hh (revision 45) @@ -1,233 +1,249 @@ #ifndef FFTJET_RECOMBINEDJET_HH_ #define FFTJET_RECOMBINEDJET_HH_ //========================================================================= // RecombinedJet.hh // // This class stores the jet information produced by the energy // recombination algorithms. It includes the clustered quantity // (e.g., a 4-vector) together with the original precluster from // the pattern recognition stage, as well as several other // quantities. This information can be retrieved using the // following methods: // // precluster() -- Returns the precluster. // // vec() -- Returns the result of the recombination // (normally, a 4-vector). // // ncells() -- The weighted number of energy discretization grid // cells contributing to this jet. Depending on the // algorithm settings (in particular, the "dataCutoff" // parameter of KernelRecombinationAlg and similar // classes), this number may or may not coincide with // the jet area. // // etSum() -- The weighted sum of Et (pt, etc.) values // contributing to this jet. // // centroidEta() -- (Pseudo)rapidity of the weighted Et centroid. // // centroidPhi() -- Azimuthal angle of the weighted Et centroid. // // etaWidth() -- Weighted eta RMS of the Et cells contributing // to this jet. // // phiWidth() -- Weighted phi RMS of the Et cells contributing // to this jet. // // etaPhiCorr() -- The weighted correlation coefficient between // eta and phi. // // fuzziness() -- This quantity characterizes how far away the // jet is from all other jets. In the "fuzzy" mode // this is the fractional Et error due to the // uncertainty in assigning the grid cells to this // jet, calculated according to the generalized // binomial distribution model. In the "crisp" mode // this quantity does not have a well-defined meaning, // but it is also a dimensionless number which // becomes close to 0 for well-separated jets. // +// isolation() -- This quantity is close to 1 for well-isolated +// jets with small constribution from background. +// +// entropy() -- Rougly, the logarithm of the effective number +// of constituents. +// // All FFTJet algorithms produce this information. User-developed // algorithms are encouraged to follow the same convention. // // In addition, there are two user-settable quantities with obvious // suggested meaning which can be retrieved using methods "pileup()" // and "uncertainty()". These quantities are not calculated by FFTJet // algorithms. Another user-settable quantity, convergenceDistance(), // could be used to indicate jet position and momentum convergence // if the jet reconstruction is performed iteratively. // // Note that, in this package, "etaWidth" and "phiWidth" are defined // with respect to the (eta, phi) centroid, so the width is the smallest // possible. If you want the width with respect to the jet direction, // add in quadrature the angular distance from the jet direction to the // direction of the centroid. // // I. Volobouev // March 2009 //========================================================================= #include "fftjet/Peak.hh" namespace fftjet { template class RecombinedJet { public: inline RecombinedJet(const Peak& peak, const VectorLike& p4, double ncells, double etSum, double centroidEta, double centroidPhi, double etaWidth, double phiWidth, double etaPhiCorr, double fuzziness, + double isolation, double entropy, double pileup=0.0, double uncertainty=0.0, double convergenceDistance=-1.0) : peak_(peak), vec_(p4), ncells_(ncells), etSum_(etSum), centroidEta_(centroidEta), centroidPhi_(centroidPhi), etaWidth_(etaWidth), phiWidth_(phiWidth), etaPhiCorr_(etaPhiCorr), - fuzziness_(fuzziness), pileup_(pileup), + fuzziness_(fuzziness), isolation_(isolation), + entropy_(entropy), pileup_(pileup), uncertainty_(uncertainty), convergenceDistance_(convergenceDistance) {} // The following functions simply return the arguments // provided in the class constructor inline const Peak& precluster() const {return peak_;} inline const VectorLike& vec() const {return vec_;} inline double ncells() const {return ncells_;} inline double etSum() const {return etSum_;} inline double centroidEta() const {return centroidEta_;} inline double centroidPhi() const {return centroidPhi_;} inline double etaWidth() const {return etaWidth_;} inline double phiWidth() const {return phiWidth_;} inline double etaPhiCorr() const {return etaPhiCorr_;} inline double fuzziness() const {return fuzziness_;} + inline double isolation() const {return isolation_;} + inline double entropy() const {return entropy_;} // It is possible to change the precluster inline void setPrecluster(const Peak& peak) {peak_ = peak;} // Modifiers for user-defined quantities inline void setPileup(const double d) {pileup_ = d;} inline void setUncertainty(const double d) {uncertainty_ = d;} inline void setConvergenceDistance(const double d) {convergenceDistance_ = d;} // Corresponding accessors inline double pileup() const {return pileup_;} inline double uncertainty() const {return uncertainty_;} inline double convergenceDistance() const {return convergenceDistance_;} // Some useful utilities bool operator==(const RecombinedJet& r) const; bool operator!=(const RecombinedJet& r) const; // We will use the 4-vector pt as the "magnitude" double magnitude() const; // Sorting by the magnitude bool operator<(const RecombinedJet& r) const; bool operator>(const RecombinedJet& r) const; // We need to provide functions which allow this class to serve // together with "Peak" in various templated code. The angular // functions are not provided in this manner: the user might want // to use either rapidity or pseudorapidity as "eta", and we can't // make this choice ahead of time. Nevertheless, the corresponding // peak functions are forwarded with a changed method name. inline double peakEta() const {return peak_.eta();} inline double peakPhi() const {return peak_.phi();} inline double peakMagnitude() const {return peak_.magnitude();} inline double driftSpeed() const {return peak_.driftSpeed();} inline double magSpeed() const {return peak_.magSpeed();} inline double lifetime() const {return peak_.lifetime();} inline double splitTime() const {return peak_.splitTime();} inline double mergeTime() const {return peak_.mergeTime();} inline double scale() const {return peak_.scale();} inline double nearestNeighborDistance() const {return peak_.nearestNeighborDistance();} inline double membershipFactor() const {return peak_.membershipFactor();} inline double recoScale() const {return peak_.recoScale();} inline double recoScaleRatio() const {return peak_.recoScaleRatio();} inline double clusterRadius() const {return peak_.clusterRadius();} inline double clusterSeparation() const {return peak_.clusterSeparation();} inline void hessian(double hessianArray[3]) const {peak_.hessian(hessianArray);} + inline double relativeScaledHDetDeriv() const + {return peak_.relativeScaledHDetDeriv();} inline int code() const {return peak_.code();} inline int status() const {return peak_.status();} inline ScaleSpaceKernel* membershipFunction() const {return peak_.membershipFunction();} // Wrappers for Peak methods related to blob detection inline double hessianDeterminant() const {return peak_.hessianDeterminant();} inline double scaledHDeterminant() const {return peak_.scaledHDeterminant();} inline double relativeScaledHDeterminant() const {return peak_.relativeScaledHDeterminant();} inline double laplacian() const {return peak_.laplacian();} inline double scaledLaplacian(double bwEta, double bwPhi) const {return peak_.scaledLaplacian(bwEta, bwPhi);} inline double relativeScaledLaplacian(double bwEta, double bwPhi) const {return peak_.relativeScaledLaplacian(bwEta, bwPhi);} // Various other modifiers inline void setPeakEtaPhi(const double eta, const double phi) {peak_.setEtaPhi(eta, phi);} inline void setPeakMagnitude(const double d) {peak_.setMagnitude(d);} inline void setDriftSpeed(const double d) {peak_.setDriftSpeed(d);} inline void setMagSpeed(const double d) {peak_.setMagSpeed(d);} inline void setLifetime(const double d) {peak_.setLifetime(d);} inline void setSplitTime(const double d) {peak_.setSplitTime(d);} inline void setMergeTime(const double d) {peak_.setMergeTime(d);} inline void setScale(const double d) {peak_.setScale(d);} inline void setNearestNeighborDistance(const double d) {peak_.setNearestNeighborDistance(d);} inline void setMembershipFactor(const double f) {peak_.setMembershipFactor(f);} inline void setRecoScale(const double s) {peak_.setRecoScale(s);} inline void setRecoScaleRatio(const double d) {peak_.setRecoScaleRatio(d);} inline void setClusterRadius(const double d) {peak_.setClusterRadius(d);} inline void setClusterSeparation(const double d) {peak_.setClusterSeparation(d);} inline void setCode(const int c) {peak_.setCode(c);} inline void setStatus(const int s) {peak_.setStatus(s);} inline void setHessian(const double hessian[3]) {peak_.setHessian(hessian);} + inline void setRelativeScaledHDetDeriv(const double s) + {peak_.setRelativeScaledHDetDeriv(s);} inline void setMembershipFunction(ScaleSpaceKernel* f) {peak_.setMembershipFunction(f);} // Sometimes there is a real need to create a dummy jet... static RecombinedJet dummy(); private: RecombinedJet(); Peak peak_; VectorLike vec_; double ncells_; double etSum_; double centroidEta_; double centroidPhi_; double etaWidth_; double phiWidth_; double etaPhiCorr_; double fuzziness_; + double isolation_; + double entropy_; double pileup_; double uncertainty_; double convergenceDistance_; }; } #include "fftjet/RecombinedJet.icc" #endif // FFTJET_RECOMBINEDJET_HH_ Index: trunk/fftjet/FasterKernelRecombinationAlg.icc =================================================================== --- trunk/fftjet/FasterKernelRecombinationAlg.icc (revision 44) +++ trunk/fftjet/FasterKernelRecombinationAlg.icc (revision 45) @@ -1,347 +1,352 @@ #include +#include namespace fftjet { template FasterKernelRecombinationAlg:: FasterKernelRecombinationAlg( ScaleSpaceKernel* ikernel, const Functor2* bgWeight, const double iunlikelyBgWeight, const double idataCutoff, const bool iwinnerTakesAll, const bool ibuildCorrelationMatrix, const bool ibuildClusterMask, const unsigned ietaBinMin, const unsigned ietaBinMax) : KernelRecombinationAlg::KernelRecombinationAlg( ikernel, bgWeight, iunlikelyBgWeight, idataCutoff, iwinnerTakesAll, ibuildCorrelationMatrix, ibuildClusterMask, ietaBinMin, ietaBinMax), // The following initializer will make sure that ikernel is an AbsKernel2d kernelScans(dynamic_cast(ikernel)), scanTable(0), scanTableLen(0) { } template FasterKernelRecombinationAlg:: ~FasterKernelRecombinationAlg() { delete [] scanTable; } template void FasterKernelRecombinationAlg:: setupTableBuf(const unsigned njets) { if (njets > scanTableLen) { delete [] scanTable; scanTable = new const Real*[njets]; scanTableLen = njets; } } template int FasterKernelRecombinationAlg::run( const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, const unsigned nBgEta, const unsigned nBgPhi, std::vector >* outputJets, VectorLike* unclustered, double* unused) { if (this->performPreliminaryProcessing( peaks, eventData, bgData, nBgEta, nBgPhi, outputJets, unclustered, unused)) return 0; // Make sure individual peak functors are not set for (unsigned ijet=0; ijetmembershipFunction() == 0); // Get memory for the kernel scan pointers setupTableBuf(B::lastNJets); // Set up kernel scan pointers for (unsigned ijet=0; ijetrecombine4Vectors(); // Go over all grid points and figure out to which jet // that point belongs. for (unsigned ieta=B::etaBinMin; ietaeta(); // We will try to get rid of each point in the inner // cycle as fast as possible. The following few lines // of code depend, unfortunately, on the internal layout // of the Grid2d object. On the positive side, we can // get rid of a point below the cutoff without making // the "uncheckedAt(ieta, iphi)" function call. const Real* databuf = eventData.data() + ieta*B::lastNPhi; const BgData* bgbuf; switch (B::bgDim) { case 0: bgbuf = bgData; break; case 1: bgbuf = bgData + ieta; break; case 2: bgbuf = bgData + ieta*B::lastNPhi; break; default: assert(0); } for (unsigned iphi=0; iphivMaker(gridValue, eta, phi)); const double* dPhi = B::dphiBuffer + iphi*B::lastNJets; const double bgW((*B::bgWeightCalc)(gridValue, bgbuf[iphi*bgStride])); // Cycle over jets. This is the most speed-critical // cycle in the whole algorithm. double sum(0.0), maxW(0.0); unsigned maxWJet(B::lastNJets); for (unsigned ijet=0; ijet maxW) { maxW = f; maxWJet = ijet; } } } if (maxWJet == B::lastNJets) { // No contributing cluster, so pure unclustered *unclustered += gridVec; *unused += gridValue; } else { const double totalSum = sum + bgW; if (B::buildClusterMask && maxW > bgW) B::clusterMask[ieta*B::lastNPhi + iphi] = maxWJet + 1; if (B::winnerTakesAll) { // Hard clustering. All energy goes // to the cluster with the highest weight. if (maxW > bgW) { if (use4Vec) B::jets[maxWJet] += gridVec; B::weightSum[maxWJet] += 1.0; B::energySum[maxWJet] += gridValue; + B::entropySum[maxWJet] += gridValue*log(gridValue); + B::doubleWeightedEnergy[maxWJet] += gridValue*maxW/totalSum; B::etaSum[maxWJet] += gridValue*eta; B::etaEnergySum[maxWJet] += gridValue*B::dEta[maxWJet]*B::dEta[maxWJet]; B::phiEnergySum[maxWJet] += gridValue*dPhi[maxWJet]*dPhi[maxWJet]; B::energyVariance[maxWJet] += gridValue*gridValue*(1.0 - maxW/totalSum); B::phiSum[maxWJet] += gridValue*dPhi[maxWJet]; } else { *unclustered += gridVec; *unused += gridValue; } } else { // Fuzzy clustering. Distribute the cell energy // among all clusters according to their weights. // Unclustered energy fraction const double bgfrac = bgW/totalSum; if (bgfrac > 0.0) { *unclustered += gridVec*bgfrac; *unused += gridValue*bgfrac; } // Fractions which belong to each jet for (unsigned ijet=0; ijet 0.0) { const double ePart = gridValue*w; if (use4Vec) B::jets[ijet] += gridVec*w; B::weightSum[ijet] += w; B::energySum[ijet] += ePart; + B::entropySum[ijet] += ePart*log(ePart); + B::doubleWeightedEnergy[ijet] += ePart*w; B::etaEnergySum[ijet] += ePart*B::dEta[ijet]*B::dEta[ijet]; B::etaSum[ijet] += ePart*eta; B::phiEnergySum[ijet] += ePart*dPhi[ijet]*dPhi[ijet]; B::energyVariance[ijet] += gridValue*ePart*(1.0 - w); B::phiSum[ijet] += ePart*dPhi[ijet]; } } } } } } this->buildOutput(outputJets, true); return 0; } template FasterKernelRecombinationAlg:: KernelScanCollection::KernelScanCollection(const AbsKernel2d* ker) : kernel(ker), mappedEta(0), mapBufLen(0) { assert(kernel); } template FasterKernelRecombinationAlg:: KernelScanCollection::~KernelScanCollection() { for (typename std::map*>::iterator it = kernelScans.begin(); it != kernelScans.end(); ++it) delete it->second; delete [] mappedEta; } template void FasterKernelRecombinationAlg:: KernelScanCollection::mapCoords( const Grid2d& eventData, const unsigned etaBinMin, const unsigned etaBinMax, const Peak** peakPositions, const unsigned lastNJets, unsigned **etaBuf, unsigned **phiBuf) { assert(peakPositions); assert(etaBuf); assert(phiBuf); // Allocate enough memory const unsigned lastNEta = eventData.nEta(); const unsigned lastNPhi = eventData.nPhi(); const unsigned mapLen = lastNJets*(lastNPhi + lastNEta); if (mapLen > mapBufLen) { delete [] mappedEta; mappedEta = new unsigned[mapLen]; mapBufLen = mapLen; } unsigned* mappedPhi = mappedEta + lastNJets*lastNEta; // Scan eta const int halfEta = lastNEta/2; for (unsigned ijet=0; ijeteta()); for (int ieta=etaBinMin; ieta < static_cast(etaBinMax); ++ieta) { const int ipos = ieta - ijetEta + halfEta; if (ipos >= 0 && ipos < static_cast(lastNEta)) mappedEta[ieta*lastNJets + ijet] = ipos; else mappedEta[ieta*lastNJets + ijet] = lastNEta; } } // Scan phi const unsigned halfPhi = lastNPhi/2; for (unsigned ijet=0; ijetphi()); for (unsigned iphi=0; iphi const Real* FasterKernelRecombinationAlg:: KernelScanCollection::getKernelScan( const Grid2d& eventData, const double scale) { typename std::map*>::iterator it = kernelScans.find(scale); if (it == kernelScans.end()) { // This should not happen very often. // We need to perform a new kernel scan. Grid2d* scan = new Grid2d(eventData); const unsigned nEta = scan->nEta(); const unsigned nPhi = scan->nPhi(); const Real x0 = scan->etaBinCenter(nEta/2); const Real y0 = scan->phiBinCenter(nPhi/2); long double sum = 0.0L; for (unsigned ieta=0; ietaetaBinCenter(ieta) - x0; for (unsigned iphi=0; iphiphiBinCenter(iphi) - y0; while (dphi > M_PI) dphi -= (2.0*M_PI); while (dphi < -M_PI) dphi += (2.0*M_PI); const double f = (*kernel)(dx, dphi, scale); sum += f; scan->uncheckedSetBin(ieta, iphi, f); } } const double integ = sum*scan->etaBinWidth()*scan->phiBinWidth(); assert(integ); scan->scaleData(1.0/integ); it = kernelScans.insert(std::make_pair(scale,scan)).first; } if (!it->second->isCompatible(eventData)) { // This is not the intended mode of operation for this class assert(!"Please maintain consistency of the data binning!"); } return it->second->data(); } } Index: trunk/fftjet/AbsClusteringTree.icc =================================================================== --- trunk/fftjet/AbsClusteringTree.icc (revision 44) +++ trunk/fftjet/AbsClusteringTree.icc (revision 45) @@ -1,930 +1,947 @@ #include #include #include #include #include #include #include "fftjet/StatAccumulator.hh" #include "fftjet/SimplePredicates.hh" namespace fftjet { template inline AbsClusteringTree::AbsClusteringTree() : badId(UINT_MAX,UINT_MAX) { } template inline AbsClusteringTree::AbsClusteringTree( const AbsClusteringTree& /* r */) : badId(UINT_MAX,UINT_MAX) { } template inline AbsClusteringTree& AbsClusteringTree::operator=(const AbsClusteringTree&) { return *this; } template inline double AbsClusteringTree::distanceToParent( const TreeNodeId& id) const { return clusterDistance(id, parentId(id)); } template double AbsClusteringTree::minScale() const { const unsigned nlevels = this->nLevels(); if (nlevels > 1) return this->getScale(nlevels - 1); else return 0.0; } template double AbsClusteringTree::maxScale() const { const unsigned nlevels = this->nLevels(); if (nlevels > 1) return this->getScale(1); else return 0.0; } template unsigned AbsClusteringTree::nClusters() const { unsigned nclu = 0; const unsigned nlevels = this->nLevels(); for (unsigned lev=0; levnClusters(lev); return nclu; } template template TreeNodeId AbsClusteringTree::closestNeighbor( const TreeNodeId& idIn, const BooleanPredicate& pred, double* distance) const { assert(distance); *distance = DBL_MAX; TreeNodeId bestNeighbor(badId); TreeNodeId id(idIn.first,0); // The following will verify that idIn.first is legal const unsigned nclus = nClusters(idIn.first); // Make sure that idIn.second is legal as well assert(idIn.second < nclus); // Cycle over clusters on this level for (; id.second template void AbsClusteringTree::closestDauOrSelf( const TreeNodeId& target, const TreeNodeId& id, const BooleanPredicate& pred, double* bestDistanceSoFar, TreeNodeId* bestIdSoFar) const { // Prune the recursion as soon as possible const double d = uncheckedClusterDistance(id, target); if (d - uncheckedNodeRadius(id) > *bestDistanceSoFar) return; // Check if this cluster provides the best distance if (pred(uncheckedCluster(id))) if (d < *bestDistanceSoFar) { *bestDistanceSoFar = d; *bestIdSoFar = id; } // Go over the daughters. It is more likely that more // far away daughter will end up closer to the target. const int ndaus = nDaughters(id); for (int i=ndaus-1; i>=0; --i) closestDauOrSelf(target, daughterId(id,i), pred, bestDistanceSoFar, bestIdSoFar); } template template void AbsClusteringTree::farthestDauOrSelf( const TreeNodeId& target, const TreeNodeId& id, const BooleanPredicate& pred, double* bestDistanceSoFar, TreeNodeId* bestIdSoFar) const { // Prune the recursion as soon as possible const double d = uncheckedClusterDistance(id, target); if (d + uncheckedNodeRadius(id) <= *bestDistanceSoFar) return; // Check if this cluster provides the best distance if (pred(uncheckedCluster(id))) if (d > *bestDistanceSoFar) { *bestDistanceSoFar = d; *bestIdSoFar = id; } // Go over the daughters. It is more likely that more // far away daughter will end up further away from the target. const int ndaus = nDaughters(id); for (int i=ndaus-1; i>=0; --i) farthestDauOrSelf(target, daughterId(id,i), pred, bestDistanceSoFar, bestIdSoFar); } template template TreeNodeId AbsClusteringTree::closestNonDescendant( const TreeNodeId& idIn, const BooleanPredicate& pred, double* distance) const { assert(distance); *distance = DBL_MAX; TreeNodeId bestId(badId); // The following will verify that idIn.first is legal const unsigned nclus = nClusters(idIn.first); // Make sure that idIn.second is legal as well assert(idIn.second < nclus); if (nclus == 2) { TreeNodeId id(idIn.first, 1U - idIn.second); closestDauOrSelf(idIn, id, pred, distance, &bestId); } else if (nclus > 2) { neigbors_.clear(); neigbors_.reserve(nclus - 1); TreeNodeId id(idIn.first,0); for (; id.second template double AbsClusteringTree::conditionalNodeRadius( const TreeNodeId& id, const BooleanPredicate& pred) const { TreeNodeId dummy(id); double distance = 0.0; const int ndaus = nDaughters(id); for (int i=ndaus-1; i>=0; --i) farthestDauOrSelf(id, daughterId(id,i), pred, &distance, &dummy); return distance; } template TreeNodeId AbsClusteringTree::closestDaughter( const TreeNodeId& id) const { if (nDaughters(id)) return daughterId(id, 0); else return badId; } template bool AbsClusteringTree::operator==( const AbsClusteringTree& r) const { const unsigned nlevels = r.nLevels(); if (nlevels != nLevels()) return false; for (unsigned lev=0; lev 0.0); const double scale = getScale(lev); assert(scale > 0.0); if (rscale != scale) return false; if (getLevel(scale) != r.getLevel(scale)) return false; const unsigned nclus = r.nClusters(lev); if (nclus != nClusters(lev)) return false; // Compare the level info, but not for level 0 if (lev) if (getLevelInfo(lev) != r.getLevelInfo(lev)) return false; // Compare individual clusters TreeNodeId id(lev,0); for (; id.second set1, set2; for (unsigned idau=0; idau inline bool AbsClusteringTree::operator!=( const AbsClusteringTree& r) const { return !(*this == r); } template inline void AbsClusteringTree::clear() { for (unsigned i = nLevels()-1; i > 0; --i) remove(i); } template unsigned AbsClusteringTree::nParentsWithDaus( const unsigned level, const unsigned nMinDau) const { assert(level < nLevels()); unsigned count = 0; TreeNodeId id(level,0); const unsigned nclus = nClusters(level); for (; id.second= nMinDau) ++count; return count; } template unsigned AbsClusteringTree::nRootDaughters( const unsigned level) const { assert(level < nLevels()); unsigned nRootDaus = 0; TreeNodeId id(level,0); const TreeNodeId rootId(0,0); const unsigned nclus = nClusters(level); for (; id.second double AbsClusteringTree::levelMatchDistance( const unsigned daughterLevel) const { assert(daughterLevel < nLevels()); if (daughterLevel < 2) return 0.0; const unsigned parentLevel = daughterLevel - 1; const unsigned nClusInParent = nClusters(parentLevel); const unsigned nClusInDau = nClusters(daughterLevel); const unsigned nSterileClus = nClusInParent - nParentsWithDaus(parentLevel, 1); if (nClusInParent == nClusInDau && nSterileClus == 0) { // There is a one-to-one match between the clusters // on these levels. This should be good enough. return 0.0; } // Figure out the number of "strange" clusters. "Strange" // clusters are: // 1. Dau level: directly connected to the top node // 2. Parent level: daughterless (sterile) // 3. Parent level: more than 2 daughters const double nStrange = nRootDaughters(daughterLevel) + nSterileClus + nParentsWithDaus(parentLevel, 3); // Return the ratio of the number of "strange" clusters // to the total number of clusters under consideration return nStrange/(nClusInParent + nClusInDau); } template double AbsClusteringTree::nextBestScale( const double minRatioLog) const { double worstDistance = 0.0; unsigned worstLevel = 0; const unsigned nlevels = nLevels(); double previousScale = getScale(1); for (unsigned ilev=2; ilev minRatioLog) { const double d = levelMatchDistance(ilev); assert(d >= 0.0); if (d > worstDistance) { worstDistance = d; worstLevel = ilev; } } previousScale = scale; } if (worstLevel) return sqrt(getScale(worstLevel-1)*getScale(worstLevel)); else return 0.0; } template bool AbsClusteringTree::isScaleUsed( const double scale) const { const unsigned nlevels = nLevels(); for (unsigned ilev = 1; ilev < nlevels; ++ilev) if (getScale(ilev) == scale) return true; return false; } template void AbsClusteringTree::getAllScales( std::vector* scales, const bool increasingOrder) const { assert(scales); scales->clear(); const unsigned nlevels = nLevels(); if (nlevels > 1) { scales->reserve(nlevels - 1); if (increasingOrder) for (unsigned ilev = nlevels - 1; ilev; --ilev) scales->push_back(getScale(ilev)); else for (unsigned ilev = 1; ilev < nlevels; ++ilev) scales->push_back(getScale(ilev)); } } template unsigned AbsClusteringTree::traceLifetime( const TreeNodeId& id, const unsigned topLevel) { const unsigned ndaus = nDaughters(id); if (ndaus == 0) return id.first; // Calculate lifetimes for all daughters other than the best. // This should not happen very often: ndaus == 1 most of the times. const TreeNodeId bestDau(closestDaughter(id)); if (ndaus > 1) for (unsigned idau=0; idau(uncheckedCluster(dauId))).setLifetime( log(hiscale/loscale)); } } // Calculate the lifetime for the best daughter const unsigned lolev = traceLifetime(bestDau, topLevel); const double hiscale = getScale(topLevel); const double loscale = getScale(lolev); (const_cast(uncheckedCluster(bestDau))).setLifetime( log(hiscale/loscale)); return lolev; } template void AbsClusteringTree::calculateLifetimes() { // Root node does not have the "best" daughter const TreeNodeId rootId(0,0); const unsigned ndaus = nDaughters(rootId); for (unsigned idau=0; idau(uncheckedCluster(dauId))).setLifetime( log(hiscale/loscale)); } } template void AbsClusteringTree::calculateDriftSpeeds() { const double invalidSpeed = -2.0; const unsigned nlevels = nLevels(); const TreeNodeId rootId(0,0); if (nlevels < 2) { // The tree is empty, we have the root level only return; } else if (nlevels == 2) { // Only one level of cluster data. // Can't create any useful drift speed info. const unsigned nclus = nClusters(1); TreeNodeId id(1,0); for (; id.second(uncheckedCluster(id))).setDriftSpeed( invalidSpeed); } else { // More than one level of cluster data. // Process normally. for (unsigned ilev=1; ilev 0.0) { assert(loscale > 0.0); assert(hiscale > loscale); speed = d/log(hiscale/loscale); } (const_cast(uncheckedCluster(id))).setDriftSpeed( speed); } } } } template void AbsClusteringTree::calculateMagnitudeSpeeds() { const double invalidSpeed = -6.0; + const double invalidDeriv = -11.0; const unsigned nlevels = nLevels(); const TreeNodeId rootId(0,0); if (nlevels < 2) { // The tree is empty, we have the root level only return; } else if (nlevels == 2) { // Only one level of cluster data. // Can't create any useful drift speed info. const unsigned nclus = nClusters(1); TreeNodeId id(1,0); for (; id.second(uncheckedCluster(id))).setMagSpeed( - invalidSpeed); + { + Cluster& cl = const_cast(uncheckedCluster(id)); + cl.setMagSpeed(invalidSpeed); + cl.setRelativeScaledHDetDeriv(invalidDeriv); + } } else { // More than one level of cluster data. // Process normally. for (unsigned ilev=1; ilev(uncheckedCluster(id)); if (hiscale > 0.0) { assert(loscale > 0.0); assert(hiscale > loscale); + const double d = hiCl->magnitude()/loCl->magnitude(); assert(d > 0.0); - speed = log(d)/log(hiscale/loscale); + const double logDelta = log(hiscale/loscale); + speed = log(d)/logDelta; + + const double dD = hiCl->relativeScaledHDeterminant()/ + loCl->relativeScaledHDeterminant(); + assert(dD > 0.0); + const double rshd = cl.relativeScaledHDeterminant(); + deriv = rshd*log(dD)/logDelta/levelScale; } - (const_cast(uncheckedCluster(id))).setMagSpeed( - speed); + + cl.setMagSpeed(speed); + cl.setRelativeScaledHDetDeriv(deriv); } } } } template void AbsClusteringTree::updateClusterRadiusInfo() { const unsigned maxlevel = nLevels() - 1U; for (unsigned ilev=1; ilev( uncheckedCluster(id))).setClusterRadius(nodeRadius(id)); } } template void AbsClusteringTree::updateClusterSeparationInfo( const double failDistance) { double distance; Always allPass; const unsigned maxlevel = nLevels() - 1U; for (unsigned ilev=1; ilev( uncheckedCluster(id))).setClusterSeparation(distance); } } } template void AbsClusteringTree::calculateNearestNeighbors( const double noNeighborDistance) { double distance; Always allPass; const unsigned nlevels = nLevels(); for (unsigned ilev=1; ilev( uncheckedCluster(id))).setNearestNeighborDistance( distance); } } } template template void AbsClusteringTree::getPassingNodes( const unsigned level, const BooleanPredicate& pred, std::vector* nodesToFill) const { assert(level < nLevels()); nodesToFill->clear(); TreeNodeId id(level,0); const unsigned nclus = nClusters(level); for (; id.secondpush_back(id); } template template unsigned AbsClusteringTree::clusterCount( const unsigned level, const BooleanPredicate& pred) const { assert(level < nLevels()); const unsigned nclus = nClusters(level); unsigned count = 0; TreeNodeId id(level,0); for (; id.second template void AbsClusteringTree::occupancyInScaleSpace( const BooleanPredicate& pred, std::vector* occupancy) const { assert(occupancy); occupancy->clear(); const unsigned n = nLevels(); occupancy->reserve(n); occupancy->push_back(1U); for (unsigned i=1; ipush_back(clusterCount(i, pred)); } template template int AbsClusteringTree::heritageLine( const TreeNodeId& idIn, const BooleanPredicate& pred, const double maxDistance, std::vector* clustersToFill) const { clustersToFill->clear(); if (!pred(getCluster(idIn))) return -1; // Find the topmost parent which satisfies the criteria // and to which we have an access through the first daughter // sequence unsigned nup = 0; TreeNodeId topmost = idIn; for (TreeNodeId parId = parentId(idIn); parId.first != 0 && parId != badId; parId = parentId(parId)) { const double d = uncheckedClusterDistance(parId, idIn); if (d > maxDistance) break; if (closestDaughter(parId) != topmost) break; if (!pred(uncheckedCluster(parId))) break; topmost = parId; ++nup; } // Find the lowest daughter which satisfies the criteria // and to which we have an access through the first daughter // sequence unsigned ndown = 0; TreeNodeId lowest = idIn; for (TreeNodeId dauId = closestDaughter(idIn); dauId != badId; dauId = closestDaughter(dauId)) { const double d = uncheckedClusterDistance(dauId, idIn); if (d > maxDistance) break; if (!pred(uncheckedCluster(dauId))) break; lowest = dauId; ++ndown; } // Form the sequence of clusters clustersToFill->reserve(nup + ndown + 1); if (nup || ndown) { TreeNodeId id = topmost; clustersToFill->push_back(uncheckedCluster(id)); do { id = closestDaughter(id); clustersToFill->push_back(uncheckedCluster(id)); } while (id != lowest); } else clustersToFill->push_back(uncheckedCluster(idIn)); return nup; } template template void AbsClusteringTree::clusterStats( const unsigned level, const PropertySelector& property, const BooleanPredicate& pred, StatAccumulator* acc) const { assert(level < nLevels()); acc->reset(); const unsigned nclus = nClusters(level); TreeNodeId id(level,0); for (; id.secondaccumulate(property(clu)); } } template template bool AbsClusteringTree::clusterCountLevels( const unsigned requestedCount, const BooleanPredicate& pred, unsigned* minLevel, unsigned* maxLevel) const { assert(minLevel); assert(maxLevel); *minLevel = 0; *maxLevel = 0; bool status = false; const unsigned nl = nLevels(); for (unsigned level = 1; level < nl; ++level) if (clusterCount(level, pred) == requestedCount) { if (*minLevel == 0) { // This is the first time the requested number // of clusters is encountered *minLevel = level; status = true; } else { // Not the first time. We need to reset // the status to "false" if the sequence // is not continuous. if (*maxLevel != level - 1) status = false; } *maxLevel = level; } return status; } template template unsigned AbsClusteringTree::stableClusterCount( const BooleanPredicate& pred, unsigned* minLevel, unsigned* maxLevel, const double alpha, const unsigned minStartingLevel, const unsigned maxStartingLevel) const { assert(minLevel); assert(maxLevel); *minLevel = 0; *maxLevel = 0; unsigned bestN = 0; double bestStability = 0.0; const unsigned lmin = minStartingLevel > 0 ? minStartingLevel : 1; const unsigned nm1 = nLevels() - 1; const unsigned lmax = maxStartingLevel < nm1 ? maxStartingLevel : nm1; unsigned oldN = clusterCount(lmin, pred); double oldScale = getScale(lmin); unsigned oldMinLevel = lmin; for (unsigned level = lmin+1; level <= lmax; ++level) { const double newScale = getScale(level); const unsigned N = clusterCount(level, pred); if (N == oldN) { const double stability = N ? pow(N,alpha)*fabs( log(oldScale/newScale)) : 0.0; if (stability > bestStability) { bestN = N; bestStability = stability; *minLevel = oldMinLevel; *maxLevel = level; } } else { oldN = N; oldScale = newScale; oldMinLevel = level; } } return bestN; } template bool AbsClusteringTree::isAncestor( const TreeNodeId& id1, const TreeNodeId& id2) const { const unsigned level1 = id1.first; const unsigned index1 = id1.second; for (TreeNodeId parent = parentId(id2); parent.first >= level1 && parent.first != UINT_MAX; parent = parentId(parent)) if (parent.first == level1 && parent.second == index1) return true; return false; } } Index: trunk/fftjet/KernelRecombinationAlg.icc =================================================================== --- trunk/fftjet/KernelRecombinationAlg.icc (revision 44) +++ trunk/fftjet/KernelRecombinationAlg.icc (revision 45) @@ -1,753 +1,767 @@ #include -#include +#include namespace fftjet { template KernelRecombinationAlg::KernelRecombinationAlg( ScaleSpaceKernel* ikernel, const Functor2* bgWeight, const double iunlikelyBgWeight, const double idataCutoff, const bool iwinnerTakesAll, const bool ibuildCorrelationMatrix, const bool ibuildClusterMask, const unsigned ietaBinMin, const unsigned ietaBinMax) : defaultKernel(ikernel), bgWeightCalc(bgWeight), unlikelyBgWeight(iunlikelyBgWeight), dataCutoff(static_cast(idataCutoff)), winnerTakesAll(iwinnerTakesAll), buildCorrelationMatrix(ibuildCorrelationMatrix), buildClusterMask(ibuildClusterMask), etaBinMin(ietaBinMin), etaBinMax0(ietaBinMax), vMaker(), weights(0), kernelScales(0), clusterScales(0), clusterScaleRatios(0), dEta(0), energySum(0), + doubleWeightedEnergy(0), energyVariance(0), weightSum(0), + entropySum(0), etaEnergySum(0), phiEnergySum(0), etaPhiESum(0), etaSum(0), phiSum(0), clusterMap(0), peakPositions(0), memFcns2d(0), memFcns3d(0), jets(0), nWeights(0), clusterMask(0), maskMemSize(0), dphiBuffer(0), dphiBufSize(0), use3d(true), lastNJets(0), lastNEta(0), lastNPhi(0), passedJets(0), bgDim(1000) { assert(bgWeightCalc); } template KernelRecombinationAlg::~KernelRecombinationAlg() { delete [] dphiBuffer; delete [] clusterMask; delete [] jets; delete [] memFcns2d; delete [] memFcns3d; delete [] peakPositions; delete [] clusterMap; delete [] weights; } template void KernelRecombinationAlg::processKernelScales( const std::vector& peaks) { unsigned npass = 0; use3d = true; for (unsigned ijet=0; ijetmembershipFactor()); assert(member >= 0.0); if (member > 0.0) { peakPositions[npass] = peak; kernelScales[npass] = member; // Figure out the type of the cluster membership function. // First, try to set up 3d handling. AbsMembershipFunction* m3d = dynamic_cast(peak->membershipFunction()); if (m3d == 0) m3d = dynamic_cast(defaultKernel); if (m3d) { memFcns3d[npass] = m3d; memFcns2d[npass] = 0; } else { // Failed to set up 3d membership function. Try 2d. AbsKernel2d* m2d = dynamic_cast(peak->membershipFunction()); if (m2d == 0) m2d = dynamic_cast(defaultKernel); if (m2d == 0) { // Failed to set up any membership function assert(!"Unspecified recombination membership function"); } else { use3d = false; memFcns3d[npass] = 0; memFcns2d[npass] = m2d; } } // Figure out the recombination scale double recoScale = peak->recoScale(); if (recoScale == 0.0) recoScale = peak->scale(); assert(recoScale > 0.0); clusterScales[npass] = recoScale; clusterScaleRatios[npass] = peak->recoScaleRatio(); ++npass; } } // Verify the consistency of the membership function type. // If "use3d" is still true, we know for sure that everything // is 3d. However, if "use3d" is not true, there may be // a mixture of 3d and 2d cases. if (!use3d) for (unsigned ijet=0; ijet void KernelRecombinationAlg::calculateDphi( const Grid2d& eventData) { const unsigned bufsiz = lastNJets*lastNPhi; if (bufsiz > dphiBufSize) { delete [] dphiBuffer; dphiBuffer = new double[bufsiz]; dphiBufSize = bufsiz; } for (unsigned iphi=0; iphiphi(); while (dphi > M_PI) dphi -= (2.0*M_PI); while (dphi < -M_PI) dphi += (2.0*M_PI); dPhi[ijet] = dphi; } } } template void KernelRecombinationAlg::setupBuffers( const unsigned njets, const unsigned nEta, const unsigned nPhi) { const bool setup4vec = recombine4Vectors(); if (njets > nWeights) { delete [] weights; - weights = new double[13*njets]; + weights = new double[15*njets]; kernelScales = weights + njets; clusterScales = kernelScales + njets; clusterScaleRatios = clusterScales + njets; dEta = clusterScaleRatios + njets; energySum = dEta + njets; - energyVariance = energySum + njets; + doubleWeightedEnergy = energySum + njets; + energyVariance = doubleWeightedEnergy + njets; weightSum = energyVariance + njets; - etaEnergySum = weightSum + njets; + entropySum = weightSum + njets; + etaEnergySum = entropySum + njets; phiEnergySum = etaEnergySum + njets; etaPhiESum = phiEnergySum + njets; etaSum = etaPhiESum + njets; phiSum = etaSum + njets; if (buildClusterMask) { delete [] clusterMap; clusterMap = new unsigned[njets+1]; } delete [] peakPositions; peakPositions = new const Peak*[njets]; delete [] memFcns2d; memFcns2d = new AbsKernel2d*[njets]; delete [] memFcns3d; memFcns3d = new AbsMembershipFunction*[njets]; if (setup4vec) { delete [] jets; jets = new VectorLike[njets]; } nWeights = njets; } if (setup4vec) { const VectorLike zero = VectorLike(); for (unsigned i=0; i maskMemSize) { delete [] clusterMask; clusterMask = new unsigned[nbins]; maskMemSize = nbins; } for (unsigned i=0; i void KernelRecombinationAlg::allUnclustered( const Grid2d& eventData, VectorLike* unclustered, double* unused) { assert(lastNJets == 0); const unsigned etaBinMax = etaBinMax0 < lastNEta ? etaBinMax0 : lastNEta; for (unsigned ieta=etaBinMin; ieta dataCutoff) { *unclustered += vMaker( gridValue, eta, eventData.phiBinCenter(iphi)); *unused += gridValue; } } } } template inline const unsigned* KernelRecombinationAlg::getClusterMask() const { return clusterMask; } template void KernelRecombinationAlg::remapClusterMask() { if (buildClusterMask) { // Calculate the number of output jets and // build the remapping array for the cluster mask clusterMap[0] = 0; for (unsigned ijet=0; ijet 0.0 && energySum[ijet] > 0.0) clusterMap[ijet+1] = ++passedJets; else clusterMap[ijet+1] = 0; } // Remap the cluster mask const unsigned nbins = lastNEta*lastNPhi; for (unsigned i=0; i 0.0 && energySum[ijet] > 0.0) ++passedJets; } } template bool KernelRecombinationAlg ::performPreliminaryProcessing( const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, const unsigned nBgEta, const unsigned nBgPhi, std::vector >* outputJets, VectorLike* unclustered, double* unused) { // Check the input assert(bgData); assert(outputJets); assert(unclustered); assert(unused); bgDim = AbsRecombinationAlg::bgDimensionality( eventData, nBgEta, nBgPhi); // Reset the output outputJets->clear(); *unclustered = VectorLike(); *unused = 0.0; // Allocate and initialize some necessary buffers setupBuffers(peaks.size(), eventData.nEta(), eventData.nPhi()); // Figure out the kernel scales from the scaler. // Get rid of the peaks suppressed by the scaler. processKernelScales(peaks); if (lastNJets) { // Figure out the "dphi" values for each phi bin // and for each jet. This must be done after calculating // "peakPositions" and kernel bandwidth values. calculateDphi(eventData); return false; } else { // Just sum up the unclustered energy and return // indicating that we are done in case there are no jets allUnclustered(eventData, unclustered, unused); return true; } } template void KernelRecombinationAlg::buildOutput( std::vector >* outputJets, const bool setRecoScaleRatiosToZero) { remapClusterMask(); const bool use4Vec = recombine4Vectors(); const bool useCentroid = useEtCentroid(); // Build the output jets outputJets->reserve(passedJets); for (unsigned ijet=0; ijet 0.0 && energySum[ijet] > 0.0) { const double et = energySum[ijet]; const double centroidEta = etaSum[ijet]/et; const double detaPeak = centroidEta - peakPositions[ijet]->eta(); const double dphiPeak = phiSum[ijet]/et; const double centroidPhi = dphiPeak + peakPositions[ijet]->phi(); const double eSigma = energyVariance[ijet] > 0.0 ? sqrt(energyVariance[ijet]) : 0.0; const double fuzziness = eSigma/et; + const double iso = doubleWeightedEnergy[ijet]/et; + const double entropy = log(et) - entropySum[ijet]/et; double etaWidth = etaEnergySum[ijet]/et - detaPeak*detaPeak; etaWidth = etaWidth > 0.0 ? sqrt(etaWidth) : 0.0; double phiWidth = phiEnergySum[ijet]/et - dphiPeak*dphiPeak; phiWidth = phiWidth > 0.0 ? sqrt(phiWidth) : 0.0; double rho = 0.0; if (etaWidth > 0.0 && phiWidth > 0.0) { double comboWidth = etaPhiESum[ijet]/et - detaPeak*dphiPeak; rho = comboWidth/etaWidth/phiWidth; if (rho > 1.0) rho = 1.0; if (rho < -1.0) rho = -1.0; } Peak inPeak(*peakPositions[ijet]); inPeak.setRecoScale(clusterScales[ijet]); if (setRecoScaleRatiosToZero) inPeak.setRecoScaleRatio(0.0); if (use4Vec) { outputJets->push_back(RecombinedJet( inPeak, jets[ijet], weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } else if (useCentroid) { outputJets->push_back(RecombinedJet( inPeak, vMaker(et, centroidEta, centroidPhi), weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } else { outputJets->push_back(RecombinedJet( inPeak, vMaker(et, inPeak.eta(), inPeak.phi()), weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } } } template int KernelRecombinationAlg::run( const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, const unsigned nBgEta, const unsigned nBgPhi, std::vector >* outputJets, VectorLike* unclustered, double* unused) { if (performPreliminaryProcessing( peaks, eventData, bgData, nBgEta, nBgPhi, outputJets, unclustered, unused)) return 0; // Some variables which will be needed below const unsigned etaBinMax = etaBinMax0 < lastNEta ? etaBinMax0 : lastNEta; const unsigned bgStride = bgDim == 2 ? 1 : 0; const bool use4Vec = recombine4Vectors(); // Go over all grid points and figure out to which jet // that point belongs. for (unsigned ieta=etaBinMin; ietaeta(); // We will try to get rid of each point in the inner // cycle as fast as possible. The following few lines // of code depend, unfortunately, on the internal layout // of the Grid2d object. On the positive side, we can // get rid of a point below the cutoff without making // the "uncheckedAt(ieta, iphi)" function call. const Real* databuf = eventData.data() + ieta*lastNPhi; const BgData* bgbuf; switch (bgDim) { case 0: bgbuf = bgData; break; case 1: bgbuf = bgData + ieta; break; case 2: bgbuf = bgData + ieta*lastNPhi; break; default: assert(0); } for (unsigned iphi=0; iphi(gridValue), bgbuf[iphi*bgStride])); const VectorLike gridVec(vMaker(gridValue, eta, phi)); // Cycle over jets. This is the most speed-critical // cycle in the whole algorithm. double sum(0.0), maxW(0.0); unsigned maxWJet(lastNJets); if (use3d) { for (unsigned ijet=0; ijet 0.0) memFcns3d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); const double f = kernelScales[ijet]* (*memFcns3d[ijet])(dEta[ijet], dPhi[ijet], gridValue, clusterScales[ijet]); sum += f; weights[ijet] = f; if (f > maxW) { maxW = f; maxWJet = ijet; } } } else { for (unsigned ijet=0; ijet 0.0) memFcns2d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); const double f = kernelScales[ijet]* (*memFcns2d[ijet])(dEta[ijet], dPhi[ijet], clusterScales[ijet]); sum += f; weights[ijet] = f; if (f > maxW) { maxW = f; maxWJet = ijet; } } } if (maxWJet == lastNJets) { // No contributing cluster const bool poorModelfit = bgW < unlikelyBgWeight; // We should attempt to rectify the situation when // the model used in this recombination algorithm // breaks down. There are several breakdown modes: // // 1) The initial guess for the jet Pt was too low, // and the cell's energy exceeds the allowed energy // range for this jet and for the delta phi and // delta eta between the cell and the jet. // // 2) More than one jet contributes energy to the // current cell, so that, again, the cell's energy // exceeds the allowed energy range for each of // these jets. // // 3) The pattern recognition stage did not find // all the jets (or the user tries to fit a certain // topology to every event). Then the cell's energy // may have been generated by a jet not taken // into account in the recombination. // // Cases 1) and 2) are treated in the code below // by splitting the cell's energy between the closest // two jets (rather than one) and the background. // Ideally, one should choose the split in such a way // that the energy assignment probability is maximized. // Unfortunately, such an approach would be extremely // CPU-intensive. Instead, we will simply allow the jets // to take as much energy as they can from this cell. // // Case 3) can not be handled here -- it is up to the // user to choose the background membership function // and the "unlikelyBgWeight" parameter in such a way // that the overall algorithm is sufficiently robust // w.r.t. this breakdown mode. // // When 2d kernels are used, energy splitting between // nearby jets becomes unnatural -- there is no reason // for the cell to have low membership function value // apart from being too far away from the jet center // in the first place. Therefore, such cells are // assigned to the background/noise anyway. // if (poorModelfit && use3d) { double maxE[2] = {0.0, 0.0}; unsigned maxEJet[2] = {UINT_MAX, UINT_MAX}; for (unsigned ijet=0; ijet 0.0) memFcns3d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); const double e = memFcns3d[ijet]->absorbableEnergy( dEta[ijet], dPhi[ijet], clusterScales[ijet]); if (e > maxE[0]) { maxE[1] = maxE[0]; maxEJet[1] = maxEJet[0]; maxE[0] = e; maxEJet[0] = ijet; } else if (e > maxE[1]) { maxE[1] = e; maxEJet[1] = ijet; } } double cellE = gridValue; for (unsigned i=0; i < 2 && cellE > 0.0; ++i) if (maxE[i] > 0.0) { const unsigned ijet = maxEJet[i]; const double absorbedE = cellE > maxE[i] ? maxE[i] : cellE; cellE -= absorbedE; if (use4Vec) jets[ijet] += vMaker(absorbedE, eta, phi); weightSum[ijet] += 1.0; energySum[ijet] += absorbedE; + entropySum[ijet] += absorbedE*log(absorbedE); + doubleWeightedEnergy[ijet] += absorbedE; etaSum[ijet] += absorbedE*eta; etaEnergySum[ijet] += absorbedE*dEta[ijet]*dEta[ijet]; phiEnergySum[ijet] += absorbedE*dPhi[ijet]*dPhi[ijet]; etaPhiESum[ijet] += absorbedE*dEta[ijet]*dPhi[ijet]; phiSum[ijet] += absorbedE*dPhi[ijet]; // // How does this affect the energy variance? // No easy way to make a meaningful estimate. // energyVariance[ijet] += 0.0; } if (cellE > 0.0) { // There is still some energy left to // assign it to the noise/background if (maxE[0] > 0.0) // Some energy was assigned to jets *unclustered += vMaker(cellE, eta, phi); else // No energy was assigned to jets. // No need to recalculate the 4-vector. *unclustered += gridVec; *unused += cellE; } } else { // Just assign the whole cell to the unclustered // energy *unclustered += gridVec; *unused += gridValue; } } else { const double totalSum = sum + bgW; if (buildClusterMask && maxW > bgW) clusterMask[ieta*lastNPhi + iphi] = maxWJet + 1; if (winnerTakesAll) { // Hard clustering. All energy goes // to the cluster with the highest weight. if (maxW > bgW) { if (use4Vec) jets[maxWJet] += gridVec; weightSum[maxWJet] += 1.0; energySum[maxWJet] += gridValue; + entropySum[maxWJet] += gridValue*log(gridValue); + doubleWeightedEnergy[maxWJet] += gridValue*maxW/totalSum; etaSum[maxWJet] += gridValue*eta; etaEnergySum[maxWJet] += gridValue*dEta[maxWJet]*dEta[maxWJet]; phiEnergySum[maxWJet] += gridValue*dPhi[maxWJet]*dPhi[maxWJet]; etaPhiESum[maxWJet] += gridValue*dEta[maxWJet]*dPhi[maxWJet]; energyVariance[maxWJet] += gridValue*gridValue*(1.0 - maxW/totalSum); phiSum[maxWJet] += gridValue*dPhi[maxWJet]; } else { *unclustered += gridVec; *unused += gridValue; } } else { // Fuzzy clustering. Distribute the cell energy // among all clusters according to their weights. // Unclustered energy fraction const double bgfrac = bgW/totalSum; if (bgfrac > 0.0) { *unclustered += gridVec*bgfrac; *unused += gridValue*bgfrac; } // Fractions which belong to each jet for (unsigned ijet=0; ijet 0.0) { const double ePart = gridValue*w; if (use4Vec) jets[ijet] += gridVec*w; weightSum[ijet] += w; energySum[ijet] += ePart; + entropySum[ijet] += ePart*log(ePart); + doubleWeightedEnergy[ijet] += ePart*w; etaEnergySum[ijet] += ePart*dEta[ijet]*dEta[ijet]; etaSum[ijet] += ePart*eta; phiEnergySum[ijet] += ePart*dPhi[ijet]*dPhi[ijet]; etaPhiESum[ijet] += ePart*dEta[ijet]*dPhi[ijet]; energyVariance[ijet] += gridValue*ePart*(1.0 - w); phiSum[ijet] += ePart*dPhi[ijet]; } } } } } } buildOutput(outputJets, false); return 0; } } Index: trunk/fftjet/Peak.hh =================================================================== --- trunk/fftjet/Peak.hh (revision 44) +++ trunk/fftjet/Peak.hh (revision 45) @@ -1,258 +1,262 @@ #ifndef FFTJET_PEAK_HH_ #define FFTJET_PEAK_HH_ //========================================================================= // Peak.hh // // Utility class which holds initial cluster (precluster) parameters. // PeakFinder, ProximityClusteringTree, ClusteringSequencer, etc. use // this class for its calculations. // // Objects of this class hold the following information: // // eta() -- Normally, rapidity or pseudorapidity of the // precluster. The exact meaning of this quantity // depends on the choices made by the user during // the energy discretization step. // // phi() -- Azimuthal angle of the precluster. // // magnitude() -- The "magnitude" of the precluster: the height // of the peak of the discretized energy // distribution convoluted with the kernel. // // scale() -- Pattern recognition resolution scale at which // this precluster was found. // // hessian(array) -- Returns the peak Hessian matrix (w.r.t. eta // and phi variables) as a 1d array. Since this // matrix is symmetric, it is enough to return // three values. The order is h[0][0], h[0][1], // h[1][1]. Note that both Hessian and Laplacian // are provided by the peak finder only when // it operates in the "subcell" resolution mode. // // hessianDeterminant() -- Returns the determinant of the peak Hessian // matrix. // // laplacian() -- Returns the peak Laplacian. // // scaledLaplacian(bwEta, bwPhi) -- Returns the scale-normalized peak // Laplacian which can be used for blob // detection in the Gaussian scale space. // // The following quantities are defined if the preclusters were arranged // in a clustering tree, and the corresponding calculations were performed // by the tree code (by default, all these quantities are calculated): // // driftSpeed() -- The speed with which the precluster location // moves in the scale space: // d distance/d log(scale). Here, "distance" // is defined in terms of the user-selected // distance function. // // magSpeed() -- The speed with which the precluster magnitude // changes in the scale space: // d log(magnitude)/d log(scale) // // lifetime() -- The "lifetime" of the precluster in the scale // space. It is computed as // log(high_scale) - log(low_scale) // where "low_scale" and "high_scale" define // the range of resolution scales for which the // precluster exists as a feature of the energy // distribution. Typically, the lifetime is traced // from the smallest scale in the clustering tree // to the scale where the precluster becomes // a part of a larger precluster. It may be useful // to check the precluster lifetime in case // the pattern recognition was performed with // a kernel which tends to produce spurious // modes. Such spurious preclusters simply // disappear at smaller scales without leaving // any descendants, which results in their low // lifetimes. // // nearestNeighborDistance() -- The distance to the nearest precluster // at the same resolution scale. // // The following quantities are optional. They become really meaningful // only if the complete event is inserted at the bottom level of the // clustering tree, and this is a computationally expensive operation. // // clusterRadius() -- The distance from the precluster location // to the most distant daughter // // clusterSeparation() -- The distance from the precluster location // to the nearest daughter of a different // precluster // // The following quantities are user-settable. They are not calculated // by FFTJet algorithms. // // code() -- The suggested use of these quantities is to // status() indicate some feature and/or status of the // user-developed pattern recognition algorithm. // // The following quantities may be used to steer the behavior of the // energy recombination codes: // // recoScale() -- The initial recombination scale which can be // set during pattern recognition stage and then // passed to the jet membership function. If this // scale is not set, the recombination stage will // use the precluster resolution scale instead. // // recoScaleRatio() -- Ratio of eta to phi recombination scale. // Can be set on per-cluster basis (in particular, // this is useful if the typical cluster shape // is changing depending on the cluster location // in the detector). // // membershipFactor() -- Can be used to multiply the jet membership // function by a precluster-dependent factor. // It may be useful to set this factor, for example, // in proportion to jet energy if the membership // function is a normalized energy profile, etc. // It can also be used for fast trimming of // preclusters between pattern recognition and // recombination stages -- just set the membership // function factor to 0 to get rid of a precluster. // // membershipFunction() -- User-settable jet membership function. // By default, this function is not set, and // all jets are recombined using the same // jet membership function provided when // the recombination algorithm object is created. // However, in pursuit of ultimate resolution, // the user may introduce precluster-specific // assumptions (for example, about jet flavor). // In this case it makes sense to use different // membership functions for each jet. Note that // all such function must have consistent // dimensionality: either all of them should // be derived from "AbsMembershipFunction" which // includes an additional energy parameter or all // of them should be derived from "AbsKernel2d" // which operates only in the eta-phi space. // I. Volobouev // May 2008 //========================================================================= #include namespace fftjet { class ScaleSpaceKernel; class Peak { public: // This object will not own the "membershipFunction" pointer Peak(double eta, double phi, double magnitude, const double hessian[3], double driftSpeed=-1.0, double magSpeed=-5.0, double lifetime=-1.0, double scale=-1.0, double nearestDistance=-1.0, double membershipFactor=1.0, double recoScale=0.0, double recoScaleRatio=0.0, double clusterRadius=-1.0, double clusterSeparation=-1.0, int code=INT_MIN, int status=-1, ScaleSpaceKernel* membershipFunction=0); // Inspectors inline double eta() const {return eta_;} inline double phi() const {return phi_;} inline double magnitude() const {return magnitude_;} inline double driftSpeed() const {return speed_;} inline double magSpeed() const {return magSpeed_;} inline double lifetime() const {return lifetime_;} inline double splitTime() const {return splitTime_;} inline double mergeTime() const {return mergeTime_;} inline double scale() const {return scale_;} inline double nearestNeighborDistance() const {return nearestD_;} inline double membershipFactor() const {return membershipFactor_;} inline double recoScale() const {return recoScale_;} inline double recoScaleRatio() const {return recoScaleRatio_;} inline double clusterRadius() const {return clusterRadius_;} inline double clusterSeparation() const {return clusterSeparation_;} void hessian(double hessianArray[3]) const; + inline double relativeScaledHDetDeriv() const {return relHDetDeriv_;} inline int code() const {return code_;} inline int status() const {return status_;} inline ScaleSpaceKernel* membershipFunction() const {return memFcn_;} // Modifiers inline void setEtaPhi(const double eta, const double phi) {eta_ = eta; phi_ = phi;} inline void setMagnitude(const double mag) {magnitude_ = mag;} inline void setDriftSpeed(const double s) {speed_ = s;} inline void setMagSpeed(const double s) {magSpeed_ = s;} inline void setLifetime(const double t) {lifetime_ = t;} inline void setSplitTime(const double t) {splitTime_ = t;} inline void setMergeTime(const double t) {mergeTime_ = t;} inline void setScale(const double scale) {scale_ = scale;} inline void setNearestNeighborDistance(const double d) {nearestD_ = d;} inline void setMembershipFactor(const double s) {membershipFactor_ = s;} inline void setRecoScale(const double d) {recoScale_ = d;} inline void setRecoScaleRatio(const double d) {recoScaleRatio_ = d;} inline void setClusterRadius(const double d) {clusterRadius_ = d;} inline void setClusterSeparation(const double d) {clusterSeparation_ = d;} void setHessian(const double hessian[3]); + inline void setRelativeScaledHDetDeriv(const double s) + {relHDetDeriv_ = s;} inline void setCode(const int c) {code_= c;} inline void setStatus(const int s) {status_= s;} inline void setMembershipFunction(ScaleSpaceKernel* f) {memFcn_ = f;} // Methods related to blob detection double hessianDeterminant() const; double scaledHDeterminant() const; double relativeScaledHDeterminant() const; inline double laplacian() const {return hessian_[0] + hessian_[2];} double scaledLaplacian(double bwEta, double bwPhi) const; double relativeScaledLaplacian(double bwEta, double bwPhi) const; // Default comparison is by magnitude times scale squared inline bool operator<(const Peak& r) const {return scale_*scale_*magnitude_ < r.scale_*r.scale_*r.magnitude_;} inline bool operator>(const Peak& r) const {return scale_*scale_*magnitude_ > r.scale_*r.scale_*r.magnitude_;} bool operator==(const Peak& r) const; inline bool operator!=(const Peak& r) const {return !(*this == r);} // Normally, the default constructor of this object // should not be used: a peak without coordinates makes // no sense. The following function can help when // there is no way to avoid building a dummy object. static Peak dummy(); private: Peak(); double eta_; double phi_; double magnitude_; double speed_; double magSpeed_; double lifetime_; double splitTime_; double mergeTime_; double scale_; double nearestD_; double membershipFactor_; double recoScale_; double recoScaleRatio_; double clusterRadius_; double clusterSeparation_; double hessian_[3]; + double relHDetDeriv_; int code_; int status_; ScaleSpaceKernel* memFcn_; }; } #endif // FFTJET_PEAK_HH_ Index: trunk/fftjet/KernelVectorRecombinationAlg.icc =================================================================== --- trunk/fftjet/KernelVectorRecombinationAlg.icc (revision 44) +++ trunk/fftjet/KernelVectorRecombinationAlg.icc (revision 45) @@ -1,667 +1,681 @@ #include -#include +#include namespace fftjet { template KernelVectorRecombinationAlg::KernelVectorRecombinationAlg( ScaleSpaceKernel* ikernel, VectorLikeMemberFunction const ietFcn, VectorLikeMemberFunction const ietaFcn, VectorLikeMemberFunction const iphiFcn, const Functor2* bgWeight, const double iunlikelyBgWeight, const bool iwinnerTakesAll, const bool ibuildCorrelationMatrix, const bool ibuildClusterMask) : defaultKernel(ikernel), etFcn(ietFcn), etaFcn(ietaFcn), phiFcn(iphiFcn), bgWeightCalc(bgWeight), unlikelyBgWeight(iunlikelyBgWeight), winnerTakesAll(iwinnerTakesAll), buildCorrelationMatrix(ibuildCorrelationMatrix), buildClusterMask(ibuildClusterMask), vMaker(), weights(0), kernelScales(0), clusterScales(0), clusterScaleRatios(0), energySum(0), + doubleWeightedEnergy(0), energyVariance(0), weightSum(0), + entropySum(0), etaEnergySum(0), phiEnergySum(0), etaPhiESum(0), etaSum(0), phiSum(0), detaBuf(0), dphiBuf(0), clusterMap(0), peakPositions(0), memFcns2d(0), memFcns3d(0), jets(0), nWeights(0), clusterMask(0), maskMemSize(0), use3d(true), lastNJets(0), lastNData(0), bgIsAVector(false) { assert(bgWeightCalc); } template KernelVectorRecombinationAlg::~KernelVectorRecombinationAlg() { delete [] clusterMask; delete [] jets; delete [] memFcns2d; delete [] memFcns3d; delete [] peakPositions; delete [] clusterMap; delete [] weights; } template void KernelVectorRecombinationAlg::processKernelScales( const std::vector& peaks) { unsigned npass = 0; use3d = true; for (unsigned ijet=0; ijetmembershipFactor()); assert(member >= 0.0); if (member > 0.0) { peakPositions[npass] = peak; kernelScales[npass] = member; // Figure out the type of the cluster membership function. // First, try to set up 3d handling. AbsMembershipFunction* m3d = dynamic_cast(peak->membershipFunction()); if (m3d == 0) m3d = dynamic_cast(defaultKernel); if (m3d) { memFcns3d[npass] = m3d; memFcns2d[npass] = 0; } else { // Failed to set up 3d membership function. Try 2d. AbsKernel2d* m2d = dynamic_cast(peak->membershipFunction()); if (m2d == 0) m2d = dynamic_cast(defaultKernel); if (m2d == 0) { // Failed to set up any membership function assert(!"Unspecified recombination membership function"); } else { use3d = false; memFcns3d[npass] = 0; memFcns2d[npass] = m2d; } } // Figure out the recombination scale double recoScale = peak->recoScale(); if (recoScale == 0.0) recoScale = peak->scale(); assert(recoScale > 0.0); clusterScales[npass] = recoScale; clusterScaleRatios[npass] = peak->recoScaleRatio(); ++npass; } } // Verify the consistency of the membership function type. // If "use3d" is still true, we know for sure that everything // is 3d. However, if "use3d" is not true, there may be // a mixture of 3d and 2d cases. if (!use3d) for (unsigned ijet=0; ijet void KernelVectorRecombinationAlg::setupBuffers( const unsigned njets, const unsigned ndata) { const bool setup4vec = recombine4Vectors(); if (njets > nWeights) { delete [] weights; - weights = new double[14*njets]; + weights = new double[16*njets]; kernelScales = weights + njets; clusterScales = kernelScales + njets; clusterScaleRatios = clusterScales + njets; energySum = clusterScaleRatios + njets; - energyVariance = energySum + njets; + doubleWeightedEnergy = energySum + njets; + energyVariance = doubleWeightedEnergy + njets; weightSum = energyVariance + njets; - etaEnergySum = weightSum + njets; + entropySum = weightSum + njets; + etaEnergySum = entropySum + njets; phiEnergySum = etaEnergySum + njets; etaPhiESum = phiEnergySum + njets; etaSum = etaPhiESum + njets; phiSum = etaSum + njets; detaBuf = phiSum + njets; dphiBuf = detaBuf + njets; if (buildClusterMask) { delete [] clusterMap; clusterMap = new unsigned[njets+1]; } delete [] peakPositions; peakPositions = new const Peak*[njets]; delete [] memFcns2d; memFcns2d = new AbsKernel2d*[njets]; delete [] memFcns3d; memFcns3d = new AbsMembershipFunction*[njets]; if (setup4vec) { delete [] jets; jets = new VectorLike[njets]; } nWeights = njets; } if (setup4vec) { const VectorLike zero = VectorLike(); for (unsigned i=0; i maskMemSize) { delete [] clusterMask; clusterMask = new unsigned[ndata]; maskMemSize = ndata; } for (unsigned i=0; i void KernelVectorRecombinationAlg::allUnclustered( const std::vector& eventData, VectorLike* unclustered, double* unused) { assert(lastNJets == 0); for (unsigned i=0; i inline const unsigned* KernelVectorRecombinationAlg::getClusterMask() const { return clusterMask; } template void KernelVectorRecombinationAlg::remapClusterMask() { if (buildClusterMask) { // Calculate the number of output jets and // build the remapping array for the cluster mask clusterMap[0] = 0; for (unsigned ijet=0; ijet 0.0 && energySum[ijet] > 0.0) clusterMap[ijet+1] = ++passedJets; else clusterMap[ijet+1] = 0; } // Remap the cluster mask for (unsigned i=0; i 0.0 && energySum[ijet] > 0.0) ++passedJets; } } template bool KernelVectorRecombinationAlg ::performPreliminaryProcessing( const std::vector& peaks, const std::vector& eventData, const BgData* bgData, const unsigned bgDataLength, std::vector >* outputJets, VectorLike* unclustered, double* unused) { // Check the input assert(bgData); assert(outputJets); assert(unclustered); assert(unused); // Make sure background length is compatible with data length lastNData = eventData.size(); bgIsAVector = (bgDataLength == lastNData); assert(bgIsAVector || bgDataLength == 1); // Reset the output outputJets->clear(); *unclustered = VectorLike(); *unused = 0.0; // Allocate and initialize some necessary buffers setupBuffers(peaks.size(), lastNData); // Figure out the kernel scales from the scaler. // Get rid of the peaks suppressed by the scaler. processKernelScales(peaks); if (lastNJets && lastNData) { // Will need to run the main data processing code return false; } else { // Just sum up the unclustered energy and return // indicating that we are done in case there are no jets allUnclustered(eventData, unclustered, unused); return true; } } template void KernelVectorRecombinationAlg::buildOutput( std::vector >* outputJets) { remapClusterMask(); const bool use4Vec = recombine4Vectors(); const bool useCentroid = useEtCentroid(); // Build the output jets outputJets->reserve(passedJets); for (unsigned ijet=0; ijet 0.0 && energySum[ijet] > 0.0) { const double et = energySum[ijet]; const double centroidEta = etaSum[ijet]/et; const double detaPeak = centroidEta - peakPositions[ijet]->eta(); const double dphiPeak = phiSum[ijet]/et; const double centroidPhi = dphiPeak + peakPositions[ijet]->phi(); const double eSigma = energyVariance[ijet] > 0.0 ? sqrt(energyVariance[ijet]) : 0.0; const double fuzziness = eSigma/et; + const double iso = doubleWeightedEnergy[ijet]/et; + const double entropy = log(et) - entropySum[ijet]/et; double etaWidth = etaEnergySum[ijet]/et - detaPeak*detaPeak; etaWidth = etaWidth > 0.0 ? sqrt(etaWidth) : 0.0; double phiWidth = phiEnergySum[ijet]/et - dphiPeak*dphiPeak; phiWidth = phiWidth > 0.0 ? sqrt(phiWidth) : 0.0; double rho = 0.0; if (etaWidth > 0.0 && phiWidth > 0.0) { double comboWidth = etaPhiESum[ijet]/et - detaPeak*dphiPeak; rho = comboWidth/etaWidth/phiWidth; if (rho > 1.0) rho = 1.0; if (rho < -1.0) rho = -1.0; } Peak inPeak(*peakPositions[ijet]); inPeak.setRecoScale(clusterScales[ijet]); if (use4Vec) { outputJets->push_back(RecombinedJet( inPeak, jets[ijet], weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } else if (useCentroid) { outputJets->push_back(RecombinedJet( inPeak, vMaker(et, centroidEta, centroidPhi), weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } else { outputJets->push_back(RecombinedJet( inPeak, vMaker(et, inPeak.eta(), inPeak.phi()), weightSum[ijet], et, centroidEta, centroidPhi, - etaWidth, phiWidth, rho, fuzziness)); + etaWidth, phiWidth, rho, fuzziness, iso, entropy)); } } } template int KernelVectorRecombinationAlg::run( const std::vector& peaks, const std::vector& eventData, const BgData* bgData, const unsigned bgDataLength, std::vector >* outputJets, VectorLike* unclustered, double* unused) { if (performPreliminaryProcessing( peaks, eventData, bgData, bgDataLength, outputJets, unclustered, unused)) return 0; // Some variables which will be needed below const bool use4Vec = recombine4Vectors(); const unsigned bgStride = bgIsAVector ? 1 : 0; // Go over all input data points and figure out to which jet // that point belongs for (unsigned ipt=0; ipteta(); double dphi = phi - peakPositions[ijet]->phi(); while (dphi > M_PI) dphi -= (2.0*M_PI); while (dphi < -M_PI) dphi += (2.0*M_PI); dphiBuf[ijet] = dphi; double f = 0.0; if (use3d) { if (clusterScaleRatios[ijet] > 0.0) memFcns3d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); f = kernelScales[ijet]* (*memFcns3d[ijet])(detaBuf[ijet], dphi, gridValue, clusterScales[ijet]); } else { if (clusterScaleRatios[ijet] > 0.0) memFcns2d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); f = kernelScales[ijet]* (*memFcns2d[ijet])(detaBuf[ijet], dphi, clusterScales[ijet]); } sum += f; weights[ijet] = f; if (f > maxW) { maxW = f; maxWJet = ijet; } } if (maxWJet == lastNJets) { // No contributing cluster const bool poorModelfit = bgW < unlikelyBgWeight; // We should attempt to rectify the situation when // the model used in this recombination algorithm // breaks down. There are several breakdown modes: // // 1) The initial guess for the jet Pt was too low, // and the cell's energy exceeds the allowed energy // range for this jet and for the delta phi and // delta eta between the cell and the jet. // // 2) More than one jet contributes energy to the // current cell, so that, again, the cell's energy // exceeds the allowed energy range for each of // these jets. // // 3) The pattern recognition stage did not find // all the jets (or the user tries to fit a certain // topology to every event). Then the cell's energy // may have been generated by a jet not taken // into account in the recombination. // // Cases 1) and 2) are treated in the code below // by splitting the cell's energy between the closest // two jets (rather than one) and the background. // Ideally, one should choose the split in such a way // that the energy assignment probability is maximized. // Unfortunately, such an approach would be extremely // CPU-intensive. Instead, we will simply allow the jets // to take as much energy as they can from this cell. // // Case 3) can not be handled here -- it is up to the // user to choose the background membership function // and the "unlikelyBgWeight" parameter in such a way // that the overall algorithm is sufficiently robust // w.r.t. this breakdown mode. // // When 2d kernels are used, energy splitting between // nearby jets becomes unnatural -- there is no reason // for the cell to have low membership function value // apart from being too far away from the jet center // in the first place. Therefore, such cells are // assigned to the background/noise anyway. // if (poorModelfit && use3d) { double maxE[2] = {0.0, 0.0}; unsigned maxEJet[2] = {UINT_MAX, UINT_MAX}; for (unsigned ijet=0; ijet 0.0) memFcns3d[ijet]->setScaleRatio( clusterScaleRatios[ijet]); const double e = memFcns3d[ijet]->absorbableEnergy( detaBuf[ijet], dphiBuf[ijet], clusterScales[ijet]); if (e > maxE[0]) { maxE[1] = maxE[0]; maxEJet[1] = maxEJet[0]; maxE[0] = e; maxEJet[0] = ijet; } else if (e > maxE[1]) { maxE[1] = e; maxEJet[1] = ijet; } } double cellE = gridValue; for (unsigned i=0; i < 2 && cellE > 0.0; ++i) if (maxE[i] > 0.0) { const unsigned ijet = maxEJet[i]; const double absorbedE = cellE > maxE[i] ? maxE[i] : cellE; cellE -= absorbedE; if (use4Vec) jets[ijet] += vMaker(absorbedE, eta, phi); weightSum[ijet] += 1.0; energySum[ijet] += absorbedE; + entropySum[ijet] += absorbedE*log(absorbedE); + doubleWeightedEnergy[ijet] += absorbedE; etaSum[ijet] += absorbedE*eta; etaEnergySum[ijet] += absorbedE*detaBuf[ijet]*detaBuf[ijet]; phiEnergySum[ijet] += absorbedE*dphiBuf[ijet]*dphiBuf[ijet]; etaPhiESum[ijet] += absorbedE*detaBuf[ijet]*dphiBuf[ijet]; phiSum[ijet] += absorbedE*dphiBuf[ijet]; // // How does this affect the energy variance? // No easy way to make a meaningful estimate. // energyVariance[ijet] += 0.0; } if (cellE > 0.0) { // There is still some energy left to // assign it to the noise/background if (maxE[0] > 0.0) // Some energy was assigned to jets *unclustered += vMaker(cellE, eta, phi); else // No energy was assigned to jets. // No need to recalculate the 4-vector. *unclustered += gridVec; *unused += cellE; } } else { // Just assign the whole cell to the unclustered // energy *unclustered += gridVec; *unused += gridValue; } } else { const double totalSum = sum + bgW; if (buildClusterMask && maxW > bgW) clusterMask[ipt] = maxWJet + 1; if (winnerTakesAll) { // Hard clustering. All energy goes // to the cluster with the highest weight. if (maxW > bgW) { if (use4Vec) jets[maxWJet] += gridVec; weightSum[maxWJet] += 1.0; energySum[maxWJet] += gridValue; + entropySum[maxWJet] += gridValue*log(gridValue); + doubleWeightedEnergy[maxWJet] += gridValue*maxW/totalSum; etaSum[maxWJet] += gridValue*eta; etaEnergySum[maxWJet] += gridValue*detaBuf[maxWJet]*detaBuf[maxWJet]; phiEnergySum[maxWJet] += gridValue*dphiBuf[maxWJet]*dphiBuf[maxWJet]; etaPhiESum[maxWJet] += gridValue*detaBuf[maxWJet]*dphiBuf[maxWJet]; energyVariance[maxWJet] += gridValue*gridValue*(1.0 - maxW/totalSum); phiSum[maxWJet] += gridValue*dphiBuf[maxWJet]; } else { *unclustered += gridVec; *unused += gridValue; } } else { // Fuzzy clustering. Distribute the cell energy // among all clusters according to their weights. // Unclustered energy fraction const double bgfrac = bgW/totalSum; if (bgfrac > 0.0) { *unclustered += gridVec*bgfrac; *unused += gridValue*bgfrac; } // Fractions which belong to each jet for (unsigned ijet=0; ijet 0.0) { const double ePart = gridValue*w; if (use4Vec) jets[ijet] += gridVec*w; weightSum[ijet] += w; energySum[ijet] += ePart; + entropySum[ijet] += ePart*log(ePart); + doubleWeightedEnergy[ijet] += ePart*w; etaEnergySum[ijet] += ePart*detaBuf[ijet]*detaBuf[ijet]; etaSum[ijet] += ePart*eta; phiEnergySum[ijet] += ePart*dphiBuf[ijet]*dphiBuf[ijet]; etaPhiESum[ijet] += ePart*detaBuf[ijet]*dphiBuf[ijet]; energyVariance[ijet] += gridValue*ePart*(1.0 - w); phiSum[ijet] += ePart*dphiBuf[ijet]; } } } } } buildOutput(outputJets); return 0; } } Index: trunk/fftjet/PeakFinder.icc =================================================================== --- trunk/fftjet/PeakFinder.icc (revision 44) +++ trunk/fftjet/PeakFinder.icc (revision 45) @@ -1,295 +1,296 @@ +#include #include namespace fftjet { inline void PeakFinder::setPeakHeightCutoff(const double cutoff) { peakHeightCutoff_ = cutoff; } // The function below works by expanding the grid // values into truncated orthogonal polynomial series template bool PeakFinder::find3by3(const Real g[3][3], Real* x, Real* y, double hessian[3]) { const Real xm1(g[0][0] + g[0][1] + g[0][2]); const Real x0(g[1][0] + g[1][1] + g[1][2]); const Real xp1(g[2][0] + g[2][1] + g[2][2]); const Real ym1(g[0][0] + g[1][0] + g[2][0]); const Real y0(g[0][1] + g[1][1] + g[2][1]); const Real yp1(g[0][2] + g[1][2] + g[2][2]); const Real cx((xp1 - xm1)/6); const Real cy((yp1 - ym1)/6); const Real cxy((g[0][0] - g[0][2] - g[2][0] + g[2][2])/4); const Real cxsq((xm1 + xp1 - 2*x0)/3); const Real cysq((ym1 + yp1 - 2*y0)/3); const Real det(cxsq*cysq - cxy*cxy); if (det) { *x = (cxy*cy - cx*cysq)/det; *y = (cx*cxy - cxsq*cy)/det; } else { *x = 0; *y = 0; } hessian[0] = cxsq; hessian[1] = cxy; hessian[2] = cysq; return cxsq <= 0 && cysq <= 0 && det >= 0; } // The meaning of filled mask values is as follows: // 0 -- Definitely not a peak (there is a higher value // neighbor or the value is not above the cutoff). // 1 -- Plateau: all neighbors have the same values. // 2 -- Peak seed: there are no higher neighbors and // there is at least one lower neighbor // 3 -- Definite peak: all neighbors are lower. // 4 -- This value can be set only by the mask // processing function. It means that the // peak position is exactly at that point, // and that subcell resolution should not // be applied. This happens, for example, // in the rare case of a 3x3 plateau. An attempt // to determine the resolution with subcell // precision would in this case result in // a near-zero determinant and a very unstable // solution. template void PeakFinder::fillMask(const Real* data, const unsigned nEta, const unsigned nPhi, const unsigned minEtaBin, const unsigned maxEtaBin, unsigned *nseeds, unsigned *npeaks) { clearMaskBuffer(nEta, nPhi); *nseeds = 0; *npeaks = 0; for (unsigned ieta = minEtaBin; ieta < maxEtaBin; ++ieta) { const Real* etam1 = data + (ieta - 1)*nPhi; const Real* etadata = data + ieta*nPhi; const Real* etap1 = data + (ieta + 1)*nPhi; int* m = mask + ieta*nPhi; for (unsigned iphi = 0; iphi < nPhi; ++iphi) { const Real value(etadata[iphi]); if (value <= peakHeightCutoff_) continue; const unsigned iphim1 = iphi ? iphi - 1 : nPhi - 1; const unsigned iphip1 = (iphi + 1) % nPhi; if (value >= etadata[iphim1] && value >= etadata[iphip1] && value >= etam1[iphim1] && value >= etam1[iphi] && value >= etam1[iphip1] && value >= etap1[iphim1] && value >= etap1[iphi] && value >= etap1[iphip1]) { m[iphi] = 1; if (value > etadata[iphim1] || value > etadata[iphip1] || value > etam1[iphim1] || value > etam1[iphi] || value > etam1[iphip1] || value > etap1[iphim1] || value > etap1[iphi] || value > etap1[iphip1]) { m[iphi] = 2; ++*nseeds; if (value > etadata[iphim1] && value > etadata[iphip1] && value > etam1[iphim1] && value > etam1[iphi] && value > etam1[iphip1] && value > etap1[iphim1] && value > etap1[iphi] && value > etap1[iphip1]) { m[iphi] = 3; ++*npeaks; } } } } } } template void PeakFinder::print3by3(const Real g[3][3]) const { std::cout << '{'; for (unsigned i=0; i<3; ++i) { if (i) std::cout << ' '; std::cout << '{'; for (unsigned j=0; j<3; ++j) { if (j) std::cout << ' '; std::cout << g[i][j]; } std::cout << '}'; } std::cout << '}' << std::endl; } template int PeakFinder::find(const Real* data, const unsigned nEta, const unsigned nPhi, std::vector* result) { assert(data); assert(nEta); assert(nPhi); assert(result); int status = 0; result->clear(); const unsigned minEtaBin = minEtaBin_ ? minEtaBin_ : 1; const unsigned maxEtaBin = maxEtaBin_ > nEta-1 ? nEta-1 : maxEtaBin_; // Take care of the mask memory allocateMaskBuffer(nEta, nPhi); // Mask potential peak positions unsigned nMaskSeeds, nMaskPeaks; fillMask(data, nEta, nPhi, minEtaBin, maxEtaBin, &nMaskSeeds, &nMaskPeaks); if (nMaskSeeds == 0) return status; if (nMaskSeeds != nMaskPeaks) // Perform some pattern recognition on the mask processMask(nEta, nPhi, minEtaBin, maxEtaBin); // Fill out the vector of results double hessian[3] = {0.0, 0.0, 0.0}; if (subCellResolution_) { Real peakEta, peakPhi, g[3][3]; for (unsigned ieta = minEtaBin; ieta < maxEtaBin; ++ieta) { const int* m = mask + ieta*nPhi; for (unsigned iphi = 0; iphi < nPhi; ++iphi) { if (m[iphi] == 4) { hessian[0] = 0.0; hessian[1] = 0.0; hessian[2] = 0.0; result->push_back( Peak(ieta, iphi, data[ieta*nPhi + iphi], hessian)); } else if (m[iphi] == 3) { const Real* etam1 = data + (ieta - 1)*nPhi; const Real* etadata = data + ieta*nPhi; const Real* etap1 = data + (ieta + 1)*nPhi; const unsigned iphim1 = iphi ? iphi - 1 : nPhi - 1; const unsigned iphip1 = (iphi + 1) % nPhi; g[0][0] = etam1[iphim1]; g[0][1] = etam1[iphi]; g[0][2] = etam1[iphip1]; g[1][0] = etadata[iphim1]; g[1][1] = etadata[iphi]; g[1][2] = etadata[iphip1]; g[2][0] = etap1[iphim1]; g[2][1] = etap1[iphi]; g[2][2] = etap1[iphip1]; if (find3by3(g, &peakEta, &peakPhi, hessian)) { // We should keep it between -1 and 1. // Otherwise we get into various troubles, // incuding the problem that we could // obtain a negative bin number for // the peak position in eta. if (peakEta >= -1.0 && peakEta <= 1.0 && peakPhi >= -1.0 && peakPhi <= 1.0) { result->push_back(Peak( static_cast(ieta) + peakEta, static_cast(iphi) + peakPhi, etadata[iphi], hessian)); } else { ++status; // The fit result is not likely // to be reliable for extrapolation. // Append the coordinates of the // center cell instead. result->push_back( Peak(ieta, iphi, etadata[iphi], hessian)); if (printFitWarnings_) { // Print out a message about the problem std::cout << "WARNING! In fftjet::PeakFinder : " << "fitted peak position is out of " << "reticle." << std::endl; std::cout << "ieta = " << ieta << ", iphi = " << iphi << ", peak is at {" << peakEta << ' ' << peakPhi << '}' << " from this data:" << std::endl; print3by3(g); } } } else { ++status; // The fit tells us that this is not // a real peak. We will follow this // suggestion and ignore the candidate. if (printFitWarnings_) { // Print out a message about the problem std::cout << "WARNING! In fftjet::PeakFinder : " << "peak candidate rejected." << std::endl; std::cout << "ieta = " << ieta << ", iphi = " << iphi << ", candidate is at {" << peakEta << ' ' << peakPhi << '}' << " from this data:" << std::endl; print3by3(g); } } } } } } else { for (unsigned ieta = minEtaBin; ieta < maxEtaBin; ++ieta) { const int* m = mask + ieta*nPhi; for (unsigned iphi = 0; iphi < nPhi; ++iphi) if (m[iphi] >= 3) result->push_back( Peak(ieta, iphi, data[ieta*nPhi + iphi], hessian)); } } return status; } } Index: trunk/fftjet/KernelVectorRecombinationAlg.hh =================================================================== --- trunk/fftjet/KernelVectorRecombinationAlg.hh (revision 44) +++ trunk/fftjet/KernelVectorRecombinationAlg.hh (revision 45) @@ -1,226 +1,228 @@ #ifndef FFTJET_KERNELVECTORRECOMBINATIONALG_HH_ #define FFTJET_KERNELVECTORRECOMBINATIONALG_HH_ //========================================================================= // KernelVectorRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // // It is not intended for this class to be constructed and destroyed // often -- it does too many allocations/deallocations of memory // buffers to work efficiently in this mode. Instead, create one // instance of this class at the beginning of your event processing // loop and call the "run" function for each event. // // The "VBuilder" functor on which this class is templated must // implement a method with the following prototype: // // VectorLike operator()(double energyLike, double eta, double phi) const; // // This method builds VectorLike objects (e.g., 4-vectors). There is no // abstract base class for VBuilder because it is used inside a pretty // tight loop where execution speed is important. // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // I. Volobouev // June 2009 //========================================================================= #include #include "fftjet/AbsVectorRecombinationAlg.hh" #include "fftjet/AbsKernel2d.hh" #include "fftjet/AbsMembershipFunction.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template class KernelVectorRecombinationAlg : public AbsVectorRecombinationAlg { public: typedef double (VectorLike::* VectorLikeMemberFunction)() const; // The meaning of the constructor arguments is as follows: // // kernel -- The default kernel functor used to represent // the jet membership function centered at each // peak position. This default functor will be // used whenever the functor associated with a Peak // object is not set. Here, we can use an instance // of any class derived either from AbsKernel2d or // from AbsMembershipFunction. The default functor // can be NULL in which case all individual Peak // functors must be provided. // // etFcn, -- Member functions of the "VectorLike" class // etaFcn, which return Et (or pt), eta, and phi for // phiFcn the input data. // // bgWeight -- Background/noise membership functor. Must // implement the operator // // "double operator()(const double& et, // const BgData& bg) const" // // which returns the function value. "et" argument // will be set to the result of the "etFcn" member // function call. "bg" argument will be set to the // corresponding background description. // // unlikelyBgWeight -- If the membership function values // for each jet are 0 and the background // weight is smaller than this parameter, // the cell is a poor fit to the probability // model used. For membership functions with // explicit cell Et argument, this condition // triggers an alternative handling of cell // energy: the algorithm attempts to split // the energy between several sources. // // winnerTakesAll -- If "true", there will be one-to-one // mapping of grid cells to clustered jets or // to background. That is, each cell will be // assigned to one jet only ("crisp" clustering). // If "false", each grid cell will be assigned // to multiple jets and background using weights // ("fuzzy" clustering). // // buildCorrelationMatrix -- This parameter is reserved for // future use (currently ignored). // // buildClusterMask -- If "true", the code will remember how // input vectors are assigned to clustered jets // for the last processed event. The assignments // are calculated as if the "winnerTakesAll" // parameter is "true" no matter what its actual // value is. The mask can be later retrieved // using the "getClusterMask" function. // // This class does not assume ownership of the jet and background // membership functors. It is a responsibility of the user of this // class to make sure that the "run" method is called only when // the membership functor objects are still alive. // KernelVectorRecombinationAlg( ScaleSpaceKernel* kernel, VectorLikeMemberFunction etFcn, VectorLikeMemberFunction etaFcn, VectorLikeMemberFunction phiFcn, const Functor2* bgWeight, double unlikelyBgWeight, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask); virtual ~KernelVectorRecombinationAlg(); // In this particular algorithm, the "etaWidth" and "phiWidth" of // the output jets will be estimated with respect to the Et centroid // direction (minimal width), not the jet direction. If you need it // w.r.t. the jet direction, add in quadrature the angular distance // from the jet direction to the centroid. virtual int run(const std::vector& peaks, const std::vector& eventData, const BgData* bgData, unsigned bgDataLength, std::vector >* jets, VectorLike* unclustered, double* unused); virtual unsigned getLastNData() const {return lastNData;} virtual unsigned getLastNJets() const {return passedJets;} virtual const unsigned* getClusterMask() const; protected: // Override the following function to tell the // algorithm whether it should build jets using 4-vectors inline virtual bool recombine4Vectors() {return true;} // Override the following function to tell the // algorithm whether it should build 4-vectors // out of Et centroids (this is only used when // the 4-vectors are not in use). inline virtual bool useEtCentroid() {return true;} // The following function builds the output 4-vectors void buildOutput(std::vector >* outputJets); // If the following function returns "true", we are done // with this event bool performPreliminaryProcessing( const std::vector& peaks, const std::vector& eventData, const BgData* bgData, unsigned bgDataLength, std::vector >* outJets, VectorLike* unclustered, double* unused); // Input parameters ScaleSpaceKernel* const defaultKernel; VectorLikeMemberFunction const etFcn; VectorLikeMemberFunction const etaFcn; VectorLikeMemberFunction const phiFcn; const Functor2* const bgWeightCalc; const double unlikelyBgWeight; const bool winnerTakesAll; const bool buildCorrelationMatrix; const bool buildClusterMask; // Vector builder VBuilder vMaker; // Various useful buffers with the size // dependent on the number of input peaks double* weights; double* kernelScales; double* clusterScales; double* clusterScaleRatios; double* energySum; + double* doubleWeightedEnergy; double* energyVariance; double* weightSum; + double* entropySum; double* etaEnergySum; double* phiEnergySum; double* etaPhiESum; double* etaSum; double* phiSum; double* detaBuf; double* dphiBuf; unsigned* clusterMap; const Peak** peakPositions; AbsKernel2d** memFcns2d; AbsMembershipFunction** memFcns3d; VectorLike* jets; unsigned nWeights; // Cluster mask buffer and its size unsigned* clusterMask; unsigned maskMemSize; // Using 2d or 3d membership functions? bool use3d; // Various variables from the last event unsigned lastNJets; unsigned lastNData; unsigned passedJets; // Background data is a vector? bool bgIsAVector; private: KernelVectorRecombinationAlg(); KernelVectorRecombinationAlg(const KernelVectorRecombinationAlg&); KernelVectorRecombinationAlg& operator=(const KernelVectorRecombinationAlg&); void remapClusterMask(); void setupBuffers(unsigned njets, unsigned ndata); void allUnclustered(const std::vector& eventData, VectorLike* unclus, double* unused); void processKernelScales(const std::vector& peaks); }; } #include "fftjet/KernelVectorRecombinationAlg.icc" #endif // FFTJET_KERNELVECTORRECOMBINATIONALG_HH_ Index: trunk/fftjet/Peak.cc =================================================================== --- trunk/fftjet/Peak.cc (revision 44) +++ trunk/fftjet/Peak.cc (revision 45) @@ -1,125 +1,127 @@ #include #include "fftjet/Peak.hh" namespace fftjet { Peak::Peak(const double eta, const double phi, const double mag, const double hessian[3], const double drsp, const double magsp, const double lifetime, const double scale, const double nearestDistance, const double memScale, const double recoScale, const double recoScaleRatio, const double clusterRadius, const double clusterSeparation, const int code, const int status, ScaleSpaceKernel* membershipFunction) : eta_(eta), phi_(phi), magnitude_(mag), speed_(drsp), magSpeed_(magsp), lifetime_(lifetime), splitTime_(-1.0), mergeTime_(-1.0), scale_(scale), nearestD_(nearestDistance), membershipFactor_(memScale), recoScale_(recoScale), recoScaleRatio_(recoScaleRatio), clusterRadius_(clusterRadius), clusterSeparation_(clusterSeparation), + relHDetDeriv_(-10.0), code_(code), status_(status), memFcn_(membershipFunction) { hessian_[0] = hessian[0]; hessian_[1] = hessian[1]; hessian_[2] = hessian[2]; } void Peak::hessian(double hessianArray[3]) const { hessianArray[0] = hessian_[0]; hessianArray[1] = hessian_[1]; hessianArray[2] = hessian_[2]; } void Peak::setHessian(const double hessianArray[3]) { hessian_[0] = hessianArray[0]; hessian_[1] = hessianArray[1]; hessian_[2] = hessianArray[2]; } double Peak::hessianDeterminant() const { return hessian_[0]*hessian_[2] - hessian_[1]*hessian_[1]; } double Peak::scaledHDeterminant() const { const double s2 = scale_*scale_; return s2*s2*(hessian_[0]*hessian_[2] - hessian_[1]*hessian_[1]); } double Peak::relativeScaledHDeterminant() const { assert(magnitude_ > 0.0); const double tmp = scale_*scale_/magnitude_; return tmp*tmp*(hessian_[0]*hessian_[2] - hessian_[1]*hessian_[1]); } double Peak::scaledLaplacian(const double bwEta, const double bwPhi) const { // Note that both scale_*scale_*hessian_[0] and // scale_*scale_*hessian_[2] should simultaneously peak // at the correct scale. The addition of the scale factors here // simply makes the peak magnitudes comparable (of course, // only if the bandwidth ratio used for clustering correctly // reflects the bandwidth ratio in the data). return scale_*scale_*(hessian_[0]*bwEta*bwEta + hessian_[2]*bwPhi*bwPhi); } double Peak::relativeScaledLaplacian( const double bwEta, const double bwPhi) const { assert(magnitude_ > 0.0); return scale_*scale_/magnitude_*(hessian_[0]*bwEta*bwEta + hessian_[2]*bwPhi*bwPhi); } bool Peak::operator==(const Peak& r) const { const bool hessians_equal = hessian_[0] == r.hessian_[0] && hessian_[1] == r.hessian_[1] && hessian_[2] == r.hessian_[2]; return hessians_equal && eta_ == r.eta_ && phi_ == r.phi_ && magnitude_ == r.magnitude_ && speed_ == r.speed_ && magSpeed_ == r.magSpeed_ && lifetime_ == r.lifetime_ && splitTime_ == r.splitTime_ && mergeTime_ == r.mergeTime_ && scale_ == r.scale_ && nearestD_ == r.nearestD_ && membershipFactor_ == r.membershipFactor_ && recoScale_ == r.recoScale_ && recoScaleRatio_ == r.recoScaleRatio_ && clusterRadius_ == r.clusterRadius_ && clusterSeparation_ == r.clusterSeparation_ && + relHDetDeriv_ == r.relHDetDeriv_ && code_ == r.code_ && status_ == r.status_; // There isn't any good way to compare membership functions // for equality. Comparing their pointers for identity does // not make a lot of sense here. } Peak Peak::dummy() { double hess[3] = {0.0, 0.0, 0.0}; return Peak(0.0, 0.0, 0.0, hess); } } Index: trunk/fftjet/KernelRecombinationAlg.hh =================================================================== --- trunk/fftjet/KernelRecombinationAlg.hh (revision 44) +++ trunk/fftjet/KernelRecombinationAlg.hh (revision 45) @@ -1,253 +1,255 @@ #ifndef FFTJET_KERNELRECOMBINATIONALG_HH_ #define FFTJET_KERNELRECOMBINATIONALG_HH_ //========================================================================= // KernelRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // // It is not intended for this class to be constructed and destroyed // often -- it does too many allocations/deallocations of memory // buffers to work efficiently in this mode. Instead, create one // instance of this class at the beginning of your event processing // loop and call the "run" function for each event. // // The "VBuilder" functor on which this class is templated must // implement a method with the following prototype: // // VectorLike operator()(Real energyLike, Real eta, Real phi) const; // // This method builds VectorLike objects (e.g., 4-vectors) out of grid // points. These objects are later agglomerated by the recombination // algorithm. There is no abstract base class for VBuilder because it is // used inside a pretty tight loop where execution speed is important. // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // I. Volobouev // March 2008 //========================================================================= #include #include "fftjet/AbsRecombinationAlg.hh" #include "fftjet/AbsKernel2d.hh" #include "fftjet/AbsMembershipFunction.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class KernelRecombinationAlg : public AbsRecombinationAlg { public: // The meaning of the constructor arguments is as follows: // // kernel -- The default kernel functor used to represent // the jet membership function centered at each // peak position. This default functor will be // used whenever the functor associated with a Peak // object is not set. Here, we can use an instance // of any class derived either from AbsKernel2d or // from AbsMembershipFunction. The default functor // can be NULL in which case all individual Peak // functors must be provided. // // bgWeight -- Background/noise membership functor. Must // implement the operator // // "double operator()(const double& et, // const BgData& bg) const" // // which returns the function value. "et" argument // will be set to the cell value from the event // energy discretization grid. "bg" argument will be // set to the corresponding background description. // // unlikelyBgWeight -- If the membership function values // for each jet are 0 and the background // weight is smaller than this parameter, // the cell is a poor fit to the probability // model used. For membership functions with // explicit cell Et argument, this condition // triggers an alternative handling of cell // energy: the algorithm attempts to split // the energy between several sources. // // dataCutoff -- Only the data points with values above // this cutoff will contribute to the final // clusters. The cutoff is not intended // for noise suppression. Instead, it // should be used to skip a large number of // zero-level points in the data grid which // are present in case the data is sparse. // In this case the cutoff should be set to 0.0. // If the data is not sparse, it is best to set // this cutoff to a negative number of large // magnitude. // // winnerTakesAll -- If "true", there will be one-to-one // mapping of grid cells to clustered jets or // to background. That is, each cell will be // assigned to one jet only ("crisp" clustering). // If "false", each grid cell will be assigned // to multiple jets and background using weights // ("fuzzy" clustering). // // buildCorrelationMatrix -- This parameter is reserved for // future use (currently ignored). // // buildClusterMask -- If "true", the code will remember how // grid cells are assigned to clustered jets // for the last processed event. The assignments // are calculated as if the "winnerTakesAll" // parameter is "true" no matter what its actual // value is. The mask can be later retrieved // using the "getClusterMask" function. // // etaBinMin -- The minimum grid bin number in eta. The code // will not attempt to cluster the cells below // that bin number. // // etaBinMax -- The maximum grid bin number in eta. The code // will not attempt to cluster the cells for that // bin number or larger. If this parameter is // larger than the grid size then the grid size // will be used internally. // // This class does not assume ownership of the jet and background // membership functors. It is a responsibility of the user of this // class to make sure that the "run" method is called only when // the membership functor objects are still alive. // KernelRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX); virtual ~KernelRecombinationAlg(); // In this particular algorithm, the "etaWidth" and "phiWidth" of // the output jets will be estimated with respect to the Et centroid // direction (minimal width), not the jet direction. If you need it // w.r.t. the jet direction, add in quadrature the angular distance // from the jet direction to the centroid. virtual int run(const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, unsigned nBgEta, unsigned nBgPhi, std::vector >* outputJets, VectorLike* unclustered, double* unused); virtual inline unsigned getLastNEta() const {return lastNEta;} virtual inline unsigned getLastNPhi() const {return lastNPhi;} virtual inline unsigned getLastNJets() const {return passedJets;} virtual const unsigned* getClusterMask() const; protected: // Override the following function to tell the // algorithm whether it should build jets using 4-vectors inline virtual bool recombine4Vectors() {return true;} // Override the following function to tell the // algorithm whether it should build 4-vectors // out of Et centroids (this is only used when // the 4-vectors are not in use). inline virtual bool useEtCentroid() {return true;} // The following function builds the output 4-vectors void buildOutput(std::vector >* outputJets, bool setRecoScaleRatiosToZero); // If the following function returns "true", we are done // with this event bool performPreliminaryProcessing( const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, unsigned nBgEta, unsigned nBgPhi, std::vector >* outJets, VectorLike* unclustered, double* unused); // Input parameters ScaleSpaceKernel* const defaultKernel; const Functor2* const bgWeightCalc; const double unlikelyBgWeight; const Real dataCutoff; const bool winnerTakesAll; const bool buildCorrelationMatrix; const bool buildClusterMask; const unsigned etaBinMin; const unsigned etaBinMax0; // Vector builder VBuilder vMaker; // Various useful buffers with the size // dependent on the number of input peaks double* weights; double* kernelScales; double* clusterScales; double* clusterScaleRatios; double* dEta; double* energySum; + double* doubleWeightedEnergy; double* energyVariance; double* weightSum; + double* entropySum; double* etaEnergySum; double* phiEnergySum; double* etaPhiESum; double* etaSum; double* phiSum; unsigned* clusterMap; const Peak** peakPositions; AbsKernel2d** memFcns2d; AbsMembershipFunction** memFcns3d; VectorLike* jets; unsigned nWeights; // Cluster mask buffer and its size unsigned* clusterMask; unsigned maskMemSize; // Buffer for dphi values and its size double* dphiBuffer; unsigned dphiBufSize; // Using 2d or 3d membership functions? bool use3d; // Various variables from the last event unsigned lastNJets; unsigned lastNEta; unsigned lastNPhi; unsigned passedJets; unsigned bgDim; private: KernelRecombinationAlg(); KernelRecombinationAlg(const KernelRecombinationAlg&); KernelRecombinationAlg& operator=(const KernelRecombinationAlg&); void remapClusterMask(); void setupBuffers(unsigned njets, unsigned nEta, unsigned nPhi); void calculateDphi(const Grid2d& eventData); void allUnclustered(const Grid2d& evData, VectorLike* unclus, double* unused); void processKernelScales(const std::vector& peaks); }; } #include "fftjet/KernelRecombinationAlg.icc" #endif // FFTJET_KERNELRECOMBINATIONALG_HH_