[Rivet-svn] r2134 - in trunk: . src/Analyses src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Fri Dec  4 15:53:25 2009
New Revision: 2134

Log:
Correcting norms of DELPHI_1996 and tidying Hemispheres implementation: hemisphere masses are still not right for OPAL_2004 on the Z pole, but *are* correct for DELPHI_1996 on the Z... it uses M^2/Evis^2 rather than OPAL's (I think) M/sqrt(s). Hmm. Can someone else please take a look and see if I'm missing something obvious from the OPAL observable definition? The shape is now ~right, but squeezed to smaller values of M/sqrt(s) than they should be.

Modified:
   trunk/ChangeLog
   trunk/src/Analyses/DELPHI_1996_S3430090.cc
   trunk/src/Analyses/OPAL_2004_S6132243.cc
   trunk/src/Projections/Hemispheres.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Fri Dec  4 14:14:36 2009	(r2133)
+++ trunk/ChangeLog	Fri Dec  4 15:53:25 2009	(r2134)
@@ -1,5 +1,12 @@
 2009-12-04  Andy Buckley  <andy at insectnation.org>
 
+	* Use scaling rather than normalising in DELPHI_1996: this is
+	generally desirable, since normalizing to 1 for 1/sig dsig/dx
+	observables isn't correct if any events fall outside the histo
+	bounds.
+
+	* Many fixes to OPAL_2004.
+
 	* Improved Hemispheres interface to remove unnecessary consts on
 	returned doubles, and to also return non-squared versions
 	of (scaled) hemisphere masses.

Modified: trunk/src/Analyses/DELPHI_1996_S3430090.cc
==============================================================================
--- trunk/src/Analyses/DELPHI_1996_S3430090.cc	Fri Dec  4 14:14:36 2009	(r2133)
+++ trunk/src/Analyses/DELPHI_1996_S3430090.cc	Fri Dec  4 15:53:25 2009	(r2134)
@@ -40,8 +40,10 @@
       : Analysis("DELPHI_1996_S3430090")
     {
       setBeams(ELECTRON, POSITRON);
-      _weightedTotalPartNum = 0;
-      _passedCutWeightSum = 0;
+      _weightedTotalPartNum = 0.0;
+      _passedCutWeightSum = 0.0;
+      _passedCut3WeightSum = 0.0;
+      _passedCut4WeightSum = 0.0;
     }
  
  
@@ -173,16 +175,20 @@
    
       // Jets
       const FastJets& durjet = applyProjection<FastJets>(e, "DurhamJets");
-      if (durjet.clusterSeq()) {
-        _histDiffRate2Durham->fill(durjet.clusterSeq()->exclusive_ymerge(2), weight);
-        _histDiffRate3Durham->fill(durjet.clusterSeq()->exclusive_ymerge(3), weight);
-        _histDiffRate4Durham->fill(durjet.clusterSeq()->exclusive_ymerge(4), weight);
-      }
       const FastJets& jadejet = applyProjection<FastJets>(e, "JadeJets");
-      if (jadejet.clusterSeq()) {
+      if (durjet.clusterSeq() && jadejet.clusterSeq()) {
+        _histDiffRate2Durham->fill(durjet.clusterSeq()->exclusive_ymerge(2), weight);
         _histDiffRate2Jade->fill(jadejet.clusterSeq()->exclusive_ymerge(2), weight);
-        _histDiffRate3Jade->fill(jadejet.clusterSeq()->exclusive_ymerge(3), weight);
-        _histDiffRate4Jade->fill(jadejet.clusterSeq()->exclusive_ymerge(4), weight);
+        if (numParticles >= 3) {
+          _passedCut3WeightSum += weight;
+          _histDiffRate3Durham->fill(durjet.clusterSeq()->exclusive_ymerge(3), weight);
+          _histDiffRate3Jade->fill(jadejet.clusterSeq()->exclusive_ymerge(3), weight);
+        }
+        if (numParticles >= 4) {
+          _passedCut4WeightSum += weight;
+          _histDiffRate4Durham->fill(durjet.clusterSeq()->exclusive_ymerge(4), weight);
+          _histDiffRate4Jade->fill(jadejet.clusterSeq()->exclusive_ymerge(4), weight);
+        }
       }
    
       // Sphericities
@@ -410,33 +416,33 @@
       scale(_histMultiXi1530_0, 1.0/_passedCutWeightSum);
       scale(_histMultiLambdaB0, 1.0/_passedCutWeightSum);
 
-      normalize(_hist1MinusT);
-      normalize(_histTMajor);
-      normalize(_histTMinor);
-      normalize(_histOblateness);
-
-      normalize(_histSphericity);
-      normalize(_histAplanarity);
-      normalize(_histPlanarity);
-
-      normalize(_histHemiMassD);
-      normalize(_histHemiMassH);
-      normalize(_histHemiMassL);
-
-      normalize(_histHemiBroadW);
-      normalize(_histHemiBroadN);
-      normalize(_histHemiBroadT);
-      normalize(_histHemiBroadD);
-
-      normalize(_histCParam);
-      normalize(_histDParam);
-
-      normalize(_histDiffRate2Durham);
-      normalize(_histDiffRate2Jade);
-      normalize(_histDiffRate3Durham);
-      normalize(_histDiffRate3Jade);
-      normalize(_histDiffRate4Durham);
-      normalize(_histDiffRate4Jade);
+      scale(_hist1MinusT, 1.0/_passedCutWeightSum);
+      scale(_histTMajor, 1.0/_passedCutWeightSum);
+      scale(_histTMinor, 1.0/_passedCutWeightSum);
+      scale(_histOblateness, 1.0/_passedCutWeightSum);
+
+      scale(_histSphericity, 1.0/_passedCutWeightSum);
+      scale(_histAplanarity, 1.0/_passedCutWeightSum);
+      scale(_histPlanarity, 1.0/_passedCutWeightSum);
+
+      scale(_histHemiMassD, 1.0/_passedCutWeightSum);
+      scale(_histHemiMassH, 1.0/_passedCutWeightSum);
+      scale(_histHemiMassL, 1.0/_passedCutWeightSum);
+
+      scale(_histHemiBroadW, 1.0/_passedCutWeightSum);
+      scale(_histHemiBroadN, 1.0/_passedCutWeightSum);
+      scale(_histHemiBroadT, 1.0/_passedCutWeightSum);
+      scale(_histHemiBroadD, 1.0/_passedCutWeightSum);
+
+      scale(_histCParam, 1.0/_passedCutWeightSum);
+      scale(_histDParam, 1.0/_passedCutWeightSum);
+
+      scale(_histDiffRate2Durham, 1.0/_passedCutWeightSum);
+      scale(_histDiffRate2Jade, 1.0/_passedCutWeightSum);
+      scale(_histDiffRate3Durham, 1.0/_passedCut3WeightSum);
+      scale(_histDiffRate3Jade, 1.0/_passedCut3WeightSum);
+      scale(_histDiffRate4Durham, 1.0/_passedCut4WeightSum);
+      scale(_histDiffRate4Jade, 1.0/_passedCut4WeightSum);
     }
 
     //@}
