[Rivet-svn] r1844 - in trunk: . include/Rivet/Projections src/Analyses src/Projections

blackhole at projects.hepforge.org blackhole at projects.hepforge.org
Thu Oct 1 18:21:30 BST 2009


Author: buckley
Date: Thu Oct  1 18:21:30 2009
New Revision: 1844

Log:
Tidying up, adding sorted particle accessors to FinalState cf. JetAlg (sorry, header change forces big rebuild) and adding an untested SUSY validation analysis like the one we've been using in Atlas for Herwig++ vs. HERWIG validation.

Added:
   trunk/src/Analyses/MC_LHC_SUSY.cc
      - copied, changed from r1836, trunk/src/Analyses/MC_LHC_LEADINGJETS.cc
Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Projections/ChargedLeptons.hh
   trunk/include/Rivet/Projections/FinalState.hh
   trunk/include/Rivet/Projections/JetAlg.hh
   trunk/include/Rivet/Projections/LeadingParticlesFinalState.hh
   trunk/src/Analyses/MC_LHC_DIJET.cc
   trunk/src/Analyses/MC_LHC_LEADINGJETS.cc
   trunk/src/Analyses/MC_LHC_WANALYSIS.cc
   trunk/src/Analyses/MC_LHC_ZANALYSIS.cc
   trunk/src/Analyses/Makefile.am
   trunk/src/Projections/ChargedLeptons.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/ChangeLog	Thu Oct  1 18:21:30 2009	(r1844)
@@ -1,3 +1,13 @@
+2009-10-01  Andy Buckley  <andy.buckley at cern.ch>
+
+	* Moving projection setup and registration into init() method for
+	analyses from ALEPH, CDF and the MC_ group.
+
+	* Adding generic SUSY validation analysis, based on plots used in
+	ATLAS Herwig++ validation.
+
+	* Adding sorted particle accessors to FinalState (cf. JetAlg).
+
 2009-09-29  Andy Buckley  <andy at insectnation.org>
 
 	* Adding optional use of args as regex match expressions with

Modified: trunk/include/Rivet/Projections/ChargedLeptons.hh
==============================================================================
--- trunk/include/Rivet/Projections/ChargedLeptons.hh	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/include/Rivet/Projections/ChargedLeptons.hh	Thu Oct  1 18:21:30 2009	(r1844)
@@ -3,7 +3,7 @@
 #define RIVET_ChargedLeptons_HH
 
 #include "Rivet/Projection.hh"
-#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 
@@ -20,7 +20,7 @@
     ChargedLeptons(const FinalState& fsp)
     { 
       setName("ChargedLeptons");
-      addProjection(fsp, "FS");
+      addProjection(ChargedFinalState(fsp), "ChFS");
     }
     
     /// Clone on the heap.
@@ -31,7 +31,7 @@
   protected:
     
     /// Apply the projection to the event.
-    void project(const Event& e);
+    void project(const Event& evt);
     
     /// Compare projections.
     int compare(const Projection& other) const;

Modified: trunk/include/Rivet/Projections/FinalState.hh
==============================================================================
--- trunk/include/Rivet/Projections/FinalState.hh	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/include/Rivet/Projections/FinalState.hh	Thu Oct  1 18:21:30 2009	(r1844)
@@ -35,9 +35,32 @@
     //@}
 
 
-    /// Access the projected final-state particles.
+    /// Get the final-state particles.
     virtual const ParticleVector& particles() const { return _theParticles; }
 
+    /// Get the final-state particles, ordered by supplied sorting function object.
+    template <typename F>
+    const ParticleVector& particles(F sorter) const {
+      std::sort(_theParticles.begin(), _theParticles.end(), sorter);
+      return _theParticles;
+    }
+
+    /// Get the final-state particles, ordered by \f$ p_T \f$.
+    const ParticleVector& particlesByPt() const {
+      return particles(cmpParticleByPt);
+    }
+
+    /// Get the final-state particles, ordered by \f$ E \f$.
+    const ParticleVector& particlesByE() const {
+      return particles(cmpParticleByE);
+    }
+
+    /// Get the final-state particles, ordered by \f$ E_T \f$.
+    const ParticleVector& particlesByEt() const {
+      return particles(cmpParticleByEt);
+    }
+
+
     /// Access the projected final-state particles.
     virtual const size_t size() const { return _theParticles.size(); }
 
