Index: trunk/fftjet/scanFrequencyKernel.hh =================================================================== --- trunk/fftjet/scanFrequencyKernel.hh (revision 43) +++ trunk/fftjet/scanFrequencyKernel.hh (revision 44) @@ -1,118 +1,118 @@ +#ifndef FFTJET_SCANFREQUENCYKERNEL_HH_ +#define FFTJET_SCANFREQUENCYKERNEL_HH_ + //======================================================================== // scanFrequencyKernel.hh // // A standard way to scan frequency kernels // // I. Volobouev // May 2009 //======================================================================== -#ifndef FFTJET_SCANFREQUENCYKERNEL_HH_ -#define FFTJET_SCANFREQUENCYKERNEL_HH_ - #include #include "fftjet/AbsFFTEngine.hh" #include "fftjet/DiscreteGauss1d.hh" #include "fftjet/DiscreteGauss2d.hh" namespace fftjet { template void scanFrequencyKernel( const AbsFFTEngine* const engine, const AbsFrequencyKernel* kernel, const double scale, Complex* to) { assert(engine); assert(kernel); assert(to); // The "DiscreteGauss2d" kernel is too important // to allow for random screw-ups. Force binning // compatibility here. { const DiscreteGauss2d* dg = dynamic_cast(kernel); if (dg) assert(dg->isCompatible(engine->nEta(), engine->nPhi())); } const int nEta = static_cast(engine->nEta()); const int nPhi = static_cast(engine->nPhi()); for (int ieta=0; ietasetTransformPoint(to, ieta, iphi, (*kernel)(ie, ip, scale)); } } } template void scanFrequencyKernel1d( const AbsFFTEngine* const engine, const AbsFrequencyKernel1d* kernel, const double scale, Complex* to) { assert(engine); assert(kernel); assert(to); // The "DiscreteGauss1d" kernel is too important // to allow for random screw-ups. Force binning // compatibility here. { const DiscreteGauss1d* dg = dynamic_cast(kernel); if (dg) assert(dg->nPhi() == engine->nPhi()); } const int nPhi = static_cast(engine->nPhi()); for (int iphi=0; iphisetTransformPoint(to, 0, iphi, (*kernel)(ip, scale)); } } // The following function can be used to normalize // the Fourier image, so that its power corresponds // to the power of the given spatial image treated // as a density template Real normalizeFrequencyAsDensity( const AbsFFTEngine* const engine, Complex* imageBuf, const Real* timeRep) { assert(engine); assert(imageBuf); assert(timeRep); const unsigned nEta = engine->nEta(); const unsigned nPhi = engine->nPhi(); const unsigned nbins = nEta*nPhi; long double sum = 0.0L, sumsq = 0.0L; for (unsigned i=0; itotalPower(imageBuf); assert(imagePower > 0.0); // The image power should be equal to timePower*nEta*nPhi engine->scaleTransform(imageBuf, sqrt(timePower*nbins/imagePower), imageBuf); return static_cast(sum); } } #endif // FFTJET_SCANFREQUENCYKERNEL_HH_ Index: trunk/fftjet/PeakFinder.hh =================================================================== --- trunk/fftjet/PeakFinder.hh (revision 43) +++ trunk/fftjet/PeakFinder.hh (revision 44) @@ -1,150 +1,150 @@ +#ifndef FFTJET_PEAKFINDER_HH_ +#define FFTJET_PEAKFINDER_HH_ + //========================================================================= // PeakFinder.hh // // This class implements an algorithm for locating precluster positions // // I. Volobouev // March 2008 //========================================================================= -#ifndef FFTJET_PEAKFINDER_HH_ -#define FFTJET_PEAKFINDER_HH_ - #include #include #include "fftjet/Peak.hh" namespace fftjet { class PeakFinder { public: // // The "peakHeightCutoff" parameter of the constructor should be // used in case large plateaus are expected at the value of // "peakHeightCutoff" or just below. A typical use for this // parameter would be to raise the peak threshold above the // round-off noise of the FFT algorithm. If such plateaus are not // expected, it is best to set the "peakHeightCutoff" parameter to // some negative number of large magnitude. // // Other constructor arguments have the following meaning: // // subCellResolution -- Set this to "true" in order to improve // the peak position determination by fitting // the 3x3 region with the peak in the center // using a 2-d quadratic polynomial (and then // finding the peak of that polynomial). // // minEtaBin -- Minimum eta bin number for peak searching. // Meaningful values of this parameter are 1 and higher. // Value of 0 will be internally turned into 1. // // maxEtaBin -- Maximum eta bin number for peak searching. // Meaningful values of this parameter are nEta-1 // and lower, where "nEta" is the size of the grid // with which the "find" function will be called. // // printFitWarnings -- If "true", the peak finder will be verbose // about problems encountered during the peak // position fitting. This parameter is // meaningful only when the "subCellResolution" // value is "true". // explicit PeakFinder(double peakHeightCutoff, bool subCellResolution=true, unsigned minEtaBin=1, unsigned maxEtaBin=UINT_MAX, bool printFitWarnings=false); PeakFinder(const PeakFinder&); ~PeakFinder(); // Change the cutoff for the peak height void setPeakHeightCutoff(double cutoff); // Trivial inspectors of the contents inline double peakHeightCutoff() const {return peakHeightCutoff_;} inline bool subCellResolution() const {return subCellResolution_;} inline unsigned smallestEtaBin() const {return minEtaBin_;} inline unsigned largestEtaBin() const {return maxEtaBin_;} // The following function will not look for peaks at // etaBin = 0 and etaBin = nEta-1. The peaks are searched // for only inside the grid, not on the boundaries. // The searched range can be limited even more (and the // code speed improved) by making the eta range smaller // with the "minEtaBin" and "maxEtaBin" parameters // in the constructor. // // It is assumed that the phi coordinate is periodic // (it has no boundary). The "data" array should be laid // out so that the fastest changing coordinate is phi. // // The peak coordinates (eta and phi) are returned // in the units of grid cell size. For example, if the peak // is in the middle of the input cell data[5*nPhi + 7] then // the returned peak coordinate will be (5, 7). // // The function returns 0 in case everything is OK. When peaks // with subcell resultion are requested, it returns the number of // "unstable" peaks -- that is, peaks whose subcell position are // outside of the initial highest cell, the peaks for which the // fitted second derivatives are not negative, etc. Significant // number of such peaks usually indicates that the round-off errors // are too large and that the value of the "peakHeightCutoff" // constructor parameter should be increased. // template int find(const Real* data, unsigned nEta, unsigned nPhi, std::vector* result); // Helper function for fitting the peak position in a more precise // manner when curve values are known on a 3x3 grid. The values // are fitted to a quadratic polynomial in 2d and the extremum of // that polynomial is found. This code assumes that the grid // locations are at -1, 0, and 1 in both directions. Correct // shifting and scaling is up to the user of this function. // // The function returns "true" if the extremum found is actually // a maximum, "false" otherwise. The extremum coordinates are // placed into the locations pointed to by x and y. // template static bool find3by3(const Real gridValues[3][3], Real* x, Real* y, double hessian[3]); private: PeakFinder(); PeakFinder& operator=(const PeakFinder&); void allocateMaskBuffer(unsigned nEta, unsigned nPhi); void clearMaskBuffer(unsigned nEta, unsigned nPhi); template void fillMask(const Real* data, unsigned nEta, unsigned nPhi, unsigned minEtaBin, unsigned maxEtaBin, unsigned *nseeds, unsigned *npeaks); void processMask(unsigned nEta, unsigned nPhi, unsigned minEtaBin, unsigned maxEtaBin); bool match3by3(unsigned iEta, unsigned iPhi, unsigned nPhi, const int pattern[3][3]); void tag3by3(unsigned iEta, unsigned iPhi, unsigned nPhi, int tagValue); template void print3by3(const Real gridValues[3][3]) const; double peakHeightCutoff_; const unsigned minEtaBin_; const unsigned maxEtaBin_; int *mask; unsigned nbins; const bool subCellResolution_; const bool printFitWarnings_; static const int plateau3x3[3][3]; static const int ridges3x3[4][3][3]; }; } #include "fftjet/PeakFinder.icc" #endif // FFTJET_PEAKFINDER_HH_ Index: trunk/fftjet/SimpleFunctors.hh =================================================================== --- trunk/fftjet/SimpleFunctors.hh (revision 43) +++ trunk/fftjet/SimpleFunctors.hh (revision 44) @@ -1,72 +1,72 @@ +#ifndef FFTJET_SIMPLEFUNCTORS_HH_ +#define FFTJET_SIMPLEFUNCTORS_HH_ + //========================================================================= // SimpleFunctors.hh // // Base classes for a variety of functor-based calculations // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_SIMPLEFUNCTORS_HH_ -#define FFTJET_SIMPLEFUNCTORS_HH_ - namespace fftjet { template struct Functor1 { typedef Arg1 argument_type; typedef Result result_type; virtual ~Functor1() {} virtual Result operator()(const Arg1&) const = 0; }; template struct Functor2 { typedef Arg1 first_argument_type; typedef Arg2 second_argument_type; typedef Result result_type; virtual ~Functor2() {} virtual Result operator()(const Arg1&, const Arg2&) const = 0; }; // A simple functor which returns its argument template struct Same : public Functor1 { inline Result operator()(const Result& a) const {return a;} }; // Adaptation for using single-argument functors with simple // cmath-style functions. Do not use this struct as a base class. template struct FcnFunctor1 : public Functor1 { inline explicit FcnFunctor1(Result (*fcn)(Arg1)) : fcn_(fcn) {} inline Result operator()(const Arg1& a) const {return fcn_(a);} private: FcnFunctor1(); Result (*fcn_)(Arg1); }; // Adaptation for using two-argument functors with simple // cmath-style functions. Do not use this struct as a base class. template struct FcnFunctor2 : public Functor2 { inline explicit FcnFunctor2(Result (*fcn)(Arg1, Arg2)) : fcn_(fcn) {} inline Result operator()(const Arg1& x, const Arg2& y) const {return fcn_(x, y);} private: FcnFunctor2(); Result (*fcn_)(Arg1, Arg2); }; } #endif // FFTJET_SIMPLEFUNCTORS_HH_ Index: trunk/fftjet/AbsScalableKernel1d.hh =================================================================== --- trunk/fftjet/AbsScalableKernel1d.hh (revision 43) +++ trunk/fftjet/AbsScalableKernel1d.hh (revision 44) @@ -1,117 +1,117 @@ +#ifndef FFTJET_ABSSCALABLEKERNEL1D_HH_ +#define FFTJET_ABSSCALABLEKERNEL1D_HH_ + //========================================================================= // AbsScalableKernel1d.hh // // Interface class for scalable 1d kernel functions. Derive your kernel // from this class if only the width of your kernel depends on scale but // not the shape. // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_ABSSCALABLEKERNEL1D_HH_ -#define FFTJET_ABSSCALABLEKERNEL1D_HH_ - #include #include #include "fftjet/AbsKernel1d.hh" namespace fftjet { class AbsScalableKernel1d : public AbsKernel1d { public: // xScaleFactor converts the scale into bandwidth value: // x_bandwidth = xScaleFactor*pow(scale, scalePower), inline AbsScalableKernel1d(const double xScaleFactor, const int scalePower) : xs_(xScaleFactor), scalePower_(scalePower) { assert(xs_); } virtual ~AbsScalableKernel1d() {} // Trivial accessors inline double xScaleFactor() const {return xs_;} inline int scalePower() const {return scalePower_;} // The function value inline double operator()(const double x, const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwx(xs_*pscale); return calculate(x/bwx)/bwx; } // The support interval inline KernelSupportInterval supportInterval(const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } return unscaledInterval() * (xs_*pscale); } // The randomizer inline double random(const double r1, const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } return xs_*pscale*unscaledRandom(r1); } private: AbsScalableKernel1d(); virtual double calculate(double x) const = 0; virtual KernelSupportInterval unscaledInterval() const = 0; virtual double unscaledRandom(double r1) const = 0; const double xs_; const int scalePower_; }; } #endif // FFTJET_ABSSCALABLEKERNEL1D_HH_ Index: trunk/fftjet/FasterKernelRecombinationAlg.hh =================================================================== --- trunk/fftjet/FasterKernelRecombinationAlg.hh (revision 43) +++ trunk/fftjet/FasterKernelRecombinationAlg.hh (revision 44) @@ -1,186 +1,186 @@ +#ifndef FFTJET_FASTERKERNELRECOMBINATIONALG_HH_ +#define FFTJET_FASTERKERNELRECOMBINATIONALG_HH_ + //========================================================================= // FasterKernelRecombinationAlg.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. // // This code usually runs faster than "KernelRecombinationAlg" because // it maintains a collection of kernel scans internally. Because of this, // it does not have to reevaluate the membership function for each // jet, instead it looks it up in a table of precomputed values. This // internal cacheing results in a useful acceleration if the kernel // function is expensive to evaluate (for simple kernel functions // the timing of the algorithm is dominated by other code). The // catch is that the precluster positions are assumed to fall in the // centers of the grid cells. This results in a small additional // uncertainty due to binning (which can often be ignored). // // Note that this recombination algorithm can not meaningfully use // kernels whose eta to phi bandwidth ratio is supposed to change // from one peak to another. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_FASTERKERNELRECOMBINATIONALG_HH_ -#define FFTJET_FASTERKERNELRECOMBINATIONALG_HH_ - #include "fftjet/KernelRecombinationAlg.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class FasterKernelRecombinationAlg : public KernelRecombinationAlg { public: // The meaning of the constructor arguments is as follows: // // kernel -- The kernel functor used to represent the cluster // membership function centered at each peak position. // The functor associated with each peak must not // be set, and the kernel itself must be derived // from "AbsKernel2d". // // bgWeight -- Background/pileup membership functor. Must // implement the operator // "double operator()(const double& et, // const BgData& bg) const" // which calculates the background weight. // "et" argument will be set to the cell value // from the event data grid. "bg" argument will be // set to the corresponding background description. // // unlikelyBgWeight -- Reserved for future use (currently ignored) // // dataCutoff -- Only the data points with values above // this cutoff will contribute to the final // clusters. The cutoff is not intended // for background 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/background using weights. // // 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 kernel object. // It is a responsibility of the user of this class to make // sure that the "run" method is called only when the kernel // object is still alive. // FasterKernelRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX); virtual ~FasterKernelRecombinationAlg(); virtual int run(const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, unsigned nBgEta, unsigned nBgPhi, std::vector >* outJets, VectorLike* unclustered, double* unused); private: typedef KernelRecombinationAlg B; // Nested class to manage kernel scan collections class KernelScanCollection { public: explicit KernelScanCollection(const AbsKernel2d* kernel); virtual ~KernelScanCollection(); void mapCoords(const Grid2d& eventData, unsigned etaBinMin, unsigned etaBinMax, const Peak** peakPositions, unsigned njets, unsigned **etaBuf, unsigned **phiBuf); const Real* getKernelScan(const Grid2d& eventData, double scale); private: KernelScanCollection(); KernelScanCollection(const KernelScanCollection&); KernelScanCollection& operator=(const KernelScanCollection&); const AbsKernel2d* const kernel; unsigned *mappedEta; unsigned mapBufLen; std::map*> kernelScans; }; FasterKernelRecombinationAlg(); FasterKernelRecombinationAlg(const FasterKernelRecombinationAlg&); FasterKernelRecombinationAlg& operator=( const FasterKernelRecombinationAlg&); void setupTableBuf(unsigned njets); KernelScanCollection kernelScans; const Real** scanTable; unsigned scanTableLen; }; } #include "fftjet/FasterKernelRecombinationAlg.icc" #endif // FFTJET_FASTERKERNELRECOMBINATIONALG_HH_ Index: trunk/fftjet/JetMagnitudeMapper.hh =================================================================== --- trunk/fftjet/JetMagnitudeMapper.hh (revision 43) +++ trunk/fftjet/JetMagnitudeMapper.hh (revision 44) @@ -1,59 +1,59 @@ +#ifndef FFTJET_JETMAGNITUDEMAPPER_HH_ +#define FFTJET_JETMAGNITUDEMAPPER_HH_ + //========================================================================= // JetMagnitudeMapper.hh // // Simple jet corrector used to invert a jet response curve // which is independent from either location or scale // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_JETMAGNITUDEMAPPER_HH_ -#define FFTJET_JETMAGNITUDEMAPPER_HH_ - #include "fftjet/AbsJetCorrector.hh" #include "fftjet/LinearInterpolator1d.hh" namespace fftjet { template class JetMagnitudeMapper : public AbsJetCorrector { public: // The constructor arguments are as follows: // // f is the jet response curve (ratio of // reconstructed/actual jet magnitude). Should // define "double operator()(double) const". // // maxMagnitude is the maximum value of _actual_ // magnitude for which the jet response // curve is build. Beyond this magnitude // the jet response is assumed to be constant // and equal to its value at maxMagnitude // // npoints number of points to use for modeling // the inverse curve between 0 and the // reconstructed value which corresponds // to maxMagnitude // template JetMagnitudeMapper(const Functor& f, double maxMagnitude, unsigned npoints); inline virtual ~JetMagnitudeMapper() {delete interp;} inline double operator()(const double& magnitude, const Jet&) const {return magnitude*(*interp)(magnitude);} private: JetMagnitudeMapper(); JetMagnitudeMapper(const JetMagnitudeMapper&); JetMagnitudeMapper& operator=(const JetMagnitudeMapper&); LinearInterpolator1d* interp; }; } #include "fftjet/JetMagnitudeMapper.icc" #endif // FFTJET_JETMAGNITUDEMAPPER_HH_ Index: trunk/fftjet/PeakEtaPhiDistance.hh =================================================================== --- trunk/fftjet/PeakEtaPhiDistance.hh (revision 43) +++ trunk/fftjet/PeakEtaPhiDistance.hh (revision 44) @@ -1,43 +1,43 @@ +#ifndef FFTJET_PEAKETAPHIDISTANCE_HH_ +#define FFTJET_PEAKETAPHIDISTANCE_HH_ + //========================================================================= // PeakEtaPhiDistance.hh // // A simple implementation of the AbsDistanceCalculator class for use with // ProximityClusteringTree. Distance is set to the delta R // between peaks in the eta-phi space, adjusted for a non-unit eta/phi // bandwidth ratio if necessary. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_PEAKETAPHIDISTANCE_HH_ -#define FFTJET_PEAKETAPHIDISTANCE_HH_ - #include #include #include "fftjet/AbsDistanceCalculator.hh" #include "fftjet/Peak.hh" namespace fftjet { class PeakEtaPhiDistance : public AbsDistanceCalculator { public: inline explicit PeakEtaPhiDistance(const double etaToPhiBandwidthRatio=1.0) : bwEta(sqrt(etaToPhiBandwidthRatio)), bwPhi(1.0/bwEta) { assert(etaToPhiBandwidthRatio > 0.0); } virtual ~PeakEtaPhiDistance() {} private: double distanceBetweenClusters( double, const Peak& smallerJet, double, const Peak& biggerJet) const; const double bwEta; const double bwPhi; }; } #endif // FFTJET_PEAKETAPHIDISTANCE_HH_ Index: trunk/fftjet/AbsFFTEngine.hh =================================================================== --- trunk/fftjet/AbsFFTEngine.hh (revision 43) +++ trunk/fftjet/AbsFFTEngine.hh (revision 44) @@ -1,125 +1,125 @@ +#ifndef FFTJET_ABSFFTENGINE_HH_ +#define FFTJET_ABSFFTENGINE_HH_ + //========================================================================= // AbsFFTEngine.hh // // Abstract FFT engine for jet reconstruction. User code should normally // avoid dependence on concrete implementations of the engine and use // this header instead. // // I. Volobouev // January 2008 //========================================================================= -#ifndef FFTJET_ABSFFTENGINE_HH_ -#define FFTJET_ABSFFTENGINE_HH_ - #include #include namespace fftjet { // Normally, each derived engine should have a constructor like this: // // ClassName(unsigned nEta, unsigned nPhi); // // so that it can call the base constructor with appropriate arguments. // If the number of eta points specified is 1, the engine should perform // 1-d transforms with "nPhi" points. // // Normally, "Real" type should be one of "float", "double", or // "long double". The "Complex" type should represent complex numbers // at the same level of precision. // template class AbsFFTEngine { public: inline AbsFFTEngine(const unsigned nEta, const unsigned nPhi) : nEta_(nEta), nPhi_(nPhi) { assert(nEta_); assert(nPhi_); } virtual ~AbsFFTEngine() {} // It is assumed that all complex buffers on which member // functions operate are created using "allocateComplex()" // method. virtual void transformForward(const Real* from, Complex* to) const=0; virtual void transformBack(const Complex* from, Real* to) const=0; // Arithmetic operations with complex 2d arrays. Note that // implementations of all of these operations should allow for // "result" to point into the same location as either "a" or "b". virtual void multiplyTransforms(const Complex* a, const Complex* b, Complex* result) const=0; virtual void addTransforms(const Complex* a, const Complex* b, Complex* result) const=0; // The following function calculates a/b virtual void divideTransforms(const Complex* a, const Complex* b, Complex* result) const=0; // The following function calculates a/b, but it does not // complain in case both a and b are complex zeros. In this // case the result is set to complex zero. virtual void deconvolutionRatio(const Complex* a, const Complex* b, Complex* result) const=0; // The following function calculates a - b virtual void subtractTransforms(const Complex* a, const Complex* b, Complex* result) const=0; // The following function multiplies a transform by a real number. // Note that "a" and "result" should be allowed to point at the same // location. virtual void scaleTransform(const Complex* a, Real scale, Complex* result) const=0; // The following function should allow "a" and "result" // to point at the same location virtual void amplitudeSquared(const Complex* a, Complex* result) const=0; // The following function computes the sum of absoule // values squared virtual double totalPower(const Complex* a) const=0; // The following function sets all array elements to Re=value, Im=0 virtual void setToReal(Complex* a, Real value) const=0; // The following function sets all array elements to complex zero. // Override if a faster implementation is available. inline virtual void zeroOut(Complex* a) const { setToReal(a, static_cast(0)); } // The following functions can be used to get/set transform points // using the standard C++ complex class which, in general, will // not be the class used by the engine to represent complex numbers virtual void setTransformPoint(Complex* points, unsigned iEta, unsigned iPhi, const std::complex& p) const=0; virtual std::complex getTransformPoint(const Complex* points, unsigned iEta, unsigned iPhi) const=0; // Allocation functions should generate 2d arrays of proper // size and alignment for the particular engine used. User code // should manage creation and destruction of such arrays by // calling these engine functions. It should be safe to call // "destroyComplex" on a NULL pointer. virtual Complex* allocateComplex() const=0; virtual void destroyComplex(Complex*) const=0; // Dimensionality of the real arrays. Note that the dimensionality // of the complex arrays is likely to be different, and it will // depend on particular engine implementation. inline unsigned nEta() const {return nEta_;} inline unsigned nPhi() const {return nPhi_;} protected: const unsigned nEta_; const unsigned nPhi_; private: // Disable the default constructor so that the derived classes // are forced to set "nEta_" and "nPhi_" correctly AbsFFTEngine(); }; } #endif // FFTJET_ABSFFTENGINE_HH_ Index: trunk/fftjet/AbsConvolverBase.hh =================================================================== --- trunk/fftjet/AbsConvolverBase.hh (revision 43) +++ trunk/fftjet/AbsConvolverBase.hh (revision 44) @@ -1,53 +1,53 @@ +#ifndef FFTJET_ABSCONVOLVERBASE_HH_ +#define FFTJET_ABSCONVOLVERBASE_HH_ + //========================================================================= // AbsConvolverBase.hh // // Abstract base class for performing kernel convolutions // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_ABSCONVOLVERBASE_HH_ -#define FFTJET_ABSCONVOLVERBASE_HH_ - #include namespace fftjet { template class AbsConvolverBase { public: virtual ~AbsConvolverBase() {} // Some inspectors of the functionality virtual bool fixingEfficiency() const = 0; virtual unsigned minFixedBin() const = 0; virtual unsigned maxFixedBin() const = 0; virtual unsigned nEtaFFT() const = 0; virtual unsigned nPhiFFT() const = 0; // In the next two methods, the "nEta" and "nPhi" arguments // which define array dimensionalities must correspond to the // grid dimensions of the FFT engine used to construct this object. // "setEventData" must be called at least once before any call // to "convolveWithKernel". virtual void setEventData(const Real* data, unsigned nEta, unsigned nPhi) = 0; virtual void convolveWithKernel(double scale, Real* result, unsigned nEta, unsigned nPhi) = 0; // Functions related to scale management virtual void processScale(double scale) = 0; // Find out how many scales were processed so far virtual unsigned nProcessedScales() const = 0; // Return the processed scales in the reverse sorted order virtual void getProcessedScales(std::vector* scales) const = 0; // Check whether a certain scale is already processed virtual bool isScaleProcessed(double scale) const = 0; }; } #endif // FFTJET_ABSCONVOLVERBASE_HH_ Index: trunk/fftjet/CompositeKernel.hh =================================================================== --- trunk/fftjet/CompositeKernel.hh (revision 43) +++ trunk/fftjet/CompositeKernel.hh (revision 44) @@ -1,75 +1,75 @@ +#ifndef FFTJET_COMPOSITEKERNEL_HH_ +#define FFTJET_COMPOSITEKERNEL_HH_ + //========================================================================= // CompositeKernel.hh // // A linear combination of kernels. The first member of the pair // is the factor for a kernel component in the combination. The // pointer to the kernel component itself is the second member // of the pair. Note that, if the factors do not sum up to 1, // the resulting kernel is not likely to be a density. // // This class is mainly intended for modeling different smearing // for charged and neutral parts of a jet (e.g., widening by // a magnetic field). It is up to the user of this class to make // sure that the combined result makes sense. In particular, // use of negative component factors should normally be avoided. // // If the "takePointerOwnership" argument is true, the object // will call the "delete" operator on the component kernel // pointers in the destructor -- that is, the object will own // the pointers. You have to decide on the ownership policy // at the construction time and follow it throughout the use // of the object. // // Note that there is no way to use different bandwidth values // for the component kernels. // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_COMPOSITEKERNEL_HH_ -#define FFTJET_COMPOSITEKERNEL_HH_ - #include #include #include "fftjet/AbsKernel2d.hh" namespace fftjet { class CompositeKernel : public AbsKernel2d, public std::vector > { public: // An empty CompositeKernel can be created. // The components can be subsequently added // by using the "push_back" function. explicit CompositeKernel(bool takePointerOwnership=false); CompositeKernel(const std::vector >&, bool takePointerOwnership=false); virtual ~CompositeKernel(); void setScaleRatio(double r); bool isDensity() const; double operator()(double x, double y, double scale) const; void supportRectangle( double scale, KernelSupportRectangle *r) const; double rectangleAverage(double x, double y, double scale, double dx, double dy) const; // The following function assumes that all component // kernels are normalized void random(double r1, double r2, double scale, double* px, double* py) const; inline bool takesPointerOwnership() const {return destroyKernelSet_;} private: CompositeKernel(const CompositeKernel&); CompositeKernel& operator=(const CompositeKernel&); const bool destroyKernelSet_; mutable std::vector cdf_; }; } #endif // FFTJET_COMPOSITEKERNEL_HH_ Index: trunk/fftjet/LinearInterpolator2d.hh =================================================================== --- trunk/fftjet/LinearInterpolator2d.hh (revision 43) +++ trunk/fftjet/LinearInterpolator2d.hh (revision 44) @@ -1,134 +1,134 @@ +#ifndef FFTJET_LINEARINTERPOLATOR2D_HH_ +#define FFTJET_LINEARINTERPOLATOR2D_HH_ + //========================================================================= // LinearInterpolator2d.hh // // This class can be used to interpolate histograms or similar regularly // spaced data in two dimensions // // I. Volobouev // August 2008 //========================================================================= -#ifndef FFTJET_LINEARINTERPOLATOR2D_HH_ -#define FFTJET_LINEARINTERPOLATOR2D_HH_ - #include #include #include namespace fftjet { class LinearInterpolator2d { public: // Constructor from a regularly spaced data. // Parameter "fOutside" can be used to specify the function // value in case the argument is outside the grid rectangle. // Parameter "useFOutside" tells whether, indeed, "fOutside" // should be used to do this or whether the function should // return the grid value at the closest rectangle edge. // // The "data" pointer can be NULL. In this case all internal // interpolator data will be set to 0. It is assumed that "setData" // method will be later used to provide the actual data. // template LinearInterpolator2d(const Real* data, unsigned nxpoints, double xmin, double xmax, unsigned nypoints, double ymin, double ymax, double fOutside=0.0, bool useFOutside=true); // Constructor which builds a function returning the given constant explicit LinearInterpolator2d(double c); // Copy constructor LinearInterpolator2d(const LinearInterpolator2d&); ~LinearInterpolator2d(); LinearInterpolator2d& operator=(const LinearInterpolator2d&); bool operator==(const LinearInterpolator2d& r) const; inline bool operator!=(const LinearInterpolator2d& r) const {return !(*this == r);} inline unsigned nx() const {return nxpoints;} inline unsigned ny() const {return nypoints;} inline double xMin() const {return xmin;} inline double xMax() const {return xmin + nxpoints*xbinwidth;} inline double yMin() const {return ymin;} inline double yMax() const {return ymin + nypoints*ybinwidth;} inline double outsideValue() const {return fOutside;} inline bool usingOutsideValue() const {return useFOutside;} inline const double* getData() const {return data;} double operator()(const double& x, const double& y) const; // The following function tells whether the function // values are non-negative (and, therefore, the function // can be treated as a density) bool isNonNegative() const; // The following function returns the function integral // inside the rectangle double integral() const; // The following function can be used to generate // random numbers according to the function represented // by the interpolator, assuming that the interpolated // function can be used as a density void random(double rnd1, double rnd2, double *p1, double *p2) const; // The following function can be used to change the // interpolated data "on the fly". The number of points // in the input array should be at least as large as nx()*ny(). template void setData(const Real* data); // The following function can be used to normalize the // function integral inside the rectangle to the given value void normalize(double value); // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static LinearInterpolator2d* read(std::istream& in); // Helper function for generating random numbers // inside a single cell static void singleCellRandom(double xlo, double ylo, double bwx, double bwy, double z00, double z01, double z10, double z11, double rnd1, double rnd2, double* p1, double* p2); private: LinearInterpolator2d(); void buildCdf(); double cornerIntegral(unsigned idx) const; double* data; double* cdf; unsigned nxpoints; unsigned nypoints; double xmin; double xbinwidth; double ymin; double ybinwidth; double fOutside; bool useFOutside; // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 5;} }; } #include "fftjet/LinearInterpolator2d.icc" #endif // FFTJET_LINEARINTERPOLATOR2D_HH_ Index: trunk/fftjet/Grid2d.hh =================================================================== --- trunk/fftjet/Grid2d.hh (revision 43) +++ trunk/fftjet/Grid2d.hh (revision 44) @@ -1,269 +1,269 @@ +#ifndef FFTJET_GRID2D_HH_ +#define FFTJET_GRID2D_HH_ + //========================================================================= // Grid2d.hh // // 2d grid (histogram) for deposited energy (Pt, etc.) with cylindrical // topology. The coverage in phi is 2 Pi. // // Note that the data from this grid will be fed directly into // the FFT engine. This means that the number of bins in this grid // should be appropriate for FFT and that there should be some // additional space in eta to avoid energy leakage between high // and low eta regions during convolution. It is better to add this // additional space on both sides, as the peaks will not be found // in the first and the last eta bin. // // This class is used inside some pretty tight loops. Don't create // any virtual functions in it. // // I. Volobouev // January 2008 //========================================================================= -#ifndef FFTJET_GRID2D_HH_ -#define FFTJET_GRID2D_HH_ - #include #include namespace fftjet { class StatAccumulator; template class Grid2d { public: // Constructors Grid2d(unsigned nEtaBins, Real etaMin, Real etaMax, unsigned nPhiBins, Real phiBin0Edge, const char* title = ""); Grid2d(const Grid2d&); // Destructor ~Grid2d(); // Assignment operator Grid2d& operator=(const Grid2d&); // Change title (could be useful after copying) void setTitle(const char* newtitle); // Comparison for equality bool operator==(const Grid2d& r) const; bool operator!=(const Grid2d& r) const; // Basic accessors unsigned nEta() const; unsigned nPhi() const; Real etaMin() const; Real etaMax() const; Real phiBin0Edge() const; const char* title() const; // Bin width quantities Real etaBinWidth() const; Real phiBinWidth() const; // Translators from coordinates to bin numbers. // Eta bin number can be negative or very large // which means that eta is out of histogram range. int getEtaBin(Real eta) const; unsigned getPhiBin(Real phi) const; // Translators from bin numbers to coordinates. Real etaBinCenter(int etaBin) const; Real phiBinCenter(unsigned phiBin) const; // Translators from real-valued bin numbers to coordinates. // These function perform the linear mapping from the coordinates // in the units of bin width to the normal coordinates. Real etaFromBinNum(Real etaBin) const; Real phiFromBinNum(Real phiBin) const; // The following function checks grid compatibility // for various calculations bool isCompatible(const Grid2d& g) const; // The following function checks whether all data // enstries are non-negative bool isNonNegative() const; // Finding minimum/maximum of the data values void findMinimum(unsigned* etaBin, unsigned* phiBin, Real* min) const; void findMaximum(unsigned* etaBin, unsigned* phiBin, Real* max) const; // Pointer to the data. To avoid any future incompatibilities, // do not manipulate the data "by hand". FFT engine needs // the direct data access for speed. It knows how the data // is laid out internally (and you are not supposed to). const Real* data() const; // Access to data values Real binValue(unsigned etaBin, unsigned phiBin) const; Real coordValue(Real eta, Real phi) const; // Values using unchecked array bounds Real uncheckedAt(unsigned etaBin, unsigned phiBin) const; // Sum of all data bins Real sum() const; // Integral of the data (sum of the data bins times bin area) Real integral() const; // Number of grid cells above some threshold value unsigned countAboveThreshold(Real threshold) const; // Integral of the data above some threshold value Real integralAboveThreshold(Real threshold) const; // Fuzzy integral of the data using the given threshold value. // Each datum is partitioned between the "signal" integral // (returned by this function) and "background" integral // according to the weights determined by datum/(threshold+datum) // and threshold/(threshold+datum), respectively. Note that this // partitioning makes sense only when both datum and threshold // are non-negative. The uncertainty of the integral is calculated // assuming that all cells are independent and that each cell // belongs completely to either signal or background. void fuzzyIntegral(Real thresh, Real *value, Real *uncertainty) const; // Accumulate statistics over the grid values (just the values, // not the coordinates). The accumulator will not be reset, so // the statistics can be accumulated from several grids (or from // multiple events). void accumulateDataStats(StatAccumulator* accumulator) const; // Fill function by coordinate. Overflows are ignored. void fill(Real eta, Real phi, Real energy); // The standard "fill" function fills the nearest four bins // in such a way that the center of mass of these fills // coincides with the coordinate given. Instead, the functions // below fill or set grid bins as in a histogram. Use these // functions (they are much faster) if you do not care about // loss of resolution due to binning. void fillFast(Real eta, Real phi, Real energy); void setFast(Real eta, Real phi, Real value); // Fill/set functions by bin number. Bin number // out of range will result in an assertion error. void fillBin(unsigned etaBin, unsigned phiBin, Real energy); void setBin(unsigned etaBin, unsigned phiBin, Real value); // Unchecked fill/set functions by bin number. // Can be used to improve the code speed in case // you are sure about your array bounds. void uncheckedFillBin(unsigned etaBin, unsigned phiBin, Real energy); void uncheckedSetBin(unsigned etaBin, unsigned phiBin, Real value); // Set everything. Assume that the data layout is correct. void blockSet(const Real *from, unsigned nEtaFrom, unsigned nPhiFrom); // Fill everything (add to internal data). // This function assumes that the data layout is correct. void blockFill(const Real *from, unsigned nEtaFrom, unsigned nPhiFrom); // Fill everything from another grid. The binning of the other // grid must be compatible Grid2d& operator+=(const Grid2d& r); // Subtract another grid. The binning of the other grid // must be compatible. Grid2d& operator-=(const Grid2d& r); // The following function sets all data values to the given constant void reset(Real value=0); // The following function multiplies all data values // by the given constant void scaleData(Real scaleFactor); // The following function multiplies data values by an eta-dependent // scale factor void scaleData(const Real* scaleFactorArray, unsigned arrayLength); // The following function sets all data below or equal to // the given threshold to 0 void applyThreshold(Real thresh); // The following function sets to 0 all data below or equal to // the ratio of the given threshold to cosh(eta). This is useful // for suppressing noise whose energy (not Et) distribution is // independent from eta while the grid itself is filled using Et. void applyCoshThreshold(Real thresh); // The following function sets all data below or equal to // the given threshold to 0 and above the given threshold to 1 void applyLogicalThreshold(Real thresh); // It may be interesting to cluster pow(Et, alpha) instead of // just Et. Apply the "power" function after filling the grid // to modify the data accordingly. Note that the "alpha" argument // cannot be negative, and that all grid entries less or equal // to 0 are set to 0 (even if "alpha" is 0). // // Doing this should be considered as an attempt to modify // the metric of the angular space (energies in real calorimeter // cells are added linearly no matter what we do here). One should // keep in mind that the occupancy away from the jet axis decays // approximately as r^-2 and that the average energy scales roughly // as r^-1, so that the angular energy profile decays as r^-3. // This means that, for example, using the "power" function with // argument 2 will cause the angular Et^2 profile to decay as r^-4. // The same asymptotic behavior of the Et profile can be achieved // by working with the variable r' = pow(r, k), k = 1/2. This is // because the occupancy behaves as r'^-2 no matter what k is // (logarithmic infrared divergence in the number of soft particles), // while average energy now decays as r'^(-1/k). In a similar vein, // setting "alpha" to 1/2 results in k = 2. // // It is not clear how important these arguments are in the presence // of calorimeter thresholds. void power(double alpha); // The following function clears some eta range. // Can be useful in order to prepare the data for FFT. // Bin number which corresponds to minEta is cleared // while bin number which corresponds to maxEta // is not cleared (but the last bin is cleared // in case maxEta is beyound the grid range). void clearEtaRange(Real minEta, Real maxEta); // Same thing but indexing by bin number. Bin number // etaBinMin is cleared and etaBinMax is not. void clearEtaBins(unsigned etaBinMin, unsigned etaBinMax); // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static Grid2d* read(std::istream& in); private: Grid2d(); void fillTwo(unsigned etaBin, Real phi, Real energy); Real* data_; std::string title_; Real phiBinWidth_; Real etaMin_; Real etaMax_; Real etaRange_; Real phiBin0Edge_; unsigned nEta_; unsigned nPhi_; // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 7;} }; } #include "fftjet/Grid2d.icc" #endif // FFTJET_GRID2D_HH_ Index: trunk/fftjet/JetMagnitudeMapper2d.hh =================================================================== --- trunk/fftjet/JetMagnitudeMapper2d.hh (revision 43) +++ trunk/fftjet/JetMagnitudeMapper2d.hh (revision 44) @@ -1,83 +1,83 @@ +#ifndef FFTJET_JETMAGNITUDEMAPPER2D_HH_ +#define FFTJET_JETMAGNITUDEMAPPER2D_HH_ + //========================================================================= // JetMagnitudeMapper2d.hh // // This jet corrector can used to invert a jet response curve // which depends on one of the jet properties (for example, scale or eta) // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_JETMAGNITUDEMAPPER2D_HH_ -#define FFTJET_JETMAGNITUDEMAPPER2D_HH_ - #include "fftjet/AbsJetCorrector.hh" #include "fftjet/LinearInterpolator2d.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template class JetMagnitudeMapper2d : public AbsJetCorrector { public: // The constructor arguments are as follows: // // f is the jet response curve (ratio of // reconstructed/actual jet magnitude). Should // define "double operator()(double, double) const". // The first variable is assumed to be the // predictor on which the response depends // (such as scale or eta) and the second // variable is the jet magnitude. // // predictorCalc Pointer to the functor which calculates // the predictor variable using the jet. // // ownPredictor Set this to "true" if the predictor should be // deleted in the destructor. // // minPredictor, Together, these three variables describe // maxPredictor, how to create the grid in the response // nPredPoints predictor dimension. // // maxMagnitude is the maximum value of _actual_ // magnitude for which the jet response // curve is build. Beyond this magnitude // the jet response is assumed to be constant // and equal to its value at maxMagnitude // (for each predictor value separately). // // nMagPoints number of reconstructed magnitude points // to use for modeling the inverse curve between 0 // and the reconstructed value which corresponds // to maxMagnitude // template JetMagnitudeMapper2d(const Functor2d& f, const Functor1* predictorCalc, bool ownPredictor, double minPredictor, double maxPredictor, unsigned nPredPoints, double maxMagnitude, unsigned nMagPoints); virtual ~JetMagnitudeMapper2d(); inline double operator()(const double& magnitude, const Jet& j) const {return magnitude*(*interp)((*fcn)(j), magnitude);} private: JetMagnitudeMapper2d(); JetMagnitudeMapper2d(const JetMagnitudeMapper2d&); JetMagnitudeMapper2d& operator=(const JetMagnitudeMapper2d&); LinearInterpolator2d* interp; const Functor1* fcn; const bool ownFcn; }; } #include "fftjet/JetMagnitudeMapper2d.icc" #endif // FFTJET_JETMAGNITUDEMAPPER2D_HH_ Index: trunk/fftjet/OpenDXPeakTree.hh =================================================================== --- trunk/fftjet/OpenDXPeakTree.hh (revision 43) +++ trunk/fftjet/OpenDXPeakTree.hh (revision 44) @@ -1,70 +1,70 @@ +#ifndef FFTJET_OPENDXPEAKTREE_HH_ +#define FFTJET_OPENDXPEAKTREE_HH_ + //========================================================================= // OpenDXPeakTree.hh // // This class provides a way to write out AbsClusteringTree // and SparseClusteringTree information in a form suitable // for visualization by OpenDX. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_OPENDXPEAKTREE_HH_ -#define FFTJET_OPENDXPEAKTREE_HH_ - #include #include "fftjet/AbsOpenDXTreeFormatter.hh" #include "fftjet/Peak.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template < typename LevelInfo, template class Tree > class OpenDXPeakTree : public AbsOpenDXTreeFormatter { public: inline OpenDXPeakTree( const Functor1* glyphSizeFunctor, const Functor1* glyphColorFunctor, const Tree& tree, const unsigned long runNum, const unsigned long evNum, const double etaRange = 0.0, const double scaleRange = 0.0) : AbsOpenDXTreeFormatter( tree, runNum, evNum, etaRange, scaleRange), glyphSizeFunctor_(glyphSizeFunctor), glyphColorFunctor_(glyphColorFunctor) {} inline OpenDXPeakTree( const Functor1* glyphSizeFunctor, const Functor1* glyphColorFunctor, const double etaRange = 0.0, const double scaleRange = 0.0) : AbsOpenDXTreeFormatter(etaRange,scaleRange), glyphSizeFunctor_(glyphSizeFunctor), glyphColorFunctor_(glyphColorFunctor) {} inline virtual ~OpenDXPeakTree() {} private: OpenDXPeakTree(); const Functor1* const glyphSizeFunctor_; const Functor1* const glyphColorFunctor_; inline virtual DXGlyphInfo buildGlyph(const Peak& peak, const int parent) const { return DXGlyphInfo(peak.eta(), peak.phi(), log(peak.scale()), (*glyphSizeFunctor_)(peak), (*glyphColorFunctor_)(peak), parent); } }; } #endif // FFTJET_OPENDXPEAKTREE_HH_ Index: trunk/fftjet/EquidistantSequence.hh =================================================================== --- trunk/fftjet/EquidistantSequence.hh (revision 43) +++ trunk/fftjet/EquidistantSequence.hh (revision 44) @@ -1,39 +1,39 @@ +#ifndef FFTJET_EQUIDISTANTSEQUENCE_HH_ +#define FFTJET_EQUIDISTANTSEQUENCE_HH_ + //======================================================================== // EquidistantSequence.hh // // A sequence of points equidistant in linear or log space // // I. Volobouev // March 2009 //======================================================================== -#ifndef FFTJET_EQUIDISTANTSEQUENCE_HH_ -#define FFTJET_EQUIDISTANTSEQUENCE_HH_ - #include namespace fftjet { class EquidistantInLinearSpace : public std::vector { public: EquidistantInLinearSpace(double minScale, double maxScale, unsigned nScales); virtual ~EquidistantInLinearSpace() {} private: EquidistantInLinearSpace(); }; class EquidistantInLogSpace : public std::vector { public: EquidistantInLogSpace(double minScale, double maxScale, unsigned nScales); virtual ~EquidistantInLogSpace() {} private: EquidistantInLogSpace(); }; } #endif // FFTJET_EQUIDISTANTSEQUENCE_HH_ Index: trunk/fftjet/JetProperty.hh =================================================================== --- trunk/fftjet/JetProperty.hh (revision 43) +++ trunk/fftjet/JetProperty.hh (revision 44) @@ -1,110 +1,110 @@ +#ifndef FFTJET_JETPROPERTY_HH_ +#define FFTJET_JETPROPERTY_HH_ + //========================================================================= // JetProperty.hh // // Selectors of properties for Peak/RecombinedJet classes. For use // with OpenDXTreeFormatter and in other similar contexts. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_JETPROPERTY_HH_ -#define FFTJET_JETPROPERTY_HH_ - #include #include #include "fftjet/SimpleFunctors.hh" namespace fftjet { // Generic wrapper for all properties which can be // expressed as jet member functions with no arguments // returning double template class JetProperty : public Functor1 { public: typedef double (Jet::* JetMemberFunction)() const; inline explicit JetProperty(JetMemberFunction fn) : fn_(fn) {} inline double operator()(const Jet& j) const {return (j.*fn_)();} private: JetProperty(); JetMemberFunction fn_; }; // Generic wrapper for taking logarithms of other properties template class LogProperty : public Functor1 { public: inline LogProperty(const Functor1* f, bool own=false) : func_(f), ownsPointer_(own) {} inline ~LogProperty() {if (ownsPointer_) delete func_;} inline double operator()(const Jet& j) const { const double value = (*func_)(j); assert(value > 0.0); return log(value); } private: LogProperty(); const Functor1* func_; const bool ownsPointer_; }; // We need couple extra constants to calculate the Laplacian template struct MinusScaledLaplacian : public Functor1 { inline MinusScaledLaplacian(const double bwx, const double bwy) : bwx_(bwx), bwy_(bwy) {} inline double operator()(const Jet& c) const { // The following normally makes the Laplacian non-negative return -c.scaledLaplacian(bwx_, bwy_); } private: MinusScaledLaplacian(); const double bwx_; const double bwy_; }; // Some potentially useful combinations of jet variables template struct ScaledHessianDet : public Functor1 { inline double operator()(const Jet& c) const { const double s = c.scale(); const double ssquared = s*s; return ssquared*ssquared*c.hessianDeterminant(); } }; template struct ScaledMagnitude : public Functor1 { inline double operator()(const Jet& c) const { return c.scale()*c.magnitude(); } }; template struct ScaledMagnitude2 : public Functor1 { inline double operator()(const Jet& c) const { const double s = c.scale(); return s*s*c.magnitude(); } }; } #endif // FFTJET_JETPROPERTY_HH_ Index: trunk/fftjet/FFTWEngine.hh =================================================================== --- trunk/fftjet/FFTWEngine.hh (revision 43) +++ trunk/fftjet/FFTWEngine.hh (revision 44) @@ -1,78 +1,78 @@ +#ifndef FFTJET_FFTWENGINE_HH_ +#define FFTJET_FFTWENGINE_HH_ + //========================================================================= // FFTWEngine.hh // // FFT engine for jet reconstruction which uses the FFTW library // from http://www.fftw.org. User code can not instantiate this // class -- instead, use FFTWDoubleEngine and FFTWFloatEngine. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_FFTWENGINE_HH_ -#define FFTJET_FFTWENGINE_HH_ - #include "fftjet/AbsFFTEngine.hh" namespace fftjet { template class FFTWEngine : public AbsFFTEngine { public: virtual ~FFTWEngine() {} void transformForward(const Real* from, Complex* to) const; void transformBack(const Complex* from, Real* to) const; void multiplyTransforms(const Complex* a, const Complex* b, Complex* result) const; void addTransforms(const Complex* a, const Complex* b, Complex* result) const; void divideTransforms(const Complex* a, const Complex* b, Complex* result) const; void deconvolutionRatio(const Complex* a, const Complex* b, Complex* result) const; void subtractTransforms(const Complex* a, const Complex* b, Complex* result) const; void scaleTransform(const Complex* a, Real scale, Complex* result) const; void amplitudeSquared(const Complex* a, Complex* result) const; double totalPower(const Complex* a) const; void setToReal(Complex* a, Real value) const; void setTransformPoint(Complex* points, unsigned iEta, unsigned iPhi, const std::complex& value) const; std::complex getTransformPoint(const Complex* points, unsigned iEta, unsigned iPhi) const; Complex* allocateComplex() const; void destroyComplex(Complex*) const; protected: FFTWEngine(unsigned nEta, unsigned nPhi); Real *in; Complex *out; Plan pf; Plan pb; void (*engineExecute)(const Plan); void* (*engineMalloc)(size_t); void (*engineFree)(void *); private: FFTWEngine(const FFTWEngine&); FFTWEngine& operator=(const FFTWEngine&); typedef AbsFFTEngine B; const unsigned nReal_; const unsigned nComplex_; }; } #include "fftjet/FFTWEngine.icc" #endif // FFTJET_FFTWENGINE_HH_ Index: trunk/fftjet/InterpolatedMembershipFcn.hh =================================================================== --- trunk/fftjet/InterpolatedMembershipFcn.hh (revision 43) +++ trunk/fftjet/InterpolatedMembershipFcn.hh (revision 44) @@ -1,146 +1,146 @@ +#ifndef FFTJET_INTERPOLATEDMEMBERSHIPFCN_HH_ +#define FFTJET_INTERPOLATEDMEMBERSHIPFCN_HH_ + //========================================================================= // InterpolatedMembershipFcn.hh // // Jet membership function whose data is tabulated on a 4-d grid // (partially zero suppressed). In between, the data is interpolated // linearly. For scales below the range of scales covered, the slice // with the smallest scale is used. For scales above the range the slice // with the largest scale is used. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_INTERPOLATEDMEMBERSHIPFCN_HH_ -#define FFTJET_INTERPOLATEDMEMBERSHIPFCN_HH_ - #include #include #include "fftjet/AbsMembershipFunction.hh" #include "fftjet/EtaPhiEtInterpolator.hh" namespace fftjet { template class InterpolatedMembershipFcn : public AbsMembershipFunction { public: // The "scales" vector in the constructor must represent // the scales arranged in the increasing order (not // necessarily equidistant). // // The "useLogSpaceForScale" argument specifies whether // the code should use scale or log(scale) as the variable // in which the interpolation is linear. // // The "energyFractionToIgnore" argument specifies the // energy fraction of the whole distribution which can be // ignored in order to improve the compression. Numbers // about 1.0e-3 or 1.0e-4 seem to work well. // // "nEtaBins", "etaMin", and "etaMax" arguments specify // the binning in the eta direction. It is assumed that // the first bin is centered at // etaMin + (etaMax - etaMin)/(2 nEtaBins), just like // in a histogram. Same thing for the energy variable. // // The phi variable is assumed to have 2 Pi range, // with the first bin edge at "phiBin0Edge". // InterpolatedMembershipFcn( const std::vector& scales, bool useLogSpaceForScale, Real energyFractionToIgnore, unsigned nEtaBins, Real etaMin, Real etaMax, unsigned nPhiBins, Real phiBin0Edge, unsigned nEnergyBins, Real energyMin, Real energyMax); virtual ~InterpolatedMembershipFcn(); // Comparison (useful for testing) bool operator==(const InterpolatedMembershipFcn& r) const; inline bool operator!=(const InterpolatedMembershipFcn& r) const {return !(*this == r);} // Basic accessors unsigned nEta() const; unsigned nPhi() const; unsigned nEnergy() const; Real ignoredFraction() const; Real etaMin() const; Real etaMax() const; Real phiBin0Edge() const; Real energyMin() const; Real energyMax() const; bool usesLogSpace() const; const std::vector& scales() const; // Set the data for one of the scales. The "data" array should // have nEtaBins*nPhiBins*nEnergyBins elements. It is assumed // that the array index associated with the energy variable // changes most often, while the array index associated with the // eta variable changes least often. template void setScaleData(unsigned scaleBinNumber, const InputReal* data, double dataScaleFactor=1.0); // The following function must be overriden from the base inline void setScaleRatio(double) {} // The following function returns "true" if the data were provided // for all the scales inline bool isReady() const {return ready_;} // Info for one of the scales Real totalEnergy(unsigned scaleBin) const; Real multiplicity(unsigned scaleBin) const; // AbsMembershipFunction methods to implement double operator()(double eta, double phi, double et, double scale) const; double absorbableEnergy(double eta, double phi, double scale) const; // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static InterpolatedMembershipFcn* read(std::istream& in); private: InterpolatedMembershipFcn(); InterpolatedMembershipFcn(const InterpolatedMembershipFcn&); InterpolatedMembershipFcn& operator=(const InterpolatedMembershipFcn&); const std::vector scales_; const Real ignoredFraction_; const Real etaMin_; const Real etaMax_; const Real phiBin0Edge_; const Real energyMin_; const Real energyMax_; const unsigned nEta_; const unsigned nPhi_; const unsigned nEnergy_; const bool useLogSpace_; bool ready_; std::vector*> interpols_; // The "setScale" function for use with I/O operators. // This class will assume ownership of the "interp" object. void setScale(unsigned scaleBin, EtaPhiEtInterpolator* interp); // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 2;} }; } #include "fftjet/InterpolatedMembershipFcn.icc" #endif // FFTJET_INTERPOLATEDMEMBERSHIPFCN_HH_ Index: trunk/fftjet/EtCentroidRecombinationAlg.hh =================================================================== --- trunk/fftjet/EtCentroidRecombinationAlg.hh (revision 43) +++ trunk/fftjet/EtCentroidRecombinationAlg.hh (revision 44) @@ -1,66 +1,66 @@ +#ifndef FFTJET_ETCENTROIDRECOMBINATIONALG_HH_ +#define FFTJET_ETCENTROIDRECOMBINATIONALG_HH_ + //========================================================================= // EtCentroidRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // This particular class assumes that the data grid is populated with // transverse energy values. // // 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 Et, Real eta, Real phi) const; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the energy centroid. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ETCENTROIDRECOMBINATIONALG_HH_ -#define FFTJET_ETCENTROIDRECOMBINATIONALG_HH_ - #include "fftjet/KernelRecombinationAlg.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class EtCentroidRecombinationAlg : public KernelRecombinationAlg { public: inline EtCentroidRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) : KernelRecombinationAlg( kernel, bgWeight, unlikelyBgWeight, dataCutoff, winnerTakesAll, buildCorrelationMatrix, buildClusterMask, etaBinMin, etaBinMax) {} inline virtual ~EtCentroidRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} }; } #endif // FFTJET_ETCENTROIDRECOMBINATIONALG_HH_ Index: trunk/fftjet/Kernel2dFactory.hh =================================================================== --- trunk/fftjet/Kernel2dFactory.hh (revision 43) +++ trunk/fftjet/Kernel2dFactory.hh (revision 44) @@ -1,85 +1,85 @@ +#ifndef FFTJET_KERNEL2DFACTORY_HH_ +#define FFTJET_KERNEL2DFACTORY_HH_ + //========================================================================= // Kernel2dFactory.hh // // Factories for 2d kernel functions // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_KERNEL2DFACTORY_HH_ -#define FFTJET_KERNEL2DFACTORY_HH_ - #include #include #include #include #include "fftjet/AbsKernel2d.hh" namespace fftjet { // Abstract factory for Kernel2d types class AbsKernel2dFactory { public: virtual ~AbsKernel2dFactory() {} virtual AbsKernel2d* create(double sx, double sy, int scalePower, const std::vector&) const = 0; virtual int nParameters() const = 0; }; // Concrete factory for Kernel2d types template class Kernel2dFactory : public AbsKernel2dFactory { public: virtual ~Kernel2dFactory() {} inline int nParameters() const {return T::nParameters();} T* create(double sx, double sy, int scalePower, const std::vector& params) const { if (nParameters() >= 0) assert(params.size() == static_cast(nParameters())); return new T(sx, sy, scalePower, params); } }; // Default kernel factory. It is intended for use inside some // interpretive language or testing framework and can iterate // over kernel classes. Usage is like this: // // DefaultKernel2dFactory factory; // AbsKernel2d *kernel = factory[userName]->create(sx, sy, p, // userParameters); // // where "userName" is one of the known kernel names, sx and sy are // doubles, p is an integer which determines the dependence of the // kernel bandwidth values on the scale, and "userParameters" is // a vector of doubles. At some point after this one has to call // "delete" operator on the kernel pointer. Known kernel names // correspond to the names of classes derived from "AbsKernel2d". // For example, to create the Gaussian kernel, one can use // // AbsKernel2d *kernel = factory["Gauss2d"]->create(1, 1, 1, // std::vector()); // // You can add your own kernel to this factory like this: // // factory["YourClass"] = new Kernel2dFactory; // class DefaultKernel2dFactory : public std::map { public: DefaultKernel2dFactory(); virtual ~DefaultKernel2dFactory(); private: DefaultKernel2dFactory(const DefaultKernel2dFactory&); DefaultKernel2dFactory& operator=(const DefaultKernel2dFactory&); }; } #endif // FFTJET_KERNEL2DFACTORY_HH_ Index: trunk/fftjet/FrequencyKernelConvolver.hh =================================================================== --- trunk/fftjet/FrequencyKernelConvolver.hh (revision 43) +++ trunk/fftjet/FrequencyKernelConvolver.hh (revision 44) @@ -1,46 +1,46 @@ +#ifndef FFTJET_FREQUENCYKERNELCONVOLVER_HH_ +#define FFTJET_FREQUENCYKERNELCONVOLVER_HH_ + //========================================================================= // FrequencyKernelConvolver.hh // // Class for performing kernel convolutions by FFT using a single kernel // whose representation is given in the frequency domain // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_FREQUENCYKERNELCONVOLVER_HH_ -#define FFTJET_FREQUENCYKERNELCONVOLVER_HH_ - #include "fftjet/AbsKernelConvolver.hh" #include "fftjet/AbsFrequencyKernel.hh" namespace fftjet { template class FrequencyKernelConvolver : public AbsKernelConvolver { public: // This FrequencyKernelConvolver will not own the AbsFFTEngine // and AbsFrequencyKernel objects. It is the responsibility // of the user of this class to ensure that these objects // are not destroyed before the KernelConvolver itself. // // If the "minFixBin" and "maxFixBin" arguments are // changed from their default values, it means that // the convolver will try to fix the bump reconstruction // efficiency (will apply weights to the data to simulate // the effects of boundary kernel). FrequencyKernelConvolver(const AbsFFTEngine* fftEngine, const AbsFrequencyKernel* kernel, unsigned minFixBin=0, unsigned maxFixBin=0); inline virtual ~FrequencyKernelConvolver() {} protected: virtual KernelData buildKernelImage(double scale); const AbsFrequencyKernel* const kernel; }; } #include "fftjet/FrequencyKernelConvolver.icc" #endif // FFTJET_FREQUENCYKERNELCONVOLVER_HH_ Index: trunk/fftjet/FasterEtCentroidRecombinationAlg.hh =================================================================== --- trunk/fftjet/FasterEtCentroidRecombinationAlg.hh (revision 43) +++ trunk/fftjet/FasterEtCentroidRecombinationAlg.hh (revision 44) @@ -1,66 +1,66 @@ +#ifndef FFTJET_FASTERETCENTROIDRECOMBINATIONALG_HH_ +#define FFTJET_FASTERETCENTROIDRECOMBINATIONALG_HH_ + //========================================================================= // FasterEtCentroidRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // This particular class assumes that the data grid is populated with // transverse energy values. // // 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 Et, Real eta, Real phi) const; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the energy centroid. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_FASTERETCENTROIDRECOMBINATIONALG_HH_ -#define FFTJET_FASTERETCENTROIDRECOMBINATIONALG_HH_ - #include "fftjet/FasterKernelRecombinationAlg.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class FasterEtCentroidRecombinationAlg : public FasterKernelRecombinationAlg { public: inline FasterEtCentroidRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) : FasterKernelRecombinationAlg( kernel, bgWeight, unlikelyBgWeight, dataCutoff, winnerTakesAll, buildCorrelationMatrix, buildClusterMask, etaBinMin, etaBinMax) {} inline virtual ~FasterEtCentroidRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} }; } #endif // FFTJET_FASTERETCENTROIDRECOMBINATIONALG_HH_ Index: trunk/fftjet/AbsKernelConvolver.hh =================================================================== --- trunk/fftjet/AbsKernelConvolver.hh (revision 43) +++ trunk/fftjet/AbsKernelConvolver.hh (revision 44) @@ -1,142 +1,142 @@ +#ifndef FFTJET_ABSKERNELCONVOLVER_HH_ +#define FFTJET_ABSKERNELCONVOLVER_HH_ + //========================================================================= // AbsKernelConvolver.hh // // Base class for performing kernel convolutions by FFT. This class manages // a collection of kernel scans, images and, if requested, the data // normalization curves. The actual building of kernel images is up to // derived classes. // // Normally, all scans should be performed using the same kernel or // set of kernels, and with the same grid dimensions. The scans should // differ only by the scale. // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_ABSKERNELCONVOLVER_HH_ -#define FFTJET_ABSKERNELCONVOLVER_HH_ - #include #include "fftjet/AbsConvolverBase.hh" #include "fftjet/AbsFFTEngine.hh" #include "fftjet/KernelData.hh" namespace fftjet { // Scan information for multiple scales. // This class will destroy all KernelData // data buffers in its own destructor. template class AbsKernelConvolver : public AbsConvolverBase { public: // This AbsKernelConvolver object will not own the AbsFFTEngine. // It is the responsibility of the user of this class to ensure // that this object is not destroyed before the AbsKernelConvolver // itself is destroyed. AbsKernelConvolver(const AbsFFTEngine* fftEngine, unsigned minFixBin, unsigned maxFixBin); virtual ~AbsKernelConvolver(); // Some trivial inspectors inline bool fixingEfficiency() const {return fixEfficiency;} inline unsigned minFixedBin() const {return minFixBin;} inline unsigned maxFixedBin() const {return maxFixBin;} inline unsigned nEtaFFT() const {return fftEngine->nEta();} inline unsigned nPhiFFT() const {return fftEngine->nPhi();} // In the next two methods, the "nEta" and "nPhi" arguments // which define array dimensionalities must correspond to the // grid dimensions of the FFT engine used to construct this object. // "setEventData" must be called at least once before any call // to "convolveWithKernel". void setEventData(const Real* data, unsigned nEta, unsigned nPhi); void convolveWithKernel(double scale, Real* result, unsigned nEta, unsigned nPhi); // There is no accessor method to get the event data back. // This is because there is no particular reason for this // class to keep the original event data in case the efficiency // needs no fixing. In this case only the Fourier image // of the data is stored internally. // The next method simply requests construction of kernel data // (scan, image, and, possibly, normalization) for the given scale. // After this, "isScaleProcessed" method will return "true" for // the given scale, and scan/image/normalization accessors will // return valid pointers. It is not necessary to call this function // explicitly if all you want is to convolve the data with // the kernel at the given scale. void processScale(double scale); // Find out how many scales were processed so far unsigned nProcessedScales() const; // Return the processed scales in the reverse sorted order void getProcessedScales(std::vector* scales) const; // Check whether a certain scale is already processed bool isScaleProcessed(double scale) const; // The following methods allow the user to examine the kernel // scans/images/normalization curves for the given scale. // These methods return NULL if the scale was never processed. // // Note that the pointer to KernelData can be // invalidated whenever a new scale is added to the system, // but pointers to the scan, image, and normalization curve // will remain valid. // const KernelData* getKernelData(double scale) const; const Real* getKernelScan(double scale) const; const Complex* getKernelImage(double scale) const; const Real* getNormalization(double scale) const; // The following function returns 0.0 in case the scale // was never processed double getScanArea(double scale) const; protected: // The following helper function can be used to build // the normalization curve from the scanned kernel Real* buildEffNorm(const Real* scan) const; // The following function must be implemented by the derived classes. // It must build the kernel scan, image and, if requested, the inverse // efficiency curve for data normalization. The image and scan buffers // must be allocated on the heap using fftEngine->allocateComplex() // and "new Real[]", respectively. Construction of the normalization // curve can be performed by the "buildEffNorm" function if the kernel // image is available in the normal space (not Fourier transformed). virtual KernelData buildKernelImage(double scale) = 0; // Input arguments const AbsFFTEngine* const fftEngine; const unsigned minFixBin; const unsigned maxFixBin; const bool fixEfficiency; // Complex buffer which can be used by the derived classes // inside the "buildKernelImage" function. It will be allocated // before "buildKernelImage" is called. Complex* complexBuffer; private: AbsKernelConvolver(); AbsKernelConvolver(const AbsKernelConvolver&); AbsKernelConvolver& operator=(const AbsKernelConvolver&); // Collection of kernel images and normalization curves std::map > kernelProperties; // Various memory buffers Complex* dataImage; const Real* dataBuffer; mutable Real* effBuf; }; } #include "fftjet/AbsKernelConvolver.icc" #endif // FFTJET_ABSKERNELCONVOLVER_HH_ Index: trunk/fftjet/cdfCellNumber.hh =================================================================== --- trunk/fftjet/cdfCellNumber.hh (revision 43) +++ trunk/fftjet/cdfCellNumber.hh (revision 44) @@ -1,24 +1,24 @@ +#ifndef FFTJET_CDFCELLNUMBER_HH_ +#define FFTJET_CDFCELLNUMBER_HH_ + //========================================================================= // cdfCellNumber.hh // // Utility functions for building inverse cumulative distributions // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_CDFCELLNUMBER_HH_ -#define FFTJET_CDFCELLNUMBER_HH_ - namespace fftjet { // A utility function for looking up the cumulative distribution // function cell number given the cdf value unsigned cdfCellNumber(double cdf, const double* cdfdata, unsigned ncells, unsigned stride=1); // Look for the solution of the equation a x^2 + b x == c // on the x interval from 0 to 1 double invertQuadraticCdf(double a, double b, double c); } #endif // FFTJET_CDFCELLNUMBER_HH_ Index: trunk/fftjet/EtaPhiEtInterpolator.hh =================================================================== --- trunk/fftjet/EtaPhiEtInterpolator.hh (revision 43) +++ trunk/fftjet/EtaPhiEtInterpolator.hh (revision 44) @@ -1,167 +1,167 @@ +#ifndef FFTJET_ETAPHIETINTERPOLATOR_HH_ +#define FFTJET_ETAPHIETINTERPOLATOR_HH_ + //========================================================================= // EtaPhiEtInterpolator.hh // // A class for an efficient storage of eta-phi-energy histograms. It is // assumed that for most eta-phi values only a small number of the lowest // energy bins are populated. This class is useful for representing // jet shapes. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ETAPHIETINTERPOLATOR_HH_ -#define FFTJET_ETAPHIETINTERPOLATOR_HH_ - #include #include namespace fftjet { template class EtaPhiEtInterpolator { public: // Constructor from a regularly spaced data. // Note that we will be using it on some noisy // data (produced by FFT). Therefore, we need // to ignore bins which contribute very little // to the overall energy. The "energyFractionToIgnore" // argument tells us how much energy we are allowed // to ignore. The size of the "data" array should be // at least nEtaBins*nEtaBins*nEnergyBins, and the // ordering of the dimensions should be [eta][phi][energy]. template EtaPhiEtInterpolator(const InputReal* data, Real energyFractionToIgnore, unsigned nEtaBins, Real etaMin, Real etaMax, unsigned nPhiBins, Real phiBin0Edge, unsigned nEnergyBins, Real energyMin, Real energyMax, const char* title = ""); // Copy constructor EtaPhiEtInterpolator(const EtaPhiEtInterpolator&); // Destructor ~EtaPhiEtInterpolator(); // Assignment operator EtaPhiEtInterpolator& operator=(const EtaPhiEtInterpolator&); // Comparison (useful for testing) bool operator==(const EtaPhiEtInterpolator& r) const; inline bool operator!=(const EtaPhiEtInterpolator& r) const {return !(*this == r);} // Change title (could be useful after copying) void setTitle(const char* newtitle); // Basic accessors unsigned nEta() const; unsigned nPhi() const; unsigned nEnergy() const; Real ignoredFraction() const; Real etaMin() const; Real etaMax() const; Real phiBin0Edge() const; Real energyMin() const; Real energyMax() const; const char* title() const; // The 3-d bin size Real binSize() const; // The length of the zero-suppressed data buffer unsigned zeroSuppressedBufLen() const; // Integral of (bin contents times by energy), // after removing the ignored fraction Real totalEnergy() const; // Integral of all data bins (after removing the ignored fraction) Real multiplicity() const; // Data lookup functions Real binValue(unsigned etaBin, unsigned phiBin, unsigned eBin) const; Real uncheckedAt(unsigned etaBin, unsigned phiBin, unsigned eBin) const; // Linearly interpolated value of the energy density Real operator()(Real eta, Real phi, Real et) const; // Maximum energy value for which the density // is not set to 0. The value will be linearly // interpolated from the four closest energy bins. Real maxNon0Energy(Real eta, Real phi) const; // Bin number lookups. Eta and energy bin numbers // can be negative or very large which means that // the input argument is out of histogram range. int getEtaBin(Real eta) const; unsigned getPhiBin(Real phi) const; int getEnergyBin(Real et) const; // Translators from bin numbers to coordinates Real etaBinCenter(int etaBin) const; Real phiBinCenter(unsigned phiBin) const; Real energyBinCenter(int eBin) const; // The following function multiplies all data values // by the given constant void scaleData(Real scaleFactor); // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static EtaPhiEtInterpolator* read(std::istream& in); private: EtaPhiEtInterpolator(); // Constructor from already sparsified data // (for use in read/write procedures). Assumes // ownership of "sparseData" and "binBounds". EtaPhiEtInterpolator(Real* sparseData, unsigned* binBounds, Real energyFractionToIgnore, Real totalEnergy, Real totalCount, unsigned nEtaBins, Real etaMin, Real etaMax, unsigned nPhiBins, Real phiBin0Edge, unsigned nEnergyBins, Real energyMin, Real energyMax, const char* title); Real* data_; unsigned* binBounds_; std::string title_; Real ignoredFraction_; Real totalEnergy_; Real totalCount_; Real phiBinWidth_; Real etaMin_; Real etaMax_; Real etaRange_; Real phiBin0Edge_; Real energyMin_; Real energyMax_; Real energyRange_; unsigned nEta_; unsigned nPhi_; unsigned nEnergy_; template static InputReal determineCutoff(const InputReal* data, unsigned angleBins, unsigned eBins, Real energyMin, Real energyMax, Real dataFractionToIgnore); // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 1;} }; } #include "fftjet/EtaPhiEtInterpolator.icc" #endif // FFTJET_ETAPHIETINTERPOLATOR_HH_ Index: trunk/fftjet/SpecialFunctions.hh =================================================================== --- trunk/fftjet/SpecialFunctions.hh (revision 43) +++ trunk/fftjet/SpecialFunctions.hh (revision 44) @@ -1,36 +1,36 @@ +#ifndef FFTJET_SPECIALFUNCTIONS_HH_ +#define FFTJET_SPECIALFUNCTIONS_HH_ + //========================================================================= // SpecialFunctions.hh // // A few mathematical special functions needed by this package. // Note that these are not optimized in any way, and may be slow. // They are not used in speed-critical parts of the package code. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_SPECIALFUNCTIONS_HH_ -#define FFTJET_SPECIALFUNCTIONS_HH_ - namespace fftjet { // Inverse cumulative distribition function for 1d Gaussian double inverseGaussCdf(double cdf); // Regularized incomplete beta function double incompleteBeta(double a, double b, double x); // Inverse regularized incomplete beta function double inverseIncompleteBeta(double a, double b, double x); // The gamma function for positive real arguments double Gamma(double x); // Incomplete gamma ratio and its complemented function double incompleteGamma(double a, double x); double incompleteGammaC(double a, double x); // Inverse incomplete gamma ratio double inverseIncompleteGamma(double a, double x); } #endif // FFTJET_SPECIALFUNCTIONS_HH_ Index: trunk/fftjet/StatAccumulator.hh =================================================================== --- trunk/fftjet/StatAccumulator.hh (revision 43) +++ trunk/fftjet/StatAccumulator.hh (revision 44) @@ -1,45 +1,45 @@ +#ifndef FFTJET_STATACCUMULATOR_HH_ +#define FFTJET_STATACCUMULATOR_HH_ + //========================================================================= // StatAccumulator.hh // // Simple single-pass (which means imprecise) accumulator of statistical // summary for some set of numbers. Calculates minimum value, maximum // value, mean, and standard deviation. "mean" and "stdev" functions // will cause an assertion error in case "accumulate" function was never // called after the object was created (or after it was reset). // // Use with care. // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_STATACCUMULATOR_HH_ -#define FFTJET_STATACCUMULATOR_HH_ - namespace fftjet { class StatAccumulator { public: StatAccumulator(); inline unsigned long count() const {return count_;} inline double min() const {return min_;} inline double max() const {return max_;} double mean() const; double stdev() const; void accumulate(double value); void accumulate(const StatAccumulator&); void reset(); private: long double sum_; long double sumsq_; double min_; double max_; unsigned long count_; }; } #endif // FFTJET_STATACCUMULATOR_HH_ Index: trunk/fftjet/AbsFrequencyKernel1d.hh =================================================================== --- trunk/fftjet/AbsFrequencyKernel1d.hh (revision 43) +++ trunk/fftjet/AbsFrequencyKernel1d.hh (revision 44) @@ -1,31 +1,31 @@ +#ifndef FFTJET_ABSFREQUENCYKERNEL1D_HH_ +#define FFTJET_ABSFREQUENCYKERNEL1D_HH_ + //========================================================================= // AbsFrequencyKernel1d.hh // // Interface class for 1d functions in the FFT transform space. // Can be used to represent filters. // // This interface is much simpler than the interface used for real // kernels. For the intended range of applications, we do not anticipate // the need to represent delta functions in the frequency space and to // generate random harmonics. // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_ABSFREQUENCYKERNEL1D_HH_ -#define FFTJET_ABSFREQUENCYKERNEL1D_HH_ - #include namespace fftjet { struct AbsFrequencyKernel1d { virtual ~AbsFrequencyKernel1d() {} // i is the harmonic number (can be both positive and negative) virtual std::complex operator()(int i, double scale) const = 0; }; } #endif // FFTJET_ABSFREQUENCYKERNEL1D_HH_ Index: trunk/fftjet/Peak.hh =================================================================== --- trunk/fftjet/Peak.hh (revision 43) +++ trunk/fftjet/Peak.hh (revision 44) @@ -1,257 +1,258 @@ +#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 farmost daughter +// 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 //========================================================================= -#ifndef FFTJET_PEAK_HH_ -#define FFTJET_PEAK_HH_ - -#include -#include #include -#include "fftjet/ScaleSpaceKernel.hh" - namespace fftjet { + class ScaleSpaceKernel; + class Peak { public: // This object will not own the "membershipFunction" pointer - Peak(double eta, double phi, double mag, + 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 int code() const {return code_;} inline int status() const {return status_;} - void hessian(double hessianArray[3]) const; - double hessianDeterminant() const; - inline double laplacian() const {return hessian_[0] + hessian_[2];} - double scaledLaplacian(double bwEta, double bwPhi) const; - inline ScaleSpaceKernel* membershipFunction() const - {return memFcn_;} + 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 setCode(const int c) {code_= c;} inline void setStatus(const int s) {status_= s;} - void setHessian(const double hessian[3]); - inline void setMembershipFunction(ScaleSpaceKernel* f) - {memFcn_ = f;} + 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]; int code_; int status_; ScaleSpaceKernel* memFcn_; }; } #endif // FFTJET_PEAK_HH_ Index: trunk/fftjet/EtCentroidVectorRecombinationAlg.hh =================================================================== --- trunk/fftjet/EtCentroidVectorRecombinationAlg.hh (revision 43) +++ trunk/fftjet/EtCentroidVectorRecombinationAlg.hh (revision 44) @@ -1,61 +1,61 @@ +#ifndef FFTJET_ETCENTROIDVECTORRECOMBINATIONALG_HH_ +#define FFTJET_ETCENTROIDVECTORRECOMBINATIONALG_HH_ + //========================================================================= // EtCentroidVectorRecombinationAlg.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; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the energy centroid. // // I. Volobouev // June 2009 //========================================================================= -#ifndef FFTJET_ETCENTROIDVECTORRECOMBINATIONALG_HH_ -#define FFTJET_ETCENTROIDVECTORRECOMBINATIONALG_HH_ - #include "fftjet/KernelVectorRecombinationAlg.hh" namespace fftjet { template class EtCentroidVectorRecombinationAlg : public KernelVectorRecombinationAlg { public: typedef double (VectorLike::* VectorLikeMemberFunction)() const; inline EtCentroidVectorRecombinationAlg( ScaleSpaceKernel* kernel, VectorLikeMemberFunction etFcn, VectorLikeMemberFunction etaFcn, VectorLikeMemberFunction phiFcn, const Functor2* bgWeight, double unlikelyBgWeight, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask) : KernelVectorRecombinationAlg( kernel, etFcn, etaFcn, phiFcn, bgWeight, unlikelyBgWeight, winnerTakesAll, buildCorrelationMatrix, buildClusterMask) {} inline virtual ~EtCentroidVectorRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} }; } #endif // FFTJET_ETCENTROIDVECTORRECOMBINATIONALG_HH_ Index: trunk/fftjet/MultiKernelConvolver.hh =================================================================== --- trunk/fftjet/MultiKernelConvolver.hh (revision 43) +++ trunk/fftjet/MultiKernelConvolver.hh (revision 44) @@ -1,76 +1,76 @@ +#ifndef FFTJET_MULTIKERNELCONVOLVER_HH_ +#define FFTJET_MULTIKERNELCONVOLVER_HH_ + //========================================================================= // MultiKernelConvolver.hh // // Class for performing kernel convolutions by FFT using a set of kernels // // I. Volobouev // April 2008 //========================================================================= -#ifndef FFTJET_MULTIKERNELCONVOLVER_HH_ -#define FFTJET_MULTIKERNELCONVOLVER_HH_ - #include #include "fftjet/AbsKernelConvolver.hh" #include "fftjet/AbsKernel2d.hh" #include "fftjet/AbsFrequencyKernel.hh" namespace fftjet { class KernelSet { public: explicit KernelSet(bool ownsPointers, double regularizationFraction=0.0); ~KernelSet(); inline bool ownsPointers() const {return ownsPointers_;} inline double regularizationFraction() const {return regFraction_;} bool isEmpty() const; // The regularization fraction (fraction of frequencies to zero out) // should be between 0 and 1. void setRegularizationFraction(double fraction); std::vector filter; std::vector numerator; std::vector denominator; std::vector denoiser; private: KernelSet(); KernelSet(const KernelSet&); KernelSet& operator=(const KernelSet&); double regFraction_; const bool ownsPointers_; }; template class MultiKernelConvolver : public AbsKernelConvolver { public: // This MultiKernelConvolver will not own the AbsFFTEngine // and KernelSet objects. It is the responsibility of // the user of this class to ensure that these objects // are not destroyed before the MultiKernelConvolver itself. MultiKernelConvolver(const AbsFFTEngine* fftEngine, const KernelSet* kernelSequence, unsigned minFixBin=0, unsigned maxFixBin=0); virtual ~MultiKernelConvolver(); protected: virtual KernelData buildKernelImage(double scale); const KernelSet* kernelSequence; private: void regularizeDenominator(Complex* imageBuf) const; Complex* denomBuf; }; } #include "fftjet/MultiKernelConvolver.icc" #endif // FFTJET_MULTIKERNELCONVOLVER_HH_ Index: trunk/fftjet/AbsDistanceCalculator.hh =================================================================== --- trunk/fftjet/AbsDistanceCalculator.hh (revision 43) +++ trunk/fftjet/AbsDistanceCalculator.hh (revision 44) @@ -1,39 +1,39 @@ +#ifndef FFTJET_ABSDISTANCECALCULATOR_HH_ +#define FFTJET_ABSDISTANCECALCULATOR_HH_ + //========================================================================= // AbsDistanceCalculator.hh // // Interface class for the distance calculator in the proximity-matched // clustering tree // // I. Volobouev // April 2008 //========================================================================= -#ifndef FFTJET_ABSDISTANCECALCULATOR_HH_ -#define FFTJET_ABSDISTANCECALCULATOR_HH_ - namespace fftjet { template class AbsDistanceCalculator { public: virtual ~AbsDistanceCalculator() {} inline double operator()(const double scale1, const Cluster& jet1, const double scale2, const Cluster& jet2) const { if (scale1 < scale2) return distanceBetweenClusters(scale1, jet1, scale2, jet2); else return distanceBetweenClusters(scale2, jet2, scale1, jet1); } private: virtual double distanceBetweenClusters( double smallerScale, const Cluster& smallerScaleCluster, double biggerScale, const Cluster& biggerScaleCluster) const = 0; }; } #endif // FFTJET_ABSDISTANCECALCULATOR_HH_ Index: trunk/fftjet/InterpolatedKernel.hh =================================================================== --- trunk/fftjet/InterpolatedKernel.hh (revision 43) +++ trunk/fftjet/InterpolatedKernel.hh (revision 44) @@ -1,75 +1,75 @@ +#ifndef FFTJET_INTERPOLATEDKERNEL_HH_ +#define FFTJET_INTERPOLATEDKERNEL_HH_ + //========================================================================= // InterpolatedKernel.hh // // Kernel whose data is tabulated on a 2d grid (histogram bin centers). // In between, the data is interpolated linearly. // // I. Volobouev // August 2008 //========================================================================= -#ifndef FFTJET_INTERPOLATEDKERNEL_HH_ -#define FFTJET_INTERPOLATEDKERNEL_HH_ - #include #include #include "fftjet/AbsScalableKernel.hh" #include "fftjet/LinearInterpolator2d.hh" namespace fftjet { class InterpolatedKernel : public AbsScalableKernel { public: template inline InterpolatedKernel(double xScaleFactor, double yScaleFactor, int scalePow, const Real* data, unsigned nx, double xmin, double xmax, unsigned ny, double ymin, double ymax) : AbsScalableKernel(xScaleFactor, yScaleFactor, scalePow), in(data, nx, xmin, xmax, ny, ymin, ymax) { assert(data); in.normalize(1.0); isDensity_ = in.isNonNegative(); } virtual ~InterpolatedKernel() {} bool operator==(const InterpolatedKernel& r) const; inline bool operator!=(const InterpolatedKernel& r) const {return !(*this == r);} inline bool isDensity() const {return isDensity_;} // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static InterpolatedKernel* read(std::istream& in); private: // Constructor for the binary IO inline InterpolatedKernel(double xScaleFactor, double yScaleFactor, int scalePow, const LinearInterpolator2d& i) : AbsScalableKernel(xScaleFactor, yScaleFactor, scalePow), in(i), isDensity_(in.isNonNegative()) {} inline double calculate(double x, double y) const {return in(x, y);} void unscaledRectangle(KernelSupportRectangle *r) const; void unscaledRandom(double r1, double r2, double* px, double* py) const; LinearInterpolator2d in; bool isDensity_; // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 4;} }; } #endif // FFTJET_INTERPOLATEDKERNEL_HH_ Index: trunk/fftjet/DXGlyphInfo.hh =================================================================== --- trunk/fftjet/DXGlyphInfo.hh (revision 43) +++ trunk/fftjet/DXGlyphInfo.hh (revision 44) @@ -1,67 +1,67 @@ +#ifndef FFTJET_DXGLYPHINFO_HH_ +#define FFTJET_DXGLYPHINFO_HH_ + //========================================================================= // DXGlyphInfo.hh // // A little helper class which colects all the info necessary to build // OpenDX glyphs. Also, this header declares several stand-alone // glyph-related utility functions. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_DXGLYPHINFO_HH_ -#define FFTJET_DXGLYPHINFO_HH_ - #include #include namespace fftjet { class DXGlyphInfo { public: inline DXGlyphInfo() : eta_(0.0), phi_(0.0), scale_(0.0), size_(0.0), color_(0.0), parent_(-1) {} inline DXGlyphInfo(double eta, double phi, double scale, double size, double color, int parent) : eta_(eta), phi_(phi), scale_(scale), size_(size), color_(color), parent_(parent) {} // Accessors inline double eta() const {return eta_;} inline double phi() const {return phi_;} inline double scale() const {return scale_;} inline double size() const {return size_;} inline double color() const {return color_;} inline int parent() const {return parent_;} // Necessary modifiers inline void setParent(const int par) {parent_ = par;} inline void setScale(const double s) {scale_ = s;} private: double eta_; double phi_; double scale_; double size_; double color_; int parent_; }; // The following function modifies the collection of DXGlyphInfo // objects in order to produce short connections through the // phi = 0 = 2*pi section. void fixDXCylindricalGeometry(std::vector* glyphs, double etaRange); // The following function remaps the scales. It does nothing // if the "range" parameter is 0.0. void remapDXscales(std::vector* glyphs, double range); // The following function writes out a vector of DXGlyphInfo objects void writeDXGlyphVector(const std::vector& glyphs, unsigned long runNumber, unsigned long eventNumber, std::ostream& os); } #endif // FFTJET_DXGLYPHINFO_HH_ Index: trunk/fftjet/DiscreteGauss2d.hh =================================================================== --- trunk/fftjet/DiscreteGauss2d.hh (revision 43) +++ trunk/fftjet/DiscreteGauss2d.hh (revision 44) @@ -1,58 +1,58 @@ +#ifndef FFTJET_DISCRETEGAUSS2D_HH_ +#define FFTJET_DISCRETEGAUSS2D_HH_ + //========================================================================= // DiscreteGauss2d.hh // // The Fourier transform of the Gaussian kernel, corrected for the grid // discretization effects // // I. Volobouev // May 2009 //========================================================================= -#ifndef FFTJET_DISCRETEGAUSS2D_HH_ -#define FFTJET_DISCRETEGAUSS2D_HH_ - #include "fftjet/AbsFrequencyKernel.hh" namespace fftjet { class DiscreteGauss2d : public AbsFrequencyKernel { public: // Parameters "sEta" and "sPhi" are scale factors for // eta and phi directions, respectively. They have the same // meaning as the corresponding "sx", "sy" parameters for // the "Gauss2d" kernel. // // Parameters "nEta" and "nPhi" represent the number of cells // in the eta-phi discretization grid with which this kernel // will be used. // DiscreteGauss2d(double sEta, double sPhi, unsigned nEta, unsigned nPhi); inline ~DiscreteGauss2d() {} inline double sEta() const {return sEta_;} inline double sPhi() const {return sPhi_;} inline unsigned nEta() const {return nEta_;} inline unsigned nPhi() const {return nPhi_;} // The following function checks whether a grid with // given numbers of eta and phi bins is compatible // with this kernel inline bool isCompatible(const unsigned etaBins, const unsigned phiBins) const {return nEta_ == etaBins && nPhi_ == phiBins;} // The following inherited function must be overriden std::complex operator()(int ix, int iy, double scale) const; private: DiscreteGauss2d(); const double sEta_; const double sPhi_; const unsigned nEta_; const unsigned nPhi_; }; } #endif // FFTJET_DISCRETEGAUSS2D_HH_ Index: trunk/fftjet/AbsScalablePhiKernel.hh =================================================================== --- trunk/fftjet/AbsScalablePhiKernel.hh (revision 43) +++ trunk/fftjet/AbsScalablePhiKernel.hh (revision 44) @@ -1,113 +1,113 @@ +#ifndef FFTJET_ABSSCALABLEPHIKERNEL_HH_ +#define FFTJET_ABSSCALABLEPHIKERNEL_HH_ + //========================================================================= // AbsScalablePhiKernel.hh // // Interface class for kernel functions which are products of delta // function in eta and some reasonable distribution in phi. The shape // of the distribution in phi is assumed to be independent of the scale. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSSCALABLEPHIKERNEL_HH_ -#define FFTJET_ABSSCALABLEPHIKERNEL_HH_ - #include #include "fftjet/AbsPhiKernel.hh" namespace fftjet { class AbsScalablePhiKernel : public AbsPhiKernel { public: inline AbsScalablePhiKernel(const double scaleFactor, const int scalePower) : phis_(scaleFactor), scalePower_(scalePower) { assert(phis_); } virtual ~AbsScalablePhiKernel() {} private: AbsScalablePhiKernel(); inline double phiFcn(const double y, const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwy(phis_*pscale); return unscaledPhiFcn(y/bwy)/bwy; } inline void phiSupport(const double scale, double *phimin, double *phimax) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwy(phis_*pscale); unscaledPhiSupport(phimin, phimax); *phimin *= bwy; *phimax *= bwy; } inline double phiRandom(const double rnd, const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwy(phis_*pscale); return bwy*unscaledPhiRandom(rnd); } virtual double unscaledPhiFcn(double phi) const = 0; virtual void unscaledPhiSupport(double *phimin, double *phimax) const = 0; virtual double unscaledPhiRandom(double rnd) const = 0; const double phis_; const int scalePower_; }; } #endif // FFTJET_ABSSCALABLEPHIKERNEL_HH_ Index: trunk/fftjet/AbsSymmetricKernel.hh =================================================================== --- trunk/fftjet/AbsSymmetricKernel.hh (revision 43) +++ trunk/fftjet/AbsSymmetricKernel.hh (revision 44) @@ -1,122 +1,123 @@ +#ifndef FFTJET_ABSSYMMETRICKERNEL_HH_ +#define FFTJET_ABSSYMMETRICKERNEL_HH_ + //========================================================================= // AbsSymmetricKernel.hh // // Interface class for scalable 2d kernel functions with rotational // symmetry // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSSYMMETRICKERNEL_HH_ -#define FFTJET_ABSSYMMETRICKERNEL_HH_ - #include + #include "fftjet/AbsScalableKernel.hh" namespace fftjet { class AbsSymmetricKernel : public AbsScalableKernel { public: inline AbsSymmetricKernel(const double xScaleFactor, const double yScaleFactor, const int scalePower) : AbsScalableKernel(xScaleFactor, yScaleFactor, scalePower), halfW_(-1.0), reso_(-1.0) {} virtual ~AbsSymmetricKernel() {} // The following function returns the radius of the circle // enclosing one half of the unscaled kernel mass. // This quantity is useful for selecting compatible scales // for different kernels. This function is just a memoized // version of "calculateHalfWeightRadius()". inline double halfWeightRadius() const { if (halfW_ < 0.0) halfW_ = calculateHalfWeightRadius(); return halfW_; } // The following function returns the unscaled kernel mass // within the circle of radius r. The "nsteps" parameter can be // used to specify the number of numerical integration steps. // Override if a better than default implementation is available. virtual double axialWeight(double r, unsigned nsteps=1000) const; // The following function returns the minimum value of r0 // for which convolution of the unscaled kernel with the function // delta(r - r0)/(2 Pi r) results in a local minimum at the center // of coordinates (0, 0). Here, "delta" is the Dirac delta function. // This function is just a memoized version of "calculateResolution()". inline double circleResolution() const { if (reso_ < 0.0) reso_ = calculateResolution(); return reso_; } protected: // The following function is used to implement the default // behavior of the "calculateResolution()" function. The scan // of the radial distance is performed from the lowest value // "rmin" to some reasonably selected large value by multiplying // the distance with "scaleFactor" on every cycle. "scaleFactor" // argument must be slightly larger than 1.0. If "rmin" argument // is 0.0 or negative, a reasonable starting radius will be // selected automatically. "nkernels" kernels are placed // around the circle at this increasing radial distance. // The cycling stops (it is switched to root finding) when // a minimum is found at the circle center. double scanForResolution(unsigned nkernels, double rmin=0.0, double scaleFactor=1.01, double eps=1.e-6) const; private: mutable double halfW_; mutable double reso_; inline double calculate(const double x, const double y) const { return eval(x*x + y*y); } inline void unscaledRectangle(KernelSupportRectangle *r) const { const double d(supportDistance()); r->xmin = -d; r->xmax = d; r->ymin = -d; r->ymax = d; } inline void unscaledRandom(const double r1, const double r2, double* px, double* py) const { const double r(randomRadius(r1)); const double phi(2.0*M_PI*r2); *px = r*cos(phi); *py = r*sin(phi); } // The following functions must be implemented // by the derived classes virtual double eval(double rsquared) const = 0; virtual double supportDistance() const = 0; virtual double randomRadius(double rnd) const = 0; // Override the following functions if better than default // implementations are available inline virtual double calculateHalfWeightRadius() const { return randomRadius(0.5); } inline virtual double calculateResolution() const { return scanForResolution(6); } // Helper function for calculating the sum value of n kernels // placed on a circle with radius "a" double circularResponse(double x, double y, unsigned n, double a) const; }; } #endif // FFTJET_ABSSYMMETRICKERNEL_HH_ Index: trunk/fftjet/PeakSelectors.hh =================================================================== --- trunk/fftjet/PeakSelectors.hh (revision 43) +++ trunk/fftjet/PeakSelectors.hh (revision 44) @@ -1,79 +1,79 @@ +#ifndef FFTJET_PEAKSELECTORS_HH_ +#define FFTJET_PEAKSELECTORS_HH_ + //========================================================================= // PeakSelectors.hh // // Some simple functors for peak selection // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_PEAKSELECTORS_HH_ -#define FFTJET_PEAKSELECTORS_HH_ - #include #include "fftjet/SimpleFunctors.hh" namespace fftjet { // Simplistic, scale and data-independent peak selector class SimplePeakSelector : public Functor1 { public: explicit SimplePeakSelector(double magCut, double driftSpeedCut=DBL_MAX, double magSpeedCut=-DBL_MAX, double lifeTimeCut=-DBL_MAX, double NNDCut=-DBL_MAX, double etaCut=DBL_MAX, double splitTimeCut=-DBL_MAX, double mergeTimeCut=-DBL_MAX); virtual ~SimplePeakSelector() {} virtual bool operator()(const Peak& peak) const; private: SimplePeakSelector(); double magCut_; double driftSpeedCut_; double magSpeedCut_; double lifeTimeCut_; double NNDCut_; double etaCut_; double splitTimeCut_; double mergeTimeCut_; }; // Selector based on scale-dependent magnitude. Intended for use // before clustering tree is constructed. The magnitude cut formula // looks like this: // // a*pow(scale,p) + b // // Note that, normally, "a" should be positive, "b" should be // non-negative, and "p" should be a negative number smaller // or equal to 2 in magnitude. If these conditions are not // satisfied, the peaks produced by stand-alone energy depositions // will "pop up" with increasing scale rather than disappear. // class ScalePowerPeakSelector : public Functor1 { public: ScalePowerPeakSelector(double a, double p, double b, double etaCut=DBL_MAX); virtual ~ScalePowerPeakSelector() {} virtual bool operator()(const Peak& peak) const; private: ScalePowerPeakSelector(); double a_; double p_; double b_; double etaCut_; }; } #include "fftjet/PeakSelectors.icc" #endif // FFTJET_PEAKSELECTORS_HH_ Index: trunk/fftjet/AbsTreeFormatter.hh =================================================================== --- trunk/fftjet/AbsTreeFormatter.hh (revision 43) +++ trunk/fftjet/AbsTreeFormatter.hh (revision 44) @@ -1,75 +1,75 @@ +#ifndef FFTJET_ABSTREEFORMATTER_HH_ +#define FFTJET_ABSTREEFORMATTER_HH_ + //========================================================================= // AbsTreeFormatter.hh // // This class provides an interface for writing out AbsClusteringTree // and SparseClusteringTree information for subsequent visualization. // // The derived classes must implement the "write" function. After this, // formatters can be written out to an ostream, and it is not necessary // to implement separate operator<< for them. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSTREEFORMATTER_HH_ -#define FFTJET_ABSTREEFORMATTER_HH_ - #include #include namespace fftjet { template < typename Cluster, typename LevelInfo, template class AbsTree > class AbsTreeFormatter { public: inline AbsTreeFormatter() : tree_(0), run_(0), n_(0) {} // The "runNum" and "evNum" arguments are some kind of // an identifier for the tree, typically run and event numbers inline AbsTreeFormatter( const AbsTree& tree, const unsigned long runNum, const unsigned long evNum) : tree_(&tree), run_(runNum), n_(evNum) {} inline virtual ~AbsTreeFormatter() {} inline void setTree(const AbsTree& tree, unsigned long runNum, unsigned long evNum) {tree_ = &tree; run_ = runNum; n_ = evNum;} inline void writeTree(std::ostream& os) const {assert(tree_); this->write(*tree_, run_, n_, os);} private: const AbsTree* tree_; unsigned long run_; unsigned long n_; virtual void write(const AbsTree& tree, unsigned long runNum, unsigned long evNum, std::ostream& os) const = 0; }; } template < typename Cluster, typename LevelInfo, template class AbsTree > inline std::ostream& operator<<( std::ostream& strm, const fftjet::AbsTreeFormatter& f) { f.writeTree(strm); return strm; } #endif // FFTJET_ABSTREEFORMATTER_HH_ Index: trunk/fftjet/RecombinedJet.hh =================================================================== --- trunk/fftjet/RecombinedJet.hh (revision 43) +++ trunk/fftjet/RecombinedJet.hh (revision 44) @@ -1,224 +1,233 @@ +#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. // // 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 //========================================================================= -#ifndef FFTJET_RECOMBINEDJET_HH_ -#define FFTJET_RECOMBINEDJET_HH_ - #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 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), 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_;} // 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 int code() const {return peak_.code();} - inline int status() const {return peak_.status();} inline void hessian(double hessianArray[3]) const {peak_.hessian(hessianArray);} + 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 ScaleSpaceKernel* membershipFunction() const - {return peak_.membershipFunction();} + 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 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 pileup_; double uncertainty_; double convergenceDistance_; }; } #include "fftjet/RecombinedJet.icc" #endif // FFTJET_RECOMBINEDJET_HH_ Index: trunk/fftjet/GaussianNoiseMembershipFcn.hh =================================================================== --- trunk/fftjet/GaussianNoiseMembershipFcn.hh (revision 43) +++ trunk/fftjet/GaussianNoiseMembershipFcn.hh (revision 44) @@ -1,55 +1,55 @@ +#ifndef FFTJET_GAUSSIANNOISEMEMBERSHIPFCN_HH_ +#define FFTJET_GAUSSIANNOISEMEMBERSHIPFCN_HH_ + //========================================================================= // GaussianNoiseMembershipFcn.hh // // Noise membership functor for a simple model in which the noise // in each digitization cell is Gaussian, pedestal is 0, and energy // spectrum of the signal does not change substantially near 0. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_GAUSSIANNOISEMEMBERSHIPFCN_HH_ -#define FFTJET_GAUSSIANNOISEMEMBERSHIPFCN_HH_ - #include #include #include "fftjet/SimpleFunctors.hh" namespace fftjet { class GaussianNoiseMembershipFcn : public Functor2 { public: // The "prior" paramerer should be inversely proportional // to an apriori estimate of the signal occupancy. inline GaussianNoiseMembershipFcn(const double minWeight, const double prior) : minWeight_(minWeight), prior_(prior), sqrt2pi_(sqrt(2.0*M_PI)) { assert(minWeight_ >= 0.0); assert(prior_ >= 0.0); } inline double operator()(const double& et, const double& sigma) const { if (prior_) { assert(et >= 0.0 && sigma > 0.0); const double x = et/sigma; return minWeight_ + prior_/sqrt2pi_/sigma*exp(-x*x/2.0); } else return minWeight_; } private: GaussianNoiseMembershipFcn(); const double minWeight_; const double prior_; const double sqrt2pi_; }; } #endif // FFTJET_GAUSSIANNOISEMEMBERSHIPFCN_HH_ Index: trunk/fftjet/SparseClusteringTree.hh =================================================================== --- trunk/fftjet/SparseClusteringTree.hh (revision 43) +++ trunk/fftjet/SparseClusteringTree.hh (revision 44) @@ -1,348 +1,348 @@ +#ifndef FFTJET_SPARSECLUSTERINGTREE_HH_ +#define FFTJET_SPARSECLUSTERINGTREE_HH_ + //========================================================================= // SparseClusteringTree.hh // // Hierarchical clustering tree without complete level-by-level structure. // // Do not create this object every event -- this will lead to many // unnecessary memory allocations and deallocations. Instead, create it // at the start of your program and reuse. // // I. Volobouev // June 2010 //========================================================================= -#ifndef FFTJET_SPARSECLUSTERINGTREE_HH_ -#define FFTJET_SPARSECLUSTERINGTREE_HH_ - #include #include #include #include "fftjet/SmallVector.hh" namespace fftjet { // Same template parameters as in AbsClusteringTree. The second // parameter is currently unused, but included here anyway for // uniformity with AbsClusteringTree (this is useful for generic // code in the I/O classes). template class SparseClusteringTree { public: typedef Cluster cluster_type; typedef LevelInfo info_type; typedef unsigned NodeId; // Node mask bits enum { MIN_LEVEL = 1U, // level with the largest scale MAX_LEVEL = 2U, // level with the smallest scale USER_LEVEL = 4U, // user-specified scale to preserve SPLIT_NODE = 8U, // more than one daughter TERMINAL_NODE = 16U, // no daughters MIDRANGE_NODE = 32U, // node in the middle of the branch SLOWEST_DRIFT = 64U, // node with the slowest drift speed RESERVED_LOCK = 128U, // reserved for internal use RESERVED_CNT = 256U, // reserved for internal use RESERVED_TAG = 512U, // reserved for internal use RESERVED_PASS = 1024U // reserved for internal use }; class Node { public: Node(const Cluster& jet, unsigned level, unsigned mask); inline const Cluster& getCluster() const {return jet_;} inline unsigned originalLevel() const {return originalLevel_;} inline unsigned mask() const {return nodeMask_;} inline unsigned parent() const {return parent_;} inline const SmallVector& daus() const {return daus_;} inline unsigned nDaus() const {return daus_.size();} // Non-const access to the payload inline Cluster& mutableCluster() {return jet_;} bool operator==(const Node& r) const; inline bool operator!=(const Node& r) const {return !(*this == r);} private: Node(); // Node payload Cluster jet_; // Original level number helps to speed-up various algorithms unsigned originalLevel_; // The reason why this node was included in this tree. // Should be constructed by bitwise "OR" using the enums // declared at the beginning of the public section. // Also used by several SparseClusteringTree algorithms. mutable unsigned nodeMask_; // Tree navigation. It is assumed that the node with number 0 // is the top node. It will have a dummy payload. unsigned parent_; SmallVector daus_; friend class SparseClusteringTree; }; // Default constructor builds an empty tree with a top node SparseClusteringTree(); SparseClusteringTree(const SparseClusteringTree&); virtual ~SparseClusteringTree() {} SparseClusteringTree& operator=(const SparseClusteringTree&); // The following function returns the node index which // can be used for quick lookup of node contents. Note // that the parent index must be correct. Direct daughters // of top node must have parent index 0. // // It is assumed that this tree is produced from another, // larger tree via some kind of a pruning process. We do not // anticipate the need to shrink this tree even further, // and this is why there is no "removeNode" method. If this // becomes really necessary, this tree should be cleared // completely and recreated (or use another, more appropriate // data structure). // NodeId addNode(const Node& node, NodeId parent); // Remove all nodes except the top one and forget all scales void clear(); // If you know in advance how many nodes/scales you are going // to make, use the following functions inline void reserveNodes(const unsigned n) {nodes_.reserve(n);} inline void reserveScales(const unsigned n) {scales_.reserve(n);} // Add a scale. For correct operation, all scales should be added, // starting with the level 0. void addScale(double scale); // Sort the contents of the tree by the product of cluster // magnitude and scale squared. If necessary, this sorting // should be performed after all nodes have been added. void sortNodes(); inline bool isSorted() const {return sorted_;} // Functions for retrieving the maximum level number stored inline unsigned maxStoredLevel() const {return maxLevel_;} // The following functions have the same meaning as in // AbsClusteringTree. It is assumed that the scale info // was correctly transferred from AbsClusteringTree into // this object. inline unsigned nLevels() const {return scales_.size();} inline double getScale(const unsigned level) const {return scales_.at(level);} unsigned getLevel(double scale) const; // Minimum and maximum scales. // 0 is returned in case the scales are not populated. double minScale() const; double maxScale() const; // Scales for all levels except level 0 void getAllScales(std::vector* scales, bool increasingOrder=false) const; // How many nodes are there (including the top node) inline unsigned size() const {return nodes_.size();} // Look up the node. Index 0 corresponds to the top node. inline const Node& getNode(const NodeId index) const {return nodes_.at(index);} inline const Node& uncheckedNode(const NodeId index) const {return nodes_[index];} // Non-const access to the nodes inline Node& mutableNode(const NodeId index) {return nodes_.at(index);} inline Node& uncheckedMutableNode(const NodeId index) {return nodes_[index];} // Some useful utilities bool operator==(const SparseClusteringTree& r) const; inline bool operator!=(const SparseClusteringTree& r) const {return !(*this == r);} // Number of clusters on the given level. Note that // this operation is a lot slower than the similar method // of the "AbsClusteringTree" class: the code here actually // has to drill down the branch structure and count branches // which intersect the given level. unsigned nClusters(unsigned level) const; // Number of (slowest drift) clusters on the branches which // intersect the given level and which satisfy a certain // predicate. The predicate must define the operator // "bool operator()(const Cluster&) const" which should return // "true" in case the cluster is to be included in the count. template unsigned clusterCount(unsigned level, const BooleanPredicate&) const; // Fill a vector of (slowest drift) nodes on the branches // which intersect the given level void getLevelNodes(unsigned level, std::vector* nodes) const; // Fill a vector of nodes whose associated clusters satisfy // a certain predicate. The predicate must define the operator // "bool operator()(const Cluster&) const" which should return // "true" in case the node is to be included in the output. template void getPassingNodes(unsigned level, const BooleanPredicate& pred, std::vector* nodes) const; // The following method returns the number of tree branches // which satisfy a certain predicate as a function of // level number. The slowest drift speed node on each branch // will be tested by the predicate to determine if the branch // satisfies the predicate. The occupancy vector will have // maxStoredLevel()+1 elements. By convention, occupancy of // level 0 will be set to 1. template void occupancyInScaleSpace(const BooleanPredicate& pred, std::vector* occupancy) const; // The following function finds minimum and maximum tree levels // for which the number of clusters satisfying a certain predicate // equals the desired count. If no levels have the requested number // of clusters, zeros are returned for both min and max levels. // This function returns "true" if all the levels between // minimum and maximum have the same number of clusters // (equal to requested), otherwise the function returns "false". template bool clusterCountLevels(unsigned desiredCount, const BooleanPredicate&, unsigned* minLevel, unsigned* maxLevel) const; // The following function has the same meaning as the function // with the same name in the "AbsClusteringTree" class template unsigned stableClusterCount(const BooleanPredicate& pred, unsigned* minLevel, unsigned* maxLevel, double alpha = 0.0, unsigned minStartingLevel=0, unsigned maxStartingLevel=UINT_MAX) const; // The following function determines the largest number of // (slowest drift) clusters which can be constructed in such // a way that the value of cluster magnitude times the scale // squared, MagS2, is above a certain threshold, and the // clusters satisfy the given predicate. The threshold values // corresponding to the given cluster count are returned in // the "thresholds" vector. For example, thresholds[N] == t // means that we can not construct more than N clusters with // MagS2 >= t, no matter what the branch selection is. The // parameter "maxClusterN" tells the function to stop when the // result for thresholds[maxClusterN] is obtained. The size // of the output vector will be maxClusterN + 1 or less. // // The tree must be sorted before this method is called. // template void nClustersAboveMagS2Threshold( const BooleanPredicate& pred, std::vector* thresholds, unsigned maxClusterN = UINT_MAX-1U) const; // The following function returns the set of node numbers // optimized in such a manner that, when clusters are arranged // in the decreasing order using the MagS2 variable, the cluster // with number "nOpt" has the largest possible value of MagS2. // If the tree has nodes satisfying the predicate beyond nOpt, // they will be taken from the requested level. Note that the // levels (and, therefore, the resolution scales) for the first // nOpt nodes will be chosen automatically. // // The function can also fill the "thresholds" vector, in a manner // identical to "nClustersAboveMagS2Threshold" method. The size of // the filled vector will normally be nOpt + 2. // // The tree must be sorted before this method is called. // template void getMagS2OptimalNodes(const BooleanPredicate& pred, unsigned nOpt, unsigned level, std::vector* nodes, std::vector* thresholds = 0) const; // The following function says whether node with id // id1 is an ancestor of the node with id id2 bool isAncestor(NodeId id1, NodeId id2) const; // Bad id will be returned for the parent index of the top node const NodeId badId; protected: std::vector nodes_; std::vector scales_; unsigned maxLevel_; private: bool sorted_; // Branches are identified by the index of the parent node // and by daughter number in the parent. The parent node itself // does not belong to the branch. void getLevelNodesOnBranch(unsigned level, std::vector* nodes, unsigned parentIndex, unsigned dauNumber, unsigned vetoBits = 0U) const; unsigned countLevelNodesOnBranch(unsigned level, unsigned parentIndex, unsigned dauNumber) const; void scaledMagnitudeOnBranch(unsigned parentIndex, unsigned dauNumber, double currentValue); void lockDaus(unsigned parentIndex, unsigned dauNumber, bool lock) const; unsigned relockDausAboveId(unsigned parentIndex, unsigned dauNumber, NodeId id) const; // Helper methods for magS2-related calculations template void magS2Init(const BooleanPredicate& pred, std::vector* thresholds, unsigned maxClusterN) const; // In the function below, nodeNumber is the in/out parameter. // The function returns the actual number of fills made (or // would be made if thresholds is NULL). The first "dummy" // fill is not counted. unsigned magS2Cycle(std::vector* thresholds, unsigned* nodeNumber, unsigned maxFills) const; // Check whether the daughter with number "dau" of the given parent // and/or at least one of the daughter descendants satisfy the // following conditions: // 1) "slowest drift speed" // 2) not an ancestor of the node with id "id" // 3) id number is smaller than "id" bool slowDausAboveId(unsigned parent, unsigned dau, NodeId id) const; template void branchOccupancy(const BooleanPredicate& pred, std::vector* occupancy, unsigned parentIndex, unsigned dauNumber) const; // These are essentially memory buffers for various methods std::vector swapped_; mutable std::vector inverseMap_; // The following buffer should not be modified anywhere outside // the "sortNodes" (and "scaledMagnitudeOnBranch" called by the // "sortNodes") method. This is because its contents are used // by other functions when the tree is in the sorted state. std::vector > moveMap_; }; } #include "fftjet/SparseClusteringTree.icc" #endif // FFTJET_SPARSECLUSTERINGTREE_HH_ Index: trunk/fftjet/KernelConvolver.hh =================================================================== --- trunk/fftjet/KernelConvolver.hh (revision 43) +++ trunk/fftjet/KernelConvolver.hh (revision 44) @@ -1,45 +1,45 @@ +#ifndef FFTJET_KERNELCONVOLVER_HH_ +#define FFTJET_KERNELCONVOLVER_HH_ + //========================================================================= // KernelConvolver.hh // // Class for performing kernel convolutions by FFT using a single kernel // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_KERNELCONVOLVER_HH_ -#define FFTJET_KERNELCONVOLVER_HH_ - #include "fftjet/AbsKernelConvolver.hh" #include "fftjet/AbsKernel2d.hh" namespace fftjet { template class KernelConvolver : public AbsKernelConvolver { public: // This KernelConvolver will not own the AbsFFTEngine // and AbsKernel2d objects. It is the responsibility of // the user of this class to ensure that these objects // are not destroyed before the KernelConvolver itself. // // If the "minFixBin" and "maxFixBin" arguments are // changed from their default values, it means that // the convolver will try to fix the bump reconstruction // efficiency (will apply weights to the data to simulate // the effects of boundary kernel) KernelConvolver(const AbsFFTEngine* fftEngine, const AbsKernel2d* kernel, unsigned minFixBin=0, unsigned maxFixBin=0); inline virtual ~KernelConvolver() {} protected: virtual KernelData buildKernelImage(double scale); const AbsKernel2d* const kernel; }; } #include "fftjet/KernelConvolver.icc" #endif // FFTJET_KERNELCONVOLVER_HH_ Index: trunk/fftjet/InterpolatedKernel3d.hh =================================================================== --- trunk/fftjet/InterpolatedKernel3d.hh (revision 43) +++ trunk/fftjet/InterpolatedKernel3d.hh (revision 44) @@ -1,119 +1,119 @@ +#ifndef FFTJET_INTERPOLATEDKERNEL3D_HH_ +#define FFTJET_INTERPOLATEDKERNEL3D_HH_ + //========================================================================= // InterpolatedKernel3d.hh // // Kernel whose data is tabulated on a 3d grid. In between, the data // is interpolated linearly. For scales below the range of scales covered, // the slice with the smallest scale is used. For scales above the range // the slice with the largest scale is used. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_INTERPOLATEDKERNEL3D_HH_ -#define FFTJET_INTERPOLATEDKERNEL3D_HH_ - #include #include #include "fftjet/AbsKernel2d.hh" #include "fftjet/LinearInterpolator2d.hh" namespace fftjet { class InterpolatedKernel3d : public AbsKernel2d { public: // The "scales" vector in the constructor must represent // the scales arranged in the increasing order // (not necessarily equidistant). The "data" array // can be NULL. In that case it is expected that the // data for each scale will be provided by calling the // "setScaleData" function. If the "data" array in the // constructor is not NULL, the number of elements there // should be at least scales.size()*nx*ny. It is assumed // that the array index which corresponds to the y variable // changes most often and that the index which corresponds // to the scale variable changes least often. // // The "useLogSpaceForScale" argument specifies whether // the code should use scale or log(scale) as the variable // in which the interpolation is linear. // // "nx", "xmin", and "xmax" arguments specify the binning // in the x direction. It is assumed that the first bin // is centered at xmin + (xmax - xmin)/(2 nx), just like // in a histogram. Same thing for y. // template InterpolatedKernel3d(const Real* data, const std::vector& scales, bool useLogSpaceForScale, unsigned nx, double xmin, double xmax, unsigned ny, double ymin, double ymax); virtual ~InterpolatedKernel3d() {} bool operator==(const InterpolatedKernel3d& r) const; inline bool operator!=(const InterpolatedKernel3d& r) const {return !(*this == r);} inline double xMin() const {return xmin_;} inline double xMax() const {return xmax_;} inline double yMin() const {return ymin_;} inline double yMax() const {return ymax_;} inline unsigned nx() const {return nxpoints_;} inline unsigned ny() const {return nypoints_;} inline bool usesLogSpace() const {return useLogSpace_;} inline const std::vector& scales() const {return scales_;} template void setScaleData(unsigned scaleBinNumber, const Real* data); // AbsKernel2d functions which must be implemented void setScaleRatio(double r); bool isDensity() const; double operator()(double x, double y, double scale) const; void supportRectangle(double scale, KernelSupportRectangle *r) const; void random(double r1, double r2, double scale, double* px, double* py) const; // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static InterpolatedKernel3d* read(std::istream& in); private: InterpolatedKernel3d(); void buildRandomizer(unsigned lowerBound, double scale) const; const double xmin_; const double xmax_; const double ymin_; const double ymax_; const unsigned nxpoints_; const unsigned nypoints_; const bool useLogSpace_; const std::vector scales_; std::vector interpols_; mutable LinearInterpolator2d randomizer_; mutable double randomizerScale_; // isDensity_: 1 - yes, 0 - unknown, -1 - no mutable int isDensity_; // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 3;} }; } #include "fftjet/InterpolatedKernel3d.icc" #endif // FFTJET_INTERPOLATEDKERNEL3D_HH_ Index: trunk/fftjet/ProximityClusteringTree.hh =================================================================== --- trunk/fftjet/ProximityClusteringTree.hh (revision 43) +++ trunk/fftjet/ProximityClusteringTree.hh (revision 44) @@ -1,176 +1,176 @@ +#ifndef FFTJET_PROXIMITYCLUSTERINGTREE_HH_ +#define FFTJET_PROXIMITYCLUSTERINGTREE_HH_ + //========================================================================= // ProximityClusteringTree.hh // // Clustering tree based on simple proximity matching in the eta-phi space // from the clusters at lower scale to clusters at higher scale // // I. Volobouev // April 2008 //========================================================================= -#ifndef FFTJET_PROXIMITYCLUSTERINGTREE_HH_ -#define FFTJET_PROXIMITYCLUSTERINGTREE_HH_ - #include "fftjet/SmallVector.hh" #include "fftjet/AbsClusteringTree.hh" #include "fftjet/AbsDistanceCalculator.hh" namespace fftjet { template class ProximityClusteringTree : public AbsClusteringTree { public: explicit ProximityClusteringTree( const AbsDistanceCalculator* calc, bool treeOwnsCalc=false); virtual ~ProximityClusteringTree(); // Copy constructor ProximityClusteringTree(const ProximityClusteringTree&); // Assignment operator ProximityClusteringTree& operator=(const ProximityClusteringTree&); // Inserts a new level and returns its number unsigned insert(double scale, const std::vector& clusters, const LevelInfo& levelInfo); // Change the level info void setLevelInfo(unsigned level, const LevelInfo& info); // Remove the given level. An attempt to remove level 0 // should be ignored. An attempt to remove level out of // range should result in a run-time error. void remove(unsigned level); // Reset the whole tree leaving only the root node void clear(); // Number of scale levels in the tree. Level 0 corresponds // to the root level (it will always be there). unsigned nLevels() const; // Level corresponding to the given scale unsigned getLevel(double scale) const; // Return the data which was provided during the level // creation. Calling this function on level 0 should // result in a run-time error. void getLevelData(unsigned level, double* scale, std::vector* clustersToFill, LevelInfo* levelInfo) const; // Number of clusters on the given level unsigned nClusters(unsigned level) const; // Scale corresponding to the given level double getScale(unsigned level) const; // Level-wide information (e.g., unclustered energy) provided // during the level creation const LevelInfo& getLevelInfo(unsigned level) const; // Cluster data for the given tree node. The "unchecked" // function should work faster by assuming that the cluster // with the given id exists. const Cluster& getCluster(const TreeNodeId&) const; const Cluster& uncheckedCluster(const TreeNodeId& id) const; // "Distance" between the two clusters with given ids double clusterDistance(const TreeNodeId& id1, const TreeNodeId& id2) const; // Same thing, but should work faster by assuming that // the clusters with the given ids exist double uncheckedClusterDistance(const TreeNodeId& id1, const TreeNodeId& id2) const; // The extent of the node descendants (think balltree) double nodeRadius(const TreeNodeId& id) const; double uncheckedNodeRadius(const TreeNodeId& id) const; // "Distance" to parent for the given cluster. // The smaller the distance, the better the match // between the daughter and its parent. double distanceToParent(const TreeNodeId& id) const; // Tree navigation. Note that insert and remove operations // can invalidate all node ids obtained earlier. TreeNodeId parentId(const TreeNodeId& id) const; unsigned nDaughters(const TreeNodeId& id) const; TreeNodeId daughterId(const TreeNodeId& id, unsigned idau) const; private: ProximityClusteringTree(); struct TreeNode; typedef TreeNode* TreeNodePtr; struct TreeNode { TreeNode(const TreeNodeId& id, const Cluster& jet); // Balltree construction void updateRadius(const ProximityClusteringTree& tree); // Add daughter in the right order void addDaughter(TreeNodePtr dau); // Sorting will be according to the cluster magnitude inline bool operator<(const TreeNode& r) const {return jet < r.jet;} inline bool operator>(const TreeNode& r) const {return jet > r.jet;} // node payload Cluster jet; // tree navigation TreeNodeId id; TreeNodePtr parent; double distanceToParent; double radius; SmallVector daus; private: TreeNode(); // Largest distance from this node and from all // of its subnodes to some other node double largestToNode(const ProximityClusteringTree& tree, const TreeNode* other) const; }; struct TreeLevel { TreeLevel(const TreeLevel&); TreeLevel(double scale, const LevelInfo& unclus, unsigned njets); const double scale; LevelInfo unclustered; std::vector nodes; private: TreeLevel(); TreeLevel& operator=(const TreeLevel&); }; void renumberLevels(unsigned startingLevel); void connectToParent(TreeLevel* dau, TreeLevel* parent); void disconnectFromParent(TreeLevel* dau); TreeNode* getNodeById(const TreeNodeId& id) const; void reestablishParenthood(); void calculateRadii(); std::vector levels; const AbsDistanceCalculator* distance; bool treeOwnsCalc; mutable bool recalculateRadii; }; } #include "fftjet/ProximityClusteringTree.icc" #endif // FFTJET_PROXIMITYCLUSTERINGTREE_HH_ Index: trunk/fftjet/AbsVectorRecombinationAlg.hh =================================================================== --- trunk/fftjet/AbsVectorRecombinationAlg.hh (revision 43) +++ trunk/fftjet/AbsVectorRecombinationAlg.hh (revision 44) @@ -1,113 +1,113 @@ +#ifndef FFTJET_ABSVECTORRECOMBINATIONALG_HH_ +#define FFTJET_ABSVECTORRECOMBINATIONALG_HH_ + //========================================================================= // AbsVectorRecombinationAlg.hh // // Interface class for fuzzy/crisp jet recombination algorithms in which // event energy flow is described by a collection of 4-vectors // // I. Volobouev // June 2009 //========================================================================= -#ifndef FFTJET_ABSVECTORRECOMBINATIONALG_HH_ -#define FFTJET_ABSVECTORRECOMBINATIONALG_HH_ - #include #include "fftjet/RecombinedJet.hh" namespace fftjet { template class AbsVectorRecombinationAlg { public: virtual ~AbsVectorRecombinationAlg() {} // Prototype for the main recombination algorithm function. // The derived classes must provide an implementation. Arguments // are as follows: // // peaks -- Collection of preclusters from the pattern recognition // stage (produced by the "PeakFinder" or similar class). // Implementations of this algorithm should ignore peaks // for which "membershipScaleFactor()" member function // returns 0.0. // // eventData -- A vector of elements which describe the // energy flow in the event. "VectorLike" is // a type (typically, a 4-vector) from which // it is possible to extract Et (or pt), // eta, and phi information. // // bgData -- Some representation of the average pile-up energy, // noise, background, etc. Exact interpretation of // this argument is up to concrete implementations // of this class. // // bgDataLength -- Dimensionality of the "bgData" array. Normally, // the derived classes should support the following // interpretation of this variable: // // bgDataLength == 1. This means that "bgData" is // a scalar, and the same background description // will be used for every data point in the event. // // bgDataLength == eventData.size(). This means // that "bgData" is a 1-d array. Each data point // has its own associated background. // // All other bgDataLength values should be // deemed illegal. // // jets -- Output jets. // // unclustered -- Vector sum of data elements not included // in any particular jet (produced on output). // // unused -- Scalar sum of Et (or pt) values in the data not // included in any particular jet (produced on output). // // The function should return a status word. 0 means everything // is fine. The meaning of other codes is up to concrete // implementations. // virtual int run(const std::vector& peaks, const std::vector& eventData, const BgData* bgData, unsigned bgDataLength, std::vector >* jets, VectorLike* unclustered, double* unused) = 0; // Get some info about the data last processed. // Useful for retrieving the energy correlation // matrix and the cell-to-jet assignment mask. virtual unsigned getLastNData() const = 0; virtual unsigned getLastNJets() const = 0; // The following function should return the jet energy // covariance matrix if it is available. In case it is not // available, the function should return the null pointer. // The data should remain valid at least until the next // "run" call (or until this object is destroyed). // The returned array should be dimensioned [njets+1][njets+1], // where "njets" is the number of jets returned by the // previous "run" call. C storage convention should be used // (row-major). Index 0 in this matrix should be reserved // for the unclustered energy. virtual const double* getEnergyCovarianceMatrix() const {return 0;} // The following call should return the assignments of the // input vectors to jets. In case these assignments are not // available, the function should return the null pointer. // The result should remain valid at least until the next // "run" call. The returned array should have "NData" // elements consistent with the data length from the previous // "run" call. // // clusterMask[i] should be set to the jet number // plus one. Jet numbering should be performed according to // the jet sequence returned by "run". Number 0 is reserved // for the unclustered energy. virtual const unsigned* getClusterMask() const {return 0;} }; } #endif // FFTJET_ABSVECTORRECOMBINATIONALG_HH_ Index: trunk/fftjet/AbsMembershipFunction.hh =================================================================== --- trunk/fftjet/AbsMembershipFunction.hh (revision 43) +++ trunk/fftjet/AbsMembershipFunction.hh (revision 44) @@ -1,39 +1,39 @@ +#ifndef FFTJET_ABSMEMBERSHIPFUNCTION_HH_ +#define FFTJET_ABSMEMBERSHIPFUNCTION_HH_ + //========================================================================= // AbsMembershipFunction.hh // // Interface class for jet recombination membership functions // which depend on both angular and energy variables // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ABSMEMBERSHIPFUNCTION_HH_ -#define FFTJET_ABSMEMBERSHIPFUNCTION_HH_ - #include "fftjet/ScaleSpaceKernel.hh" namespace fftjet { struct AbsMembershipFunction : public ScaleSpaceKernel { virtual ~AbsMembershipFunction() {} // The following member sets eta to phi (or x to y) scale ratio virtual void setScaleRatio(double r) = 0; // The function value virtual double operator()(double eta, double phi, double et, double scale) const = 0; // Maximum energy (Et, Pt, etc) which can be absorbed by // the jet at the given angle and scale. Override this function // if your jet recombination algorithm can meaningfully employ // the model in which the cell energy can be shared (actually // split, not just use probabilities that the cell belongs // to one jet or another) between several jets. virtual double absorbableEnergy(double /* eta */, double /* phi */, double /* scale */) const {return 0.0;} }; } #endif // FFTJET_ABSMEMBERSHIPFUNCTION_HH_ Index: trunk/fftjet/ClusteringTreeSparsifier.hh =================================================================== --- trunk/fftjet/ClusteringTreeSparsifier.hh (revision 43) +++ trunk/fftjet/ClusteringTreeSparsifier.hh (revision 44) @@ -1,86 +1,86 @@ +#ifndef FFTJET_CLUSTERINGTREESPARSIFIER_HH_ +#define FFTJET_CLUSTERINGTREESPARSIFIER_HH_ + //========================================================================= // ClusteringTreeSparsifier.hh // // This class converts a dense clustering tree (with all levels) into // a sparse clustering tree by selecting peaks with slow coordinate // drift in the scale space (thus avoiding bifurcation points). // // I. Volobouev // June 2010 //========================================================================= -#ifndef FFTJET_CLUSTERINGTREESPARSIFIER_HH_ -#define FFTJET_CLUSTERINGTREESPARSIFIER_HH_ - #include "fftjet/AbsClusteringTree.hh" #include "fftjet/SparseClusteringTree.hh" namespace fftjet { template class ClusteringTreeSparsifier { public: // Constructor arguments are as follows: // // maxLevelNumber -- Maximum level to write into the sparsified // tree. If this number is 0 or negative, this // number will be added to the maximum level // number in the input tree in order to obtain // the sparsified maximum level number. Note that // the trees we process will not necessarily // have the same number of levels (e.g., adaptive // trees). Also, the level with the highest // number may actually contain the complete set // of energy flow objects (however, this typically // does not affect the drift speed calculated for // the next-to-last level). This is why the default // value is -1 and not 0. // // filterMask -- Some bits of this mask can be set to 0. // Corresponding reasons to write out the cluster // will be ignored (per enum at the beginning of the // public section of the "SparseClusteringTree" // class). Do not mask out either MIDRANGE_NODE or // SLOWEST_DRIFT, or you will be killing the very // cluster you want to look at! // // userScales -- Scales for which all peaks should be written // out. The code will pick the nearest level. // Note, again, that these are not level numbers // because the number of levels in the input // trees is not necessarily fixed. // // nUserScales -- Number of scales in the "userScales" array. // ClusteringTreeSparsifier(int maxLevelNumber = -1, unsigned filterMask = 0xffffffffU, const double* userScales = 0, unsigned nUserScales = 0); inline virtual ~ClusteringTreeSparsifier() {} virtual void sparsify( const AbsClusteringTree& in, SparseClusteringTree* output) const; protected: void transferScaleInfo( const AbsClusteringTree& in, SparseClusteringTree* output) const; virtual void descendTreeBranch( const AbsClusteringTree& input, unsigned maxLevel, const TreeNodeId& topId, unsigned parentIndex, SparseClusteringTree* output) const; std::vector userScales_; mutable std::vector userLevels_; int maxLevelNumber_; unsigned filterMask_; }; } #include "fftjet/ClusteringTreeSparsifier.icc" #endif // FFTJET_CLUSTERINGTREESPARSIFIER_HH_ Index: trunk/fftjet/invertJetResponse.hh =================================================================== --- trunk/fftjet/invertJetResponse.hh (revision 43) +++ trunk/fftjet/invertJetResponse.hh (revision 44) @@ -1,33 +1,33 @@ +#ifndef FFTJET_INVERTJETRESPONSE_HH_ +#define FFTJET_INVERTJETRESPONSE_HH_ + //========================================================================= // invertJetResponse.hh // // The intended typical use of these functions is figuring out jet // corrections from the jet response in the form of median (mean, etc) // ratio between reconstructed and actual jet energy. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_INVERTJETRESPONSE_HH_ -#define FFTJET_INVERTJETRESPONSE_HH_ - namespace fftjet { // The following function solves numerically the equation x*f(x) = y // for unknown x. It is assumed that x*f(x) is a monotonously increasing // function, x*f(x) -> 0 when x -> 0, f(x) is positive, and that both // x and y are non-negative. template double invertJetResponse(const Functor& f, double y); // The following function solves numerically the equation x*f(a,x) = y // for unknown x, where "a" is some fixed parameter. It is assumed that // x*f(a,x) is a monotonously increasing function of x, x*f(a,x) -> 0 when // x -> 0, f(a,x) is positive, and that both x and y are non-negative. template double invertJetResponse2d(const Functor2d& f, double a, double y); } #include "fftjet/invertJetResponse.icc" #endif // FFTJET_INVERTJETRESPONSE_HH_ Index: trunk/fftjet/AbsKernel1d.hh =================================================================== --- trunk/fftjet/AbsKernel1d.hh (revision 43) +++ trunk/fftjet/AbsKernel1d.hh (revision 44) @@ -1,135 +1,135 @@ +#ifndef FFTJET_ABSKERNEL1D_HH_ +#define FFTJET_ABSKERNEL1D_HH_ + //========================================================================= // AbsKernel1d.hh // // Interface class for 1d kernel functions in scale space // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_ABSKERNEL1D_HH_ -#define FFTJET_ABSKERNEL1D_HH_ - namespace fftjet { // Utility class for kernel-related calculations class KernelSupportInterval { public: inline KernelSupportInterval(double minval, double maxval) : xmin_(minval), xmax_(maxval) {} inline double xmin() const {return xmin_;} inline double xmax() const {return xmax_;} // Multiplication by a scale factor inline KernelSupportInterval operator*(const double& scale) const { return KernelSupportInterval(xmin_*scale, xmax_*scale); } private: KernelSupportInterval(); double xmin_; double xmax_; }; // Interface class for 1d kernel functions. In all calculations, // it should be assumed that the "scale" parameter is positive. // // Subsequently, derived classes will be used with a templated // factory. In order to play nice with that factory, each kernel // class must have a constructor which looks like // // ClassName(double s, int n, const std::vector& v); // // "s" is typically (but not necessarily) the scale factor // which maps the scale argument into bandwidth values. n is normally // the power of the scale in the bandwidth calculations: 1 means // bandwidth is proportional to the scale, 0 means the bandwidth // is independent from the scale, and -1 means that bandwidth is // inversely proportional to the scale. Vector "v" is a vector of // parameters for the kernels. In addition, each class must have // a static function // // static int nParameters(); // // which defines the number of parameters the class expects to get // in "v". If this function returns a negative number, it means that // the number of parameters is arbitrary (for example, this would be // the case for kernels defined by interpolation tables). // // Of course, more complicated constructors may be needed for some // kernels. Such constructors can not be used with the common factory. // These kernels will, consequently, require special parsers in // a configuration file or in an interpretive language interface. // class AbsKernel1d { public: virtual ~AbsKernel1d() {} // The function value virtual double operator()(double x, double scale) const = 0; // The following function should return the smallest interval // which bounds the region in which the kernel value is not 0. virtual KernelSupportInterval supportInterval(double scale) const = 0; // An average over a small interval centered at x // with side dx. This is more precise (but slower) // than calculating the value of the function at // the interval center. // // Degenerate kernels (e.g., delta function) must override // the default implementation of this function. virtual double intervalAverage(double x, double sc, double dx) const; // The following function tells us whether the kernel // can be treated as a probability density. Normally, // it means that the kernel must be non-negative everywhere // and that its integral over the whole space must be positive. virtual bool isDensity() const = 0; // The following function returns random numbers distributed // according to the kernel treated as a probability density. // The argument "rnd" should be a random number between 0 and 1. // Calling this function on a kernel which is not a density // should result in a run-time error. virtual double random(double rnd, double scale) const = 0; // The following function scans the kernel for use on intervals. // The function returns the area under the scan before normalization. // It is assumed that the coordinate of the leftmost point in the // grid is xmin + xstep/2, where xstep = (xmax - xmin)/nx. // x0 is the position of the kernel center. template double scanInterval(Real* to, unsigned nx, double x0, double scale, double xmin, double xmax, double normalizationArea=1.0, bool normalize=true) const; // The following function scans the kernel for use with // periodic topologies. The function returns the // area under the scan before normalization. template double scanCircle(Real* to, unsigned nPhi, double phi0, double scale, double phiBin0Edge, double normalizationArea=1.0, bool normalize=true) const; // The following function assumes periodic 1d region // with the 2*Pi period. to[0][0] corresponds to kernel // argument (0, 0). The function returns the area under // the scan before normalization. template double scanFFT(Real* to, unsigned nPhi, double scale, double normalizationArea=1.0, bool normalize=true) const; }; } #include "fftjet/AbsKernel1d.icc" #endif // FFTJET_ABSKERNEL1D_HH_ Index: trunk/fftjet/convertPeakCoords.hh =================================================================== --- trunk/fftjet/convertPeakCoords.hh (revision 43) +++ trunk/fftjet/convertPeakCoords.hh (revision 44) @@ -1,58 +1,58 @@ +#ifndef FFTJET_CONVERTPEAKCOORDS_HH_ +#define FFTJET_CONVERTPEAKCOORDS_HH_ + //======================================================================== // convertPeakCoords.hh // // The following function will convert the peak positions from the units // of cell size into absolute eta/phi units. For input grids with // reasonable phiBin0Edge values the phi angle will end up in the // interval [0, 2*Pi). // // I. Volobouev // May 2008 //======================================================================== -#ifndef FFTJET_CONVERTPEAKCOORDS_HH_ -#define FFTJET_CONVERTPEAKCOORDS_HH_ - #include #include #include "fftjet/Grid2d.hh" #include "fftjet/Peak.hh" namespace fftjet { template void convertPeakCoords(const Grid2d& grid, std::vector* peaks) { const double twopi = 2.0*M_PI; const unsigned npeaks = peaks->size(); const double etaWidth = grid.etaBinWidth(); const double phiWidth = grid.phiBinWidth(); const double h0 = etaWidth*etaWidth; const double h1 = etaWidth*phiWidth; const double h2 = phiWidth*phiWidth; double hessian[3]; for (unsigned ijet=0; ijet= twopi) phi -= twopi; peak.setEtaPhi(grid.etaFromBinNum(peak.eta()), phi); // Convert peak Hessian peak.hessian(hessian); hessian[0] /= h0; hessian[1] /= h1; hessian[2] /= h2; peak.setHessian(hessian); } } } #endif // FFTJET_CONVERTPEAKCOORDS_HH_ Index: trunk/fftjet/AbsScalableKernel.hh =================================================================== --- trunk/fftjet/AbsScalableKernel.hh (revision 43) +++ trunk/fftjet/AbsScalableKernel.hh (revision 44) @@ -1,149 +1,149 @@ +#ifndef FFTJET_ABSSCALABLEKERNEL_HH_ +#define FFTJET_ABSSCALABLEKERNEL_HH_ + //========================================================================= // AbsScalableKernel.hh // // Interface class for scalable 2d kernel functions. Derive your kernel // from this class if only the width of your kernel depends on scale but // not the shape. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSSCALABLEKERNEL_HH_ -#define FFTJET_ABSSCALABLEKERNEL_HH_ - #include #include #include "fftjet/AbsKernel2d.hh" namespace fftjet { class AbsScalableKernel : public AbsKernel2d { public: // xScaleFactor and yScaleFactor convert the scale // into x and y bandwidth values, respectively: // x_bandwidth = xScaleFactor*pow(scale, scalePower), // y_bandwidth = yScaleFactor*pow(scale, scalePower) inline AbsScalableKernel(const double xScaleFactor, const double yScaleFactor, const int scalePower) : xs_(xScaleFactor), ys_(yScaleFactor), geomMean_(sqrt(xScaleFactor*yScaleFactor)), scalePower_(scalePower) { assert(xs_ > 0.0); assert(ys_ > 0.0); } virtual ~AbsScalableKernel() {} // The following member sets eta to phi (or x to y) scale ratio inline void setScaleRatio(const double r) { assert(r > 0.0); const double sqr = sqrt(r); xs_ = geomMean_*sqr; ys_ = geomMean_/sqr; } // Trivial accessors inline double xScaleFactor() const {return xs_;} inline double yScaleFactor() const {return ys_;} inline int scalePower() const {return scalePower_;} // The function value inline double operator()(const double x, const double y, const double scale) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwx(xs_*pscale); const double bwy(ys_*pscale); return calculate(x/bwx, y/bwy)/bwx/bwy; } // The support rectangle inline void supportRectangle(const double scale, KernelSupportRectangle *r) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } const double bwx(xs_*pscale); const double bwy(ys_*pscale); unscaledRectangle(r); r->xmin *= bwx; r->xmax *= bwx; r->ymin *= bwy; r->ymax *= bwy; } // The randomizer inline void random(const double r1, const double r2, const double scale, double* px, double* py) const { double pscale; switch (scalePower_) { case -1: pscale = 1.0/scale; break; case 0: pscale = 1.0; break; case 1: pscale = scale; break; default: pscale = pow(scale, scalePower_); break; } unscaledRandom(r1, r2, px, py); *px *= (xs_*pscale); *py *= (ys_*pscale); } private: AbsScalableKernel(); virtual double calculate(double x, double y) const = 0; virtual void unscaledRectangle(KernelSupportRectangle *r) const = 0; virtual void unscaledRandom(double r1, double r2, double* px, double* py) const = 0; double xs_; double ys_; const double geomMean_; const int scalePower_; }; } #endif // FFTJET_ABSSCALABLEKERNEL_HH_ Index: trunk/fftjet/SimplePredicates.hh =================================================================== --- trunk/fftjet/SimplePredicates.hh (revision 43) +++ trunk/fftjet/SimplePredicates.hh (revision 44) @@ -1,25 +1,25 @@ +#ifndef FFTJET_SIMPLEPREDICATES_HH_ +#define FFTJET_SIMPLEPREDICATES_HH_ + //========================================================================= // SimplePredicates.hh // // Some trivial predicates for use with AbsClusteringTree // and other classes. Usually, "SimplePeakSelector" class will // be more useful for the purpose of cluster filtering. // // I. Volobouev // February 2009 //========================================================================= -#ifndef FFTJET_SIMPLEPREDICATES_HH_ -#define FFTJET_SIMPLEPREDICATES_HH_ - #include "fftjet/SimpleFunctors.hh" namespace fftjet { template struct Always : public Functor1 { inline bool operator()(const Argument&) const {return result;} }; } #endif // FFTJET_SIMPLEPREDICATES_HH_ Index: trunk/fftjet/quartic_lib.hh =================================================================== --- trunk/fftjet/quartic_lib.hh (revision 43) +++ trunk/fftjet/quartic_lib.hh (revision 44) @@ -1,42 +1,42 @@ +#ifndef FFTJET_QUARTIC_LIB_HH_ +#define FFTJET_QUARTIC_LIB_HH_ + //========================================================================= // quartic_lib.hh // // Numerically sound routines for solving algebraic equations. // Original code by Don Herbison-Evans, with minimal adaptations for // this package by igv. See the following webpage for a description // of methods used: // // http://linus.socs.uts.edu.au/~don/pubs/solving.html // // I. Volobouev // November 2008 //========================================================================= -#ifndef FFTJET_QUARTIC_LIB_HH_ -#define FFTJET_QUARTIC_LIB_HH_ - namespace fftjet { // Numerically stable code to solve quadratic equations: // x**2 + b*x + c == 0 // Returns the number of real roots found. The roots are placed // into *x1 and *x2. The case of 0 discriminant is treated as if // there are two equal roots. int quadratic(double b, double c, double *x1, double *x2); // Find the real roots of the cubic: // x**3 + p*x**2 + q*x + r = 0 // The number of real roots is returned, and the roots are placed // into the "v3" array. int cubic(double p, double q, double r, double v3[3]); // Solve quartic equation using, hopefully, a numerically stable // scheme (see the web page for details). The number of real roots // is returned, and the roots are placed into the "rts" array. int quartic(double a, double b, double c, double d, double rts[4]); // Set the debug level for cubic and quartic. // <1 for lots of diagnostics, >5 for none. void set_quartic_debug(int level); } #endif /* FFTJET_QUARTIC_LIB_HH_ */ Index: trunk/fftjet/PeakDiffusionDistance.hh =================================================================== --- trunk/fftjet/PeakDiffusionDistance.hh (revision 43) +++ trunk/fftjet/PeakDiffusionDistance.hh (revision 44) @@ -1,42 +1,42 @@ +#ifndef FFTJET_PEAKDIFFUSIONDISTANCE_HH_ +#define FFTJET_PEAKDIFFUSIONDISTANCE_HH_ + //========================================================================= // PeakDiffusionDistance.hh // // An implementation of the AbsDistanceCalculator class for use with // ProximityClusteringTree. Employs diffusion distance described // in http://iopscience.iop.org/article/10.1088/1742-6596/404/1/012054/pdf // // I. Volobouev // September 2018 //========================================================================= -#ifndef FFTJET_PEAKDIFFUSIONDISTANCE_HH_ -#define FFTJET_PEAKDIFFUSIONDISTANCE_HH_ - #include #include #include "fftjet/AbsDistanceCalculator.hh" #include "fftjet/Peak.hh" namespace fftjet { class PeakDiffusionDistance : public AbsDistanceCalculator { public: inline explicit PeakDiffusionDistance(const double etaToPhiBandwidthRatio=1.0) : bwEta(sqrt(etaToPhiBandwidthRatio)), bwPhi(1.0/bwEta) { assert(etaToPhiBandwidthRatio > 0.0); } virtual ~PeakDiffusionDistance() {} private: double distanceBetweenClusters( double smallerScale, const Peak& smallerJet, double biggerScale, const Peak& biggerJet) const; const double bwEta; const double bwPhi; }; } #endif // FFTJET_PEAKDIFFUSIONDISTANCE_HH_ Index: trunk/fftjet/Kernels.hh =================================================================== --- trunk/fftjet/Kernels.hh (revision 43) +++ trunk/fftjet/Kernels.hh (revision 44) @@ -1,222 +1,223 @@ +#ifndef FFTJET_KERNELS_HH_ +#define FFTJET_KERNELS_HH_ + //========================================================================= // Kernels.hh // // Simple 2d kernel functions. More complicated kernels are implemented // in their own separate files. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_KERNELS_HH_ -#define FFTJET_KERNELS_HH_ - #include + #include "fftjet/AbsSymmetricKernel.hh" namespace fftjet { // Symmetric beta kernel. Constant kernel can be obtained // by using symmetric beta kernel with power 0. Epanechnikov // kernel is the symmetric beta kernel with power 1, etc. class SymmetricBeta : public AbsSymmetricKernel { public: inline SymmetricBeta(double sx, double sy, int scalePow, double power) : AbsSymmetricKernel(sx, sy, scalePow), n_(power) {} SymmetricBeta(double sx, double sy, int scalePow, const std::vector& params); virtual ~SymmetricBeta() {} inline bool isDensity() const {return true;} inline double power() const {return n_;} inline static int nParameters() {return 1;} private: double eval(double rsquared) const; double randomRadius(double rnd) const; inline double supportDistance() const {return 1.0;} double calculateResolution() const; double n_; }; // Linear kernel. The radial profile descends linearly // from a constant value at 0 to zero value at 1. class Linear2d : public AbsSymmetricKernel { public: inline Linear2d(double sx, double sy, int scalePow) : AbsSymmetricKernel(sx, sy, scalePow) {} inline Linear2d(double sx, double sy, int scalePow, const std::vector&) : AbsSymmetricKernel(sx, sy, scalePow) {} inline bool isDensity() const {return true;} inline static int nParameters() {return 0;} private: double eval(double rsquared) const; double randomRadius(double rnd) const; inline double supportDistance() const {return 1.0;} inline double calculateHalfWeightRadius() const {return 0.5;} inline double calculateResolution() const {return 1.0;} }; // Gaussian kernel. Standard deviation is 1 in any direction. class Gauss2d : public AbsSymmetricKernel { public: inline Gauss2d(double sx, double sy, int scalePow) : AbsSymmetricKernel(sx, sy, scalePow) {} inline Gauss2d(double sx, double sy, int scalePow, const std::vector&) : AbsSymmetricKernel(sx, sy, scalePow) {} inline bool isDensity() const {return true;} double axialWeight(double r, unsigned nsteps=0) const; inline static int nParameters() {return 0;} private: double eval(double rsquared) const; double randomRadius(double rnd) const; double supportDistance() const; inline double calculateHalfWeightRadius() const {return 1.177410022515474691;} inline double calculateResolution() const {return 1.414213562373095;} }; // The following functor, together with RealFrequencyKernel class, // can be used to calculate stable rotationally symmetric Sub-Gaussian // distributions. If you need to know what they are see, for example, // section 6 (page 40) of // http://academic2.american.edu/~jpnolan/stable/DataAnalysis.ps class SubGauss : public AbsSymmetricKernel { public: SubGauss(double sx, double sy, int scalePow, double alpha); SubGauss(double sx, double sy, int scalePow, const std::vector& params); virtual ~SubGauss() {} inline bool isDensity() const {return true;} inline double alpha() const {return alpha_;} inline static int nParameters() {return 1;} private: double eval(double rsquared) const; double supportDistance() const; double randomRadius(double rnd) const; double alpha_; double normfactor_; static double calculateNormfactor(double alpha); }; // Huber kernel. Weight of the tail must be >= 0 and < 1. class Huber2d : public AbsSymmetricKernel { public: Huber2d(double sx, double sy, int scalePow, double tailWeight); Huber2d(double sx, double sy, int scalePow, const std::vector& params); virtual ~Huber2d() {} inline bool isDensity() const {return true;} inline double tailWeight() const {return tWeight_;} inline double transition() const {return a_;} inline static int nParameters() {return 1;} private: double eval(double rsquared) const; double randomRadius(double rnd) const; double supportDistance() const; double calculateHalfWeightRadius() const; void setup(); double weight_(double asquared) const; double tWeight_; double a_; double normfactor_; }; class DeltaFunctionKernel : public AbsKernel2d { public: inline explicit DeltaFunctionKernel(double volume) : value_(volume) {} DeltaFunctionKernel(double sx, double sy, int scalePow, const std::vector& params); virtual ~DeltaFunctionKernel() {} inline void setScaleRatio(double) {} inline bool isDensity() const {return true;} inline double volume() const {return value_;} double operator()(double x, double y, double scale) const; void supportRectangle(double scale, KernelSupportRectangle *r) const; double rectangleAverage(double x, double y, double scale, double dx, double dy) const; void random(double r1, double r2, double scale, double* px, double* py) const; inline static int nParameters() {return 1;} private: double value_; }; // Inverse power kernel proportional to 1/(1 + (r^2)^n). // n must be >= 1.05. For n <= 1 the integral diverges, // and for n between 1 and 1.05 we encounter various // numerical difficulties. class InvPower2d : public AbsSymmetricKernel { public: InvPower2d(double sx, double sy, int scalePow, double n); InvPower2d(double sx, double sy, int scalePow, const std::vector& params); inline virtual ~InvPower2d() {delete [] cdf_;} inline bool isDensity() const {return true;} inline double power() const {return n_;} double axialWeight(double r, unsigned nsteps=0) const; inline static int nParameters() {return 1;} private: InvPower2d(const InvPower2d&); InvPower2d& operator=(const InvPower2d&); inline double calculateResolution() const {return 1.0;} double eval(double rsquared) const; double randomRadius(double rnd) const; double supportDistance() const; void setup(); void buildCdf(); double tailIntegral(double r) const; double headIntegral(double r) const; double annulusIntegral(double middleR, double width) const; double calculateHalfWeightRadius() const; double n_; double normfactor_; double lowbound_; double lowweight_; double upbound_; double upweight_; double* cdf_; }; } #endif // FFTJET_KERNELS_HH_ Index: trunk/fftjet/KernelData.hh =================================================================== --- trunk/fftjet/KernelData.hh (revision 43) +++ trunk/fftjet/KernelData.hh (revision 44) @@ -1,46 +1,46 @@ +#ifndef FFTJET_KERNELDATA_HH_ +#define FFTJET_KERNELDATA_HH_ + //========================================================================= // KernelData.hh // // Scan information for a single scale. This class will not own // the pointers. "normalization()" method returns the normalization // curve (or NULL pointer in case there is no normalization data). // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_KERNELDATA_HH_ -#define FFTJET_KERNELDATA_HH_ - namespace fftjet { template class KernelData { public: // The "fftScanArea" parameter is not the area under the scan // curve (which is usually normalized). Instead, it should // represent the area returned by the "scanFFT" function // of the kernel. If it is very different from 1.0, it means // there are significant distortions of the kernel shape // due to binning. KernelData(const Real* scan, const Complex* image, const Real* norm, double fftScanArea); inline const Real* scan() const {return scan_;} inline const Complex* image() const {return image_;} inline const Real* normalization() const {return normalization_;} inline double fftScanArea() const {return scanArea_;} private: KernelData(); const Real* scan_; const Complex* image_; const Real* normalization_; double scanArea_; }; } #include "fftjet/KernelData.icc" #endif // FFTJET_KERNELDATA_HH_ Index: trunk/fftjet/PeakEtaDependentDistance.hh =================================================================== --- trunk/fftjet/PeakEtaDependentDistance.hh (revision 43) +++ trunk/fftjet/PeakEtaDependentDistance.hh (revision 44) @@ -1,42 +1,42 @@ +#ifndef FFTJET_PEAKETADEPENDENTDISTANCE_HH_ +#define FFTJET_PEAKETADEPENDENTDISTANCE_HH_ + //========================================================================= // PeakEtaDependentDistance.hh // // An implementation of the AbsDistanceCalculator class similar to // PeakEtaPhiDistance but in which the ratio between eta and phi // bandwidth values is given by an eta-dependent interpolation table. // Note that it is possible to (mis)specify the interpolation table // in such a way that this function is no longer a pseudo-metric. // Therefore, use at your own risk. // // I. Volobouev // June 2010 //========================================================================= -#ifndef FFTJET_PEAKETADEPENDENTDISTANCE_HH_ -#define FFTJET_PEAKETADEPENDENTDISTANCE_HH_ - #include "fftjet/AbsDistanceCalculator.hh" #include "fftjet/LinearInterpolator1d.hh" #include "fftjet/Peak.hh" namespace fftjet { class PeakEtaDependentDistance : public AbsDistanceCalculator { public: explicit PeakEtaDependentDistance( const LinearInterpolator1d& tableOfEtaToPhiBandwidthRatios); virtual ~PeakEtaDependentDistance() {} private: PeakEtaDependentDistance(); double distanceBetweenClusters( double smallerScale, const Peak& smallerJet, double biggerScale, const Peak& biggerJet) const; LinearInterpolator1d table_; }; } #endif // FFTJET_PEAKETADEPENDENTDISTANCE_HH_ Index: trunk/fftjet/KernelRecombinationAlg.hh =================================================================== --- trunk/fftjet/KernelRecombinationAlg.hh (revision 43) +++ trunk/fftjet/KernelRecombinationAlg.hh (revision 44) @@ -1,253 +1,253 @@ +#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 //========================================================================= -#ifndef FFTJET_KERNELRECOMBINATIONALG_HH_ -#define FFTJET_KERNELRECOMBINATIONALG_HH_ - #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* energyVariance; double* weightSum; 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_ Index: trunk/fftjet/CUFFTFloatEngine.hh =================================================================== --- trunk/fftjet/CUFFTFloatEngine.hh (revision 43) +++ trunk/fftjet/CUFFTFloatEngine.hh (revision 44) @@ -1,68 +1,68 @@ +#ifndef FFTJET_CUFFTFLOATENGINE_HH_ +#define FFTJET_CUFFTFLOATENGINE_HH_ + //========================================================================= // CUFFTFloatEngine.hh // // Single precision FFT engine for jet reconstruction which uses // CUDA-based FFT library from NVIDIA. For more details, see // http://en.wikipedia.org/wiki/CUDA // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_CUFFTFLOATENGINE_HH_ -#define FFTJET_CUFFTFLOATENGINE_HH_ - #include "fftjet/AbsFFTEngine.hh" #include "cufft.h" namespace fftjet { class CUFFTFloatEngine : public AbsFFTEngine { public: CUFFTFloatEngine(unsigned nEta, unsigned nPhi); virtual ~CUFFTFloatEngine(); void transformForward(const float* from, cufftComplex* to) const; void transformBack(const cufftComplex* from, float* to) const; void multiplyTransforms(const cufftComplex* a, const cufftComplex* b, cufftComplex* result) const; void addTransforms(const cufftComplex* a, const cufftComplex* b, cufftComplex* result) const; void divideTransforms(const cufftComplex* a, const cufftComplex* b, cufftComplex* result) const; void deconvolutionRatio(const cufftComplex* a, const cufftComplex* b, cufftComplex* result) const; void subtractTransforms(const cufftComplex* a, const cufftComplex* b, cufftComplex* result) const; void scaleTransform(const cufftComplex* a, float scale, cufftComplex* result) const; void amplitudeSquared(const cufftComplex* a, cufftComplex* result) const; double totalPower(const cufftComplex* a) const; void setToReal(cufftComplex* a, float value) const; void setTransformPoint(cufftComplex* points, unsigned iEta, unsigned iPhi, const std::complex& value) const; std::complex getTransformPoint(const cufftComplex* points, unsigned iEta, unsigned iPhi) const; cufftComplex* allocateComplex() const; void destroyComplex(cufftComplex*) const; private: CUFFTFloatEngine(const CUFFTFloatEngine&); CUFFTFloatEngine& operator=(const CUFFTFloatEngine&); void checkFFTStatus(cufftResult status) const; float *in; cufftComplex *out; cufftHandle pf; cufftHandle pb; }; } #endif // FFTJET_CUFFTFLOATENGINE_HH_ Index: trunk/fftjet/PythiaKernel_30_100_v0.hh =================================================================== --- trunk/fftjet/PythiaKernel_30_100_v0.hh (revision 43) +++ trunk/fftjet/PythiaKernel_30_100_v0.hh (revision 44) @@ -1,30 +1,30 @@ +#ifndef FFTJET_PYTHIAKERNEL_30_100_V0_HH_ +#define FFTJET_PYTHIAKERNEL_30_100_V0_HH_ + //========================================================================= // PythiaKernel_30_100_v0.hh // // Angular energy distribution for jets produced by light quarks // with momenta between 30 and 100 GeV/c, generated using Pythia // jet gun. The natural scale for this kernel is 1.0/(parton Pt), // (if scalePower == 1) where parton Pt is in units GeV/c. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_PYTHIAKERNEL_30_100_V0_HH_ -#define FFTJET_PYTHIAKERNEL_30_100_V0_HH_ - #include "fftjet/LogProfileKernel.hh" namespace fftjet { class PythiaKernel_30_100_v0 : public LogProfileKernel { public: explicit PythiaKernel_30_100_v0(int scalePower = 1); PythiaKernel_30_100_v0(double sx, double sy, int scalePower, const std::vector& params); inline static int nParameters() {return 0;} }; } #endif // FFTJET_PYTHIAKERNEL_30_100_V0_HH_ Index: trunk/fftjet/FFTWDoubleEngine.hh =================================================================== --- trunk/fftjet/FFTWDoubleEngine.hh (revision 43) +++ trunk/fftjet/FFTWDoubleEngine.hh (revision 44) @@ -1,29 +1,29 @@ +#ifndef FFTJET_FFTWDOUBLEENGINE_HH_ +#define FFTJET_FFTWDOUBLEENGINE_HH_ + //========================================================================= // FFTWDoubleEngine.hh // // Double precision FFT engine for jet reconstruction which uses // the FFTW library from http://www.fftw.org // // I. Volobouev // January 2008 //========================================================================= -#ifndef FFTJET_FFTWDOUBLEENGINE_HH_ -#define FFTJET_FFTWDOUBLEENGINE_HH_ - #include "fftjet/FFTWEngine.hh" #include "fftw3.h" namespace fftjet { class FFTWDoubleEngine : public FFTWEngine { public: FFTWDoubleEngine(unsigned nEta, unsigned nPhi, unsigned optimization = FFTW_PATIENT); virtual ~FFTWDoubleEngine(); }; } #include "fftjet/FFTWDoubleEngine.icc" #endif // FFTJET_FFTWDOUBLEENGINE_HH_ Index: trunk/fftjet/Kernels1d.hh =================================================================== --- trunk/fftjet/Kernels1d.hh (revision 43) +++ trunk/fftjet/Kernels1d.hh (revision 44) @@ -1,84 +1,85 @@ +#ifndef FFTJET_KERNELS1D_HH_ +#define FFTJET_KERNELS1D_HH_ + //========================================================================= // Kernels1d.hh // // Simple 1d kernel functions. More complicated kernels are implemented // in their own separate files. // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_KERNELS1D_HH_ -#define FFTJET_KERNELS1D_HH_ - #include + #include "fftjet/AbsScalableKernel1d.hh" namespace fftjet { // Gaussian kernel class Gauss1d : public AbsScalableKernel1d { public: inline Gauss1d(double sx, int scalePow) : AbsScalableKernel1d(sx, scalePow) {} inline Gauss1d(double sx, int scalePow, const std::vector&) : AbsScalableKernel1d(sx, scalePow) {} inline bool isDensity() const {return true;} inline static int nParameters() {return 0;} private: double calculate(double x) const; KernelSupportInterval unscaledInterval() const; double unscaledRandom(double r1) const; }; // Symmetric beta kernel class SymmetricBeta1d : public AbsScalableKernel1d { public: SymmetricBeta1d(double sx, int scalePow, double power); SymmetricBeta1d(double sx, int scalePow, const std::vector& p); virtual ~SymmetricBeta1d() {} inline bool isDensity() const {return true;} inline double power() const {return n_;} inline static int nParameters() {return 1;} private: double calculate(double x) const; KernelSupportInterval unscaledInterval() const; double unscaledRandom(double r1) const; double calculateNorm() const; double n_; double norm_; }; class DeltaFunction1d : public AbsKernel1d { public: inline explicit DeltaFunction1d(double area) : value_(area) {} DeltaFunction1d(double sx, int scalePow, const std::vector& params); virtual ~DeltaFunction1d() {} inline double area() const {return value_;} double operator()(double x, double scale) const; KernelSupportInterval supportInterval(double scale) const; double intervalAverage(double x, double sc, double dx) const; inline bool isDensity() const {return true;} double random(double rnd, double scale) const; inline static int nParameters() {return 1;} private: double value_; }; } #endif // FFTJET_KERNELS1D_HH_ Index: trunk/fftjet/peakEtLifetime.hh =================================================================== --- trunk/fftjet/peakEtLifetime.hh (revision 43) +++ trunk/fftjet/peakEtLifetime.hh (revision 44) @@ -1,31 +1,31 @@ +#ifndef FFTJET_PEAKETLIFETIME_HH_ +#define FFTJET_PEAKETLIFETIME_HH_ + //========================================================================= // peakEtLifetime.hh // // Cluster lifetime weighted by the local cluster Et fraction. // Intended for use with sparse clustering trees. // // I. Volobouev // May 2013 //========================================================================= -#ifndef FFTJET_PEAKETLIFETIME_HH_ -#define FFTJET_PEAKETLIFETIME_HH_ - namespace fftjet { template double peakEtSplitTime(const SparseTree& tree, typename SparseTree::NodeId id, double minScale); template double peakEtMergeTime(const SparseTree& tree, typename SparseTree::NodeId id, double maxScale); template void updateSplitMergeTimes(SparseTree& tree, double minScale, double maxScale); } #include "fftjet/peakEtLifetime.icc" #endif // FFTJET_PEAKETLIFETIME_HH_ Index: trunk/fftjet/AbsJetCorrector.hh =================================================================== --- trunk/fftjet/AbsJetCorrector.hh (revision 43) +++ trunk/fftjet/AbsJetCorrector.hh (revision 44) @@ -1,51 +1,51 @@ +#ifndef FFTJET_ABSJETCORRECTOR_HH_ +#define FFTJET_ABSJETCORRECTOR_HH_ + //========================================================================= // AbsJetCorrector.hh // // Interface class for jet correction-type calculations. It is assumed // that we do not want to adjust the jet location, only its magnitude. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ABSJETCORRECTOR_HH_ -#define FFTJET_ABSJETCORRECTOR_HH_ - #include namespace fftjet { template struct AbsJetCorrector { virtual ~AbsJetCorrector() {} // It is assumed that "magnitude" is some kind of a guess // for the jet energy-like variable (Et, Pt, etc) which we // want to correct. The Jet class should provide all // other info about the jet (eta, phi, scale, etc). // The function should return adjusted magnitude. virtual double operator()(const double& magnitude, const Jet& jet) const = 0; }; // The following class does not perform any correction template struct IdleJetCorrector : public AbsJetCorrector { inline double operator()(const double& magnitude, const Jet&) const {return magnitude;} }; // The following can sometimes be used for debugging purposes template struct InvalidJetCorrector : public AbsJetCorrector { inline double operator()(const double&, const Jet&) const { assert(!"InvalidJetCorrector::operator() triggered"); return 0.0; } }; } #endif // FFTJET_ABSJETCORRECTOR_HH_ Index: trunk/fftjet/RealFrequencyKernel1d.hh =================================================================== --- trunk/fftjet/RealFrequencyKernel1d.hh (revision 43) +++ trunk/fftjet/RealFrequencyKernel1d.hh (revision 44) @@ -1,50 +1,50 @@ +#ifndef FFTJET_REALFREQUENCYKERNEL1D_HH_ +#define FFTJET_REALFREQUENCYKERNEL1D_HH_ + //========================================================================= // RealFrequencyKernel1d.hh // // 1d functions in the Fourier transform space whose transform // is pure real. Implemented via real kernels. // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_REALFREQUENCYKERNEL1D_HH_ -#define FFTJET_REALFREQUENCYKERNEL1D_HH_ - #include #include "fftjet/AbsFrequencyKernel1d.hh" #include "fftjet/AbsKernel1d.hh" namespace fftjet { class RealFrequencyKernel1d : public AbsFrequencyKernel1d { public: inline RealFrequencyKernel1d(const AbsKernel1d* realKernel, const bool assumePointerOwnership=false) : kern_(realKernel), owns_(assumePointerOwnership) {assert(kern_);} inline virtual ~RealFrequencyKernel1d() { if (owns_) delete const_cast(kern_); } inline bool ownsPointer() const {return owns_;} inline std::complex operator()( const int ix, const double scale) const { // It is natural to invert the scale in the frequency space return std::complex((*kern_)(ix, 1.0/scale), 0.0); } private: RealFrequencyKernel1d(); RealFrequencyKernel1d(const RealFrequencyKernel1d&); RealFrequencyKernel1d& operator=(const RealFrequencyKernel1d&); const AbsKernel1d* const kern_; const bool owns_; }; } #endif // FFTJET_REALFREQUENCYKERNEL1D_HH_ Index: trunk/fftjet/KernelVectorRecombinationAlg.hh =================================================================== --- trunk/fftjet/KernelVectorRecombinationAlg.hh (revision 43) +++ trunk/fftjet/KernelVectorRecombinationAlg.hh (revision 44) @@ -1,226 +1,226 @@ +#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 //========================================================================= -#ifndef FFTJET_KERNELVECTORRECOMBINATIONALG_HH_ -#define FFTJET_KERNELVECTORRECOMBINATIONALG_HH_ - #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* energyVariance; double* weightSum; 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/LinearInterpolator1d.hh =================================================================== --- trunk/fftjet/LinearInterpolator1d.hh (revision 43) +++ trunk/fftjet/LinearInterpolator1d.hh (revision 44) @@ -1,122 +1,122 @@ +#ifndef FFTJET_LINEARINTERPOLATOR1D_HH_ +#define FFTJET_LINEARINTERPOLATOR1D_HH_ + //========================================================================= // LinearInterpolator1d.hh // // This class can be used to interpolate histograms or similar regularly // spaced data in one dimension // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_LINEARINTERPOLATOR1D_HH_ -#define FFTJET_LINEARINTERPOLATOR1D_HH_ - #include #include #include #include "fftjet/SimpleFunctors.hh" namespace fftjet { class LinearInterpolator1d : public Functor1 { public: // Constructor from a regularly spaced data. "f_low" // is the function value below x_min, and "f_hi" is // the function value above "x_max". template LinearInterpolator1d(const Real* data, unsigned npoints, double x_min, double x_max, double f_low=0.0, double f_hi=0.0); // Constructor from a list of points, not necessarily regularly // spaced. The first member of the pair is the x coordinate, the // second is the tabulated function value. The input list will be // interpolated to "npoints" internal points linearly. LinearInterpolator1d(const std::vector >& v, unsigned npoints); // Constructor which builds a function returning the given constant explicit LinearInterpolator1d(double c); // Copy constructor LinearInterpolator1d(const LinearInterpolator1d&); virtual ~LinearInterpolator1d(); LinearInterpolator1d& operator=(const LinearInterpolator1d&); bool operator==(const LinearInterpolator1d& r) const; inline bool operator!=(const LinearInterpolator1d& r) const {return !(*this == r);} inline unsigned nx() const {return npoints;} inline double xMin() const {return xmin;} inline double xMax() const {return xmin + npoints*binwidth;} inline double fLow() const {return flow;} inline double fHigh() const {return fhi;} inline const double* getData() const {return data;} virtual double operator()(const double& x) const; // The following function tells whether the function // values are non-negative (and, therefore, the function // can be treated as a density) bool isNonNegative() const; // The following function returns the function integral // inside the interval double integral() const; // The following function returns the value of // the integral as if the interpolator is treated // as a density between xmin and xmax double getCdf(double x) const; // The following function can be used to generate // random numbers according to the function represented // by the interpolator, assuming that the interpolated // function can be used as a density double random(double rnd) const; // The following function can be used to normalize the // function integral inside the interval to the given value void normalize(double value); // The "write" function returns "true" if the object // is successfully written out. bool write(std::ostream& of) const; // The following function reads the object in. // Returns NULL in case of a problem. static LinearInterpolator1d* read(std::istream& in); // Helper function for generating random numbers // inside a single cell static double singleCellRandom(double xlo, double bwx, double z00, double z01, double rnd); private: LinearInterpolator1d(); void buildCdf(); double* data; double* cdf; double xmin; double flow; double fhi; double binwidth; unsigned npoints; // Version number for the I/O static inline unsigned version() {return 1;} // Class id for the I/O static inline unsigned classId() {return 6;} }; } #include "fftjet/LinearInterpolator1d.icc" #endif // FFTJET_LINEARINTERPOLATOR1D_HH_ Index: trunk/fftjet/AbsRecombinationAlg.hh =================================================================== --- trunk/fftjet/AbsRecombinationAlg.hh (revision 43) +++ trunk/fftjet/AbsRecombinationAlg.hh (revision 44) @@ -1,130 +1,130 @@ +#ifndef FFTJET_ABSRECOMBINATIONALG_HH_ +#define FFTJET_ABSRECOMBINATIONALG_HH_ + //========================================================================= // AbsRecombinationAlg.hh // // Interface class for fuzzy/crisp jet recombination algorithms // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSRECOMBINATIONALG_HH_ -#define FFTJET_ABSRECOMBINATIONALG_HH_ - #include #include "fftjet/Grid2d.hh" #include "fftjet/RecombinedJet.hh" namespace fftjet { template class AbsRecombinationAlg { public: virtual ~AbsRecombinationAlg() {} // Prototype for the main recombination algorithm function. // The derived classes must provide an implementation. Arguments // are as follows: // // peaks -- Collection of preclusters from the pattern recognition // stage (produced by the "PeakFinder" or similar class). // Implementations of this algorithm should ignore peaks // for which "membershipScaleFactor()" member function // returns 0.0. // // eventData -- Actual discretized data from the current event. // Usually this is the same Grid2d object as the one // used for pattern recognition. // // bgData -- Some representation of the average pile-up energy, // noise, background, etc. Exact interpretation of // this argument is up to concrete implementations // of this class. // // nBgEta, -- Dimensionalities of the "bgData" array. Normally, // nBgPhi the derived classes should support the following // interpretation of these variables: // // nBgEta == 1, nBgPhi == 1. This means that "bgData" is // a scalar, and the same background description will be // used for each grid cell. // // nBgEta == eventData.nEta(), nBgPhi == 1. This means // that "bgData" is a 1-d array. The same background // description will be used for each phi cell number // in each given eta, but backgrounds are different // for different eta. // // nBgEta == eventData.nEta(), nBgPhi == eventData.nPhi() // This means that "bgData" is a 2-d array, and there is // a separate background description object for every // grid cell. // // All other nBgEta, nBgPhi combinations should be // deemed illegal. The "bgDimensionality" method // can be used to enforce this policy. // // jets -- Output jets. // // unclustered -- Vector sum of grid cells not included // in any particular jet (produced on output). // // unused -- Scalar sum of values in the event grid not included // in any particular jet (produced on output). // // The function should return a status word. 0 means everything // is fine. The meaning of other codes is up to concrete // implementations. // virtual int run(const std::vector& peaks, const Grid2d& eventData, const BgData* bgData, unsigned nBgEta, unsigned nBgPhi, std::vector >* jets, VectorLike* unclustered, double* unused) = 0; // Get some info about the data last processed. // Useful for retrieving the energy correlation // matrix and the cell-to-jet assignment mask. virtual unsigned getLastNEta() const = 0; virtual unsigned getLastNPhi() const = 0; virtual unsigned getLastNJets() const = 0; // The following function should return the jet energy // covariance matrix if it is available. In case it is not // available, the function should return the null pointer. // The data should remain valid at least until the next // "run" call (or until this object is destroyed). // The returned array should be dimensioned [njets+1][njets+1], // where "njets" is the number of jets returned by the // previous "run" call. C storage convention should be used // (row-major). Index 0 in this matrix should be reserved // for the unclustered energy. virtual const double* getEnergyCovarianceMatrix() const {return 0;} // The following call should return the cell-to-jet mask // if it is available. In case it is not available, // the function should return the null pointer. The data // should remain valid at least until the next "run" call. // The returned array should be dimensioned [nEta][nPhi] // with sizes corresponding to the data grid size from // the previous "run" call. // // clusterMask[ieta][iphi] should be set to the jet number // plus one. Jet numbering should be performed according to // the jet sequence returned by "run". Number 0 is reserved // for the unclustered energy. virtual const unsigned* getClusterMask() const {return 0;} // The following function will checks the "nBgEta" and // "nBgPhi" arguments and will return the effective // dimensionality of the "bgData" array (0, 1, or 2). // A run-time error will be generated in case nBgEta, nBgPhi // arguments are inconsistent with the Grid2d dimensions. static unsigned bgDimensionality(const Grid2d& eventData, unsigned nBgEta, unsigned nBgPhi); }; } #include "fftjet/AbsRecombinationAlg.icc" #endif // FFTJET_ABSRECOMBINATIONALG_HH_ Index: trunk/fftjet/PhiKernels.hh =================================================================== --- trunk/fftjet/PhiKernels.hh (revision 43) +++ trunk/fftjet/PhiKernels.hh (revision 44) @@ -1,71 +1,71 @@ +#ifndef FFTJET_PHIKERNELS_HH_ +#define FFTJET_PHIKERNELS_HH_ + //========================================================================= // PhiKernels.hh // // Kernel functions in phi // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_PHIKERNELS_HH_ -#define FFTJET_PHIKERNELS_HH_ - #include #include "fftjet/AbsScalablePhiKernel.hh" namespace fftjet { class PhiGauss : public AbsScalablePhiKernel { public: inline PhiGauss(const double sy, const int scalePow) : AbsScalablePhiKernel(sy, scalePow) {} inline PhiGauss(double, double sy, int scalePow, const std::vector&) : AbsScalablePhiKernel(sy, scalePow) {} inline bool isDensity() const {return true;} inline static int nParameters() {return 0;} private: double unscaledPhiFcn(double phi) const; void unscaledPhiSupport(double *phimin, double *phimax) const; double unscaledPhiRandom(double rnd) const; }; // Interpolation kernel for smearing in the phi direction only. // It is assumed that the kernel is symmetric w.r.t. the change of // sign of the phi angle, and the values are calculated on the phi // interval [0, Pi/2] with uniform step in phi. The first value // corresponds to r = 0.0, and the last one corresponds to // Pi/2*(1.0 - 1.0/profile.size()). Each profile value must be // non-negative. class PhiProfileKernel : public AbsScalablePhiKernel { public: PhiProfileKernel(double sx, double sy, int scalePower, const std::vector& profile); inline virtual ~PhiProfileKernel() {delete [] cdf_;} inline bool isDensity() const {return true;} inline static int nParameters() {return -1;} private: PhiProfileKernel(const PhiProfileKernel&); PhiProfileKernel& operator=(const PhiProfileKernel&); double unscaledPhiFcn(double phi) const; void unscaledPhiSupport(double *phimin, double *phimax) const; double unscaledPhiRandom(double rnd) const; unsigned rCellNumber(const double r, double *diff) const; double randomDistance(double rnd) const; const std::vector params_; double* cdf_; double normfactor_; }; } #endif // FFTJET_PHIKERNELS_HH_ Index: trunk/fftjet/Matrix.hh =================================================================== --- trunk/fftjet/Matrix.hh (revision 43) +++ trunk/fftjet/Matrix.hh (revision 44) @@ -1,103 +1,103 @@ +#ifndef FFTJET_MATRIX_HH_ +#define FFTJET_MATRIX_HH_ + //========================================================================= // Matrix.hh // // A simple helper class for matrix manipulations. Depending on how much // space is provided with the "Len" parameter, the data will be placed // either on the stack or on the heap. // // C storage convention is used for internal data. In the element access // operations, array bounds are not checked. // // Note that this class is slower and less complete than dedicated matrix // classes based on expression templates (such as those in boost uBLAS // or in Blitz++). Don't use it for any calculations in which speed is // really important. // // I. Volobouev // November 2008 //========================================================================= -#ifndef FFTJET_MATRIX_HH_ -#define FFTJET_MATRIX_HH_ - #include namespace fftjet { template class Matrix { public: // Default constructor creates an unitialized matrix // which can be assigned from other matrix Matrix(); // The following constructor creates an unitialized matrix // which can be assigned element-by-element or from another // matrix with the same dimensions Matrix(unsigned nrows, unsigned ncols); // The following constructor initializes the matrix as follows: // initCode = 0 All elements are initialized to 0 // initCode = 1 Matrix must be square; diagonal elements // are initialized to 1 Matrix(unsigned nrows, unsigned ncols, int initCode); // The following constructor initializes the matrix from the // given 1d array using C storage conventions Matrix(unsigned nrows, unsigned ncols, const Numeric* data); Matrix(const Matrix&); ~Matrix(); Matrix& operator=(const Matrix&); inline unsigned nrows() const {return nrows_;} inline unsigned ncols() const {return ncols_;} inline Numeric* data() const {return data_;} // The following function resets the object to an unintialized state void reset(); // The following function changes the object dimensions. // All data is lost in the process. void resize(unsigned nrows, unsigned ncols); // The following function sets all elements to 0 void zeroOut(); bool operator==(const Matrix&) const; bool operator!=(const Matrix&) const; Numeric* operator[](unsigned) const; Matrix operator*(const Matrix& r) const; Matrix operator*(Numeric r) const; Matrix operator/(Numeric r) const; Matrix operator+(const Matrix& r) const; Matrix operator-(const Matrix& r) const; Matrix operator+() const; Matrix operator-() const; Matrix& operator*=(Numeric r); Matrix& operator/=(Numeric r); Matrix& operator+=(const Matrix& r); Matrix& operator-=(const Matrix& r); Matrix T() const; Matrix symmetrize() const; private: Numeric local_[Len]; Numeric* data_; unsigned nrows_; unsigned ncols_; unsigned len_; }; } template std::ostream& operator<<(std::ostream& os, const fftjet::Matrix& m); #include "fftjet/Matrix.icc" #endif // FFTJET_MATRIX_HH_ Index: trunk/fftjet/Peak.cc =================================================================== --- trunk/fftjet/Peak.cc (revision 43) +++ trunk/fftjet/Peak.cc (revision 44) @@ -1,102 +1,125 @@ +#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), 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_ && 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 sense here. + // 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/EtSumVectorRecombinationAlg.hh =================================================================== --- trunk/fftjet/EtSumVectorRecombinationAlg.hh (revision 43) +++ trunk/fftjet/EtSumVectorRecombinationAlg.hh (revision 44) @@ -1,62 +1,62 @@ +#ifndef FFTJET_ETSUMVECTORRECOMBINATIONALG_HH_ +#define FFTJET_ETSUMVECTORRECOMBINATIONALG_HH_ + //========================================================================= // EtSumVectorRecombinationAlg.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; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the precluster. // // I. Volobouev // June 2009 //========================================================================= -#ifndef FFTJET_ETSUMVECTORRECOMBINATIONALG_HH_ -#define FFTJET_ETSUMVECTORRECOMBINATIONALG_HH_ - #include "fftjet/KernelVectorRecombinationAlg.hh" namespace fftjet { template class EtSumVectorRecombinationAlg : public KernelVectorRecombinationAlg { public: typedef double (VectorLike::* VectorLikeMemberFunction)() const; inline EtSumVectorRecombinationAlg( ScaleSpaceKernel* kernel, VectorLikeMemberFunction etFcn, VectorLikeMemberFunction etaFcn, VectorLikeMemberFunction phiFcn, const Functor2* bgWeight, double unlikelyBgWeight, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask) : KernelVectorRecombinationAlg( kernel, etFcn, etaFcn, phiFcn, bgWeight, unlikelyBgWeight, winnerTakesAll, buildCorrelationMatrix, buildClusterMask) {} inline virtual ~EtSumVectorRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} inline virtual bool useEtCentroid() {return false;} }; } #endif // FFTJET_ETSUMVECTORRECOMBINATIONALG_HH_ Index: trunk/fftjet/AbsClusteringTree.hh =================================================================== --- trunk/fftjet/AbsClusteringTree.hh (revision 43) +++ trunk/fftjet/AbsClusteringTree.hh (revision 44) @@ -1,386 +1,386 @@ +#ifndef FFTJET_ABSCLUSTERINGTREE_HH_ +#define FFTJET_ABSCLUSTERINGTREE_HH_ + //========================================================================= // AbsClusteringTree.hh // // Interface class for the clustering tree // // I. Volobouev // May 2008 //========================================================================= -#ifndef FFTJET_ABSCLUSTERINGTREE_HH_ -#define FFTJET_ABSCLUSTERINGTREE_HH_ - #include #include #include namespace fftjet { // Forward declarations class StatAccumulator; // The tree nodes are identified by their scale space level // and their sequential number within the level typedef std::pair TreeNodeId; template class AbsClusteringTree { public: typedef Cluster cluster_type; typedef LevelInfo info_type; typedef TreeNodeId NodeId; AbsClusteringTree(); AbsClusteringTree(const AbsClusteringTree&); virtual ~AbsClusteringTree() {} // Assignment operator AbsClusteringTree& operator=(const AbsClusteringTree&); // This function inserts a new level and returns its number. // The concrete implementations of this class should normally // arrange levels in such a way that the levels with the // lagest scales have smallest numbers. virtual unsigned insert(double scale, const std::vector& clusters, const LevelInfo& levelInfo) = 0; // Change the level info virtual void setLevelInfo(unsigned level, const LevelInfo& info) = 0; // Remove the given level. An attempt to remove level 0 // should be ignored. An attempt to remove level out of // range should result in a run-time error. virtual void remove(unsigned level) = 0; // Reset the whole tree leaving only the root node virtual void clear(); // Number of scale levels in the tree. Level 0 corresponds // to the root level (it will always be there). virtual unsigned nLevels() const = 0; // Level corresponding to the given scale virtual unsigned getLevel(double scale) const = 0; // Check whether a certain scale value was used in building the tree virtual bool isScaleUsed(double scale) const; // Return the data which was provided during the level // creation. Calling this function on level 0 should // result in a run-time error. virtual void getLevelData(unsigned level, double* scale, std::vector* clustersToFill, LevelInfo* levelInfo) const = 0; // Return the list of nodes whose associated clusters satisfy // a certain predicate. The predicate must define the operator // "bool operator()(const Cluster&) const" which should return // "true" in case the node is to be included in the output. template void getPassingNodes(unsigned level, const BooleanPredicate& pred, std::vector* nodesToFill) const; // Number of clusters on the given level virtual unsigned nClusters(unsigned level) const = 0; // Total number of clusters in the tree. This count should // include the top-level cluster on the root level. virtual unsigned nClusters() const; // Number of clusters on the given level satisfying // a certain predicate. The predicate must define the operator // "bool operator()(const Cluster&) const" which should return // "true" in case the cluster is to be included in the count. template unsigned clusterCount(unsigned level, const BooleanPredicate&) const; // Same thing for all levels. The answers are placed in a vector. // By convention, occupancy of level 0 will be set to 1. template void occupancyInScaleSpace(const BooleanPredicate& pred, std::vector* occupancy) const; // Statistics for clusters on the given level satisfying // a certain predicate. The predicate must define the operator // "bool operator()(const Cluster&) const". // The "PropertySelector" functor must define the operator // "double operator()(const Cluster&) const". template void clusterStats(unsigned level, const PropertySelector&, const BooleanPredicate&, StatAccumulator*) const; // The following function finds minimum and maximum tree levels // for which the number of clusters satisfying a certain predicate // equals the desired count. If no levels have the requested number // of clusters, zeros are returned for both min and max levels. // This function returns "true" if all the levels between // minimum and maximum have the same number of clusters // (equal to requested), otherwise the function returns "false". template bool clusterCountLevels(unsigned desiredCount, const BooleanPredicate&, unsigned* minLevel, unsigned* maxLevel) const; // The following function returns the number of clusters N for which // stability of the configuration is the highest. Here, "stability" // is defined as pow(N,alpha)*log(scaleMax/scaleMin), where "scaleMax" // is the maximum scale with exactly N clusters satisfying the // predicate, "scaleMin" is the minimum scale with N such clusters, // and all tree levels with scales between scaleMax and scaleMin // have the same number of clusters, N. // // The input variables "minStartingLevel" and "maxStartingLevel" // can be used to limit the searched range of scales. The scale // which corresponds to "maxStartingLevel" is used in the search, // and by default the full range of tree scales is used. The levels // which correspond to the bounds of the most stable configuration // are returned in *minLevel and *maxLevel. // // If there are no stable configurations (that is, the number of // clusters is different on each level) the function returns 0 // and sets *minLevel and *maxLevel to 0. // // Note that this function is going to be useful only when the // clustering kernel (as well as the input predicate) are // "reasonable". In particular, it is highly desirable that // the number of clusters decreases monotonously as the scale // increases. // template unsigned stableClusterCount(const BooleanPredicate& pred, unsigned* minLevel, unsigned* maxLevel, double alpha = 0.0, unsigned minStartingLevel=0, unsigned maxStartingLevel=UINT_MAX) const; // Scale corresponding to the given level virtual double getScale(unsigned level) const = 0; // Minimum and maximum scales. The default implementations // of these methods assume that level 1 has the largest scale. // 0 is returned in case the tree is not populated. virtual double minScale() const; virtual double maxScale() const; // Scales for all levels except level 0 virtual void getAllScales(std::vector* scales, bool increasingOrder=false) const; // Level-wide information (e.g., unclustered energy) provided // during the level creation. Note that the reference can be // invalidated by insert/remove/clear operations. virtual const LevelInfo& getLevelInfo(unsigned level) const = 0; // Cluster data for the given tree node. The "unchecked" // function should work faster by assuming that the cluster // with the given id exists. virtual const Cluster& getCluster(const NodeId& id) const=0; virtual const Cluster& uncheckedCluster(const NodeId& id) const=0; // "Distance" between the two clusters with given ids. // The exact meaning of this distance is up to derived classes. // It is very important, however, that it is symmetric and that // it satisfies the triangle inequality. virtual double clusterDistance(const NodeId& id1, const NodeId& id2) const = 0; // Same thing, but should work faster by assuming that // the clusters with the given ids exist virtual double uncheckedClusterDistance(const NodeId& id1, const NodeId& id2) const=0; // The extent of the node descendants (think balltree). The size // of the whole tree is undefined because the toplevel cluster // does not have a meaningful location. virtual double nodeRadius(const NodeId& id) const = 0; virtual double uncheckedNodeRadius(const NodeId& id) const = 0; // "Distance" to the parent for the given cluster. // The smaller the distance, the better the match between // the daughter and its parent. The derived classes are // encouraged to override the default implementation because // distances to parents are likely to be calculated during // the tree construction anyway. virtual double distanceToParent(const NodeId& id) const; // The following function returns the id of the closest cluster // at the same tree level which satisfies the predicate, together // with the distance to that cluster. Bad node id is returned // in case there are no neighbor clusters satisfying the predicate. template NodeId closestNeighbor(const NodeId& id, const BooleanPredicate& pred, double* distance) const; // The following function returns the id of the closest cluster // at the same or lower tree level which satisfies the predicate // and which is not one of the descendants of this cluster. // Bad node id is returned in case there are no such clusters. // // The use of this function should be limited to predicates // which can be satisfied by a resonably large fraction of // clusters, otherwise the code may end up searching almost the // entire tree at and below the level of the node with given id. template NodeId closestNonDescendant(const NodeId& id, const BooleanPredicate& pred, double* distance) const; // The following function returns the node radius using // only the clusters satisfying the predicate template double conditionalNodeRadius(const NodeId& id, const BooleanPredicate& pred) const; // Tree navigation. Note that insert and remove operations // can invalidate all node ids and numbers of daughters // obtained earlier. // // The daughters must be arranged in the order of increasing // distance to parent. badId should be returned as a parent // of the root node. // virtual NodeId parentId(const NodeId& id) const = 0; virtual unsigned nDaughters(const NodeId& id) const = 0; virtual NodeId daughterId(const NodeId& id, unsigned idau) const = 0; // The following function should be implemented if the tree // itself knows how to grow optimally. It should return // the scale which would improve the level matching in the best // possible place (plug the biggest hole). The function should // return 0.0 when there is nothing more to gain from adding // intermediate levels. // // Note that the tree should be initialized with at least two // data levels (with the largest and the smallest scales) for // the default implementation of this function to work. // // The "minRatioLog" parameter prevents the code from trying // to request extremely close scale values in pathological // configurations. If log(parentLevelScale/daughterLevelScale) // is less than "minRatioLog" then the function should not // request a new scale between these existing scales. virtual double nextBestScale(double minRatioLog) const; // The default implementation of the nextBestScale() function // uses the following function to figure out how well a level // is matched to its parent level. The function should return 0 // for a good match. If you want to improve on the default tree // growing algorithm but do not want to create a fancy multi-step // optimization strategy then it should be sufficient to override // just this function and rely on the default nextBestScale(). virtual double levelMatchDistance(unsigned daughterLevel) const; // Various useful utilities virtual NodeId closestDaughter(const NodeId& id) const; bool operator==(const AbsClusteringTree& r) const; bool operator!=(const AbsClusteringTree& r) const; // The following function says whether node with id // id1 is an ancestor of the node with id id2 virtual bool isAncestor(const NodeId& id1, const NodeId& id2) const; // The following function fills a continuous sequence of // clusters (in the order of increasing level number) which are // all related to each other through the parent/closest daughter // relationship, satisfy the given predicate, and are all // within the given distance from the given initial cluster. // The function returns the position of the initial cluster // inside the sequence. // // -1 is returned (and the sequence is cleared) if the initial // cluster does not satisfy the predicate // template int heritageLine(const NodeId& initialId, const BooleanPredicate& pred, double maxDistance, std::vector* clusterSequence) const; // Number of clusters on the given level with the // number of daughters equal to or above the given // minimum number virtual unsigned nParentsWithDaus(unsigned level, unsigned nMinDau) const; // Number of clusters connected directly to the root node. // All clusters are like that at the level with the // largest scale. virtual unsigned nRootDaughters(unsigned level) const; // The following function calculates and remembers // the drift speed for all clusters. Call this function // as needed after completing the tree construction. virtual void calculateDriftSpeeds(); // The following function calculates and remembers // the magnitude change speed for all clusters. Call this // function as needed after completing the tree construction. virtual void calculateMagnitudeSpeeds(); // The following function calculates and remembers // the lifetime for all clusters. Call this function // as needed after completing the tree construction. virtual void calculateLifetimes(); // The following function calculates and remembers the cluster // distance to the nearest neighbor cluster at the same scale. // Call this function as needed after completing the tree // construction. The "noNeighborDistance" argument specifies // what number to use when the cluster has no neighbors. virtual void calculateNearestNeighbors(double noNeighborDistance=-2.0); // The following functions transfer the information from the tree // nodes to the clusters. The "updateClusterRadiusInfo" method uses // "nodeRadius" method to evaluate the radii, and the // "updateClusterSeparationInfo" method uses "closestNonDescendant" // with an "all pass" predicate. The "failDistance" argument tells // us how to set the result if the quantity can not be calculated // (for example, there is only one cluster in the whole tree). // The results are not produced for the top level and for the level // with the smallest scale (for these levels they are meaningless). virtual void updateClusterRadiusInfo(); virtual void updateClusterSeparationInfo(double failDistance=-2.0); // A convenience function which invokes the default sequence // of the tree post-processing steps. If you do not like // this sequence, just call the functions you need directly. virtual inline void postProcess() { calculateDriftSpeeds(); calculateMagnitudeSpeeds(); calculateLifetimes(); calculateNearestNeighbors(); } // Bad node id. Should be returned when the user requests // the parent id of the root node or the closest daughter // of a daughterless node. const NodeId badId; private: virtual unsigned traceLifetime(const TreeNodeId& id, unsigned topLev); template void closestDauOrSelf(const TreeNodeId& target, const TreeNodeId& id, const BooleanPredicate& pred, double* bestDistanceSoFar, TreeNodeId* bestIdSoFar) const; template void farthestDauOrSelf(const TreeNodeId& target, const TreeNodeId& id, const BooleanPredicate& pred, double* bestDistanceSoFar, TreeNodeId* bestIdSoFar) const; mutable std::vector > neigbors_; }; } #include "fftjet/AbsClusteringTree.icc" #endif // FFTJET_ABSCLUSTERINGTREE_HH_ Index: trunk/fftjet/interpolate.hh =================================================================== --- trunk/fftjet/interpolate.hh (revision 43) +++ trunk/fftjet/interpolate.hh (revision 44) @@ -1,43 +1,43 @@ +#ifndef FFTJET_INTERPOLATE_HH_ +#define FFTJET_INTERPOLATE_HH_ + //======================================================================== // interpolate.hh // // A set of simple interpolating functions // // I. Volobouev // March 2009 //======================================================================== -#ifndef FFTJET_INTERPOLATE_HH_ -#define FFTJET_INTERPOLATE_HH_ - #include namespace fftjet { inline double interpolate_quadratic( const double xmin, const double xmax, const double v0, const double v1, const double v2, const double x) { const double h = (xmax - xmin)/2.0; const double d1 = (v2 - v0)/(xmax - xmin); const double d2 = (v2 + v0 - 2.0*v1)/h/h; const double dx = x - (xmax + xmin)/2.0; return v1 + dx*(d1 + d2*dx/2.0); } inline double lin_interpolate_1d(const double x0, const double x1, const double z0, const double z1, const double x) { return z0 + (z1 - z0)*((x - x0)/(x1 - x0)); } inline double log_interpolate_1d(const double x0, const double x1, const double z0, const double z1, const double x) { return z0 + (z1 - z0)*(log(x/x0)/log(x1/x0)); } } #endif // FFTJET_INTERPOLATE_HH_ Index: trunk/fftjet/FrequencySequentialConvolver.hh =================================================================== --- trunk/fftjet/FrequencySequentialConvolver.hh (revision 43) +++ trunk/fftjet/FrequencySequentialConvolver.hh (revision 44) @@ -1,55 +1,55 @@ +#ifndef FFTJET_FREQUENCYSEQUENTIALCONVOLVER_HH_ +#define FFTJET_FREQUENCYSEQUENTIALCONVOLVER_HH_ + //========================================================================= // FrequencySequentialConvolver.hh // // Class for performing sequential convolutions by FFT using 1d kernels // whose representations are given in the frequency domain // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_FREQUENCYSEQUENTIALCONVOLVER_HH_ -#define FFTJET_FREQUENCYSEQUENTIALCONVOLVER_HH_ - #include #include "fftjet/AbsSequentialConvolver.hh" #include "fftjet/AbsFrequencyKernel1d.hh" namespace fftjet { template class FrequencySequentialConvolver : public AbsSequentialConvolver { public: FrequencySequentialConvolver( const AbsFFTEngine* etaEngine, const AbsFFTEngine* phiEngine, const AbsFrequencyKernel1d* etaKernel, const AbsFrequencyKernel1d* phiKernel, const std::vector& phiScales, unsigned minEtaBin=0, unsigned maxEtaBin=UINT_MAX, bool fixEfficiency=false); inline virtual ~FrequencySequentialConvolver() {} protected: virtual KernelData buildKernelImageEta( double scale, const AbsFFTEngine* engine, unsigned minFixBin, unsigned maxFixBin); virtual KernelData buildKernelImagePhi( double scale, const AbsFFTEngine* engine); private: KernelData buildImage( double scale, const AbsFFTEngine* engine, const AbsFrequencyKernel1d* kernel, unsigned minFixBin, unsigned maxFixBin); const AbsFrequencyKernel1d* const etaKernel; const AbsFrequencyKernel1d* const phiKernel; }; } #include "fftjet/FrequencySequentialConvolver.icc" #endif // FFTJET_FREQUENCYSEQUENTIALCONVOLVER_HH_ Index: trunk/fftjet/rectangleOverlap.hh =================================================================== --- trunk/fftjet/rectangleOverlap.hh (revision 43) +++ trunk/fftjet/rectangleOverlap.hh (revision 44) @@ -1,57 +1,57 @@ +#ifndef FFTJET_RECTANGLEOVERLAP_HH_ +#define FFTJET_RECTANGLEOVERLAP_HH_ + //========================================================================= // rectangleOverlap.hh // // Helper functions for the signal discretization process // // I. Volobouev // January 2008 //========================================================================= -#ifndef FFTJET_RECTANGLEOVERLAP_HH_ -#define FFTJET_RECTANGLEOVERLAP_HH_ - #include #include namespace fftjet { // The "intervalOverlap" function returns the length of the overlap // of two intervals. If x1_min > x1_max or x2_min > x2_max then // the arguments will be swapped internally. double intervalOverlap(double x1_min, double x1_max, double x2_min, double x2_max); // The "rectangleRectangleOverlap" function returns the area // of the overlap of two rectangles whose sides are parallel // to the coordinate axes. // // If x1_min > x1_max then arguments will be swapped internally. // Same is true for all other such pairs. // double rectangleRectangleOverlap(double x1_min, double y1_min, double x1_max, double y1_max, double x2_min, double y2_min, double x2_max, double y2_max); // The "rectanglePolygonOverlap" function returns the area // of the overlap of a rectangle whose sides are parallel // to the coordinate axes with a convex polygon. Polygon // is defined by its sequence of vertices. Each vertex is // represented by std::pair where the first // element of the pair defines the x coordinate and the second // defines the y. The vertex sequence is represented by // std::vector of such pairs. // // If xrect_min > xrect_max then arguments will be swapped // internally. Same is true for y. // // Don't use this function in those cases where the // "rectangleRectangleOverlap" function can be used. // "rectangleRectangleOverlap" runs significantly faster. // double rectanglePolygonOverlap( double xrect_min, double yrect_min, double xrect_max, double yrect_max, const std::vector >& polygon); } #endif // FFTJET_RECTANGLEOVERLAP_HH_ Index: trunk/fftjet/AbsSequentialConvolver.hh =================================================================== --- trunk/fftjet/AbsSequentialConvolver.hh (revision 43) +++ trunk/fftjet/AbsSequentialConvolver.hh (revision 44) @@ -1,185 +1,185 @@ +#ifndef FFTJET_ABSSEQUENTIALCONVOLVER_HH_ +#define FFTJET_ABSSEQUENTIALCONVOLVER_HH_ + //========================================================================= // AbsSequentialConvolver.hh // // Base class for performing kernel convolutions sequentially (first, // convolve each eta bin, then each phi bin). This allows for different // kernel scales in different eta bins. // // This class manages a collection of kernel scans, images and, if // requested, the data normalization curves. The actual building of // kernel images is up to derived classes. // // Normally, all scans should be performed using the same kernel or // set of kernels, and with the same grid dimensions. The scans should // differ only by scale. // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_ABSSEQUENTIALCONVOLVER_HH_ -#define FFTJET_ABSSEQUENTIALCONVOLVER_HH_ - #include #include "fftjet/AbsConvolverBase.hh" #include "fftjet/AbsFFTEngine.hh" #include "fftjet/KernelData.hh" namespace fftjet { // There is one KernelData element for each eta bin // and one eta kernel for all phi bins template class SequentialKernelData : public std::vector > { public: inline SequentialKernelData(const KernelData& etaKernel) : etaKernel_(etaKernel) {} inline const KernelData& etaKernelData() const {return etaKernel_;} private: SequentialKernelData(); KernelData etaKernel_; }; // Scan information for multiple scales. // This class will destroy all KernelData // data buffers in its own destructor. template class AbsSequentialConvolver : public AbsConvolverBase { public: // This AbsSequentialConvolver object will not own the AbsFFTEngine // objects. It is the responsibility of the user of this class // to ensure that these objects are not destroyed before the // AbsSequentialConvolver itself is destroyed. // // Note that in both engines the nEta() function should return 1 // and nPhi() should return the number of bins. The number of // elements in the "phiScales" vector should be equal to nPhi() // of the "etaEngine". The scale for each eta bin will be obtained // by multiplying the overall scale by the corresponding bin scale. // AbsSequentialConvolver(const AbsFFTEngine* etaEngine, const AbsFFTEngine* phiEngine, const std::vector& phiScales, unsigned minEtaBin, unsigned maxEtaBin, bool fixEfficiency); virtual ~AbsSequentialConvolver(); // Some trivial inspectors inline bool fixingEfficiency() const {return fixEfficiency;} inline unsigned minFixedBin() const {return minEtaBin;} inline unsigned maxFixedBin() const {return maxEtaBin;} inline unsigned nEtaFFT() const {return etaEngine->nPhi();} inline unsigned nPhiFFT() const {return phiEngine->nPhi();} // In the next two methods, the "nEta" and "nPhi" arguments // which define array dimensionalities must correspond to the // grid dimensions of the FFT engine used to construct this object. // "setEventData" must be called at least once before any call // to "convolveWithKernel". void setEventData(const Real* data, unsigned nEta, unsigned nPhi); void convolveWithKernel(double scale, Real* result, unsigned nEta, unsigned nPhi); // There is no accessor method to get the event data back. // This is because there is no particular reason for this // class to keep the original event data. Only the Fourier // image of the data is stored internally. // The next method simply requests construction of kernel data // (scan, image, and, possibly, normalization) for the given scale. // After this, "isScaleProcessed" method will return "true" for // the given scale, and scan/image/normalization accessors will // return valid pointers. It is not necessary to call this function // explicitly if all you want is to convolve the data with // the kernel at the given scale. void processScale(double scale); // Find out how many scales were processed so far unsigned nProcessedScales() const; // Return the processed scales in the reverse sorted order void getProcessedScales(std::vector* scales) const; // Check whether a certain scale is already processed bool isScaleProcessed(double scale) const; // The following methods allow the user to examine the kernel // scans/images/normalization curves for the given scale. // These methods return NULL if the scale was never processed. // // Note that the pointer to SequentialKernelData can be // invalidated whenever a new scale is added to the system, // but pointers to the scan and image curves will remain valid. // const SequentialKernelData* getKernelData( double scale) const; const Real* getKernelScanPhi(double scale, unsigned etaBin) const; const Complex* getKernelImagePhi(double scale, unsigned etaBin) const; double getScanAreaPhi(double scale, unsigned etaBin) const; const Real* getKernelScanEta(double scale) const; const Complex* getKernelImageEta(double scale) const; const Real* getNormalization(double scale) const; double getScanAreaEta(double scale) const; protected: // The following helper function can be used to build // the normalization curve from the scanned kernel Real* buildEffNorm(const Real* scan, unsigned scanLength) const; // The following functions must be implemented by the derived classes. // It must build the kernel scan, image and, if requested, the inverse // efficiency curve for data normalization (the latter must be done // in case minFixBin < maxFixBin). virtual KernelData buildKernelImageEta( double scale, const AbsFFTEngine* engine, unsigned minFixBin, unsigned maxFixBin) = 0; virtual KernelData buildKernelImagePhi( double scale, const AbsFFTEngine* engine) = 0; // Complex buffer which can be used by the derived classes. // The buffer will be sufficiently large to accomodate the image // of either phi or eta kernel. Complex* complexBuffer; private: AbsSequentialConvolver(); AbsSequentialConvolver(const AbsSequentialConvolver&); AbsSequentialConvolver& operator=(const AbsSequentialConvolver&); // Collection of kernel images and normalization curves std::map > kernelProperties; typedef typename std::map >::iterator MapIterator; typedef typename std::map >::const_iterator ConstMapIterator; // Unique kernel images for the phi convolutions std::map > uniquePhiKernels; // Collection of buffers for the phi images in each eta bin std::vector dataImages; // Work buffer for eta convolutions Real* workbuf; // Input arguments const AbsFFTEngine* const etaEngine; const AbsFFTEngine* const phiEngine; const std::vector phiScales; const unsigned minEtaBin; const unsigned maxEtaBin; const bool fixEfficiency; // See if we have data bool haveData; }; } #include "fftjet/AbsSequentialConvolver.icc" #endif // FFTJET_ABSSEQUENTIALCONVOLVER_HH_ Index: trunk/fftjet/FasterEtSumRecombinationAlg.hh =================================================================== --- trunk/fftjet/FasterEtSumRecombinationAlg.hh (revision 43) +++ trunk/fftjet/FasterEtSumRecombinationAlg.hh (revision 44) @@ -1,67 +1,67 @@ +#ifndef FFTJET_FASTERETSUMRECOMBINATIONALG_HH_ +#define FFTJET_FASTERETSUMRECOMBINATIONALG_HH_ + //========================================================================= // FasterEtSumRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // This particular class assumes that the data grid is populated with // transverse energy values. // // 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 Et, Real eta, Real phi) const; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the precluster. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_FASTERETSUMRECOMBINATIONALG_HH_ -#define FFTJET_FASTERETSUMRECOMBINATIONALG_HH_ - #include "fftjet/FasterKernelRecombinationAlg.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class FasterEtSumRecombinationAlg : public FasterKernelRecombinationAlg { public: inline FasterEtSumRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) : FasterKernelRecombinationAlg( kernel, bgWeight, unlikelyBgWeight, dataCutoff, winnerTakesAll, buildCorrelationMatrix, buildClusterMask, etaBinMin, etaBinMax) {} inline virtual ~FasterEtSumRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} inline virtual bool useEtCentroid() {return false;} }; } #endif // FFTJET_FASTERETSUMRECOMBINATIONALG_HH_ Index: trunk/fftjet/LogProfileKernel.hh =================================================================== --- trunk/fftjet/LogProfileKernel.hh (revision 43) +++ trunk/fftjet/LogProfileKernel.hh (revision 44) @@ -1,60 +1,61 @@ +#ifndef FFTJET_LOGPROFILEKERNEL_HH_ +#define FFTJET_LOGPROFILEKERNEL_HH_ + //========================================================================= // LogProfileKernel.hh // // Interpolation kernel using data-defined radial profile // curve in the log(density) space. In the log profile, // the points are assumed to be _in the center_ of uniformly // spaced bins between 0 and rmax. Beyond rmax the kernel // is decaying exponentially, with the decay time determined // from the last two points in the profile. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_LOGPROFILEKERNEL_HH_ -#define FFTJET_LOGPROFILEKERNEL_HH_ - #include + #include "fftjet/AbsSymmetricKernel.hh" namespace fftjet { class LogProfileKernel : public AbsSymmetricKernel { public: // The "profile" vector must have at least two points, // and the point before last must be larger than the // last (to ensure exponential decay at large distances) LogProfileKernel(double sx, double sy, int scalePower, const std::vector& profile, double rmax=1.0); inline virtual ~LogProfileKernel() {delete [] cdf_;} inline bool isDensity() const {return true;} double axialWeight(double r, unsigned nsteps=0) const; inline static int nParameters() {return -1;} private: LogProfileKernel(const LogProfileKernel&); LogProfileKernel& operator=(const LogProfileKernel&); virtual double eval(double rsquared) const; virtual double randomRadius(double rnd) const; inline double supportDistance() const {return support_;} double annulusIntegral(double middleR, double width) const; double estimateSupport() const; void normalize(); void buildCdf(); const std::vector params_; const double maxRadius_; double normfactor_; double weightBeforeMR_; double support_; double* cdf_; }; } #endif // FFTJET_LOGPROFILEKERNEL_HH_ Index: trunk/fftjet/RealFrequencyKernel.hh =================================================================== --- trunk/fftjet/RealFrequencyKernel.hh (revision 43) +++ trunk/fftjet/RealFrequencyKernel.hh (revision 44) @@ -1,50 +1,50 @@ +#ifndef FFTJET_REALFREQUENCYKERNEL_HH_ +#define FFTJET_REALFREQUENCYKERNEL_HH_ + //========================================================================= // RealFrequencyKernel.hh // // 2d functions in the Fourier transform space whose transform // is pure real. Implemented via real kernels. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_REALFREQUENCYKERNEL_HH_ -#define FFTJET_REALFREQUENCYKERNEL_HH_ - #include #include "fftjet/AbsFrequencyKernel.hh" #include "fftjet/AbsKernel2d.hh" namespace fftjet { class RealFrequencyKernel : public AbsFrequencyKernel { public: inline RealFrequencyKernel(const AbsKernel2d* realKernel, const bool assumePointerOwnership=false) : kern_(realKernel), owns_(assumePointerOwnership) {assert(kern_);} inline virtual ~RealFrequencyKernel() { if (owns_) delete const_cast(kern_); } inline bool ownsPointer() const {return owns_;} inline std::complex operator()( const int ix, const int iy, const double scale) const { // It is natural to invert the scale in the frequency space return std::complex((*kern_)(ix, iy, 1.0/scale), 0.0); } private: RealFrequencyKernel(); RealFrequencyKernel(const RealFrequencyKernel&); RealFrequencyKernel& operator=(const RealFrequencyKernel&); const AbsKernel2d* const kern_; const bool owns_; }; } #endif // FFTJET_REALFREQUENCYKERNEL_HH_ Index: trunk/fftjet/AbsOpenDXTreeFormatter.hh =================================================================== --- trunk/fftjet/AbsOpenDXTreeFormatter.hh (revision 43) +++ trunk/fftjet/AbsOpenDXTreeFormatter.hh (revision 44) @@ -1,124 +1,124 @@ +#ifndef FFTJET_ABSOPENDXTREEFORMATTER_HH_ +#define FFTJET_ABSOPENDXTREEFORMATTER_HH_ + //========================================================================= // AbsOpenDXTreeFormatter.hh // // This class provides a way to write out AbsClusteringTree and // SparseClusteringTree information in a form suitable for visualization // by OpenDX ( www.opendx.org ). // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ABSOPENDXTREEFORMATTER_HH_ -#define FFTJET_ABSOPENDXTREEFORMATTER_HH_ - #include "fftjet/AbsTreeFormatter.hh" #include "fftjet/DXGlyphInfo.hh" namespace fftjet { // Forward declaration template < typename Cluster, typename LevelInfo, template class AbsTree > class AbsOpenDXTreeFormatter; // We need to specialize AbsOpenDXTreeFormatter class slightly // so that it works a little bit differently for AbsClusteringTree // and for SparseClusteringTree. Unfortunately, C++ does not allow // us to partially specialize member (or non-member) functions. // This is why the following work-around class is needed. namespace Private { template < typename Cluster, typename LevelInfo, template class AbsTree > struct AbsOpenDXTreeFormatterWriter { void operator()( const AbsOpenDXTreeFormatter&, const AbsTree& tree, unsigned long runNum, unsigned long evNum, std::ostream& os); }; } // The following class performs all the work necessary // to construct the OpenDX "field" object except building // the glyphs themselves. The task of building glyphs // is delegated to derived classes. template < typename Cluster, typename LevelInfo, template class AbsTree > class AbsOpenDXTreeFormatter : public AbsTreeFormatter { public: // // The "etaRange" parameter can be useful to specify whether // you would like the trees to be displayed in the eta // range at least from -etaRange to etaRange. If it is // specified as 0.0 (default) then the range of the tree // itself will be used as the range of the plot. This may be // inconvenient if you plan on looking at many plots. // // The "scaleRange" parameter can be used to change the range // of scales produced by the code. If this parameter // is left at its default value, 0.0, then the code will // simply keep the calculated scale values unchanged. // If this parameter is not 0.0 then these values will be // remapped linearly to the range from 0.0 to "scaleRange". // This may be useful to get a better picture in OpenDX // (however, the plotted scales will be somewhat arbitrary). // If you want to do this remapping, a good value of // "scaleRange" is around 6 or so -- then the scale range // in the visualization will be similar to the phi range. // inline AbsOpenDXTreeFormatter( const AbsTree& tree, const unsigned long runNum, const unsigned long evNum, const double etaRange = 0.0, const double scaleRange = 0.0) : AbsTreeFormatter(tree, runNum, evNum), etaRange_(fabs(etaRange)), scaleRange_(scaleRange) { } inline explicit AbsOpenDXTreeFormatter( const double etaRange = 0.0, const double scaleRange = 0.0) : AbsTreeFormatter(), etaRange_(fabs(etaRange)), scaleRange_(scaleRange) {} inline virtual ~AbsOpenDXTreeFormatter() {} inline double getEtaRange() const {return etaRange_;} inline double getScaleRange() const {return scaleRange_;} private: friend struct Private::AbsOpenDXTreeFormatterWriter< Cluster,LevelInfo,AbsTree>; AbsOpenDXTreeFormatter(); const double etaRange_; const double scaleRange_; void write(const AbsTree& tree, unsigned long runNum, unsigned long evNum, std::ostream& os) const; virtual DXGlyphInfo buildGlyph(const Cluster& clus,int parent) const=0; }; } #include "fftjet/AbsOpenDXTreeFormatter.icc" #endif // FFTJET_ABSOPENDXTREEFORMATTER_HH_ Index: trunk/fftjet/AbsPatternRecognitionAlg.hh =================================================================== --- trunk/fftjet/AbsPatternRecognitionAlg.hh (revision 43) +++ trunk/fftjet/AbsPatternRecognitionAlg.hh (revision 44) @@ -1,48 +1,48 @@ +#ifndef FFTJET_ABSPATTERNRECOGNITIONALG_HH_ +#define FFTJET_ABSPATTERNRECOGNITIONALG_HH_ + //========================================================================= // AbsPatternRecognitionAlg.hh // // Interface class for multiscale pattern recognition algorithms. // // The utility of this class is not expected to be great -- it is likely // that different pattern recognition strategies will need different // concepts to describe their success or failure (the "PatternInfo" type // on which this class is templated). In such a case an abstract base // class no longer serves its main purpose very well. Anyway, it may be // still useful to derive concrete pattern recognition strategies from // this class. This standardizes data flow through the sequence of // algorithm steps and makes the code somewhat easier to read and // understand (and, of course, it makes sure that different algorithms // using the same PatternInfo will be interchangeable). // // I. Volobouev // May 2009 //========================================================================= -#ifndef FFTJET_ABSPATTERNRECOGNITIONALG_HH_ -#define FFTJET_ABSPATTERNRECOGNITIONALG_HH_ - #include "fftjet/AbsClusteringTree.hh" #include "fftjet/Peak.hh" namespace fftjet { template struct AbsPatternRecognitionAlg { virtual ~AbsPatternRecognitionAlg() {} // The "run" function should examine the given clustering // tree and fill out the collection of preclusters (peaks). // The "PatternInfo" class should provide the necessary // information about the "quality" (stability in the scale // space, etc) of the precluster collection found. // // The function should return a status word. 0 means everything // is fine. The meaning of other codes is up to concrete // implementations. // virtual int run(const AbsClusteringTree& tree, std::vector* peaks, PatternInfo* info) = 0; }; } #endif // FFTJET_ABSPATTERNRECOGNITIONALG_HH_ Index: trunk/fftjet/ClusteringSequencer.hh =================================================================== --- trunk/fftjet/ClusteringSequencer.hh (revision 43) +++ trunk/fftjet/ClusteringSequencer.hh (revision 44) @@ -1,129 +1,129 @@ +#ifndef FFTJET_CLUSTERINGSEQUENCER_HH_ +#define FFTJET_CLUSTERINGSEQUENCER_HH_ + //========================================================================= // ClusteringSequencer.hh // // Generic sequencer of steps for constructing the clustering tree // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_CLUSTERINGSEQUENCER_HH_ -#define FFTJET_CLUSTERINGSEQUENCER_HH_ - #include #include "fftjet/AbsConvolverBase.hh" #include "fftjet/AbsClusteringTree.hh" #include "fftjet/SimpleFunctors.hh" #include "fftjet/PeakFinder.hh" #include "fftjet/Grid2d.hh" namespace fftjet { template class ClusteringSequencer { public: // The meaning of the constructor arguments is as follows: // // convolver -- Object which calculates and manages Fourier // transforms and convolutions. // // selector -- Object which defines whether the peak // should be included in the clustering tree. // The code will try to cast the selector // dynamically to type AbsPeakSelector. // If the cast succeeds, the "setEventData" // method of the selector will be called. // // peakFinder -- Object which performs peak finding // // initialScales -- Initial set of scales for clustering // // maxAdaptiveScales -- Maximum number of scales to use for // growing the tree on top of the initial // set of scales // // minRatioLog -- This argument will be passed to the // "nextBestScale" function of the clustering // tree. Used only when "maxAdaptiveScales" // is not 0. // // This class does not takes ownership of the "convolver" and // "selector" objects. It is a responsibility of the user of this // class to make sure that the "processEvent" method is called // only when these objects are alive. // ClusteringSequencer( AbsConvolverBase* convolver, Functor1* selector, const PeakFinder& peakFinder, const std::vector& initialScales, unsigned maxAdaptiveScales=0, double minRatioLog=0.01); virtual ~ClusteringSequencer(); // The main data processing function. It returns a status word. // 0 means everything is OK. virtual int run(const Grid2d& eventData, AbsClusteringTree* outputTree); // The following function inserts the complete event into // the clustering tree (but only those grid points whose // energy, pt, etc. is above the "dataCutoff" and also within // the eta bin limits of the peak finder). Status word is // returned (0 mean OK). The "scale" argument must be positive // but smaller than any of the scales already present in the tree. // This function will be useful in case the clustering tree // will be employed in the data analysis as a balltree. // // Note that Laplacian, Hessian, lifetime, etc. properties // will not be set in any meaningful way for this tree layer // because raw cells do not necessarily correspond to any peaks. // // This function causes calculation of the radii and separations // for the clusters in the tree if the flag "updateRadii" is true. // virtual int insertCompleteEvent( double scale, const Grid2d& eventData, AbsClusteringTree* outputTree, double dataCutoff=0.0, bool updateRadii=true); // Various inspectors inline unsigned maxAdaptiveScales() const {return maxAdaptiveScales_;} inline double minRatioLog() const {return minRatioLog_;} inline unsigned nInitialScales() const {return initialScales.size();} inline const std::vector& getInitialScales() const {return initialScales;} const PeakFinder& getPeakFinder() const {return peakFinder;} protected: virtual int processScale(double scale, const Grid2d& eventData, AbsClusteringTree* outTree); void runPeakSelector(double scale); AbsConvolverBase* const convolver; Functor1* const peakSelector; PeakFinder peakFinder; std::vector initialScales; const unsigned maxAdaptiveScales_; const double minRatioLog_; const unsigned nEta; const unsigned nPhi; std::vector peaks; std::vector selected; Real* convolvedData; private: ClusteringSequencer(); ClusteringSequencer(const ClusteringSequencer&); ClusteringSequencer& operator=(const ClusteringSequencer&); }; } #include "fftjet/ClusteringSequencer.icc" #endif // FFTJET_CLUSTERINGSEQUENCER_HH_ Index: trunk/fftjet/DiscreteGauss1d.hh =================================================================== --- trunk/fftjet/DiscreteGauss1d.hh (revision 43) +++ trunk/fftjet/DiscreteGauss1d.hh (revision 44) @@ -1,45 +1,45 @@ +#ifndef FFTJET_DISCRETEGAUSS1D_HH_ +#define FFTJET_DISCRETEGAUSS1D_HH_ + //========================================================================= // DiscreteGauss1d.hh // // The Fourier transform of the Gaussian kernel, corrected for the grid // discretization effects // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_DISCRETEGAUSS1D_HH_ -#define FFTJET_DISCRETEGAUSS1D_HH_ - #include "fftjet/AbsFrequencyKernel1d.hh" namespace fftjet { class DiscreteGauss1d : public AbsFrequencyKernel1d { public: // Parameter "sPhi" is the scale factors which has // the same meaning as the corresponding "sx" parameter // for the "Gauss1d" kernel. // // Parameter "nPhi" represents the number of cells // in the discretization grid with which this kernel // will be used. // DiscreteGauss1d(double sPhi, unsigned nPhi); inline ~DiscreteGauss1d() {} inline double sPhi() const {return sPhi_;} inline unsigned nPhi() const {return nPhi_;} // The following inherited function must be overriden std::complex operator()(int ix, double scale) const; private: DiscreteGauss1d(); const double sPhi_; const unsigned nPhi_; }; } #endif // FFTJET_DISCRETEGAUSS1D_HH_ Index: trunk/fftjet/EtSumRecombinationAlg.hh =================================================================== --- trunk/fftjet/EtSumRecombinationAlg.hh (revision 43) +++ trunk/fftjet/EtSumRecombinationAlg.hh (revision 44) @@ -1,67 +1,67 @@ +#ifndef FFTJET_ETSUMRECOMBINATIONALG_HH_ +#define FFTJET_ETSUMRECOMBINATIONALG_HH_ + //========================================================================= // EtSumRecombinationAlg.hh // // Class for fuzzy/crisp recombination algorithms in which // the proximity to the peak is determined by a kernel function. // This particular class assumes that the data grid is populated with // transverse energy values. // // 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 Et, Real eta, Real phi) const; // // The "BgData" class should contain all the info necessary for // calculating the background weight. // // In this algorithm, the eta and phi positions of the jet are taken // from the eta and phi of the precluster. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_ETSUMRECOMBINATIONALG_HH_ -#define FFTJET_ETSUMRECOMBINATIONALG_HH_ - #include "fftjet/KernelRecombinationAlg.hh" namespace fftjet { template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class EtSumRecombinationAlg : public KernelRecombinationAlg { public: inline EtSumRecombinationAlg(ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) : KernelRecombinationAlg( kernel, bgWeight, unlikelyBgWeight, dataCutoff, winnerTakesAll, buildCorrelationMatrix, buildClusterMask, etaBinMin, etaBinMax) {} inline virtual ~EtSumRecombinationAlg() {} protected: inline virtual bool recombine4Vectors() {return false;} inline virtual bool useEtCentroid() {return false;} }; } #endif // FFTJET_ETSUMRECOMBINATIONALG_HH_ Index: trunk/fftjet/SmallVector.hh =================================================================== --- trunk/fftjet/SmallVector.hh (revision 43) +++ trunk/fftjet/SmallVector.hh (revision 44) @@ -1,100 +1,100 @@ +#ifndef FFTJET_SMALLVECTOR_HH_ +#define FFTJET_SMALLVECTOR_HH_ + //========================================================================= // SmallVector.hh // // A vector class which maintains its data on the stack until its size // exceeds the template parameter, after which it migrates to the heap. // This is not a full-blown implementation, just something to keep daugher // pointers in various tree nodes mainly on the stack. The contained type // must possess the default constructor and the assignment operator // (or to be a simple built-in type). // // I. Volobouev // June 2010 //========================================================================= -#ifndef FFTJET_SMALLVECTOR_HH_ -#define FFTJET_SMALLVECTOR_HH_ - #include namespace fftjet { template class SmallVector { public: // Default constructor SmallVector(); // Constructor which makes a vector with the given size SmallVector(unsigned n); // The copy constructor SmallVector(const SmallVector&); // Converting constructor. It looks more general than the copy // constructor, but the actual copy constructor has to be created // anyway -- otherwise the compiler will generate an incorrect // default copy constructor. template SmallVector(const SmallVector&); // Destructor ~SmallVector(); // Assignment operator SmallVector& operator=(const SmallVector&); // Converting assignment operator template SmallVector& operator=(const SmallVector&); // Comparison for equality template bool operator==(const SmallVector& r) const; template inline bool operator!=(const SmallVector& r) const {return !(*this == r);} // Subscripting inline T& operator[](const unsigned i) {return data_[i];} inline const T& operator[](const unsigned i) const {return data_[i];} inline T& at(const unsigned i) {assert(i < size_); return data_[i];} inline const T& at(const unsigned i) const {assert(i < size_); return data_[i];} // Other members with obvious meaning inline unsigned size() const {return size_;} inline unsigned capacity() const {return capacity_;} inline bool empty() const {return !size_;} void clear(); void reserve(unsigned nElements); void resize(unsigned nElements); void push_back(const T& elem); void insert(unsigned position, const T& elem); void erase(unsigned position); // The following method returns the vector size if the element // is not found, otherwise it returns the element's position unsigned find(const T& elem) const; private: template friend class SmallVector; T* makeBuffer(unsigned sizeNeeded); T localData_[Len]; T* data_; unsigned capacity_; unsigned size_; }; } #include "fftjet/SmallVector.icc" #endif // FFTJET_SMALLVECTOR_HH_ Index: trunk/fftjet/Kernel1dFactory.hh =================================================================== --- trunk/fftjet/Kernel1dFactory.hh (revision 43) +++ trunk/fftjet/Kernel1dFactory.hh (revision 44) @@ -1,84 +1,84 @@ +#ifndef FFTJET_KERNEL1DFACTORY_HH_ +#define FFTJET_KERNEL1DFACTORY_HH_ + //========================================================================= // Kernel1dFactory.hh // // Factories for 1d kernel functions // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_KERNEL1DFACTORY_HH_ -#define FFTJET_KERNEL1DFACTORY_HH_ - #include #include #include #include #include "fftjet/AbsKernel1d.hh" namespace fftjet { // Abstract factory for Kernel1d types class AbsKernel1dFactory { public: virtual ~AbsKernel1dFactory() {} virtual AbsKernel1d* create(double sx, int scalePower, const std::vector&) const = 0; virtual int nParameters() const = 0; }; // Concrete factory for Kernel1d types template class Kernel1dFactory : public AbsKernel1dFactory { public: virtual ~Kernel1dFactory() {} inline int nParameters() const {return T::nParameters();} T* create(double sx, int scalePower, const std::vector& params) const { if (nParameters() >= 0) assert(params.size() == static_cast(nParameters())); return new T(sx, scalePower, params); } }; // Default kernel factory. It is intended for use inside some // interpretive language or testing framework and can iterate // over kernel classes. Usage is like this: // // DefaultKernel1dFactory factory; // AbsKernel1d *kernel = factory[userName]->create(sx, p, userParameters); // // where "userName" is one of the known kernel names, sx is // double, p is an integer which determines the dependence of the // kernel bandwidth values on the scale, and "userParameters" is // a vector of doubles. At some point after this one has to call // "delete" operator on the kernel pointer. Known kernel names // correspond to the names of classes derived from "AbsKernel1d". // For example, to create the Gaussian kernel, one can use // // AbsKernel1d *kernel = factory["Gauss1d"]->create(1, 1, // std::vector()); // // You can add your own kernel to this factory like this: // // factory["YourClass"] = new Kernel1dFactory(); // class DefaultKernel1dFactory : public std::map { public: DefaultKernel1dFactory(); virtual ~DefaultKernel1dFactory(); private: DefaultKernel1dFactory(const DefaultKernel1dFactory&); DefaultKernel1dFactory& operator=(const DefaultKernel1dFactory&); }; } #endif // FFTJET_KERNEL1DFACTORY_HH_ Index: trunk/fftjet/ConstScaleReconstruction.hh =================================================================== --- trunk/fftjet/ConstScaleReconstruction.hh (revision 43) +++ trunk/fftjet/ConstScaleReconstruction.hh (revision 44) @@ -1,77 +1,77 @@ +#ifndef FFTJET_CONSTSCALERECONSTRUCTION_HH_ +#define FFTJET_CONSTSCALERECONSTRUCTION_HH_ + //========================================================================= // ConstScaleReconstruction.hh // // Driver class for running both pattern recognition and jet recombination // at the same constant scale throughout the whole data grid // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_CONSTSCALERECONSTRUCTION_HH_ -#define FFTJET_CONSTSCALERECONSTRUCTION_HH_ - #include #include "fftjet/AbsConvolverBase.hh" #include "fftjet/SimpleFunctors.hh" #include "fftjet/PeakFinder.hh" #include "fftjet/AbsRecombinationAlg.hh" namespace fftjet { template class ConstScaleReconstruction { public: // This object will not own "convolver", "selector", and // "recombiner" pointers. It is a responsibility of the user // of this class to make sure those objects are still // alive when the "run" method is called. // // The code will attempt to cast the selector dynamically // to type AbsPeakSelector. If the cast succeeds, // the "setEventData" method of the selector will be called. // // It is possible to use other powers of the energy-like // variable besides 1 for pattern recognition (argument "etPower"). // See the comment to the "power" function inside the "Grid2d.hh" // header file. // ConstScaleReconstruction( AbsConvolverBase* convolver, Functor1* selector, const PeakFinder& peakFinder, AbsRecombinationAlg* recombiner, double etPower = 1.0); virtual ~ConstScaleReconstruction(); // The return value from the following function is the status // produced by the AbsRecombinationAlg "run" function virtual int run(double scale, const Grid2d& eventData, const BgData* bgData, unsigned nBgEta, unsigned nBgPhi, std::vector >* outputJets, VectorLike* unclustered, double* unused); private: ConstScaleReconstruction(); ConstScaleReconstruction(const ConstScaleReconstruction&); ConstScaleReconstruction& operator=(const ConstScaleReconstruction&); AbsConvolverBase* const convolver_; Functor1* const selector_; PeakFinder peakFinder_; AbsRecombinationAlg* const recombiner_; const unsigned nEta_; const unsigned nPhi_; Real* const convolvedData_; const double etPower_; std::vector peaks; Grid2d* copyGrid; }; } #include "fftjet/ConstScaleReconstruction.icc" #endif // FFTJET_CONSTSCALERECONSTRUCTION_HH_ Index: trunk/fftjet/extractSubTree.hh =================================================================== --- trunk/fftjet/extractSubTree.hh (revision 43) +++ trunk/fftjet/extractSubTree.hh (revision 44) @@ -1,55 +1,55 @@ +#ifndef FFTJET_EXTRACTSUBTREE_HH_ +#define FFTJET_EXTRACTSUBTREE_HH_ + //========================================================================= // extractSubTree.hh // // Extract all daughters of some cluster from SparseClusteringTree // into another SparseClusteringTree, as if that cluster was the // top node of the filled tree. // // I. Volobouev // May 2016 //========================================================================= -#ifndef FFTJET_EXTRACTSUBTREE_HH_ -#define FFTJET_EXTRACTSUBTREE_HH_ - #include namespace fftjet { namespace Private { template void fillSubTree(const Tree& fromTree, const typename Tree::Node& parent, Tree* toTree, const typename Tree::NodeId toParentNode) { typedef typename Tree::Node Node; typedef typename Tree::NodeId NodeId; const unsigned ndaus = parent.nDaus(); for (unsigned idau = 0; idau < ndaus; ++idau) { const Node& dau = fromTree.getNode(parent.daus()[idau]); const NodeId added = toTree->addNode( Node(dau.getCluster(), dau.originalLevel(), dau.mask()), toParentNode); fillSubTree(fromTree, dau, toTree, added); } } } template void extractSubTree(const Tree& fromTree, const typename Tree::NodeId fromParentNode, Tree* toTree) { assert(toTree); toTree->clear(); const unsigned nLevels = fromTree.nLevels(); toTree->reserveScales(nLevels); for (unsigned i=0; iaddScale(fromTree.getScale(i)); Private::fillSubTree(fromTree, fromTree.getNode(fromParentNode), toTree, 0); } } #endif // FFTJET_EXTRACTSUBTREE_HH_ Index: trunk/fftjet/ScaleSpaceKernel.hh =================================================================== --- trunk/fftjet/ScaleSpaceKernel.hh (revision 43) +++ trunk/fftjet/ScaleSpaceKernel.hh (revision 44) @@ -1,35 +1,35 @@ +#ifndef FFTJET_SCALESPACEKERNEL_HH_ +#define FFTJET_SCALESPACEKERNEL_HH_ + //========================================================================= // ScaleSpaceKernel.hh // // Generic base class for kernel functions in scale space. // Application code is not supposed to derive directly from this class. // Instead, use "AbsKernel2d" or "AbsMembershipFunction" as bases. // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_SCALESPACEKERNEL_HH_ -#define FFTJET_SCALESPACEKERNEL_HH_ - namespace fftjet { class AbsKernel2d; struct AbsMembershipFunction; class ScaleSpaceKernel { public: virtual ~ScaleSpaceKernel() {} // The following member sets eta to phi (or x to y) scale ratio virtual void setScaleRatio(double r) = 0; private: friend class AbsKernel2d; friend struct AbsMembershipFunction; inline ScaleSpaceKernel() {} }; } #endif // FFTJET_SCALESPACEKERNEL_KERNEL Index: trunk/fftjet/AbsPeakSelector.hh =================================================================== --- trunk/fftjet/AbsPeakSelector.hh (revision 43) +++ trunk/fftjet/AbsPeakSelector.hh (revision 44) @@ -1,44 +1,44 @@ +#ifndef FFTJET_ABSPEAKSELECTOR_HH_ +#define FFTJET_ABSPEAKSELECTOR_HH_ + //========================================================================= // AbsPeakSelector.hh // // Interface for peak selectors in ClusteringSequencer and // ConstScaleReconstruction. Derive you peak selector from this class // if your selector needs to do some calculations using the discretized // event data before the selection is made. Note that predicates used // by clustering trees do not get the event data -- they should normally // be derived from Functor1. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSPEAKSELECTOR_HH_ -#define FFTJET_ABSPEAKSELECTOR_HH_ - #include "fftjet/Peak.hh" #include "fftjet/Grid2d.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template struct AbsPeakSelector : public Functor1 { virtual ~AbsPeakSelector() {} // The following function will be called once per event, // before the peak selection must be made virtual void setEventData(const Grid2d& eventData) = 0; // The following operator should return "true" if the peak // is to be included in the subsequent processing virtual bool operator()(const Peak& peak) const = 0; }; // A simple selector which selects everything struct AllPeaksPass : public Functor1 { inline bool operator()(const Peak&) const {return true;} }; } #endif // FFTJET_ABSPEAKSELECTOR_HH_ Index: trunk/fftjet/AbsPhiKernel.hh =================================================================== --- trunk/fftjet/AbsPhiKernel.hh (revision 43) +++ trunk/fftjet/AbsPhiKernel.hh (revision 44) @@ -1,78 +1,78 @@ +#ifndef FFTJET_ABSPHIKERNEL_HH_ +#define FFTJET_ABSPHIKERNEL_HH_ + //========================================================================= // AbsPhiKernel.hh // // Interface class for kernel functions which are products of delta // function in eta and some reasonable distribution in phi // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSPHIKERNEL_HH_ -#define FFTJET_ABSPHIKERNEL_HH_ - #include #include "fftjet/AbsKernel2d.hh" namespace fftjet { class AbsPhiKernel : public AbsKernel2d { public: virtual ~AbsPhiKernel() {} // Attempt to set the scale ratio is ignored inline void setScaleRatio(double) {} // The function value inline double operator()(const double x, double, double) const { if (x == 0.0) assert(!"Attempt to evaluate delta function at 0"); return 0.0; } // The support rectangle inline void supportRectangle(const double scale, KernelSupportRectangle *r) const { phiSupport(scale, &r->ymin, &r->ymax); const double eps = 1.0/sqrt(DBL_MAX); r->xmin = -eps; r->xmax = eps; } virtual inline double rectangleAverage( const double x, const double y, const double scale, const double dx, const double dy) const { const double halfx(fabs(dx/2.0)); if (x - halfx <= 0.0 && x + halfx > 0.0) { // 0.28867513459481288225 is 1/(2 sqrt(3)) const double yshift = 0.28867513459481288225*dy; return (phiFcn(y - yshift, scale) + phiFcn(y + yshift, scale))/2.0/dx; } else return 0.0; } // The randomizer inline void random(double, const double r2, const double scale, double* px, double* py) const { *px = 0.0; *py = phiRandom(r2, scale); } private: virtual double phiFcn(double phi, double scale) const = 0; virtual void phiSupport(double scale, double *phimin, double *phimax) const = 0; virtual double phiRandom(double rnd, double scale) const = 0; }; } #endif // FFTJET_ABSPHIKERNEL_HH_ Index: trunk/fftjet/SequentialConvolver.hh =================================================================== --- trunk/fftjet/SequentialConvolver.hh (revision 43) +++ trunk/fftjet/SequentialConvolver.hh (revision 44) @@ -1,63 +1,63 @@ +#ifndef FFTJET_SEQUENTIALCONVOLVER_HH_ +#define FFTJET_SEQUENTIALCONVOLVER_HH_ + //========================================================================= // SequentialConvolver.hh // // Class for performing kernel convolutions by FFT using separate kernels // in eta and phi // // I. Volobouev // November 2009 //========================================================================= -#ifndef FFTJET_SEQUENTIALCONVOLVER_HH_ -#define FFTJET_SEQUENTIALCONVOLVER_HH_ - #include #include "fftjet/AbsSequentialConvolver.hh" #include "fftjet/AbsKernel1d.hh" namespace fftjet { template class SequentialConvolver : public AbsSequentialConvolver { public: // This SequentialConvolver will not own the AbsFFTEngine // and AbsKernel1d objects. It is the responsibility of // the user of this class to ensure that these objects // are not destroyed before the SequentialConvolver itself. // // The convolver will try to fix the bump reconstruction // efficiency (will apply weights to the data to simulate // the effects of boundary kernel) in case "fixEfficiency" // parameter is "true". // SequentialConvolver(const AbsFFTEngine* etaEngine, const AbsFFTEngine* phiEngine, const AbsKernel1d* etaKernel, const AbsKernel1d* phiKernel, const std::vector& phiScales, unsigned minEtaBin=0, unsigned maxEtaBin=UINT_MAX, bool fixEfficiency=false); inline virtual ~SequentialConvolver() {} protected: virtual KernelData buildKernelImageEta( double scale, const AbsFFTEngine* engine, unsigned minFixBin, unsigned maxFixBin); virtual KernelData buildKernelImagePhi( double scale, const AbsFFTEngine* engine); private: KernelData buildImage( double scale, const AbsFFTEngine* engine, const AbsKernel1d* kernel, unsigned minFixBin, unsigned maxFixBin); const AbsKernel1d* const etaKernel; const AbsKernel1d* const phiKernel; }; } #include "fftjet/SequentialConvolver.icc" #endif // FFTJET_SEQUENTIALCONVOLVER_HH_ Index: trunk/fftjet/binaryIO.hh =================================================================== --- trunk/fftjet/binaryIO.hh (revision 43) +++ trunk/fftjet/binaryIO.hh (revision 44) @@ -1,54 +1,54 @@ +#ifndef FFTJET_BINARYIO_HH_ +#define FFTJET_BINARYIO_HH_ + //========================================================================= // binaryIO.hh // // Utility functions for reading/writing "plain old data" in binary form // (not platform independent yet) // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_BINARYIO_HH_ -#define FFTJET_BINARYIO_HH_ - #include #include namespace fftjet { // The following functions perform binary I/O of built-in types template inline void write_pod(std::ostream& of, const T& pod) { of.write(reinterpret_cast(&pod), sizeof(T)); } template inline void read_pod(std::istream& in, T* pod) { assert(pod); in.read(reinterpret_cast(pod), sizeof(T)); } template inline void write_pod_array(std::ostream& of, const T* pod, const unsigned len) { if (len) { assert(pod); of.write(reinterpret_cast(pod), len*sizeof(T)); } } template inline void read_pod_array(std::istream& in, T* pod, const unsigned len) { if (len) { assert(pod); in.read(reinterpret_cast(pod), len*sizeof(T)); } } } #endif // FFTJET_BINARYIO_HH_ Index: trunk/fftjet/MagneticSmearingKernel.hh =================================================================== --- trunk/fftjet/MagneticSmearingKernel.hh (revision 43) +++ trunk/fftjet/MagneticSmearingKernel.hh (revision 44) @@ -1,131 +1,131 @@ +#ifndef FFTJET_MAGNETICSMEARINGKERNEL_HH_ +#define FFTJET_MAGNETICSMEARINGKERNEL_HH_ + //========================================================================= // MagneticSmearingKernel.hh // // This function represents smearing of a jet in phi due to the // presence of magnetic field parallel to the beam axis, and it is // really a 1-d function (no smearing in eta). // // The scale for this kernel should be set to 1/(jet Pt). // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_MAGNETICSMEARINGKERNEL_HH_ -#define FFTJET_MAGNETICSMEARINGKERNEL_HH_ - #include "fftjet/AbsPhiKernel.hh" #include "fftjet/LinearInterpolator1d.hh" namespace fftjet { // We will need the jet fragmentation function in order to build // the kernel. We will assume that this function is implemented // as a functor which has "double operator()(double x) const". // We will also assume that the lifetime of the fragmentation // function is longer than the lifetime of the kernel. // template class MagneticSmearingKernel : public AbsPhiKernel { public: // Charged particle deviation in phi due to the presence // of magnetic field can be calculated as follows: // // numeratorConst // sin(phi) = -------------- // Pt // // Where "numeratorConst" depends on the units in which Pt, // magnetic field, and calorimeter radius are measured. // The "nominal" value of this constant is determined from // the relationship // // sin(phi) = R/(2*a) // a = Pt*1.0e9/(c*B), so that numeratorConst = c*B*R/2.0e9 // // where // // Pt -- transverse particle momentum, in GeV/c // R -- "typical" calorimeter radius, in m // a -- gyration radius, in m, for a particle with given Pt // and charge equal in magnitude to the charge of // the electron (that is, |Z| == 1). // c -- speed of light in m/s (c = 299792458.0) // B -- magnetic field, in Tesla // // Of course, if particle Pt is below "numeratorConst" then // it simply never reaches the calorimeter. // // The sign of phi depends on the sign of charge and direction // of the magnetic field. We will assume that the number of // negatively charged particles is the same, on average, as // the number of positively charged particles, and that the // smearing is symmetric in phi. We will further assume that // the fragmentation function for particles with |Z| == 1 // is the same as for |Z| == 2. // // The constructor arguments are as follows: // // fcn -- the fragmentation function functor // // numeratorConst -- explained above // // charge1Fraction -- average energy fraction deposited by // particles with |Z| == 1 // // charge2Fraction -- average energy fraction deposited by // particles with |Z| == 2 // // samplesPerBin -- number of points per each phi bin // in which the fragmentation function // gets evaluated (and then integrated over) // // objectOwnsFcn -- if "true", the fragmentation function functor // will be deleted in the destructor // // It is assumed that the energy fraction deposited by neutral // particles is (1 - charge1Fraction - charge2Fraction). // MagneticSmearingKernel(const FragmentationFunction* fcn, double numeratorConst, double charge1Fraction, double charge2Fraction, unsigned samplesPerBin=100, bool objectOwnsFcn=false); virtual ~MagneticSmearingKernel(); inline double numeratorConst() const {return numeratorConst_;} inline double charge1Fraction() const {return charge1Fraction_;} inline double charge2Fraction() const {return charge2Fraction_;} inline unsigned samplesPerBin() const {return samplesPerBin_;} inline bool isDensity() const {return true;} double rectangleAverage(double x, double y, double scale, double dx, double dy) const; private: MagneticSmearingKernel(); MagneticSmearingKernel(const MagneticSmearingKernel&); MagneticSmearingKernel& operator=(const MagneticSmearingKernel&); double phiFcn(double phi, double scale) const; void phiSupport(double scale, double *phimin, double *phimax) const; double phiRandom(double rnd, double scale) const; double deltaFunAverage(double x, double y, double dx, double dy) const; void buildRandomizer(); const FragmentationFunction* fcn_; const double numeratorConst_; const double charge0Fraction_; const double charge1Fraction_; const double charge2Fraction_; const unsigned samplesPerBin_; const bool objectOwnsFcn_; LinearInterpolator1d* randomizer_; }; } #include "fftjet/MagneticSmearingKernel.icc" #endif // FFTJET_MAGNETICSMEARINGKERNEL_HH_ Index: trunk/fftjet/FFTWFloatEngine.hh =================================================================== --- trunk/fftjet/FFTWFloatEngine.hh (revision 43) +++ trunk/fftjet/FFTWFloatEngine.hh (revision 44) @@ -1,29 +1,29 @@ +#ifndef FFTJET_FFTWFLOATENGINE_HH_ +#define FFTJET_FFTWFLOATENGINE_HH_ + //========================================================================= // FFTWFloatEngine.hh // // Single precision FFT engine for jet reconstruction which uses // the FFTW library from http://www.fftw.org // // I. Volobouev // January 2008 //========================================================================= -#ifndef FFTJET_FFTWFLOATENGINE_HH_ -#define FFTJET_FFTWFLOATENGINE_HH_ - -#include "fftjet/FFTWEngine.hh" #include "fftw3.h" +#include "fftjet/FFTWEngine.hh" namespace fftjet { class FFTWFloatEngine : public FFTWEngine { public: FFTWFloatEngine(unsigned nEta, unsigned nPhi, unsigned optimization = FFTW_PATIENT); virtual ~FFTWFloatEngine(); }; } #include "fftjet/FFTWFloatEngine.icc" #endif // FFTJET_FFTWFLOATENGINE_HH_ Index: trunk/fftjet/VectorRecombinationAlgFactory.hh =================================================================== --- trunk/fftjet/VectorRecombinationAlgFactory.hh (revision 43) +++ trunk/fftjet/VectorRecombinationAlgFactory.hh (revision 44) @@ -1,130 +1,130 @@ +#ifndef FFTJET_VECTORRECOMBINATIONALGFACTORY_HH_ +#define FFTJET_VECTORRECOMBINATIONALGFACTORY_HH_ + //========================================================================= // VectorRecombinationAlgFactory.hh // // Factories for algorithms used to recombine collections of 4-vectors // // Intended usage is like this: once you know which classes to use // in place of VectorLike, BgData, and VBuilder, instantiate // the "DefaultVectorRecombinationAlgFactory". You can then choose // which algorithm to use dynamically, by name, as follows: // // std::string algoName; // .... Get the algoName somehow .... // AbsVectorRecombinationAlg* alg = 0; // DefaultVectorRecombinationAlgFactory<...> aFactory; // if (aFactory[algoName] == NULL) // { // .... Process invalid algorithm name error ..... // } // else // alg = aFactory[algoName]->create(kernel, ...); // .... Perform some work with alg .... // delete alg; // // The correspondence between valid algorithm names and classes is // // algoName: Actual class name: // // "Kernel" KernelVectorRecombinationAlg // "EtCentroid" EtCentroidVectorRecombinationAlg // "EtSum" EtSumVectorRecombinationAlg // // I. Volobouev // June 2009 //========================================================================= -#ifndef FFTJET_VECTORRECOMBINATIONALGFACTORY_HH_ -#define FFTJET_VECTORRECOMBINATIONALGFACTORY_HH_ - -#include #include +#include #include "fftjet/ScaleSpaceKernel.hh" #include "fftjet/AbsVectorRecombinationAlg.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template class AbsVectorRecombinationAlgFactory { public: typedef double (VectorLike::* VectorLikeMemberFunction)() const; virtual ~AbsVectorRecombinationAlgFactory() {} virtual AbsVectorRecombinationAlg* create( ScaleSpaceKernel* kernel, VectorLikeMemberFunction etFcn, VectorLikeMemberFunction etaFcn, VectorLikeMemberFunction phiFcn, const Functor2* bgWeight, double unlikelyBgWeight, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask) const = 0; }; template < template < typename VectorLike, typename BgData, typename VBuilder > class VectorRecombinationAlg, typename VectorLike, typename BgData, typename VBuilder > class VectorRecombinationAlgFactory : public AbsVectorRecombinationAlgFactory { public: typedef double (VectorLike::* VectorLikeMemberFunction)() const; virtual ~VectorRecombinationAlgFactory() {} VectorRecombinationAlg* create( ScaleSpaceKernel* kernel, VectorLikeMemberFunction etFcn, VectorLikeMemberFunction etaFcn, VectorLikeMemberFunction phiFcn, const Functor2* bgWeight, double unlikelyBgWeight, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask) const { return new VectorRecombinationAlg( kernel, etFcn, etaFcn, phiFcn, bgWeight, unlikelyBgWeight, winnerTakesAll, buildCorrelationMatrix, buildClusterMask); } }; template < typename VectorLike, typename BgData, typename VBuilder > class DefaultVectorRecombinationAlgFactory : public std::map*> { public: DefaultVectorRecombinationAlgFactory(); virtual ~DefaultVectorRecombinationAlgFactory(); private: typedef AbsVectorRecombinationAlgFactory BaseFactory; DefaultVectorRecombinationAlgFactory( const DefaultVectorRecombinationAlgFactory&); DefaultVectorRecombinationAlgFactory& operator=( const DefaultVectorRecombinationAlgFactory&); }; } #include "fftjet/VectorRecombinationAlgFactory.icc" #endif // FFTJET_VECTORRECOMBINATIONALGFACTORY_HH_ Index: trunk/fftjet/AbsFrequencyKernel.hh =================================================================== --- trunk/fftjet/AbsFrequencyKernel.hh (revision 43) +++ trunk/fftjet/AbsFrequencyKernel.hh (revision 44) @@ -1,33 +1,33 @@ +#ifndef FFTJET_ABSFREQUENCYKERNEL_HH_ +#define FFTJET_ABSFREQUENCYKERNEL_HH_ + //========================================================================= // AbsFrequencyKernel.hh // // Interface class for 2d functions in the FFT transform space. // Can be used to represent filters. // // This interface is much simpler than the interface used for real // kernels. For the intended range of applications, we do not anticipate // the need to represent delta functions in the frequency space and to // generate random harmonics. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSFREQUENCYKERNEL_HH_ -#define FFTJET_ABSFREQUENCYKERNEL_HH_ - #include namespace fftjet { struct AbsFrequencyKernel { virtual ~AbsFrequencyKernel() {} // ix and iy are harmonics. The convention used here // is that they can be both positive and negative. virtual std::complex operator()(int ix, int iy, double scale) const = 0; }; } #endif // FFTJET_ABSFREQUENCYKERNEL_HH_ Index: trunk/fftjet/AbsKernel2d.hh =================================================================== --- trunk/fftjet/AbsKernel2d.hh (revision 43) +++ trunk/fftjet/AbsKernel2d.hh (revision 44) @@ -1,166 +1,166 @@ +#ifndef FFTJET_ABSKERNEL2D_HH_ +#define FFTJET_ABSKERNEL2D_HH_ + //========================================================================= // AbsKernel2d.hh // // Interface class for 2d kernel functions in scale space // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_ABSKERNEL2D_HH_ -#define FFTJET_ABSKERNEL2D_HH_ - #include "fftjet/ScaleSpaceKernel.hh" namespace fftjet { // Utility class for kernel-related calculations struct KernelSupportRectangle { double xmin; double xmax; double ymin; double ymax; }; // Interface class for 2d kernel functions. In all calculations, // it should be assumed that the "scale" parameter is positive. // // Subsequently, derived classes will be used with a templated // factory. In order to play nice with that factory, each kernel // class must have a constructor which looks like // // ClassName(double sx, double sy, int n, const std::vector& v); // // "sx" and "sy" are typically (but not necessarily) scale factors // for the x and y directions which map the scale argument into // bandwidth values. n is normally the power of the scale in the // bandwidth calculations: 1 means bandwidth is proportional to // the scale, 0 means the bandwidth is independent from the scale, and // -1 means that bandwidth is inversely proportional to the scale. // Vector "v" is a vector of parameters for the kernels. In addition, // each class must have a static function // // static int nParameters(); // // which defines the number of parameters the class expects to get // in "v". If this function returns a negative number, it means that // the number of parameters is arbitrary (for example, this would be // the case for kernels defined by interpolation tables). // // Of course, more complicated constructors may be needed for some // kernels. Such constructors can not be used with the common factory. // These kernels will, consequently, require special parsers in // a configuration file or in an interpretive language interface. // class AbsKernel2d : public ScaleSpaceKernel { public: virtual ~AbsKernel2d() {} // The following member sets eta to phi (or x to y) scale ratio virtual void setScaleRatio(double r) = 0; // The function value virtual double operator()(double x, double y, double scale) const = 0; // The following function should return the smallest rectangle // which bounds the region in which the kernel value is not 0. // // The size of the rectangle can strongly affect the speed // of various scanning algorithms. Therefore, resist the // temptation to return DBL_MAX for kernels with infinite support // and do something more computationally meaningful instead. // For example, the Gaussian kernel, while having infinite // support theoretically, turns to 0 in double precision // representation outside the square with corners at (-39,-39) // and (39,39). virtual void supportRectangle(double scale, KernelSupportRectangle *r) const = 0; // An average over a small rectangle centered at (x, y) // with sides dx and dy. This is more precise (but slower) // than calculating the value of the function at the rectangle // center. // // Degenerate kernels (e.g., delta functions in one or two // dimensions) must override the default implementation // of this function. virtual double rectangleAverage(double x, double y, double scale, double dx, double dy) const; // The following function tells us whether the kernel // can be treated as a probability density. Normally, // it means that the kernel must be non-negative everywhere // and that its integral over the whole space must be positive. virtual bool isDensity() const = 0; // The following function returns random numbers distributed // according to the kernel treated as a probability density. The // arguments "r1" and "r2" should be random numbers between 0 and 1. // Calling this function on a kernel which is not a density // should result in a run-time error. virtual void random(double r1, double r2, double scale, double* px, double* py) const = 0; // The following function scans the kernel for use with planar // geometries. The function returns the area under the scan before // normalization. It is assumed that the coordinate of the // bottom left point in the grid is (xmin + xstep/2, ymin + ystep/2), // where xstep = (xmax - xmin)/nx and ystep = (ymax - ymin)/ny. template double scanPlane(Real* to, unsigned nx, unsigned ny, double x0, double y0, double scale, double xmin, double xmax, double ymin, double ymax, double normalizationArea=1.0, bool normalize=true) const; // The following function scans the kernel for use with // cylindrical geometries. The function returns the // area under the scan before normalization. template double scanCylinder(Real* to, unsigned nEta, unsigned nPhi, double eta0, double phi0, double scale, double etaMin, double etaMax, double phiBin0Edge, double normalizationArea=1.0, bool normalize=true) const; // The following function assumes periodic 2d region // with the 2*Pi periods in each direction. to[0][0] // corresponds to kernel argument (0, 0). The function // returns the area under the scan before normalization. template double scanFFT(Real* to, unsigned nx, unsigned ny, double scale, double normalizationArea=1.0, bool normalize=true) const; private: // Various utility functions follow template long double fillFFTy_(Real* to, unsigned nx, unsigned ny, double x, double scale) const; void optimalLineScan(unsigned nx, double xmin, double xmax, double rmin, double rmax, unsigned* nxmin, unsigned* nxmax) const; // Clear out the whole scanned area if the following functions // return "true". In this case some optimal scan may be // available. bool optimalPlaneScan(unsigned nx, unsigned ny, double x0, double y0, double scale, double xmin, double xmax, double ymin, double ymax, unsigned* nxmin, unsigned* nxmax, unsigned* nymin, unsigned* nymax) const; bool optimalCylinderScan(unsigned nEta, unsigned nPhi, double eta0, double phi0, double scale, double etaMin, double etaMax, double phiBin0Edge, unsigned* netaMin, unsigned* netaMax, unsigned* nphiMin, unsigned* nphiMax) const; }; } #include "fftjet/AbsKernel2d.icc" #endif // FFTJET_ABSKERNEL2D_HH_ Index: trunk/fftjet/ProfileKernel.hh =================================================================== --- trunk/fftjet/ProfileKernel.hh (revision 43) +++ trunk/fftjet/ProfileKernel.hh (revision 44) @@ -1,56 +1,56 @@ +#ifndef FFTJET_PROFILEKERNEL_HH_ +#define FFTJET_PROFILEKERNEL_HH_ + //========================================================================= // ProfileKernel.hh // // Interpolation kernel using data-defined radial profile // curve. It is assumed that profile values are calculated // on the radius interval [0, 1] with uniform step in radius. // The first value corresponds to r = 0.0, and the last one // corresponds to 1.0 - 1.0/profile.size(). Each profile value // must be non-negative. // // I. Volobouev // March 2009 //========================================================================= -#ifndef FFTJET_PROFILEKERNEL_HH_ -#define FFTJET_PROFILEKERNEL_HH_ - #include #include "fftjet/AbsSymmetricKernel.hh" namespace fftjet { class ProfileKernel : public AbsSymmetricKernel { public: ProfileKernel(double sx, double sy, int scalePower, const std::vector& profile); inline virtual ~ProfileKernel() {delete [] cdf_;} inline bool isDensity() const {return true;} double axialWeight(double r, unsigned nsteps=0) const; inline static int nParameters() {return -1;} protected: void normalizeProfile(); double rawInterpolatedValue(double r) const; double gaussIntegral() const; private: ProfileKernel(const ProfileKernel&); ProfileKernel& operator=(const ProfileKernel&); virtual double eval(double rsquared) const; virtual double randomRadius(double rnd) const; inline virtual double supportDistance() const {return 1.0;} void buildCdf(); unsigned rCellNumber(double r, double *delta) const; const std::vector params_; double normfactor_; double* cdf_; }; } #endif // FFTJET_PROFILEKERNEL_HH_ Index: trunk/fftjet/RecombinationAlgFactory.hh =================================================================== --- trunk/fftjet/RecombinationAlgFactory.hh (revision 43) +++ trunk/fftjet/RecombinationAlgFactory.hh (revision 44) @@ -1,133 +1,133 @@ +#ifndef FFTJET_RECOMBINATIONALGFACTORY_HH_ +#define FFTJET_RECOMBINATIONALGFACTORY_HH_ + //========================================================================= // RecombinationAlgFactory.hh // // Factories for recombination algorithms. // // Intended usage is like this: once you know which classes to use // in place of Real, VectorLike, BgData, and VBuilder, instantiate // the "DefaultRecombinationAlgFactory". You can then choose which // algorithm to use dynamically, by name, as follows: // // std::string algoName; // .... Get the algoName somehow .... // AbsRecombinationAlg* alg = 0; // DefaultRecombinationAlgFactory<...> aFactory; // if (aFactory[algoName] == NULL) // { // .... Process invalid algorithm name error ..... // } // else // alg = aFactory[algoName]->create(kernel, ...); // .... Perform some work with alg .... // delete alg; // // The correspondence between valid algorithm names and classes is // // algoName: Actual class name: // // "Kernel" KernelRecombinationAlg // "EtCentroid" EtCentroidRecombinationAlg // "EtSum" EtSumRecombinationAlg // "FasterKernel" FasterKernelRecombinationAlg // "FasterEtCentroid" FasterEtCentroidRecombinationAlg // "FaterEtSum" FasterEtSumRecombinationAlg // // I. Volobouev // April 2009 //========================================================================= -#ifndef FFTJET_RECOMBINATIONALGFACTORY_HH_ -#define FFTJET_RECOMBINATIONALGFACTORY_HH_ - #include #include #include "fftjet/ScaleSpaceKernel.hh" #include "fftjet/AbsRecombinationAlg.hh" #include "fftjet/SimpleFunctors.hh" namespace fftjet { template class AbsRecombinationAlgFactory { public: virtual ~AbsRecombinationAlgFactory() {} virtual AbsRecombinationAlg* create( ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) const = 0; }; template < template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class RecombinationAlg, typename Real, typename VectorLike, typename BgData, typename VBuilder > class RecombinationAlgFactory : public AbsRecombinationAlgFactory { public: virtual ~RecombinationAlgFactory() {} RecombinationAlg* create( ScaleSpaceKernel* kernel, const Functor2* bgWeight, double unlikelyBgWeight, double dataCutoff, bool winnerTakesAll, bool buildCorrelationMatrix, bool buildClusterMask, unsigned etaBinMin=0, unsigned etaBinMax=UINT_MAX) const { return new RecombinationAlg( kernel, bgWeight, unlikelyBgWeight, dataCutoff, winnerTakesAll, buildCorrelationMatrix, buildClusterMask, etaBinMin, etaBinMax); } }; template < typename Real, typename VectorLike, typename BgData, typename VBuilder > class DefaultRecombinationAlgFactory : public std::map*> { public: DefaultRecombinationAlgFactory(); virtual ~DefaultRecombinationAlgFactory(); private: typedef AbsRecombinationAlgFactory BaseFactory; DefaultRecombinationAlgFactory( const DefaultRecombinationAlgFactory&); DefaultRecombinationAlgFactory& operator=( const DefaultRecombinationAlgFactory&); }; } #include "fftjet/RecombinationAlgFactory.icc" #endif // FFTJET_RECOMBINATIONALGFACTORY_HH_