implement new way for constraint using the same jets as the fit

DiJetReader.cc

@@ -1,6 +1,6 @@
 //
 //    first version: Hartmut Stadie 2008/12/12
-//    $Id: DiJetReader.cc,v 1.36 2010/01/25 17:35:20 stadie Exp $
+//    $Id: DiJetReader.cc,v 1.37 2010/01/25 17:51:35 stadie Exp $
 //   
 #include "DiJetReader.h"
 
@@ -12,6 +12,7 @@
 #include "Jet.h"
 #include "JetWithTowers.h"
 #include "TwoJetsPtBalanceEvent.h"
+#include "JetConstraintEvent.h"
 #include "NJetSel.h"
 #include "CorFactors.h"
 #include "CorFactorsFactory.h"
@@ -395,6 +396,10 @@ int DiJetReader::createJetTruthEvents(std::vector<Event*>& data)
     JetTruthEvent* jte = new JetTruthEvent(jet,nJet_->GenJetColEt[genJetIdx],1.);//nJet_->Weight);
     data.push_back(jte);
     ++njets;
+    //add jet to constraints
+    for(unsigned int i = 0 ; i < constraints_.size() ; ++i) {
+      constraints_[i]->addJet(jet,new Function(&Parametrization::correctedJetEt,0,0,0,0,cp_));
+    }
   }     
   delete [] terr;
   return njets;

EventReader.cc

@@ -1,21 +1,24 @@
 //
-// $Id: EventReader.cc,v 1.6 2010/01/08 18:23:28 mschrode Exp $
+// $Id: EventReader.cc,v 1.7 2010/01/28 16:06:22 stadie Exp $
 //
 #include "EventReader.h"
 
 #include "ConfigFile.h"
 #include "Parameters.h" 
+#include "Parametrization.h"
 #include "CorFactorsFactory.h"
+#include "JetConstraintEvent.h" 
 #include "TChain.h"
 #include "ToyMC.h"
 #include "TTree.h"
 
 #include <dlfcn.h>
 
 unsigned int EventReader::numberOfEventReaders_ = 0;
+std::vector<JetConstraintEvent*> EventReader::constraints_;
 
 EventReader::EventReader(const std::string& configfile, TParameters* param) 
-  : config_(0),par_(param),corFactorsFactory_(0)
+  : config_(0),par_(param),corFactorsFactory_(0),cp_(new ConstParametrization())
 {
   numberOfEventReaders_++;
 
@@ -84,12 +87,32 @@ EventReader::EventReader(const std::string& configfile, TParameters* param)
   correctToL3_ = config_->read<bool>("correct jets to L3",false);
   if(correctToL3_) {
     std::cout << "Jets will be corrected to Level3." << std::endl;
+  } 
+  if(! constraints_.size() ) {
+    std::vector<double> jet_constraint = bag_of<double>(config_->read<std::string>( "jet constraints",""));
+    if(jet_constraint.size() % 5 == 0) {
+      for(unsigned int i = 0 ; i < jet_constraint.size() ; i += 5) {
+	constraints_.push_back(new JetConstraintEvent(jet_constraint[i],jet_constraint[i+1],jet_constraint[i+2],jet_constraint[i+3],jet_constraint[i+4]));
+      } 
+    } else if(jet_constraint.size() > 1) {
+      std::cout << "wrong number of arguments for jet constraint:" << jet_constraint.size() << '\n';
+    }
+    for(unsigned int i = 0 ; i < constraints_.size() ; ++i) {
+      const JetConstraintEvent* jce = constraints_[i];
+      std::cout << "adding constraint for jets with " << jce->minEta() << " < |eta| <  " 
+		<< jce->maxEta() << " and " << jce->minPt() << " < pt < " << jce->maxPt() 
+		<< " with weight " << jce->weight() << "\n";
+    }
   }
 }
 
 EventReader::~EventReader()
 {
   delete config_;
+  for(unsigned int i = 0 ; i < constraints_.size() ; ++i) {
+    delete constraints_[i];
+  }
+  constraints_.clear();
 }
 
 