@@ -77,7 +100,7 @@
     double _ptmin;
     
     /// The final-state particles.
-    ParticleVector _theParticles;
+    mutable ParticleVector _theParticles;
     
   };
 

Modified: trunk/include/Rivet/Projections/JetAlg.hh
==============================================================================
--- trunk/include/Rivet/Projections/JetAlg.hh	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/include/Rivet/Projections/JetAlg.hh	Thu Oct  1 18:21:30 2009	(r1844)
@@ -73,7 +73,7 @@
 
     /// Get the jets, ordered by supplied sorting function object.
     template <typename F>
-      Jets jets(F sorter, double ptmin=0.0) const {
+    Jets jets(F sorter, double ptmin=0.0) const {
       Jets js = jets(ptmin);
       if (sorter != 0) {
         std::sort(js.begin(), js.end(), sorter);

Modified: trunk/include/Rivet/Projections/LeadingParticlesFinalState.hh
==============================================================================
--- trunk/include/Rivet/Projections/LeadingParticlesFinalState.hh	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/include/Rivet/Projections/LeadingParticlesFinalState.hh	Thu Oct  1 18:21:30 2009	(r1844)
@@ -17,8 +17,9 @@
 
     /// Constructor: the supplied FinalState projection is assumed to live through the run.
     /// @todo Why specify the rap & pT cuts again?
-    LeadingParticlesFinalState(FinalState & fsp, double mineta = -MAXRAPIDITY, double maxeta = MAXRAPIDITY, double minpt = 0.0 * GeV)
-  :  FinalState(mineta, maxeta, minpt) {
+    LeadingParticlesFinalState(const FinalState& fsp, double mineta=-MAXRAPIDITY, double maxeta=MAXRAPIDITY, double minpt=0.0*GeV)
+      :  FinalState(mineta, maxeta, minpt) 
+    {
       setName("LeadingParticlesFinalState");
       addProjection(fsp, "FS");
     }
@@ -40,22 +41,29 @@
       _ids.insert(-id);
       return *this;
     } 
+
+    // /// Check if a particle of a particular ID was found in the current event
+    // bool hasParticleId(const PdgId pid) const;
+
+    // /// Get a particle of a particular ID (check it exists first)
+    // bool get(const PdgId pid) const;
+
     
   protected:
 
     /// Apply the projection on the supplied event.
-    void project(const Event & e);
+    void project(const Event& e);
 
     /// Compare projections.
-    int compare(const Projection & p) const;
+    int compare(const Projection& p) const;
 
-    /// check if the particle's id is in the list 
-    bool inList(const Particle & particle) const;
+    /// Check if the particle's ID is in the list
+    bool inList(const Particle& particle) const;
 
   private:
 
     /// IDs of the leading particles to be selected
-    std::set < long >_ids;
+    std::set<long>_ids;
 
   };
 

Modified: trunk/src/Analyses/MC_LHC_DIJET.cc
==============================================================================
--- trunk/src/Analyses/MC_LHC_DIJET.cc	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Analyses/MC_LHC_DIJET.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -14,12 +14,6 @@
     /// Default constructor
     MC_LHC_DIJET() 
       : Analysis("MC_LHC_DIJET") {
-      FinalState fs;
-      ChargedFinalState cfs;
-      addProjection(fs, "FS");
-      addProjection(cfs, "CFS");
-      addProjection(FastJets(fs, FastJets::KT, 0.7), "Jets");
-      addProjection(FastJets(cfs, FastJets::KT, 0.7), "ChargedJets");
     }
     
     
@@ -27,6 +21,13 @@
     //@{
 
     void init() { 
+      FinalState fs(-4, 4, 0.5*GeV);
+      ChargedFinalState cfs(fs);
+      addProjection(fs, "FS");
+      addProjection(cfs, "CFS");
+      addProjection(FastJets(fs, FastJets::ANTIKT, 0.7), "Jets");
+      addProjection(FastJets(cfs, FastJets::ANTIKT, 0.7), "ChargedJets");
+
       _hist_jetcount = bookHistogram1D("d01-x01-y01", 5, 0., 10.);
       _hist_jetpt = bookHistogram1D("d02-x01-y01", 30, 30.,100.);
       _hist_jetptlog = bookHistogram1D("d03-x01-y01", 20, 0.,8.);

Modified: trunk/src/Analyses/MC_LHC_LEADINGJETS.cc
==============================================================================
--- trunk/src/Analyses/MC_LHC_LEADINGJETS.cc	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Analyses/MC_LHC_LEADINGJETS.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -20,15 +20,6 @@
       : Analysis("MC_LHC_LEADINGJETS")
     { 
       setBeams(PROTON, PROTON);
-      
-      // Final state for the jet finding
-      const FinalState fsj(-4.0, 4.0, 0.0*GeV);
-      addProjection(fsj, "FSJ");
-      addProjection(FastJets(fsj, FastJets::KT, 0.7), "Jets");
-      
-      // Charged final state for the distributions
-      const ChargedFinalState cfs(-1.0, 1.0, 0.5*GeV);
-      addProjection(cfs, "CFS");
     }
     
     
@@ -37,6 +28,15 @@
     
     // Book histograms
     void init() {
+      // Final state for the jet finding
+      const FinalState fsj(-4.0, 4.0, 0.0*GeV);
+      addProjection(fsj, "FSJ");
+      addProjection(FastJets(fsj, FastJets::KT, 0.7), "Jets");
+      
+      // Charged final state for the distributions
+      const ChargedFinalState cfs(-1.0, 1.0, 0.5*GeV);
+      addProjection(cfs, "CFS");
+
       const double maxpt1 = 500.0/GeV;
       _hist_pnchg      = bookProfile1D("trans-nchg", 50, 0.0, maxpt1);
       _hist_pmaxnchg   = bookProfile1D("trans-maxnchg", 50, 0.0, maxpt1);
@@ -61,7 +61,7 @@
       const Jets jets = jetpro.jetsByPt();
       getLog() << Log::DEBUG << "Jet multiplicity = " << jets.size() << endl;
 
-      // We require the leading jet to be within |eta| < 2
+      // Require the leading jet to be within |eta| < 2
       if (jets.size() < 1 || fabs(jets[0].momentum().pseudorapidity()) > 2) {
         getLog() << Log::DEBUG << "Failed jet cut" << endl;
         vetoEvent;

Copied and modified: trunk/src/Analyses/MC_LHC_SUSY.cc (from r1836, trunk/src/Analyses/MC_LHC_LEADINGJETS.cc)
==============================================================================
--- trunk/src/Analyses/MC_LHC_LEADINGJETS.cc	Thu Sep 24 16:04:52 2009	(r1836, copy source)
+++ trunk/src/Analyses/MC_LHC_SUSY.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -5,30 +5,24 @@
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
+#include "Rivet/Projections/IdentifiedFinalState.hh"
+#include "Rivet/Projections/TotalVisibleMomentum.hh"
+#include "Rivet/Projections/LeadingParticlesFinalState.hh"
 
 namespace Rivet {
 
 
-  /* Underlying event in leading jet, extended to the LHC
+  /* Basic SUSY type validation analysis for the LHC
    * @author Andy Buckley
    */ 
-  class MC_LHC_LEADINGJETS : public Analysis {
+  class MC_LHC_SUSY : public Analysis {
   public:
     
     /// Constructor
-    MC_LHC_LEADINGJETS()
-      : Analysis("MC_LHC_LEADINGJETS")
+    MC_LHC_SUSY()
+      : Analysis("MC_LHC_SUSY")
     { 
       setBeams(PROTON, PROTON);
-      
-      // Final state for the jet finding
-      const FinalState fsj(-4.0, 4.0, 0.0*GeV);
-      addProjection(fsj, "FSJ");
-      addProjection(FastJets(fsj, FastJets::KT, 0.7), "Jets");
-      
-      // Charged final state for the distributions
-      const ChargedFinalState cfs(-1.0, 1.0, 0.5*GeV);
-      addProjection(cfs, "CFS");
     }
     
     
@@ -37,139 +31,209 @@
     
     // Book histograms
     void init() {
-      const double maxpt1 = 500.0/GeV;
-      _hist_pnchg      = bookProfile1D("trans-nchg", 50, 0.0, maxpt1);
-      _hist_pmaxnchg   = bookProfile1D("trans-maxnchg", 50, 0.0, maxpt1);
-      _hist_pminnchg   = bookProfile1D("trans-minnchg", 50, 0.0, maxpt1);
-      _hist_pcptsum    = bookProfile1D("trans-ptsum", 50, 0.0, maxpt1);
-      _hist_pmaxcptsum = bookProfile1D("trans-maxptsum", 50, 0.0, maxpt1);
-      _hist_pmincptsum = bookProfile1D("trans-minptsum", 50, 0.0, maxpt1);
-      _hist_pcptave    = bookProfile1D("trans-ptavg", 50, 0.0, maxpt1);
+      // Basic final state
+      const FinalState fs(-4.0, 4.0, 0.5*GeV);
+
+      // Tracks and jets
+      addProjection(ChargedFinalState(fs), "Tracks");
+      addProjection(FastJets(fs, FastJets::ANTIKT, 0.7), "Jets");
+
+      IdentifiedFinalState photonfs(fs);
+      photonfs.acceptId(PHOTON);
+      addProjection(photonfs, "AllPhotons");
+      /// @todo Isolated photons
+
+      IdentifiedFinalState efs(fs);
+      efs.acceptIdPair(ELECTRON);
+      addProjection(efs, "Electrons");
+
+      IdentifiedFinalState mufs(fs);
+      efs.acceptIdPair(MUON);
+      addProjection(mufs, "Muons");
+
+      addProjection(TotalVisibleMomentum(fs), "MET");
+      LeadingParticlesFinalState lpfs(fs);
+      lpfs.addParticleIdPair(ELECTRON);
+      lpfs.addParticleIdPair(MUON);
+      addProjection(lpfs, "LeadingParticles");
+
+      _hist_n_trk   = bookHistogram1D("n-trk", 50, 0.5, 300.5);
+      _hist_phi_trk = bookHistogram1D("phi-trk", 50, -PI, PI);
+      _hist_eta_trk = bookHistogram1D("eta-trk", 50, -4, 4);
+      _hist_pt_trk  = bookHistogram1D("pt-trk", 100, 0.0, 1500);
+
+      _hist_n_jet   = bookHistogram1D("n-jet", 21, -0.5, 20.5);
+      _hist_phi_jet = bookHistogram1D("phi-jet", 50, -PI, PI);
+      _hist_eta_jet = bookHistogram1D("eta-jet", 50, -4, 4);
+      _hist_pt_jet  = bookHistogram1D("pt-jet", 100, 0.0, 1500);
+
+      _hist_n_e   = bookHistogram1D("n-e", 11, -0.5, 10.5);
+      _hist_phi_e = bookHistogram1D("phi-e", 50, -PI, PI);
+      _hist_eta_e = bookHistogram1D("eta-e", 50, -4, 4);
+      _hist_pt_e  = bookHistogram1D("pt-e", 100, 0.0, 500);
+
+      _hist_n_mu   = bookHistogram1D("n-mu", 11, -0.5, 10.5);
+      _hist_phi_mu = bookHistogram1D("phi-mu", 50, -PI, PI);
+      _hist_eta_mu = bookHistogram1D("eta-mu", 50, -4, 4);
+      _hist_pt_mu  = bookHistogram1D("pt-mu", 100, 0.0, 500);
+
+      _hist_n_gamma   = bookHistogram1D("n-gamma", 11, -0.5, 10.5);
+      _hist_phi_gamma = bookHistogram1D("phi-gamma", 50, -PI, PI);
+      _hist_eta_gamma = bookHistogram1D("eta-gamma", 50, -4, 4);
+      _hist_pt_gamma  = bookHistogram1D("pt-gamma", 100, 0.0, 500);
+
+      /// @todo Isolated photons
+
+      _hist_met = bookHistogram1D("Etmiss", 100, 0.0, 1500);
+
+      _hist_mll_ossf_ee   = bookHistogram1D("mll-ossf-ee", 50, 0.0, 500);
+      _hist_mll_ossf_mumu = bookHistogram1D("mll-ossf-mumu", 50, 0.0, 500);
+      _hist_mll_osof_emu  = bookHistogram1D("mll-osof-mumu", 50, 0.0, 500);
+
+      // LSP eta, pT, phi, mass
     }
 
 
     // Do the analysis
-    void analyze(const Event& e) {
-
-      const FinalState& fsj = applyProjection<FinalState>(e, "FSJ");
-      if (fsj.particles().empty()) {
+    void analyze(const Event& evt) {
+      const FinalState& tracks = applyProjection<FinalState>(evt, "Tracks");
+      if (tracks.particles().empty()) {
         getLog() << Log::DEBUG << "Failed multiplicity cut" << endl;
         vetoEvent;
       }
 
-      const FastJets& jetpro = applyProjection<FastJets>(e, "Jets");
+      // Get event weight
+      const double weight = evt.weight();
+
+      // Fill track histos
+      _hist_n_trk->fill(tracks.size(), weight);
+      foreach (const Particle& t, tracks.particles()) {
+        const FourMomentum& p = t.momentum();
+        _hist_phi_trk->fill(p.phi(), weight);
+        _hist_eta_trk->fill(p.eta(), weight);
+        _hist_pt_trk->fill(p.pT()/GeV, weight);
+      }
+
+      // Get jets and fill jet histos
+      const FastJets& jetpro = applyProjection<FastJets>(evt, "Jets");
       const Jets jets = jetpro.jetsByPt();
       getLog() << Log::DEBUG << "Jet multiplicity = " << jets.size() << endl;
+      _hist_n_jet->fill(jets.size(), weight);
+      foreach (const Jet& j, jets) {
+        const FourMomentum& pj = j.momentum();
+        _hist_phi_jet->fill(pj.phi(), weight);
+        _hist_eta_jet->fill(pj.eta(), weight);
+        _hist_pt_jet->fill(pj.pT()/GeV, weight);
+      }
 
-      // We require the leading jet to be within |eta| < 2
-      if (jets.size() < 1 || fabs(jets[0].momentum().pseudorapidity()) > 2) {
-        getLog() << Log::DEBUG << "Failed jet cut" << endl;
-        vetoEvent;
+      /// @todo Resum photons around electrons
+
+      // Fill final state electron/positron histos
+      const FinalState& efs = applyProjection<FinalState>(evt, "Electrons");
+      _hist_n_e->fill(efs.size(), weight);
+      foreach (const Particle& e, efs.particles()) {
+        const FourMomentum& p = e.momentum();
+        _hist_phi_e->fill(p.phi(), weight);
+        _hist_eta_e->fill(p.eta(), weight);
+        _hist_pt_e->fill(p.pT()/GeV, weight);
       }
 
-      const double jetphi = jets[0].momentum().phi();
-      const double jetpT  = jets[0].momentum().pT();
-      getLog() << Log::DEBUG << "Leading jet: pT = " << jetpT
-               << ", eta = " << jets[0].momentum().pseudorapidity()
-               << ", phi = " << jetphi << endl;
-
-      // Get the event weight
-      const double weight = e.weight();
-
-      // Get the final states to work with for filling the distributions
-      const FinalState& cfs = applyProjection<ChargedFinalState>(e, "CFS");
-
-      size_t   numOverall(0),     numToward(0),     numTrans1(0),     numTrans2(0),     numAway(0);
-      double ptSumOverall(0.0), ptSumToward(0.0), ptSumTrans1(0.0), ptSumTrans2(0.0), ptSumAway(0.0);
-      double ptMaxOverall(0.0), ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);
-
-      // Calculate all the charged stuff
-      foreach (const Particle& p, cfs.particles()) {
-        const double dPhi = deltaPhi(p.momentum().phi(), jetphi);
-        const double pT = p.momentum().pT();
-        const double phi = p.momentum().azimuthalAngle();
-        const double rotatedphi = phi - jetphi;
-
-        ptSumOverall += pT;
-        ++numOverall;
-        if (pT > ptMaxOverall) ptMaxOverall = pT;
-
-        if (dPhi < PI/3.0) {
-          ptSumToward += pT;
-          ++numToward;
-          if (pT > ptMaxToward) ptMaxToward = pT;
-        }
-        else if (dPhi < 2*PI/3.0) {
-          if (rotatedphi <= PI) {
-            ptSumTrans1 += pT;
-            ++numTrans1;
-            if (pT > ptMaxTrans1) {
-              ptMaxTrans1 = pT;
-            } else {
-              ptSumTrans2 += pT;
-              ++numTrans2;
-              if (pT > ptMaxTrans2) ptMaxTrans2 = pT;
-            }
-          }
-        }
-        else {
-          ptSumAway += pT;
-          ++numAway;
-          if (pT > ptMaxAway) ptMaxAway = pT;
+      /// @todo Resum photons around muons
+
+      // Fill final state muon/antimuon histos
+      const FinalState& mufs = applyProjection<FinalState>(evt, "Muons");
+      _hist_n_e->fill(mufs.size(), weight);
+      foreach (const Particle& mu, efs.particles()) {
+        const FourMomentum& p = mu.momentum();
+        _hist_phi_mu->fill(p.phi(), weight);
+        _hist_eta_mu->fill(p.eta(), weight);
+        _hist_pt_mu->fill(p.pT()/GeV, weight);
+      }
+
+      // Fill final state non-isolated photon histos
+      const FinalState& allphotonfs = applyProjection<FinalState>(evt, "AllPhotons");
+      _hist_n_e->fill(allphotonfs.size(), weight);
+      foreach (const Particle& ph, allphotonfs.particles()) {
+        const FourMomentum& p = ph.momentum();
+        _hist_phi_mu->fill(p.phi(), weight);
+        _hist_eta_mu->fill(p.eta(), weight);
+        _hist_pt_mu->fill(p.pT()/GeV, weight);
+      }
+
+      /// @todo Isolated photons
+
+      // Calculate and fill missing Et histos
+      const TotalVisibleMomentum& met = applyProjection<TotalVisibleMomentum>(evt, "MET");
+      _hist_met->fill(met.scalarET()/GeV);
+
+      // Choose highest-pT leptons of each sign and flavour for dilepton mass edges
+      const FinalState& lpfs = applyProjection<FinalState>(evt, "LeadingParticles");
+      bool eplus_ok, eminus_ok, muplus_ok, muminus_ok;
+      FourMomentum peplus, peminus, pmuplus, pmuminus;
+      foreach (const Particle& p, lpfs.particles()) {
+        const PdgId pid = p.pdgId();
+        if (pid == ELECTRON) {
+          eminus_ok = true;
+          peminus = p.momentum();
+        } else if (pid == POSITRON) {
+          eplus_ok = true;
+          peplus = p.momentum();
+        } else if (pid == MUON) {
+          muminus_ok = true;
+          pmuminus = p.momentum();
+        } else if (pid == ANTIMUON) {
+          muplus_ok = true;
+          pmuplus = p.momentum();
+        } else {
+          throw Error("Unexpected particle type in leading particles FS!");
         }
       }
-      
-      
-      // Fill the histograms
-      //_hist_tnchg->fill(jetpT, numToward/(4*PI/3), weight);
-      _hist_pnchg->fill(jetpT, (numTrans1+numTrans2)/(4*PI/3), weight);
-      _hist_pmaxnchg->fill(jetpT, (numTrans1>numTrans2 ? numTrans1 : numTrans2)/(2*PI/3), weight);
-      _hist_pminnchg->fill(jetpT, (numTrans1<numTrans2 ? numTrans1 : numTrans2)/(2*PI/3), weight);
-      //_hist_pdifnchg->fill(jetpT, abs(numTrans1-numTrans2)/(2*PI/3), weight);
-      //_hist_anchg->fill(jetpT, numAway/(4*PI/3), weight);
-      
-      //_hist_tcptsum->fill(jetpT, ptSumToward/(4*PI/3), weight);
-      _hist_pcptsum->fill(jetpT, (ptSumTrans1+ptSumTrans2)/(4*PI/3), weight);
-      _hist_pmaxcptsum->fill(jetpT, (ptSumTrans1>ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/(2*PI/3), weight);
-      _hist_pmincptsum->fill(jetpT, (ptSumTrans1<ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/(2*PI/3), weight);
-      //_hist_pdifcptsum->fill(jetpT, fabs(ptSumTrans1-ptSumTrans2)/(2*PI/3), weight);
-      //_hist_acptsum->fill(jetpT, ptSumAway/(4*PI/3), weight);
-      
-      //if (numToward > 0) {
-      //  _hist_tcptave->fill(jetpT, ptSumToward/numToward, weight);
-      //  _hist_tcptmax->fill(jetpT, ptMaxToward, weight);
-      //}
-      if ((numTrans1+numTrans2) > 0) {
-        _hist_pcptave->fill(jetpT, (ptSumTrans1+ptSumTrans2)/(numTrans1+numTrans2), weight);
-        //_hist_pcptmax->fill(jetpT, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2), weight);
-      }
-      //if (numAway > 0) {
-      //  _hist_acptave->fill(jetpT, ptSumAway/numAway, weight);
-      //  _hist_acptmax->fill(jetpT, ptMaxAway, weight);
-      //}
+      // m_ee
+      if (eminus_ok && eplus_ok) {
+        const double m_ee = FourMomentum(peplus + peminus).mass();
+        _hist_mll_ossf_ee->fill(m_ee/GeV, weight);
+      }
+      // m_mumu
+      if (muminus_ok && muplus_ok) {
+        const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
+        _hist_mll_ossf_mumu->fill(m_mumu/GeV, weight);
+      }
+      // m_emu (both configurations)
+      if (eminus_ok && muplus_ok) {
+        const double m_emu = FourMomentum(pmuplus + peminus).mass();
+        _hist_mll_osof_emu->fill(m_emu/GeV, weight);
+      }
+      if (muminus_ok && eplus_ok) {
+        const double m_mue = FourMomentum(peplus + pmuminus).mass();
+        _hist_mll_osof_emu->fill(m_mue/GeV, weight);
+      }
+
     }
     
     
     void finalize() {  
-      //
+      /// @todo Normalisations
     }
+
+    //@}    
     
-    
+
   private:
     
-    AIDA::IProfile1D *_hist_pnchg;
-    AIDA::IProfile1D *_hist_pmaxnchg;
-    AIDA::IProfile1D *_hist_pminnchg;
-    AIDA::IProfile1D *_hist_pcptsum;
-    AIDA::IProfile1D *_hist_pmaxcptsum;
-    AIDA::IProfile1D *_hist_pmincptsum;
-    AIDA::IProfile1D *_hist_pcptave;  
+    AIDA::IHistogram1D *_hist_n_trk, *_hist_phi_trk, *_hist_eta_trk, *_hist_pt_trk;
+    AIDA::IHistogram1D *_hist_n_jet, *_hist_phi_jet, *_hist_eta_jet, *_hist_pt_jet;
+    AIDA::IHistogram1D *_hist_n_e, *_hist_phi_e, *_hist_eta_e, *_hist_pt_e;
+    AIDA::IHistogram1D *_hist_n_mu, *_hist_phi_mu, *_hist_eta_mu, *_hist_pt_mu;
+    AIDA::IHistogram1D *_hist_n_gamma, *_hist_phi_gamma, *_hist_eta_gamma, *_hist_pt_gamma;
+    /// @todo Isolated photons
+    AIDA::IHistogram1D *_hist_met;
+    AIDA::IHistogram1D *_hist_mll_ossf_ee, *_hist_mll_ossf_mumu, *_hist_mll_osof_emu;
     
   };
   
   
   
   // This global object acts as a hook for the plugin system
-  AnalysisBuilder<MC_LHC_LEADINGJETS> plugin_MC_LHC_LEADINGJETS;
+  AnalysisBuilder<MC_LHC_SUSY> plugin_MC_LHC_SUSY;
   
 }

Modified: trunk/src/Analyses/MC_LHC_WANALYSIS.cc
==============================================================================
--- trunk/src/Analyses/MC_LHC_WANALYSIS.cc	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Analyses/MC_LHC_WANALYSIS.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -16,13 +16,6 @@
     MC_LHC_WANALYSIS()
       : Analysis("MC_LHC_WANALYSIS") 
     {
-      const ChargedFinalState cfs;
-      addProjection(cfs, "CFS");
-      /// @todo Handle muon-decay Ws as well
-      const WFinder wf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 30.0*GeV, 110.0*GeV, 0.2);
-      addProjection(wf, "WF");
-      FastJets fastjets(wf.remainingFinalState(), FastJets::KT, 0.7);
-      addProjection(fastjets, "Jets");
     }
     
 
@@ -31,6 +24,14 @@
     //@{
 
     void init() { 
+      const ChargedFinalState cfs;
+      addProjection(cfs, "CFS");
+      /// @todo Handle muon-decay Ws as well
+      const WFinder wf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 30.0*GeV, 110.0*GeV, 0.2);
+      addProjection(wf, "WF");
+      FastJets fastjets(wf.remainingFinalState(), FastJets::KT, 0.7);
+      addProjection(fastjets, "Jets");
+
       _hist_chargemulti = bookHistogram1D("d01-x01-y01", 30, 0.5, 250.5);
       _hist_chargept = bookHistogram1D("d02-x01-y01", 32, 0., 25.);
       _hist_chargemeanpt = bookHistogram1D("d03-x01-y01", 25, 0., 10.);

Modified: trunk/src/Analyses/MC_LHC_ZANALYSIS.cc
==============================================================================
--- trunk/src/Analyses/MC_LHC_ZANALYSIS.cc	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Analyses/MC_LHC_ZANALYSIS.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -16,13 +16,6 @@
     MC_LHC_ZANALYSIS()
       : Analysis("MC_LHC_ZANALYSIS") 
     {
-      const ChargedFinalState cfs;
-      addProjection(cfs, "CFS");
-      /// @todo Handle muon-decay Zs as well
-      const ZFinder zf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 30.0*GeV, 115.0*GeV, 0.2);
-      addProjection(zf, "ZF");
-      FastJets fastjets(zf.remainingFinalState(), FastJets::KT, 0.7);
-      addProjection(fastjets, "Jets");
     }
     
 
@@ -31,6 +24,14 @@
     //@{
 
     void init() { 
+      const ChargedFinalState cfs;
+      addProjection(cfs, "CFS");
+      /// @todo Handle muon-decay Zs as well
+      const ZFinder zf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 30.0*GeV, 115.0*GeV, 0.2);
+      addProjection(zf, "ZF");
+      FastJets fastjets(zf.remainingFinalState(), FastJets::KT, 0.7);
+      addProjection(fastjets, "Jets");
+
       _hist_chargemulti = bookHistogram1D("d01-x01-y01", 30, 0.5, 250.5);
       _hist_chargept = bookHistogram1D("d02-x01-y01", 32, 0., 25.);
       _hist_chargemeanpt = bookHistogram1D("d03-x01-y01", 25, 0., 10.);

Modified: trunk/src/Analyses/Makefile.am
==============================================================================
--- trunk/src/Analyses/Makefile.am	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Analyses/Makefile.am	Thu Oct  1 18:21:30 2009	(r1844)
@@ -109,4 +109,5 @@
     MC_LHC_LEADINGJETS.cc \
     MC_LHC_DIJET.cc \
     MC_LHC_WANALYSIS.cc \
-    MC_LHC_ZANALYSIS.cc
+    MC_LHC_ZANALYSIS.cc \
+    MC_LHC_SUSY.cc

Modified: trunk/src/Projections/ChargedLeptons.cc
==============================================================================
--- trunk/src/Projections/ChargedLeptons.cc	Thu Oct  1 17:04:43 2009	(r1843)
+++ trunk/src/Projections/ChargedLeptons.cc	Thu Oct  1 18:21:30 2009	(r1844)
@@ -8,25 +8,21 @@
 
 
   int ChargedLeptons::compare(const Projection& other) const {
-    return mkNamedPCmp(other, "FS");
+    return mkNamedPCmp(other, "ChFS");
   }
 
 
-  void ChargedLeptons::project(const Event& e) {
+  void ChargedLeptons::project(const Event& evt) {
     // Reset result
     _theChargedLeptons.clear();
 
-    // Get hadron and charge info for each particle, and fill counters appropriately
-    const FinalState& fs = applyProjection<FinalState>(e, "FS");
+    // Loop over charged particles and fill vector with leptons
+    const FinalState& fs = applyProjection<FinalState>(evt, "ChFS");
     foreach (const Particle& p, fs.particles()) {
       if (PID::isLepton(p.pdgId())) {
-        if (PID::threeCharge(p.pdgId()) != 0) {
-          // Put it into the C.L. vector
-          _theChargedLeptons.push_back(Particle(p));
-        }
+        _theChargedLeptons += Particle(p);
       }
     }
-    getLog() << Log::DEBUG << "Done" << endl;
   }
 
 


More information about the Rivet-svn mailing list