|
[Rivet-svn] r4211 - in trunk: . include/Rivet src src/Analyses src/Core src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Mar 7 17:20:02 GMT 2013
Author: buckley Date: Thu Mar 7 17:20:02 2013 New Revision: 4211 Log: Converting Event to use a different implementation of original and modified GenParticles, and to manage the memory in a more future-proof way. Event::genParticle() now returns a const pointer rather than a reference, to signal that the user is leaving the happy pastures of 'normal' Rivet behind. Converting Particle::genParticle() to return a const pointer rather than a reference, for the same reason as below (+ consistency within Rivet and with the HepMC pointer-centric coding design). Adding new Particle(const GenParticle*) constructor. Updating util classes, projections, and analyses to deal with the HepMC return value changes. Modified: trunk/ChangeLog trunk/include/Rivet/Event.hh trunk/include/Rivet/Particle.hh trunk/src/Analyses/ARGUS_1993_S2653028.cc trunk/src/Analyses/ARGUS_1993_S2669951.cc trunk/src/Analyses/ARGUS_1993_S2789213.cc trunk/src/Analyses/ATLAS_2011_I9035664.cc trunk/src/Analyses/ATLAS_2011_I944826.cc trunk/src/Analyses/BABAR_2003_I593379.cc trunk/src/Analyses/BABAR_2005_S6181155.cc trunk/src/Analyses/BABAR_2007_S7266081.cc trunk/src/Analyses/BELLE_2001_S4598261.cc trunk/src/Analyses/BELLE_2006_S6265367.cc trunk/src/Analyses/CDF_2006_S6653332.cc trunk/src/Analyses/CDF_2008_S7540469.cc trunk/src/Analyses/D0_2008_S7863608.cc trunk/src/Analyses/D0_2009_S8349509.cc trunk/src/Analyses/D0_2010_S8570965.cc trunk/src/Analyses/H1_1994_S2919893.cc trunk/src/Analyses/H1_2000_S4129130.cc trunk/src/Analyses/LHCB_2010_S8758301.cc trunk/src/Analyses/LHCB_2011_I917009.cc trunk/src/Analyses/LHCB_2012_I1119400.cc trunk/src/Analyses/MC_PDFS.cc trunk/src/Analyses/MC_PRINTEVENT.cc trunk/src/Analyses/MC_XS.cc trunk/src/Analyses/STAR_2006_S6870392.cc trunk/src/Analyses/STAR_2008_S7993412.cc trunk/src/Core/Event.cc trunk/src/Core/Jet.cc trunk/src/Core/Particle.cc trunk/src/Makefile.am trunk/src/Projections/Beam.cc trunk/src/Projections/DISFinalState.cc trunk/src/Projections/DISLepton.cc trunk/src/Projections/FinalState.cc trunk/src/Projections/InvMassFinalState.cc trunk/src/Projections/LeadingParticlesFinalState.cc trunk/src/Projections/MergedFinalState.cc trunk/src/Projections/VetoedFinalState.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/ChangeLog Thu Mar 7 17:20:02 2013 (r4211) @@ -1,5 +1,21 @@ 2013-03-07 Andy Buckley <andy.buckley at cern.ch> + * Updating util classes, projections, and analyses to deal with + the HepMC return value changes. + + * Adding new Particle(const GenParticle*) constructor. + + * Converting Particle::genParticle() to return a const pointer + rather than a reference, for the same reason as below (+ + consistency within Rivet and with the HepMC pointer-centric coding + design). + + * Converting Event to use a different implementation of original + and modified GenParticles, and to manage the memory in a more + future-proof way. Event::genParticle() now returns a const pointer + rather than a reference, to signal that the user is leaving the + happy pastures of 'normal' Rivet behind. + * Adding a Particles typedef by analogy to Jets, and in preference to the cumbersome ParticleVector. Modified: trunk/include/Rivet/Event.hh ============================================================================== --- trunk/include/Rivet/Event.hh Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/include/Rivet/Event.hh Thu Mar 7 17:20:02 2013 (r4211) @@ -25,37 +25,49 @@ /// @name Standard constructors and destructors. //@{ - /// The default constructor. - Event(const GenEvent& ge); - - /// The copy constructor. - Event(const Event& e); - - /// The destructor - ~Event(); + /// Constructor from a HepMC GenEvent reference + Event(const GenEvent& ge) + : _originalGenEvent(&ge), _genEvent(ge) + { _init(ge); } + + /// Constructor from a HepMC GenEvent pointer + Event(const GenEvent* ge) + : _originalGenEvent(ge), _genEvent(*ge) + { assert(ge); _init(*ge); } + + /// Copy constructor + Event(const Event& e) + : _originalGenEvent(e._originalGenEvent), _genEvent(e._genEvent) + { } //@} public: - /// Return the generated event obtained from an external event generator. - const GenEvent& genEvent() const; + /// The generated event obtained from an external event generator + const GenEvent* genEvent() const { + return &_genEvent; + } - /// The weight associated with the event. + /// @brief The generation weight associated with the event + /// + /// @todo This needs to be revisited when we finally add the mechanism to + /// support NLO counter-events and weight vectors. double weight() const { - return _weight; + return (!_genEvent.weights().empty()) ? _genEvent.weights()[0] : 1.0; } public: - /// Add a projection \a p to this Event. If an equivalent Projection - /// has been applied before, the Projection::project(const Event &) - /// of \a p is not called and a reference to the previous equivalent - /// projection is returned. If no previous Projection was found, the - /// Projection::project(const Event &) of \a p is called and a - /// reference to p is returned. + /// @brief Add a projection @a p to this Event. + /// + /// If an equivalent Projection has been applied before, the + /// Projection::project(const Event&) of @a p is not called and a reference + /// to the previous equivalent projection is returned. If no previous + /// Projection was found, the Projection::project(const Event&) of @a p is + /// called and a reference to @a p is returned. template <typename PROJ> const PROJ& applyProjection(PROJ& p) const { const Projection* cpp(&p); @@ -81,22 +93,41 @@ private: - void _geNormAlignment(); + /// @brief Actual (shared) implementation of the constructors from GenEvents + /// + /// @todo When we can require C++11, constructors can call each other and + /// this can be removed. + void _init(const GenEvent& ge); + + /// @brief Convert the GenEvent to use conventional alignment + /// + /// For example, FHerwig only produces DIS events in the unconventional + /// hadron-lepton orientation and has to be corrected for DIS analysis + /// portability. + void _geNormAlignment(); - /// @brief The generated event, obtained from an external generator. - /// Note that it is only externally accessible via a const reference. - GenEvent const& _genEvent; - - GenEvent* _modGenEvent; + /// @brief The generated event, as obtained from an external generator. + /// + /// This is the original GenEvent. In practise the version seen by users + /// will often/always be a modified one. + /// + /// @todo Provide access to this via an Event::originalGenEvent() method? If requested... + const GenEvent* _originalGenEvent; + + /// @brief The GenEvent used by Rivet analysis projections etc. + /// + /// This version may be rotated to a "normal" alignment, have + /// generator-specific particles stripped out, etc. If an analysis is + /// affected by these modifications, it is probably an unphysical analysis! + /// + /// Stored as a non-pointer since it may get overwritten, and memory for + /// copying and cleanup is much neater this way. + GenEvent _genEvent; /// The set of Projection objects applied so far. mutable std::set<ConstProjectionPtr> _projections; - /// @brief The generation weight associated with the event. - /// Usually 1.0. Only copied from the HepMC event once, at construction time. - double _weight; - }; Modified: trunk/include/Rivet/Particle.hh ============================================================================== --- trunk/include/Rivet/Particle.hh Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/include/Rivet/Particle.hh Thu Mar 7 17:20:02 2013 (r4211) @@ -36,19 +36,19 @@ _momentum(gp.momentum()) { } + /// Constructor from a HepMC GenParticle pointer. + Particle(const GenParticle* gp) + : ParticleBase(), + _original(gp), _id(gp->pdg_id()), + _momentum(gp->momentum()) + { } + public: /// Get a const reference to the original GenParticle. - const GenParticle& genParticle() const { - assert(_original); - return *_original; - } - - - /// Check if the particle corresponds to a GenParticle. - bool hasGenParticle() const { - return bool(_original); + const GenParticle* genParticle() const { + return _original; } @@ -79,12 +79,13 @@ return momentum().mass(); } + /// @todo Enable? // /// The charge of this Particle. // double charge() const { // return PID::charge(*this); // } - // /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). + /// @todo Enable? // /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). // int threeCharge() const { // return PID::threeCharge(*this); // } Modified: trunk/src/Analyses/ARGUS_1993_S2653028.cc ============================================================================== --- trunk/src/Analyses/ARGUS_1993_S2653028.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/ARGUS_1993_S2653028.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -23,17 +23,17 @@ void analyze(const Event& e) { const double weight = e.weight(); - const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); - - // find the upsilons + // Find the upsilons ParticleVector upsilons; - // first in unstable final state - foreach (const Particle& p, ufs.particles()) - if(p.pdgId()==300553) upsilons.push_back(p); - // then in whole event if fails - if(upsilons.empty()) { + // First in unstable final state + const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); + foreach (const Particle& p, ufs.particles()) { + if (p.pdgId() == 300553) upsilons.push_back(p); + } + // Then in whole event if that failed + if (upsilons.empty()) { foreach (GenParticle* p, Rivet::particles(e.genEvent())) { - if(p->pdg_id()!=300553) continue; + if (p->pdg_id() != 300553) continue; const GenVertex* pv = p->production_vertex(); bool passed = true; if (pv) { @@ -45,17 +45,16 @@ } } } - if(passed) upsilons.push_back(Particle(*p)); + if (passed) upsilons.push_back(Particle(*p)); } } - // find an upsilons + // Find an upsilon foreach (const Particle& p, upsilons) { _weightSum += weight; vector<GenParticle *> pionsA,pionsB,protonsA,protonsB,kaons; - // find the decay products we want - findDecayProducts(p.genParticle(),pionsA,pionsB, - protonsA,protonsB,kaons); + // Find the decay products we want + findDecayProducts(p.genParticle(), pionsA, pionsB, protonsA, protonsB, kaons); LorentzTransform cms_boost; if(p.momentum().vector3().mod()>0.001) cms_boost = LorentzTransform(-p.momentum().boostVector()); @@ -96,7 +95,7 @@ } // analyze void finalize() { - if(_weightSum>0.) { + if (_weightSum > 0.) { scale(_histPiA, 1./_weightSum); scale(_histPiB, 1./_weightSum); scale(_histKA , 1./_weightSum); @@ -133,7 +132,7 @@ private: //@{ - // count of weights + /// Count of weights double _weightSum; // Histograms // spectra @@ -151,28 +150,26 @@ Histo1DPtr _multpB; //@} - void findDecayProducts(const GenParticle & p, - vector<GenParticle *> & pionsA, - vector<GenParticle *> & pionsB, - vector<GenParticle *> & protonsA, - vector<GenParticle *> & protonsB, - vector<GenParticle *> & kaons) { - int parentId = p.pdg_id(); - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + void findDecayProducts(const GenParticle* p, + vector<GenParticle*>& pionsA, vector<GenParticle*>& pionsB, + vector<GenParticle*>& protonsA, vector<GenParticle*>& protonsB, + vector<GenParticle*>& kaons) + { + int parentId = p->pdg_id(); + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = abs((*pp)->pdg_id()); - if(id==PIPLUS) { - if(parentId != LAMBDA && parentId != K0S) { + if (id==PIPLUS) { + if (parentId != LAMBDA && parentId != K0S) { pionsA.push_back(*pp); pionsB.push_back(*pp); } else pionsB.push_back(*pp); } - else if(id==PROTON) { - if(parentId != LAMBDA && parentId != K0S) { + else if (id==PROTON) { + if (parentId != LAMBDA && parentId != K0S) { protonsA.push_back(*pp); protonsB.push_back(*pp); } @@ -183,8 +180,7 @@ kaons.push_back(*pp); } else if((*pp)->end_vertex()) - findDecayProducts(**pp,pionsA,pionsB, - protonsA,protonsB,kaons); + findDecayProducts(*pp, pionsA, pionsB, protonsA, protonsB, kaons); } } }; Modified: trunk/src/Analyses/ARGUS_1993_S2669951.cc ============================================================================== --- trunk/src/Analyses/ARGUS_1993_S2669951.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/ARGUS_1993_S2669951.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -34,9 +34,9 @@ ParticleVector upsilons; // first in unstable final state foreach (const Particle& p, ufs.particles()) - if(p.pdgId()==553 || p.pdgId()==100553 ) upsilons.push_back(p); + if (p.pdgId()==553 || p.pdgId()==100553 ) upsilons.push_back(p); // then in whole event if fails - if(upsilons.empty()) { + if (upsilons.empty()) { foreach (GenParticle* p, Rivet::particles(e.genEvent())) { if( p->pdg_id() != 553 && p->pdg_id() != 100553 ) continue; const GenVertex* pv = p->production_vertex(); @@ -85,7 +85,7 @@ _weightSum_Ups2 += weight; ParticleVector unstable; // find the decay products we want - findDecayProducts(ups.genParticle(),unstable); + findDecayProducts(ups.genParticle(), unstable); LorentzTransform cms_boost; if(ups.momentum().vector3().mod()>0.001) cms_boost = LorentzTransform(-ups.momentum().boostVector()); @@ -96,17 +96,17 @@ FourMomentum p2 = cms_boost.transform(p.momentum()); double xp = 2.*p2.t()/mass; double beta = p2.vector3().mod()/p2.t(); - if(id==9010221) { + if (id==9010221) { if(parentId==553) _hist_Ups1_f0->fill(xp,weight/beta); else _hist_Ups2_f0->fill(xp,weight/beta); ++nf0; } - else if(id==331) { + else if (id==331) { if(xp>0.35) ++nEtaA; ++nEtaB; } } - if(parentId==553) { + if (parentId==553) { _count_f0[0] += nf0*weight; _count_etaPrime_highZ[0] += nEtaA*weight; _count_etaPrime_allZ[0] += nEtaB*weight; @@ -176,9 +176,10 @@ _hist_cont_f0 = bookHisto1D( 2,1,1); _hist_Ups1_f0 = bookHisto1D( 3,1,1); _hist_Ups2_f0 = bookHisto1D( 4,1,1); - } // init + } - private: + + private: //@{ vector<double> _count_etaPrime_highZ; @@ -189,26 +190,26 @@ Histo1DPtr _hist_Ups1_f0; Histo1DPtr _hist_Ups2_f0; - // count of weights double _weightSum_cont,_weightSum_Ups1,_weightSum_Ups2; //@} - void findDecayProducts(const GenParticle & p, + + void findDecayProducts(const GenParticle* p, ParticleVector & unstable) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = abs((*pp)->pdg_id()); - if(id == 331 || id == 9010221 ) { - unstable.push_back(Particle(**pp)); - } - else if((*pp)->end_vertex()) - findDecayProducts(**pp,unstable); + if (id == 331 || id == 9010221) { + unstable.push_back(Particle(*pp)); + } else if ((*pp)->end_vertex()) + findDecayProducts(*pp, unstable); } } + }; + // The hook for the plugin system DECLARE_RIVET_PLUGIN(ARGUS_1993_S2669951); Modified: trunk/src/Analyses/ARGUS_1993_S2789213.cc ============================================================================== --- trunk/src/Analyses/ARGUS_1993_S2789213.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/ARGUS_1993_S2789213.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -27,35 +27,34 @@ const Beam beamproj = applyProjection<Beam>(e, "Beams"); const double s = sqr(beamproj.sqrtS()); const double roots = sqrt(s); - const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); - // find the upsilons + // Find the upsilons ParticleVector upsilons; - // first in unstable final state + // First in unstable final state + const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); foreach (const Particle& p, ufs.particles()) - if(p.pdgId()==300553 || p.pdgId()==553 ) upsilons.push_back(p); - // then in whole event if fails - if(upsilons.empty()) { + if (p.pdgId() == 300553 || p.pdgId() == 553) upsilons.push_back(p); + // Then in whole event if that failed + if (upsilons.empty()) { foreach (GenParticle* p, Rivet::particles(e.genEvent())) { - if( p->pdg_id() != 300553 && p->pdg_id() != 553 ) continue; + if (p->pdg_id() != 300553 && p->pdg_id() != 553) continue; const GenVertex* pv = p->production_vertex(); bool passed = true; if (pv) { - for (GenVertex::particles_in_const_iterator - pp = pv->particles_in_const_begin() ; - pp != pv->particles_in_const_end() ; ++pp) { + /// @todo Use better looping + for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ; pp != pv->particles_in_const_end() ; ++pp) { if ( p->pdg_id() == (*pp)->pdg_id() ) { passed = false; break; } } } - if(passed) upsilons.push_back(Particle(*p)); + if (passed) upsilons.push_back(Particle(*p)); } } // continuum - if(upsilons.empty()) { + if (upsilons.empty()) { _weightSum_cont += weight; unsigned int nOmega(0),nRho0(0),nKStar0(0),nKStarPlus(0),nPhi(0); foreach (const Particle& p, ufs.particles()) { @@ -261,27 +260,27 @@ Histo1DPtr _hist_cont_Omega ; Histo1DPtr _hist_Ups1_Omega ; - // count of weights double _weightSum_cont,_weightSum_Ups1,_weightSum_Ups4; //@} - void findDecayProducts(const GenParticle & p, - ParticleVector & unstable) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + + void findDecayProducts(const GenParticle* p, ParticleVector& unstable) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = abs((*pp)->pdg_id()); - if(id == 113 || id == 313 || id == 323 || - id == 333 || id == 223 ) { - unstable.push_back(Particle(**pp)); + if (id == 113 || id == 313 || id == 323 || + id == 333 || id == 223 ) { + unstable.push_back(Particle(*pp)); } - else if((*pp)->end_vertex()) - findDecayProducts(**pp,unstable); + else if ((*pp)->end_vertex()) + findDecayProducts(*pp, unstable); } } + }; + // The hook for the plugin system DECLARE_RIVET_PLUGIN(ARGUS_1993_S2789213); Modified: trunk/src/Analyses/ATLAS_2011_I9035664.cc ============================================================================== --- trunk/src/Analyses/ATLAS_2011_I9035664.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/ATLAS_2011_I9035664.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -10,7 +10,7 @@ namespace Rivet { - /// @brief J/Psi production at ATLAS + /// @brief J/psi production at ATLAS class ATLAS_2011_I9035664: public Analysis { public: @@ -50,8 +50,8 @@ const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); foreach (const Particle& p, ufs.particles()) { - if (abs(p.pdgId())!=443) continue; - HepMC::GenVertex* gv = p.genParticle().production_vertex(); + if (abs(p.pdgId()) != 443) continue; + HepMC::GenVertex* gv = p.genParticle()->production_vertex(); bool nonPrompt = false; if (gv) { foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { @@ -80,7 +80,7 @@ else if (!nonPrompt) _PrRapMedLow->fill(xp, weight); _IncRapMedLow->fill(xp, weight); } - + else if (rapidity<=0.75) { if (nonPrompt) _nonPrRapLow->fill(xp, weight); else if (!nonPrompt) _PrRapLow->fill(xp, weight); Modified: trunk/src/Analyses/ATLAS_2011_I944826.cc ============================================================================== --- trunk/src/Analyses/ATLAS_2011_I944826.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/ATLAS_2011_I944826.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -6,16 +6,10 @@ #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/UnstableFinalState.hh" -//#include "LWH/Histogram1D.h" -#include "HepMC/GenParticle.h" -#include "HepMC/GenVertex.h" - -//using namespace std; -//using namespace Rivet; - namespace Rivet { + class ATLAS_2011_I944826 : public Analysis { public: @@ -36,7 +30,7 @@ UnstableFinalState ufs(-MAXRAPIDITY, MAXRAPIDITY, 100*MeV); addProjection(ufs, "UFS"); - + vector<pair<double, double> > etaRanges; etaRanges.push_back(make_pair(-3.84, -2.09)); etaRanges.push_back(make_pair(2.09, 3.84)); @@ -65,7 +59,7 @@ //_temp_lambdabar_v_y.reset( new LWH::Histogram1D(10, 0.0, 2.5)); //_temp_lambda_v_pT.reset( new LWH::Histogram1D(18, 0.5, 4.1)); //_temp_lambdabar_v_pT.reset(new LWH::Histogram1D(18, 0.5, 4.1)); - } + } else if (fuzzyEquals(sqrtS()*GeV, 900, 1E-3)) { _hist_Ks_pT = bookHisto1D(4,1,1); _hist_Ks_y = bookHisto1D(5,1,1); @@ -73,7 +67,7 @@ _hist_L_pT = bookHisto1D(10,1,1); _hist_L_y = bookHisto1D(11,1,1); _hist_L_mult = bookHisto1D(12,1,1); - _hist_Ratio_v_y = bookScatter2D(15,1,1); + _hist_Ratio_v_y = bookScatter2D(15,1,1); _hist_Ratio_v_pT = bookScatter2D(16,1,1); /// @todo YODA //_temp_lambda_v_y.reset( new LWH::Histogram1D(5, 0.0, 2.5)); @@ -85,24 +79,24 @@ // This function is required to impose the flight time cuts on Kaons and Lambdas inline double getPerpFlightDistance(const Rivet::Particle& p) { - const HepMC::GenParticle& genp = p.genParticle(); - HepMC::GenVertex* prodV = genp.production_vertex(); - HepMC::GenVertex* decV = genp.end_vertex(); + const HepMC::GenParticle* genp = p.genParticle(); + HepMC::GenVertex* prodV = genp->production_vertex(); + HepMC::GenVertex* decV = genp->end_vertex(); const HepMC::ThreeVector prodPos = prodV->point3d(); if (decV) { - const HepMC::ThreeVector decPos = decV->point3d(); + const HepMC::ThreeVector decPos = decV->point3d(); double dy = prodPos.y() - decPos.y(); double dx = prodPos.x() - decPos.x(); - return sqrt(dx*dx + dy*dy); + return add_quad(dx, dy); } else return 9999999.; } inline bool daughtersSurviveCuts(const Rivet::Particle& p) { - // We require the Kshort or Lambda to decay into two charged + // We require the Kshort or Lambda to decay into two charged // particles with at least 100MeV pT inside acceptance region - const HepMC::GenParticle& genp = p.genParticle(); - HepMC::GenVertex* decV = genp.end_vertex(); + const HepMC::GenParticle* genp = p.genParticle(); + HepMC::GenVertex* decV = genp->end_vertex(); bool decision = true; if (!decV) return false; @@ -157,14 +151,14 @@ // This ufs holds all the Kaons and Lambdas const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(event, "UFS"); - + // Some conters - int n_KS0 = 0; - int n_LAMBDA = 0; + int n_KS0 = 0; + int n_LAMBDA = 0; // Particle loop foreach (const Particle& p, ufs.particles()) { - + // General particle quantities const double pT = p.momentum().pT()*GeV; const double y = p.momentum().rapidity(); @@ -174,7 +168,7 @@ // Look for Kaons, Lambdas switch (pid) { - + case K0S: flightd = getPerpFlightDistance(p); if (!inRange(flightd, 4., 450.) ){ @@ -184,7 +178,7 @@ if (daughtersSurviveCuts(p) ) { _hist_Ks_y ->fill(y, weight); _hist_Ks_pT->fill(pT, weight); - + _sum_w_ks += weight; n_KS0++; } @@ -193,7 +187,7 @@ case LAMBDA: if (pT < 0.5) { // Lambdas have an additional pT cut of 500 MeV MSG_DEBUG("Lambda failed pT cut:" << pT); - break; + break; } flightd = getPerpFlightDistance(p); if (!inRange(flightd, 17., 450.)) { @@ -205,14 +199,14 @@ /// @todo YODA //_temp_lambda_v_y ->fill(fabs(y), weight); //_temp_lambda_v_pT ->fill(pT, weight); - + _hist_L_y->fill( y, weight); _hist_L_pT->fill(pT, weight); - + _sum_w_lambda += weight; n_LAMBDA++; } - + else if (p.pdgId() == -3122) { /// @todo YODA //_temp_lambdabar_v_y ->fill(fabs(y), weight); @@ -220,16 +214,16 @@ } } break; - + } // End of switch - + }// End of particle loop // Fill multiplicity histos _hist_Ks_mult->fill(n_KS0, weight); _hist_L_mult->fill(n_LAMBDA, weight); } - + /// Normalise histograms etc., after the run @@ -245,15 +239,15 @@ scale(_hist_L_pT, 1.0/_sum_w_lambda); scale(_hist_L_y, 1.0/_sum_w_lambda); scale(_hist_L_mult, 1.0/_sum_w_passed); - + /// @todo YODA //// Division of histograms to obtain lambdabar/lambda ratios //if (fuzzyEquals(sqrtS()*GeV, 7000, 1E-3)) { // histogramFactory().divide(histoPath("d13-x01-y01"), *_temp_lambdabar_v_y, *_temp_lambda_v_y ); // histogramFactory().divide(histoPath("d14-x01-y01"), *_temp_lambdabar_v_pT, *_temp_lambda_v_pT); - //} - //else if (fuzzyEquals(sqrtS()*GeV, 900, 1E-3)) { + //} + //else if (fuzzyEquals(sqrtS()*GeV, 900, 1E-3)) { // histogramFactory().divide(histoPath("d15-x01-y01"), *_temp_lambdabar_v_y, *_temp_lambda_v_y ); // histogramFactory().divide(histoPath("d16-x01-y01"), *_temp_lambdabar_v_pT, *_temp_lambda_v_pT); //} @@ -279,7 +273,7 @@ Scatter2DPtr _hist_Ratio_v_pT; Scatter2DPtr _hist_Ratio_v_y; - + /// @todo YODA //shared_ptr<LWH::Histogram1D> _temp_lambda_v_y, _temp_lambdabar_v_y; //shared_ptr<LWH::Histogram1D> _temp_lambda_v_pT, _temp_lambdabar_v_pT; Modified: trunk/src/Analyses/BABAR_2003_I593379.cc ============================================================================== --- trunk/src/Analyses/BABAR_2003_I593379.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/BABAR_2003_I593379.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -9,7 +9,7 @@ namespace Rivet { - /// @brief Babar chamronium spectra + /// @brief Babar charmonium spectra /// @author Peter Richardson class BABAR_2003_I593379 : public Analysis { public: @@ -21,21 +21,21 @@ void analyze(const Event& e) { const double weight = e.weight(); - const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); - // find the upsilons + + // Find the charmonia ParticleVector upsilons; - // first in unstable final state + // First in unstable final state + const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); foreach (const Particle& p, ufs.particles()) if (p.pdgId()==300553) upsilons.push_back(p); - // then in whole event if fails + // Then in whole event if fails if (upsilons.empty()) { foreach (GenParticle* p, Rivet::particles(e.genEvent())) { if (p->pdg_id()!=300553) continue; const GenVertex* pv = p->production_vertex(); bool passed = true; if (pv) { - for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ; - pp != pv->particles_in_const_end() ; ++pp) { + for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin(); pp != pv->particles_in_const_end(); ++pp) { if ( p->pdg_id() == (*pp)->pdg_id() ) { passed = false; break; @@ -55,12 +55,12 @@ findDecayProducts(p.genParticle(), allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2); LorentzTransform cms_boost(-p.momentum().boostVector()); - for (size_t i=0; i<allJpsi.size(); i++) { + for (size_t i = 0; i < allJpsi.size(); i++) { double pcm = cms_boost.transform(FourMomentum(allJpsi[i]->momentum())).vector3().mod(); _hist_all_Jpsi->fill(pcm, weight); } _mult_JPsi->fill(10.58, weight*double(allJpsi.size())); - for (size_t i=0; i<primaryJpsi.size(); i++) { + for (size_t i = 0; i < primaryJpsi.size(); i++) { double pcm = cms_boost.transform(FourMomentum(primaryJpsi[i]->momentum())).vector3().mod(); _hist_primary_Jpsi->fill(pcm, weight); } @@ -70,13 +70,13 @@ _hist_Psi_prime->fill(pcm, weight); } _mult_Psi2S->fill(10.58, weight*double(Psiprime.size())); - for (size_t i=0; i<all_chi_c1.size(); i++) { + for (size_t i = 0; i < all_chi_c1.size(); i++) { double pcm = cms_boost.transform(FourMomentum(all_chi_c1[i]->momentum())).vector3().mod(); _hist_chi_c1->fill(pcm, weight); } _mult_chi_c1->fill(10.58, weight*double(all_chi_c1.size())); _mult_chi_c1_direct->fill(10.58, weight*double(primary_chi_c1.size())); - for (size_t i=0; i<all_chi_c2.size(); i++) { + for (size_t i = 0; i < all_chi_c2.size(); i++) { double pcm = cms_boost.transform(FourMomentum(all_chi_c2[i]->momentum())).vector3().mod(); _hist_chi_c2->fill(pcm, weight); } @@ -85,8 +85,8 @@ } } // analyze - void finalize() { + void finalize() { scale(_hist_all_Jpsi , 0.5*0.1/_weightSum); scale(_hist_chi_c1 , 0.5*0.1/_weightSum); scale(_hist_chi_c2 , 0.5*0.1/_weightSum); @@ -140,27 +140,24 @@ Histo1DPtr _mult_Psi2S; //@} - void findDecayProducts(const GenParticle & p, - vector<GenParticle *> & allJpsi, - vector<GenParticle *> & primaryJpsi, - vector<GenParticle *> & Psiprime, - vector<GenParticle *> & all_chi_c1, - vector<GenParticle *> & all_chi_c2, - vector<GenParticle *> & primary_chi_c1, - vector<GenParticle *> & primary_chi_c2) { - const GenVertex* dv = p.end_vertex(); - bool isOnium(false); - for (GenVertex::particles_in_const_iterator pp = dv->particles_in_const_begin() ; - pp != dv->particles_in_const_end() ; ++pp) { + void findDecayProducts(const GenParticle* p, + vector<GenParticle*>& allJpsi, + vector<GenParticle*>& primaryJpsi, + vector<GenParticle*>& Psiprime, + vector<GenParticle*>& all_chi_c1, vector<GenParticle*>& all_chi_c2, + vector<GenParticle*>& primary_chi_c1, vector<GenParticle*>& primary_chi_c2) { + const GenVertex* dv = p->end_vertex(); + bool isOnium = false; + /// @todo Use better looping + for (GenVertex::particles_in_const_iterator pp = dv->particles_in_const_begin() ; pp != dv->particles_in_const_end() ; ++pp) { int id = (*pp)->pdg_id(); id = id%1000; id -= id%10; id /= 10; if (id==44) isOnium = true; } - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); if (id==100443) { Psiprime.push_back(*pp); @@ -178,11 +175,11 @@ if (!isOnium) primaryJpsi.push_back(*pp); } if ((*pp)->end_vertex()) { - findDecayProducts(**pp, allJpsi, primaryJpsi, Psiprime, - all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2); + findDecayProducts(*pp, allJpsi, primaryJpsi, Psiprime, all_chi_c1, all_chi_c2, primary_chi_c1, primary_chi_c2); } } } + }; Modified: trunk/src/Analyses/BABAR_2005_S6181155.cc ============================================================================== --- trunk/src/Analyses/BABAR_2005_S6181155.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/BABAR_2005_S6181155.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -98,52 +98,50 @@ Histo1DPtr _histOffResonance_norm; //@} - bool checkDecay(const GenParticle & p) { - unsigned int nstable=0,npip=0,npim=0; - unsigned int nXim=0,nXip=0; - findDecayProducts(p,nstable,npip,npim, - nXip,nXim); - int id = p.pdg_id(); + bool checkDecay(const GenParticle* p) { + unsigned int nstable = 0, npip = 0, npim = 0; + unsigned int nXim = 0, nXip = 0; + findDecayProducts(p, nstable, npip, npim, nXip, nXim); + int id = p->pdg_id(); // Xi_c - if(id==4132) { - if(nstable==2&&nXim==1&&npip==1) return true; + if (id==4132) { + if (nstable==2 && nXim==1 && npip==1) return true; } - else if(id==-4132) { - if(nstable==2&&nXip==1&&npim==1) return true; + else if (id==-4132) { + if (nstable==2 && nXip==1 && npim==1) return true; } return false; } - void findDecayProducts(const GenParticle & p, - unsigned int & nstable, - unsigned int & npip, unsigned int & npim, - unsigned int & nXip, unsigned int & nXim) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + void findDecayProducts(const GenParticle* p, + unsigned int& nstable, + unsigned int& npip, unsigned int& npim, + unsigned int& nXip, unsigned int& nXim) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); - if(id==3312) { + if (id==3312) { ++nXim; ++nstable; - } - else if(id==-3312) { + } else if (id==-3312) { ++nXip; ++nstable; - } - else if(id==111||id==221) + } else if(id==111||id==221) { ++nstable; - else if((*pp)->end_vertex()) - findDecayProducts(**pp,nstable,npip,npim,nXip,nXim); - else { - if(id!=22) ++nstable; + } else if ((*pp)->end_vertex()) { + findDecayProducts(*pp, nstable, npip, npim, nXip, nXim); + } else { + if (id != 22) ++nstable; if (id == 211) ++npip; else if(id == -211) ++npim; } } } + }; + // The hook for the plugin system DECLARE_RIVET_PLUGIN(BABAR_2005_S6181155); Modified: trunk/src/Analyses/BABAR_2007_S7266081.cc ============================================================================== --- trunk/src/Analyses/BABAR_2007_S7266081.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/BABAR_2007_S7266081.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -22,10 +22,9 @@ void analyze(const Event& e) { - const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); - - // find the taus + // Find the taus ParticleVector taus; + const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); foreach (const Particle& p, ufs.particles()) { if(abs(p.pdgId())!=15) continue; _weight_total += 1.; @@ -36,14 +35,14 @@ if(p.momentum().vector3().mod()>0.001) cms_boost = LorentzTransform(-p.momentum().boostVector()); // find the decay products we want - findDecayProducts(p.genParticle(),nstable,pip,pim,Kp,Km); - if(p.pdgId()<0) { + findDecayProducts(p.genParticle(), nstable, pip, pim, Kp, Km); + if (p.pdgId()<0) { swap(pip,pim); swap(Kp ,Km ); } - if(nstable!=4) continue; + if (nstable!=4) continue; // pipipi - if(pim.size()==2&&pip.size()==1) { + if (pim.size()==2&&pip.size()==1) { _weight_pipippi += 1.; _hist_pipipi_pipipi-> fill((pip[0].momentum()+pim[0].momentum()+pim[1].momentum()).mass(),1.); @@ -52,7 +51,7 @@ _hist_pipipi_pipi-> fill((pip[0].momentum()+pim[1].momentum()).mass(),1.); } - else if(pim.size()==1&&pip.size()==1&&Km.size()==1) { + else if (pim.size()==1&&pip.size()==1&&Km.size()==1) { _weight_Kpipi += 1.; _hist_Kpipi_Kpipi-> fill((pim[0].momentum()+pip[0].momentum()+Km[0].momentum()).mass(),1.); @@ -61,7 +60,7 @@ _hist_Kpipi_pipi-> fill((pim[0].momentum()+pip[0].momentum()).mass(),1.); } - else if(Kp.size()==1&&Km.size()==1&&pim.size()==1) { + else if (Kp.size()==1&&Km.size()==1&&pim.size()==1) { _weight_KpiK += 1.; _hist_KpiK_KpiK-> fill((Kp[0].momentum()+Km[0].momentum()+pim[0].momentum()).mass(),1.); @@ -70,7 +69,7 @@ _hist_KpiK_piK-> fill((Kp[0].momentum()+pim[0].momentum()).mass(),1.); } - else if(Kp.size()==1&&Km.size()==2) { + else if (Kp.size()==1&&Km.size()==2) { _weight_KKK += 1.; _hist_KKK_KKK-> fill((Kp[0].momentum()+Km[0].momentum()+Km[1].momentum()).mass(),1.); @@ -82,6 +81,7 @@ } } // analyze + void finalize() { if(_weight_pipippi>0.) { scale(_hist_pipipi_pipipi, 1./_weight_pipippi); @@ -120,6 +120,7 @@ //br_KKK->point(0)->coordinate(1)->setErrorMinus( 100.*sqrt(_weight_KKK)/_weight_total); } // finalize + void init() { addProjection(UnstableFinalState(), "UFS"); @@ -136,6 +137,7 @@ } // init + private: //@{ @@ -154,13 +156,13 @@ double _weight_total,_weight_pipippi,_weight_Kpipi,_weight_KpiK,_weight_KKK; //@} - void findDecayProducts(const GenParticle & p, unsigned int & nstable, - ParticleVector & pip, ParticleVector & pim, - ParticleVector & Kp, ParticleVector & Km) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + void findDecayProducts(const GenParticle* p, + unsigned int & nstable, + ParticleVector& pip, ParticleVector& pim, + ParticleVector& Kp, ParticleVector& Km) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); if( id == PI0 ) ++nstable; @@ -183,14 +185,17 @@ ++nstable; } else if((*pp)->end_vertex()) { - findDecayProducts(**pp,nstable,pip,pim,Kp,Km); + findDecayProducts(*pp, nstable, pip, pim, Kp, Km); } else ++nstable; } } + + }; + // The hook for the plugin system DECLARE_RIVET_PLUGIN(BABAR_2007_S7266081); Modified: trunk/src/Analyses/BELLE_2001_S4598261.cc ============================================================================== --- trunk/src/Analyses/BELLE_2001_S4598261.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/BELLE_2001_S4598261.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -22,29 +22,29 @@ void analyze(const Event& e) { const double weight = e.weight(); - const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); // find the upsilons ParticleVector upsilons; // first in unstable final state + const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS"); foreach (const Particle& p, ufs.particles()) - if(p.pdgId()==300553) upsilons.push_back(p); + if (p.pdgId()==300553) upsilons.push_back(p); // then in whole event if fails - if(upsilons.empty()) { - foreach (GenParticle* p, Rivet::particles(e.genEvent())) { - if(p->pdg_id()!=300553) continue; + if (upsilons.empty()) { + foreach (const GenParticle* p, Rivet::particles(e.genEvent())) { + if (p->pdg_id() != 300553) continue; const GenVertex* pv = p->production_vertex(); bool passed = true; if (pv) { - for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ; - pp != pv->particles_in_const_end() ; ++pp) { + /// @todo Use better looping + for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ; pp != pv->particles_in_const_end() ; ++pp) { if ( p->pdg_id() == (*pp)->pdg_id() ) { passed = false; break; } } } - if(passed) upsilons.push_back(Particle(*p)); + if (passed) upsilons.push_back(Particle(p)); } } @@ -53,7 +53,7 @@ _weightSum += weight; // find the neutral pions from the decay vector<GenParticle *> pions; - findDecayProducts(p.genParticle(),pions); + findDecayProducts(p.genParticle(), pions); LorentzTransform cms_boost(-p.momentum().boostVector()); for(unsigned int ix=0;ix<pions.size();++ix) { double pcm = @@ -90,20 +90,19 @@ Histo1DPtr _histMult; //@} - void findDecayProducts(const GenParticle & p, - vector<GenParticle *> & pions) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + void findDecayProducts(const GenParticle* p, vector<GenParticle*>& pions) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); - if(id==111) { + if (id==111) { pions.push_back(*pp); - } - else if((*pp)->end_vertex()) - findDecayProducts(**pp,pions); + } else if((*pp)->end_vertex()) + findDecayProducts(*pp, pions); } } + + }; Modified: trunk/src/Analyses/BELLE_2006_S6265367.cc ============================================================================== --- trunk/src/Analyses/BELLE_2006_S6265367.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/BELLE_2006_S6265367.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -122,8 +122,8 @@ _sigmaDStarPlusB->fill(10.6,weight); _sigmaDStarPlusC->fill(10.6,weight); } - const GenParticle * Dmeson = &p.genParticle(); - const GenVertex* dv = p.genParticle().end_vertex(); + const GenParticle* Dmeson = p.genParticle(); + const GenVertex* dv = p.genParticle()->end_vertex(); bool D0decay(false), Pi0decay(false), Piplusdecay(false), Dplusdecay(false); for (GenVertex::particles_out_const_iterator @@ -141,13 +141,13 @@ Piplusdecay = true; } } - if (D0decay && Piplusdecay && checkDecay(*Dmeson)) { + if (D0decay && Piplusdecay && checkDecay(Dmeson)) { if (onresonance) _histXpDstarplus2D0_R->fill(xp, s*weight); else _histXpDstarplus2D0_C->fill(xp, s*weight); } - else if (Dplusdecay && Pi0decay && checkDecay(*Dmeson)) { + else if (Dplusdecay && Pi0decay && checkDecay(Dmeson)) { if (onresonance) _histXpDstarplus2Dplus_R->fill(xp, s*weight); else @@ -173,12 +173,11 @@ xp = mom/sqrt(s/4.0 - mH2); if(!onresonance) _sigmaDStar0 ->fill(10.6,weight); MSG_DEBUG("xp = " << xp); - const GenParticle * Dmeson = &p.genParticle(); - const GenVertex* dv = p.genParticle().end_vertex(); + const GenParticle* Dmeson = p.genParticle(); + const GenVertex* dv = p.genParticle()->end_vertex(); bool D0decay(false), Pi0decay(false); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { if (abs((*pp)->pdg_id()) == 421) { Dmeson = *pp; D0decay = true; @@ -187,7 +186,7 @@ Pi0decay = true; } } - if (D0decay && Pi0decay && checkDecay(*Dmeson)) { + if (D0decay && Pi0decay && checkDecay(Dmeson)) { if (onresonance) _histXpDstar0_R->fill(xp, s*weight); else { @@ -349,12 +348,11 @@ Histo1DPtr _histXpDstar0_R_N; //@} - bool checkDecay(const GenParticle & p) { + bool checkDecay(const GenParticle* p) { unsigned int nstable=0,npip=0,npim=0; unsigned int np=0,nap=0,nKp=0,nKm=0,nPhi=0; - findDecayProducts(p,nstable,npip,npim, - np,nap,nKp,nKm,nPhi); - int id = p.pdg_id(); + findDecayProducts(p, nstable, npip, npim, np, nap, nKp, nKm, nPhi); + int id = p->pdg_id(); //D0 if(id==421) { if(nstable==2&&nKm==1&&npip==1) return true; @@ -390,24 +388,23 @@ return false; } - void findDecayProducts(const GenParticle & p, + void findDecayProducts(const GenParticle* p, unsigned int & nstable, unsigned int & npip, unsigned int & npim , unsigned int & np, unsigned int & nap , unsigned int & nKp, unsigned int & nKm , unsigned int & nPhi) { - const GenVertex* dv = p.end_vertex(); - for (GenVertex::particles_out_const_iterator - pp = dv->particles_out_const_begin(); - pp != dv->particles_out_const_end(); ++pp) { + const GenVertex* dv = p->end_vertex(); + /// @todo Use better looping + for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { int id = (*pp)->pdg_id(); if(id==333) ++nPhi; else if(id==111||id==221) ++nstable; else if((*pp)->end_vertex()) - findDecayProducts(**pp,nstable,npip,npim,np,nap,nKp,nKm,nPhi); + findDecayProducts(*pp, nstable, npip, npim, np, nap, nKp, nKm, nPhi); else { - if(id!=22) ++nstable; + if (id != 22) ++nstable; if (id == 211) ++npip; else if(id == -211) ++npim; else if(id == 2212) ++np; @@ -417,6 +414,7 @@ } } } + }; Modified: trunk/src/Analyses/CDF_2006_S6653332.cc ============================================================================== --- trunk/src/Analyses/CDF_2006_S6653332.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/CDF_2006_S6653332.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -91,8 +91,7 @@ /// @todo Use jet contents rather than accessing quarks directly ParticleVector bquarks; /// @todo Use nicer looping - for (GenEvent::particle_const_iterator p = event.genEvent().particles_begin(); - p != event.genEvent().particles_end(); ++p) { + for (GenEvent::particle_const_iterator p = event.genEvent()->particles_begin(); p != event.genEvent()->particles_end(); ++p) { if ( fabs((*p)->pdg_id()) == BQUARK ) { bquarks.push_back(Particle(**p)); } Modified: trunk/src/Analyses/CDF_2008_S7540469.cc ============================================================================== --- trunk/src/Analyses/CDF_2008_S7540469.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/CDF_2008_S7540469.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -49,8 +49,7 @@ // Skip if the event is empty const FinalState& fs = applyProjection<FinalState>(event, "FS"); if (fs.empty()) { - MSG_DEBUG("Skipping event " << event.genEvent().event_number() - << " because no final state pair found "); + MSG_DEBUG("Skipping event " << numEvents() << " because no final state pair found"); vetoEvent; } @@ -61,22 +60,22 @@ for (size_t i=0; i<all_els.size(); ++i) { for (size_t j=i+1; j<all_els.size(); ++j) { bool candidate=true; - double mZ=FourMomentum(all_els[i].momentum()+all_els[j].momentum()).mass()/GeV; - if (mZ<66.0 || mZ>116.0) { - candidate=false; - } - double abs_eta_0=fabs(all_els[i].momentum().pseudorapidity()); - double abs_eta_1=fabs(all_els[j].momentum().pseudorapidity()); - if (abs_eta_1<abs_eta_0) { - double tmp=abs_eta_0; - abs_eta_0=abs_eta_1; - abs_eta_1=tmp; + double mZ = FourMomentum(all_els[i].momentum()+all_els[j].momentum()).mass()/GeV; + if (mZ < 66.0 || mZ > 116.0) { + candidate = false; + } + double abs_eta_0 = fabs(all_els[i].momentum().pseudorapidity()); + double abs_eta_1 = fabs(all_els[j].momentum().pseudorapidity()); + if (abs_eta_1 < abs_eta_0) { + double tmp = abs_eta_0; + abs_eta_0 = abs_eta_1; + abs_eta_1 = tmp; } - if (abs_eta_0>1.0) { - candidate=false; + if (abs_eta_0 > 1.0) { + candidate = false; } - if (!(abs_eta_1<1.0 || (abs_eta_1>1.2 && abs_eta_1<2.8))) { - candidate=false; + if (!(abs_eta_1 < 1.0 || (inRange(abs_eta_1, 1.2, 2.8)))) { + candidate = false; } if (candidate) { Z_candidates.push_back(make_pair(all_els[i], all_els[j])); @@ -84,8 +83,7 @@ } } if (Z_candidates.size() != 1) { - MSG_DEBUG("Skipping event " << event.genEvent().event_number() - << " because no unique electron pair found "); + MSG_DEBUG("Skipping event " << numEvents() << " because no unique electron pair found "); vetoEvent; } @@ -107,10 +105,10 @@ copy = false; } } else { - if (p.genParticle().barcode()==Z_candidates[0].first.genParticle().barcode()) { + if (p.genParticle()->barcode() == Z_candidates[0].first.genParticle()->barcode()) { copy = false; } - if (p.genParticle().barcode()==Z_candidates[0].second.genParticle().barcode()) { + if (p.genParticle()->barcode() == Z_candidates[0].second.genParticle()->barcode()) { copy = false; } } Modified: trunk/src/Analyses/D0_2008_S7863608.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7863608.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/D0_2008_S7863608.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -68,8 +68,7 @@ // Return if there are no jets: if(jets_cut.size()<1) { - MSG_DEBUG("Skipping event " << e.genEvent().event_number() - << " because no jets pass cuts "); + MSG_DEBUG("Skipping event " << numEvents() << " because no jets pass cuts "); vetoEvent; } Modified: trunk/src/Analyses/D0_2009_S8349509.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8349509.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/D0_2009_S8349509.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -17,8 +17,9 @@ //@{ /// Constructor - D0_2009_S8349509() : Analysis("D0_2009_S8349509"), - _inclusive_Z_sumofweights(0.0) + D0_2009_S8349509() + : Analysis("D0_2009_S8349509"), + _inclusive_Z_sumofweights(0.0) { } //@} @@ -81,8 +82,7 @@ // Return if there are no jets: if (jets.size() < 1) { - MSG_DEBUG("Skipping event " << event.genEvent().event_number() - << " because no jets pass cuts "); + MSG_DEBUG("Skipping event " << numEvents() << " because no jets pass cuts "); vetoEvent; } Modified: trunk/src/Analyses/D0_2010_S8570965.cc ============================================================================== --- trunk/src/Analyses/D0_2010_S8570965.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/D0_2010_S8570965.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -68,12 +68,12 @@ double phi_P = photon.momentum().phi(); double Etsum=0.0; foreach (const Particle& p, fs) { - if (p.genParticle().barcode()!=photon.genParticle().barcode() && + if (p.genParticle()->barcode() != photon.genParticle()->barcode() && deltaR(eta_P, phi_P, p.momentum().eta(), p.momentum().phi()) < 0.4) { Etsum += p.momentum().Et(); } } - if (Etsum<2.5*GeV) { + if (Etsum < 2.5*GeV) { isolated_photons.push_back(photon); } } Modified: trunk/src/Analyses/H1_1994_S2919893.cc ============================================================================== --- trunk/src/Analyses/H1_1994_S2919893.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/H1_1994_S2919893.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -52,10 +52,10 @@ // Extract the particles other than the lepton ParticleVector particles; particles.reserve(fs.particles().size()); - const GenParticle& dislepGP = dl.out().genParticle(); + const GenParticle* dislepGP = dl.out().genParticle(); foreach (const Particle& p, fs.particles()) { - const GenParticle& loopGP = p.genParticle(); - if (&loopGP == &dislepGP) continue; + const GenParticle* loopGP = p.genParticle(); + if (loopGP == dislepGP) continue; particles.push_back(p); } Modified: trunk/src/Analyses/H1_2000_S4129130.cc ============================================================================== --- trunk/src/Analyses/H1_2000_S4129130.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/H1_2000_S4129130.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -45,12 +45,11 @@ // Extract the particles other than the lepton ParticleVector particles; particles.reserve(fs.particles().size()); - const GenParticle& dislepGP = dl.out().genParticle(); - for (ParticleVector::const_iterator p = fs.particles().begin(); - p != fs.particles().end(); ++p) { - const GenParticle& loopGP = p->genParticle(); - if (&loopGP == &dislepGP) continue; - particles.push_back(*p); + const GenParticle* dislepGP = dl.out().genParticle(); + foreach (const Particle& p, fs.particles()) { + const GenParticle* loopGP = p.genParticle(); + if (loopGP == dislepGP) continue; + particles.push_back(p); } // Cut on the forward energy Modified: trunk/src/Analyses/LHCB_2010_S8758301.cc ============================================================================== --- trunk/src/Analyses/LHCB_2010_S8758301.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/LHCB_2010_S8758301.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -171,8 +171,8 @@ const GenParticle* getLongestLivedAncestor(const Particle& p, double& lifeTime) { const GenParticle* ret = NULL; lifeTime = 1.; - if (!p.hasGenParticle()) return NULL; - const GenParticle* pmother = &(p.genParticle()); + if (p.genParticle() == NULL) return NULL; + const GenParticle* pmother = p.genParticle(); double longest_ctau = 0.; double mother_ctau; int mother_pid, n_inparts; Modified: trunk/src/Analyses/LHCB_2011_I917009.cc ============================================================================== --- trunk/src/Analyses/LHCB_2011_I917009.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/LHCB_2011_I917009.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -168,10 +168,10 @@ // Data members like post-cuts event weight counters go here const double getMotherLifeTimeSum(const Particle& p) { - if ( !p.hasGenParticle() ) return -1.; + if (p.genParticle() == NULL) return -1.; double lftSum = 0.; double plft = 0.; - const GenParticle* part = &(p.genParticle()); + const GenParticle* part = p.genParticle(); GenVertex* ivtx = part->production_vertex(); while(ivtx) { Modified: trunk/src/Analyses/LHCB_2012_I1119400.cc ============================================================================== --- trunk/src/Analyses/LHCB_2012_I1119400.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/LHCB_2012_I1119400.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -200,10 +200,10 @@ // Data members like post-cuts event weight counters go here const double getMotherLifeTimeSum(const Particle& p) { - if ( !p.hasGenParticle() ) return -1.; + if (p.genParticle() == NULL) return -1.; double lftSum = 0.; double plft = 0.; - const GenParticle* part = &(p.genParticle()); + const GenParticle* part = p.genParticle(); GenVertex* ivtx = part->production_vertex(); while(ivtx) { Modified: trunk/src/Analyses/MC_PDFS.cc ============================================================================== --- trunk/src/Analyses/MC_PDFS.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/MC_PDFS.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -43,8 +43,8 @@ const double weight = event.weight(); // This analysis needs a valid HepMC PDF info object to do anything - if (event.genEvent().pdf_info() == 0) vetoEvent; - HepMC::PdfInfo pdfi = *event.genEvent().pdf_info(); + if (event.genEvent()->pdf_info() == 0) vetoEvent; + HepMC::PdfInfo pdfi = *(event.genEvent()->pdf_info()); MSG_DEBUG("PDF Q = " << pdfi.scalePDF() << " for (id, x) = " << "(" << pdfi.id1() << ", " << pdfi.x1() << ") " Modified: trunk/src/Analyses/MC_PRINTEVENT.cc ============================================================================== --- trunk/src/Analyses/MC_PRINTEVENT.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/MC_PRINTEVENT.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -153,30 +153,31 @@ /// Perform the per-event analysis void analyze(const Event& event) { - //printEvent(event.genEvent()); + /// @todo Wouldn't this be nice... if HepMC::IO_AsciiParticles was sane :-/ + // printEvent(event.genEvent()); - const GenEvent& evt = event.genEvent(); + const GenEvent* evt = event.genEvent(); cout << string(120, '=') << "\n" << endl; // Weights - cout << "Weights(" << evt.weights().size() << ")="; - for (HepMC::WeightContainer::const_iterator wgt = evt.weights().begin(); wgt != evt.weights().end(); ++wgt) { + cout << "Weights(" << evt->weights().size() << ")="; + for (HepMC::WeightContainer::const_iterator wgt = evt->weights().begin(); wgt != evt->weights().end(); ++wgt) { cout << *wgt << " "; } cout << "\n" - << "EventScale " << evt.event_scale() - << " [energy] \t alphaQCD=" << evt.alphaQCD() - << "\t alphaQED=" << evt.alphaQED() << endl; - - if (evt.pdf_info()) { - cout << "PdfInfo: id1=" << evt.pdf_info()->id1() - << " id2=" << evt.pdf_info()->id2() - << " x1=" << evt.pdf_info()->x1() - << " x2=" << evt.pdf_info()->x2() - << " q=" << evt.pdf_info()->scalePDF() - << " xpdf1=" << evt.pdf_info()->pdf1() - << " xpdf2=" << evt.pdf_info()->pdf2() + << "EventScale " << evt->event_scale() + << " [energy] \t alphaQCD=" << evt->alphaQCD() + << "\t alphaQED=" << evt->alphaQED() << endl; + + if (evt->pdf_info()) { + cout << "PdfInfo: id1=" << evt->pdf_info()->id1() + << " id2=" << evt->pdf_info()->id2() + << " x1=" << evt->pdf_info()->x1() + << " x2=" << evt->pdf_info()->x2() + << " q=" << evt->pdf_info()->scalePDF() + << " xpdf1=" << evt->pdf_info()->pdf1() + << " xpdf2=" << evt->pdf_info()->pdf2() << endl; } else { cout << "PdfInfo: EMPTY"; @@ -196,7 +197,7 @@ // Print all particles // const HepPDT::ParticleDataTable* pdt = m_ppsvc->PDT(); - for (HepMC::GenEvent::particle_const_iterator p = evt.particles_begin(); p != evt.particles_end(); ++p) { + for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) { int p_bcode = (*p)->barcode(); int p_pdg_id = (*p)->pdg_id(); double p_px = (*p)->momentum().px(); Modified: trunk/src/Analyses/MC_XS.cc ============================================================================== --- trunk/src/Analyses/MC_XS.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/MC_XS.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -40,24 +40,23 @@ /// Perform the per-event analysis void analyze(const Event& event) { _h_N->fill(0.5,1.); - _h_pmXS->fill(0.5*(event.weight()>0?1.:-1),abs(event.weight())); - _h_pmN->fill(0.5*(event.weight()>0?1.:-1),1.); -#ifdef HEPMC_HAS_CROSS_SECTION - _mc_xs = event.genEvent().cross_section()->cross_section(); - _mc_error = event.genEvent().cross_section()->cross_section_error(); -#endif + _h_pmXS->fill(0.5*(event.weight()> 0 ? 1. : -1), abs(event.weight())); + _h_pmN->fill(0.5*(event.weight() > 0 ? 1. : -1), 1.); + #ifdef HEPMC_HAS_CROSS_SECTION + _mc_xs = event.genEvent()->cross_section()->cross_section(); + _mc_error = event.genEvent()->cross_section()->cross_section_error(); + #endif } /// Normalise histograms etc., after the run void finalize() { scale(_h_pmXS, crossSection()/sumOfWeights()); -#ifndef HEPMC_HAS_CROSS_SECTION - _mc_xs=crossSection(); - _mc_error=0.0; -#endif - _h_XS->addPoint(0,_mc_xs, - 0,_mc_error); + #ifndef HEPMC_HAS_CROSS_SECTION + _mc_xs = crossSection(); + _mc_error = 0.0; + #endif + _h_XS->addPoint(0, _mc_xs, 0, _mc_error); } //@} @@ -74,11 +73,9 @@ double _mc_xs, _mc_error; //@} - }; - // The hook for the plugin system DECLARE_RIVET_PLUGIN(MC_XS); Modified: trunk/src/Analyses/STAR_2006_S6870392.cc ============================================================================== --- trunk/src/Analyses/STAR_2006_S6870392.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/STAR_2006_S6870392.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -1,9 +1,9 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" +#include "Rivet/RivetYODA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" -#include "Rivet/RivetYODA.hh" namespace Rivet { @@ -40,8 +40,7 @@ // Skip if the event is empty const FinalState& fs = applyProjection<FinalState>(event, "FS"); if (fs.empty()) { - MSG_DEBUG("Skipping event " << event.genEvent().event_number() - << " because no final state found "); + MSG_DEBUG("Skipping event " << numEvents() << " because no final state found "); vetoEvent; } Modified: trunk/src/Analyses/STAR_2008_S7993412.cc ============================================================================== --- trunk/src/Analyses/STAR_2008_S7993412.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Analyses/STAR_2008_S7993412.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -34,8 +34,7 @@ // Skip if the event is empty const FinalState& fs = applyProjection<FinalState>(event, "FS"); if (fs.empty()) { - MSG_DEBUG("Skipping event " << event.genEvent().event_number() - << " because no final state found "); + MSG_DEBUG("Skipping event " << numEvents() << " because no final state found "); vetoEvent; } Modified: trunk/src/Core/Event.cc ============================================================================== --- trunk/src/Core/Event.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Core/Event.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -7,37 +7,57 @@ namespace Rivet { - void _geRot180x(GenEvent& ge) { - for (HepMC::GenEvent::particle_iterator ip = ge.particles_begin(); ip != ge.particles_end(); ++ip) { - const HepMC::FourVector& mom = (*ip)->momentum(); - (*ip)->set_momentum(HepMC::FourVector(mom.px(), -mom.py(), -mom.pz(), mom.e())); - } - for (HepMC::GenEvent::vertex_iterator iv = ge.vertices_begin(); iv != ge.vertices_end(); ++iv) { - const HepMC::FourVector& pos = (*iv)->position(); - (*iv)->set_position(HepMC::FourVector(pos.x(), -pos.y(), -pos.z(), pos.t())); + void Event::_init(const GenEvent& ge) { + // Use Rivet's preferred units if possible + #ifdef HEPMC_HAS_UNITS + _genEvent.use_units(HepMC::Units::GEV, HepMC::Units::MM); + #endif + + // Use the conventional alignment + _geNormAlignment(); + + // @todo Filter the event to remove generator-specific particles: optional + // behaviour? Maybe disableable in an inconvenient way, e.g. with an env + // var, to communicate the appropriate distaste for this sort of truth + // analysis ;-) + + // Debug printout to check that copying/mangling has worked + /// @todo Enable this when HepMC has been fixed to allow printing to a stream like the Rivet logger. + //_genEvent.print(); + } + + + namespace { // unnamed namespace for hiding + + void _geRot180x(GenEvent& ge) { + /// @todo Use nicer iterators over HepMC particles + for (HepMC::GenEvent::particle_iterator ip = ge.particles_begin(); ip != ge.particles_end(); ++ip) { + const HepMC::FourVector& mom = (*ip)->momentum(); + (*ip)->set_momentum(HepMC::FourVector(mom.px(), -mom.py(), -mom.pz(), mom.e())); + } + /// @todo Use nicer iterators over HepMC vertices + for (HepMC::GenEvent::vertex_iterator iv = ge.vertices_begin(); iv != ge.vertices_end(); ++iv) { + const HepMC::FourVector& pos = (*iv)->position(); + (*iv)->set_position(HepMC::FourVector(pos.x(), -pos.y(), -pos.z(), pos.t())); + } } + } - // Convert the GenEvent to use conventional alignment - // (proton or electron on +ve z-axis?) - // For example, FHerwig only produces DIS events in the - // unconventional orientation and has to be corrected void Event::_geNormAlignment() { if (!_genEvent.valid_beam_particles()) return; typedef pair<HepMC::GenParticle*, HepMC::GenParticle*> GPPair; GPPair bps = _genEvent.beam_particles(); - //const PdgIdPair beamids = make_pdgid_pair(bps.first->pdg_id(), bps.second->pdg_id()); - //Log::getLog("Rivet.Event") << Log::TRACE << "Beam IDs: " << beamids << endl; - const HepMC::GenParticle* plusgp = 0; - bool rot = false; // Rotate e+- p and ppbar to put p along +z - /// @todo e+ e- convention? B-factories different from LEP? + /// @todo Is there an e+ e- convention for longitudinal boosting, e.g. at B-factories? Different from LEP? // if (compatible(beamids, make_pdgid_pair(ELECTRON, PROTON)) || // compatible(beamids, make_pdgid_pair(POSITRON, PROTON)) || // compatible(beamids, make_pdgid_pair(ANTIPROTON, PROTON)) ) { // Log::getLog("Rivet.Event") << Log::TRACE << "May need to rotate event..." << endl; + bool rot = false; + const HepMC::GenParticle* plusgp = 0; if (bps.first->pdg_id() != PROTON || bps.second->pdg_id() != PROTON) { if (bps.first->pdg_id() == PROTON) { plusgp = bps.first; @@ -52,58 +72,14 @@ // Do the rotation if (rot) { if (Log::getLog("Rivet.Event").isActive(Log::TRACE)) { - Log::getLog("Rivet.Event") << Log::TRACE << "Rotating event" << endl; - Log::getLog("Rivet.Event") << Log::TRACE << "Before rotation: " + Log::getLog("Rivet.Event") << Log::TRACE << "Rotating event\n" + << "Before rotation: " << bps.first->pdg_id() << "@pz=" << bps.first->momentum().pz()/GeV << ", " << bps.second->pdg_id() << "@pz=" << bps.second->momentum().pz()/GeV << endl; } - if (!_modGenEvent) _modGenEvent = new GenEvent(_genEvent); - _geRot180x(*_modGenEvent); + _geRot180x(_genEvent); } } - Event::Event(const GenEvent& ge) - : _genEvent(ge), _modGenEvent(NULL), _weight(1.0) - { - // Set the weight if there is one, otherwise default to 1.0 - if (!_genEvent.weights().empty()) { - _weight = ge.weights()[0]; - } - - // Use Rivet's preferred units if possible - #ifdef HEPMC_HAS_UNITS - if (_genEvent.momentum_unit() != HepMC::Units::GEV || - _genEvent.length_unit() != HepMC::Units::MM) { - if (!_modGenEvent) _modGenEvent = new GenEvent(ge); - _modGenEvent->use_units(HepMC::Units::GEV, HepMC::Units::MM); - } - #endif - - // Use the conventional alignment - _geNormAlignment(); - - // Debug printout to check that copying/mangling has worked - //_genEvent.print(); - } - - - Event::Event(const Event& e) - : _genEvent(e._genEvent), _modGenEvent(e._modGenEvent), - _weight(e._weight) - { - // - } - - - Event::~Event() { - if (_modGenEvent) delete _modGenEvent; - } - - const GenEvent& Event::genEvent() const { - if (_modGenEvent) return *_modGenEvent; - return _genEvent; - } - - } Modified: trunk/src/Core/Jet.cc ============================================================================== --- trunk/src/Core/Jet.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Core/Jet.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -56,9 +56,9 @@ bool Jet::containsParticle(const Particle& particle) const { - const int barcode = particle.genParticle().barcode(); + const int barcode = particle.genParticle()->barcode(); foreach (const Particle& p, particles()) { - if (p.genParticle().barcode() == barcode) return true; + if (p.genParticle()->barcode() == barcode) return true; } return false; } @@ -114,7 +114,7 @@ const PdgId pid = p.pdgId(); if (abs(pid) == CQUARK) return true; if (PID::isHadron(pid) && PID::hasCharm(pid)) return true; - HepMC::GenVertex* gv = p.genParticle().production_vertex(); + HepMC::GenVertex* gv = p.genParticle()->production_vertex(); if (gv) { foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { const PdgId pid2 = pi->pdg_id(); @@ -131,7 +131,7 @@ const PdgId pid = p.pdgId(); if (abs(pid) == BQUARK) return true; if (PID::isHadron(pid) && PID::hasBottom(pid)) return true; - HepMC::GenVertex* gv = p.genParticle().production_vertex(); + HepMC::GenVertex* gv = p.genParticle()->production_vertex(); if (gv) { foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { const PdgId pid2 = pi->pdg_id(); Modified: trunk/src/Core/Particle.cc ============================================================================== --- trunk/src/Core/Particle.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Core/Particle.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -5,7 +5,7 @@ bool Particle::hasAncestor(PdgId pdg_id) const { - GenVertex* prodVtx = genParticle().production_vertex(); + GenVertex* prodVtx = genParticle()->production_vertex(); if (prodVtx == 0) return false; foreach (const GenParticle* ancestor, particles(prodVtx, HepMC::ancestors)) { if (ancestor->pdg_id() == pdg_id) return true; Modified: trunk/src/Makefile.am ============================================================================== --- trunk/src/Makefile.am Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Makefile.am Thu Mar 7 17:20:02 2013 (r4211) @@ -1,4 +1,4 @@ -SUBDIRS = Core Projections Tools +SUBDIRS = Core Tools Projections if ENABLE_ANALYSES SUBDIRS += Analyses endif Modified: trunk/src/Projections/Beam.cc ============================================================================== --- trunk/src/Projections/Beam.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/Beam.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -48,15 +48,15 @@ void Beam::project(const Event& e) { - assert(e.genEvent().particles_size() >= 2); - if (e.genEvent().valid_beam_particles()) { - pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = e.genEvent().beam_particles(); + assert(e.genEvent()->particles_size() >= 2); + if (e.genEvent()->valid_beam_particles()) { + pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = e.genEvent()->beam_particles(); assert(beams.first && beams.second); _theBeams.first = *(beams.first); _theBeams.second = *(beams.second); - } else if(e.genEvent().barcode_to_particle(1) && e.genEvent().barcode_to_particle(2)) { - _theBeams.first = *(e.genEvent().barcode_to_particle(1)); - _theBeams.second = *(e.genEvent().barcode_to_particle(2)); + } else if(e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) { + _theBeams.first = *(e.genEvent()->barcode_to_particle(1)); + _theBeams.second = *(e.genEvent()->barcode_to_particle(2)); } else { _theBeams.first = Particle(ANY, FourMomentum()); Modified: trunk/src/Projections/DISFinalState.cc ============================================================================== --- trunk/src/Projections/DISFinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/DISFinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -15,9 +15,9 @@ // lepton, with momenta boosted into the appropriate frame. _theParticles.clear(); _theParticles.reserve(fs.particles().size()-1); - const GenParticle& dislepGP = dislep.out().genParticle(); + const GenParticle* dislepGP = dislep.out().genParticle(); foreach (const Particle& p, fs.particles()) { - if (&p.genParticle() != &dislepGP) { //< Ensure that we skip the DIS lepton + if (p.genParticle() != dislepGP) { //< Ensure that we skip the DIS lepton Particle temp(p); temp.setMomentum(hcmboost.transform(temp.momentum())); _theParticles.push_back(temp); Modified: trunk/src/Projections/DISLepton.cc ============================================================================== --- trunk/src/Projections/DISLepton.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/DISLepton.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -6,6 +6,7 @@ namespace Rivet { + int DISLepton::compare(const Projection& p) const { const DISLepton& other = pcast<DISLepton>(p); return @@ -16,12 +17,12 @@ void DISLepton::project(const Event& e) { const ParticlePair& inc = applyProjection<Beam>(e, "Beam").beams(); - + Particle inLepton; - + bool firstIsLepton = PID::isLepton(inc.first.pdgId()); bool secondIsLepton = PID::isLepton(inc.second.pdgId()); - + if(firstIsLepton && !secondIsLepton){ _incoming = inc.first; }else if(!firstIsLepton && secondIsLepton){ @@ -30,12 +31,12 @@ //eek! throw Error("DISLepton projector could not find the correct beam. "); } - + _sign = (_incoming.momentum().pz() > 0.0)? 1.0: -1.0; long id = _incoming.pdgId(); - + double pzMax = -1000000000.0; - + const FinalState& fs = applyProjection<FinalState>(e, "FS"); foreach (const Particle& p, fs.particles()) { double pz = _sign * p.momentum().pz(); @@ -44,9 +45,11 @@ pzMax = pz; } } - - if (!_outgoing.hasGenParticle()) { + + if (_outgoing.genParticle() == NULL) { throw Error("DISLepton projector could not find the scattered lepton."); } } + + } Modified: trunk/src/Projections/FinalState.cc ============================================================================== --- trunk/src/Projections/FinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/FinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -94,7 +94,7 @@ /// Decide if a particle is to be accepted or not. bool FinalState::accept(const Particle& p) const { // Not having s.c. == 1 should never happen! - assert(!p.hasGenParticle() || p.genParticle().status() == 1); + assert(p.genParticle() == NULL || p.genParticle()->status() == 1); // Check pT cut if (_ptmin > 0.0) { Modified: trunk/src/Projections/InvMassFinalState.cc ============================================================================== --- trunk/src/Projections/InvMassFinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/InvMassFinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -177,7 +177,7 @@ MSG_DEBUG("Selected " << _theParticles.size() << " particles " << "(" << _particlePairs.size() << " pairs)"); if (getLog().isActive(Log::TRACE)) { foreach (const Particle& p, _theParticles) { - MSG_TRACE("ID: " << p.pdgId() << ", barcode: " << p.genParticle().barcode()); + MSG_TRACE("ID: " << p.pdgId() << ", barcode: " << p.genParticle()->barcode()); } } } Modified: trunk/src/Projections/LeadingParticlesFinalState.cc ============================================================================== --- trunk/src/Projections/LeadingParticlesFinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/LeadingParticlesFinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -35,7 +35,7 @@ MSG_DEBUG("Original final state particles size " << particles.size()); ParticleVector::const_iterator ifs; for (ifs = particles.begin(); ifs != particles.end(); ++ifs) { - if (inList(*ifs) && FinalState::accept(ifs->genParticle())) { + if (inList(*ifs) && FinalState::accept(*ifs->genParticle())) { // Look for an existing particle in tmp container map < long, ParticleVector::const_iterator >::const_iterator itmp = tmp.find(ifs->pdgId()); if (itmp != tmp.end()) { // if a particle with this type has been already selected @@ -49,7 +49,7 @@ } } } - + // Loop on the tmp container and fill _theParticles map<long, ParticleVector::const_iterator>::const_iterator i; for (i = tmp.begin(); i != tmp.end(); ++i) { Modified: trunk/src/Projections/MergedFinalState.cc ============================================================================== --- trunk/src/Projections/MergedFinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/MergedFinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -22,10 +22,10 @@ _theParticles.push_back(pa); } foreach (const Particle& pb, fsb.particles()){ - const GenParticle* originalb = &(pb.genParticle()); + const GenParticle* originalb = pb.genParticle(); bool notfound = true; foreach (const Particle& pa, fsa.particles()){ - const GenParticle* originala = &(pa.genParticle()); + const GenParticle* originala = pa.genParticle(); if (originala == originalb) { notfound = false; break; Modified: trunk/src/Projections/VetoedFinalState.cc ============================================================================== --- trunk/src/Projections/VetoedFinalState.cc Thu Mar 7 12:59:39 2013 (r4210) +++ trunk/src/Projections/VetoedFinalState.cc Thu Mar 7 17:20:02 2013 (r4211) @@ -110,9 +110,10 @@ /// @todo Improve! for (ParentVetos::const_iterator vIt = _parentVetoes.begin(); vIt != _parentVetoes.end(); ++vIt) { for (ParticleVector::iterator p = _theParticles.begin(); p != _theParticles.end(); ++p) { - GenVertex* startVtx = ((*p).genParticle()).production_vertex(); + GenVertex* startVtx = p->genParticle()->production_vertex(); bool veto = false; if (startVtx!=0) { + /// @todo Use better HepMC iteration for (GenVertex::particle_iterator pIt = startVtx->particles_begin(HepMC::ancestors); pIt != startVtx->particles_end(HepMC::ancestors) && !veto; ++pIt) { if (*vIt == (*pIt)->pdg_id()) { @@ -130,13 +131,13 @@ const FinalState& vfs = applyProjection<FinalState>(e, ifs); const ParticleVector& vfsp = vfs.particles(); for (ParticleVector::iterator icheck = _theParticles.begin(); icheck != _theParticles.end(); ++icheck) { - if (!icheck->hasGenParticle()) continue; + if (icheck->genParticle() == NULL) continue; bool found = false; for (ParticleVector::const_iterator ipart = vfsp.begin(); ipart != vfsp.end(); ++ipart){ - if (!ipart->hasGenParticle()) continue; - MSG_TRACE("Comparing barcode " << icheck->genParticle().barcode() - << " with veto particle " << ipart->genParticle().barcode()); - if (ipart->genParticle().barcode() == icheck->genParticle().barcode()){ + if (ipart->genParticle() == NULL) continue; + MSG_TRACE("Comparing barcode " << icheck->genParticle()->barcode() + << " with veto particle " << ipart->genParticle()->barcode()); + if (ipart->genParticle()->barcode() == icheck->genParticle()->barcode()){ found = true; break; }
More information about the Rivet-svn mailing list |