@@ -167,3 +190,16 @@ TTree * EventReader::createTree(const std::string &dataType) const {
   return tree;
 }
 
+int EventReader::addConstraints(std::vector<Event*>& data) {
+  unsigned int n = constraints_.size();
+  for(unsigned int i = 0 ; i < n ; ++i) { 
+    JetConstraintEvent* jce = constraints_[i];
+    std::cout << "added constraint for jets with " << jce->minEta() << " < |eta| <  " 
+	      << jce->maxEta() << " and " << jce->minPt() << " < pt < " << jce->maxPt() 
+	      << " with weight " << jce->weight() << " and " << jce->nJets() 
+	      << " jets " << "\n";
+    data.push_back(jce);
+  }
+  constraints_.clear();
+  return n;
+}

EventReader.h

@@ -1,5 +1,5 @@
 //
-// $Id: EventReader.h,v 1.8 2009/11/27 15:28:12 stadie Exp $
+// $Id: EventReader.h,v 1.9 2010/01/08 18:23:28 mschrode Exp $
 //
 #ifndef EVENTREADER_H
 #define EVENTREADER_H
@@ -11,10 +11,13 @@ class Measurement;
 class CorFactors;
 class CorFactorsFactory;
 class TTree;
+class JetConstraintEvent;
+class Parametrization;
 
 #include <vector>
 #include <string>
 
+
 class EventReader
 {
  public:
@@ -24,6 +27,8 @@ class EventReader
   virtual ~EventReader();
   virtual int readEvents(std::vector<Event*>& data) = 0;
 
+  static int addConstraints(std::vector<Event*>& data);
+
  protected:
   //! Read CorFactors from Ntuple
   virtual CorFactors* createCorFactors(int jetid) const { return 0;}
@@ -38,7 +43,10 @@ class EventReader
 
   double (*tower_error_param)(const double *x, const Measurement *xorig, double err);
   double (*jet_error_param)  (const double *x, const Measurement *xorig, double err);
-  double (*track_error_param)(const double *x, const Measurement *xorig, double err);
+  double (*track_error_param)(const double *x, const Measurement *xorig, double err);  
+
+  static std::vector<JetConstraintEvent*> constraints_;
+  Parametrization *cp_;
 };
 
 

Function.h

@@ -2,7 +2,7 @@
 //    Class representing a correction function
 //
 //    first version: Hartmut Stadie 2008/12/14
-//    $Id: JetTruthEvent.cc,v 1.17 2009/07/23 15:51:17 stadie Exp $
+//    $Id: Function.h,v 1.4 2009/11/24 16:52:59 stadie Exp $
 //   
 #ifndef FUNCTION_H
 #define FUNCTION_H
