Skip to content

Commit

Permalink
first version of class for jet withs towers.
Browse files Browse the repository at this point in the history
  • Loading branch information
stadie committed Dec 27, 2008
1 parent 85fa378 commit 9ca80f1
Show file tree
Hide file tree
Showing 5 changed files with 255 additions and 71 deletions.
74 changes: 17 additions & 57 deletions JetTruthEvent.cc
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.1 2008/12/16 15:21:27 stadie Exp $
// $Id: JetTruthEvent.cc,v 1.2 2008/12/17 09:38:36 stadie Exp $
//

#include "JetTruthEvent.h"
Expand All @@ -19,72 +19,32 @@ double JetTruthEvent::chi2_fast(double * temp_derivative1, double * temp_derivat
double const epsilon) const
{
//find expected measurement of jet Et
double expectedEt = 0;
{
//find root of truth - jet->correctedEt(expectedEt)
int i = 0;
double x1 = truth * 0.5,x2 = 2 * truth;
double y1 = truth - jet->correctedEt(x1), y2 = truth - jet->correctedEt(x2);
//bag the root
while(((y2 > 0)||(y1 < 0))&& i < 5) {
if(y2 > 0) x2 *= 2;
else x1 /= 2;
++i;
}
double xm,ym;
if(i == 5) {
i = 99;
x1 = x2 = truth;
}
else i = 0;
while(std::abs(x1 - x2) > 0.1) {
if(i < 20) break;
xm = 0.5 * (x1 + x2);
ym = jet->correctedEt(xm);
if(ym < 0) {
y2 = ym;
x2 = xm;
} else {
x1 = xm;
y1 = ym;
}
++i;
}
//std::cout << "I:" << i << ", " << y1 << "," << y2 << '\n';
expectedEt = 0.5 * (x1 + x2);
}
//expectedEt = jet->Et();
double scale = truth;
double expectedEt = jet->expectedEt(truth,scale);
//calculate chi2
double err2inv = jet->Error();
err2inv *= err2inv;
err2inv = 1/err2inv;
double et = jet->correctedEt(expectedEt);
double hadf = jet->HadEt()/jet->Et();
double s = et/expectedEt;
et = s * jet->Et();
double s2inv =hadf/s;
s2inv *= s2inv;
double chi2 = et - truth;
chi2 *= chi2 * err2inv * s2inv;
double chi2 = jet->Et() - expectedEt;
chi2 *= chi2 * err2inv;
chi2 = weight * TData::ScaleResidual(chi2);
//calculate chi2 for derivatives
int npar = jet->nPar();
double et1,et2,temp1,temp2;
for(int i = 0 ; i < npar ; ++i) {
int parid = jet->varyPar(i,epsilon,expectedEt,et2,et1);
s = et1/expectedEt;
et1 = s * jet->Et();
temp1 = et1 - truth;
//s2inv = hadf/s;
//s2inv *= s2inv;
temp1 *= temp1 * err2inv * s2inv;
int parid = jet->varyPar(i,epsilon,truth,scale,et2,et1);
//std::cout << truth << ", " << expectedEt << "," << et1 << ", " << et2 << std::endl;
if((std::abs((et1 - expectedEt)/scale) > 0.05) ||
(std::abs((et2 - expectedEt)/scale) > 0.05)) {
std::cout << "strange extrapolation result:" << expectedEt << "; "
<< et1 << "; " << et2 << " jet:" << jet->Et()
<< " truth:" << truth << std::endl;
}
temp1 = et1 - jet->Et();
temp1 *= temp1 * err2inv;
temp1 = weight * TData::ScaleResidual(temp1);
s = et2/expectedEt;
et2 = s * jet->Et();
temp2 = et2 - truth;
//s2inv = hadf/s;
//s2inv *= s2inv;
temp2 *= temp2 * err2inv * s2inv;
temp2 = et2 - jet->Et();
temp2 *= temp2 * err2inv;
temp2 = weight * TData::ScaleResidual(temp2);
temp_derivative1[parid] += (temp2 - temp1); // for 1st derivative
temp_derivative2[parid] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative
Expand Down
123 changes: 123 additions & 0 deletions JetWithTowers.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
//
// Class for jets with towers
//
// first version: Hartmut Stadie 2008/12/25
// $Id: Jet.h,v 1.1 2008/12/16 15:21:27 stadie Exp $
//
#include"JetWithTowers.h"

#include "TLorentzVector.h"

JetWithTowers::JetWithTowers(double Et, double EmEt, double HadEt,
double OutEt, double E,double eta,double phi,
Flavor flavor,
double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars)
: Jet(Et,EmEt,HadEt,OutEt,E,eta,phi,flavor,func,err,firstpar,id,npars),
njetpars(npars),ntowerpars(0)
{
}

JetWithTowers::~JetWithTowers()
{
for(TowerCollIter i = towers.begin() ; i != towers.end() ; ++i) {
delete *i;
}
}

void JetWithTowers::ChangeParAddress(double* oldpar, double* newpar)
{
Jet::ChangeParAddress(oldpar,newpar);
for(TowerCollIter i = towers.begin() ; i != towers.end() ; ++i) {
(*i)->ChangeParAddress(oldpar,newpar);
towerpars[(*i)->FirstPar()] = (*i)->Par();
}
}

double JetWithTowers::correctedEt(double Et) 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);
}
//std::cout << "scale ET:" << Et << " cor. sum of tower:" << ccet << "\n";
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 << "alternating par:" << id + towpar << " = " << p[towpar] << std::endl;
double orig = p[towpar];
p[towpar] += eps;
upperEt = expectedEt(Et,scale,true);
p[towpar] = orig - eps;
lowerEt = expectedEt(Et,scale,true);
p[towpar] = orig;
return id + towpar;
}



