LimitSetting/RooPDF_BKG.cxx at main · IMSA-CMS/LimitSetting · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
/*****************************************************************************
 * Project: RooFit                                                           *
 *                                                                           *
 * This code was autogenerated by RooClassFactory                            *
 *****************************************************************************/

// Your description goes here...


#include "Riostream.h"

#include "RooPDF_BKG.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"
#include <RooGenericPdf.h>

// Indices in bkg_params for each type of background
static const int I_TTBAR = 1;
static const int I_ZZ = 2;
static const int I_DY = 0;
static const int I_TTZ = 3;


ClassImp(RooPDF_BKG);

 RooPDF_BKG::RooPDF_BKG(const char *name, const char *title,
                        RooAbsReal& _x,
                        const std::vector<std::vector<double>>& _bkg_params,
                        const std::string _channel_name) :
   RooAbsPdf(name,title),
   x("x","x",this,_x),
   bkg_params(_bkg_params),
   channel_name(_channel_name)
 {
  // std::cout << "Starting constructor\n";
  // for (auto vec : bkg_params)
  // {
  //   for (auto num : vec)
  //   {
  //     std::cout << "Entry: " << num << '\n';
  //   }
  // }
 }


 RooPDF_BKG::RooPDF_BKG(const RooPDF_BKG& other, const char* name) :
   RooAbsPdf(other,name),
   x("x",this,other.x),
   bkg_params(other.bkg_params),
   channel_name(other.channel_name)
 {
 }

 RooGenericPdf RooPDF_BKG::ttbar()
 {
  std::vector<double> bkg_type = bkg_params[I_TTBAR];
  std::string bkg_Function = std::to_string(bkg_type[0]) + "* (@0-" + std::to_string(bkg_type[1]) + ")^" + std::to_string(bkg_type[2]);
  RooGenericPdf bkg_ttbar((channel_name + "_bkg_ttbar").c_str(), bkg_Function.c_str(), RooArgSet(*x.absArg()));
  return bkg_ttbar;
 }

 RooGenericPdf RooPDF_BKG::ZZ()
 {
  std::vector<double> bkg_type = bkg_params[I_ZZ];
  std::string bkg_Function = std::to_string(bkg_type[0]) + "* (@0-" + std::to_string(bkg_type[1]) + ")^" + std::to_string(bkg_type[2]);
  RooGenericPdf bkg_ZZ((channel_name + "_bkg_ZZ").c_str(), bkg_Function.c_str(), RooArgSet(*x.absArg()));
  return bkg_ZZ;
 }

 RooGenericPdf RooPDF_BKG::DY()
 {
  std::vector<double> bkg_type = bkg_params[I_DY];
  std::string bkg_Function = std::to_string(bkg_type[0]) + "* (@0-" + std::to_string(bkg_type[1]) + ")^" + std::to_string(bkg_type[2]);
  RooGenericPdf bkg_DY((channel_name + "_bkg_DY").c_str(), bkg_Function.c_str(), RooArgSet(*x.absArg()));
  return bkg_DY;
 }

 RooGenericPdf RooPDF_BKG::TTZ()
 {
  std::vector<double> bkg_type = bkg_params[I_TTZ];
  std::string bkg_Function = std::to_string(bkg_type[0]) + "* (@0-" + std::to_string(bkg_type[1]) + ")^" + std::to_string(bkg_type[2]);
  RooGenericPdf bkg_TTZ((channel_name + "_bkg_TTZ").c_str(), bkg_Function.c_str(), RooArgSet(*x.absArg()));
  return bkg_TTZ;
 }


 Double_t RooPDF_BKG::evaluate() const
 {
  //std::cout << "Starting evaluate()\n";
  std::vector<double> component_functions;
  double function_sum = 0;
	for (std::vector<double> bkg_type : bkg_params)
	{
    if (bkg_type.size() < 2)
    {
      continue;
    }
    //std::cerr << "Next loop iteration\n";
    //std::cerr << bkg_type[0] << '\n' << bkg_type[1] << '\n' << bkg_type[2] <<'\n';

    //std::cout << "BKG Components: " << bkg_type[0] << " " << bkg_type[1] << " " << bkg_type[2] << "\n";


	  double bkg_component = bkg_type[0] * std::pow((x - bkg_type[1]),bkg_type[2]);
    bool isNaN = std::isnan(bkg_component);
    if (isNaN)
    {
      //std::cout << "NaN value from evaluating " << bkg_type[0] << " * " << "(" << x << "-" << bkg_type[1] << ")^" << bkg_type[2] << "\n";
    }
    // else
    // {
    //   std::cout << "Is a number" << "\n";
    // }
    //bool secondInvalidValueTest = bkg_type[2] <= 0;
    //std::cout << "Test statement: " << invalidValueTest << "\n";
		component_functions.push_back(bkg_component);
    function_sum += bkg_component;
	}
  //std::cout << GetName() << " bkg function_sum: " << function_sum << "\n";
  //std::cout << "\n \n \n";
  return function_sum;
 }