Incorporating Kinematic Fit implementation for ZGamma

README.md

@@ -37,7 +37,7 @@ you can use the below commands for the setup.
 This also sets up the environment for the batch system at UCSB
 
 ~~~~bash
-git clone --recurse-submodules [email protected]:richstu/draw_pico.git
+git clone --recursive [email protected]:richstu/draw_pico.git
 cd draw_pico
 source set_env.sh
 ./compile.py or scons

inc/zgamma/KinZfitter.h

@@ -0,0 +1,121 @@
+/*************************************************************************
+*  Authors:   Tongguang CHeng(IHEP, Beijing) Hualin Mei(UF)
+*************************************************************************/
+#ifndef KinZfitter_h
+#define KinZfitter_h
+
+// C++ includes
+#include <iostream>
+#include <complex>
+#include <cstdlib>
+#include <vector>
+#include <fstream>
+#include <sstream>
+#include <cmath>
+#include <map>
+
+// ROOT includes
+#include "TString.h"
+#include "TLorentzVector.h"
+
+// ROOFIT
+#include "RooRealVar.h"
+#include "RooArgSet.h"
+#include "RooGaussian.h"
+#include "RooBreitWigner.h"
+#include "RooProdPdf.h"
+#include "RooDataSet.h"
+#include "RooGlobalFunc.h"
+#include "RooDataHist.h"
+#include "RooHistPdf.h"
+#include "RooCBShape.h"
+#include "RooMinuit.h"
+#include "RooFormulaVar.h"
+#include "RooAddPdf.h"
+#include "RooGenericPdf.h"
+#include "RooFFTConvPdf.h"
+#include "RooFitResult.h"
+
+// fit result covariance matrix
+#include <TMatrixDSym.h>
+
+using namespace std;
+
+class KinZfitter {
+ public:
+
+  KinZfitter(bool isData);
+
+  /// Kinematic fit of lepton momenta
+  void Setup(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhotons, std::map<unsigned int, double> errorLeptons);
+
+  ///
+  void KinRefitZ1();
+
+  int  PerZ1Likelihood(double & l1, double & l2, double & lph1, double & lph2);
+  void SetZ1Result(double l1, double l2, double lph1, double lph2);
+  double pterr(TLorentzVector fsrPhoton);
+  double masserror(std::vector<TLorentzVector> p4s, std::vector<double> pTErrs);
+  double masserrorFullCov(std::vector<TLorentzVector> p4s, TMatrixDSym covMatrix);
+
+  // result wrappers
+  double GetRefitM4l();
+  double GetM4l();
+  double GetRefitMZ1();
+
+  double GetMZ1Err();
+  double GetRefitM4lErr();
+  double GetM4lErr();
+  double GetRefitM4lErrFullCov();
+
+  std::vector<TLorentzVector> GetRefitP4s();
+  std::vector<TLorentzVector> GetP4s();
+
+  ////////////////
+
+  // Need to check in deep whether this function is doable
+  // RooFormulaVar p1DOTp2(RooRealVar pT1, RooRealVar theta1, RooRealVar phi1, RooRealVar m1, TString index1, RooRealVar pT2, RooRealVar theta2, RooRealVar phi2, RooRealVar m2, TString index2);
+
+ private:
+
+  /// True mZ/mZ1 shape, final states
+  TString PDFName_, fs_;
+
+  /// debug flag
+  bool debug_;
+
+  /// whether use correction for pT error
+  bool isCorrPTerr_;
+  /// whether use data or mc correction
+  bool isData_;
+
+  void initZs(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhoton, std::map<unsigned int, double> errorLeptons);
+
+  /// lepton ids for Z1
+  std::vector<int> idsZ1_;
+  /// lepton ids that fsr photon associated to
+  std::vector<int> idsFsrZ1_;
+  /// (Four) TLorentzVectors that form the Higgs Candidate
+  std::vector<TLorentzVector> p4sZ1_, p4sZ1ph_;
+  std::vector<TLorentzVector> p4sZ1REFIT_, p4sZ1phREFIT_;
+
+  /// pTerr vector
+  std::vector<double> pTerrsZ1_, pTerrsZ1ph_;
+  std::vector<double> pTerrsZ1REFIT_, pTerrsZ1phREFIT_;
+
+  // covariance matrix
+  // what directly coming from Refit
+  TMatrixDSym covMatrixZ1_;
+
+  // refit energy scale with respect to reco pT
+  double lZ1_l1_, lZ1_l2_;
+  double lZ1_ph1_, lZ1_ph2_;
+
+  // True mZ1 shape parameters
+  // double sgVal_, aVal_, nVal_, fVal_, meanVal_, sigmaVal_, f1Val_;
+  double meanCB_, sigmaCB_, alphaCB_, nCB_, meanGauss1_, sigmaGauss1_, f1_, meanGauss2_, sigmaGauss2_, f2_, meanGauss3_, sigmaGauss3_, f3_;
+  // double aVal_, meanVal_, nVal_, sgVal_, bwsigVal_;
+
+};
+
+#endif

inc/zgamma/zg_utilities.hpp

@@ -32,5 +32,9 @@ namespace ZgUtilities {
   double cos_theta(const Baby &b, bool gen = false);
   double cos_Theta(const Baby &b, bool gen = false);
   double Getphi(const Baby &b, bool gen = false);
+  double KinRefit(const Baby &b);
+  std::vector<TLorentzVector> RefitP4(const Baby &b);
+  double AssignL1Error(const Baby &b);
+  double AssignL2Error(const Baby &b);
 }
 #endif