@@ -30,8 +30,8 @@ class Function {
   bool hasInverse() { return invfunc_;}
   double inverse(const Measurement* x) const { return invfunc_ ? (param_->*invfunc_)(x,firstpar_) : 0;}
  private:
-  const ParametrizationFunction func_;
-  const ParametrizationFunction invfunc_;
+  ParametrizationFunction func_;
+  ParametrizationFunction invfunc_;
   double *firstpar_;
   int parindex_;
   int npars_;

Jet.cc

@@ -2,7 +2,7 @@
 //    Class for basic jets 
 //
 //    first version: Hartmut Stadie 2008/12/14
-//    $Id: Jet.cc,v 1.36 2009/11/27 15:28:12 stadie Exp $
+//    $Id: Jet.cc,v 1.37 2010/01/25 17:35:20 stadie Exp $
 //   
 #include "Jet.h"  
 
@@ -24,6 +24,15 @@ Jet::Jet(double Et, double EmEt, double HadEt ,double OutEt, double E,
   varcoll.resize(f.nPars() + gf.nPars());
 }
 
+Jet::Jet(const Jet& j) 
+  : Measurement(j), flavor_(j.flavor_), genPt_(j.genPt_),dR_(j.dR_), 
+    ptHat_(j.ptHat_),corFactors_(new CorFactors(*(j.corFactors_))),f(j.f),
+    gf(j.gf),errf(j.errf),etmin(j.etmin),EoverPt(j.EoverPt),gsl_impl(this)
+{
+  temp = *this;
+  varcoll.resize(f.nPars() + gf.nPars());
+}
+
 Jet::~Jet()
 {
   delete corFactors_;

Jet.h

@@ -4,7 +4,7 @@
 //!
 //!    \date 2008/12/14
 //!
-//!    $Id: Jet.h,v 1.29 2010/01/12 19:24:48 mschrode Exp $
+//!    $Id: Jet.h,v 1.30 2010/01/25 17:35:20 stadie Exp $
 #ifndef JET_H
 #define JET_H
 
@@ -137,10 +137,13 @@ class Jet : public Measurement
   static void printInversionStats();  //!< Print some info on inversion
 
   int parIndex() const { return f.parIndex(); }
+
+  virtual Jet* clone() const { return new Jet(*this);} //!< Clone this jet
+  void setGlobalFunction(const Function& ngf) { gf = ngf;} //!< Set global correction function, needed for constraints
  protected:
   mutable VariationColl varcoll;
   virtual double expectedEt(double truth, double start, bool fast = false);
-
+  Jet(const Jet&j); //!< disallow copies!
  private: 
   Flavor flavor_;           //!< The jet's Flavor
   double genPt_;            //!< The genjet pt

JetConstraintEvent.cc

@@ -2,7 +2,7 @@
 //    Class for constraints on the jet correction
 //
 //    first version: Hartmut Stadie 2009/07/23
-//    $Id: JetConstraintEvent.cc,v 1.2 2009/07/24 09:38:34 stadie Exp $
+//    $Id: JetConstraintEvent.cc,v 1.3 2010/01/28 16:07:27 stadie Exp $
 //   
 
 
@@ -12,75 +12,85 @@
 //interface to Data
 JetConstraintEvent::~JetConstraintEvent()
 {
-  for(unsigned int i = 0, njets = jets.size() ; i < njets ; ++i) {
-    delete jets[i];
+  for(unsigned int i = 0, njets = jets_.size() ; i < njets ; ++i) {
+    delete jets_[i];
   }
-  jets.clear();
+  jets_.clear();
 }
 
-void JetConstraintEvent::addJet(Jet* j) 
+void JetConstraintEvent::addJet(const Jet* j, const Function* globalFunc) 
 {
-  jets.push_back(j); 
-  trusum += j->Et();
+  if((j->pt() > maxpt_) || (j->pt() < minpt_) || 
+     (std::abs(j->eta()) > maxeta_) || (std::abs(j->eta()) < mineta_)) return; 
+
+  Jet* jet = j->clone();
+  if(globalFunc) jet->setGlobalFunction(*globalFunc);
+  jets_.push_back(jet); 
+  trusum_ += jet->Et();
   //add parameter ids to variation  map
-  const Jet::VariationColl& varcoll = j->varyParsDirectly(0.001);
+  const Jet::VariationColl& varcoll = jet->varyParsDirectly(0.001);
   for(Jet::VariationCollIter i = varcoll.begin() ; i != varcoll.end() ; ++i) {
-    varmap[i->parid] = Variation();
+    varmap_[i->parid] = Variation();
   }
 }
 
 
 void JetConstraintEvent::ChangeParAddress(double* oldpar, double* newpar) { 
-  for(unsigned int i = 0, njets = jets.size() ; i < njets ; ++i) {
-    jets[i]->ChangeParAddress(oldpar,newpar);
+  for(unsigned int i = 0, njets = jets_.size() ; i < njets ; ++i) {
+    jets_[i]->ChangeParAddress(oldpar,newpar);
   }
 }
 
 double JetConstraintEvent::chi2() const
 {
   double sumet = 0;
-  unsigned int njets = jets.size();
+  unsigned int njets = jets_.size();
+  if(! njets) return 0;
   for(unsigned int i = 0 ; i < njets ; ++i) {
-    sumet += jets[i]->correctedEt(jets[i]->Et());
+    sumet += jets_[i]->correctedEt(jets_[i]->Et());
   }
-  double chi2 = sumet / trusum -1;
-  return weight * chi2 * chi2;
+  double chi2 = sumet / trusum_ -1;
+  return weight_ * chi2 * chi2;
 }
 
 double JetConstraintEvent::chi2_fast(double * temp_derivative1, double * temp_derivative2, double const epsilon) const
 {
   double sumet = 0;
-  unsigned int njets = jets.size();
+  unsigned int njets = jets_.size();
+  if(! njets) {
+    chi2plots_ = 0;
+    return 0;
+  }
   for(unsigned int i = 0 ; i < njets ; ++i) {
-    sumet += jets[i]->correctedEt(jets[i]->Et());
+    sumet += jets_[i]->correctedEt(jets_[i]->Et());
   }
-  double chi2 = sumet / trusum - 1;
+  double chi2 = sumet / trusum_ - 1;
   chi2 *= chi2;
-  chi2 *= weight;
+  chi2 *= weight_;
   //std::cout << "JetConstraintEvent::chi2_fast:  Et:" << jets[0]->Et() 
   //   	    << " chi2:" << chi2 << " sum Et:" << sumet << std::endl;
   //derivative....
-  for(VarMap::iterator i = varmap.begin() ; i != varmap.end() ; ++i) {
-    i->second.uppersum = sumet;
-    i->second.lowersum = sumet;
+  for(VarMap::iterator i = varmap_.begin() ; i != varmap_.end() ; ++i) {
+    i->second.uppersum_ = sumet;
+    i->second.lowersum_ = sumet;
   }
   for(unsigned int i = 0 ; i < njets ; ++i) {
-    const Jet::VariationColl& varcoll = jets[i]->varyParsDirectly(epsilon);
-    double et = jets[i]->correctedEt(jets[i]->Et());
+    const Jet::VariationColl& varcoll = jets_[i]->varyParsDirectly(epsilon);
+    double et = jets_[i]->correctedEt(jets_[i]->Et());
     for(Jet::VariationCollIter j = varcoll.begin() ; j != varcoll.end() ; ++j) {
-      VarMap::iterator k = varmap.find(j->parid);
-      assert(k != varmap.end());
-      k->second.uppersum += j->upperEt - et;
-      k->second.lowersum += j->lowerEt - et;
+      VarMap::iterator k = varmap_.find(j->parid);
+      assert(k != varmap_.end());
+      k->second.uppersum_ += j->upperEt - et;
+      k->second.lowersum_ += j->lowerEt - et;
     }
   }
-  for(VarMap::iterator i = varmap.begin() ; i != varmap.end() ; ++i) {
-    double temp1 = i->second.lowersum / trusum - 1;
+  for(VarMap::iterator i = varmap_.begin() ; i != varmap_.end() ; ++i) {
+    double temp1 = i->second.lowersum_ / trusum_ - 1;
     temp1 *= temp1;
-    temp1 *= weight;
-    double temp2 = i->second.uppersum / trusum - 1;
+    temp1 *= weight_;
+    double temp2 = i->second.uppersum_ / trusum_ - 1;
     temp2 *= temp2;
-    temp2 *= weight;
+    temp2 *= weight_;
 //     std::cout << "JetConstraintEvent::chi2_fast:  Et:" << jets[0]->Et() 
 // 	      << " temp1:" << temp1 << " sum Et:" <<  i->second.lowersum 
 // 	      << std::endl;
@@ -90,7 +100,7 @@ double JetConstraintEvent::chi2_fast(double * temp_derivative1, double * temp_de
     temp_derivative1[i->first] += (temp2 - temp1); // for 1st derivative
     temp_derivative2[i->first] += (temp2 + temp1 - 2 * chi2); // for 2nd derivative   
   }
-  chi2plots = chi2;
+  chi2plots_ = chi2;
   return chi2;
 }