Index: trunk/code/PL_JETMASS.plot
===================================================================
--- trunk/code/PL_JETMASS.plot	(revision 0)
+++ trunk/code/PL_JETMASS.plot	(revision 473)
@@ -0,0 +1,7 @@
+# BEGIN PLOT /PL_JETMASS/jetmass
+LogY=0
+# END PLOT
+
+# BEGIN PLOT /PL_JETMASS/jetmass_norec
+LogY=0
+# END PLOT
Index: trunk/code/PL_JETMASS.info
===================================================================
--- trunk/code/PL_JETMASS.info	(revision 0)
+++ trunk/code/PL_JETMASS.info	(revision 473)
@@ -0,0 +1,34 @@
+Name: PL_JETMASS
+Year: <Insert year of publication>
+Summary: <Insert short PL_JETMASS description>
+Experiment: <Insert the experiment name>
+Collider: <Insert the collider name>
+SpiresID: <Insert the Inspire ID>
+Status: UNVALIDATED
+Authors:
+ - Your Name <your@email.address>
+#References:
+# - <Example: Phys.Lett.B639:151-158,2006, Erratum-ibid.B658:285-289,2008>
+# - <Example: doi:10.1016/j.physletb.2006.04.048>
+# - <Example: arXiv:hep-ex/0511054 (plus erratum)>
+RunInfo:
+  <Insert event types (not gen-specific), energy, any kinematic
+  efficiency cut(s) that may be needed; essentially any details needed to set
+  up a generator to reproduce the data.>
+NumEvents: 1000000
+NeedCrossSection: no
+#Beams: <Insert beam pair(s), e.g. [p-, p+] or [[p-, e-], [p-, e+]]>
+#Energies: <Insert list of run energies or beam energy pairs in GeV,
+#           e.g. [1960] or [[8.0, 3.5]] or [630, 1800]. Order pairs to match "Beams">
+#PtCuts: <Insert list of kinematic pT cuts in GeV, e.g. [0, 20]>
+#NeedCrossSection: True
+Description:
+  '<Insert a fairly long description, including what is measured
+  and if possible what it is useful for in terms of MC validation
+  and tuning. Use LaTeX for maths like $\pT > 50\;\GeV$.
+  Use single quotes around the block if required (see http://www.yaml.org)>'
+BibKey: 
+BibTeX: ''
+ToDo:
+ - Implement the analysis, test it, remove this ToDo, and mark as VALIDATED :-)
+
Index: trunk/code/PL_JETMASS.cc
===================================================================
--- trunk/code/PL_JETMASS.cc	(revision 0)
+++ trunk/code/PL_JETMASS.cc	(revision 473)
@@ -0,0 +1,178 @@
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/ClusterSequence.hh"
+
+namespace Rivet {
+
+
+  using namespace fastjet;
+
+  class PL_JETMASS : public Analysis {
+  public:
+
+    /// Constructor
+    PL_JETMASS()
+      : Analysis("PL_JETMASS")
+    {    }
+
+
+    /// Book histograms and initialise projections before the run
+    void init() {
+
+	    _h_jetmass = bookHisto1D("jetmass", 25, 0., 50.);
+	    _h_jetmass_norec = bookHisto1D("jetmass_norec", 25, 0., 50.);
+	    _h_jetmass_unsub = bookHisto1D("jetmass_unsub", 25, 0., 50.);
+
+    }
+
+
+    /// 4-momentum subtraction
+    FourMomentum SubtractJetMom(PseudoJet pjet, const vector<PseudoJet> * scatcenters) {
+        if (scatcenters->empty()) return Get4Mom(pjet);
+        FourMomentum sub(0.,0.,0.,0.);
+        foreach (PseudoJet pj, pjet.constituents()) {
+            if (pj.E() < 1e-5 && pj.E() > 1e-7) {
+	     FourMomentum test(Get4Mom(pj));
+	     test *= 10000.;
+	     //cout<<"Found dummy :"<<Get4Mom(pj)<<endl;
+	     double dRmin(10.);
+	     PseudoJet matchscatcen;
+	     foreach (PseudoJet sc, *scatcenters) {
+			double dR = deltaR(test.eta(), test.phi(ZERO_2PI), Get4Mom(sc).eta(),Get4Mom(sc).phi(ZERO_2PI));
+			if (dR < dRmin) {
+			    dRmin = dR;
+			    matchscatcen = sc;
+			}
+	     }
+	     if (dRmin < 1e-5) {
+		    //cout<<"Found matching scattering centre :"<<Get4Mom(matchscatcen)<<endl;
+		    sub += Get4Mom(matchscatcen);
+	     }
+	     else cout<<"Error: did not find scattering centre matching dummy.\n";
+	     } 
+        }
+        //cout<<"size of map: "<<scmap.size()<<endl;
+        FourMomentum jetmom(Get4Mom(pjet));
+        //cout<<"Original momentum: "<<jetmom<<endl;
+        jetmom -= sub;
+        //cout<<"subtracted momentum: "<<jetmom<<endl;
+        //cout<<"Subtracted pt: "<<subpt<<"    "<<pjet.pt()<<endl;
+        return jetmom;
+    }
+    
+
+     
+        
+    /// 4-momentum adapter
+    FourMomentum Get4Mom(PseudoJet pjet){
+        FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz());
+        return mom;
+    }
+
+    
+    
+    
+    /// Perform the per-event analysis
+    void analyze(const Event& event) {
+      const double weight = event.weight();
+
+      PseudoJets pevent, pjets, pjetsnorec;
+      PseudoJets evnorec;      
+      list<double> ptsubs;
+        PseudoJets scatcens; 
+        vector<FourMomentum> scafter;
+        foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
+            FourMomentum mom(p->momentum());
+            if (fabs(mom.eta()) < 4.) {
+               if(p->status() == 1) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+                   pevent.push_back(pjet);
+                   evnorec.push_back(pjet);
+               }
+               if (p->status() == 3) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+                   scatcens.push_back(pjet);
+               }
+               if(p->status() == 4) {
+                   PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E());
+                   pevent.push_back(pjet);
+                   scafter.push_back(mom);
+               }
+            }
+        }
+      
+      JetDefinition jet_def(antikt_algorithm, 1.0);
+      ClusterSequence cs(pevent, jet_def);
+      pjets = sorted_by_pt(cs.inclusive_jets(50.0));
+
+      foreach (PseudoJet pjet, pjets) {
+         if (fabs(pjet.eta()) < 2.) {
+		  
+	      double mass, mass2;	  
+	      FourMomentum punsub(Get4Mom(pjet));
+	      if (punsub.pT() > 60.) {
+		mass2 = punsub.mass2();
+		mass = mass2>0.?sqrt(mass2):0.;
+		_h_jetmass_unsub->fill(mass, weight);
+	      }
+	      
+             FourMomentum jetsub  = SubtractJetMom(pjet, &scatcens);
+             double ptsub;
+             if (jetsub.E() < 0. || deltaR(jetsub.eta(), jetsub.phi(ZERO_2PI), Get4Mom(pjet).eta(), Get4Mom(pjet).phi(ZERO_2PI)) > 1.) ptsub = 0.;
+             else ptsub = jetsub.pT();
+
+             if (ptsub < 60.) continue;
+	      
+             mass2 = jetsub.mass2();
+             mass = mass2>0.?sqrt(mass2):0.;
+	      _h_jetmass->fill(mass, weight);
+	  }
+      }
+
+      ClusterSequence csnorec(evnorec, jet_def);
+      pjetsnorec = sorted_by_pt(csnorec.inclusive_jets(60.0));
+      
+      foreach (PseudoJet pjet, pjetsnorec) {
+         if (fabs(pjet.eta()) < 2.) {
+      
+	      FourMomentum jetmom(Get4Mom(pjet));
+		  
+             double mass2(jetmom.mass2());
+             double mass(mass2>0.?sqrt(mass2):0.);
+	      _h_jetmass_norec->fill(mass, weight);
+	  }
+      }
+      
+	      
+    }
+
+
+    /// Normalise histograms etc., after the run
+    void finalize() {
+
+       normalize(_h_jetmass);
+       normalize(_h_jetmass_norec);
+       normalize(_h_jetmass_unsub);
+
+    }
+
+
+  private:
+
+    // Data members like post-cuts event weight counters go here
+
+
+    Histo1DPtr _h_jetmass, _h_jetmass_norec, _h_jetmass_unsub;
+
+
+  };
+
+
+
+  // The hook for the plugin system
+  DECLARE_RIVET_PLUGIN(PL_JETMASS);
+
+
+}