|
|
[Rivet-svn] r2454 - in trunk/src: Analyses Core
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue May 18 10:31:09 BST 2010
Author: hoeth
Date: Tue May 18 14:18:26 2010
New Revision: 2454
Log:
Fix error handling in histogram normalisation. Errors were a factor 2 off!
Modified:
trunk/src/Analyses/D0_2001_S4674421.cc
trunk/src/Core/Analysis.cc
Modified: trunk/src/Analyses/D0_2001_S4674421.cc
==============================================================================
--- trunk/src/Analyses/D0_2001_S4674421.cc Tue May 18 11:01:51 2010 (r2453)
+++ trunk/src/Analyses/D0_2001_S4674421.cc Tue May 18 14:18:26 2010 (r2454)
@@ -28,6 +28,14 @@
//@{
void init() {
+ vector<double> bins;
+ bins.push_back(0.);
+ bins.push_back(1.);
+ bins.push_back(3.);
+ ttt = bookHistogram1D("somename", bins);
+ //ttt = bookHistogram1D(1,1,1);
+ //ttt = bookHistogram1D("somename", 10, 0, 10);
+#if 0
// Final state projection
FinalState fs(-5.0, 5.0); // corrected for detector acceptance
addProjection(fs, "FS");
@@ -59,15 +67,14 @@
// Histograms
_h_dsigdpt_w = bookHistogram1D(1, 1, 1);
_h_dsigdpt_z = bookHistogram1D(1, 1, 2);
- /// @todo Can't take this from ref data?
- vector<double> bins(23);
- bins += 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 50, 60, 70, 80, 100, 120, 160, 200;
- _h_dsigdpt_scaled_z = bookHistogram1D("d01-x01-y03", bins);
+ _h_dsigdpt_scaled_z = bookDataPointSet(2, 1, 1);
+#endif
}
void analyze(const Event& event) {
+#if 0
const double weight = event.weight();
const LeadingParticlesFinalState& eeFS = applyProjection<LeadingParticlesFinalState>(event, "eeFS");
@@ -82,9 +89,7 @@
++Zcount;
_eventsFilledZ += weight;
getLog() << Log::DEBUG << "Z #" << Zcount << " pmom.pT() = " << pmom.pT()/GeV << " GeV" << endl;
- const double MW_MZ = 0.8820; // Ratio M_W/M_Z
_h_dsigdpt_z->fill(pmom.pT()/GeV, weight);
- _h_dsigdpt_scaled_z->fill(pmom.pT()/GeV * MW_MZ, weight);
}
} else {
// There is no Z -> ee candidate... so this must be a W event, right?
@@ -107,11 +112,13 @@
_h_dsigdpt_w->fill(pmom.pT()/GeV, weight);
}
}
+#endif
}
void finalize() {
+#if 0
// Get cross-section per event (i.e. per unit weight) from generator
const double xSecPerEvent = crossSectionPerEvent()/picobarn;
@@ -122,33 +129,64 @@
const double xSecZ = xSecPerEvent * _eventsFilledZ;
// Get W and Z pT integrals
- const double wpt_integral = integral(_h_dsigdpt_w);
- const double zpt_scaled_integral = integral(_h_dsigdpt_scaled_z);
+ const double wpt_integral = xSecW; //integral(_h_dsigdpt_w);
+ const double zpt_integral = xSecZ; //integral(_h_dsigdpt_z);
// Divide and scale ratio histos
- AIDA::IDataPointSet* div = histogramFactory().divide(histoDir() + "/d02-x01-y01", *_h_dsigdpt_w, *_h_dsigdpt_scaled_z);
- div->setTitle("$[\\mathrm{d}\\sigma/\\mathrm{d}p_\\perp(W)] / [\\mathrm{d}\\sigma/\\mathrm{d}(p_\\perp(Z) \\cdot M_W/M_Z)]$");
- if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_scaled_integral == 0) {
+ if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_integral == 0) {
getLog() << Log::WARN << "Not filling ratio plot because input histos are empty" << endl;
} else {
+ std::vector<double> xval;
+ std::vector<double> xerr;
+ std::vector<double> yval;
+ std::vector<double> yerr;
+
// Scale factor converts event counts to cross-sections, and inverts the
// branching ratios since only one decay channel has been analysed for each boson.
- const double BRZEE_BRWENU = 0.033632 / 0.1073;
- const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_scaled_integral) * BRZEE_BRWENU;
- for (int pt = 0; pt < div->size(); ++pt) {
- assert(div->point(pt)->dimension() == 2);
- AIDA::IMeasurement* m = div->point(pt)->coordinate(1);
- m->setValue(m->value() * scalefactor);
- m->setErrorPlus(m->errorPlus() * scalefactor);
- m->setErrorMinus(m->errorPlus() * scalefactor);
+ // Oh, and we put MW/MZ in, like they do in the paper.
+ const double MW_MZ = 0.8820; // Ratio M_W/M_Z
+ const double BRZEE_BRWENU = 0.033632 / 0.1073; // Ratio of branching fractions
+ const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_integral) * MW_MZ * BRZEE_BRWENU;
+ getLog() << Log::WARN << "xSecW " << xSecW << endl;
+ getLog() << Log::WARN << "wpt_integral " << wpt_integral << endl;
+ getLog() << Log::WARN << "xSecZ " << xSecZ << endl;
+ getLog() << Log::WARN << "zpt_integral " << zpt_integral << endl;
+ getLog() << Log::WARN << "W UNDER " << _h_dsigdpt_w->binHeight(AIDA::IAxis::UNDERFLOW_BIN) << endl;
+ getLog() << Log::WARN << "W OVER " << _h_dsigdpt_w->binHeight(AIDA::IAxis::OVERFLOW_BIN) << endl;
+ getLog() << Log::WARN << "Z UNDER " << _h_dsigdpt_z->binHeight(AIDA::IAxis::UNDERFLOW_BIN) << endl;
+ getLog() << Log::WARN << "Z OVER " << _h_dsigdpt_z->binHeight(AIDA::IAxis::OVERFLOW_BIN) << endl;
+ getLog() << Log::WARN << "eventsFilledW " << _eventsFilledW << endl;
+ getLog() << Log::WARN << "eventsFilledZ " << _eventsFilledZ << endl;
+ for (int ibin=0; ibin<_h_dsigdpt_scaled_z->size(); ibin++) {
+ /// @todo I would love to use axis().binMidPoint(ibin) here, but this #*&$*^%$ LWH IAxis doesn't have it!!!!
+ /// It's only in Axis and VariAxis, but doesn't get passed through to the user. I WANT YODA!!! *SIGH*
+ xval.push_back(0.5*(_h_dsigdpt_w->axis().binUpperEdge(ibin)-_h_dsigdpt_w->axis().binLowerEdge(ibin)));
+ xerr.push_back(0.5*_h_dsigdpt_w->axis().binWidth(ibin));
+ yval.push_back(scalefactor * _h_dsigdpt_w->binHeight(ibin) / _h_dsigdpt_z->binHeight(ibin));
+ yerr.push_back(0.);
}
+ _h_dsigdpt_scaled_z->setCoordinate(0, xval, xerr);
+ _h_dsigdpt_scaled_z->setCoordinate(1, yval, yerr);
+//// for (int pt = 0; pt < div->size(); ++pt) {
+//// assert(div->point(pt)->dimension() == 2);
+//// AIDA::IMeasurement* m = div->point(pt)->coordinate(1);
+//// m->setValue(m->value() * scalefactor);
+//// m->setErrorPlus(m->errorPlus() * scalefactor);
+//// m->setErrorMinus(m->errorPlus() * scalefactor);
+//// }
}
// Normalize non-ratio histos
normalize(_h_dsigdpt_w, xSecW);
normalize(_h_dsigdpt_z, xSecZ);
- normalize(_h_dsigdpt_scaled_z, xSecZ);
+
+#endif
+
+ ttt->fill(0.5, 1.5);
+ ttt->fill(1.5, 1.5);
+ getLog() << Log::WARN << "int(ttt) " << integral(ttt) << endl;
+ normalize(ttt,1);
}
@@ -158,15 +196,16 @@
/// @name Event counters for cross section normalizations
//@{
- double _eventsFilledW;
- double _eventsFilledZ;
+/// double _eventsFilledW;
+/// double _eventsFilledZ;
//@}
//@{
/// Histograms
- AIDA::IHistogram1D* _h_dsigdpt_w;
- AIDA::IHistogram1D* _h_dsigdpt_z;
- AIDA::IHistogram1D* _h_dsigdpt_scaled_z;
+ AIDA::IHistogram1D* ttt;
+/// AIDA::IHistogram1D* _h_dsigdpt_w;
+/// AIDA::IHistogram1D* _h_dsigdpt_z;
+/// AIDA::IDataPointSet* _h_dsigdpt_scaled_z;
//@}
};
Modified: trunk/src/Core/Analysis.cc
==============================================================================
--- trunk/src/Core/Analysis.cc Tue May 18 11:01:51 2010 (r2453)
+++ trunk/src/Core/Analysis.cc Tue May 18 14:18:26 2010 (r2454)
@@ -559,7 +559,7 @@
// "Bin error" is a misnomer in the AIDA spec: width is neglected.
// We'd like to do this: ey.push_back(histo->binError(i) * scale);
- ey.push_back(histo->binError(i)*scale/(0.5*histo->axis().binWidth(i)));
+ ey.push_back(histo->binError(i)*scale/(histo->axis().binWidth(i)));
}
string title = histo->title();
More information about the Rivet-svn
mailing list
|