diff --git a/JetTruthEvent.cc b/JetTruthEvent.cc index edd3ca0..2d8b30f 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.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" @@ -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 diff --git a/JetWithTowers.cc b/JetWithTowers.cc new file mode 100644 index 0000000..74cbfd3 --- /dev/null +++ b/JetWithTowers.cc @@ -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::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; +} + diff --git a/JetWithTowers.h b/JetWithTowers.h new file mode 100644 index 0000000..1c1c40b --- /dev/null +++ b/JetWithTowers.h @@ -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 +#include + +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 TowerColl; + typedef TowerColl::iterator TowerCollIter; + typedef TowerColl::const_iterator TowerCollConstIter; + TowerColl towers; + int njetpars,ntowerpars; + std::map towerpars; +}; + +#endif diff --git a/Makefile b/Makefile index a1804a1..d7f7c40 100644 --- a/Makefile +++ b/Makefile @@ -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 $< @@ -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 @@ -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 diff --git a/PhotonJetReader.cc b/PhotonJetReader.cc index 9b49a6f..0688847 100644 --- a/PhotonJetReader.cc +++ b/PhotonJetReader.cc @@ -4,13 +4,14 @@ // This class reads events according fo the GammaJetSel // // first version: Hartmut Stadie 2008/12/12 -// $Id: PhotonJetReader.cc,v 1.4 2008/12/16 15:21:26 stadie Exp $ +// $Id: PhotonJetReader.cc,v 1.5 2008/12/17 08:16:04 stadie Exp $ // #include "PhotonJetReader.h" #include "CalibData.h" #include "JetTruthEvent.h" #include "Jet.h" +#include "JetWithTowers.h" #include "ConfigFile.h" #include "ToyMC.h" #include "Parameters.h" @@ -36,7 +37,7 @@ PhotonJetReader::PhotonJetReader(const std::string& configfile, TParameters* p) Rel_cut_on_gamma = config->read("Relative Rest Jet Cut",0.2); dataClass = config->read("Gamma-Jet data class", 0); - if((dataClass != 0) && (dataClass != 1)) dataClass = 0; + if((dataClass < 0) || (dataClass > 2)) dataClass = 0; TTree* tchain_gammajet; vector input_gammajet = @@ -93,7 +94,7 @@ int PhotonJetReader::readEvents(std::vector& data) TData* ev = 0; if(dataClass == 0) ev = createTruthMultMessEvent(); - else if(dataClass == 1) ev = createJetTruthEvent(); + else if(dataClass > 0) ev = createJetTruthEvent(); if(ev) { @@ -113,7 +114,7 @@ TData* PhotonJetReader::createJetTruthEvent() double out = 0; double err2 = 0; TMeasurement tower; - double terr; + double* terr = new double[gammajet.NobjTowCal]; for(int n = 0; n < gammajet.NobjTowCal; ++n) { em += gammajet.TowEm[n]; had += gammajet.TowHad[n]; @@ -126,9 +127,10 @@ TData* PhotonJetReader::createJetTruthEvent() tower.eta = gammajet.TowEta[n]; tower.phi = gammajet.TowPhi[n]; tower.E = gammajet.TowE[n]; - terr = sqrt(1.3 * 1.3/gammajet.TowHad[n] + 0.056 * 0.056) * tower.HadF; - double err = tower_error_param(&tower.pt,&tower,terr); - err2 += err * err; + terr[n] = sqrt(1.3 * 1.3/gammajet.TowHad[n] + 0.056 * 0.056) * tower.HadF; + double err = tower_error_param(&tower.pt,&tower,terr[n]); + terr[n] = err * err; + err2 += terr[n]; } //calc jet error double factor = gammajet.JetCalEt / gammajet.JetCalE; @@ -145,11 +147,36 @@ TData* PhotonJetReader::createJetTruthEvent() int jet_index = p->GetJetBin(p->GetJetEtaBin(gammajet.TowId_eta[0]), p->GetJetPhiBin(gammajet.TowId_phi[0])); double* firstpar = p->GetJetParRef(jet_index); - Jet *j = new Jet(gammajet.JetCalEt,em * factor,had * factor,out * factor,gammajet.JetCalEta, - gammajet.JetCalPhi,gammajet.JetCalE,TJet::uds,p->jet_parametrization,sqrt(err2), - firstpar,firstpar - p->GetPars(),p->GetNumberOfJetParametersPerBin()); - + Jet *j; + if(dataClass == 2) { + JetWithTowers *jt = + new JetWithTowers(gammajet.JetCalEt,em * factor,had * factor, + out * factor,gammajet.JetCalE,gammajet.JetCalEta, + gammajet.JetCalPhi,TJet::uds, + p->jet_parametrization,sqrt(err2), + firstpar,firstpar - p->GetPars(), + p->GetNumberOfJetParametersPerBin()); + for(int i = 0; i < gammajet.NobjTowCal; ++i) { + double scale = gammajet.TowEt[i]/gammajet.TowE[i]; + int id = p->GetBin(p->GetEtaBin(gammajet.TowId_eta[i]), + p->GetPhiBin(gammajet.TowId_phi[i])); + jt->addTower(gammajet.TowEt[i],gammajet.TowEm[i]*scale, + gammajet.TowHad[i]*scale,gammajet.TowOE[i]*scale, + gammajet.TowE[i],gammajet.TowEta[i],gammajet.TowPhi[i], + p->tower_parametrization,terr[i],p->GetTowerParRef(id), + id,p->GetNumberOfTowerParametersPerBin()); + } + j = jt; + } + else { + j = new Jet(gammajet.JetCalEt,em * factor,had * factor,out * factor, + gammajet.JetCalE,gammajet.JetCalEta,gammajet.JetCalPhi, + TJet::uds,p->jet_parametrization,sqrt(err2), + firstpar,firstpar - p->GetPars(), + p->GetNumberOfJetParametersPerBin()); + } JetTruthEvent* jte = new JetTruthEvent(j,gammajet.PhotonEt,gammajet.EventWeight); + delete [] terr; return jte; }