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);
+
+
+}