Skip to content

Commit

Permalink
compute fraction of tower to jet Et based on the hadronic Et, use old…
Browse files Browse the repository at this point in the history 
… chi2 computation for events that fail the inversion
  • Loading branch information
stadie committed Jan 13, 2009
1 parent 3972333 commit 5a74398
Show file tree
Hide file tree
Showing 5 changed files with 108 additions and 68 deletions.
77 changes: 45 additions & 32 deletions Jet.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Class for basic jets
//
// first version: Hartmut Stadie 2008/12/14
// $Id: Jet.cc,v 1.5 2009/01/06 13:42:36 stadie Exp $
// $Id: Jet.cc,v 1.6 2009/01/09 18:09:58 stadie Exp $
//
#include "Jet.h"

Expand Down Expand Up @@ -53,7 +53,7 @@ const Jet::VariationColl& Jet::varyPars(double eps, double Et, double scale)
}


double Jet::correctedEt(double Et) const {
double Jet::correctedEt(double Et, bool fast) const {
//assume that only the hadronic energy gets modified!
temp.pt = Et;
temp.HadF = Et - OutF - EMF;
Expand All @@ -71,30 +71,32 @@ double Jet::expectedEt(double truth, double& scale, bool extrapolate)
// y: truth - jet->correctedEt(expectedEt)
// f: jet->correctedEt(expectedEt)
double x1 = scale;
double f1 = correctedEt(x1);
double f1 = correctedEt(x1,extrapolate);
//get second point assuming a constant correction factor
double x2 = (truth - EMF) * (x1 - EMF)/(f1 - EMF) + EMF;
double x2 = (truth - EMF - OutF) * (x1 - EMF - OutF)/(f1 - EMF - OutF) + EMF;
if((x2 > up )||(x2 < low)) x2 = scale;
double f2 = correctedEt(x2);
double f2 = correctedEt(x2,true);
double y2 = truth - f2;
//std::cout << "truth:" << truth << " scale:" << scale << " f1:" << f1 << " x2:" << x2
// << " f2:" << f2 << '\n';
if(extrapolate || (std::abs(y2) < eps)) {
//std::cout << "extrapolated:" << x2 << ", " << y2 << " at scale " << scale << std::endl;
return x2;
}
x2 = secant(truth,x1,x2,eps);
//x2 = falseposition(truth,x1,x2,eps);
if(! secant(truth,x1,x2,eps)) return -1;
scale = x2;
f1 = correctedEt(scale);
x2 = (truth - EMF) * (scale - EMF)/(f1 - EMF) + EMF;
f1 = correctedEt(scale,true);
x2 = (truth - EMF - OutF) * (scale - EMF - OutF)/(f1 - EMF - OutF) + EMF + OutF;
//std::cout << i << ": scale:" << scale << ", expected:" << (truth - EMF) * (scale - EMF)/(f1 - EMF) + EMF << " dist for scale:" << truth - f1 << "\n";
//assert(std::abs(correctedEt(x2)-truth)/truth < eps);
return ((x2 < up )&&(x2 > low)) ? x2 : scale;
}

double Jet::falseposition(double truth, double x1, double x2,double eps)
bool Jet::falseposition(double truth, double& x1, double& x2,double eps)
{
double f1 = correctedEt(x1);
double f2 = correctedEt(x2);
//x2 is the best estimate!
double f1 = correctedEt(x1,true);
double f2 = correctedEt(x2,true);
double temp;
double step = 0.1 * truth;
++ncalls;
Expand All @@ -116,13 +118,13 @@ double Jet::falseposition(double truth, double x1, double x2,double eps)
// << " y1,2:" << y1 << ", " << y2 << std::endl;
if(f1 > truth) {
x1 -= step;
f1 = correctedEt(x1);
f1 = correctedEt(x1,true);
y1 = truth - f1;
++ntries;
}
if(f2 < truth) {
x2 += step;
f2 = correctedEt(x2);
f2 = correctedEt(x2,true);
y2 = truth - f2;
++ntries;
}
Expand All @@ -133,7 +135,7 @@ double Jet::falseposition(double truth, double x1, double x2,double eps)
while(std::abs((x2-x1)/x1) > eps) {
//std::cout << i << ":" << x1 << ", " << x2 << " : " << y1 << ", " << y2 << std::endl;
double x3 = x1 + y1 * (x2-x1)/(f2 - f1);
double f3 = correctedEt(x3);
double f3 = correctedEt(x3,true);
double y3 = truth - f3;
++i;
if(y1 * y3 < 0) {
Expand All @@ -147,25 +149,29 @@ double Jet::falseposition(double truth, double x1, double x2,double eps)
}
if(i > 100) {
++nfails;
break;
ntries += i;
return false;
}
}
ntries += i;
return 0.5*(x1 + x2);
x2 = 0.5*(x1 + x2);
return true;
}

double Jet::secant(double truth, double x1, double x2,double eps)

bool Jet::secant(double truth, double& x1, double& x2,double eps)
{
//x2 is the best estimate!
const double up = 4 * truth;
const double low = 0.2 * truth;
double f1 = correctedEt(x1);
double f2 = correctedEt(x2);
double f1 = correctedEt(x1,true);
double f2 = correctedEt(x2,true);
double y2 = truth - f2;
double y1 = truth - f1;
++ncalls;
int i = 0;
double dx = std::abs(x1-x2), dx1 = truth, dx2 = truth;
while(std::abs(dx/x1) > eps) {
while((dx/x1 > eps)&&(i < 100)) {
//std::cout << i << ":" << x1 << ", " << x2 << " : " << y1 << ", " << y2 << std::endl;
double x3 = x1 + y1 * (x2-x1)/(f2 - f1);
if((x3 > up )||(x3 < low)) {
Expand All @@ -174,7 +180,7 @@ double Jet::secant(double truth, double x1, double x2,double eps)
dx2 = dx1;
dx1 = dx;
dx = std::abs(x2 - x3);
if((i > 0) && (dx2 < dx)) {
if(dx2 < dx) {
//std::cout << "Warning: fit alternating!\n";
//std::cout << i << ": last three intervall sizes " << dx << ", "
// << dx1 << ", " << dx2 << std::endl;
Expand All @@ -188,7 +194,7 @@ double Jet::secant(double truth, double x1, double x2,double eps)
++i;
continue;
}
double f3 = correctedEt(x3);
double f3 = correctedEt(x3,true);
double y3 = truth - f3;
//use false position if root is bracketed
if(y1 * y3 < 0) {
Expand All @@ -205,16 +211,23 @@ double Jet::secant(double truth, double x1, double x2,double eps)
}
++i;
//std::cout << i << ":" << x1 << ", " << x2 << ":" << truth - f1 << "; " << truth - f2 << "\n";
if(i > 2000) {
//std::cout << "failed to find good root\n";
//std::cout << i << ":" << x1 << ", " << x2 << ":" << truth - f2 << "\n";
x2 = 0.5 * (x2+x1);
++nfails;
break;
}
}
// if(i > 100) {
// //std::cout << "failed to find good root\n";
// //std::cout << i << ":" << x1 << ", " << x2 << ":" << truth - f2 << "\n";
// x2 = 0.5 * (x2+x1);
// ++nfails;
// ntries += i;
// return false;
// }
}
ntries += i;
return x2;
if(std::abs(y2) > 0.001 * truth) {
//std::cout << "failed to find good root\n";
//std::cout << i << ":" << x1 << ", " << x2 << ":" << truth - f2 << "\n";
++nfails;
return false;
}
return true;
}

int Jet::ncalls = 0;
Expand Down
8 changes: 4 additions & 4 deletions Jet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Class for basic jets
//
// first version: Hartmut Stadie 2008/12/14
// $Id: Jet.h,v 1.3 2009/01/04 16:21:06 stadie Exp $
// $Id: Jet.h,v 1.4 2009/01/09 18:09:58 stadie Exp $
//
#ifndef JET_H
#define JET_H
Expand All @@ -26,7 +26,7 @@ class Jet : public TJet
double phi() const {return TMeasurement::phi;}
Flavor flavor() const {return TJet::flavor;}
virtual void ChangeParAddress(double* oldpar, double* newpar) {par += newpar - oldpar;}
virtual double correctedEt(double Et) const;
virtual double correctedEt(double Et, bool fast = false) const;
double expectedEt(double truth, double& scale, bool extrapolate = false);
double Error() const {return error;}
virtual int nPar() const {return npar;}
Expand Down Expand Up @@ -54,8 +54,8 @@ class Jet : public TJet
mutable VariationColl varcoll;
private:
mutable TMeasurement temp;
double secant(double truth, double x1, double x2, double eps);
double falseposition(double truth, double x1, double x2, double eps);
bool secant(double truth, double& x1, double& x2, double eps);
bool falseposition(double truth, double& x1, double& x2, double eps);
static int ncalls, ntries, nfails, nwarns;
};

Expand Down
30 changes: 23 additions & 7 deletions JetTruthEvent.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Class for all events with one jet and truth informatio
//
// first version: Hartmut Stadie 2008/12/14
// $Id: JetTruthEvent.cc,v 1.3 2008/12/27 16:39:02 stadie Exp $
// $Id: JetTruthEvent.cc,v 1.4 2009/01/04 16:21:06 stadie Exp $
//

#include "JetTruthEvent.h"
Expand All @@ -21,22 +21,38 @@ double JetTruthEvent::chi2_fast(double * temp_derivative1, double * temp_derivat
//find expected measurement of jet Et
double scale = truth;
double expectedEt = jet->expectedEt(truth,scale);
if(expectedEt < 0) {//return 1.5;
//return chi2 from modified measurement instead!
double et = jet->correctedEt(jet->Et());
double s= et/jet->Et();
double err2inv = s * jet->Error();
err2inv *= err2inv;
err2inv = 1/err2inv;
double chi2 = truth - et;
chi2 *= chi2 * err2inv;
chi2 = weight * TData::ScaleResidual(chi2);
return chi2;
}
//calculate chi2
double err2inv = jet->Error();
err2inv *= err2inv;
err2inv = 1/err2inv;
double chi2 = jet->Et() - expectedEt;
chi2 *= chi2 * err2inv;
chi2 = weight * TData::ScaleResidual(chi2);
//if(chi2 != chi2) {//check for NAN
// std::cout << truth << ", " << scale << ", " << expectedEt << ", " << chi2 << '\n';
//}
//calculate chi2 for derivatives
double temp1,temp2;
const Jet::VariationColl& varcoll = jet->varyPars(epsilon,truth,scale);
for(Jet::VariationCollIter i = varcoll.begin() ; i != varcoll.end() ; ++i) {
if((std::abs((i->lowerEt - expectedEt)/scale) > 0.05) ||
(std::abs((i->upperEt - expectedEt)/scale) > 0.05)) {
std::cout << "strange extrapolation result:" << expectedEt << "; "
<< i->lowerEt << "; " << i->upperEt << " jet:" << jet->Et()
<< " truth:" << truth << std::endl;
if((std::abs((i->lowerEt - expectedEt)/scale) > 0.01) ||
(std::abs((i->upperEt - expectedEt)/scale) > 0.01)) {
std::cout << "strange extrapolation result modifying par:" << i->parid << ":"
<< expectedEt << " - " << i->lowerEt << " + " << i->upperEt
<< " uncor jet Et:" << jet->Et() << " truth:" << truth << std::endl;
continue;
}
temp1 = i->lowerEt - jet->Et();
temp1 *= temp1 * err2inv;
Expand All @@ -46,6 +62,6 @@ double JetTruthEvent::chi2_fast(double * temp_derivative1, double * temp_derivat
temp2 = weight * TData::ScaleResidual(temp2);
temp_derivative1[i->parid] += (temp2 - temp1); // for 1st derivative
temp_derivative2[i->parid] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative
}
}
return chi2;
}
48 changes: 28 additions & 20 deletions JetWithTowers.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Class for jets with towers
//
// first version: Hartmut Stadie 2008/12/25
// $Id: JetWithTowers.cc,v 1.2 2009/01/04 16:21:06 stadie Exp $
// $Id: JetWithTowers.cc,v 1.3 2009/01/09 18:09:58 stadie Exp $
//
#include"JetWithTowers.h"

Expand Down Expand Up @@ -37,17 +37,24 @@ void JetWithTowers::ChangeParAddress(double* oldpar, double* newpar)
}
}

double JetWithTowers::correctedEt(double Et) const
double JetWithTowers::correctedEt(double Et,bool fast) const
{
double cet = 0;
for(TowerCollConstIter i = towers.begin() ; i != towers.end() ; ++i) {
cet += (*i)->projectionToJetAxis() * (*i)->correctedEt((*i)->Et());
}
//std::cout << "jet ET:" << pt << " sum of tower:" << cet << "\n";
double ccet = 0;
for(TowerCollConstIter i = towers.begin() ; i != towers.end() ; ++i) {
ccet += (*i)->projectionToJetAxis() *
(*i)->correctedEt((*i)->lastCorrectedEt()/cet * Et);
double ccet = EmEt() + OutEt();
double HadEt = Et - ccet;
if(! fast) {
double chad = 0;
for(TowerCollConstIter i = towers.begin() ; i != towers.end() ; ++i) {
chad += (*i)->projectionToJetAxis() * (*i)->correctedHadEt((*i)->HadEt());
}
//std::cout << "jet ET:" << pt << " sum of tower:" << chad + EMF + OutF << "\n";
for(TowerCollConstIter i = towers.begin() ; i != towers.end() ; ++i) {
(*i)->setFractionOfJetHadEt((*i)->lastCorrectedHadEt()/chad);
ccet += (*i)->projectionToJetAxis() * (*i)->correctedHadEt((*i)->fractionOfJetHadEt() * HadEt);
}
} else {
for(TowerCollConstIter i = towers.begin() ; i != towers.end() ; ++i) {
ccet += (*i)->projectionToJetAxis() * (*i)->correctedHadEt((*i)->fractionOfJetHadEt() * HadEt);
}
}
//std::cout << "scale ET:" << Et << " cor. sum of tower:" << ccet << "\n";
return Jet::correctedEt(ccet);
Expand All @@ -67,7 +74,7 @@ int JetWithTowers::varyPar(int i, double eps, double Et, double scale,
for(int i = 0 ; i < towid ; ++i) ++iter;
double *p = iter->second;
int id = iter->first;
//std::cout << "alternating par:" << id + towpar << " = " << p[towpar] << std::endl;
//std::cout << "truth: " << Et << "alternating par:" << id + towpar << " = " << p[towpar] << std::endl;
double orig = p[towpar];
p[towpar] += eps;
double s = scale;
Expand All @@ -92,7 +99,8 @@ const Jet::VariationColl& JetWithTowers::varyPars(double eps, double Et, double
int id = iter->first;
//std::cout << p << "," << id << " tower[0]:" << towers[0]->Par() << '\n';
for(int towpar = 0 ; towpar < ntowerpars ; ++towpar) {
//std::cout << "alternating par:" << id + towpar << " = " << p[towpar] << std::endl;
//std::cout << "truth: " << Et << " alternating par:" << id + towpar << " = " << p[towpar]
// << std::endl;
double orig = p[towpar];
p[towpar] += eps;
double s = scale;
Expand Down Expand Up @@ -139,19 +147,19 @@ JetWithTowers::Tower::Tower(double Et, double EmEt, double HadEt ,
double alpha,double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars)
: TMeasurement(Et,EmEt,HadEt,OutEt,E,eta,phi), alpha(alpha), par(firstpar),
npar(npars), parid(id), error(err),f(func)
npar(npars), parid(id), error(err), lastCorHadEt(0), fraction(0), f(func)
{
temp = *this;
}

double JetWithTowers::Tower::correctedEt(double Et) const
double JetWithTowers::Tower::correctedHadEt(double HadEt) const
{
//assume that only the hadronic energy gets modified!
temp.pt = Et;
temp.HadF = Et - OutF - EMF;
temp.E = TMeasurement::E * Et/pt;
lastCorEt = f(&temp,par);
temp.pt = HadEt + OutF + EMF;
temp.HadF = HadEt;
temp.E = TMeasurement::E * temp.pt/pt;
lastCorHadEt = f(&temp,par) - OutF - EMF;
//std::cout << pt << ", " << Et << ":" << lastCorEt << " par:" << par[0] << std::endl;
return lastCorEt;
return lastCorHadEt;
}

13 changes: 8 additions & 5 deletions JetWithTowers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Class for jets with towers
//
// first version: Hartmut Stadie 2008/12/25
// $Id: JetWithTowers.h,v 1.2 2009/01/04 16:21:06 stadie Exp $
// $Id: JetWithTowers.h,v 1.3 2009/01/09 18:09:58 stadie Exp $
//
#ifndef JETWITHTOWERS_H
#define JETWITHTOWERS_H
Expand All @@ -22,7 +22,7 @@ class JetWithTowers : public Jet
virtual ~JetWithTowers();
virtual int nPar() const {return njetpars + towerpars.size() * ntowerpars;}
virtual void ChangeParAddress(double* oldpar, double* newpar);
virtual double correctedEt(double Et) const;
virtual double correctedEt(double Et,bool fast = false) const;
//varies the i'th parameter for this jet by eps and returns its overall
// parameter id and sets the Et for the par + eps and par - eps result
virtual int varyPar(int i, double eps, double Et, double scale, double& upperEt, double& lowerEt);
Expand All @@ -48,9 +48,11 @@ class JetWithTowers : public Jet
double eta() const {return TMeasurement::eta;}
double phi() const {return TMeasurement::phi;}
double projectionToJetAxis() const {return alpha;}
double fractionOfJetHadEt() const { return fraction;}
void setFractionOfJetHadEt(double frac) const { fraction = frac;}
void ChangeParAddress(double* oldpar, double* newpar) {par += newpar - oldpar;}
double correctedEt(double Et) const;
double lastCorrectedEt() const { return lastCorEt;}
double correctedHadEt(double HadEt) const;
double lastCorrectedHadEt() const { return lastCorHadEt;}
double Error() const {return error;}
int nPar() const {return npar;}
int FirstPar() const {return parid;}
Expand All @@ -61,7 +63,8 @@ class JetWithTowers : public Jet
int npar,parid;
double error;
mutable TMeasurement temp;
mutable double lastCorEt;
mutable double lastCorHadEt;
mutable double fraction;
double const(*f)(TMeasurement *const x, double *const par);
};
typedef std::vector<Tower*> TowerColl;
Expand Down

0 comments on commit 5a74398

Please sign in to comment.