|
[Rivet-svn] r1934 - in trunk: . include/Rivet include/Rivet/Projections src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgMon Oct 19 15:09:01 BST 2009
Author: hoeth Date: Mon Oct 19 15:09:01 2009 New Revision: 1934 Log: new MergedFinalState projection. As in the NeutralFinalState, the compare method is not yet working properly. Added: trunk/include/Rivet/Projections/MergedFinalState.hh trunk/src/Projections/MergedFinalState.cc Modified: trunk/ChangeLog trunk/include/Rivet/Makefile.am trunk/include/Rivet/Particle.hh trunk/include/Rivet/ParticleBase.hh trunk/include/Rivet/Projections/NeutralFinalState.hh trunk/src/Analyses/STAR_2009_UE_HELEN.cc trunk/src/Projections/Makefile.am trunk/src/Projections/NeutralFinalState.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/ChangeLog Mon Oct 19 15:09:01 2009 (r1934) @@ -10,6 +10,10 @@ state takes E_T instead of p_T as argument (makes more sense for neutral particles). The compare() method does not yet work as expected (E_T comparison still missing). + * Adding new MergedFinalState projection. This merges two final + states, removing duplicate particles. Duplicates are identified + by looking at the genParticle(), so users need to take care of + any manually added particles themselves. 2009-10-17 Andy Buckley <andy at insectnation.org> Modified: trunk/include/Rivet/Makefile.am ============================================================================== --- trunk/include/Rivet/Makefile.am Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/include/Rivet/Makefile.am Mon Oct 19 15:09:01 2009 (r1934) @@ -61,6 +61,7 @@ Projections/KtJets.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ + Projections/MergedFinalState.hh \ Projections/Multiplicity.hh \ Projections/NeutralFinalState.hh \ Projections/ParisiTensor.hh \ Modified: trunk/include/Rivet/Particle.hh ============================================================================== --- trunk/include/Rivet/Particle.hh Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/include/Rivet/Particle.hh Mon Oct 19 15:09:01 2009 (r1934) @@ -9,17 +9,17 @@ namespace Rivet { - + /// Representation of particles from a HepMC::GenEvent. class Particle : public ParticleBase { - public: + public: /// Default constructor. Particle() : ParticleBase(), _original(0), _id(0), _momentum(), _mass(0.0) { } - + /// Constructor from a HepMC GenParticle. Particle(const GenParticle& gp) : ParticleBase(), _original(&gp), _id(gp.pdg_id()), @@ -29,16 +29,16 @@ public: /// Get a const reference to the original GenParticle. - const GenParticle& genParticle() const { - assert(_original); - return *_original; + const GenParticle& genParticle() const { + assert(_original); + return *_original; } - + /// Check if the particle corresponds to a GenParticle. bool hasGenParticle() const { return bool(_original); } - + /// The PDG ID code for this Particle. const long pdgId() const { return _id; } @@ -54,11 +54,11 @@ /// Set the momentum of this Particle. Particle& setMomentum(const FourMomentum& momentum) { _momentum = momentum; return *this; } - + /// The mass of this Particle. const double mass() const { return _mass; } - + bool hasAncestor(long pdg_id) const { GenVertex* prodVtx = genParticle().production_vertex(); if (prodVtx == 0) return false; @@ -68,24 +68,24 @@ if ((*ancestor)->pdg_id() == pdg_id) { return true; } - } + } return false; } - + private: /// A pointer to the original GenParticle from which this Particle is projected. const GenParticle* _original; - + /// The PDG ID code for this Particle. long _id; - + /// The momentum of this projection of the Particle. FourMomentum _momentum; - + /// The mass of this Particle, stored for numerical hygiene. - double _mass; + double _mass; }; @@ -135,7 +135,7 @@ - + } #endif Modified: trunk/include/Rivet/ParticleBase.hh ============================================================================== --- trunk/include/Rivet/ParticleBase.hh Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/include/Rivet/ParticleBase.hh Mon Oct 19 15:09:01 2009 (r1934) @@ -13,100 +13,100 @@ #include "Rivet/Math/Vectors.hh" namespace Rivet{ - + class ParticleBase{ - + public: - + ParticleBase(){ }; - + virtual ~ParticleBase(){}; - + virtual FourMomentum &momentum()=0;//{return _momentum;}; - + virtual const FourMomentum &momentum() const =0;//{return _momentum;}; - - /// struct for sorting by increasing transverse momentum in stl set, sort etc. - + + /// struct for sorting by increasing transverse momentum in stl set, sort etc. + struct byPTAscending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ double pt2left = left.momentum().pT2(); double pt2right = right.momentum().pT2(); return pt2left < pt2right; } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - - /// struct for sorting by decreasing transverse momentum in stl set, sort etc. - + + /// struct for sorting by decreasing transverse momentum in stl set, sort etc. + struct byPTDescending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ return byPTAscending()(right, left); } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - + /// struct for sorting by increasing transverse energy - + struct byETAscending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ double pt2left = left.momentum().Et2(); double pt2right = right.momentum().Et2(); return pt2left < pt2right; } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - + /// struct for sorting by decreasing transverse energy - + struct byETDescending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ return byETAscending()(right, left); } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - + /// struct for sorting by increaing energy - + struct byEAscending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ double pt2left = left.momentum().E(); double pt2right = right.momentum().E(); return pt2left < pt2right; } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - + /// struct for sorting by decreasing energy - + struct byEDescending{ bool operator()(const ParticleBase &left, const ParticleBase &right) const{ return byEAscending()(right, left); } - + bool operator()(const ParticleBase *left, const ParticleBase *right) const{ return (*this)(*left, *right); } }; - + protected: - + }; } Added: trunk/include/Rivet/Projections/MergedFinalState.hh ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/include/Rivet/Projections/MergedFinalState.hh Mon Oct 19 15:09:01 2009 (r1934) @@ -0,0 +1,47 @@ +// -*- C++ -*- +#ifndef RIVET_MergedFinalState_HH +#define RIVET_MergedFinalState_HH + +#include "Rivet/Tools/Logging.hh" +#include "Rivet/Rivet.hh" +#include "Rivet/Particle.hh" +#include "Rivet/Event.hh" +#include "Rivet/Projection.hh" +#include "Rivet/Projections/FinalState.hh" + +namespace Rivet { + + + /// Project only merged final state particles. + class MergedFinalState : public FinalState { + + public: + + /// @name Constructors + //@{ + MergedFinalState(const FinalState& fspa, const FinalState& fspb) { + setName("MergedFinalState"); + addProjection(fspa, "FSA"); + addProjection(fspb, "FSB"); + } + + /// Clone on the heap. + virtual const Projection* clone() const { + return new MergedFinalState(*this); + } + //@} + + protected: + + /// Apply the projection on the supplied event. + void project(const Event& e); + + /// Compare projections. + int compare(const Projection& p) const; + }; + + +} + + +#endif Modified: trunk/include/Rivet/Projections/NeutralFinalState.hh ============================================================================== --- trunk/include/Rivet/Projections/NeutralFinalState.hh Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/include/Rivet/Projections/NeutralFinalState.hh Mon Oct 19 15:09:01 2009 (r1934) @@ -16,18 +16,18 @@ class NeutralFinalState : public FinalState { public: - + /// @name Constructors //@{ NeutralFinalState(const FinalState& fsp) : _Etmin(0.0*GeV) { setName("NeutralFinalState"); addProjection(fsp, "FS"); } - + NeutralFinalState(double mineta = -MAXRAPIDITY, double maxeta = MAXRAPIDITY, double minEt = 0.0*GeV) : _Etmin(minEt) - { + { setName("NeutralFinalState"); addProjection(FinalState(mineta, maxeta, 0.0*GeV), "FS"); } @@ -39,10 +39,10 @@ //@} protected: - + /// Apply the projection on the supplied event. void project(const Event& e); - + /// The minimum allowed transverse energy. double _Etmin; @@ -50,7 +50,7 @@ int compare(const Projection& p) const; }; - + } Modified: trunk/src/Analyses/STAR_2009_UE_HELEN.cc ============================================================================== --- trunk/src/Analyses/STAR_2009_UE_HELEN.cc Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/src/Analyses/STAR_2009_UE_HELEN.cc Mon Oct 19 15:09:01 2009 (r1934) @@ -5,6 +5,7 @@ #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" +#include "fastjet/SISConePlugin.hh" namespace Rivet { @@ -28,10 +29,11 @@ void init() { // Final state for the jet finding - const FinalState fsj(-4.0, 4.0, 0.0*GeV); + const FinalState fsj(-1, 1, 0.2*GeV); addProjection(fsj, "FSJ"); - /// @todo Check which merging parameter they've been using - addProjection(FastJets(fsj, FastJets::SISCONE, 0.7), "AllJets"); + // Split-merge is 0.75, so we need to initialize the plugin manually: + // R = 0.7, overlap_threshold = 0.75 + addProjection(FastJets(fsj, fastjet::SISConePlugin(0.7, 0.75)), "AllJets"); // Charged final state for the distributions const ChargedFinalState cfs(-1.0, 1.0, 0.2*GeV); @@ -61,8 +63,8 @@ jets.push_back(jet); } - // We require the leading jet to be within |eta|<2 - if (jets.size() < 1 || fabs(jets[0].momentum().eta()) >= 2) { + // We require the leading jet to be within |eta|<(1-R)=0.3 + if (jets.size() < 1 || fabs(jets[0].momentum().eta()) >= 0.3) { getLog() << Log::DEBUG << "Failed jet cut" << endl; vetoEvent; } Modified: trunk/src/Projections/Makefile.am ============================================================================== --- trunk/src/Projections/Makefile.am Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/src/Projections/Makefile.am Mon Oct 19 15:09:01 2009 (r1934) @@ -20,6 +20,7 @@ JetShape.cc \ LeadingParticlesFinalState.cc \ LossyFinalState.cc \ + MergedFinalState.cc \ Multiplicity.cc \ NeutralFinalState.cc \ ParisiTensor.cc \ Added: trunk/src/Projections/MergedFinalState.cc ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/src/Projections/MergedFinalState.cc Mon Oct 19 15:09:01 2009 (r1934) @@ -0,0 +1,48 @@ +// -*- C++ -*- +#include "Rivet/Rivet.hh" +#include "Rivet/Projections/MergedFinalState.hh" +#include "Rivet/Tools/ParticleIdUtils.hh" +#include "Rivet/Cmp.hh" +#include <algorithm> + +namespace Rivet { + + + int MergedFinalState::compare(const Projection& p) const { + /// @todo: This needs to be fixed!! We need to + // - check if the size matches (if it doesn't, they are not equal) + // - for equal sizes compare the elements (yuck!) + return mkNamedPCmp(p, "FSA"); + } + + + void MergedFinalState::project(const Event& e) { + const FinalState& fsa = applyProjection<FinalState>(e, "FSA"); + const FinalState& fsb = applyProjection<FinalState>(e, "FSB"); + _theParticles.clear(); + foreach (const Particle& pa, fsa.particles()){ + _theParticles.push_back(pa); + } + foreach (const Particle& pb, fsb.particles()){ + const GenParticle* originalb = &(pb.genParticle()); + bool notfound = true; + foreach (const Particle& pa, fsa.particles()){ + const GenParticle* originala = &(pa.genParticle()); + if (originala = originalb) { + notfound = false; + break; + } + } + if (notfound) { + _theParticles.push_back(pb); + } + } + getLog() << Log::DEBUG << "Number of particles in the two final states to be merged: = \n" + << " 1st final state = " << fsa.particles().size() << endl + << " 2nd final state = " << fsb.particles().size() << endl; + getLog() << Log::DEBUG << "Number of merged final-state particles = " + << _theParticles.size() << endl; + } + + +} Modified: trunk/src/Projections/NeutralFinalState.cc ============================================================================== --- trunk/src/Projections/NeutralFinalState.cc Mon Oct 19 15:04:43 2009 (r1933) +++ trunk/src/Projections/NeutralFinalState.cc Mon Oct 19 15:09:01 2009 (r1934) @@ -22,14 +22,14 @@ _theParticles.push_back(p); if (getLog().isActive(Log::TRACE)) { getLog() << Log::TRACE - << "Selected: ID = " << p.pdgId() - << ", Et = " << p.momentum().Et() - << ", eta = " << p.momentum().eta() + << "Selected: ID = " << p.pdgId() + << ", Et = " << p.momentum().Et() + << ", eta = " << p.momentum().eta() << ", charge = " << PID::threeCharge(p.pdgId())/3.0 << endl; } } } - getLog() << Log::DEBUG << "Number of neutral final-state particles = " + getLog() << Log::DEBUG << "Number of neutral final-state particles = " << _theParticles.size() << endl; }
More information about the Rivet-svn mailing list |