void JetWithTowers::addTower(double Et, double EmEt, double HadEt ,
double OutEt, double E,double eta,double phi,
double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars)
{
TLorentzVector jet, towp;
jet.SetPtEtaPhiM(TMeasurement::pt,TMeasurement::eta,TMeasurement::phi,0);
towp.SetPtEtaPhiM(Et,eta,phi,0);
jet -= towp;
double projection = (TMeasurement::pt-jet.Et())/Et;
projection = 1.0;
//std::cout << "Jet:" << jet.Eta() << ", " << jet.Phi()
// << " tower:" << towp.Eta() << ", " << towp.Phi() << " :"
// << projection << std::endl;
towers.push_back(new Tower(Et,EmEt,HadEt,OutEt,E,eta,phi,projection,func,
err,firstpar,id,npars));
ntowerpars = npars;
towerpars[id] = firstpar;
Jet::npar = njetpars + towerpars.size() * ntowerpars;
//std::cout << "tower:" << Et << " projected:" << projection * Et << std::endl;
}



JetWithTowers::Tower::Tower(double Et, double EmEt, double HadEt ,
double OutEt, double E,double eta,double phi,
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)
{
temp = *this;
}

double JetWithTowers::Tower::correctedEt(double Et) 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);
//std::cout << pt << ", " << Et << ":" << lastCorEt << " par:" << par[0] << std::endl;
return lastCorEt;
}

71 changes: 71 additions & 0 deletions JetWithTowers.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
//
// Class for jets with towers
//
// first version: Hartmut Stadie 2008/12/25
// $Id: Jet.h,v 1.1 2008/12/16 15:21:27 stadie Exp $
//
#ifndef JETWITHTOWERS_H
#define JETWITHTOWERS_H

#include"Jet.h"

#include <vector>
#include <map>

class JetWithTowers : public Jet
{
public:
JetWithTowers(double Et, double EmEt, double HadEt ,double OutEt, double E,
double eta,double phi, Flavor flavor,
double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars);
virtual ~JetWithTowers();
virtual void ChangeParAddress(double* oldpar, double* newpar);
virtual double correctedEt(double Et) 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);
void addTower(double Et, double EmEt, double HadEt ,double OutEt, double E,
double eta,double phi,
double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars);
protected:
class Tower : public TMeasurement {
public:
Tower(double Et, double EmEt, double HadEt ,double OutEt, double E,
double eta,double phi, double alpha,
double const(*func)(TMeasurement *const x, double *const par),
double err, double* firstpar, int id, int npars);
double Et() const {return pt;}
double EmEt() const {return EMF;}
double HadEt() const {return HadF;}
double OutEt() const {return OutF;}
double E() const {return TMeasurement::E;}
double eta() const {return TMeasurement::eta;}
double phi() const {return TMeasurement::phi;}
double projectionToJetAxis() const {return alpha;}
void ChangeParAddress(double* oldpar, double* newpar) {par += newpar - oldpar;}
double correctedEt(double Et) const;
double lastCorrectedEt() const { return lastCorEt;}
double Error() const {return error;}
int nPar() const {return npar;}
int FirstPar() const {return parid;}
double *Par() const {return par;}
private:
double alpha;
double* par;//address to first parameter for this jet
int npar,parid;
double error;
mutable TMeasurement temp;
mutable double lastCorEt;
double const(*f)(TMeasurement *const x, double *const par);
};
typedef std::vector<Tower*> TowerColl;
typedef TowerColl::iterator TowerCollIter;
typedef TowerColl::const_iterator TowerCollConstIter;
TowerColl towers;
int njetpars,ntowerpars;
std::map<int,double*> towerpars;
};