@@ -449,7 +455,12 @@
     /// inclusive single particle distributions' normalisations.
     double _weightedTotalPartNum;
 
+    /// @name Sums of weights past various cuts
+    //@{
     double _passedCutWeightSum;
+    double _passedCut3WeightSum;
+    double _passedCut4WeightSum;
+    //@}
 
     /// @name Histograms
     //@{

Modified: trunk/src/Analyses/OPAL_2004_S6132243.cc
==============================================================================
--- trunk/src/Analyses/OPAL_2004_S6132243.cc	Fri Dec  4 14:14:36 2009	(r2133)
+++ trunk/src/Analyses/OPAL_2004_S6132243.cc	Fri Dec  4 15:53:25 2009	(r2134)
@@ -152,7 +152,9 @@
    
       // Hemispheres
       const Hemispheres& hemi = applyProjection<Hemispheres>(event, "Hemispheres");
+      /// @todo Something wrong here: values are too low
       const double hemi_mh = hemi.Mhigh()/sqrtS();
+      /// @todo Something wrong here: values are too low
       const double hemi_ml = hemi.Mlow()/sqrtS();
       const double hemi_bmax = hemi.Bmax();
       const double hemi_bmin = hemi.Bmin();

Modified: trunk/src/Projections/Hemispheres.cc
==============================================================================
--- trunk/src/Projections/Hemispheres.cc	Fri Dec  4 14:14:36 2009	(r2133)
+++ trunk/src/Projections/Hemispheres.cc	Fri Dec  4 15:53:25 2009	(r2134)
@@ -13,7 +13,7 @@
     FourMomentum p4With, p4Against;
     double Evis(0), broadWith(0), broadAgainst(0), broadDenom(0);
     const FinalState& fs = applyProjection<FinalState>(e, ax.getProjection("FS"));
-    const ParticleVector particles = fs.particles();
+    const ParticleVector& particles = fs.particles();
     getLog() << Log::DEBUG << "number of particles = " << particles.size() << endl;
     foreach (const Particle& p, particles) {
       const FourMomentum p4 = p.momentum();
@@ -48,21 +48,21 @@
     _E2vis = Evis * Evis;
 
     // Calculate masses.
-    const double mass2With = p4With.invariant();
-    const double mass2Against = p4Against.invariant();
-    const bool withIsMaxMass = (mass2With > mass2Against);
-    _M2high = (withIsMaxMass) ? mass2With : mass2Against;
-    _M2low = (withIsMaxMass) ? mass2Against : mass2With;
+    const double mass2With = p4With.mass2();
+    const double mass2Against = p4Against.mass2();
+    _M2high = max(mass2With, mass2Against);
+    _M2low = min(mass2With, mass2Against);
 
     // Calculate broadenings.
     broadWith /= broadDenom;
     broadAgainst /= broadDenom;
-    const bool withIsMaxBroad = (broadWith > broadAgainst);
-    _Bmax = (withIsMaxBroad) ? broadWith : broadAgainst;
-    _Bmin = (withIsMaxBroad) ? broadAgainst : broadWith;
+    _Bmax = max(broadWith, broadAgainst);
+    _Bmin = min(broadWith, broadAgainst);
 
     // Calculate high-max correlation flag.
-    _highMassEqMaxBroad = ((withIsMaxMass && withIsMaxBroad) || (!withIsMaxMass && !withIsMaxBroad));
+    const int maxMassID = (mass2With >= mass2Against);
+    const int maxBroadID = (broadWith >= broadAgainst);
+    _highMassEqMaxBroad = (maxMassID == maxBroadID);
   }