|
|
[Rivet-svn] r2302 - trunk/src/Analyses
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Wed Mar 3 10:24:55 GMT 2010
Author: fsiegert
Date: Wed Mar 3 10:24:55 2010
New Revision: 2302
Log:
Bugfix in CDF_1996_S3418421.
Improvements in CDF_1996_S3349578 and CDF_1997_S3541940: After boosting the
jets around the numerics is not stable enough to guarantee mass2()>-1e8. Use
_safeMass() to avoid running into that assert(). Also avoid any boosting of
the average beam vector in case the jet system rapidity is \approx 0 (where
no boosting is necessary anymore). This should make that part numerically more
stable.
Modified:
trunk/src/Analyses/CDF_1996_S3349578.cc
trunk/src/Analyses/CDF_1996_S3418421.cc
trunk/src/Analyses/CDF_1997_S3541940.cc
Modified: trunk/src/Analyses/CDF_1996_S3349578.cc
==============================================================================
--- trunk/src/Analyses/CDF_1996_S3349578.cc Tue Mar 2 22:58:09 2010 (r2301)
+++ trunk/src/Analyses/CDF_1996_S3349578.cc Wed Mar 3 10:24:55 2010 (r2302)
@@ -127,7 +127,7 @@
}
if (sumEt < 420.0*GeV) return;
- const double m3J = jetsystem.mass();
+ const double m3J = _safeMass(jetsystem);
if (m3J<600*GeV) {
return;
}
@@ -155,9 +155,9 @@
const double X4 = 2.0*p4.E()/m3J;
const double psi3 = _psi(p3, pAV, p4, p5);
- const double f3 = p3.mass()/m3J;
- const double f4 = p4.mass()/m3J;
- const double f5 = p5.mass()/m3J;
+ const double f3 = _safeMass(p3)/m3J;
+ const double f4 = _safeMass(p4)/m3J;
+ const double f5 = _safeMass(p5)/m3J;
_h_3_mNJ->fill(m3J, weight);
_h_3_X3->fill(X3, weight);
@@ -184,7 +184,7 @@
}
if (sumEt < 420.0*GeV) return;
- const double m4J = jetsystem.mass();
+ const double m4J = _safeMass(jetsystem);
if (m4J < 650*GeV) return;
LorentzTransform cms_boost(-jetsystem.boostVector());
@@ -215,11 +215,11 @@
// fill histograms
const double X4 = 2.0*p4.E()/m4J;
const double psi3 = _psi(p3, pAV, p4, p5);
- const double f3 = p3.mass()/m4J;
- const double f4 = p4.mass()/m4J;
- const double f5 = p5.mass()/m4J;
- const double fA = pA.mass()/m4J;
- const double fB = pB.mass()/m4J;
+ const double f3 = _safeMass(p3)/m4J;
+ const double f4 = _safeMass(p4)/m4J;
+ const double f5 = _safeMass(p5)/m4J;
+ const double fA = _safeMass(pA)/m4J;
+ const double fB = _safeMass(pB)/m4J;
const double XA = pA.E()/(pA.E()+pB.E());
const double psiAB = _psi(pA, pB, pA+pB, pAV);
@@ -251,7 +251,7 @@
}
if (sumEt < 420.0*GeV) return;
- const double m5J = jetsystem.mass();
+ const double m5J = _safeMass(jetsystem);
if (m5J < 750*GeV) return;
LorentzTransform cms_boost(-jetsystem.boostVector());
@@ -278,15 +278,15 @@
const double X3 = 2.0*p3.E()/m5J;
const double X4 = 2.0*p4.E()/m5J;
const double psi3 = _psi(p3, pAV, p4, p5);
- const double f3 = p3.mass()/m5J;
- const double f4 = p4.mass()/m5J;
- const double f5 = p5.mass()/m5J;
- const double fA = pA.mass()/m5J;
- const double fB = pB.mass()/m5J;
+ const double f3 = _safeMass(p3)/m5J;
+ const double f4 = _safeMass(p4)/m5J;
+ const double f5 = _safeMass(p5)/m5J;
+ const double fA = _safeMass(pA)/m5J;
+ const double fB = _safeMass(pB)/m5J;
const double XA = pA.E()/(pA.E()+pB.E());
const double psiAB = _psi(pA, pB, pA+pB, pAV);
- const double fC = pC.mass()/m5J;
- const double fD = pD.mass()/m5J;
+ const double fC = _safeMass(pC)/m5J;
+ const double fD = _safeMass(pD)/m5J;
const double XC = pC.E()/(pC.E()+pD.E());
const double psiCD = _psi(pC, pD, pC+pD, pAV);
@@ -383,14 +383,16 @@
}
FourMomentum _avg_beam_in_lab(const double& m, const double& y) {
- FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
- LorentzTransform cms_boost(-boostvec.boostVector());
- cms_boost = cms_boost.inverse();
const double mt = m/2.0;
FourMomentum beam1(mt, 0, 0, mt);
FourMomentum beam2(mt, 0, 0, -mt);
- beam1=cms_boost.transform(beam1);
- beam2=cms_boost.transform(beam2);
+ if (fabs(y)>1e-3) {
+ FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
+ LorentzTransform cms_boost(-boostvec.boostVector());
+ cms_boost = cms_boost.inverse();
+ beam1=cms_boost.transform(beam1);
+ beam2=cms_boost.transform(beam2);
+ }
if (beam1.E()>beam2.E()) {
return beam1-beam2;
}
@@ -405,6 +407,16 @@
Vector3 p3xp4 = p3.vector3().cross(p4.vector3());
return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit())));
}
+
+ double _safeMass(const FourMomentum& p) {
+ double mass2=p.mass2();
+ if (mass2>0.0) return sqrt(mass2);
+ else if (mass2<-1.0e-5) {
+ getLog() << Log::WARNING << "m2 = " << m2 << ". Assuming m2=0." << endl;
+ return 0.0;
+ }
+ else return 0.0;
+ }
private:
Modified: trunk/src/Analyses/CDF_1996_S3418421.cc
==============================================================================
--- trunk/src/Analyses/CDF_1996_S3418421.cc Tue Mar 2 22:58:09 2010 (r2301)
+++ trunk/src/Analyses/CDF_1996_S3418421.cc Wed Mar 3 10:24:55 2010 (r2302)
@@ -51,16 +51,16 @@
void analyze(const Event& event) {
const double weight = event.weight();
- Jets jets = applyProjection<FastJets>(event, "Jets").jetsByEt(50.0*GeV);
+ Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(50.0*GeV);
if (jets.size()<2) {
vetoEvent;
}
FourMomentum jet1 = jets[0].momentum();
FourMomentum jet2 = jets[1].momentum();
- double eta1 = fabs(jet1.eta());
- double eta2 = fabs(jet2.eta());
+ double eta1 = jet1.eta();
+ double eta2 = jet2.eta();
double chi = exp(fabs(eta1-eta2));
- if (eta2>2.0 || eta1>2.0 || chi>5.0) {
+ if (fabs(eta2)>2.0 || fabs(eta1)>2.0 || chi>5.0) {
vetoEvent;
}
Modified: trunk/src/Analyses/CDF_1997_S3541940.cc
==============================================================================
--- trunk/src/Analyses/CDF_1997_S3541940.cc Tue Mar 2 22:58:09 2010 (r2301)
+++ trunk/src/Analyses/CDF_1997_S3541940.cc Wed Mar 3 10:24:55 2010 (r2302)
@@ -81,7 +81,7 @@
vetoEvent;
}
- double m6J=jetsystem.mass();
+ double m6J=_safeMass(jetsystem);
if (m6J<520.0*GeV) {
vetoEvent;
}
@@ -126,27 +126,27 @@
_h_costheta3ppp->fill(costheta3ppp, weight);
double psi3ppp=_psi(p3ppp, pAV, p4ppp, p5ppp);
_h_psi3ppp->fill(psi3ppp, weight);
- _h_f3ppp->fill(p3ppp.mass()/m6J, weight);
- _h_f4ppp->fill(p4ppp.mass()/m6J, weight);
- _h_f5ppp->fill(p5ppp.mass()/m6J, weight);
+ _h_f3ppp->fill(_safeMass(p3ppp)/m6J, weight);
+ _h_f4ppp->fill(_safeMass(p4ppp)/m6J, weight);
+ _h_f5ppp->fill(_safeMass(p5ppp)/m6J, weight);
// 4 -> 3 jet variables
- _h_fApp->fill(pApp.mass()/m6J, weight);
- _h_fBpp->fill(pApp.mass()/m6J, weight);
+ _h_fApp->fill(_safeMass(pApp)/m6J, weight);
+ _h_fBpp->fill(_safeMass(pApp)/m6J, weight);
_h_XApp->fill(pApp.E()/(pApp.E()+pBpp.E()), weight);
double psiAppBpp=_psi(pApp, pBpp, pApp+pBpp, pAV);
_h_psiAppBpp->fill(psiAppBpp, weight);
// 5 -> 4 jet variables
- _h_fCp->fill(pCp.mass()/m6J, weight);
- _h_fDp->fill(pDp.mass()/m6J, weight);
+ _h_fCp->fill(_safeMass(pCp)/m6J, weight);
+ _h_fDp->fill(_safeMass(pDp)/m6J, weight);
_h_XCp->fill(pCp.E()/(pCp.E()+pDp.E()), weight);
double psiCpDp=_psi(pCp, pDp, pCp+pDp, pAV);
_h_psiCpDp->fill(psiCpDp, weight);
// 6 -> 5 jet variables
- _h_fE->fill(pE.mass()/m6J, weight);
- _h_fF->fill(pF.mass()/m6J, weight);
+ _h_fE->fill(_safeMass(pE)/m6J, weight);
+ _h_fF->fill(_safeMass(pF)/m6J, weight);
_h_XE->fill(pE.E()/(pE.E()+pF.E()), weight);
double psiEF=_psi(pE, pF, pE+pF, pAV);
_h_psiEF->fill(psiEF, weight);
@@ -207,14 +207,16 @@
}
FourMomentum _avg_beam_in_lab(const double& m, const double& y) {
- FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
- LorentzTransform cms_boost(-boostvec.boostVector());
- cms_boost = cms_boost.inverse();
const double mt = m/2.0;
FourMomentum beam1(mt, 0, 0, mt);
FourMomentum beam2(mt, 0, 0, -mt);
- beam1=cms_boost.transform(beam1);
- beam2=cms_boost.transform(beam2);
+ if (fabs(y)>1e-3) {
+ FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
+ LorentzTransform cms_boost(-boostvec.boostVector());
+ cms_boost = cms_boost.inverse();
+ beam1=cms_boost.transform(beam1);
+ beam2=cms_boost.transform(beam2);
+ }
if (beam1.E()>beam2.E()) {
return beam1-beam2;
}
@@ -230,6 +232,16 @@
return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit())));
}
+ double _safeMass(const FourMomentum& p) {
+ double mass2=p.mass2();
+ if (mass2>0.0) return sqrt(mass2);
+ else if (mass2<-1.0e-5) {
+ getLog() << Log::WARNING << "m2 = " << m2 << ". Assuming m2=0." << endl;
+ return 0.0;
+ }
+ else return 0.0;
+ }
+
private:
More information about the Rivet-svn
mailing list
|