|
|
[Rivet-svn] r4276 - in trunk/src: Analyses Core
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue May 14 18:54:14 BST 2013
Author: buckley
Date: Tue May 14 18:54:14 2013
New Revision: 4276
Log:
Another histogramming TODO done
Modified:
trunk/src/Analyses/ATLAS_2011_S9128077.cc
trunk/src/Core/Analysis.cc
Modified: trunk/src/Analyses/ATLAS_2011_S9128077.cc
==============================================================================
--- trunk/src/Analyses/ATLAS_2011_S9128077.cc Tue May 14 17:52:45 2013 (r4275)
+++ trunk/src/Analyses/ATLAS_2011_S9128077.cc Tue May 14 18:54:14 2013 (r4276)
@@ -11,18 +11,11 @@
class ATLAS_2011_S9128077 : public Analysis {
public:
- /// @name Constructors etc.
- //@{
-
/// Constructor
ATLAS_2011_S9128077()
: Analysis("ATLAS_2011_S9128077")
{ }
- //@}
-
-
- public:
/// @name Analysis methods
//@{
@@ -30,16 +23,16 @@
/// Book histograms and initialise projections before the run
void init() {
+ // Projections
const FinalState fs;
-
FastJets j4(fs, FastJets::ANTIKT, 0.4);
j4.useInvisibles();
addProjection(j4, "AntiKtJets04");
-
FastJets j6(fs, FastJets::ANTIKT, 0.6);
j6.useInvisibles();
addProjection(j6, "AntiKtJets06");
+ // Persistent histograms
_h_jet_multi_inclusive = bookHisto1D(1, 1, 1);
_h_jet_multi_ratio = bookScatter2D(2, 1, 1);
_h_jet_pT.resize(4);
@@ -50,28 +43,8 @@
_h_HT_2 = bookHisto1D(7, 1, 1);
_h_HT_3 = bookHisto1D(8, 1, 1);
_h_HT_4 = bookHisto1D(9, 1, 1);
-
- /// temporary histograms which will be divided in the end for the dsigma3/dsigma2 ratios
- _h_tmp_pTlead_R06_60_2 = bookHisto1D(10, 1, 1, "tmp1");
- _h_tmp_pTlead_R06_80_2 = bookHisto1D(11, 1, 1, "tmp2");
- _h_tmp_pTlead_R06_110_2 = bookHisto1D(12, 1, 1, "tmp3");
- _h_tmp_pTlead_R06_60_3 = bookHisto1D(10, 1, 1, "tmp4");
- _h_tmp_pTlead_R06_80_3 = bookHisto1D(11, 1, 1, "tmp5");
- _h_tmp_pTlead_R06_110_3 = bookHisto1D(12, 1, 1, "tmp6");
-
- _h_tmp_pTlead_R04_60_2 = bookHisto1D(13, 1, 1, "tmp7");
- _h_tmp_pTlead_R04_80_2 = bookHisto1D(14, 1, 1, "tmp8");
- _h_tmp_pTlead_R04_110_2 = bookHisto1D(15, 1, 1, "tmp9");
- _h_tmp_pTlead_R04_60_3 = bookHisto1D(13, 1, 1, "tmp10");
- _h_tmp_pTlead_R04_80_3 = bookHisto1D(14, 1, 1, "tmp11");
- _h_tmp_pTlead_R04_110_3 = bookHisto1D(15, 1, 1, "tmp12");
-
- _h_tmp_HT2_R06_2 = bookHisto1D(16, 1, 1, "tmp13");
- _h_tmp_HT2_R06_3 = bookHisto1D(16, 1, 1, "tmp14");
- _h_tmp_HT2_R04_2 = bookHisto1D(17, 1, 1, "tmp15");
- _h_tmp_HT2_R04_3 = bookHisto1D(17, 1, 1, "tmp16");
-
- _h_pTlead_R06_60_ratio = bookScatter2D(10, 1, 1);
+ //
+ _h_pTlead_R06_60_ratio = bookScatter2D(10, 1, 1);
_h_pTlead_R06_80_ratio = bookScatter2D(11, 1, 1);
_h_pTlead_R06_110_ratio = bookScatter2D(12, 1, 1);
_h_pTlead_R04_60_ratio = bookScatter2D(13, 1, 1);
@@ -79,6 +52,26 @@
_h_pTlead_R04_110_ratio = bookScatter2D(15, 1, 1);
_h_HT2_R06_ratio = bookScatter2D(16, 1, 1);
_h_HT2_R04_ratio = bookScatter2D(17, 1, 1);
+
+ // Temporary histograms to be divided for the dsigma3/dsigma2 ratios
+ _h_tmp_pTlead_R06_60_2 = Histo1D(refData(10, 1, 1));
+ _h_tmp_pTlead_R06_80_2 = Histo1D(refData(11, 1, 1));
+ _h_tmp_pTlead_R06_110_2 = Histo1D(refData(12, 1, 1));
+ _h_tmp_pTlead_R06_60_3 = Histo1D(refData(10, 1, 1));
+ _h_tmp_pTlead_R06_80_3 = Histo1D(refData(11, 1, 1));
+ _h_tmp_pTlead_R06_110_3 = Histo1D(refData(12, 1, 1));
+ //
+ _h_tmp_pTlead_R04_60_2 = Histo1D(refData(13, 1, 1));
+ _h_tmp_pTlead_R04_80_2 = Histo1D(refData(14, 1, 1));
+ _h_tmp_pTlead_R04_110_2 = Histo1D(refData(15, 1, 1));
+ _h_tmp_pTlead_R04_60_3 = Histo1D(refData(13, 1, 1));
+ _h_tmp_pTlead_R04_80_3 = Histo1D(refData(14, 1, 1));
+ _h_tmp_pTlead_R04_110_3 = Histo1D(refData(15, 1, 1));
+ //
+ _h_tmp_HT2_R06_2 = Histo1D(refData(16, 1, 1));
+ _h_tmp_HT2_R06_3 = Histo1D(refData(16, 1, 1));
+ _h_tmp_HT2_R04_2 = Histo1D(refData(17, 1, 1));
+ _h_tmp_HT2_R04_3 = Histo1D(refData(17, 1, 1));
}
@@ -88,64 +81,65 @@
vector<FourMomentum> jets04;
foreach (const Jet& jet, applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(60.0*GeV)) {
- if (fabs(jet.momentum().eta())<2.8) {
+ if (fabs(jet.momentum().eta()) < 2.8) {
jets04.push_back(jet.momentum());
}
}
- if (jets04.size()>1 && jets04[0].pT()>80.0*GeV) {
- for (size_t i=2; i<=jets04.size(); ++i) {
+ if (jets04.size() > 1 && jets04[0].pT() > 80.0*GeV) {
+ for (size_t i = 2; i <= jets04.size(); ++i) {
_h_jet_multi_inclusive->fill(i, weight);
}
- double HT=0.0;
- for (size_t i=0; i<jets04.size(); ++i) {
- if (i<_h_jet_pT.size()) _h_jet_pT[i]->fill(jets04[i].pT(), weight);
+ double HT = 0.0;
+ for (size_t i = 0; i < jets04.size(); ++i) {
+ if (i < _h_jet_pT.size()) _h_jet_pT[i]->fill(jets04[i].pT(), weight);
HT += jets04[i].pT();
}
- if (jets04.size()>=2) _h_HT_2->fill(HT, weight);
- if (jets04.size()>=3) _h_HT_3->fill(HT, weight);
- if (jets04.size()>=4) _h_HT_4->fill(HT, weight);
+ if (jets04.size() >= 2) _h_HT_2->fill(HT, weight);
+ if (jets04.size() >= 3) _h_HT_3->fill(HT, weight);
+ if (jets04.size() >= 4) _h_HT_4->fill(HT, weight);
double pT1(jets04[0].pT()), pT2(jets04[1].pT());
- double HT2=pT1+pT2;
- if (jets04.size()>=2) {
- _h_tmp_HT2_R04_2->fill(HT2, weight);
- _h_tmp_pTlead_R04_60_2->fill(pT1, weight);
- if (pT2>80.0*GeV) _h_tmp_pTlead_R04_80_2->fill(pT1, weight);
- if (pT2>110.0*GeV) _h_tmp_pTlead_R04_110_2->fill(pT1, weight);
+ double HT2 = pT1 + pT2;
+ if (jets04.size() >= 2) {
+ _h_tmp_HT2_R04_2.fill(HT2, weight);
+ _h_tmp_pTlead_R04_60_2.fill(pT1, weight);
+ if (pT2 > 80.0*GeV) _h_tmp_pTlead_R04_80_2.fill(pT1, weight);
+ if (pT2 > 110.0*GeV) _h_tmp_pTlead_R04_110_2.fill(pT1, weight);
}
- if (jets04.size()>=3) {
+ if (jets04.size() >= 3) {
double pT3(jets04[2].pT());
- _h_tmp_HT2_R04_3->fill(HT2, weight);
- _h_tmp_pTlead_R04_60_3->fill(pT1, weight);
- if (pT3>80.0*GeV) _h_tmp_pTlead_R04_80_3->fill(pT1, weight);
- if (pT3>110.0*GeV) _h_tmp_pTlead_R04_110_3->fill(pT1, weight);
+ _h_tmp_HT2_R04_3.fill(HT2, weight);
+ _h_tmp_pTlead_R04_60_3.fill(pT1, weight);
+ if (pT3 > 80.0*GeV) _h_tmp_pTlead_R04_80_3.fill(pT1, weight);
+ if (pT3 > 110.0*GeV) _h_tmp_pTlead_R04_110_3.fill(pT1, weight);
}
}
+ /// @todo It'd be better to avoid duplicating 95% of the code!
vector<FourMomentum> jets06;
foreach (const Jet& jet, applyProjection<FastJets>(event, "AntiKtJets06").jetsByPt(60.0*GeV)) {
- if (fabs(jet.momentum().eta())<2.8) {
+ if (fabs(jet.momentum().eta()) < 2.8) {
jets06.push_back(jet.momentum());
}
}
- if (jets06.size()>1 && jets06[0].pT()>80.0*GeV) {
+ if (jets06.size() > 1 && jets06[0].pT() > 80.0*GeV) {
double pT1(jets06[0].pT()), pT2(jets06[1].pT());
- double HT2=pT1+pT2;
- if (jets06.size()>=2) {
- _h_tmp_HT2_R06_2->fill(HT2, weight);
- _h_tmp_pTlead_R06_60_2->fill(pT1, weight);
- if (pT2>80.0*GeV) _h_tmp_pTlead_R06_80_2->fill(pT1, weight);
- if (pT2>110.0*GeV) _h_tmp_pTlead_R06_110_2->fill(pT1, weight);
+ double HT2 = pT1 + pT2;
+ if (jets06.size() >= 2) {
+ _h_tmp_HT2_R06_2.fill(HT2, weight);
+ _h_tmp_pTlead_R06_60_2.fill(pT1, weight);
+ if (pT2 > 80.0*GeV) _h_tmp_pTlead_R06_80_2.fill(pT1, weight);
+ if (pT2 > 110.0*GeV) _h_tmp_pTlead_R06_110_2.fill(pT1, weight);
}
- if (jets06.size()>=3) {
+ if (jets06.size() >= 3) {
double pT3(jets06[2].pT());
- _h_tmp_HT2_R06_3->fill(HT2, weight);
- _h_tmp_pTlead_R06_60_3->fill(pT1, weight);
- if (pT3>80.0*GeV) _h_tmp_pTlead_R06_80_3->fill(pT1, weight);
- if (pT3>110.0*GeV) _h_tmp_pTlead_R06_110_3->fill(pT1, weight);
+ _h_tmp_HT2_R06_3.fill(HT2, weight);
+ _h_tmp_pTlead_R06_60_3.fill(pT1, weight);
+ if (pT3 > 80.0*GeV) _h_tmp_pTlead_R06_80_3.fill(pT1, weight);
+ if (pT3 > 110.0*GeV) _h_tmp_pTlead_R06_110_3.fill(pT1, weight);
}
}
@@ -155,21 +149,16 @@
/// Normalise histograms etc., after the run
void finalize() {
- // Fill inclusive jet multi ratio
- int Nbins = _h_jet_multi_inclusive->numBins();
- std::vector<double> ratio(Nbins-1, 0.0);
- std::vector<double> err(Nbins-1, 0.0);
- for (int i = 0; i < Nbins-1; ++i) {
- if (_h_jet_multi_inclusive->bin(i).area() > 0.0 && _h_jet_multi_inclusive->bin(i+1).area() > 0.0) {
- ratio[i] = _h_jet_multi_inclusive->bin(i+1).area()/_h_jet_multi_inclusive->bin(i).area();
- double relerr_i = _h_jet_multi_inclusive->bin(i).areaErr()/_h_jet_multi_inclusive->bin(i).area();
- double relerr_j = _h_jet_multi_inclusive->bin(i+1).areaErr()/_h_jet_multi_inclusive->bin(i+1).area();
- err[i] = ratio[i] * (relerr_i + relerr_j);
+ // Fill inclusive jet multiplicity ratio
+ for (size_t i = 0; i < _h_jet_multi_ratio->numPoints(); ++i) {
+ if (_h_jet_multi_inclusive->bin(i).sumW() != 0) {
+ const double val = _h_jet_multi_inclusive->bin(i+1).sumW() / _h_jet_multi_inclusive->bin(i).sumW();
+ const double err = ( _h_jet_multi_inclusive->bin(i+1).relErr() + _h_jet_multi_inclusive->bin(i).relErr() ) * val;
+ _h_jet_multi_ratio->point(i).setY(val, err);
}
}
- //\todo YODA
- //_h_jet_multi_ratio->setCoordinate(1, ratio, err);
+ // Normalize std histos
scale(_h_jet_multi_inclusive, crossSectionPerEvent());
scale(_h_jet_pT[0], crossSectionPerEvent());
scale(_h_jet_pT[1], crossSectionPerEvent());
@@ -179,31 +168,15 @@
scale(_h_HT_3, crossSectionPerEvent());
scale(_h_HT_4, crossSectionPerEvent());
- /// create ratio histograms
- divide(_h_tmp_pTlead_R06_60_3,_h_tmp_pTlead_R06_60_2,
- _h_pTlead_R06_60_ratio);
-
- divide(_h_tmp_pTlead_R06_80_3,_h_tmp_pTlead_R06_80_2,
- _h_pTlead_R06_80_ratio);
-
- divide(_h_tmp_pTlead_R06_110_3,_h_tmp_pTlead_R06_110_2,
- _h_pTlead_R06_110_ratio);
-
- divide(_h_tmp_pTlead_R04_60_3,_h_tmp_pTlead_R04_60_2,
- _h_pTlead_R04_60_ratio);
-
- divide(_h_tmp_pTlead_R04_80_3,_h_tmp_pTlead_R04_80_2,
- _h_pTlead_R04_80_ratio);
-
- divide(_h_tmp_pTlead_R04_110_3,_h_tmp_pTlead_R04_110_2,
- _h_pTlead_R04_110_ratio);
-
- divide(_h_tmp_HT2_R06_3,_h_tmp_HT2_R06_2,
- _h_HT2_R06_ratio);
-
- divide(_h_tmp_HT2_R04_3,_h_tmp_HT2_R04_2,
- _h_HT2_R04_ratio);
-
+ /// Create ratio histograms
+ divide(_h_tmp_pTlead_R06_60_3,_h_tmp_pTlead_R06_60_2, _h_pTlead_R06_60_ratio);
+ divide(_h_tmp_pTlead_R06_80_3,_h_tmp_pTlead_R06_80_2, _h_pTlead_R06_80_ratio);
+ divide(_h_tmp_pTlead_R06_110_3,_h_tmp_pTlead_R06_110_2, _h_pTlead_R06_110_ratio);
+ divide(_h_tmp_pTlead_R04_60_3,_h_tmp_pTlead_R04_60_2, _h_pTlead_R04_60_ratio);
+ divide(_h_tmp_pTlead_R04_80_3,_h_tmp_pTlead_R04_80_2, _h_pTlead_R04_80_ratio);
+ divide(_h_tmp_pTlead_R04_110_3,_h_tmp_pTlead_R04_110_2, _h_pTlead_R04_110_ratio);
+ divide(_h_tmp_HT2_R06_3,_h_tmp_HT2_R06_2, _h_HT2_R06_ratio);
+ divide(_h_tmp_HT2_R04_3,_h_tmp_HT2_R04_2, _h_HT2_R04_ratio);
}
//@}
@@ -211,11 +184,6 @@
private:
- // Data members like post-cuts event weight counters go here
-
-
- private:
-
/// @name Histograms
//@{
Histo1DPtr _h_jet_multi_inclusive;
@@ -224,36 +192,22 @@
Histo1DPtr _h_HT_2;
Histo1DPtr _h_HT_3;
Histo1DPtr _h_HT_4;
+ //@}
- /// temporary histograms which will be divided in the end for the dsigma3/dsigma2 ratios
- Histo1DPtr _h_tmp_pTlead_R06_60_2;
- Histo1DPtr _h_tmp_pTlead_R06_80_2;
- Histo1DPtr _h_tmp_pTlead_R06_110_2;
- Histo1DPtr _h_tmp_pTlead_R06_60_3;
- Histo1DPtr _h_tmp_pTlead_R06_80_3;
- Histo1DPtr _h_tmp_pTlead_R06_110_3;
-
- Histo1DPtr _h_tmp_pTlead_R04_60_2;
- Histo1DPtr _h_tmp_pTlead_R04_80_2;
- Histo1DPtr _h_tmp_pTlead_R04_110_2;
- Histo1DPtr _h_tmp_pTlead_R04_60_3;
- Histo1DPtr _h_tmp_pTlead_R04_80_3;
- Histo1DPtr _h_tmp_pTlead_R04_110_3;
-
- Histo1DPtr _h_tmp_HT2_R06_2;
- Histo1DPtr _h_tmp_HT2_R06_3;
- Histo1DPtr _h_tmp_HT2_R04_2;
- Histo1DPtr _h_tmp_HT2_R04_3;
-
- Scatter2DPtr _h_pTlead_R06_60_ratio;
- Scatter2DPtr _h_pTlead_R06_80_ratio;
- Scatter2DPtr _h_pTlead_R06_110_ratio;
- Scatter2DPtr _h_pTlead_R04_60_ratio;
- Scatter2DPtr _h_pTlead_R04_80_ratio;
- Scatter2DPtr _h_pTlead_R04_110_ratio;
- Scatter2DPtr _h_HT2_R06_ratio;
- Scatter2DPtr _h_HT2_R04_ratio;
+ /// @name Ratio histograms
+ //@{
+ Scatter2DPtr _h_pTlead_R06_60_ratio, _h_pTlead_R06_80_ratio, _h_pTlead_R06_110_ratio;
+ Scatter2DPtr _h_pTlead_R04_60_ratio, _h_pTlead_R04_80_ratio, _h_pTlead_R04_110_ratio;
+ Scatter2DPtr _h_HT2_R06_ratio, _h_HT2_R04_ratio;
+ //@}
+ /// @name Temporary histograms to be divided for the dsigma3/dsigma2 ratios
+ //@{
+ Histo1D _h_tmp_pTlead_R06_60_2, _h_tmp_pTlead_R06_80_2, _h_tmp_pTlead_R06_110_2;
+ Histo1D _h_tmp_pTlead_R06_60_3, _h_tmp_pTlead_R06_80_3, _h_tmp_pTlead_R06_110_3;
+ Histo1D _h_tmp_pTlead_R04_60_2, _h_tmp_pTlead_R04_80_2, _h_tmp_pTlead_R04_110_2;
+ Histo1D _h_tmp_pTlead_R04_60_3, _h_tmp_pTlead_R04_80_3, _h_tmp_pTlead_R04_110_3;
+ Histo1D _h_tmp_HT2_R06_2, _h_tmp_HT2_R06_3, _h_tmp_HT2_R04_2, _h_tmp_HT2_R04_3;
//@}
};
Modified: trunk/src/Core/Analysis.cc
==============================================================================
--- trunk/src/Core/Analysis.cc Tue May 14 17:52:45 2013 (r4275)
+++ trunk/src/Core/Analysis.cc Tue May 14 18:54:14 2013 (r4276)
@@ -191,17 +191,20 @@
}
- Histo1DPtr Analysis::bookHisto1D(size_t datasetId, size_t xAxisId,
- size_t yAxisId, const string& title,
- const string& xtitle, const string& ytitle)
+ Histo1DPtr Analysis::bookHisto1D(size_t datasetId, size_t xAxisId, size_t yAxisId,
+ const string& title,
+ const string& xtitle,
+ const string& ytitle)
{
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookHisto1D(axisCode, title, xtitle, ytitle);
}
- Histo1DPtr Analysis::bookHisto1D(const string& hname, const string& title,
- const string& xtitle, const string& ytitle)
+ Histo1DPtr Analysis::bookHisto1D(const string& hname,
+ const string& title,
+ const string& xtitle,
+ const string& ytitle)
{
const Scatter2D & refdata = refData(hname);
const string path = histoPath(hname);
@@ -218,7 +221,8 @@
Histo1DPtr Analysis::bookHisto1D(const string& hname,
size_t nbins, double lower, double upper,
const string& title,
- const string& xtitle, const string& ytitle) {
+ const string& xtitle,
+ const string& ytitle) {
const string path = histoPath(hname);
Histo1DPtr hist( new Histo1D(nbins, lower, upper, path, title) );
addPlot(hist);
More information about the Rivet-svn
mailing list
|