[Rivet-svn] r2455 - trunk/src/Analyses

[blackhole at projects.hepforge.org]

Author: hoeth
Date: Tue May 18 14:20:54 2010
New Revision: 2455

Log:
revert last commit for D0_2001_S4674421.

Modified:
   trunk/src/Analyses/D0_2001_S4674421.cc

Modified: trunk/src/Analyses/D0_2001_S4674421.cc
==============================================================================
--- trunk/src/Analyses/D0_2001_S4674421.cc	Tue May 18 14:18:26 2010	(r2454)
+++ trunk/src/Analyses/D0_2001_S4674421.cc	Tue May 18 14:20:54 2010	(r2455)
@@ -28,14 +28,6 @@
     //@{
  
     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");
@@ -67,14 +59,15 @@
       // Histograms
       _h_dsigdpt_w = bookHistogram1D(1, 1, 1);
       _h_dsigdpt_z = bookHistogram1D(1, 1, 2);
-      _h_dsigdpt_scaled_z = bookDataPointSet(2, 1, 1);
-#endif
+      /// @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);
     }
 
 
 
     void analyze(const Event& event) {
-#if 0
       const double weight = event.weight();
 
       const LeadingParticlesFinalState& eeFS = applyProjection<LeadingParticlesFinalState>(event, "eeFS");
@@ -89,7 +82,9 @@
           ++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?
@@ -112,13 +107,11 @@
           _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;
 
@@ -129,64 +122,33 @@
       const double xSecZ = xSecPerEvent * _eventsFilledZ;
 
       // Get W and Z pT integrals
-      const double wpt_integral = xSecW; //integral(_h_dsigdpt_w);
-      const double zpt_integral = xSecZ; //integral(_h_dsigdpt_z);
+      const double wpt_integral = integral(_h_dsigdpt_w);
+      const double zpt_scaled_integral = integral(_h_dsigdpt_scaled_z);
 
       // Divide and scale ratio histos
-      if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_integral == 0) {
+      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) {
         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.
-        // 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.);
+        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);
         }
-        _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);
     }
 
 
@@ -196,16 +158,15 @@
 
     /// @name Event counters for cross section normalizations
     //@{
-///    double _eventsFilledW;
-///    double _eventsFilledZ;
+    double _eventsFilledW;
+    double _eventsFilledZ;
     //@}
 
     //@{
     /// Histograms
-    AIDA::IHistogram1D*  ttt;
-///    AIDA::IHistogram1D*  _h_dsigdpt_w;
-///    AIDA::IHistogram1D*  _h_dsigdpt_z;
-///    AIDA::IDataPointSet* _h_dsigdpt_scaled_z;
+    AIDA::IHistogram1D* _h_dsigdpt_w;
+    AIDA::IHistogram1D* _h_dsigdpt_z;
+    AIDA::IHistogram1D* _h_dsigdpt_scaled_z;
     //@}
 
   };