|
|
[Rivet-svn] r3225 - trunk/src/Analyses
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue Jul 19 20:10:52 BST 2011
Author: davemallows
Date: Tue Jul 19 20:10:51 2011
New Revision: 3225
Log:
Replaced MC_TTBAR
Modified:
trunk/src/Analyses/MC_TTBAR.cc
Modified: trunk/src/Analyses/MC_TTBAR.cc
==============================================================================
--- trunk/src/Analyses/MC_TTBAR.cc Tue Jul 19 19:58:02 2011 (r3224)
+++ trunk/src/Analyses/MC_TTBAR.cc Tue Jul 19 20:10:51 2011 (r3225)
@@ -1,89 +1,159 @@
-// -*- C++ -*-
#include "Rivet/Analysis.hh"
-#include "Rivet/RivetAIDA.hh"
-#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
-
-namespace Rivet {
+#include "Rivet/Projections/ChargedLeptons.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/AnalysisLoader.hh"
+#include "Rivet/RivetAIDA.hh"
- /// @brief MC validation analysis for Z + jets events
- /// @todo More! This analysis just checks the \f$ \eta \f$ distribution at the moment.
+namespace Rivet {
+
class MC_TTBAR : public Analysis {
- public:
- /// @name Constructors etc.
- //@{
+ public:
- /// Constructor
MC_TTBAR()
: Analysis("MC_TTBAR")
- {
- //setNeedsCrossSection(false);
- }
-
- //@}
-
-
- public:
-
+ { }
+
/// @name Analysis methods
//@{
-
- /// Book histograms and initialise projections before the run
void init() {
+ addProjection(ChargedFinalState(-3.5, 3.5, 0.5*GeV), "CFS");
+ addProjection(ChargedLeptons(ChargedFinalState(-3.5, 3.5, 30*GeV)), "LFS");
+ addProjection(FastJets(FinalState(-2.5, 2.5, 0*GeV), FastJets::KT, 0.5), "JETS");
+
+ _h_jet_0_pT = bookHistogram1D("jet_0_pT", 50, 0, 250);
+ _h_jet_1_pT = bookHistogram1D("jet_1_pT", 50, 0, 250);
+ _h_jet_2_pT = bookHistogram1D("jet_2_pT", 50, 0, 250);
+ _h_jet_3_pT = bookHistogram1D("jet_3_pT", 50, 0, 250);
+
+ _h_bjet_0_pT = bookHistogram1D("bjet_0_pT", 50, 0, 250);
+ _h_bjet_1_pT = bookHistogram1D("bjet_1_pT", 50, 0, 250);
+
+ _h_ljet_0_pT = bookHistogram1D("ljet_0_pT", 50, 0, 250);
+ _h_ljet_1_pT = bookHistogram1D("ljet_1_pT", 50, 0, 250);
+
+ _h_W_mass = bookHistogram1D("W_mass", 75, 30, 180);
+ _h_t_mass = bookHistogram1D("t_mass", 150, 130, 430);
+ _h_t_mass_W_cut = bookHistogram1D("t_mass_W_cut", 150, 130, 430);
+ _h_W_comb_mass = bookHistogram1D("W_comb_mass", 75, 30, 180);
+ _h_t_comb_mass = bookHistogram1D("t_comb_mass", 150, 130, 430);
+ }
+
+ void analyze(const Event& event) {
+ double weight = event.weight();
- const ChargedFinalState cfs(-5.0, 5.0);
- addProjection(cfs, "CFS");
+ const FinalState& cfs = applyProjection<FinalState>(event, "CFS");
+ getLog() << Log::DEBUG << "Total charged multiplicity = "
+ << cfs.size() << endl;
+
+ const ChargedLeptons& lfs = applyProjection<ChargedLeptons>(event, "LFS");
+ getLog() << Log::DEBUG << "Charged lepton multiplicity = "
+ << lfs.chargedLeptons().size() << endl;
+ if (lfs.chargedLeptons().size() != 1) {
+ MSG_DEBUG("Event failed lepton cut");
+ vetoEvent;
+ }
+ foreach (Particle lepton, lfs.chargedLeptons()) {
+ getLog() << Log::DEBUG << "lepton pT = " << lepton.momentum().pT() << endl;
+ }
- /// @todo Book histograms here, e.g.:
- _hist_nch_eta = bookHistogram1D("nch-eta", 20, -5.0, 5.0);
- _hist_nch_pt = bookHistogram1D("nch-pt", 100, 0.0, 200.0);
- _hist_nch_phi = bookHistogram1D("nch-phi", 100, 0.0, TWOPI);
+ const FastJets& jetpro = applyProjection<FastJets>(event, "JETS");
+ const Jets jets = jetpro.jetsByPt();
+ getLog() << Log::DEBUG << "jet multiplicity = " << jets.size() << endl;
+
+ if (jets.size() < 4) {
+ getLog() << Log::DEBUG << "Event failed jet cut" << endl;
+ vetoEvent;
+ }
- }
+ _h_jet_0_pT->fill(jets[0].momentum().pT(), weight);
+ _h_jet_1_pT->fill(jets[1].momentum().pT(), weight);
+ _h_jet_2_pT->fill(jets[2].momentum().pT(), weight);
+ _h_jet_3_pT->fill(jets[3].momentum().pT(), weight);
+
+ if (jets[3].momentum().pT() < 35) {
+ getLog() << Log::DEBUG << "Event failed jet cut" << endl;
+ vetoEvent;
+ }
+ foreach (Jet jet, jets) {
+ getLog() << Log::DEBUG << "jet pT = " << jet.momentum().pT() << endl;
+ }
- /// Perform the per-event analysis
- void analyze(const Event& event) {
- const double weight = event.weight();
- const ChargedFinalState& cfs = applyProjection<ChargedFinalState>(event, "CFS");
+ Jets bjets, ljets;
+ foreach (Jet jet, jets) {
+ if (jet.momentum().pT() < 35*GeV) continue;
+ if (jet.containsBottom())
+ bjets.push_back(jet);
+ else
+ ljets.push_back(jet);
+ }
- foreach (const Particle& p, cfs.particles()) {
- double eta = p.momentum().pseudorapidity();
- double pT = p.momentum().perp();
- double phi = p.momentum().phi();
- _hist_nch_eta->fill(eta, weight);
- _hist_nch_pt->fill(pT, weight);
- _hist_nch_phi->fill(phi, weight);
+ if (bjets.size() !=2) {
+ getLog() << Log::DEBUG << "Event failed b-tagging cut" << endl;
+ vetoEvent;
}
- }
+ _h_bjet_0_pT->fill(bjets[0].momentum().pT(), weight);
+ _h_bjet_1_pT->fill(bjets[1].momentum().pT(), weight);
+ _h_ljet_0_pT->fill(ljets[0].momentum().pT(), weight);
+ _h_ljet_1_pT->fill(ljets[1].momentum().pT(), weight);
- /// Normalise histograms etc., after the run
+ FourMomentum W = ljets[0].momentum() + ljets[1].momentum();
+ FourMomentum t1 = W + bjets[0].momentum();
+ FourMomentum t2 = W + bjets[1].momentum();
+
+ _h_W_mass->fill(mass(W), weight);
+ _h_t_mass->fill(mass(t1), weight);
+ _h_t_mass->fill(mass(t2), weight);
+ if (mass(W) > 70 && mass(W) < 90) {
+ getLog() << Log::INFO << "W found with mass " << W.mass() << endl;
+ _h_t_mass_W_cut->fill(mass(t1), weight);
+ _h_t_mass_W_cut->fill(mass(t2), weight);
+ }
+
+ _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);
+
+ _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);
+ }
+
void finalize() {
- scale(_hist_nch_eta, 1.0/sumOfWeights());
- scale(_hist_nch_pt, 1.0/sumOfWeights());
- scale(_hist_nch_phi, 1.0/sumOfWeights());
+ // No histos, so nothing to do!
}
-
-
- private:
-
- /// @name Histograms
- //@{
- AIDA::IHistogram1D* _hist_nch_eta;
- AIDA::IHistogram1D* _hist_nch_pt;
- AIDA::IHistogram1D* _hist_nch_phi;
//@}
+ private:
+ AIDA::IHistogram1D * _h_jet_0_pT;
+ AIDA::IHistogram1D * _h_jet_1_pT;
+ AIDA::IHistogram1D * _h_jet_2_pT;
+ AIDA::IHistogram1D * _h_jet_3_pT;
+
+ AIDA::IHistogram1D * _h_bjet_0_pT;
+ AIDA::IHistogram1D * _h_bjet_1_pT;
+
+ AIDA::IHistogram1D * _h_ljet_0_pT;
+ AIDA::IHistogram1D * _h_ljet_1_pT;
+
+ AIDA::IHistogram1D * _h_W_mass;
+ AIDA::IHistogram1D * _h_t_mass;
+ AIDA::IHistogram1D * _h_W_comb_mass;
+ AIDA::IHistogram1D * _h_t_comb_mass;
+ AIDA::IHistogram1D * _h_t_mass_W_cut;
};
-
-
// This global object acts as a hook for the plugin system
AnalysisBuilder<MC_TTBAR> plugin_MC_TTBAR;
-
}
More information about the Rivet-svn
mailing list
|