|
[Rivet-svn] r3413 - in branches/2011-07-aida2yoda: . include/Rivet/Projections src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Oct 6 14:16:00 BST 2011
Author: hoeth Date: Thu Oct 6 14:16:00 2011 New Revision: 3413 Log: merge r3346-3347 and r3349 from trunk Modified: branches/2011-07-aida2yoda/ChangeLog branches/2011-07-aida2yoda/include/Rivet/Projections/MissingMomentum.hh branches/2011-07-aida2yoda/src/Analyses/D0_2004_S5992206.cc branches/2011-07-aida2yoda/src/Analyses/MC_SUSY.cc branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc branches/2011-07-aida2yoda/src/Analyses/UA1_1990_S2044935.cc branches/2011-07-aida2yoda/src/Projections/MissingMomentum.cc branches/2011-07-aida2yoda/src/Projections/WFinder.cc Modified: branches/2011-07-aida2yoda/ChangeLog ============================================================================== --- branches/2011-07-aida2yoda/ChangeLog Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/ChangeLog Thu Oct 6 14:16:00 2011 (r3413) @@ -1,3 +1,35 @@ +2011-09-10 Andy Buckley <andy at insectnation.org> + + * Changing the behaviour and interface of the MissingMomentum + projection to calculate vector ET correctly. This was previously + calculated according to the common definition of -E*sin(theta) of + the summed visible 4-momentum in the event, but that is incorrect + because the timelike term grows monotonically. Instead, transverse + 2-vectors of size ET need to be constructed for each visible + particle, and vector-summed in the transverse plane. + + The rewrite of this behaviour made it opportune to make an API + improvement: the previous method names scalarET/vectorET() have + been renamed to scalar/vectorEt() to better match the Rivet + FourMomentum::Et() method, and MissingMomentum::vectorEt() now + returns a Vector3 rather than a double so that the transverse + missing Et direction is also available. + + Only one data analysis has been affected by this change in + behaviour: the D0_2004_S5992206 dijet delta(phi) analysis. It's + expected that this change will not be very significant, as it is + a *veto* on significant missing ET to reduce non-QCD + contributions. MC studies using this analysis ~always run with QCD + events only, so these contributions should be small. The analysis + efficiency may have been greatly improved, as fewer events will + now fail the missing ET veto cut. + + * Add sorting of the ParticleVector returned by the ChargedLeptons + projection. + + * configure.ac: Adding a check to make sure that no-one tries to + install into --prefix=$PWD. + 2011-09-04 Andy Buckley <andy at insectnation.org> * lighthisto fixes from Christian Roehr. Modified: branches/2011-07-aida2yoda/include/Rivet/Projections/MissingMomentum.hh ============================================================================== --- branches/2011-07-aida2yoda/include/Rivet/Projections/MissingMomentum.hh Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/include/Rivet/Projections/MissingMomentum.hh Thu Oct 6 14:16:00 2011 (r3413) @@ -47,16 +47,13 @@ public: /// The vector-summed visible four-momentum in the event. - FourMomentum& visibleMomentum() { return _momentum; } - - /// The vector-summed visible four-momentum in the event. const FourMomentum& visibleMomentum() const { return _momentum; } /// The vector-summed (in)visible transverse energy in the event - double vectorET() const { return _momentum.Et(); } + const Vector3& vectorEt() const { return _vet; } - /// The scalar-summed (in)visible transverse energy in the event. - double scalarET() const { return _set; } + /// The scalar-summed visible transverse energy in the event. + double scalarEt() const { return _set; } protected: @@ -82,6 +79,9 @@ /// Scalar transverse energy double _set; + /// Vector transverse energy + Vector3 _vet; + }; Modified: branches/2011-07-aida2yoda/src/Analyses/D0_2004_S5992206.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/D0_2004_S5992206.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Analyses/D0_2004_S5992206.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -66,11 +66,11 @@ // Analyse and print some info const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); - getLog() << Log::DEBUG << "Jet multiplicity before any pT cut = " << jetpro.size() << endl; + MSG_DEBUG("Jet multiplicity before any pT cut = " << jetpro.size()); const Jets jets = jetpro.jetsByPt(40.0*GeV); if (jets.size() >= 2) { - getLog() << Log::DEBUG << "Jet multiplicity after pT > 40 GeV cut = " << jets.size() << endl; + MSG_DEBUG("Jet multiplicity after pT > 40 GeV cut = " << jets.size()); } else { vetoEvent; } @@ -79,14 +79,14 @@ if (fabs(rap1) > 0.5 || fabs(rap2) > 0.5) { vetoEvent; } - getLog() << Log::DEBUG << "Jet eta and pT requirements fulfilled" << endl; + MSG_DEBUG("Jet eta and pT requirements fulfilled"); const double pT1 = jets[0].momentum().pT(); const MissingMomentum& caloMissEt = applyProjection<MissingMomentum>(event, "CalMET"); - getLog() << Log::DEBUG << "Missing vector Et = " << caloMissEt.vectorET()/GeV << " GeV" << endl; - if (caloMissEt.vectorET() > 0.7*pT1) { - MSG_DEBUG("Vetoing event with too much missing Et: " - << caloMissEt.vectorET()/GeV << " GeV > " + MSG_DEBUG("Missing vector Et = " << caloMissEt.vectorEt()/GeV << " GeV"); + if (caloMissEt.vectorEt().mod() > 0.7*pT1) { + MSG_DEBUG("Vetoing event with too much missing ET: " + << caloMissEt.vectorEt()/GeV << " GeV > " << 0.7*pT1/GeV << " GeV"); vetoEvent; } Modified: branches/2011-07-aida2yoda/src/Analyses/MC_SUSY.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/MC_SUSY.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Analyses/MC_SUSY.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -204,7 +204,7 @@ // Calculate and fill missing Et histos const MissingMomentum& met = applyProjection<MissingMomentum>(evt, "MET"); - _hist_met->fill(met.vectorET()/GeV); + _hist_met->fill(met.vectorEt().mod()/GeV); // Choose highest-pT leptons of each sign and flavour for dilepton mass edges const FinalState& lpfs = applyProjection<FinalState>(evt, "LeadingParticles"); Modified: branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -1,54 +1,71 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedLeptons.hh" +#include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/RivetYODA.hh" - namespace Rivet { - class MC_TTBAR : public Analysis { + class MC_TTBAR : public Analysis { public: - MC_TTBAR() - : Analysis("MC_TTBAR") - { } + /// Minimal constructor + MC_TTBAR() : Analysis("MC_TTBAR") + { + _sumwPassedLepJetMET = 0; + _sumwPassedJetID = 0; + _sumwPassedWMass = 0; + } /// @name Analysis methods //@{ + /// Set up projections and book histograms void init() { + + // A FinalState is used to select particles within |eta| < 4.2 and with pT + // > 30 GeV, out of which the ChargedLeptons projection picks only the + // electrons and muons, to be accessed later as "LFS". addProjection(ChargedLeptons(FinalState(-4.2, 4.2, 30*GeV)), "LFS"); - addProjection(FastJets(FinalState(-4.2, 4.2, 0*GeV), FastJets::ANTIKT, 0.4), "Jets"); + // A second FinalState is used to select all particles in |eta| < 4.2, + // with no pT cut. This is used to construct jets and measure missing + // transverse energy. + FinalState fs(-4.2, 4.2, 0*GeV); + addProjection(FastJets(fs, FastJets::ANTIKT, 0.6), "Jets"); + addProjection(MissingMomentum(fs), "MissingET"); + // Booking of histograms _h_jet_1_pT = bookHisto1D("jet_1_pT", 50, 0, 500); _h_jet_2_pT = bookHisto1D("jet_2_pT", 50, 0, 400); _h_jet_3_pT = bookHisto1D("jet_3_pT", 50, 0, 300); _h_jet_4_pT = bookHisto1D("jet_4_pT", 50, 0, 200); - - _h_bjet_1_pT = bookHisto1D("jetb_1_pT", 50, 0, 250); - _h_bjet_2_pT = bookHisto1D("jetb_2_pT", 50, 0, 250); - - _h_ljet_1_pT = bookHisto1D("jetl_1_pT", 50, 0, 250); - _h_ljet_2_pT = bookHisto1D("jetl_2_pT", 50, 0, 250); - + _h_jet_HT = bookHisto1D("jet_HT", logspace(0, 2000, 50)); + // + _h_bjet_1_pT = bookHisto1D("jetb_1_pT", 50, 0, 400); + _h_bjet_2_pT = bookHisto1D("jetb_2_pT", 50, 0, 300); + // + _h_ljet_1_pT = bookHisto1D("jetl_1_pT", 50, 0, 400); + _h_ljet_2_pT = bookHisto1D("jetl_2_pT", 50, 0, 300); + // _h_W_mass = bookHisto1D("W_mass", 75, 30, 180); _h_t_mass = bookHisto1D("t_mass", 150, 130, 430); _h_t_mass_W_cut = bookHisto1D("t_mass_W_cut", 150, 130, 430); - _h_W_comb_mass = bookHisto1D("W_comb_mass", 75, 30, 180); - _h_t_comb_mass = bookHisto1D("t_comb_mass", 150, 130, 430); } void analyze(const Event& event) { const double weight = event.weight(); + // Use the "LFS" projection to require at least one hard charged + // lepton. This is an experimental signature for the leptonically decaying + // W. This helps to reduce pure QCD backgrounds. const ChargedLeptons& lfs = applyProjection<ChargedLeptons>(event, "LFS"); MSG_DEBUG("Charged lepton multiplicity = " << lfs.chargedLeptons().size()); - foreach (Particle lepton, lfs.chargedLeptons()) { + foreach (const Particle& lepton, lfs.chargedLeptons()) { MSG_DEBUG("Lepton pT = " << lepton.momentum().pT()); } if (lfs.chargedLeptons().empty()) { @@ -56,28 +73,65 @@ vetoEvent; } + // Use a missing ET cut to bias toward events with a hard neutrino from + // the leptonically decaying W. This helps to reduce pure QCD backgrounds. + const MissingMomentum& met = applyProjection<MissingMomentum>(event, "MissingET"); + MSG_DEBUG("Vector ET = " << met.vectorEt().mod() << " GeV"); + if (met.vectorEt().mod() < 30*GeV) { + MSG_DEBUG("Event failed missing ET cut"); + vetoEvent; + } + + // Use the "Jets" projection to check that there are at least 4 jets of + // any pT. Getting the jets sorted by pT ensures that the first jet is the + // hardest, and so on. We apply no pT cut here only because we want to + // plot all jet pTs to help optimise our jet pT cut. const FastJets& jetpro = applyProjection<FastJets>(event, "Jets"); const Jets alljets = jetpro.jetsByPt(); if (alljets.size() < 4) { MSG_DEBUG("Event failed jet multiplicity cut"); vetoEvent; } - _h_jet_1_pT->fill(alljets[0].momentum().pT(), weight); - _h_jet_2_pT->fill(alljets[1].momentum().pT(), weight); - _h_jet_3_pT->fill(alljets[2].momentum().pT(), weight); - _h_jet_4_pT->fill(alljets[3].momentum().pT(), weight); - const Jets jets = jetpro.jetsByPt(35*GeV); - foreach (const Jet& jet, jets) { - MSG_DEBUG("Jet pT = " << jet.momentum().pT()/GeV << " GeV"); - } - if (jets.size() < 4) { - MSG_DEBUG("Event failed jet pT cut"); + // Update passed-cuts counter and fill all-jets histograms + _sumwPassedLepJetMET += weight; + _h_jet_1_pT->fill(alljets[0].momentum().pT()/GeV, weight); + _h_jet_2_pT->fill(alljets[1].momentum().pT()/GeV, weight); + _h_jet_3_pT->fill(alljets[2].momentum().pT()/GeV, weight); + _h_jet_4_pT->fill(alljets[3].momentum().pT()/GeV, weight); + + // Insist that the hardest 4 jets pass pT hardness cuts. If we don't find + // at least 4 such jets, we abandon this event. + const Jets jets = jetpro.jetsByPt(30*GeV); + double ht = 0.0; + foreach (const Jet& j, jets) { ht += j.momentum().pT(); } + _h_jet_HT->fill(ht/GeV, weight); + if (jets.size() < 4 || + jets[0].momentum().pT() < 60*GeV || + jets[1].momentum().pT() < 50*GeV || + jets[3].momentum().pT() < 30*GeV) { + MSG_DEBUG("Event failed jet cuts"); vetoEvent; } + // Sort the jets into b-jets and light jets. We expect one hard b-jet from + // each top decay, so our 4 hardest jets should include two b-jets. The + // Jet::containsBottom() method is equivalent to perfect experimental + // b-tagging, in a generator-independent way. Jets bjets, ljets; foreach (const Jet& jet, jets) { + // // Don't count jets that overlap with the hard leptons + bool isolated = true; + foreach (const Particle& lepton, lfs.chargedLeptons()) { + if (deltaR(jet.momentum(), lepton.momentum()) < 0.3) { + isolated = false; + break; + } + } + if (!isolated) { + MSG_DEBUG("Jet failed lepton isolation cut"); + break; + } if (jet.containsBottom()) { bjets.push_back(jet); } else { @@ -86,47 +140,72 @@ } MSG_DEBUG("Number of b-jets = " << bjets.size()); if (bjets.size() != 2) { - MSG_DEBUG("Event failed b-tagging cut"); + MSG_DEBUG("Event failed post-lepton-isolation b-tagging cut"); vetoEvent; } + if (ljets.size() < 2) { + MSG_DEBUG("Event failed since not enough light jets remaining after lepton-isolation"); + vetoEvent; + } + + // Plot the pTs of the identified jets. + _sumwPassedJetID += weight; _h_bjet_1_pT->fill(bjets[0].momentum().pT(), weight); _h_bjet_2_pT->fill(bjets[1].momentum().pT(), weight); _h_ljet_1_pT->fill(ljets[0].momentum().pT(), weight); _h_ljet_2_pT->fill(ljets[1].momentum().pT(), weight); - const FourMomentum W = ljets[0].momentum() + ljets[1].momentum(); + // Construct the hadronically decaying W momentum 4-vector from pairs of + // non-b-tagged jets. The pair which best matches the W mass is used. We start + // with an always terrible 4-vector estimate which should always be "beaten" by + // a real jet pair. + FourMomentum W(10*sqrtS(), 0, 0, 0); + for (size_t i = 0; i < ljets.size()-1; ++i) { + for (size_t j = i + 1; j < ljets.size(); ++j) { + const FourMomentum Wcand = ljets[i].momentum() + ljets[j].momentum(); + MSG_TRACE(i << "," << j << ": candidate W mass = " << Wcand.mass()/GeV + << " GeV, vs. incumbent candidate with " << W.mass()/GeV << " GeV"); + if (fabs(Wcand.mass() - 80.4*GeV) < fabs(W.mass() - 80.4*GeV)) { + W = Wcand; + } + } + } + MSG_DEBUG("Candidate W mass = " << W.mass() << " GeV"); + + // There are two b-jets with which this can be combined to make the + // hadronically decaying top, one of which is correct and the other is + // not... but we have no way to identify which is which, so we construct + // both possible top momenta and fill the histograms with both. const FourMomentum t1 = W + bjets[0].momentum(); const FourMomentum t2 = W + bjets[1].momentum(); - _h_W_mass->fill(W.mass(), weight); _h_t_mass->fill(t1.mass(), weight); _h_t_mass->fill(t2.mass(), weight); - if (inRange(W.mass()/GeV, 70, 90)) { + + // Placing a cut on the well-known W mass helps to reduce backgrounds + if (inRange(W.mass()/GeV, 75, 85)) { MSG_DEBUG("W found with mass " << W.mass()/GeV << " GeV"); + _sumwPassedWMass += weight; _h_t_mass_W_cut->fill(t1.mass(), weight); _h_t_mass_W_cut->fill(t2.mass(), weight); } - // All combinatoric 2-jet masses - _h_W_comb_mass->fill(mass(jets[0].momentum() + jets[1].momentum()), weight); - _h_W_comb_mass->fill(mass(jets[0].momentum() + jets[2].momentum()), weight); - _h_W_comb_mass->fill(mass(jets[0].momentum() + jets[3].momentum()), weight); - _h_W_comb_mass->fill(mass(jets[1].momentum() + jets[2].momentum()), weight); - _h_W_comb_mass->fill(mass(jets[1].momentum() + jets[3].momentum()), weight); - _h_W_comb_mass->fill(mass(jets[2].momentum() + jets[3].momentum()), weight); - - // All combinatoric 3-jet masses - _h_t_comb_mass->fill(mass(jets[0].momentum() + jets[1].momentum() + jets[2].momentum()), weight); - _h_t_comb_mass->fill(mass(jets[0].momentum() + jets[1].momentum() + jets[3].momentum()), weight); - _h_t_comb_mass->fill(mass(jets[0].momentum() + jets[2].momentum() + jets[3].momentum()), weight); - _h_t_comb_mass->fill(mass(jets[1].momentum() + jets[2].momentum() + jets[3].momentum()), weight); - - /// @todo Add reconstruction of the other top from the leptonically decaying W, using WFinder } void finalize() { - // No histos, so nothing to do! + scale(_h_jet_1_pT, 1/_sumwPassedLepJetMET); + scale(_h_jet_2_pT, 1/_sumwPassedLepJetMET); + scale(_h_jet_3_pT, 1/_sumwPassedLepJetMET); + scale(_h_jet_4_pT, 1/_sumwPassedLepJetMET); + scale(_h_jet_HT, 1/_sumwPassedLepJetMET); + scale(_h_bjet_1_pT, 1/_sumwPassedJetID); + scale(_h_bjet_2_pT, 1/_sumwPassedJetID); + scale(_h_ljet_1_pT, 1/_sumwPassedJetID); + scale(_h_ljet_2_pT, 1/_sumwPassedJetID); + scale(_h_W_mass, 1/_sumwPassedJetID); + scale(_h_t_mass, 1/_sumwPassedJetID); + scale(_h_t_mass_W_cut, 1/_sumwPassedWMass); } //@} @@ -134,14 +213,20 @@ private: + // Passed-cuts counters + double _sumwPassedLepJetMET, _sumwPassedJetID, _sumwPassedWMass; + + // @name Histogram data members + //@{ + Histo1DPtr _h_jet_1_pT, _h_jet_2_pT, _h_jet_3_pT, _h_jet_4_pT; + Histo1DPtr _h_jet_HT; Histo1DPtr _h_bjet_1_pT, _h_bjet_2_pT; Histo1DPtr _h_ljet_1_pT, _h_ljet_2_pT; Histo1DPtr _h_W_mass; - Histo1DPtr _h_t_mass; - Histo1DPtr _h_W_comb_mass; - Histo1DPtr _h_t_comb_mass; - Histo1DPtr _h_t_mass_W_cut; + Histo1DPtr _h_t_mass, _h_t_mass_W_cut; + + //@} }; Modified: branches/2011-07-aida2yoda/src/Analyses/UA1_1990_S2044935.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/UA1_1990_S2044935.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Analyses/UA1_1990_S2044935.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -78,8 +78,8 @@ // Use good central detector tracks const FinalState& cfs = applyProjection<FinalState>(event, "TrackFS"); - const double Et25 = applyProjection<MissingMomentum>(event, "MET25").scalarET(); - const double Et60 = applyProjection<MissingMomentum>(event, "MET60").scalarET(); + const double Et25 = applyProjection<MissingMomentum>(event, "MET25").scalarEt(); + const double Et60 = applyProjection<MissingMomentum>(event, "MET60").scalarEt(); const unsigned int nch = cfs.size(); // Event level histos Modified: branches/2011-07-aida2yoda/src/Projections/MissingMomentum.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Projections/MissingMomentum.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Projections/MissingMomentum.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -14,18 +14,20 @@ void MissingMomentum::clear() { _momentum = FourMomentum(); _set = 0.0; + _vet = Vector3(); } void MissingMomentum::project(const Event& e) { clear(); - + // Project into final state const FinalState& vfs = applyProjection<FinalState>(e, "VisibleFS"); foreach (const Particle& p, vfs.particles()) { const FourMomentum& mom = p.momentum(); _momentum += mom; _set += mom.Et(); + _vet += mom.Et() * mom.vector3().setZ(0.0).unit(); } } Modified: branches/2011-07-aida2yoda/src/Projections/WFinder.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Projections/WFinder.cc Thu Oct 6 14:03:02 2011 (r3412) +++ branches/2011-07-aida2yoda/src/Projections/WFinder.cc Thu Oct 6 14:16:00 2011 (r3413) @@ -155,9 +155,9 @@ // Check missing ET const MissingMomentum& vismom = applyProjection<MissingMomentum>(e, "MissingET"); /// @todo Restrict missing momentum eta range? Use vectorET()? - if (vismom.scalarET() < _etMiss) { - getLog() << Log::DEBUG << "Not enough missing ET: " << vismom.scalarET()/GeV - << " GeV vs. " << _etMiss/GeV << " GeV" << endl; + if (vismom.scalarEt() < _etMiss) { + MSG_DEBUG("Not enough missing ET: " << vismom.scalarEt()/GeV + << " GeV vs. " << _etMiss/GeV << " GeV"); return; }
More information about the Rivet-svn mailing list |