From c62ced79ca83d61262dbf166d1f1a969c02d3248 Mon Sep 17 00:00:00 2001
From: hyejin kwon <hj.kwon@cern.ch>
Date: Tue, 14 Nov 2023 14:59:19 +0100
Subject: [PATCH] add DoubleSidedGausssian pdf

---
 interface/DoubleSidedGaussian.h | 60 +++++++++++++++++++++++++++++++++
 src/DoubleSidedGaussian.cc      | 59 ++++++++++++++++++++++++++++++++
 src/classes.h                   |  1 +
 src/classes_def.xml             |  1 +
 4 files changed, 121 insertions(+)
 create mode 100644 interface/DoubleSidedGaussian.h
 create mode 100644 src/DoubleSidedGaussian.cc

diff --git a/interface/DoubleSidedGaussian.h b/interface/DoubleSidedGaussian.h
new file mode 100644
index 00000000000..4e3830afcb8
--- /dev/null
+++ b/interface/DoubleSidedGaussian.h
@@ -0,0 +1,60 @@
+/*****************************************************************************
+ * Project: RooFit                                                           *
+ *                                                                           *
+  * This code was autogenerated by RooClassFactory                            * 
+ *****************************************************************************/
+
+#ifndef MY_DOUBLE_Gaussian
+#define MY_DOUBLE_Gaussian
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+#include "RooCategoryProxy.h"
+#include "RooAbsReal.h"
+#include "RooAbsCategory.h"
+#include "RooGlobalFunc.h"
+
+#include "RooArgusBG.h"
+#include "RooRealVar.h"
+#include "RooDataSet.h"
+
+#include "RooGaussian.h"
+#include "RooCBShape.h"
+#include "RooExponential.h"
+#include "RooBreitWigner.h"
+
+#include "RooConstVar.h"
+#include "RooDataHist.h"
+#include "RooFitResult.h"
+#include "RooMinuit.h"
+#include "RooPlot.h"
+ 
+class DoubleSidedGaussian : public RooAbsPdf {
+public:
+  DoubleSidedGaussian() {} ; 
+  DoubleSidedGaussian(const char *name, const char *title,
+	      RooAbsReal& _x,
+	      RooAbsReal& _mean,
+	      RooAbsReal& _sig1,
+	      RooAbsReal& _sig2,
+              Double_t _yMax);
+  DoubleSidedGaussian(const DoubleSidedGaussian& other, const char* name=0) ;
+  virtual TObject* clone(const char* newname) const { return new DoubleSidedGaussian(*this,newname); }
+  inline virtual ~DoubleSidedGaussian() { }
+
+protected:
+
+  RooRealProxy x ;
+  RooRealProxy mean ;
+  RooRealProxy sig1 ;
+  RooRealProxy sig2 ;
+  Double_t evaluate() const ;
+
+private:
+
+  Double_t yMax;
+  ClassDef(DoubleSidedGaussian,1) // Your description goes here...
+};
+ 
+#endif
+
diff --git a/src/DoubleSidedGaussian.cc b/src/DoubleSidedGaussian.cc
new file mode 100644
index 00000000000..d0e24f10274
--- /dev/null
+++ b/src/DoubleSidedGaussian.cc
@@ -0,0 +1,59 @@
+/***************************************************************************** 
+ * Project: RooFit                                                           * 
+ *                                                                           * 
+ * This code was autogenerated by RooClassFactory                            * 
+ *****************************************************************************/ 
+
+// Your description goes here... 
+
+#include "../interface/DoubleSidedGaussian.h" 
+#include "RooAbsReal.h" 
+#include "RooAbsCategory.h" 
+#include <math.h> 
+#include "TMath.h" 
+#include "RooGaussian.h"
+
+ClassImp(DoubleSidedGaussian) 
+
+DoubleSidedGaussian::DoubleSidedGaussian(const char *name, const char *title, 
+                       RooAbsReal& _x,
+                       RooAbsReal& _mean,
+                       RooAbsReal& _sig1,
+                       RooAbsReal& _sig2,
+                       Double_t _yMax) :
+  RooAbsPdf(name,title), 
+  x("x","x",this,_x),
+  mean("mean","mean",this,_mean),
+  sig1("sig1","sig1",this,_sig1),
+  sig2("sig2","sig2",this,_sig2),
+  yMax(_yMax)
+{
+} 
+
+
+DoubleSidedGaussian::DoubleSidedGaussian(const DoubleSidedGaussian& other, const char* name) :  
+  RooAbsPdf(other,name), 
+  x("x",this,other.x),
+  mean("mean",this,other.mean),
+  sig1("sig1",this,other.sig1),
+  sig2("sig2",this,other.sig2),
+  yMax(other.yMax)
+{ 
+} 
+
+
+
+Double_t DoubleSidedGaussian::evaluate() const 
+{ 
+  // From wikipedia
+  double result = -1;
+  double sqrt2OverPi = TMath::Sqrt(2. / TMath::Pi());
+  double mode = mean; // redefine mean to mode
+  double A = sqrt2OverPi / (sig2+sig1);
+  if (x < mode)
+    result = A * TMath::Exp( -1 * (x-mode)*(x-mode) / (2*sig1*sig1) );
+  else
+    result = A * TMath::Exp( -1 * (x-mode)*(x-mode) / (2*sig2*sig2) );
+  return result;
+} 
+
diff --git a/src/classes.h b/src/classes.h
index 867ffeee1c9..d334f11d428 100644
--- a/src/classes.h
+++ b/src/classes.h
@@ -69,3 +69,4 @@
 #include "HiggsAnalysis/CombinedLimit/interface/CMSExternalMorph.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSInterferenceFunc.h"
 #include "HiggsAnalysis/CombinedLimit/interface/RooEFTScalingFunction.h"
+#include "HiggsAnalysis/CombinedLimit/interface/DoubleSidedGaussian.h"
diff --git a/src/classes_def.xml b/src/classes_def.xml
index 3bf3eb4f689..e90cb94fe85 100644
--- a/src/classes_def.xml
+++ b/src/classes_def.xml
@@ -225,4 +225,5 @@
   <class name="CMSExternalMorph" />
   <class name="CMSInterferenceFunc" />
   <class name="RooEFTScalingFunction" />
+  <class name="DoubleSidedGaussian" />
 </lcgdict>