Factorized efficiency from energy smearer.

Added di-photon smearers with systematics in StatAnalysis.
Moved vertex reweghing to di-photon smearer.

BaseSmearer.cc

@@ -7,3 +7,11 @@ BaseSmearer::~BaseSmearer()
 {}
 
 bool operator == (BaseSmearer * a, const std::string & b) { return a->name() == b; };
+
+BaseDiPhotonSmearer::BaseDiPhotonSmearer()
+{}
+
+BaseDiPhotonSmearer::~BaseDiPhotonSmearer() 
+{}
+
+bool operator == (BaseDiPhotonSmearer * a, const std::string & b) { return a->name() == b; };

BaseSmearer.h

@@ -31,10 +31,34 @@ class BaseSmearer
 
 	// ! return smeared photon informations
 	virtual bool smearPhoton( PhotonReducedInfo & pho, float & weight, float syst_shift=0. ) const = 0;
-
+
+	/// virtual bool smearDiPhoton( TLorentzVector & 4p, const TVector3 & selVtx, const TVector3 & trueVtx);
 };
 
 // ! Used to search analyzers by name 
 bool operator == (BaseSmearer * a, const std::string & b);
 
+
+class BaseDiPhotonSmearer 
+{
+public:
+	// ! C-TOR
+	BaseDiPhotonSmearer();
+
+	// ! D-TOR
+	virtual ~BaseDiPhotonSmearer();
+
+	// ! the class name
+	virtual const std::string & name() const = 0; 
+
+	// !  
+	operator const std::string & () const { return this->name(); };
+
+	virtual bool smearDiPhoton( TLorentzVector & p4, TVector3 & selVtx, float & weight, const int & category, 
+				    const int & genMassPoint, const TVector3 & trueVtx, float syst_shift=0.) const = 0 ;
+};
+
+// ! Used to search analyzers by name 
+bool operator == (BaseDiPhotonSmearer * a, const std::string & b);
+
 #endif

DiPhoEfficiencySmearer.cc

@@ -0,0 +1,128 @@
+#include "DiPhoEfficiencySmearer.h"
+#include "PhotonReducedInfo.h"
+#include "TRandom3.h"
+#include <assert.h>
+
+DiPhoEfficiencySmearer::DiPhoEfficiencySmearer(const diPhoEfficiencySmearingParameters& par) : myParameters_(par)
+{
+  rgen_ = new TRandom3(0);
+  name_="DiPhoEfficiencySmearer_"+ par.categoryType + "_" + par.parameterSetName;
+
+}
+
+DiPhoEfficiencySmearer::~DiPhoEfficiencySmearer()
+{
+  delete rgen_;
+}
+
+bool DiPhoEfficiencySmearer::smearDiPhoton( TLorentzVector & p4, TVector3 & selVtx, float & weight, const int & category, 
+					    const int & genMassPoint, const TVector3 & trueVtx, float syst_shift) const
+{
+  std::string cat=Form("cat%d", category);
+
+  if (cat == "")
+    {
+      std::cout << "No category has been found associated with this diphoton. Giving Up" << std::endl;
+      return false;
+    }
+
+  /////////////////////// changing weigh of photon according to efficiencies ///////////////////////////////////////////
+  assert( ! smearing_eff_graph_.empty() );
+  if( doVtxEff_ ) {
+    cat += (selVtx - trueVtx).Mag() > 1. ? "_pass" : "_fail"; 
+  }
+
+  weight = getWeight( p4.Pt(), cat, syst_shift );
+
+  return true;
+}
+
+
+
+bool DiPhoEfficiencySmearer::init() 
+{
+
+  // if map is not empty, yuo're initilized and happy..
+  if( !smearing_eff_graph_.empty() ){
+     std:cout << "initialization of efficiencies already done; proceed with usage. " << std::endl;
+    return true;
+  }
+  if( doVtxEff_ ) { passFailWeights_ = true; }
+
+  //otherwise, get smearing functions from file and set up map
+  std::cout << "\n>>>initializing one efficiency for photon re-weighting; " <<  std::endl;
+
+  // do basic sanity checks first
+  if( effName_.empty()){
+    std::cout << "you're initialzinfg reweighting for efficiency but effName_ is empty" << std::endl;  assert(false); }
+  if( myParameters_.efficiency_file.empty()){
+    std::cout << "you're initialzinfg reweighting for efficiency: " << effName_  << " but input file with TGraphErrors is not specified; doing nothing. " << std::endl;  assert(false); }
+
+  theDiPhoEfficiencyFile_ = new TFile(myParameters_.efficiency_file.c_str());
+
+  // initialize formulas for the di-photon categories; 
+  for( int ii=0; ii<myParameters_.n_categories; ++ii ) {
+    std::string cat = Form("cat%d", ii);
+    if( passFailWeights_ ) {
+      smearing_eff_graph_[cat+"_pass"]=(TGraphAsymmErrors*) theDiPhoEfficiencyFile_->Get((effName_+"_"+cat+"_pass").c_str())->Clone();
+      smearing_eff_graph_[cat+"_fail"]=(TGraphAsymmErrors*) theDiPhoEfficiencyFile_->Get((effName_+"_"+cat+"_fail").c_str())->Clone();
+      std::cerr << "DiPhoEfficiencySmearerc " << cat+"_pass" << std::endl;
+    } else {
+      smearing_eff_graph_[cat]=(TGraphAsymmErrors*) theDiPhoEfficiencyFile_->Get((effName_+"_"+cat).c_str())->Clone();
+    }
+  }
+
+  delete  theDiPhoEfficiencyFile_;
+  return true;
+
+}
+
+
+double DiPhoEfficiencySmearer::getWeight(double pt, std::string theCategory, float syst_shift) const
+{
+  std::map<std::string,TGraphAsymmErrors*>::const_iterator theIter = smearing_eff_graph_.find(theCategory);
+  if( theIter != smearing_eff_graph_.end()  ) {
+
+    // determine the pair of bins between which  you interpolate
+    int numPoints = ( theIter->second )->GetN();
+    double x, y;
+    int myBin = -1;
+    for (int bin=0; bin<numPoints; bin++ ){
+      ( theIter->second )->GetPoint(bin, x, y);
+      if(pt > x) {
+	myBin = bin; }
+      else break;
+    }
+    int binLow, binHigh; 
+    if(myBin == -1)                      {binHigh = 0; binLow=0;}
+    else if (myBin == (numPoints-1))     {binHigh = numPoints-1; binLow=numPoints-1;}
+    else {binLow=myBin; binHigh=myBin+1;}
+
+
+    // get hold of efficiency ratio and error at either points
+    // low-high refer to the points ; up-down refers to the errors 
+    double xLow, yLow;    double xHigh, yHigh;
+    ( theIter->second )->GetPoint(binLow, xLow, yLow);
+    ( theIter->second )->GetPoint(binHigh, xHigh, yHigh);
+
+    double errLowYup    = ( theIter->second )->GetErrorYhigh(binLow);
+    double errLowYdown  = ( theIter->second )->GetErrorYlow(binLow);
+    double errHighYup   = ( theIter->second )->GetErrorYhigh(binHigh);
+    double errHighYdown = ( theIter->second )->GetErrorYlow(binHigh);
+
+    double theErrorLow, theErrorHigh;
+    if(syst_shift>0) {theErrorLow = errLowYup;   theErrorHigh = errHighYup;}
+    else             {theErrorLow = errLowYdown; theErrorHigh = errHighYdown;}
+
+    double theWeight, theError;
+    theWeight = yLow + (yHigh-yLow) / (xHigh-xLow) * (pt-xLow);
+    theError  = theErrorLow + (theErrorHigh-theErrorLow) / (xHigh-xLow) * (pt-xLow);
+
+    return  ( theWeight + (theError*syst_shift));
+  }
+  else {
+    std::cout << "category asked: " << theCategory << " was not found - which is a problem. Returning weight 1. " << std::endl;
+    return 1.;
+  }
+
+}

