diff --git a/Jet.cc b/Jet.cc index db0bc88..97ac24b 100644 --- a/Jet.cc +++ b/Jet.cc @@ -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" @@ -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; } @@ -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) @@ -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; @@ -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); @@ -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; diff --git a/Jet.h b/Jet.h index 680da42..33452ab 100644 --- a/Jet.h +++ b/Jet.h @@ -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 @@ -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;} @@ -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; @@ -46,7 +49,7 @@ class Jet : public TJet typedef std::vector::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(); diff --git a/JetTruthEvent.cc b/JetTruthEvent.cc index 50957e5..07a17ea 100644 --- a/JetTruthEvent.cc +++ b/JetTruthEvent.cc @@ -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" @@ -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; } @@ -121,23 +130,14 @@ 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; @@ -145,35 +145,41 @@ double JetTruthEvent::chi2_fast_invert(double * temp_derivative1, 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; } diff --git a/JetWithTowers.cc b/JetWithTowers.cc index 1318a6f..d88a618 100644 --- a/JetWithTowers.cc +++ b/JetWithTowers.cc @@ -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" @@ -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::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::const_iterator iter = towerpars.begin() ; @@ -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; diff --git a/JetWithTowers.h b/JetWithTowers.h index 9196ab7..8c17d4a 100644 --- a/JetWithTowers.h +++ b/JetWithTowers.h @@ -2,7 +2,7 @@ // Class for jets with towers // // first version: Hartmut Stadie 2008/12/25 -// $Id: JetWithTowers.h,v 1.5 2009/01/16 08:46:40 stadie Exp $ +// $Id: JetWithTowers.h,v 1.6 2009/01/22 15:30:30 stadie Exp $ // #ifndef JETWITHTOWERS_H #define JETWITHTOWERS_H @@ -26,12 +26,9 @@ class JetWithTowers : public Jet virtual double correctedEt(double Et,bool fast = false) const; virtual double Error() const; virtual double expectedError(double truth) 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); // 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); void addTower(double Et, double EmEt, double HadEt ,double OutEt, double E,