Skip to content

Commit

Permalink
experimentell rewrite of inversion code, accomodate changes in TJet i…
Browse files Browse the repository at this point in the history 
…nterface
  • Loading branch information
stadie committed Feb 10, 2009
1 parent 0d473bc commit 8de290d
Show file tree
Hide file tree
Showing 5 changed files with 88 additions and 107 deletions.
77 changes: 39 additions & 38 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.9 2009/01/18 13:06:15 stadie Exp $
// $Id: Jet.cc,v 1.10 2009/01/22 15:30:30 stadie Exp $
//
#include "Jet.h"

Expand All @@ -12,43 +12,41 @@ Jet::Jet(double Et, double EmEt, double HadEt ,double OutEt, double E,
double (*func)(const TMeasurement *x, const double *par),
double (*errfunc)(const double *x, const TMeasurement *xorig, double err),
double* firstpar, int id, int npars)
: TJet(Et,EmEt,HadEt,OutEt,E,eta,phi,flavor), par(firstpar), npar(npars), parid(id),
f(func),errf(errfunc),varcoll(npars)
: TJet(Et,EmEt,HadEt,OutEt,E,eta,phi,flavor,0.0,1.0,1.0,1.0,1.0), par(firstpar),
npar(npars), parid(id), f(func),errf(errfunc),varcoll(npars)
{
temp = *this;
}

//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
int Jet::varyPar(int i, double eps, double Et, double scale, double &upperEt, double& lowerEt)
Jet::Jet(double Et, double EmEt, double HadEt ,double OutEt, double E,
double eta,double phi, Flavor flavor, double genPt, double ZSPcor,
double JPTcor, double L2cor, double L3cor,
double (*func)(const TMeasurement *x, const double *par),
double (*errfunc)(const double *x, const TMeasurement *xorig, double err),
double* firstpar, int id, int npars)
: TJet(Et,EmEt,HadEt,OutEt,E,eta,phi,flavor,genPt,ZSPcor,JPTcor,L2cor,L3cor),
par(firstpar), npar(npars), parid(id), f(func),errf(errfunc),varcoll(npars)
{
double orig = par[i];
par[i] += eps;
upperEt = expectedEt(Et,scale,true);
par[i] = orig - eps;
lowerEt = expectedEt(Et,scale,true);
par[i] = orig;
return parid + i;
temp = *this;
}

// varies all parameters for this jet by eps and returns a vector of the
// parameter id and the Et for the par + eps and par - eps variation
const Jet::VariationColl& Jet::varyPars(double eps, double Et, double scale)
const Jet::VariationColl& Jet::varyPars(double eps, double Et, double start)
{
//scale = Et;
double s = scale;
//start = Et;
for(int i = 0 ; i < npar ; ++i) {
double orig = par[i];
par[i] += eps;
varcoll[i].upperEt = expectedEt(Et,s,true);
varcoll[i].upperEt = expectedEt(Et,start);
if( varcoll[i].upperEt < 0) varcoll[i].upperEt = 0.999 * pt;
varcoll[i].upperError = expectedError(varcoll[i].upperEt);
//varcoll[i].upperEt = expectedEt(Et,s,false);
s = scale;
par[i] = orig - eps;;
varcoll[i].lowerEt = expectedEt(Et,s,true);
varcoll[i].lowerEt = expectedEt(Et,start);
if( varcoll[i].lowerEt < 0) varcoll[i].lowerEt = 1.001 * pt;
varcoll[i].lowerError = expectedError(varcoll[i].lowerEt);
//varcoll[i].lowerEt = expectedEt(Et,s,false);
s = scale;
par[i] = orig;
varcoll[i].parid = parid + i;
}
Expand Down Expand Up @@ -83,35 +81,37 @@ double Jet::correctedEt(double Et, bool fast) const {
return corEt;
}

double Jet::expectedEt(double truth, double& scale, bool extrapolate)
double Jet::expectedEt(double truth, double start, bool fast)
{
static const double eps = 1.0e-5;
const double up = 4 * truth;
const double low = 0.2 * truth;
static const double eps = 1.0e-8;
//const double up = 4 * truth;
//const double low = 0.2 * truth;
double x1 = start,x2;
//find root of truth - jet->correctedEt(expectedEt)
// x: expectedEt
// y: truth - jet->correctedEt(expectedEt)
// f: jet->correctedEt(expectedEt)
double x1 = scale;
double f1 = correctedEt(x1,extrapolate);
double f1 = correctedEt(x1,fast);
//get second point assuming a constant correction factor
double x2 = (truth - EMF - OutF) * (x1 - EMF - OutF)/(f1 - EMF - OutF) + EMF;
if((x2 > up )||(x2 < low)) x2 = scale;
double f2 = correctedEt(x2,true);
double y2 = truth - f2;
x2 = (truth - EMF - OutF) * (x1 - EMF - OutF)/(f1 - EMF - OutF) + EMF;
//if((x2 > up )||(x2 < low)) x2 = x1;
///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;
}
if(! secant(truth,x1,x2,eps)) return -1;
scale = x2;
f1 = correctedEt(scale,true);
x2 = (truth - EMF - OutF) * (scale - EMF - OutF)/(f1 - EMF - OutF) + EMF + OutF;
*/
if(! secant(truth,x2,x1,eps)) return -1;
//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;
//std::cout << "x1=" << x1 << " x2=" << x2 << '\n';
assert(std::abs(correctedEt(x2)-truth)/truth < eps);
return x2;
}