DiPhoEfficiencySmearer.h

@@ -0,0 +1,72 @@
+#ifndef __DIPHOEFFICIENCYSMEARER__
+#define __DIPHOEFFICIENCYSMEARER__
+
+#include "BaseSmearer.h"
+#include "EfficiencySmearer.h"
+#include <string>
+#include <map>
+#include "TFile.h"
+#include "TGraphAsymmErrors.h"
+
+class PhotonReducedInfo;
+class TRandom3;
+
+// ------------------------------------------------------------------------------------
+class DiPhoEfficiencySmearer : public BaseDiPhotonSmearer
+{
+public:
+
+  struct diPhoEfficiencySmearingParameters
+  {
+	  int n_categories;
+	  std::string categoryType;
+	  std::string parameterSetName;
+
+	  typedef std::map<std::string,float> parameterMap;
+	  typedef std::map<std::string,float>::iterator parameterMapIt;
+	  typedef std::map<std::string,float>::const_iterator parameterMapConstIt;
+
+	  // Scale offset and smearing error should be espressed as a relative value
+	  // Example: scale_offset["EB"]=1.002 , smearing_sigma["EB"]=0.01  
+
+	  std::map<std::string,float> scale_offset;
+	  std::map<std::string,float> scale_offset_error;
+
+	  std::map<std::string,float> smearing_sigma;
+	  std::map<std::string,float> smearing_sigma_error;
+
+          std::string efficiency_file;
+  };
+
+  DiPhoEfficiencySmearer(const diPhoEfficiencySmearingParameters& par) ;
+  virtual ~DiPhoEfficiencySmearer();
+
+  virtual const std::string & name() const { return name_; };
+
+  virtual bool smearDiPhoton( TLorentzVector & p4, TVector3 & selVtx, float & weight, const int & category, const int & genMassPoint, 
+			      const TVector3 & trueVtx, float syst_shift) const ;
+
+  void name(const std::string & x) { name_ = x; };
+
+  void setEffName(const std::string & x) { effName_ =x; };
+
+  bool init();
+  void passFailWeights(bool x) { passFailWeights_ = x; };
+  void doVtxEff(bool x) { doVtxEff_ = x; };
+
+  diPhoEfficiencySmearingParameters  myParameters_;
+
+ protected:
+
+  double getWeight(double pt, std::string theCategory, float syst_shift) const;
+
+  bool passFailWeights_, doVtxEff_;
+
+  std::string   name_;
+  TRandom3     *rgen_;
+  std::string   effName_;
+  TFile        *theDiPhoEfficiencyFile_; 
+  std::map<std::string,TGraphAsymmErrors*> smearing_eff_graph_;
+};
+
+#endif