Page MenuHomeHEPForge
Files F7877438

No OneTemporary
Actions

View Options

diff --git a/data/plotinfo/H1_1995_S3167097.plot b/data/plotinfo/H1_1995_S3167097.plot
--- a/data/plotinfo/H1_1995_S3167097.plot
+++ b/data/plotinfo/H1_1995_S3167097.plot
@@ -1,78 +1,59 @@
# BEGIN PLOT /H1_1995_S3167097/1
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 11
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/2
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 12
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
-# BEGIN PLOT /H1_1995_S3167097/21tmp
-Title=$\langle E_\perp \rangle$ vs kin. bin
-XLabel=Kinematic bin
-YLabel=$\langle \mathrm{E}_{\perp}\rangle$
-# END PLOT
-
-# BEGIN PLOT /H1_1995_S3167097/22tmp
-Title=$\langle x \rangle$ vs kin. bin
-XLabel=Kinematic bin
-YLabel=$\langle x \rangle$
-# END PLOT
-
-# BEGIN PLOT /H1_1995_S3167097/23tmp
-Title=$\langle Q^2 \rangle$ vs kin. bin
-XLabel=Kinematic bin
-YLabel=$\langle Q^{2}\rangle$
-# END PLOT
-
# BEGIN PLOT /H1_1995_S3167097/24
Title=Num events vs kin. bin
XLabel=Kinematic bin
YLabel=N
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/3
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 13
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/4
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 14
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/5
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 15
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/6
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 16
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/7
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 17
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/8
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 18
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
# BEGIN PLOT /H1_1995_S3167097/9
Title=Stat $\mathrm{d}E_\perp/\mathrm{d}[c]$ CMS bin = 19
XLabel=$\eta$
YLabel=$1/\mathrm{N} \mathrm{d}E_{\perp}/d\eta [GeV]$
# END PLOT
-
diff --git a/src/Analyses/ATLAS_2012_I1083318.cc b/src/Analyses/ATLAS_2012_I1083318.cc
--- a/src/Analyses/ATLAS_2012_I1083318.cc
+++ b/src/Analyses/ATLAS_2012_I1083318.cc
@@ -1,264 +1,262 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetYODA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/LeptonClusters.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
namespace Rivet {
class ATLAS_2012_I1083318 : public Analysis {
public:
/// @name Constructors etc.
//@{
/// Constructor
ATLAS_2012_I1083318()
: Analysis("ATLAS_2012_I1083318")
{ }
//@}
public:
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
FinalState fs;
IdentifiedFinalState allleptons;
allleptons.acceptIdPair(PID::ELECTRON);
allleptons.acceptIdPair(PID::MUON);
vector< pair<double, double> > etaRanges;
etaRanges.push_back(make_pair(-2.5, 2.5));
- LeptonClusters leptons(fs, allleptons,
- 0.1, true,
- etaRanges, 20.0*GeV);
+ LeptonClusters leptons(fs, allleptons, 0.1, true, etaRanges, 20.0*GeV);
addProjection(leptons, "leptons");
// Leading neutrinos for Etmiss
LeadingParticlesFinalState neutrinos(fs);
neutrinos.addParticleIdPair(PID::NU_E);
neutrinos.addParticleIdPair(PID::NU_MU);
neutrinos.setLeadingOnly(true);
addProjection(neutrinos, "neutrinos");
// Input for the jets: "Neutrinos, electrons, and muons from decays of the
// massive W boson were not used"
VetoedFinalState veto;
veto.addVetoOnThisFinalState(leptons);
veto.addVetoOnThisFinalState(neutrinos);
FastJets jets(veto, FastJets::ANTIKT, 0.4);
jets.useInvisibles(true);
addProjection(jets, "jets");
for (size_t i = 0; i < 2; ++i) {
_h_NjetIncl[i] = bookHisto1D(1, 1, i+1);
_h_RatioNjetIncl[i] = bookScatter2D(2, 1, i+1);
_h_FirstJetPt_1jet[i] = bookHisto1D(3, 1, i+1);
_h_FirstJetPt_2jet[i] = bookHisto1D(4, 1, i+1);
_h_FirstJetPt_3jet[i] = bookHisto1D(5, 1, i+1);
_h_FirstJetPt_4jet[i] = bookHisto1D(6, 1, i+1);
_h_SecondJetPt_2jet[i] = bookHisto1D(7, 1, i+1);
_h_SecondJetPt_3jet[i] = bookHisto1D(8, 1, i+1);
_h_SecondJetPt_4jet[i] = bookHisto1D(9, 1, i+1);
_h_ThirdJetPt_3jet[i] = bookHisto1D(10, 1, i+1);
_h_ThirdJetPt_4jet[i] = bookHisto1D(11, 1, i+1);
_h_FourthJetPt_4jet[i] = bookHisto1D(12, 1, i+1);
_h_Ht_1jet[i] = bookHisto1D(13, 1, i+1);
_h_Ht_2jet[i] = bookHisto1D(14, 1, i+1);
_h_Ht_3jet[i] = bookHisto1D(15, 1, i+1);
_h_Ht_4jet[i] = bookHisto1D(16, 1, i+1);
_h_Minv_2jet[i] = bookHisto1D(17, 1, i+1);
_h_Minv_3jet[i] = bookHisto1D(18, 1, i+1);
_h_Minv_4jet[i] = bookHisto1D(19, 1, i+1);
_h_JetRapidity[i] = bookHisto1D(20, 1, i+1);
_h_DeltaYElecJet[i] = bookHisto1D(21, 1, i+1);
_h_SumYElecJet[i] = bookHisto1D(22, 1, i+1);
_h_DeltaR_2jet[i] = bookHisto1D(23, 1, i+1);
_h_DeltaY_2jet[i] = bookHisto1D(24, 1, i+1);
_h_DeltaPhi_2jet[i] = bookHisto1D(25, 1, i+1);
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
const vector<ClusteredLepton>& leptons = applyProjection<LeptonClusters>(event, "leptons").clusteredLeptons();
Particles neutrinos = applyProjection<FinalState>(event, "neutrinos").particlesByPt();
if (leptons.size() != 1 || (neutrinos.size() == 0)) {
vetoEvent;
}
FourMomentum lepton = leptons[0].momentum();
FourMomentum p_miss = neutrinos[0].momentum();
if (p_miss.Et() < 25.0*GeV) {
vetoEvent;
}
double mT = sqrt(2.0 * lepton.pT() * p_miss.Et() * (1.0 - cos( lepton.phi()-p_miss.phi()) ) );
if (mT < 40.0*GeV) {
vetoEvent;
}
double jetcuts[] = { 30.0*GeV, 20.0*GeV };
const FastJets& jetpro = applyProjection<FastJets>(event, "jets");
for (size_t i = 0; i < 2; ++i) {
vector<FourMomentum> jets;
double HT = lepton.pT() + p_miss.pT();
foreach (const Jet& jet, jetpro.jetsByPt(jetcuts[i])) {
if (fabs(jet.momentum().rapidity()) < 4.4 && deltaR(lepton, jet.momentum()) > 0.5) {
jets.push_back(jet.momentum());
HT += jet.momentum().pT();
}
}
_h_NjetIncl[i]->fill(0.0, weight);
// Njet>=1 observables
if (jets.size() < 1) continue;
_h_NjetIncl[i]->fill(1.0, weight);
_h_FirstJetPt_1jet[i]->fill(jets[0].pT(), weight);
_h_JetRapidity[i]->fill(jets[0].rapidity(), weight);
_h_Ht_1jet[i]->fill(HT, weight);
_h_DeltaYElecJet[i]->fill(lepton.rapidity()-jets[0].rapidity(), weight);
_h_SumYElecJet[i]->fill(lepton.rapidity()+jets[0].rapidity(), weight);
// Njet>=2 observables
if (jets.size() < 2) continue;
_h_NjetIncl[i]->fill(2.0, weight);
_h_FirstJetPt_2jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_2jet[i]->fill(jets[1].pT(), weight);
_h_Ht_2jet[i]->fill(HT, weight);
double m2_2jet = FourMomentum(jets[0]+jets[1]).mass2();
_h_Minv_2jet[i]->fill(m2_2jet>0.0 ? sqrt(m2_2jet) : 0.0, weight);
_h_DeltaR_2jet[i]->fill(deltaR(jets[0], jets[1]), weight);
_h_DeltaY_2jet[i]->fill(jets[0].rapidity()-jets[1].rapidity(), weight);
_h_DeltaPhi_2jet[i]->fill(deltaPhi(jets[0], jets[1]), weight);
// Njet>=3 observables
if (jets.size() < 3) continue;
_h_NjetIncl[i]->fill(3.0, weight);
_h_FirstJetPt_3jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_3jet[i]->fill(jets[1].pT(), weight);
_h_ThirdJetPt_3jet[i]->fill(jets[2].pT(), weight);
_h_Ht_3jet[i]->fill(HT, weight);
double m2_3jet = FourMomentum(jets[0]+jets[1]+jets[2]).mass2();
_h_Minv_3jet[i]->fill(m2_3jet>0.0 ? sqrt(m2_3jet) : 0.0, weight);
// Njet>=4 observables
if (jets.size() < 4) continue;
_h_NjetIncl[i]->fill(4.0, weight);
_h_FirstJetPt_4jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_4jet[i]->fill(jets[1].pT(), weight);
_h_ThirdJetPt_4jet[i]->fill(jets[2].pT(), weight);
_h_FourthJetPt_4jet[i]->fill(jets[3].pT(), weight);
_h_Ht_4jet[i]->fill(HT, weight);
double m2_4jet = FourMomentum(jets[0]+jets[1]+jets[2]+jets[3]).mass2();
_h_Minv_4jet[i]->fill(m2_4jet>0.0 ? sqrt(m2_4jet) : 0.0, weight);
// Njet>=5 observables
if (jets.size() < 5) continue;
_h_NjetIncl[i]->fill(5.0, weight);
}
}
/// Normalise histograms etc., after the run
void finalize() {
for (size_t i = 0; i < 2; ++i) {
// Construct jet multiplicity ratio
- for (int n = 0; n < _h_NjetIncl[i]->numPoints()-1; ++n) {
- YODA::Bin& b0 = _h_NjetIncl[i]->bin(n);
- YODA::Bin& b1 = _h_NjetIncl[i]->bin(n+1);
+ for (size_t n = 0; n < _h_NjetIncl[i]->numBins()-1; ++n) {
+ YODA::HistoBin1D& b0 = _h_NjetIncl[i]->bin(n);
+ YODA::HistoBin1D& b1 = _h_NjetIncl[i]->bin(n+1);
if (b0.height() == 0.0 || b1.height() == 0.0) continue;
_h_RatioNjetIncl[i]->point(n).setY(b1.height()/b0.height());
_h_RatioNjetIncl[i]->point(n).setYErr(b1.height()/b0.height() * (b0.relErr() + b1.relErr()));
}
// Scale all histos to the cross section
const double factor = crossSection()/sumOfWeights();
scale(_h_DeltaPhi_2jet[i], factor);
scale(_h_DeltaR_2jet[i], factor);
scale(_h_DeltaY_2jet[i], factor);
scale(_h_DeltaYElecJet[i], factor);
scale(_h_FirstJetPt_1jet[i], factor);
scale(_h_FirstJetPt_2jet[i], factor);
scale(_h_FirstJetPt_3jet[i], factor);
scale(_h_FirstJetPt_4jet[i], factor);
scale(_h_FourthJetPt_4jet[i], factor);
scale(_h_Ht_1jet[i], factor);
scale(_h_Ht_2jet[i], factor);
scale(_h_Ht_3jet[i], factor);
scale(_h_Ht_4jet[i], factor);
scale(_h_JetRapidity[i], factor);
scale(_h_Minv_2jet[i], factor);
scale(_h_Minv_3jet[i], factor);
scale(_h_Minv_4jet[i], factor);
scale(_h_NjetIncl[i], factor);
scale(_h_SecondJetPt_2jet[i], factor);
scale(_h_SecondJetPt_3jet[i], factor);
scale(_h_SecondJetPt_4jet[i], factor);
scale(_h_SumYElecJet[i], factor);
scale(_h_ThirdJetPt_3jet[i], factor);
scale(_h_ThirdJetPt_4jet[i], factor);
}
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_DeltaPhi_2jet[2];
Histo1DPtr _h_DeltaR_2jet[2];
Histo1DPtr _h_DeltaY_2jet[2];
Histo1DPtr _h_DeltaYElecJet[2];
Histo1DPtr _h_FirstJetPt_1jet[2];
Histo1DPtr _h_FirstJetPt_2jet[2];
Histo1DPtr _h_FirstJetPt_3jet[2];
Histo1DPtr _h_FirstJetPt_4jet[2];
Histo1DPtr _h_FourthJetPt_4jet[2];
Histo1DPtr _h_Ht_1jet[2];
Histo1DPtr _h_Ht_2jet[2];
Histo1DPtr _h_Ht_3jet[2];
Histo1DPtr _h_Ht_4jet[2];
Histo1DPtr _h_JetRapidity[2];
Histo1DPtr _h_Minv_2jet[2];
Histo1DPtr _h_Minv_3jet[2];
Histo1DPtr _h_Minv_4jet[2];
Histo1DPtr _h_NjetIncl[2];
Scatter2DPtr _h_RatioNjetIncl[2];
Histo1DPtr _h_SecondJetPt_2jet[2];
Histo1DPtr _h_SecondJetPt_3jet[2];
Histo1DPtr _h_SecondJetPt_4jet[2];
Histo1DPtr _h_SumYElecJet[2];
Histo1DPtr _h_ThirdJetPt_3jet[2];
Histo1DPtr _h_ThirdJetPt_4jet[2];
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2012_I1083318);
}
diff --git a/src/Analyses/H1_1995_S3167097.cc b/src/Analyses/H1_1995_S3167097.cc
--- a/src/Analyses/H1_1995_S3167097.cc
+++ b/src/Analyses/H1_1995_S3167097.cc
@@ -1,136 +1,115 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetYODA.hh"
#include "Rivet/Projections/DISFinalState.hh"
#include "Rivet/Projections/CentralEtHCM.hh"
namespace Rivet {
/// H1 energy flow in DIS
///
- /// @todo Check!
+ /// @todo Make histograms match those in HepData and use autobooking
/// @author Leif Lonnblad
+ /// @author Andy Buckley
class H1_1995_S3167097 : public Analysis {
public:
/// Constructor
H1_1995_S3167097()
: Analysis("H1_1995_S3167097")
{ }
/// @name Analysis methods
//@{
void init() {
+ // Projections
const DISKinematics& diskin = addProjection(DISKinematics(), "Kinematics");
const DISFinalState& fshcm = addProjection(DISFinalState(diskin, DISFinalState::HCM), "FS");
addProjection(CentralEtHCM(fshcm), "Y1HCM");
- const size_t NBINS = 9;
- _hEtFlow = vector<Histo1DPtr>(NBINS);
- _hEtFlowStat = vector<Histo1DPtr>(NBINS);
- _nev = vector<double>(NBINS);
- /// @todo Automate this sort of thing so that the analysis code is more readable.
- for (size_t i = 0; i < NBINS; ++i) {
- string istr(1, char('1' + i));
- _hEtFlow[i] = bookHisto1D(istr, 24, -6, 6);
- _hEtFlowStat[i] = bookHisto1D(istr, 24, -6, 6);
- }
- /// @todo Replace with really temp (non-persistent) histos?
- // _hAvEt = bookScatter2D("21");
- // _hAvX = bookScatter2D("22");
- // _hAvQ2 = bookScatter2D("23");
- _hAvEt = bookHisto1D("21tmp", NBINS, 1.0, 10.0);
- _hAvX = bookHisto1D("22tmp", NBINS, 1.0, 10.0);
- _hAvQ2 = bookHisto1D("23tmp", NBINS, 1.0, 10.0);
- _hN = bookHisto1D("24", NBINS, 1.0, 10.0);
+ // Histograms
+ /// @todo Convert to use autobooking and correspond to HepData data tables
+ _sumw.resize(9);
+ _hEtFlow.resize(9);
+ for (size_t i = 0; i < 9; ++i)
+ _hEtFlow[i] = bookHisto1D(lexical_cast<string>(i), 24, -6, 6);
+ _tmphAvEt = Histo1D(9, 1.0, 10.0);
+ _tmphAvX = Histo1D(9, 1.0, 10.0);
+ _tmphAvQ2 = Histo1D(9, 1.0, 10.0);
+ _tmphN = Histo1D(9, 1.0, 10.0);
}
/// Calculate the bin number from the DISKinematics projection
size_t _getbin(const DISKinematics& dk) {
if (inRange(dk.Q2()/GeV2, 5.0, 10.0)) {
if (inRange(dk.x(), 1e-4, 2e-4)) return 0;
if (inRange(dk.x(), 2e-4, 5e-4) && dk.Q2() > 6.0*GeV2) return 1;
} else if (inRange(dk.Q2()/GeV2, 10.0, 20.0)) {
if (inRange(dk.x(), 2e-4, 5e-4)) return 2;
if (inRange(dk.x(), 5e-4, 8e-4)) return 3;
if (inRange(dk.x(), 8e-4, 1.5e-3)) return 4;
if (inRange(dk.x(), 1.5e-3, 4e-3)) return 5;
} else if (inRange(dk.Q2()/GeV2, 20.0, 50.0)) {
if (inRange(dk.x(), 5e-4, 1.4e-3)) return 6;
if (inRange(dk.x(), 1.4e-3, 3e-3)) return 7;
if (inRange(dk.x(), 3e-3, 1e-2)) return 8;
}
return -1;
}
void analyze(const Event& event) {
const FinalState& fs = applyProjection<FinalState>(event, "FS");
const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
const CentralEtHCM y1 = applyProjection<CentralEtHCM>(event, "Y1HCM");
const int ibin = _getbin(dk);
if (ibin < 0) vetoEvent;
const double weight = event.weight();
for (size_t i = 0, N = fs.particles().size(); i < N; ++i) {
const double rap = fs.particles()[i].momentum().rapidity();
const double et = fs.particles()[i].momentum().Et();
_hEtFlow[ibin]->fill(rap, weight * et/GeV);
- _hEtFlowStat[ibin]->fill(rap, weight * et/GeV);
}
- _nev[ibin] += weight;
- _hAvEt->fill(ibin + 1.5, weight * y1.sumEt()/GeV);
- _hAvX->fill(ibin + 1.5, weight * dk.x());
- _hAvQ2->fill(ibin + 1.5, weight * dk.Q2()/GeV2);
- _hN->fill(ibin + 1.5, weight);
+ _sumw[ibin] += weight;
+ _tmphAvEt.fill(ibin + 1.5, weight * y1.sumEt()/GeV);
+ _tmphAvX.fill(ibin + 1.5, weight * dk.x());
+ _tmphAvQ2.fill(ibin + 1.5, weight * dk.Q2()/GeV2);
+ _tmphN.fill(ibin + 1.5, weight);
}
void finalize() {
- for (size_t ibin = 0; ibin < 9; ++ibin) {
- scale(_hEtFlow[ibin], 0.5/_nev[ibin]);
- scale(_hEtFlowStat[ibin], 0.5/_nev[ibin]);
- }
-
- "/H1_1995_S3167097/21";
-
- divide(_tmphAvEt, *_hN, );
- h->setTitle(_hAvEt->title());
- histogramFactory().destroy(_hAvEt);
-
- h = histogramFactory().divide("/H1_1995_S3167097/22", *_hAvX, *_hN);
- h->setTitle(_hAvX->title());
- histogramFactory().destroy(_hAvX);
-
- h = histogramFactory().divide("/H1_1995_S3167097/23", *_hAvQ2, *_hN);
- h->setTitle(_hAvQ2->title());
- histogramFactory().destroy(_hAvQ2);
+ for (size_t ibin = 0; ibin < 9; ++ibin) scale(_hEtFlow[ibin], 0.5/_sumw[ibin]);
+ addPlot(Scatter2DPtr( new Scatter2D(_tmphAvEt/_tmphN, histoPath("21")) ));
+ addPlot(Scatter2DPtr( new Scatter2D(_tmphAvX/_tmphN, histoPath("22")) ));
+ addPlot(Scatter2DPtr( new Scatter2D(_tmphAvQ2/_tmphN, histoPath("23")) ));
}
//@}
private:
- /// Histograms for the \f$ E_T \f$ flows
- vector<Histo1DPtr> _hEtFlow, _hEtFlowStat;
+ /// Histograms for the \f$ E_T \f$ flow
+ vector<Histo1DPtr> _hEtFlow;
- /// Histograms for averages in different kinematical bins.
- Histo1DPtr _hAvEt, _hAvX, _hAvQ2, _hN;
+ /// Temporary histograms for averages in different kinematical bins.
+ Histo1D _tmphAvEt, _tmphAvX, _tmphAvQ2, _tmphN;
- /// Helper vector;
- vector<double> _nev;
+ /// Weights counters for each kinematic bin
+ vector<double> _sumw;
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(H1_1995_S3167097);
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:34 PM (1 d, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804970
Default Alt Text
(18 KB)

Event Timeline