#endif
9 changes: 6 additions & 3 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ RCXX=$(SPECIALFLAGS) -Wno-deprecated -Wall $(ROOTCFLAGS)
RLXX=$(LFLAGS) $(ROOTLIBS) -lboost_thread -lpthread #-lrt -lpthread # -lposix4


SRC=caliber.cc GammaJetSel.cc ZJetSel.cc TrackClusterSel.cc NJetSel.cc TopSel.cc ConfigFile.cc CalibData.cc Parameters.cc ControlPlots.cc ToyMC.cc EventReader.cc PhotonJetReader.cc DiJetReader.cc TriJetReader.cc ZJetReader.cc TopReader.cc ParameterLimitsReader.cc TowerConstraintsReader.cc TrackClusterReader.cc EventProcessor.cc Jet.cc JetTruthEvent.cc
SRC=caliber.cc GammaJetSel.cc ZJetSel.cc TrackClusterSel.cc NJetSel.cc TopSel.cc ConfigFile.cc CalibData.cc Parameters.cc ControlPlots.cc ToyMC.cc EventReader.cc PhotonJetReader.cc DiJetReader.cc TriJetReader.cc ZJetReader.cc TopReader.cc ParameterLimitsReader.cc TowerConstraintsReader.cc TrackClusterReader.cc EventProcessor.cc Jet.cc JetTruthEvent.cc JetWithTowers.cc

%.o: %.cc
$(C) $(RCXX) -c $<
Expand Down Expand Up @@ -65,7 +65,7 @@ ControlPlots.o: ControlPlots.cc ControlPlots.h CalibData.h CalibMath.h ConfigFil
EventReader.o: EventReader.h EventReader.cc ConfigFile.h
$(C) $(CFLAGS) -c EventReader.cc

PhotonJetReader.o: EventReader.h PhotonJetReader.h PhotonJetReader.cc GammaJetSel.h ToyMC.h Parameters.h ConfigFile.h Jet.h JetTruthEvent.h
PhotonJetReader.o: EventReader.h PhotonJetReader.h PhotonJetReader.cc GammaJetSel.h ToyMC.h Parameters.h ConfigFile.h Jet.h JetTruthEvent.h JetWithTowers.h
$(C) $(RCXX) -c PhotonJetReader.cc

DiJetReader.o: EventReader.h DiJetReader.h DiJetReader.cc NJetSel.h ToyMC.h Parameters.h ConfigFile.h
Expand Down Expand Up @@ -98,7 +98,10 @@ Jet.o: CalibData.h Jet.h Jet.cc Parametrization.h
JetTruthEvent.o: CalibData.h Jet.h JetTruthEvent.h JetTruthEvent.cc
$(C) $(CFLAGS) -c JetTruthEvent.cc

caliber.o: caliber.cc caliber.h CalibMath.h external.h ConfigFile.h CalibData.h Parameters.h ControlPlots.h EventReader.h DiJetReader.h TriJetReader.h ZJetReader.h TopReader.h ParameterLimitsReader.h TowerConstraintsReader.h TrackClusterReader.h EventProcessor.h
JetWithTowers.o: CalibData.h Jet.h JetWithTowers.h JetWithTowers.cc Parametrization.h
$(C) $(RCXX) -c JetWithTowers.cc

caliber.o: caliber.cc caliber.h CalibMath.h external.h ConfigFile.h CalibData.h Parameters.h ControlPlots.h EventReader.h DiJetReader.h TriJetReader.h ZJetReader.h TopReader.h ParameterLimitsReader.h TowerConstraintsReader.h TrackClusterReader.h EventProcessor.h Jet.h
$(C) $(RCXX) -I/usr/include/boost -c caliber.cc

runjunk: $(SRC:.cc=.o) lbfgs.o
Expand Down
Loading

0 comments on commit 9ca80f1

Please sign in to comment.