|
[Rivet-svn] r3145 - in trunk: . include/Rivet/Projections src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgFri Jun 10 14:06:25 BST 2011
Author: fsiegert Date: Fri Jun 10 14:06:25 2011 New Revision: 3145 Log: Bugfixes in WFinder. Modified: trunk/ChangeLog trunk/include/Rivet/Projections/InvMassFinalState.hh trunk/include/Rivet/Projections/WFinder.hh trunk/src/Analyses/D0_2008_S7837160.cc trunk/src/Analyses/MC_WJETS.cc trunk/src/Analyses/MC_WPOL.cc trunk/src/Analyses/MC_WWJETS.cc trunk/src/Analyses/Makefile.am trunk/src/Projections/InvMassFinalState.cc trunk/src/Projections/WFinder.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/ChangeLog Fri Jun 10 14:06:25 2011 (r3145) @@ -1,3 +1,7 @@ +2011-06-10 Frank Siegert <frank.siegert at cern.ch> + + * Bugfixes in WFinder. + 2011-06-10 Andy Buckley <andy at insectnation.org> * Adding \physicsxcoor and \physicsycoor treatment to make-plots. Modified: trunk/include/Rivet/Projections/InvMassFinalState.hh ============================================================================== --- trunk/include/Rivet/Projections/InvMassFinalState.hh Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/include/Rivet/Projections/InvMassFinalState.hh Fri Jun 10 14:06:25 2011 (r3145) @@ -15,13 +15,15 @@ InvMassFinalState(const FinalState& fsp, const std::pair<PdgId, PdgId>& idpair, // pair of decay products double minmass, // min inv mass - double maxmass); // max inv mass + double maxmass, // max inv mass + double masstarget=-1.0); InvMassFinalState(const FinalState& fsp, const std::vector<std::pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products double minmass, // min inv mass - double maxmass); // max inv mass + double maxmass, // max inv mass + double masstarget=-1.0); /// Clone on the heap. @@ -59,6 +61,9 @@ /// Max inv mass double _maxmass; + /// Target mass if only one pair should be returned + double _masstarget; + }; Modified: trunk/include/Rivet/Projections/WFinder.hh ============================================================================== --- trunk/include/Rivet/Projections/WFinder.hh Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/include/Rivet/Projections/WFinder.hh Fri Jun 10 14:06:25 2011 (r3145) @@ -22,14 +22,14 @@ /// @name Constructors //@{ - /// Constructor taking a FinalState and type of the charged lepton, mass window, - /// and maximum dR of photons around the charged lepton to take into account for W - /// reconstruction. - WFinder(const ChargedFinalState& fs_l, - PdgId pid, - double m2_min, double m2_max, - double missingET, - double dRmax); + // /// Constructor taking a FinalState and type of the charged lepton, mass window, + // /// and maximum dR of photons around the charged lepton to take into account for W + // /// reconstruction. + // WFinder(const ChargedFinalState& fs_l, + // PdgId pid, + // double m2_min, double m2_max, + // double missingET, + // double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); /// Constructor taking single eta/pT bounds and type of the charged lepton, mass @@ -40,7 +40,7 @@ PdgId pid, double m2_min, double m2_max, double missingET, - double dRmax); + double dRmax, bool clusterPhotons=true, bool excludePhotonsFromRFS=false); /// Constructor taking multiple eta/pT bounds and type of the charged lepton, mass @@ -51,7 +51,7 @@ PdgId pid, double m2_min, const double m2_max, double missingET, - double dRmax); + double dRmax, bool clusterPhotons=true, bool excludePhotonsFromRFS=false); /// Clone on the heap. @@ -66,15 +66,12 @@ /// (e.g. for running a jet finder on it) const FinalState& remainingFinalState() const; - /// Access to the W constituent leptons and photons - const FinalState& constituentsFinalState() const; - - /// Access to the W constituent leptons - const FinalState& constituentLeptonsFinalState() const; - - /// Access to the photons which have been clustered to the leptons - const FinalState& clusteredPhotonsFinalState() 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: @@ -98,14 +95,7 @@ double pTmin, PdgId pid, double m2_min, double m2_max, double missingET, - double dRmax); - - /// Common implementation of constructor operation, taking FS. - void _init(const ChargedFinalState& fs_l, - PdgId pid, - double m2_min, double m2_max, - double missingET, - double dRmax); + double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS); private: Modified: trunk/src/Analyses/D0_2008_S7837160.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7837160.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Analyses/D0_2008_S7837160.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -37,7 +37,7 @@ void init() { // Projections /// @todo Use separate pT and ETmiss cuts in WFinder - const WFinder wfe(-5, 5, 0.0*GeV, ELECTRON, 60.0*GeV, 100.0*GeV, 0*GeV, 0.2); + const WFinder wfe(-5, 5, 25.0*GeV, ELECTRON, 60.0*GeV, 100.0*GeV, 25.0*GeV, 0.2); addProjection(wfe, "WFe"); // Cross-section histograms @@ -62,19 +62,8 @@ // Require that leptons have Et >= 25 GeV /// @todo Use pT cut in WFinder /// @todo Any ETmiss cut? - FourMomentum p_e; - int chg_e = 0; - foreach (const Particle& l, wf.constituentLeptonsFinalState().particles()) { - const FourMomentum pl = l.momentum(); - if (abs(l.pdgId()) == ELECTRON) { - chg_e = PID::threeCharge(l.pdgId()); - p_e = pl; - } - if (pl.Et()/GeV < 25.0) { - getLog() << Log::DEBUG << l.pdgId() << " fails Et cut" << endl; - vetoEvent; - } - } + FourMomentum p_e=wf.constituentLepton().momentum(); + int chg_e = PID::threeCharge(wf.constituentLepton().pdgId()); if (p_e.eta() < 0) chg_e *= -1; assert(chg_e != 0); Modified: trunk/src/Analyses/MC_WJETS.cc ============================================================================== --- trunk/src/Analyses/MC_WJETS.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Analyses/MC_WJETS.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -67,14 +67,13 @@ _h_W_pT_peak->fill(wmom.pT(), weight); _h_W_y->fill(wmom.rapidity(), weight); _h_W_phi->fill(wmom.azimuthalAngle(), weight); - foreach (const Particle& l, wfinder.constituentLeptonsFinalState().particles()) { - _h_lepton_pT->fill(l.momentum().pT(), weight); - _h_lepton_eta->fill(l.momentum().eta(), weight); - if (PID::threeCharge(l.pdgId()) != 0) { - emom = l.momentum(); - charge3_x_eta = PID::threeCharge(l.pdgId()) * emom.eta(); - charge3 = PID::threeCharge(l.pdgId()); - } + Particle l=wfinder.constituentLepton(); + _h_lepton_pT->fill(l.momentum().pT(), weight); + _h_lepton_eta->fill(l.momentum().eta(), weight); + if (PID::threeCharge(l.pdgId()) != 0) { + emom = l.momentum(); + charge3_x_eta = PID::threeCharge(l.pdgId()) * emom.eta(); + charge3 = PID::threeCharge(l.pdgId()); } assert(charge3_x_eta != 0); assert(charge3!=0); Modified: trunk/src/Analyses/MC_WPOL.cc ============================================================================== --- trunk/src/Analyses/MC_WPOL.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Analyses/MC_WPOL.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -78,7 +78,7 @@ const ParticlePair& beams = applyProjection<Beam>(event, "Beams").beams(); FourMomentum pb1(beams.second.momentum()), pb2(beams.first.momentum()); - Particle lepton=wfinder.constituentLeptonsFinalState().particles()[0]; + Particle lepton=wfinder.constituentLepton(); FourMomentum pl(lepton.momentum()); size_t idx = (PID::threeCharge(lepton.pdgId())>0 ? 0 : 1); FourMomentum plnu(wfinder.particles()[0].momentum()); Modified: trunk/src/Analyses/MC_WWJETS.cc ============================================================================== --- trunk/src/Analyses/MC_WWJETS.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Analyses/MC_WWJETS.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -32,10 +32,8 @@ addProjection(wmnufinder, "WmnuFinder"); VetoedFinalState jetinput; jetinput - .addVetoOnThisFinalState(wenufinder.constituentsFinalState()) - .addVetoOnThisFinalState(wenufinder.clusteredPhotonsFinalState()) - .addVetoOnThisFinalState(wmnufinder.constituentsFinalState()) - .addVetoOnThisFinalState(wmnufinder.clusteredPhotonsFinalState()); + .addVetoOnThisFinalState(wenufinder.originalLeptonFinalState()) + .addVetoOnThisFinalState(wmnufinder.originalLeptonFinalState()); FastJets jetpro(jetinput, FastJets::KT, 0.7); addProjection(jetpro, "Jets"); @@ -103,24 +101,10 @@ FourMomentum wmnu(wmnufinder.particles()[0].momentum()); FourMomentum ww(wenu+wmnu); // find leptons - FourMomentum ep(0.0,0.0,0.0,0.0), enu(0.0,0.0,0.0,0.0); - if (PID::threeCharge(wenufinder.constituentsFinalState().particles()[0])>0.0) { - ep=wenufinder.constituentsFinalState().particles()[0].momentum(); - enu=wenufinder.constituentsFinalState().particles()[1].momentum(); - } - else { - ep=wenufinder.constituentsFinalState().particles()[1].momentum(); - enu=wenufinder.constituentsFinalState().particles()[0].momentum(); - } - FourMomentum mnu(0.0,0.0,0.0,0.0), mm(0.0,0.0,0.0,0.0); - if (PID::threeCharge(wmnufinder.constituentsFinalState().particles()[0])==0.0) { - mnu=wmnufinder.constituentsFinalState().particles()[0].momentum(); - mm=wmnufinder.constituentsFinalState().particles()[1].momentum(); - } - else { - mnu=wmnufinder.constituentsFinalState().particles()[1].momentum(); - mm=wmnufinder.constituentsFinalState().particles()[0].momentum(); - } + FourMomentum ep=wenufinder.constituentLepton().momentum(); + FourMomentum enu=wenufinder.constituentNeutrino().momentum(); + FourMomentum mm=wmnufinder.constituentLepton().momentum(); + FourMomentum mnu=wmnufinder.constituentNeutrino().momentum(); _h_WW_pT->fill(ww.pT(),weight); _h_WW_pT_peak->fill(ww.pT(),weight); @@ -168,12 +152,6 @@ foreach (const Jet& jet, jets) { HT+=jet.momentum().pT(); } - foreach (const Particle& p, wenufinder.clusteredPhotonsFinalState().particles()) { - HT+=p.momentum().pT(); - } - foreach (const Particle& p, wmnufinder.clusteredPhotonsFinalState().particles()) { - HT+=p.momentum().pT(); - } if (HT>0.0) _h_HT->fill(HT, weight); if (jets.size()>1) { Modified: trunk/src/Analyses/Makefile.am ============================================================================== --- trunk/src/Analyses/Makefile.am Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Analyses/Makefile.am Fri Jun 10 14:06:25 2011 (r3145) @@ -112,7 +112,6 @@ endif if ENABLE_UNVALIDATED RivetCDFAnalyses_la_SOURCES += \ - CDF_1991_S2313472.cc \ CDF_1993_S2742446.cc \ CDF_1996_S3108457.cc \ CDF_1996_S3349578.cc \ @@ -145,8 +144,7 @@ if ENABLE_UNVALIDATED RivetD0Analyses_la_SOURCES += \ D0_1996_S3214044.cc \ - D0_1996_S3324664.cc \ - D0_1998_S3711838.cc + D0_1996_S3324664.cc endif Modified: trunk/src/Projections/InvMassFinalState.cc ============================================================================== --- trunk/src/Projections/InvMassFinalState.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Projections/InvMassFinalState.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -9,8 +9,9 @@ InvMassFinalState::InvMassFinalState(const FinalState& fsp, const pair<PdgId, PdgId>& idpair, // pair of decay products double minmass, // min inv mass - double maxmass) // max inv mass - : _minmass(minmass), _maxmass(maxmass) + double maxmass, // max inv mass + double masstarget) + : _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget) { setName("InvMassFinalState"); addProjection(fsp, "FS"); @@ -21,8 +22,9 @@ InvMassFinalState::InvMassFinalState(const FinalState& fsp, const vector<pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products double minmass, // min inv mass - double maxmass) // max inv mass - : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass) + double maxmass, // max inv mass + double masstarget) + : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget) { setName("InvMassFinalState"); addProjection(fsp, "FS"); @@ -85,6 +87,7 @@ vector<const Particle*> tmp; // Now calculate the inv mass + pair<double, pair<Particle, Particle> > closestPair; foreach (const Particle* i1, type1) { foreach (const Particle* i2, type2) { // check this is actually a pair @@ -119,9 +122,23 @@ } // Store accepted particle pairs _particlePairs += make_pair(*i1, *i2); + if (_masstarget>0.0) { + double diff=fabs(v4.mass()-_masstarget); + if (diff<closestPair.first) { + closestPair.first=diff; + closestPair.second=make_pair(*i1, *i2); + } + } } } } + if (_masstarget>0.0) { + _theParticles.clear(); + _particlePairs.clear(); + _theParticles += closestPair.second.first; + _theParticles += closestPair.second.second; + _particlePairs += closestPair.second; + } getLog() << Log::DEBUG << "Selected " << _theParticles.size() << " particles " << "(" << _particlePairs.size() << " pairs)" << endl; Modified: trunk/src/Projections/WFinder.cc ============================================================================== --- trunk/src/Projections/WFinder.cc Fri Jun 10 13:03:07 2011 (r3144) +++ trunk/src/Projections/WFinder.cc Fri Jun 10 14:06:25 2011 (r3145) @@ -3,7 +3,7 @@ #include "Rivet/Projections/InvMassFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/MergedFinalState.hh" -#include "Rivet/Projections/ClusteredPhotons.hh" +#include "Rivet/Projections/LeptonClusters.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Tools/ParticleIdUtils.hh" #include "Rivet/Tools/Logging.hh" @@ -12,24 +12,16 @@ namespace Rivet { - WFinder::WFinder(const ChargedFinalState& fs_l, - PdgId pid, - double m2_min, double m2_max, - double missingET, - double dRmax) { - _init(fs_l, pid, m2_min, m2_max, missingET, dRmax); - } - - WFinder::WFinder(double etaMin, double etaMax, double pTmin, PdgId pid, double m2_min, double m2_max, double missingET, - double dRmax) { + double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { vector<pair<double, double> > etaRanges; etaRanges += std::make_pair(etaMin, etaMax); - _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, dRmax); + _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, + dRmax, clusterPhotons, excludePhotonsFromRFS); } @@ -38,8 +30,9 @@ PdgId pid, double m2_min, double m2_max, double missingET, - double dRmax) { - _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, dRmax); + double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { + _init(etaRanges, pTmin, pid, m2_min, m2_max, missingET, + dRmax, clusterPhotons, excludePhotonsFromRFS); } @@ -48,20 +41,13 @@ PdgId pid, double m2_min, double m2_max, double missingET, - double dRmax) { - ChargedFinalState fs_l(etaRanges, pTmin); - _init(fs_l, pid, m2_min, m2_max, missingET, dRmax); - } - - - void WFinder::_init(const ChargedFinalState& fs_l, - PdgId pid, - double m2_min, double m2_max, - double missingET, - double dRmax) + double dRmax, bool clusterPhotons, bool excludePhotonsFromRFS) { setName("WFinder"); + + + // Check that the arguments are legal assert(abs(pid) == ELECTRON || abs(pid) == MUON); PdgId nu_pid = abs(pid) + 1; @@ -71,38 +57,36 @@ IdentifiedFinalState fs_nu; fs_nu.acceptNeutrinos(); - // Make a merged final state projection for charged and neutral leptons - MergedFinalState mergedFS(fs_l, fs_nu); + // lepton clusters + FinalState fs; + IdentifiedFinalState bareleptons(fs); + bareleptons.acceptIdPair(pid); + LeptonClusters leptons(fs, bareleptons, dRmax, + clusterPhotons, excludePhotonsFromRFS, + etaRanges, pTmin); + addProjection(leptons, "LeptonClusters"); - // Mass range - _m2_min = m2_min; - _m2_max = m2_max; - // Set ETmiss - _etMiss = missingET; + // 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); + InvMassFinalState imfs(mergedFS, l_nu_ids, m2_min, m2_max, 80.403); addProjection(imfs, "IMFS"); - // A projection for clustering photons on to the charged lepton - ClusteredPhotons cphotons(FinalState(), imfs, dRmax); - addProjection(cphotons, "CPhotons"); - - // Projection for all signal constituents - MergedFinalState signalFS(imfs, cphotons); - addProjection(signalFS, "SignalParticles"); - // Add MissingMomentum proj to calc MET - MissingMomentum vismom(signalFS); + MissingMomentum vismom(fs); addProjection(vismom, "MissingET"); + // Set ETmiss + _etMiss = missingET; // FS for non-signal bits of the event VetoedFinalState remainingFS; - remainingFS.addVetoOnThisFinalState(signalFS); + remainingFS.addVetoOnThisFinalState(leptons.constituentsFinalState()); + remainingFS.addVetoOnThisFinalState(fs_nu); addProjection(remainingFS, "RFS"); } @@ -115,19 +99,42 @@ } - const FinalState& WFinder::constituentsFinalState() const { - return getProjection<FinalState>("SignalParticles"); + Particle WFinder::constituentLepton() const { + const InvMassFinalState& imfs = getProjection<InvMassFinalState>("IMFS"); + assert(imfs.particles().size()==2); + + Particle p1,p2; + p1 = imfs.particles()[0]; + p2 = imfs.particles()[1]; + if (p1.pdgId() == ELECTRON || p1.pdgId() == MUON) { + return p1; + } + else { + return p2; + } } - const FinalState& WFinder::constituentLeptonsFinalState() const { - return getProjection<FinalState>("IMFS"); + Particle WFinder::constituentNeutrino() const { + const InvMassFinalState& imfs = getProjection<InvMassFinalState>("IMFS"); + assert(imfs.particles().size()==2); + + Particle p1,p2; + p1 = imfs.particles()[0]; + p2 = imfs.particles()[1]; + if (p1.pdgId() == ELECTRON || p1.pdgId() == MUON) { + return p2; + } + else { + return p1; + } } - const FinalState& WFinder::clusteredPhotonsFinalState() const + const FinalState& WFinder::originalLeptonFinalState() const { - return getProjection<FinalState>("CPhotons"); + const LeptonClusters& leptons=getProjection<LeptonClusters>("LeptonClusters"); + return leptons.constituentsFinalState(); } @@ -135,9 +142,6 @@ PCmp cmp = mkNamedPCmp(p, "IMFS"); if (cmp != EQUIVALENT) return cmp; - cmp = mkNamedPCmp(p, "CPhotons"); - if (cmp != EQUIVALENT) return cmp; - return EQUIVALENT; } @@ -149,33 +153,14 @@ void WFinder::project(const Event& e) { clear(); + const InvMassFinalState& imfs = applyProjection<InvMassFinalState>(e, "IMFS"); + applyProjection<FinalState>(e, "RFS"); + if (imfs.particles().size() != 2) return; Particle p1,p2; - if(imfs.particles().size() == 2) { - p1 = imfs.particles()[0]; - p2 = imfs.particles()[1]; - } - else { - // no candiate, return - if(imfs.particles().empty()) { - getLog() << Log::DEBUG << "No W+- candidates found" << " " - << imfs.particles().size() << " " << endl; - return; - } - // more than one, pick the one nearer the W mass - double deltaM = 1e30; - for(unsigned int ix=0;ix<imfs.particlePairs().size();++ix) { - FourMomentum pW = imfs.particlePairs()[ix].first .momentum()+ - imfs.particlePairs()[ix].second.momentum(); - double mW = pW.mass(); - if(fabs(mW-80.403)<deltaM) { - deltaM = fabs(mW-80.4); - p1 = imfs.particlePairs()[ix].first ; - p2 = imfs.particlePairs()[ix].second; - } - } - } + p1 = imfs.particles()[0]; + p2 = imfs.particles()[1]; FourMomentum pW = p1.momentum() + p2.momentum(); const int w3charge = PID::threeCharge(p1) + PID::threeCharge(p2); @@ -189,14 +174,6 @@ << " " << p1.momentum() << " " << p1.pdgId() << endl << " + " << p2.momentum() << " " << p2.pdgId() << endl; - // Add in clustered photons - const FinalState& photons = applyProjection<FinalState>(e, "CPhotons"); - foreach (const Particle& photon, photons.particles()) { - msg << " + " << photon.momentum() << " " << photon.pdgId() << endl; - pW += photon.momentum(); - } - msg << " = " << pW; - // Check missing ET const MissingMomentum& vismom = applyProjection<MissingMomentum>(e, "MissingET"); /// @todo Restrict missing momentum eta range? Use vectorET()? @@ -206,10 +183,6 @@ return; } - // Check mass range again - if (!inRange(pW.mass()/GeV, _m2_min, _m2_max)) return; - getLog() << Log::DEBUG << msg.str() << endl; - // Make W Particle and insert into particles list const PdgId wpid = (wcharge == 1) ? WPLUSBOSON : WMINUSBOSON; _theParticles.push_back(Particle(wpid, pW));
More information about the Rivet-svn mailing list |