bool Jet::falseposition(double truth, double& x1, double& x2,double eps)
Expand Down Expand Up @@ -181,7 +181,7 @@ bool Jet::falseposition(double truth, double& x1, double& x2,double eps)
}


bool Jet::secant(double truth, double& x1, double& x2,double eps)
bool Jet::secant(double truth, double& x2, double& x1,double eps)
{
//x2 is the best estimate!
const double up = 4 * truth;
Expand All @@ -197,6 +197,7 @@ bool Jet::secant(double truth, double& x1, double& x2,double eps)
x2 = 1.0001 * x1;
dx = std::abs(x1-x2);
}
//std::cout << "first intervall size:" << dx/x1 << '\n';
while((dx/x1 > eps)&&(i < 100)) {
//std::cout << i << ":" << x1 << ", " << x2 << " : " << y1 << ", " << y2 << std::endl;
double x3 = x1 + y1 * (x2-x1)/(f2 - f1);
Expand Down Expand Up @@ -247,7 +248,7 @@ bool Jet::secant(double truth, double& x1, double& x2,double eps)
// }
}
ntries += i;
if(std::abs(y2) > 0.001 * truth) {
if(std::abs(y2) > eps * truth) {
//std::cout << "failed to find good root\n";
//std::cout << i << ":" << x1 << ", " << x2 << ":" << truth - f2 << "\n";
++nfails;
Expand Down
15 changes: 9 additions & 6 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.6 2009/01/16 08:46:40 stadie Exp $
// $Id: Jet.h,v 1.7 2009/01/22 15:30:30 stadie Exp $
//
#ifndef JET_H
#define JET_H
Expand All @@ -17,6 +17,12 @@ class Jet : public TJet
double (*func)(const TMeasurement *x, const double *par),
double (*errfunc)(const double *x, const TMeasurement *xorig, double err),
double* firstpar, int id, int npars);
Jet(double Et, double EmEt, double HadEt ,double OutEt, double E,
double eta,double phi, Flavor flavor, double genPt, double ZSPcor,
double JPTcor, double L2cor, double L3cor,
double (*func)(const TMeasurement *x, const double *par),
double (*errfunc)(const double *x, const TMeasurement *xorig, double err),
double* firstpar, int id, int npars);
virtual ~Jet() {};
double Et() const {return pt;}
double EmEt() const {return EMF;}
Expand All @@ -28,13 +34,10 @@ class Jet : public TJet
Flavor flavor() const {return TJet::flavor;}
virtual void ChangeParAddress(double* oldpar, double* newpar) {par += newpar - oldpar;}
virtual double correctedEt(double Et, bool fast = false) const;
double expectedEt(double truth, double& scale, bool extrapolate = false);
double expectedEt(double truth, double start, bool fast = false);
virtual double Error() const {return errf(&(TMeasurement::pt),this,0);}
virtual double expectedError(double truth) const { return errf(&truth,this,0);}
virtual int nPar() const {return npar;}
//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);
struct ParameterVariation {
int parid;
double upperEt;
Expand All @@ -46,7 +49,7 @@ class Jet : public TJet
typedef std::vector<ParameterVariation>::const_iterator VariationCollIter;
// varies all parameters for this jet by eps and returns a vector of the
// parameter id and the Et for the par + eps and par - eps variation
virtual const VariationColl& varyPars(double eps, double Et, double scale);
virtual const VariationColl& varyPars(double eps, double Et, double start);
virtual const VariationColl& varyParsDirectly(double eps);

static void printInversionStats();
Expand Down
56 changes: 31 additions & 25 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.7 2009/01/22 15:30:30 stadie Exp $
// $Id: JetTruthEvent.cc,v 1.8 2009/01/22 17:48:10 stadie Exp $
//

#include "JetTruthEvent.h"
Expand Down Expand Up @@ -51,33 +51,42 @@ double JetTruthEvent::chi2_fast_simple_scaled(double * temp_derivative1,
double const epsilon) const
{
double et = jet->correctedEt(jet->Et());
double c = (et - jet->EmEt() - jet->OutEt()) / jet->HadEt();
double c = (et - jet->EmEt() - jet->OutEt()) / jet->HadEt();
if(c == 0) c = 1.0;
double err2inv = c * jet->expectedError((truth - jet->EmEt() - jet->OutEt())/c + jet->EmEt() + jet->OutEt() );
err2inv *= err2inv;
err2inv = 1/err2inv;
double chi2 = truth - et;
chi2 *= chi2 * err2inv;
chi2 = weight * TData::ScaleResidual(chi2);
if(chi2 != chi2) {//check for NAN
std::cout << truth << ", " << et << ", " << jet->Et() << ", " << c << ", " << chi2 << '\n';
}
double temp1,temp2;
const Jet::VariationColl& varcoll = jet->varyParsDirectly(epsilon);
for(Jet::VariationCollIter i = varcoll.begin() ; i != varcoll.end() ; ++i) {
temp1 = truth - i->lowerEt;
c = (i->lowerEt - jet->EmEt() - jet->OutEt()) / jet->HadEt();
c = (i->lowerEt - jet->EmEt() - jet->OutEt()) / jet->HadEt();
if(c == 0) c = 1.0;
err2inv = c*jet->expectedError((truth - jet->EmEt() - jet->OutEt())/c + jet->EmEt() + jet->OutEt() );
err2inv *= err2inv;
err2inv = 1/err2inv;
temp1 *= temp1 * err2inv;
temp1 = weight * TData::ScaleResidual(temp1);
assert(temp1 == temp1);
temp2 = truth - i->upperEt;
c = (i->upperEt - jet->EmEt() - jet->OutEt()) / jet->HadEt();
c = (i->upperEt - jet->EmEt() - jet->OutEt()) / jet->HadEt();
if(c == 0) c = 1.0;
err2inv = c * jet->expectedError((truth - jet->EmEt() - jet->OutEt())/c + jet->EmEt() + jet->OutEt() );
err2inv *= err2inv;
err2inv = 1/err2inv;
temp2 *= temp2 * err2inv;
temp2 = weight * TData::ScaleResidual(temp2);
assert(temp2 == temp2);
temp_derivative1[i->parid] += (temp2 - temp1); // for 1st derivative
temp_derivative2[i->parid] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative
temp_derivative2[i->parid] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative
}
assert(chi2 == chi2);
return chi2;
}

Expand Down Expand Up @@ -121,59 +130,56 @@ double JetTruthEvent::chi2_fast_invert(double * temp_derivative1,
double const epsilon) const
{
//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;
double expectedEt = jet->expectedEt(truth,truth);
if(expectedEt < 0) {
return chi2_fast_simple_scaled(temp_derivative1,temp_derivative2,epsilon);
}
//calculate chi2
double err2inv = jet->expectedError(expectedEt);
//std::cout << "Jet Et:" << expectedEt << "," << jet->Et() << " error:" << err2inv << '\n';
//std::cout << "Jet Et:" << expectedEt << "," << jet->Et() << " error:" << err2inv
// << " true Et:" << truth << '\n';
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';
std::cout << truth << ", " << expectedEt << ", " << chi2 << '\n';
}

//calculate chi2 for derivatives
double temp1,temp2;
const Jet::VariationColl& varcoll = jet->varyPars(epsilon,truth,scale);
const Jet::VariationColl& varcoll = jet->varyPars(epsilon,truth,expectedEt);
for(Jet::VariationCollIter i = varcoll.begin() ; i != varcoll.end() ; ++i) {
if((std::abs((i->lowerEt - expectedEt)/scale) > 0.01) ||
(std::abs((i->upperEt - expectedEt)/scale) > 0.01)) {
if((std::abs((i->lowerEt - expectedEt)/expectedEt) > 0.01) ||
(std::abs((i->upperEt - expectedEt)/expectedEt) > 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();
err2inv = i->lowerError;
err2inv *= err2inv;
err2inv = 1/err2inv;
temp1 *= temp1 * err2inv;
temp1 = weight * TData::ScaleResidual(temp1);
assert(temp1 == temp1);
temp2 = i->upperEt - jet->Et();
err2inv = i->upperError;
err2inv *= err2inv;
err2inv = 1/err2inv;
temp2 *= temp2 * err2inv;
temp2 = weight * TData::ScaleResidual(temp2);
temp2 = weight * TData::ScaleResidual(temp2);
assert(temp2 == temp2);
temp_derivative1[i->parid] += (temp2 - temp1); // for 1st derivative
temp_derivative2[i->parid] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative
assert( temp_derivative1[i->parid] == temp_derivative1[i->parid] );
assert( temp_derivative2[i->parid] == temp_derivative2[i->parid] );
}
assert(chi2 == chi2);
//std::cout << "Jet Et:" << jet->Et() << " expected Et:" << expectedEt << " error:" << sqrt(1/err2inv)
// << " true Et:" << truth << " chi2:" << chi2 << '\n';
return chi2;
}
40 changes: 7 additions & 33 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.5 2009/01/16 08:46:40 stadie Exp $
// $Id: JetWithTowers.cc,v 1.6 2009/01/22 15:30:30 stadie Exp $
//
#include"JetWithTowers.h"

Expand Down Expand Up @@ -61,37 +61,11 @@ double JetWithTowers::correctedEt(double Et,bool fast) const
return Jet::correctedEt(ccet);
}


//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
int JetWithTowers::varyPar(int i, double eps, double Et, double scale,
double& upperEt, double& lowerEt)
{
if(i < njetpars) return Jet::varyPar(i,eps,Et,scale,upperEt,lowerEt);
int towid = (i - njetpars)/ntowerpars;
int towpar = (i - njetpars) % ntowerpars;
//std::cout << "I: " << i << " tow:" << towid << "; " << towpar << std::endl;
std::map<int,double*>::const_iterator iter = towerpars.begin();
for(int i = 0 ; i < towid ; ++i) ++iter;
double *p = iter->second;
int id = iter->first;
//std::cout << "truth: " << Et << "alternating par:" << id + towpar << " = " << p[towpar] << std::endl;
double orig = p[towpar];
p[towpar] += eps;
double s = scale;
upperEt = expectedEt(Et,s,true);
p[towpar] = orig - eps;
s = scale;
lowerEt = expectedEt(Et,s,true);
p[towpar] = orig;
return id + towpar;
}

// varies all parameters for this jet by eps and returns a vector of the
// parameter id and the Et for the par + eps and par - eps variation
const Jet::VariationColl& JetWithTowers::varyPars(double eps, double Et, double scale)
const Jet::VariationColl& JetWithTowers::varyPars(double eps, double Et, double start)
{
Jet::varyPars(eps,Et,scale);
Jet::varyPars(eps,Et,start);
int i = njetpars;

for(std::map<int,double*>::const_iterator iter = towerpars.begin() ;
Expand All @@ -104,12 +78,12 @@ const Jet::VariationColl& JetWithTowers::varyPars(double eps, double Et, double
// << std::endl;
double orig = p[towpar];
p[towpar] += eps;
double s = scale;
varcoll[i].upperEt = expectedEt(Et,s,true);
varcoll[i].upperEt = expectedEt(Et,start);
if( varcoll[i].upperEt < 0) varcoll[i].upperEt = pt;
varcoll[i].upperError = expectedError(varcoll[i].upperEt);
p[towpar] = orig - eps;
s = scale;
varcoll[i].lowerEt = expectedEt(Et,s,true);
varcoll[i].lowerEt = expectedEt(Et,start);
if( varcoll[i].lowerEt < 0) varcoll[i].lowerEt = pt;
varcoll[i].lowerError = expectedError(varcoll[i].lowerEt);
p[towpar] = orig;
varcoll[i].parid = id + towpar;
Expand Down
Loading

0 comments on commit 8de290d

Please sign in to comment.