|
[Rivet-svn] r2766 - in trunk: . include/Rivet include/Rivet/Math include/Rivet/Projections src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Nov 25 14:11:16 GMT 2010
Author: buckley Date: Thu Nov 25 14:11:16 2010 New Revision: 2766 Log: Many improvements and extensions relating to jets and jet shapes. Not quite finished or tested yet, but it compiles and should be basically sensible... too many code changes to keep on coding! Modified: trunk/ChangeLog trunk/include/Rivet/Analysis.hh trunk/include/Rivet/Cmp.fhh trunk/include/Rivet/Math/MathHeader.hh trunk/include/Rivet/Math/MathUtils.hh trunk/include/Rivet/Math/Vector4.hh trunk/include/Rivet/Projections/ClosestJetShape.hh trunk/include/Rivet/Projections/FastJets.hh trunk/include/Rivet/Projections/FinalState.hh trunk/include/Rivet/Projections/JetAlg.hh trunk/include/Rivet/Projections/JetShape.hh trunk/src/Analyses/CDF_2005_S6217184.cc trunk/src/Projections/ClosestJetShape.cc trunk/src/Projections/FastJets.cc trunk/src/Projections/JetShape.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/ChangeLog Thu Nov 25 14:11:16 2010 (r2766) @@ -1,3 +1,23 @@ +2010-11-25 Andy Buckley <andy at insectnation.org> + + * Cmp.hh: Adding ASC/DESC (and ANTISORTED) as preferred + non-EQUIVALENT enum value synonyms over misleading + SORTED/UNSORTED. + + * Change of rapidity scheme enum name to RapScheme + + * Reworking JetShape a bit further: constructor args now avoid + inconsistencies (it was previously possible to define incompatible + range-ends and interval). Internal binning implementation also + reworked to use a vector of bin edges: the bin details are + available via the interface. The general jet pT cuts can be + applied via the JetShape constructor. + + * MathUtils.hh: Adding linspace and logspace utility + functions. Useful for defining binnings. + + * Adding more general cuts on jet pT and (pseudo)rapidity. + 2010-11-11 Andy Buckley <andy at insectnation.org> * Adding special handling of FourMomentum::mass() for computed Modified: trunk/include/Rivet/Analysis.hh ============================================================================== --- trunk/include/Rivet/Analysis.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Analysis.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -17,7 +17,7 @@ /// @def vetoEvent /// Preprocessor define for vetoing events, including the log message and return. -#define vetoEvent { getLog() << Log::DEBUG << "Vetoing event on line " << __LINE__ << " of " << __FILE__ << endl; return; } +#define vetoEvent { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } namespace Rivet { Modified: trunk/include/Rivet/Cmp.fhh ============================================================================== --- trunk/include/Rivet/Cmp.fhh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Cmp.fhh Thu Nov 25 14:11:16 2010 (r2766) @@ -13,9 +13,12 @@ /// Enumerate the possible states of a Cmp object. enum CmpState { UNDEFINED = -2, //< Undefined state. + ASC = -1, //< The two corresponding objects are in ascending order. ORDERED = -1, //< The two corresponding objects are ordered. EQUIVALENT = 0, //< The two corresponding objects are equivalent. - UNORDERED = 1 //< The two corresponding objects are unordered. + DESC = 1, //< The two corresponding objects are in descending order. + UNORDERED = 1, //< The two corresponding objects are anti-ordered. + ANTIORDERED = 1 //< The two corresponding objects are anti-ordered. }; Modified: trunk/include/Rivet/Math/MathHeader.hh ============================================================================== --- trunk/include/Rivet/Math/MathHeader.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Math/MathHeader.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -51,8 +51,8 @@ /// Enum for signs of numbers. enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 }; - /// Enum for longitudinal variable to be used in calculating \f$ R \f$ - enum DeltaRScheme { PSEUDORAPIDITY, RAPIDITY }; + /// Enum for rapidity variable to be used in calculating \f$ R \f$, applying rapidity cuts, etc. + enum RapScheme { PSEUDORAPIDITY = 0, ETA = 0, RAPIDITY = 1, YRAP = 1 }; /// Enum for range of \f$ \phi \f$ to be mapped into enum PhiMapping { MINUSPI_PLUSPI, ZERO_2PI }; Modified: trunk/include/Rivet/Math/MathUtils.hh ============================================================================== --- trunk/include/Rivet/Math/MathUtils.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Math/MathUtils.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -3,6 +3,7 @@ #define RIVET_MathUtils_HH #include "Rivet/Math/MathHeader.hh" +#include "Rivet/RivetBoost.hh" #include <cassert> namespace Rivet { @@ -86,6 +87,14 @@ } + /// Alternative version of inRange which accepts a pair for the range arguments. + template<typename NUM> + inline bool inRange(NUM value, pair<NUM, NUM> lowhigh, + RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) { + return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound); + } + + /// Determine if @a value is in the range @a low to @a high, with boundary /// types defined by @a lowbound and @a highbound. /// @todo Optimise to one-line at compile time? @@ -103,6 +112,47 @@ } + /// Alternative version of @c inRange<int> which accepts a pair for the range arguments. + inline bool inRange(int value, pair<int, int> lowhigh, + RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) { + return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound); + } + + + /// Make a list of @a nbins + 1 values equally spaced between @a start and @a end inclusive. + inline vector<double> linspace(double start, double end, size_t nbins) { + assert(end >= start); + assert(nbins > 0); + vector<double> rtn; + const double interval = (end-start)/static_cast<double>(nbins); + double edge = start; + while (inRange(edge, start, end, CLOSED, CLOSED)) { + rtn.push_back(edge); + edge += interval; + } + assert(rtn.size() == nbins+1); + return rtn; + } + + + /// Make a list of @a nbins + 1 values exponentially spaced between @a start and @a end inclusive. + inline vector<double> logspace(double start, double end, size_t nbins) { + assert(end >= start); + assert(start > 0); + assert(nbins > 0); + const double logstart = std::log(start); + const double logend = std::log(end); + const vector<double> logvals = linspace(logstart, logend, nbins); + vector<double> rtn; + foreach (double logval, logvals) { + rtn.push_back(std::exp(logval)); + } + assert(rtn.size() == nbins+1); + return rtn; + } + + + /// Named number-type squaring operation. template <typename Num> inline Num sqr(Num a) { @@ -208,9 +258,9 @@ /// Calculate the distance between two points in 2D rapidity-azimuthal /// ("eta-phi") space. The phi values are given in radians. - inline double deltaR(double y1, double phi1, double y2, double phi2) { + inline double deltaR(double rap1, double phi1, double rap2, double phi2) { const double dphi = deltaPhi(phi1, phi2); - return sqrt( sqr(y1-y2) + sqr(dphi) ); + return sqrt( sqr(rap1-rap2) + sqr(dphi) ); } /// Calculate a rapidity value from the supplied energy @a E and longitudinal momentum @pz. Modified: trunk/include/Rivet/Math/Vector4.hh ============================================================================== --- trunk/include/Rivet/Math/Vector4.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Math/Vector4.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -495,7 +495,7 @@ /// case since @c RAPIDITY is only a valid option for vectors whose type is /// really the FourMomentum derived class. inline double deltaR(const FourVector& a, const FourVector& b, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR(a.vector3(), b.vector3()); @@ -517,7 +517,7 @@ inline double deltaR(const FourVector& v, double eta2, double phi2, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR(v.vector3(), eta2, phi2); @@ -538,7 +538,7 @@ inline double deltaR(double eta1, double phi1, const FourVector& v, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY : return deltaR(eta1, phi1, v.vector3()); @@ -563,7 +563,7 @@ /// rapidity (a relativistic energy-momentum quantity) is to be used: this can /// be chosen via the optional scheme parameter. inline double deltaR(const FourMomentum& a, const FourMomentum& b, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR(a.vector3(), b.vector3()); @@ -576,7 +576,7 @@ inline double deltaR(const FourMomentum& v, double eta2, double phi2, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR(v.vector3(), eta2, phi2); @@ -590,7 +590,7 @@ inline double deltaR(double eta1, double phi1, const FourMomentum& v, - DeltaRScheme scheme = PSEUDORAPIDITY) { + RapScheme scheme = PSEUDORAPIDITY) { switch (scheme) { case PSEUDORAPIDITY: return deltaR(eta1, phi1, v.vector3()); Modified: trunk/include/Rivet/Projections/ClosestJetShape.hh ============================================================================== --- trunk/include/Rivet/Projections/ClosestJetShape.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Projections/ClosestJetShape.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -58,7 +58,7 @@ /// Constructor. ClosestJetShape(const FinalState& fs, const vector<FourMomentum>& jetaxes, double rmin=0.0, double rmax=0.7, double interval=0.1, - double r1minPsi=0.3, DeltaRScheme distscheme=RAPIDITY); + double r1minPsi=0.3, RapScheme distscheme=RAPIDITY); /// Clone on the heap. virtual const Projection* clone() const { @@ -156,7 +156,7 @@ double _r1minPsi; /// Rapidity scheme - DeltaRScheme _distscheme; + RapScheme _distscheme; /// Number of bins size_t _nbins; Modified: trunk/include/Rivet/Projections/FastJets.hh ============================================================================== --- trunk/include/Rivet/Projections/FastJets.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Projections/FastJets.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -96,16 +96,7 @@ } /// Get the jets (unordered). - Jets jets(double ptmin = 0.0) const; - - /// Get the jets, ordered by \f$ p_T \f$. - Jets jetsByPt(double ptmin = 0.0) const; - - /// Get the jets, ordered by \f$ E \f$. - Jets jetsByE(double ptmin = 0.0) const; - - /// Get the jets, ordered by rapidity. - Jets jetsByRapidity(double ptmin = 0.0) const; + Jets _jets(double ptmin = 0.0) const; /// Get the pseudo jets (unordered). PseudoJets pseudoJets(double ptmin = 0.0) const; Modified: trunk/include/Rivet/Projections/FinalState.hh ============================================================================== --- trunk/include/Rivet/Projections/FinalState.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Projections/FinalState.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -70,7 +70,6 @@ return particles(cmpParticleByAscAbsPseudorapidity); } - /// Access the projected final-state particles. virtual size_t size() const { return _theParticles.size(); } @@ -79,6 +78,9 @@ /// @deprecated Is this final state empty? virtual bool isEmpty() const { return _theParticles.empty(); } + /// Minimum-\f$ p_\perp \f$ requirement. + virtual size_t ptMin() const { return _ptmin; } + public: Modified: trunk/include/Rivet/Projections/JetAlg.hh ============================================================================== --- trunk/include/Rivet/Projections/JetAlg.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Projections/JetAlg.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -78,13 +78,29 @@ _useInvisibles = useinvis; } - /// Get jets in no guaranteed order, with an optional cut on min \f$ p_\perp \f$. - /// @todo Provide an extra optional cut on ptmax? - virtual Jets jets(double ptmin=0.0) const = 0; + /// Get jets in no guaranteed order, with optional cuts on \f$ p_\perp \f$ and rapidity. + /// @todo Introduce MomentumFilter objects for pT, ET, eta, y, etc. filtering, to avoid double-arg ambiguities + virtual Jets jets(double ptmin=0.0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=PSEUDORAPIDITY) const { + const Jets rawjets = _jets(ptmin); + Jets rtn; + foreach (const Jet& j, rawjets) { + const FourMomentum pj = j.momentum(); + if (!inRange(pj.pT(), ptmin, ptmax)) continue; + if (rapscheme == PSEUDORAPIDITY && !inRange(pj.eta(), rapmin, rapmax)) continue; + if (rapscheme == RAPIDITY && !inRange(pj.rapidity(), rapmin, rapmax)) continue; + rtn += j; + } + return rtn; + } - /// Get the jets, ordered by supplied sorting function object. + /// Get the jets, ordered by supplied sorting function object, with optional cuts on \f$ p_\perp \f$ and rapidity. + /// @todo Introduce MomentumFilter objects for pT, ET, eta, y, etc. filtering, to avoid double-arg ambiguities template <typename F> - Jets jets(F sorter, double ptmin=0.0) const { + Jets jets(F sorter, double ptmin=0.0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=PSEUDORAPIDITY) const { Jets js = jets(ptmin); if (sorter != 0) { std::sort(js.begin(), js.end(), sorter); @@ -92,22 +108,40 @@ return js; } - /// Get the jets, ordered by \f$ p_T \f$. - Jets jetsByPt(double ptmin=0.0) const { - return jets(cmpJetsByPt, ptmin); + /// Get the jets, ordered by \f$ p_T \f$, with optional cuts on \f$ p_\perp \f$ and rapidity. + /// @todo Introduce MomentumFilter objects for pT, ET, eta, y, etc. filtering, to avoid double-arg ambiguities + Jets jetsByPt(double ptmin=0.0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=PSEUDORAPIDITY) const { + return jets(cmpJetsByPt, ptmin, ptmax, rapmin, rapmax, rapscheme); } - /// Get the jets, ordered by \f$ E \f$. - Jets jetsByE(double ptmin=0.0) const { - return jets(cmpJetsByE, ptmin); + /// Get the jets, ordered by \f$ E \f$, with optional cuts on \f$ p_\perp \f$ and rapidity. + /// @todo Introduce MomentumFilter objects for pT, ET, eta, y, etc. filtering, to avoid double-arg ambiguities + Jets jetsByE(double ptmin=0.0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=PSEUDORAPIDITY) const { + return jets(cmpJetsByE, ptmin, ptmax, rapmin, rapmax, rapscheme); } - /// Get the jets, ordered by \f$ E_T \f$. - Jets jetsByEt(double ptmin=0.0) const { - return jets(cmpJetsByEt, ptmin); + /// Get the jets, ordered by \f$ E_T \f$, with optional cuts on \f$ p_\perp \f$ and rapidity. + /// @todo Introduce MomentumFilter objects for pT, ET, eta, y, etc. filtering, to avoid double-arg ambiguities + Jets jetsByEt(double ptmin=0.0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=PSEUDORAPIDITY) const { + return jets(cmpJetsByEt, ptmin, ptmax, rapmin, rapmax, rapscheme); } + protected: + + /// @brief Internal pure virtual method for getting jets in no guaranteed order. + /// An optional cut on min \f$ p_\perp \f$ is applied in this function, since that is + /// directly supported by FastJet and it seems a shame to not make use of that. But + /// all other jet cuts are applied at the @c ::jets() function level. + virtual Jets _jets(double ptmin=0.0) const = 0; + + public: /// Number of jets. Modified: trunk/include/Rivet/Projections/JetShape.hh ============================================================================== --- trunk/include/Rivet/Projections/JetShape.hh Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/include/Rivet/Projections/JetShape.hh Thu Nov 25 14:11:16 2010 (r2766) @@ -50,11 +50,18 @@ /// @name Constructors etc. //@{ - /// Constructor. + /// Constructor from histo range and number of bins. JetShape(const JetAlg& jetalg, - /// @todo Provide bin edges instead -- and util functions for making lin and log spaces - double rmin=0.0, double rmax=0.7, double interval=0.1, - DeltaRScheme distscheme=RAPIDITY); + double rmin, double rmax, size_t nbins, + double ptmin=0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=RAPIDITY); + + /// Constructor from vector of bin edges. + JetShape(const JetAlg& jetalg, vector<double> binedges, + double ptmin=0, double ptmax=MAXDOUBLE, + double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE, + RapScheme rapscheme=RAPIDITY); /// Clone on the heap. virtual const Projection* clone() const { @@ -77,28 +84,61 @@ /// Number of equidistant radius bins. size_t numBins() const { - return _nbins; + return _binedges.size() - 1; } /// \f$ r_\text{min} \f$ value. double rMin() const { - return _rmin; + return _binedges.front(); } /// \f$ r_\text{max} \f$ value. double rMax() const { - return _rmax; + return _binedges.back(); + } + + /// \f$ p_\perp^\text{min} \f$ value. + double ptMin() const { + return _ptcuts.first; + } + + /// \f$ p_\perp^\text{max} \f$ value. + double ptMax() const { + return _ptcuts.second; + } + + /// Central \f$ r \f$ value for bin @a rbin. + /// @todo This external indexing thing is a bit nasty... + double rBinMin(size_t rbin) const { + assert(inRange(rbin, 0, numBins())); + return _binedges[rbin]; + } + + /// Central \f$ r \f$ value for bin @a rbin. + /// @todo This external indexing thing is a bit nasty... + double rBinMax(size_t rbin) const { + assert(inRange(rbin, 0, numBins())); + return _binedges[rbin+1]; + } + + /// Central \f$ r \f$ value for bin @a rbin. + /// @todo This external indexing thing is a bit nasty... + double rBinMid(size_t rbin) const { + assert(inRange(rbin, 0, numBins())); + return (_binedges[rbin] + _binedges[rbin+1])/2.0; } /// Return value of differential jet shape profile histo bin. - /// @todo Remove this external indexing thing + /// @todo This external indexing thing is a bit nasty... double diffJetShape(size_t rbin) const { + assert(inRange(rbin, 0, numBins())); return _diffjetshapes[rbin]; } /// Return value of integrated jet shape profile histo bin. - /// @todo Remove this external indexing thing + /// @todo This external indexing thing is a bit nasty... double intJetShape(size_t rbin) const { + assert(inRange(rbin, 0, numBins())); double rtn = 0; for (size_t i = 0; i <= rbin; ++i) { rtn += _diffjetshapes[i]; @@ -108,8 +148,6 @@ /// @todo Provide int and diff jet shapes with some sort of area normalisation? - /// @todo Support some sort of jet pT binning? - // /// Return value of \f$ \Psi \f$ (integrated jet shape) at given radius for a \f$ p_T \f$ bin. // /// @todo Remove this external indexing thing // double psi(size_t pTbin) const { @@ -128,34 +166,32 @@ private: - /// @name The projected jet shapes + /// @name Jet shape parameters //@{ - /// Jet shape histo - vector<double> _diffjetshapes; + /// Vector of radius bin edges + vector<double> _binedges; - //@} + /// Lower and upper cuts on contributing jet \f$ p_\perp \f$. + pair<double, double> _ptcuts; + /// Lower and upper cuts on contributing jet (pseudo)rapidity. + pair<double, double> _rapcuts; - /// @name Jet shape parameters - //@{ - - /// Min radius (typically r=0) - double _rmin; + /// Rapidity scheme + RapScheme _rapscheme; - /// Max radius - double _rmax; + //@} - /// Length of radius interval - double _interval; - /// Rapidity scheme - DeltaRScheme _distscheme; + /// @name The projected jet shapes + //@{ - /// Number of bins - size_t _nbins; + /// Jet shape histo + vector<double> _diffjetshapes; //@} + }; Modified: trunk/src/Analyses/CDF_2005_S6217184.cc ============================================================================== --- trunk/src/Analyses/CDF_2005_S6217184.cc Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/src/Analyses/CDF_2005_S6217184.cc Thu Nov 25 14:11:16 2010 (r2766) @@ -5,7 +5,7 @@ #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" -#include "Rivet/Projections/ClosestJetShape.hh" +#include "Rivet/Projections/JetShape.hh" namespace Rivet { @@ -31,28 +31,29 @@ // Set up projections const FinalState fs(-2.0, 2.0); addProjection(fs, "FS"); - addProjection(FastJets(fs, FastJets::CDFMIDPOINT, 0.7), "Jets"); - - // Veto (anti)neutrinos, and muons with pT above 1.0 GeV - VisibleFinalState visfs(fs); - VetoedFinalState vfs(visfs); - vfs.addVetoPairDetail(MUON, 1.0*GeV, MAXDOUBLE); - addProjection(ClosestJetShape(vfs, _jetaxes, 0.0, 0.7, 0.1, 0.3), "JetShape"); + FastJets fj(fs, FastJets::CDFMIDPOINT, 0.7); + fj.useInvisibles(); + addProjection(fj, "Jets"); // Specify pT bins - _pTbins += 37.0, 45.0, 55.0, 63.0, 73.0, 84.0, 97.0, 112.0, 128.0, - 148.0, 166.0, 186.0, 208.0, 229.0, 250.0, 277.0, 304.0, 340.0, 380.0; + _ptedges = { 37.0, 45.0, 55.0, 63.0, 73.0, 84.0, 97.0, 112.0, 128.0, + 148.0, 166.0, 186.0, 208.0, 229.0, 250.0, 277.0, 304.0, 340.0, 380.0 }; - /// Book histograms - // 18 = 6x3 pT bins, one histogram each + // Register a jet shape projection and histogram for each pT bin for (size_t i = 0; i < 6; ++i) { for (size_t j = 0; j < 3; ++j) { - size_t k = i*3 + j; + const size_t k = i*3 + j; + stringstream ss; ss << "JetShape" << k; + const string pname = ss.str(); + _jsnames_pT[k] = pname; + const JetShape jsp(fj, 0.0, 0.7, 6, _ptedges[k], _ptedges[k+1], 0.1, 0.7, RAPIDITY); + addProjection(jsp, pname); _profhistRho_pT[k] = bookProfile1D(i+1, 1, j+1); _profhistPsi_pT[k] = bookProfile1D(6+i+1, 1, j+1); } } + // Final histo _profhistPsi = bookProfile1D(13, 1, 1); } @@ -62,38 +63,28 @@ void analyze(const Event& event) { // Get jets and require at least one to pass pT and y cuts - const Jets& jets = applyProjection<FastJets>(event, "Jets").jetsByPt(); + /// @todo Need a different rapidity cut -- discontinuous in |y| + const Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(37, 380, 0.1, 0.7); getLog() << Log::DEBUG << "Jet multiplicity before cuts = " << jets.size() << endl; - - // Determine the central jet axes - _jetaxes.clear(); - foreach (const Jet& jt, jets) { - const FourMomentum pj = jt.momentum(); - if (inRange(pj.pT()/GeV, 37.0, 380.0) && inRange(fabs(pj.rapidity()), 0.1, 0.7)) { - _jetaxes.push_back(jt.momentum()); - } + if (jets.size() == 0) { + MSG_DEBUG("No jets found in required pT range"); + vetoEvent(event); } - if (_jetaxes.empty()) vetoEvent; // Calculate and histogram jet shapes const double weight = event.weight(); - const ClosestJetShape& js = applyProjection<ClosestJetShape>(event, "JetShape"); - /// @todo Use BinnedHistogram, for collections of histos each for a range of values of an extra variable - for (size_t jind = 0; jind < _jetaxes.size(); ++jind) { - for (size_t ipT = 0; ipT < 18; ++ipT) { - if (_jetaxes[jind].pT() > _pTbins[ipT] && _jetaxes[jind].pT() <= _pTbins[ipT+1]) { - for (size_t rbin = 0; rbin < js.numBins(); ++rbin) { - const double rad_Rho = js.rMin() + (rbin+0.5)*js.interval(); - _profhistRho_pT[ipT]->fill(rad_Rho/0.7, (0.7/1.0)*js.diffJetShape(jind, rbin), weight); - /// @todo Calc int histos from diff histos - const double rad_Psi = js.rMin() +(rbin+1.0)*js.interval(); - _profhistPsi_pT[ipT]->fill(rad_Psi/0.7, js.intJetShape(jind, rbin), weight); - } - /// @todo Calc int histos from diff histos - _profhistPsi->fill((_pTbins[ipT] + _pTbins[ipT+1])/2.0, js.psi(jind), weight); - } + for (size_t ipT = 0; ipT < 18; ++ipT) { + const JetShape& jsipt = applyProjection<JetShape>(event, _jsnames_pT[ipT]); + for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) { + const double r_rho = jsipt.rBinMid(rbin); + _profhistRho_pT[ipT]->fill(r_rho/0.7, (0.7/1.0)*jsipt.diffJetShape(rbin), weight); + const double r_Psi = jsipt.rBinMax(rbin); + _profhistPsi_pT[ipT]->fill(r_Psi/0.7, jsipt.intJetShape(rbin), weight); } + // Final histo is 1 - Psi(0.3/R) as a function of jet pT bin + const double ptmid = (_ptedges[ipT] + _ptedges[ipT+1])/2.0; + _profhistPsi->fill(ptmid/GeV, jsipt.intJetShape(2), weight); } } @@ -112,11 +103,12 @@ /// @name Analysis data //@{ - /// Vector of jet axes - vector<FourMomentum> _jetaxes; + /// Jet \f$ p_\perp\f$ bins. + double _ptedges[19]; + + /// JetShape projection name for each \f$p_\perp\f$ bin. + string _jsnames_pT[18]; - /// \f$p_\perp\f$ bins to be distinguished during analysis - vector<double> _pTbins; //@} Modified: trunk/src/Projections/ClosestJetShape.cc ============================================================================== --- trunk/src/Projections/ClosestJetShape.cc Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/src/Projections/ClosestJetShape.cc Thu Nov 25 14:11:16 2010 (r2766) @@ -9,7 +9,7 @@ ClosestJetShape::ClosestJetShape(const FinalState& fs, const vector<FourMomentum>& jetaxes, double rmin, double rmax, double interval, - double r1minPsi, DeltaRScheme distscheme) + double r1minPsi, RapScheme distscheme) : _jetaxes(jetaxes), _rmin(rmin), _rmax(rmax), _interval(interval), _r1minPsi(r1minPsi), _distscheme(distscheme) Modified: trunk/src/Projections/FastJets.cc ============================================================================== --- trunk/src/Projections/FastJets.cc Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/src/Projections/FastJets.cc Thu Nov 25 14:11:16 2010 (r2766) @@ -156,25 +156,25 @@ } - Jets FastJets::jets(double ptmin) const { + Jets FastJets::_jets(double ptmin) const { Jets rtn = _pseudojetsToJets(pseudoJets(ptmin)); return rtn; } - Jets FastJets::jetsByPt(double ptmin) const { - return _pseudojetsToJets(pseudoJetsByPt(ptmin)); - } + // Jets FastJets::jetsByPt(double ptmin) const { + // return _pseudojetsToJets(pseudoJetsByPt(ptmin)); + // } - Jets FastJets::jetsByE(double ptmin) const { - return _pseudojetsToJets(pseudoJetsByE(ptmin)); - } + // Jets FastJets::jetsByE(double ptmin) const { + // return _pseudojetsToJets(pseudoJetsByE(ptmin)); + // } - Jets FastJets::jetsByRapidity(double ptmin) const { - return _pseudojetsToJets(pseudoJetsByRapidity(ptmin)); - } + // Jets FastJets::jetsByRapidity(double ptmin) const { + // return _pseudojetsToJets(pseudoJetsByRapidity(ptmin)); + // } PseudoJets FastJets::pseudoJets(double ptmin) const { Modified: trunk/src/Projections/JetShape.cc ============================================================================== --- trunk/src/Projections/JetShape.cc Thu Nov 25 02:51:41 2010 (r2765) +++ trunk/src/Projections/JetShape.cc Thu Nov 25 14:11:16 2010 (r2766) @@ -5,28 +5,53 @@ namespace Rivet { - /// Constructor. + // Constructor. JetShape::JetShape(const JetAlg& jetalg, - double rmin, double rmax, double interval, - DeltaRScheme distscheme) - : _rmin(rmin), _rmax(rmax), _interval(interval), - _distscheme(distscheme) + double rmin, double rmax, size_t nbins, + double ptmin, double ptmax, + double rapmin, double rapmax, + RapScheme rapscheme) + : _rapscheme(rapscheme) { setName("JetShape"); - _nbins = int(round((rmax-rmin)/interval)); + _binedges = linspace(rmin, rmax, nbins); + _ptcuts = make_pair(ptmin, ptmax); + addProjection(jetalg, "Jets"); + } + + + // Constructor. + JetShape::JetShape(const JetAlg& jetalg, + vector<double> binedges, + double ptmin, double ptmax, + double rapmin, double rapmax, + RapScheme rapscheme) + : _binedges(binedges), _rapscheme(rapscheme) + { + setName("JetShape"); + _ptcuts = make_pair(ptmin, ptmax); addProjection(jetalg, "Jets"); } int JetShape::compare(const Projection& p) const { - PCmp jcmp = mkNamedPCmp(p, "Jets"); - return jcmp; - /// @todo Also compare bin edges. + const int jcmp = mkNamedPCmp(p, "Jets"); + if (jcmp != EQUIVALENT) return jcmp; + const JetShape& other = pcast<JetShape>(p); + const int ptcmp = cmp(ptMin(), other.ptMin()) || cmp(ptMax(), other.ptMax()); + if (ptcmp != EQUIVALENT) return ptcmp; + int bincmp = cmp(numBins(), other.numBins()); + if (bincmp != EQUIVALENT) return bincmp; + for (size_t i = 0; i < _binedges.size(); ++i) { + bincmp = cmp(_binedges[i], other._binedges[i]); + if (bincmp != EQUIVALENT) return bincmp; + } + return EQUIVALENT; } void JetShape::clear() { - _diffjetshapes = vector<double>(_nbins, 0.0); + _diffjetshapes = vector<double>(numBins(), 0.0); } @@ -35,20 +60,31 @@ foreach (const Jet& j, jets) { FourMomentum pj = j.momentum(); + // These cuts are already applied when getting the jets in ::project(), but we check again for direct calc use. + if (!inRange(pj.pT(), _ptcuts)) continue; + if (_rapscheme == PSEUDORAPIDITY && !inRange(pj.eta(), _rapcuts)) continue; + if (_rapscheme == RAPIDITY && !inRange(pj.rapidity(), _rapcuts)) continue; foreach (const Particle& p, j.particles()) { - const double dR = deltaR(pj, p.momentum()); - size_t dRindex = int(floor((dR-_rmin)/(_rmax - _rmin))); - assert(dRindex < _nbins); + const double dR = deltaR(pj, p.momentum(), _rapscheme); + if (!inRange(dR, _binedges.front(), _binedges.back())) continue; //< Out of histo range + size_t dRindex = -1; + for (size_t i = 1; i < _binedges.size(); ++i) { + if (dR < _binedges[i]) { + dRindex = i-1; + break; + } + } + assert(inRange(dRindex, 0, numBins())); _diffjetshapes[dRindex] += p.momentum().pT(); } } // Normalize to total pT double integral = 0.0; - for (size_t i = 0; i < _nbins; ++i) { + for (size_t i = 0; i < numBins(); ++i) { integral += _diffjetshapes[i]; } - for (size_t i = 0; i < _nbins; ++i) { + for (size_t i = 0; i < numBins(); ++i) { _diffjetshapes[i] /= integral; } @@ -56,7 +92,8 @@ void JetShape::project(const Event& e) { - const Jets jets = applyProjection<JetAlg>(e, "Jets").jets(); + const Jets jets = applyProjection<JetAlg>(e, "Jets").jets(_ptcuts.first, _ptcuts.second, + _rapcuts.first, _rapcuts.second, _rapscheme); calc(jets); }
More information about the Rivet-svn mailing list |