|
[Rivet-svn] r3268 - in trunk: include/Rivet/Projections src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgMon Aug 1 18:28:50 BST 2011
Author: fsiegert Date: Mon Aug 1 18:28:49 2011 New Revision: 3268 Log: Another overhaul in the never-ending [WZ]Finder saga. Now I am not aware of any inconsistencies anymore and the chaining of projections is much more readable now IMO. This does not seem to have a lot of influence on observables (given the recent fixes on trunk are used as well), a comparison to Rivet 1.5.0, 1.5.1 and latest trunk is at http://fsiegert.web.cern.ch/fsiegert/tmp/wzfinder This changes the interface significantly, so we will definitely need a 1.6.0 instead of 1.5.2 (-> rewrite the ChangeLog). Furthermore I'm not sure whether we should rename these projections, since their _theParticles vector is not the W or Z anymore, but the constituent original particles. This was the only way I saw to recover the remaining final state without the bugs. Of course one can still easily access the W/Z stored as _bosons. At least this should be mentioned prominently in any updates notes. Modified: trunk/include/Rivet/Projections/InvMassFinalState.hh trunk/include/Rivet/Projections/LeptonClusters.hh trunk/include/Rivet/Projections/WFinder.hh trunk/include/Rivet/Projections/ZFinder.hh trunk/src/Analyses/CDF_2000_S4155203.cc trunk/src/Analyses/CDF_2009_S8383952.cc trunk/src/Analyses/D0_2000_S4480767.cc trunk/src/Analyses/D0_2007_S7075677.cc trunk/src/Analyses/D0_2008_S6879055.cc trunk/src/Analyses/D0_2008_S7554427.cc trunk/src/Analyses/D0_2008_S7837160.cc trunk/src/Analyses/D0_2008_S7863608.cc trunk/src/Analyses/D0_2009_S8202443.cc trunk/src/Analyses/D0_2009_S8349509.cc trunk/src/Analyses/D0_2010_S8671338.cc trunk/src/Analyses/D0_2010_S8821313.cc trunk/src/Analyses/MC_HJETS.cc trunk/src/Analyses/MC_WJETS.cc trunk/src/Analyses/MC_WPOL.cc trunk/src/Analyses/MC_WWJETS.cc trunk/src/Analyses/MC_ZJETS.cc trunk/src/Analyses/MC_ZZJETS.cc trunk/src/Projections/InvMassFinalState.cc trunk/src/Projections/LeptonClusters.cc trunk/src/Projections/WFinder.cc trunk/src/Projections/ZFinder.cc Modified: trunk/include/Rivet/Projections/InvMassFinalState.hh ============================================================================== --- trunk/include/Rivet/Projections/InvMassFinalState.hh Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/include/Rivet/Projections/InvMassFinalState.hh Mon Aug 1 18:28:49 2011 (r3268) @@ -44,6 +44,9 @@ _useTransverseMass = usetrans; } + /// Operate on a given particle vector directly instead of through project (no caching) + void calc(const ParticleVector& inparticles); + protected: Modified: trunk/include/Rivet/Projections/LeptonClusters.hh ============================================================================== --- trunk/include/Rivet/Projections/LeptonClusters.hh Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/include/Rivet/Projections/LeptonClusters.hh Mon Aug 1 18:28:49 2011 (r3268) @@ -12,61 +12,31 @@ namespace Rivet { - /// @brief Cluster photons from fs to all charged particles (typically - /// leptons) from signal and store the original charged particles and photons - /// as particles() while the newly created clustered lepton objects are - /// accessible as clusteredLeptons() - class LeptonClustersConstituents : public FinalState { - public: - - /// @name Constructors - //@{ - LeptonClustersConstituents(const FinalState& fs, const FinalState& signal, - double dRmax, bool cluster, bool track) - : _dRmax(dRmax), _cluster(cluster), _track(track) - { - setName("LeptonClustersConstituents"); - IdentifiedFinalState photonfs(fs); - photonfs.acceptId(PHOTON); - addProjection(photonfs, "Photons"); - addProjection(signal, "Signal"); - } - - - /// Clone on the heap. - virtual const Projection* clone() const { - return new LeptonClustersConstituents(*this); - } - //@} - + class ClusteredLepton : public Particle { public: + ClusteredLepton(Particle lepton) : + Particle(lepton.pdgId(), lepton.momentum()), + _constituentLepton(lepton) {} + + void addPhoton(const Particle& p, bool cluster) { + _constituentPhotons.push_back(p); + if (cluster) setMomentum(momentum()+p.momentum()); + } - const ParticleVector& clusteredLeptons() const { return _clusteredLeptons; } - - protected: - - /// Apply the projection on the supplied event. - void project(const Event& e); - - /// Compare projections. - int compare(const Projection& p) const; - + const Particle& constituentLepton() const { return _constituentLepton; } + const ParticleVector& constituentPhotons() const { return _constituentPhotons; } private: - - /// maximum cone radius to find photons in - double _dRmax; - /// whether to actually add the photon momenta to clusteredLeptons - bool _cluster; - /// whether to add the photons to the particles() vector - bool _track; - ParticleVector _clusteredLeptons; - + ParticleVector _constituentPhotons; + Particle _constituentLepton; }; - + /// @brief Cluster photons from a given FS to all charged particles (typically + /// leptons) from signal and store the original charged particles and photons + /// as particles() while the newly created clustered lepton objects are + /// accessible as clusteredLeptons() /// @brief Use LeptonClustersConstituents projection to cluster all photons to /// leptons. This projection here makes the clustered objects available /// in a FinalState manner, i.e. as particles(). The given pT and eta cuts are @@ -76,22 +46,16 @@ public: - LeptonClusters(const FinalState& fs, const FinalState& signal, double dRmax, - bool cluster, bool track, + LeptonClusters(const FinalState& photons, const FinalState& signal, + double dRmax, bool cluster, const std::vector<std::pair<double, double> >& etaRanges, - double pTmin) : - FinalState(etaRanges, pTmin) - { - setName("LeptonClusters"); - LeptonClustersConstituents constituents(fs, signal, dRmax, cluster, track); - addProjection(constituents, "Constituents"); - } + double pTmin); virtual const Projection* clone() const { return new LeptonClusters(*this); } - const FinalState& constituentsFinalState() const; + const vector<ClusteredLepton>& clusteredLeptons() const { return _clusteredLeptons; } protected: @@ -100,9 +64,21 @@ /// Compare projections. int compare(const Projection& p) const; + + private: + + /// maximum cone radius to find photons in + double _dRmax; + /// whether to actually add the photon momenta to clusteredLeptons + bool _cluster; + + /// container which stores the clustered lepton objects + vector<ClusteredLepton> _clusteredLeptons; }; + + } Modified: trunk/include/Rivet/Projections/WFinder.hh ============================================================================== --- trunk/include/Rivet/Projections/WFinder.hh Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/include/Rivet/Projections/WFinder.hh Mon Aug 1 18:28:49 2011 (r3268) @@ -8,6 +8,7 @@ #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/LeptonClusters.hh" namespace Rivet { @@ -27,7 +28,7 @@ // /// reconstruction. // WFinder(const ChargedFinalState& fs_l, // PdgId pid, - // double m2_min, double m2_max, + // double minmass, double maxmass, // double missingET, // double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); @@ -38,9 +39,11 @@ WFinder(double etaMin, double etaMax, double pTmin, PdgId pid, - double m2_min, double m2_max, + double minmass, double maxmass, double missingET, - double dRmax, bool clusterPhotons=true, bool excludePhotonsFromRFS=false); + double dRmax, bool clusterPhotons=true, bool trackPhotons=false, + double masstarget=80.4, + bool useTransverseMass=false); /// Constructor taking multiple eta/pT bounds and type of the charged lepton, transverse mass @@ -49,9 +52,11 @@ WFinder(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, const double m2_max, + double minmass, const double maxmass, double missingET, - double dRmax, bool clusterPhotons=true, bool excludePhotonsFromRFS=false); + double dRmax, bool clusterPhotons=true, bool trackPhotons=false, + double masstarget=80.4, + bool useTransverseMass=false); /// Clone on the heap. @@ -61,18 +66,22 @@ //@} + /// Access to the found bosons (currently either 0 or 1) + const ParticleVector& bosons() const { return _bosons; } + + /// Access to the W constituent clustered leptons (currently either of + /// size 0 if no boson was found or 1 if one boson was found) + const vector<Particle>& constituentLeptons() const { return _constituentLeptons; } + + /// Access to the W constituent neutrinos (currently either of size 0 if no + /// boson was found or 1 if one boson was found) + const vector<Particle>& constituentNeutrinos() const { return _constituentNeutrinos; } + /// Access to the remaining particles, after the W and clustered photons /// have been removed from the full final state /// (e.g. for running a jet finder on it) const FinalState& remainingFinalState() const; - /// Access to the W constituent clustered lepton - /// (e.g. for more fine-grained cuts on the clustered lepton) - Particle constituentLepton() const; - Particle constituentNeutrino() const; - - const FinalState& originalLeptonFinalState() const; - protected: /// Apply the projection on the supplied event. @@ -87,6 +96,9 @@ /// Clear the projection void clear() { _theParticles.clear(); + _bosons.clear(); + _constituentLeptons.clear(); + _constituentNeutrinos.clear(); } @@ -95,19 +107,45 @@ /// Common implementation of constructor operation, taking FS params. void _init(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, double m2_max, + double minmass, double maxmass, double missingET, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget, + bool useTransverseMass); private: - // Mass range - double _m2_min, _m2_max; + /// Transverse mass cuts + double _minmass, _maxmass, _masstarget; + bool _useTransverseMass; - // Missing ET cut + /// Missing ET cut double _etMiss; + /// Switch for tracking of photons (whether to add them to _theParticles) + /// This is relevant when the ZFinder::_theParticles are to be excluded + /// from e.g. the input to a jet finder, to specify whether the clustered + /// photons are to be excluded as well. + /// (Yes, some experiments make a difference between clusterPhotons and + /// trackPhotons!) + bool _trackPhotons; + + /// Lepton flavour + PdgId _pid; + + /// Neutrino flavour + PdgId _nu_pid; + + /// list of found bosons (currently either 0 or 1) + ParticleVector _bosons; + + /// Constituent leptons (currently either 0 or 1) + ParticleVector _constituentLeptons; + + /// Constituent neutrinos (currently either 0 or 1) + ParticleVector _constituentNeutrinos; + }; Modified: trunk/include/Rivet/Projections/ZFinder.hh ============================================================================== --- trunk/include/Rivet/Projections/ZFinder.hh Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/include/Rivet/Projections/ZFinder.hh Mon Aug 1 18:28:49 2011 (r3268) @@ -8,6 +8,7 @@ #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/LeptonClusters.hh" namespace Rivet { @@ -32,8 +33,9 @@ ZFinder(double etaMin, double etaMax, double pTmin, PdgId pid, - double m2_min, double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); + double minmass, double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget=91.2*GeV); /// Constructor taking multiple eta/pT bounds and type of the leptons, mass @@ -45,8 +47,9 @@ ZFinder(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, const double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); + double minmass, const double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget=91.2*GeV); /// Clone on the heap. @@ -56,16 +59,18 @@ //@} + /// Access to the found bosons (currently either 0 or 1) + const ParticleVector& bosons() const { return _bosons; } + + /// Access to the Z constituent clustered leptons + /// (e.g. for more fine-grained cuts on the clustered leptons) + const vector<Particle>& constituents() const { return _constituents; } + /// Access to the remaining particles, after the Z and clustered photons /// have been removed from the full final state /// (e.g. for running a jet finder on it) const FinalState& remainingFinalState() const; - /// Access to the Z constituent clustered leptons final state - /// (e.g. for more fine-grained cuts on the clustered leptons) - const FinalState& constituentsFinalState() const; - - const FinalState& originalConstituentsFinalState() const; protected: @@ -76,12 +81,43 @@ int compare(const Projection& p) const; + public: + + /// Clear the projection + void clear() { + _theParticles.clear(); + _bosons.clear(); + _constituents.clear(); + } + + private: /// Common implementation of constructor operation, taking FS params. void _init(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromJets); + double minmass, double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget); + + /// Mass cuts to apply to clustered leptons (cf. InvMassFinalState) + double _minmass, _maxmass, _masstarget; + + /// Switch for tracking of photons (whether to add them to _theParticles) + /// This is relevant when the ZFinder::_theParticles are to be excluded + /// from e.g. the input to a jet finder, to specify whether the clustered + /// photons are to be excluded as well. + /// (Yes, some experiments make a difference between clusterPhotons and + /// trackPhotons!) + bool _trackPhotons; + + /// Lepton flavour + PdgId _pid; + + /// list of found bosons (currently either 0 or 1) + ParticleVector _bosons; + + /// Clustered leptons + vector<Particle> _constituents; }; @@ -89,4 +125,5 @@ } + #endif Modified: trunk/src/Analyses/CDF_2000_S4155203.cc ============================================================================== --- trunk/src/Analyses/CDF_2000_S4155203.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/CDF_2000_S4155203.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -40,12 +40,12 @@ /// Do the analysis void analyze(const Event& e) { const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size() != 1) { - MSG_DEBUG("Num e+ e- pairs found = " << zfinder.particles().size()); + if (zfinder.bosons().size() != 1) { + MSG_DEBUG("Num e+ e- pairs found = " << zfinder.bosons().size()); vetoEvent; } - FourMomentum pZ = zfinder.particles()[0].momentum(); + FourMomentum pZ = zfinder.bosons()[0].momentum(); if (pZ.mass2() < 0) { MSG_DEBUG("Negative Z mass**2 = " << pZ.mass2()/GeV2 << "!"); vetoEvent; Modified: trunk/src/Analyses/CDF_2009_S8383952.cc ============================================================================== --- trunk/src/Analyses/CDF_2009_S8383952.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/CDF_2009_S8383952.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -54,8 +54,8 @@ const double weight = event.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); - if (zfinder.particles().size() == 1) { - double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); + if (zfinder.bosons().size() == 1) { + double yZ = fabs(zfinder.bosons()[0].momentum().rapidity()); _h_yZ->fill(yZ, weight); _h_xs->fill(1960.0, weight); } Modified: trunk/src/Analyses/D0_2000_S4480767.cc ============================================================================== --- trunk/src/Analyses/D0_2000_S4480767.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2000_S4480767.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -44,9 +44,9 @@ const double weight = event.weight(); const WFinder& wf = applyProjection<WFinder>(event, "WFinder"); - if (wf.size() == 0) vetoEvent; + if (wf.bosons().size() == 0) vetoEvent; - _h_W_pT->fill(wf.particles()[0].momentum().pT()/GeV, weight); + _h_W_pT->fill(wf.bosons()[0].momentum().pT()/GeV, weight); } Modified: trunk/src/Analyses/D0_2007_S7075677.cc ============================================================================== --- trunk/src/Analyses/D0_2007_S7075677.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2007_S7075677.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -41,10 +41,10 @@ const double weight = e.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size() == 1) { - const ParticleVector& el(zfinder.constituentsFinalState().particles()); + if (zfinder.bosons().size() == 1) { + const ParticleVector& el(zfinder.constituents()); if (el[0].momentum().pT() > 25.0*GeV || el[1].momentum().pT() > 25.0*GeV) { - double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); + double yZ = fabs(zfinder.bosons()[0].momentum().rapidity()); _h_yZ->fill(yZ, weight); } } Modified: trunk/src/Analyses/D0_2008_S6879055.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S6879055.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2008_S6879055.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -46,12 +46,12 @@ const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); - if (zfinder.particles().size()!=1) { + if (zfinder.bosons().size()!=1) { vetoEvent; } - FourMomentum e0 = zfinder.constituentsFinalState().particles()[0].momentum(); - FourMomentum e1 = zfinder.constituentsFinalState().particles()[1].momentum(); + FourMomentum e0 = zfinder.constituents()[0].momentum(); + FourMomentum e1 = zfinder.constituents()[1].momentum(); const double e0eta = e0.eta(); const double e0phi = e0.phi(); const double e1eta = e1.eta(); Modified: trunk/src/Analyses/D0_2008_S7554427.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7554427.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2008_S7554427.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -44,9 +44,9 @@ const double weight = e.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size() == 1) { - double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); - double pTZ = zfinder.particles()[0].momentum().pT(); + if (zfinder.bosons().size() == 1) { + double yZ = fabs(zfinder.bosons()[0].momentum().rapidity()); + double pTZ = zfinder.bosons()[0].momentum().pT(); _h_ZpT->fill(pTZ, weight); if (yZ > 2.0) { _h_forward_ZpT->fill(pTZ, weight); Modified: trunk/src/Analyses/D0_2008_S7837160.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7837160.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2008_S7837160.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -54,7 +54,7 @@ /// Do the analysis void analyze(const Event & event) { const WFinder& wf = applyProjection<WFinder>(event, "WFe"); - if (wf.size() == 0) { + if (wf.bosons().size() == 0) { getLog() << Log::DEBUG << "No W candidates found: vetoing" << endl; vetoEvent; } @@ -62,8 +62,8 @@ // Require that leptons have Et >= 25 GeV /// @todo Use pT cut in WFinder /// @todo Any ETmiss cut? - FourMomentum p_e=wf.constituentLepton().momentum(); - int chg_e = PID::threeCharge(wf.constituentLepton().pdgId()); + FourMomentum p_e=wf.constituentLeptons()[0].momentum(); + int chg_e = PID::threeCharge(wf.constituentLeptons()[0].pdgId()); if (p_e.eta() < 0) chg_e *= -1; assert(chg_e != 0); Modified: trunk/src/Analyses/D0_2008_S7863608.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7863608.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2008_S7863608.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -56,7 +56,7 @@ const double weight = e.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()==1) { + if (zfinder.bosons().size()==1) { _sum_of_weights_inclusive += weight; const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); const Jets& jets = jetpro.jetsByPt(20.0*GeV); @@ -74,7 +74,7 @@ vetoEvent; } - const FourMomentum Zmom = zfinder.particles()[0].momentum(); + const FourMomentum Zmom = zfinder.bosons()[0].momentum(); // In jet pT _h_jet_pT_cross_section->fill( jets_cut[0].momentum().pT(), weight); Modified: trunk/src/Analyses/D0_2009_S8202443.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8202443.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2009_S8202443.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -65,7 +65,7 @@ // unconstrained electrons first const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()==1) { + if (zfinder.bosons().size()==1) { _sum_of_weights += weight; const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); const Jets& jets = jetpro.jetsByPt(20.0*GeV); @@ -93,7 +93,7 @@ // constrained electrons const ZFinder& zfinder_constrained = applyProjection<ZFinder>(e, "ZFinderConstrained"); - if (zfinder_constrained.particles().size()==1) { + if (zfinder_constrained.bosons().size()==1) { _sum_of_weights_constrained += weight; const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinderConstrained"); const Jets& jets = jetpro.jetsByPt(20.0*GeV); Modified: trunk/src/Analyses/D0_2009_S8349509.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8349509.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2009_S8349509.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -64,11 +64,11 @@ const double weight = event.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); - if (zfinder.particles().size()==1) { + if (zfinder.bosons().size()==1) { // count inclusive sum of weights for histogram normalisation _inclusive_Z_sumofweights += weight; - const FourMomentum Zmom = zfinder.particles()[0].momentum(); + const FourMomentum Zmom = zfinder.bosons()[0].momentum(); if (Zmom.pT()<25.0*GeV) { vetoEvent; } Modified: trunk/src/Analyses/D0_2010_S8671338.cc ============================================================================== --- trunk/src/Analyses/D0_2010_S8671338.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2010_S8671338.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -44,8 +44,8 @@ void analyze(const Event& e) { const double weight = e.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()==1) { - double ZpT = zfinder.particles()[0].momentum().pT()/GeV; + if (zfinder.bosons().size()==1) { + double ZpT = zfinder.bosons()[0].momentum().pT()/GeV; _h_Z_pT_normalised->fill(ZpT, weight); _h_Z_pT_xs->fill(ZpT, weight); } Modified: trunk/src/Analyses/D0_2010_S8821313.cc ============================================================================== --- trunk/src/Analyses/D0_2010_S8821313.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/D0_2010_S8821313.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -61,8 +61,9 @@ const double weight = event.weight(); const ZFinder& zfinder_ee = applyProjection<ZFinder>(event, "zfinder_ee"); - if (zfinder_ee.particles().size()==1) { - ParticleVector ee=zfinder_ee.constituentsFinalState().particlesByPt(); + if (zfinder_ee.bosons().size()==1) { + ParticleVector ee=zfinder_ee.constituents(); + std::sort(ee.begin(), ee.end(), cmpParticleByPt); FourMomentum eminus=PID::threeCharge(ee[0].pdgId())<0.0?ee[0].momentum():ee[1].momentum(); FourMomentum eplus=PID::threeCharge(ee[0].pdgId())<0.0?ee[1].momentum():ee[0].momentum(); double phi_acop=M_PI-mapAngle0ToPi(eminus.phi()-eplus.phi()); @@ -71,13 +72,14 @@ if (sin2thetastar<0.0) sin2thetastar=0.0; double phistar=tan(phi_acop/2.0)*sqrt(sin2thetastar); - FourMomentum Zmom=zfinder_ee.particles()[0].momentum(); + FourMomentum Zmom=zfinder_ee.bosons()[0].momentum(); _h_phistar_ee.fill(Zmom.rapidity(), phistar, weight); } const ZFinder& zfinder_mm = applyProjection<ZFinder>(event, "zfinder_mm"); - if (zfinder_mm.particles().size()==1) { - ParticleVector mm=zfinder_mm.constituentsFinalState().particlesByPt(); + if (zfinder_mm.bosons().size()==1) { + ParticleVector mm=zfinder_mm.constituents(); + std::sort(mm.begin(), mm.end(), cmpParticleByPt); FourMomentum mminus=PID::threeCharge(mm[0].pdgId())<0.0?mm[0].momentum():mm[1].momentum(); FourMomentum mplus=PID::threeCharge(mm[0].pdgId())<0.0?mm[1].momentum():mm[0].momentum(); double phi_acop=M_PI-mapAngle0ToPi(mminus.phi()-mplus.phi()); @@ -86,7 +88,7 @@ if (sin2thetastar<0.0) sin2thetastar=0.0; double phistar=tan(phi_acop/2.0)*sqrt(sin2thetastar); - FourMomentum Zmom=zfinder_mm.particles()[0].momentum(); + FourMomentum Zmom=zfinder_mm.bosons()[0].momentum(); _h_phistar_mm.fill(Zmom.rapidity(), phistar, weight); } } Modified: trunk/src/Analyses/MC_HJETS.cc ============================================================================== --- trunk/src/Analyses/MC_HJETS.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_HJETS.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -48,18 +48,18 @@ /// Do the analysis void analyze(const Event & e) { const ZFinder& hfinder = applyProjection<ZFinder>(e, "Hfinder"); - if (hfinder.particles().size()!=1) { + if (hfinder.bosons().size()!=1) { vetoEvent; } const double weight = e.weight(); - FourMomentum hmom(hfinder.particles()[0].momentum()); + FourMomentum hmom(hfinder.bosons()[0].momentum()); _h_H_mass->fill(hmom.mass(),weight); _h_H_pT->fill(hmom.pT(),weight); _h_H_pT_peak->fill(hmom.pT(),weight); _h_H_y->fill(hmom.rapidity(),weight); _h_H_phi->fill(hmom.azimuthalAngle(),weight); - foreach (const Particle& l, hfinder.constituentsFinalState().particles()) { + foreach (const Particle& l, hfinder.constituents()) { _h_lepton_pT->fill(l.momentum().pT(), weight); _h_lepton_eta->fill(l.momentum().eta(), weight); } Modified: trunk/src/Analyses/MC_WJETS.cc ============================================================================== --- trunk/src/Analyses/MC_WJETS.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_WJETS.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -53,7 +53,7 @@ /// Do the analysis void analyze(const Event & e) { const WFinder& wfinder = applyProjection<WFinder>(e, "WFinder"); - if (wfinder.particles().size() != 1) { + if (wfinder.bosons().size() != 1) { vetoEvent; } const double weight = e.weight(); @@ -61,13 +61,13 @@ double charge3_x_eta = 0; int charge3 = 0; FourMomentum emom; - FourMomentum wmom(wfinder.particles().front().momentum()); + FourMomentum wmom(wfinder.bosons().front().momentum()); _h_W_mass->fill(wmom.mass(), weight); _h_W_pT->fill(wmom.pT(), weight); _h_W_pT_peak->fill(wmom.pT(), weight); _h_W_y->fill(wmom.rapidity(), weight); _h_W_phi->fill(wmom.azimuthalAngle(), weight); - Particle l=wfinder.constituentLepton(); + Particle l=wfinder.constituentLeptons()[0]; _h_lepton_pT->fill(l.momentum().pT(), weight); _h_lepton_eta->fill(l.momentum().eta(), weight); if (PID::threeCharge(l.pdgId()) != 0) { Modified: trunk/src/Analyses/MC_WPOL.cc ============================================================================== --- trunk/src/Analyses/MC_WPOL.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_WPOL.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -72,16 +72,16 @@ const double weight = event.weight(); const WFinder& wfinder = applyProjection<WFinder>(event, "WFinder"); - if (wfinder.particles().size() != 1) { + if (wfinder.bosons().size() != 1) { vetoEvent; } const ParticlePair& beams = applyProjection<Beam>(event, "Beams").beams(); FourMomentum pb1(beams.second.momentum()), pb2(beams.first.momentum()); - Particle lepton=wfinder.constituentLepton(); + Particle lepton=wfinder.constituentLeptons()[0]; FourMomentum pl(lepton.momentum()); size_t idx = (PID::threeCharge(lepton.pdgId())>0 ? 0 : 1); - FourMomentum plnu(wfinder.particles()[0].momentum()); + FourMomentum plnu(wfinder.bosons()[0].momentum()); LorentzTransform cms(-plnu.boostVector()); Matrix3 zrot(plnu.vector3(), Vector3(0.0, 0.0, 1.0)); @@ -98,7 +98,7 @@ Matrix3 xrot(Vector3(xref.x(), xref.y(), 0.0), Vector3(1.0, 0.0, 0.0)); pl3=xrot*pl3; - double ptw(wfinder.particles()[0].momentum().pT()/GeV); + double ptw(wfinder.bosons()[0].momentum().pT()/GeV); double thetas(pl3.theta()), phis(pl3.phi()); double costhetas(cos(thetas)), sinthetas(sin(thetas)); double cosphis(cos(phis)), sinphis(sin(phis)); Modified: trunk/src/Analyses/MC_WWJETS.cc ============================================================================== --- trunk/src/Analyses/MC_WWJETS.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_WWJETS.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -32,8 +32,8 @@ addProjection(wmnufinder, "WmnuFinder"); VetoedFinalState jetinput; jetinput - .addVetoOnThisFinalState(wenufinder.originalLeptonFinalState()) - .addVetoOnThisFinalState(wmnufinder.originalLeptonFinalState()); + .addVetoOnThisFinalState(wenufinder) + .addVetoOnThisFinalState(wmnufinder); FastJets jetpro(jetinput, FastJets::KT, 0.7); addProjection(jetpro, "Jets"); @@ -88,23 +88,23 @@ const double weight = e.weight(); const WFinder& wenufinder = applyProjection<WFinder>(e, "WenuFinder"); - if (wenufinder.particles().size()!=1) { + if (wenufinder.bosons().size()!=1) { vetoEvent; } const WFinder& wmnufinder = applyProjection<WFinder>(e, "WmnuFinder"); - if (wmnufinder.particles().size()!=1) { + if (wmnufinder.bosons().size()!=1) { vetoEvent; } - FourMomentum wenu(wenufinder.particles()[0].momentum()); - FourMomentum wmnu(wmnufinder.particles()[0].momentum()); + FourMomentum wenu(wenufinder.bosons()[0].momentum()); + FourMomentum wmnu(wmnufinder.bosons()[0].momentum()); FourMomentum ww(wenu+wmnu); // find leptons - FourMomentum ep=wenufinder.constituentLepton().momentum(); - FourMomentum enu=wenufinder.constituentNeutrino().momentum(); - FourMomentum mm=wmnufinder.constituentLepton().momentum(); - FourMomentum mnu=wmnufinder.constituentNeutrino().momentum(); + FourMomentum ep=wenufinder.constituentLeptons()[0].momentum(); + FourMomentum enu=wenufinder.constituentNeutrinos()[0].momentum(); + FourMomentum mm=wmnufinder.constituentLeptons()[0].momentum(); + FourMomentum mnu=wmnufinder.constituentNeutrinos()[0].momentum(); _h_WW_pT->fill(ww.pT(),weight); _h_WW_pT_peak->fill(ww.pT(),weight); Modified: trunk/src/Analyses/MC_ZJETS.cc ============================================================================== --- trunk/src/Analyses/MC_ZJETS.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_ZJETS.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -48,18 +48,18 @@ /// Do the analysis void analyze(const Event & e) { const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()!=1) { + if (zfinder.bosons().size()!=1) { vetoEvent; } const double weight = e.weight(); - FourMomentum zmom(zfinder.particles()[0].momentum()); + FourMomentum zmom(zfinder.bosons()[0].momentum()); _h_Z_mass->fill(zmom.mass(),weight); _h_Z_pT->fill(zmom.pT(),weight); _h_Z_pT_peak->fill(zmom.pT(),weight); _h_Z_y->fill(zmom.rapidity(),weight); _h_Z_phi->fill(zmom.azimuthalAngle(),weight); - foreach (const Particle& l, zfinder.constituentsFinalState().particles()) { + foreach (const Particle& l, zfinder.constituents()) { _h_lepton_pT->fill(l.momentum().pT(), weight); _h_lepton_eta->fill(l.momentum().eta(), weight); } Modified: trunk/src/Analyses/MC_ZZJETS.cc ============================================================================== --- trunk/src/Analyses/MC_ZZJETS.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Analyses/MC_ZZJETS.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -32,8 +32,8 @@ addProjection(zmmfinder, "ZmmFinder"); VetoedFinalState jetinput; jetinput - .addVetoOnThisFinalState(zeefinder.originalConstituentsFinalState()) - .addVetoOnThisFinalState(zmmfinder.originalConstituentsFinalState()); + .addVetoOnThisFinalState(zeefinder) + .addVetoOnThisFinalState(zmmfinder); FastJets jetpro(jetinput, FastJets::KT, 0.7); addProjection(jetpro, "Jets"); @@ -86,36 +86,36 @@ const double weight = e.weight(); const ZFinder& zeefinder = applyProjection<ZFinder>(e, "ZeeFinder"); - if (zeefinder.particles().size()!=1) { + if (zeefinder.bosons().size()!=1) { vetoEvent; } const ZFinder& zmmfinder = applyProjection<ZFinder>(e, "ZmmFinder"); - if (zmmfinder.particles().size()!=1) { + if (zmmfinder.bosons().size()!=1) { vetoEvent; } - FourMomentum zee(zeefinder.particles()[0].momentum()); - FourMomentum zmm(zmmfinder.particles()[0].momentum()); + FourMomentum zee(zeefinder.bosons()[0].momentum()); + FourMomentum zmm(zmmfinder.bosons()[0].momentum()); FourMomentum zz(zee+zmm); // find leptons FourMomentum ep(0.0,0.0,0.0,0.0), em(0.0,0.0,0.0,0.0); - if (PID::threeCharge(zeefinder.constituentsFinalState().particles()[0])>0.0) { - ep=zeefinder.constituentsFinalState().particles()[0].momentum(); - em=zeefinder.constituentsFinalState().particles()[1].momentum(); + if (PID::threeCharge(zeefinder.constituents()[0])>0.0) { + ep=zeefinder.constituents()[0].momentum(); + em=zeefinder.constituents()[1].momentum(); } else { - ep=zeefinder.constituentsFinalState().particles()[1].momentum(); - em=zeefinder.constituentsFinalState().particles()[0].momentum(); + ep=zeefinder.constituents()[1].momentum(); + em=zeefinder.constituents()[0].momentum(); } FourMomentum mp(0.0,0.0,0.0,0.0), mm(0.0,0.0,0.0,0.0); - if (PID::threeCharge(zmmfinder.constituentsFinalState().particles()[0])>0.0) { - mp=zmmfinder.constituentsFinalState().particles()[0].momentum(); - mm=zmmfinder.constituentsFinalState().particles()[1].momentum(); + if (PID::threeCharge(zmmfinder.constituents()[0])>0.0) { + mp=zmmfinder.constituents()[0].momentum(); + mm=zmmfinder.constituents()[1].momentum(); } else { - mp=zmmfinder.constituentsFinalState().particles()[1].momentum(); - mm=zmmfinder.constituentsFinalState().particles()[0].momentum(); + mp=zmmfinder.constituents()[1].momentum(); + mm=zmmfinder.constituents()[0].momentum(); } _h_ZZ_pT->fill(zz.pT(),weight); Modified: trunk/src/Projections/InvMassFinalState.cc ============================================================================== --- trunk/src/Projections/InvMassFinalState.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Projections/InvMassFinalState.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -60,15 +60,19 @@ void InvMassFinalState::project(const Event& e) { + const FinalState& fs = applyProjection<FinalState>(e, "FS"); + calc(fs.particles()); + } + + void InvMassFinalState::calc(const ParticleVector& inparticles) { _theParticles.clear(); _particlePairs.clear(); - const FinalState& fs = applyProjection<FinalState>(e, "FS"); // Containers for the particles of type specified in the pair vector<const Particle*> type1; vector<const Particle*> type2; // Get all the particles of the type specified in the pair from the particle list - foreach (const Particle& ipart, fs.particles()) { + foreach (const Particle& ipart, inparticles) { // Loop around possible particle pairs (typedef needed to keep BOOST_FOREACH happy) foreach (const PdgIdPair& ipair, _decayids) { if (ipart.pdgId() == ipair.first) { Modified: trunk/src/Projections/LeptonClusters.cc ============================================================================== --- trunk/src/Projections/LeptonClusters.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Projections/LeptonClusters.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -6,11 +6,20 @@ namespace Rivet { - const FinalState& LeptonClusters::constituentsFinalState() const + LeptonClusters::LeptonClusters(const FinalState& photons, const FinalState& signal, + double dRmax, bool cluster, + const std::vector<std::pair<double, double> >& etaRanges, + double pTmin) : + FinalState(etaRanges, pTmin), + _dRmax(dRmax), _cluster(cluster) { - return getProjection<FinalState>("Constituents"); - } + setName("LeptonClusters"); + IdentifiedFinalState photonfs(photons); + photonfs.acceptId(PHOTON); + addProjection(photonfs, "Photons"); + addProjection(signal, "Signal"); + } int LeptonClusters::compare(const Projection& p) const { // Compare the two as final states (for pT and eta cuts) @@ -18,40 +27,17 @@ int fscmp = FinalState::compare(other); if (fscmp != EQUIVALENT) return fscmp; - fscmp = mkNamedPCmp(p, "Constituents"); - return fscmp; - } - - - void LeptonClusters::project(const Event& e) { - _theParticles.clear(); - - const LeptonClustersConstituents& constituents = - applyProjection<LeptonClustersConstituents>(e, "Constituents"); - - foreach (const Particle& p, constituents.clusteredLeptons()) { - if (accept(p)) { - _theParticles.push_back(p); - } - } - } - - - - int LeptonClustersConstituents::compare(const Projection& p) const { - const PCmp fscmp = mkNamedPCmp(p, "Photons"); - if (fscmp != EQUIVALENT) return fscmp; + const PCmp phcmp = mkNamedPCmp(p, "Photons"); + if (phcmp != EQUIVALENT) return phcmp; const PCmp sigcmp = mkNamedPCmp(p, "Signal"); if (sigcmp != EQUIVALENT) return sigcmp; - const LeptonClustersConstituents& other = dynamic_cast<const LeptonClustersConstituents&>(p); - return (cmp(_dRmax, other._dRmax) || cmp(_cluster, other._cluster) || - cmp(_track, other._track)); + return (cmp(_dRmax, other._dRmax) || cmp(_cluster, other._cluster)); } - void LeptonClustersConstituents::project(const Event& e) { + void LeptonClusters::project(const Event& e) { _theParticles.clear(); _clusteredLeptons.clear(); @@ -59,14 +45,14 @@ ParticleVector bareleptons=signal.particles(); if (bareleptons.size()==0) return; + vector<ClusteredLepton> allClusteredLeptons; for (size_t i=0; i<bareleptons.size(); ++i) { - _theParticles.push_back(bareleptons[i]); - _clusteredLeptons.push_back(Particle(bareleptons[i])); + allClusteredLeptons.push_back(ClusteredLepton(bareleptons[i])); } const FinalState& photons = applyProjection<FinalState>(e, "Photons"); - foreach (const Particle& p, photons.particles()) { - const FourMomentum p_P = p.momentum(); + foreach (const Particle& photon, photons.particles()) { + const FourMomentum p_P = photon.momentum(); double dRmin=_dRmax; int idx=-1; for (size_t i=0; i<bareleptons.size(); ++i) { @@ -81,10 +67,18 @@ } } if (idx>-1) { - if (_cluster) _clusteredLeptons[idx].setMomentum(_clusteredLeptons[idx].momentum()+p_P); - if (_track) _theParticles.push_back(p); + if (_cluster) allClusteredLeptons[idx].addPhoton(photon, _cluster); } } - } + foreach (const ClusteredLepton& lepton, allClusteredLeptons) { + if (accept(lepton)) { + _clusteredLeptons.push_back(lepton); + _theParticles.push_back(lepton.constituentLepton()); + _theParticles.insert(_theParticles.end(), + lepton.constituentPhotons().begin(), + lepton.constituentPhotons().end()); + } + } + } } Modified: trunk/src/Projections/WFinder.cc ============================================================================== --- trunk/src/Projections/WFinder.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Projections/WFinder.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -15,79 +15,76 @@ WFinder::WFinder(double etaMin, double etaMax, double pTmin, PdgId pid, - double m2_min, double m2_max, + double minmass, double maxmass, double missingET, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget, + bool useTransverseMass) { vector<pair<double, double> > etaRanges; etaRanges += std::make_pair(etaMin, etaMax); - _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, - dRmax, clusterPhotons, excludePhotonsFromRFS); + _init(etaRanges, pTmin, pid, minmass, maxmass, missingET, + dRmax, clusterPhotons, trackPhotons, masstarget, useTransverseMass); } WFinder::WFinder(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, double m2_max, + double minmass, double maxmass, double missingET, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { - _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, - dRmax, clusterPhotons, excludePhotonsFromRFS); + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget, + bool useTransverseMass) { + _init(etaRanges, pTmin, pid, minmass, maxmass, missingET, + dRmax, clusterPhotons, trackPhotons, masstarget, useTransverseMass); } void WFinder::_init(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, double m2_max, + double minmass, double maxmass, double missingET, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget, + bool useTransverseMass) { setName("WFinder"); + _minmass = minmass; + _maxmass = maxmass; + _masstarget = masstarget; + _pid = pid; + _trackPhotons = trackPhotons; + _useTransverseMass = useTransverseMass; + // Check that the arguments are legal - assert(abs(pid) == ELECTRON || abs(pid) == MUON); - PdgId nu_pid = abs(pid) + 1; - assert(abs(nu_pid) == NU_E || abs(nu_pid) == NU_MU); + assert(abs(_pid) == ELECTRON || abs(_pid) == MUON); + _nu_pid = abs(_pid) + 1; + assert(abs(_nu_pid) == NU_E || abs(_nu_pid) == NU_MU); // Don't make pT or eta cuts on the neutrino - IdentifiedFinalState fs_nu; - fs_nu.acceptNeutrinos(); + IdentifiedFinalState neutrinos; + neutrinos.acceptNeutrinos(); + addProjection(neutrinos, "Neutrinos"); // Lepton clusters FinalState fs; IdentifiedFinalState bareleptons(fs); bareleptons.acceptIdPair(pid); LeptonClusters leptons(fs, bareleptons, dRmax, - clusterPhotons, excludePhotonsFromRFS, + clusterPhotons, etaRanges, pTmin); addProjection(leptons, "LeptonClusters"); - // Make a merged final state projection for charged and neutral leptons - MergedFinalState mergedFS(leptons, fs_nu); - - // Make and register an invariant mass final state for the W decay leptons - vector<pair<PdgId, PdgId> > l_nu_ids; - l_nu_ids += make_pair(abs(pid), -abs(nu_pid)); - l_nu_ids += make_pair(-abs(pid), abs(nu_pid)); - InvMassFinalState imfs(mergedFS, l_nu_ids, m2_min, m2_max, 80.403); - - ///@todo this breaks existing analyses like D0_2008_S7837160 which use mass - /// cuts. have to make this optional. - //imfs.useTransverseMass(); - - addProjection(imfs, "IMFS"); - // Add MissingMomentum proj to calc MET MissingMomentum vismom(fs); addProjection(vismom, "MissingET"); // Set ETmiss _etMiss = missingET; - // FS for non-signal bits of the event VetoedFinalState remainingFS; - remainingFS.addVetoOnThisFinalState(leptons.constituentsFinalState()); - remainingFS.addVetoOnThisFinalState(fs_nu); + remainingFS.addVetoOnThisFinalState(*this); addProjection(remainingFS, "RFS"); } @@ -96,70 +93,52 @@ const FinalState& WFinder::remainingFinalState() const { - /// @todo This must already have been applied to be useful return getProjection<FinalState>("RFS"); } + int WFinder::compare(const Projection& p) const { + PCmp LCcmp = mkNamedPCmp(p, "LeptonClusters"); + if (LCcmp != EQUIVALENT) return LCcmp; - Particle WFinder::constituentLepton() const { - /// @todo Is this safe? This isn't applying the rest of the selection cuts - const InvMassFinalState& imfs = getProjection<InvMassFinalState>("IMFS"); - /// @todo Is this really meant to be _impossible_? - assert(imfs.particles().size()==2); - - Particle p1,p2; - p1 = imfs.particles()[0]; - p2 = imfs.particles()[1]; - if (abs(p1.pdgId()) == ELECTRON || abs(p1.pdgId()) == MUON) { - return p1; - } - else { - return p2; - } - } - - - Particle WFinder::constituentNeutrino() const { - /// @todo Is this safe? This isn't applying the rest of the selection cuts - const InvMassFinalState& imfs = getProjection<InvMassFinalState>("IMFS"); - /// @todo Is this really meant to be _impossible_? - assert(imfs.particles().size()==2); - - Particle p1,p2; - p1 = imfs.particles()[0]; - p2 = imfs.particles()[1]; - if (abs(p1.pdgId()) == ELECTRON || abs(p1.pdgId()) == MUON) { - return p2; - } - else { - return p1; - } - } - - - const FinalState& WFinder::originalLeptonFinalState() const { - const LeptonClusters& leptons=getProjection<LeptonClusters>("LeptonClusters"); - return leptons.constituentsFinalState(); - } - + const WFinder& other = dynamic_cast<const WFinder&>(p); + return (cmp(_minmass, other._minmass) || cmp(_maxmass, other._maxmass) || + cmp(_useTransverseMass, other._useTransverseMass) || + cmp(_etMiss, other._etMiss) || + cmp(_pid, other._pid) || cmp(_trackPhotons, other._trackPhotons)); - int WFinder::compare(const Projection& p) const { - return mkNamedPCmp(p, "IMFS"); } void WFinder::project(const Event& e) { clear(); - const InvMassFinalState& imfs = applyProjection<InvMassFinalState>(e, "IMFS"); - if (imfs.particles().size() != 2) return; + const LeptonClusters& leptons = applyProjection<LeptonClusters>(e, "LeptonClusters"); + const FinalState& neutrinos = applyProjection<FinalState>(e, "Neutrinos"); - /// Apply the RFS here so that it can be acquired later via the remainingFinalState method - applyProjection<FinalState>(e, "RFS"); - - Particle p1,p2; - p1 = imfs.particles()[0]; - p2 = imfs.particles()[1]; + // Make and register an invariant mass final state for the W decay leptons + vector<pair<PdgId, PdgId> > l_nu_ids; + l_nu_ids += make_pair(abs(_pid), -abs(_nu_pid)); + l_nu_ids += make_pair(-abs(_pid), abs(_nu_pid)); + InvMassFinalState imfs(FinalState(), l_nu_ids, _minmass, _maxmass, _masstarget); + imfs.useTransverseMass(_useTransverseMass); + ParticleVector tmp; + tmp.insert(tmp.end(), leptons.clusteredLeptons().begin(), leptons.clusteredLeptons().end()); + tmp.insert(tmp.end(), neutrinos.particles().begin(), neutrinos.particles().end()); + imfs.calc(tmp); + + if (imfs.particlePairs().size() < 1) return; + + ParticlePair Wconstituents(imfs.particlePairs()[0]); + Particle p1(Wconstituents.first), p2(Wconstituents.second); + + if (PID::threeCharge(p1)==0) { + _constituentLeptons += p2; + _constituentNeutrinos += p1; + } + else { + _constituentLeptons += p1; + _constituentNeutrinos += p2; + } FourMomentum pW = p1.momentum() + p2.momentum(); const int w3charge = PID::threeCharge(p1) + PID::threeCharge(p2); @@ -184,7 +163,28 @@ // Make W Particle and insert into particles list const PdgId wpid = (wcharge == 1) ? WPLUSBOSON : WMINUSBOSON; - _theParticles.push_back(Particle(wpid, pW)); + _bosons.push_back(Particle(wpid, pW)); + + // find the LeptonClusters and neutrinos which survived the IMFS cut such that we can + // extract their original particles + foreach (const Particle& p, _constituentNeutrinos) { + _theParticles.push_back(p); + } + foreach (const Particle& p, _constituentLeptons) { + foreach (const ClusteredLepton& l, leptons.clusteredLeptons()) { + if (p.pdgId()==l.pdgId() && p.momentum()==l.momentum()) { + _theParticles.push_back(l.constituentLepton()); + if (_trackPhotons) { + _theParticles.insert(_theParticles.end(), + l.constituentPhotons().begin(), + l.constituentPhotons().end()); + } + } + } + } + +// /// Apply the RFS here so that it can be acquired later via the remainingFinalState method +// applyProjection<FinalState>(e, "RFS"); } Modified: trunk/src/Projections/ZFinder.cc ============================================================================== --- trunk/src/Projections/ZFinder.cc Mon Aug 1 18:00:08 2011 (r3267) +++ trunk/src/Projections/ZFinder.cc Mon Aug 1 18:28:49 2011 (r3268) @@ -1,7 +1,6 @@ // -*- C++ -*- #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/InvMassFinalState.hh" -#include "Rivet/Projections/ClusteredPhotons.hh" #include "Rivet/Projections/LeptonClusters.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Tools/ParticleIdUtils.hh" @@ -14,42 +13,49 @@ ZFinder::ZFinder(double etaMin, double etaMax, double pTmin, PdgId pid, - double m2_min, double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { + double minmass, double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget) { vector<pair<double, double> > etaRanges; etaRanges += std::make_pair(etaMin, etaMax); - _init(etaRanges, pTmin, pid, m2_min, m2_max, dRmax, clusterPhotons, excludePhotonsFromRFS); + _init(etaRanges, pTmin, pid, minmass, maxmass, dRmax, clusterPhotons, trackPhotons, masstarget); } ZFinder::ZFinder(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, const double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { - _init(etaRanges, pTmin, pid, m2_min, m2_max, dRmax, clusterPhotons, excludePhotonsFromRFS); + double minmass, const double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget) { + _init(etaRanges, pTmin, pid, minmass, maxmass, dRmax, clusterPhotons, trackPhotons, masstarget); } void ZFinder::_init(const std::vector<std::pair<double, double> >& etaRanges, double pTmin, PdgId pid, - double m2_min, double m2_max, - double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) + double minmass, double maxmass, + double dRmax, bool clusterPhotons, bool trackPhotons, + double masstarget) { setName("ZFinder"); + _minmass = minmass; + _maxmass = maxmass; + _masstarget = masstarget; + _pid = pid; + _trackPhotons = trackPhotons; + FinalState fs; IdentifiedFinalState bareleptons(fs); bareleptons.acceptIdPair(pid); LeptonClusters leptons(fs, bareleptons, dRmax, - clusterPhotons, excludePhotonsFromRFS, + clusterPhotons, etaRanges, pTmin); addProjection(leptons, "LeptonClusters"); - InvMassFinalState imfs(leptons, std::make_pair(pid, -pid), m2_min, m2_max); - addProjection(imfs, "IMFS"); VetoedFinalState remainingFS; - remainingFS.addVetoOnThisFinalState(leptons.constituentsFinalState()); + remainingFS.addVetoOnThisFinalState(*this); addProjection(remainingFS, "RFS"); } @@ -63,46 +69,59 @@ } - const FinalState& ZFinder::constituentsFinalState() const - { - return getProjection<FinalState>("IMFS"); - } - - const FinalState& ZFinder::originalConstituentsFinalState() const - { - const LeptonClusters& leptons=getProjection<LeptonClusters>("LeptonClusters"); - return leptons.constituentsFinalState(); - } - int ZFinder::compare(const Projection& p) const { - //std::cout<<"Comparing ZFinder"<<std::endl; - PCmp cmp = mkNamedPCmp(p, "IMFS"); - //std::cout<<"Result "<<(int)cmp<<std::endl; - if (cmp != EQUIVALENT) return cmp; + PCmp LCcmp = mkNamedPCmp(p, "LeptonClusters"); + if (LCcmp != EQUIVALENT) return LCcmp; - return EQUIVALENT; + const ZFinder& other = dynamic_cast<const ZFinder&>(p); + return (cmp(_minmass, other._minmass) || cmp(_maxmass, other._maxmass) || + cmp(_pid, other._pid) || cmp(_trackPhotons, other._trackPhotons)); } void ZFinder::project(const Event& e) { - _theParticles.clear(); + clear(); + + const LeptonClusters& leptons = applyProjection<LeptonClusters>(e, "LeptonClusters"); - const FinalState& imfs=applyProjection<FinalState>(e, "IMFS"); - applyProjection<FinalState>(e, "RFS"); - if (imfs.particles().size() != 2) return; - FourMomentum pZ = imfs.particles()[0].momentum() + imfs.particles()[1].momentum(); - const int z3charge = PID::threeCharge(imfs.particles()[0].pdgId()) + PID::threeCharge(imfs.particles()[1].pdgId()); + InvMassFinalState imfs(FinalState(), std::make_pair(_pid, -_pid), _minmass, _maxmass, _masstarget); + ParticleVector tmp; + tmp.insert(tmp.end(), leptons.clusteredLeptons().begin(), leptons.clusteredLeptons().end()); + imfs.calc(tmp); + + if (imfs.particlePairs().size() < 1) return; + ParticlePair Zconstituents(imfs.particlePairs()[0]); + Particle l1(Zconstituents.first), l2(Zconstituents.second); + _constituents += l1, l2; + FourMomentum pZ = l1.momentum() + l2.momentum(); + const int z3charge = PID::threeCharge(l1.pdgId()) + PID::threeCharge(l2.pdgId()); assert(z3charge == 0); stringstream msg; msg << "Z reconstructed from: " << endl - << " " << imfs.particles()[0].momentum() << " " << imfs.particles()[0].pdgId() << endl - << " + " << imfs.particles()[1].momentum() << " " << imfs.particles()[1].pdgId() << endl; + << " " << l1.momentum() << " " << l1.pdgId() << endl + << " + " << l2.momentum() << " " << l2.pdgId() << endl; - Particle Z; - Z.setMomentum(pZ); - _theParticles.push_back(Z); + _bosons.push_back(Particle(ZBOSON, pZ)); getLog() << Log::DEBUG << name() << " found one Z." << endl; + + // find the LeptonClusters which survived the IMFS cut such that we can + // extract their original particles + foreach (const Particle& p, _constituents) { + foreach (const ClusteredLepton& l, leptons.clusteredLeptons()) { + if (p.pdgId()==l.pdgId() && p.momentum()==l.momentum()) { + _theParticles.push_back(l.constituentLepton()); + if (_trackPhotons) { + _theParticles.insert(_theParticles.end(), + l.constituentPhotons().begin(), + l.constituentPhotons().end()); + } + } + } + } + +// /// Apply the RFS here so that it can be acquired later via the remainingFinalState method +// applyProjection<FinalState>(e, "RFS"); }
More information about the Rivet-svn mailing list |