JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
SuperMELA.cc
Go to the documentation of this file.
1 #include <sstream>
2 #include <cmath>
3 #include <cassert>
4 #include "MELAStreamHelpers.hh"
5 #include "TVar.hh"
6 #include "SuperMELA.h"
7 #include "MELAHXSWidth.h"
8 #include "RooArgSet.h"
9 #include "RooArgList.h"
10 
11 
12 using namespace std;
13 using namespace RooFit;
14 using MELAStreamHelpers::MELAout;
15 using MELAStreamHelpers::MELAerr;
16 
17 
18 SuperMELA::SuperMELA(
19  double mH,
20  string channel,
21  double LHCsqrts
22  ) :
23  mHVal_(mH),
24  sqrts_(LHCsqrts)
25 {
26  sigma_CB_=nullptr;
27  mean_CB_err_=nullptr;
28  sigma_CB_err_=nullptr;
29 
30  n_CB_=nullptr;
31  alpha_CB_=nullptr;
32  n2_CB_=nullptr;
33  alpha2_CB_=nullptr;
34  mean_CB_=nullptr;
35  meanTOT_CB_=nullptr;
36 
37  mean_BW_=nullptr;
38  width_BW_=nullptr;
39 
40  corr_mean_sig=nullptr;
41  corr_sigma_sig=nullptr;
42 
43  sig_CB_=nullptr;
44  sig_BW_=nullptr;
45  sig_FFT_=nullptr;
46 
47  a0_qqZZ_=nullptr;
48  a1_qqZZ_=nullptr;
49  a2_qqZZ_=nullptr;
50  a3_qqZZ_=nullptr;
51  a4_qqZZ_=nullptr;
52  a5_qqZZ_=nullptr;
53  a6_qqZZ_=nullptr;
54  a7_qqZZ_=nullptr;
55  a8_qqZZ_=nullptr;
56  a9_qqZZ_=nullptr;
57  a10_qqZZ_=nullptr;
58  a11_qqZZ_=nullptr;
59  a12_qqZZ_=nullptr;
60  a13_qqZZ_=nullptr;
61 
62  qqZZ_pdf_=nullptr;
63 
64  verbose_=false;
65  mH_rrv_=new RooRealVar("mH", "mH", mHVal_, 0., sqrts_*1000.);
66  m4l_rrv_=nullptr;
67  strChan_=channel;
68 
69  //pathToCards_="../../../HiggsAnalysis/HZZ4L_CombinationPy/CreateDatacards/SM_inputs_8TeV/";
70  pathToCards_="../data/CombinationInputs/SM_inputs_8TeV/";
71  //init();
72 }
73 
74 SuperMELA::~SuperMELA(){
75  // Delete in reverse order wrt SuperMELA::init()
76  delete qqZZ_pdf_; qqZZ_pdf_=nullptr;
77 
78  delete a0_qqZZ_; delete a1_qqZZ_; delete a2_qqZZ_; delete a3_qqZZ_;
79  delete a4_qqZZ_; delete a5_qqZZ_; delete a6_qqZZ_; delete a7_qqZZ_;
80  delete a8_qqZZ_; delete a9_qqZZ_; delete a10_qqZZ_; delete a11_qqZZ_;
81  delete a12_qqZZ_; delete a13_qqZZ_;
82  a0_qqZZ_=nullptr; a1_qqZZ_=nullptr; a2_qqZZ_=nullptr; a3_qqZZ_=nullptr;
83  a4_qqZZ_=nullptr; a5_qqZZ_=nullptr; a6_qqZZ_=nullptr; a7_qqZZ_=nullptr;
84  a8_qqZZ_=nullptr; a9_qqZZ_=nullptr; a10_qqZZ_=nullptr; a11_qqZZ_=nullptr;
85  a12_qqZZ_=nullptr; a13_qqZZ_=nullptr;
86 
87  delete sig_FFT_; sig_FFT_=nullptr;
88  delete sig_BW_; sig_BW_=nullptr;
89  delete sig_CB_; sig_CB_=nullptr;
90 
91  delete width_BW_; width_BW_=nullptr;
92  delete mean_BW_; mean_BW_=nullptr;
93 
94  delete sigma_CB_; sigma_CB_=nullptr;
95  delete meanTOT_CB_; meanTOT_CB_=nullptr;
96  delete mean_CB_; mean_CB_=nullptr;
97 
98  delete corr_sigma_sig; corr_sigma_sig=nullptr;
99  delete corr_mean_sig; corr_mean_sig=nullptr;
100 
101  delete alpha2_CB_; alpha2_CB_=nullptr;
102  delete n2_CB_; n2_CB_=nullptr;
103  delete alpha_CB_; alpha_CB_=nullptr;
104  delete n_CB_; n_CB_=nullptr;
105 
106  delete sigma_CB_err_; sigma_CB_err_=nullptr;
107  delete mean_CB_err_; mean_CB_err_=nullptr;
108 
109  delete m4l_rrv_; m4l_rrv_=nullptr;
110  delete mH_rrv_; mH_rrv_=nullptr;
111 }
112 
113 void SuperMELA::SetDecayChannel(string myChan){
114  if (verbose_) MELAout << "SuperMELA::SetDecayChannel: Switching from " << strChan_ << " to " << myChan << endl;
115  if (myChan != strChan_){ // do nothing if it's the same as before
116  strChan_=myChan;
117  bool newChanOK = checkChannel();
118  if (verbose_) MELAout << "Setting decay channel of SuperMELA to " << strChan_.c_str() << " , re-initializing..." << endl;
119  init();
120  if (verbose_ && newChanOK) MELAout << "Decay channel set successfully to " << strChan_.c_str() << endl;
121  }
122 }
123 
124 double SuperMELA::GetSigShapeSystematic(string parName){
125  if (parName=="meanCB") return mean_CB_err_->getVal();
126  else if (parName=="sigmaCB") return sigma_CB_err_->getVal();
127  else{
128  MELAout << "Error from SuperMELA::GetSigShapeSystematic, unrecognized input: " << parName.c_str() << endl;
129  try{
130  throw 40;
131  }
132  catch (int e){
133  if (e==40) MELAout << "Exception " << e << " in SuperMELA::GetSigShapeSystematic! Unrecognized type of parameter requested: " << parName.c_str() << " . Valid options are meanCB ; sigmaCB" << endl;
134  return 0;
135  }
136  }
137  return 0;
138 }
139 
140 double SuperMELA::GetSigShapeParameter(string parName){
141  if (parName=="meanCB") return meanTOT_CB_->getVal();
142  else if (parName=="sigmaCB") return sigma_CB_->getVal();
143  else if (parName=="alphaCB") return alpha_CB_->getVal();
144  else if (parName=="nCB") return n_CB_->getVal();
145  else{
146  try{
147  throw 40;
148  }
149  catch (int e){
150  if (e==40) MELAout << "Exception " << e << " in SuperMELA::GetSigShapeParameter! Unrecognized type of parameter requested: " << parName.c_str() << " . Valid options are meanCB ; sigmaCB ; alphaCB ; nCB" << endl;
151  return 0;
152  }
153  }
154  return 0;
155 }
156 
157 
158 void SuperMELA::init(){
159  if (verbose_) MELAout << "Begin SuperMELA::init..." << endl;
160 
161  // Calculate m4l ranges for the given mH, set range of rrv
162  calc_mZZ_range(mHVal_, lowMH_, highMH_);
163  if (verbose_) MELAout << "Range width=" << highMH_ - lowMH_ << endl;
164  delete m4l_rrv_; m4l_rrv_=new RooRealVar("CMS_zz4l_mass", "CMS_zz4l_mass", mHVal_, lowMH_, highMH_);
165  m4l_rrv_->setBins(2000, "fft");
166  m4l_rrv_->setRange("shape", lowMH_, highMH_);
167 
168  //parameters for signal m4l shape systematics
169  double str_mean_CB_err_e;
170  double str_mean_CB_err_m;
171  double str_sigma_CB_err_e;
172  double str_sigma_CB_err_m;
173  readSigSystFromFile(str_mean_CB_err_e, str_mean_CB_err_m, str_sigma_CB_err_e, str_sigma_CB_err_m);
174  if (verbose_){
175  MELAout << "mean_CB systematics (ele and mu): " << str_mean_CB_err_e << " / " << str_mean_CB_err_m << endl;
176  MELAout << "sigma_CB systematics (ele and mu): " << str_sigma_CB_err_e << " / " << str_sigma_CB_err_m << endl;
177  }
178 
179  //delete old stuff before reinitialization
180  delete mean_CB_err_;
181  delete sigma_CB_err_;
182  if (strChan_=="4mu"){
183  mean_CB_err_=new RooRealVar("mean_CB_err", "mean_CB_err", str_mean_CB_err_m);
184  sigma_CB_err_=new RooRealVar("sigma_CB_err", "sigma_CB_err", str_sigma_CB_err_m);
185  }
186  else if (strChan_=="4e"){
187  mean_CB_err_=new RooRealVar("mean_CB_err", "mean_CB_err", str_mean_CB_err_e);
188  sigma_CB_err_=new RooRealVar("sigma_CB_err", "sigma_CB_err", str_sigma_CB_err_e);
189  }
190  else{//2e2mu, we should have already checked that the string of the channel is a sensible one
191  mean_CB_err_=new RooRealVar("mean_CB_err", "mean_CB_err", (str_mean_CB_err_m + str_mean_CB_err_e));
192  sigma_CB_err_=new RooRealVar("sigma_CB_err", "sigma_CB_err", sqrt(pow(str_sigma_CB_err_m, 2) + pow(str_sigma_CB_err_e, 2)));
193  }
194  if (mean_CB_err_->getVal()<0.){ MELAout << "Negative error on the m4l mean ! " << mean_CB_err_->getVal() << endl; }
195  if (sigma_CB_err_->getVal()<0.){ MELAout << "Negative error on the m4l sigma ! " << sigma_CB_err_->getVal() << endl; }
196 
197 
198  //set parameters for signal m4l shape and calculate normalization
199  string str_n_CB;
200  string str_alpha_CB;
201  string str_n2_CB;
202  string str_alpha2_CB;
203  string str_mean_CB;
204  string str_sigma_CB;
205  if (verbose_)MELAout << "Reading signal shape formulas" << endl;
206  readSigParsFromFile(str_mean_CB, str_sigma_CB, str_n_CB, str_alpha_CB, str_n2_CB, str_alpha2_CB);
207  if (verbose_){
208  MELAout << "Read from input card the following formulas: " << endl;
209  MELAout << "Mean RooFormula (string): " << str_mean_CB.c_str() << endl;
210  MELAout << "Sigma RooFormula (string): " << str_sigma_CB.c_str() << endl;
211  }
212 
213  delete n_CB_;
214  delete alpha_CB_;
215  delete n2_CB_;
216  delete alpha2_CB_;
217  char rrvName[96];
218  sprintf(rrvName, "CMS_zz4l_n_sig_%s_%d", strChan_.c_str(), int(sqrts_));
219  if (verbose_) MELAout << "SuperMELA::init: Constructing n_CB_ from formula " << str_n_CB.c_str() << endl;
220  n_CB_=new RooFormulaVar(rrvName, str_n_CB.c_str(), RooArgList(*mH_rrv_));
221  sprintf(rrvName, "CMS_zz4l_alpha_sig_%s_%d", strChan_.c_str(), int(sqrts_));
222  if (verbose_) MELAout << "SuperMELA::init: Constructing alpha_CB_ from formula " << str_alpha_CB.c_str() << endl;
223  alpha_CB_=new RooFormulaVar(rrvName, str_alpha_CB.c_str(), RooArgList(*mH_rrv_));
224  sprintf(rrvName, "CMS_zz4l_n2_sig_%s_%d", strChan_.c_str(), int(sqrts_));
225  if (verbose_) MELAout << "SuperMELA::init: Constructing n2_CB_ from formula " << str_n2_CB.c_str() << endl;
226  n2_CB_=new RooFormulaVar(rrvName, str_n2_CB.c_str(), RooArgList(*mH_rrv_));
227  sprintf(rrvName, "CMS_zz4l_alpha2_sig_%s_%d", strChan_.c_str(), int(sqrts_));
228  if (verbose_) MELAout << "SuperMELA::init: Constructing alpha2_CB_ from formula " << str_alpha2_CB.c_str() << endl;
229  alpha2_CB_=new RooFormulaVar(rrvName, str_alpha2_CB.c_str(), RooArgList(*mH_rrv_));
230 
231  corr_mean_sig = new RooRealVar("CMS_zz4l_mean_sig_corrMH", "CMS_zz4l_mean_sig_corrMH", 0., -10., 10.);
232  corr_sigma_sig = new RooRealVar("CMS_zz4l_sigma_sig_corrMH", "CMS_zz4l_sigma_sig_corrMH", 0., -10., 10.);
233  corr_mean_sig->setConstant(true);
234  corr_sigma_sig->setConstant(true);
235 
236  delete mean_CB_;
237  delete meanTOT_CB_;
238  delete sigma_CB_;
239  mean_CB_=new RooFormulaVar("CMS_zz4l_mean_m_sig", Form("(%s)+@0*@1", str_mean_CB.c_str()), RooArgList(*mH_rrv_, *corr_mean_sig));//this is normalized by mHVal_
240  meanTOT_CB_=new RooFormulaVar("CMS_zz4l_mean_sig", "(@0+@1)", RooArgList(*mH_rrv_, *mean_CB_));
241  if (verbose_){ MELAout << "Signal Mean vals -> Correction: " << corr_mean_sig->getVal() << " Mean: " << mean_CB_->getVal() << " Total: " << meanTOT_CB_->getVal() << endl; }
242  sigma_CB_=new RooFormulaVar("CMS_zz4l_sigma_m_sig", Form("(%s)*(1+@0*@1)", str_sigma_CB.c_str()), RooArgList(*mH_rrv_, *corr_sigma_sig));
243 
244  //for high-mass one also needs a gamma RooFormulaVar,
245  delete mean_BW_;
246  delete width_BW_;
247  sprintf(rrvName, "CMS_zz4l_mean_BW_sig_%s_%d", strChan_.c_str(), int(sqrts_));
248  mean_BW_=new RooRealVar(rrvName, "CMS_zz4l_mean_BW", mHVal_, 100., 1000.);
249  sprintf(rrvName, "CMS_zz4l_width_BW_sig_%s_%d", strChan_.c_str(), int(sqrts_));
250  width_BW_=new RooRealVar(rrvName, "CMS_zz4l_width_BW", 1.);
251  mean_BW_->setVal(mHVal_);
252  mean_BW_->setConstant(true);
253  width_BW_->setConstant(true);
254 
255  if (verbose_){
256  MELAout << "Signal shape parameter values: " << endl;
257  MELAout << "Mean (formula value) = " << meanTOT_CB_->getVal() << endl;
258  MELAout << "Sigma (formula value) = " << sigma_CB_->getVal() << endl;
259  MELAout << "n (formula value) = " << n_CB_->getVal() << endl;
260  MELAout << "alpha (formula value) = " << alpha_CB_->getVal() << endl;
261  MELAout << "n2 (formula value) = " << n2_CB_->getVal() << endl;
262  MELAout << "alpha2 (formula value) = " << alpha2_CB_->getVal() << endl;
263  MELAout << "Mean BW (realvar value) = " << mean_BW_->getVal() << endl;
264  MELAout << "Width BW (realvar value) = " << width_BW_->getVal() << endl;
265  }
266 
267  delete sig_FFT_;
268  delete sig_CB_;
269  delete sig_BW_;
270  sig_CB_ =new MELADoubleCB("signalCB_ggH", "signalCB_ggH", *m4l_rrv_, *meanTOT_CB_, *sigma_CB_, *alpha_CB_, *n_CB_, *alpha2_CB_, *n2_CB_);
271  sig_BW_ =new MELARelBWUFParam("signalBW_ggH", "signalBW_ggH", *m4l_rrv_, *mean_BW_, *width_BW_);
272  //sig_FFT_=new RooFFTConvPdf("signal_ggH", "BW (X) CB", *m4l_rrv_, *sig_BW_, *sig_CB_, 2);
273  //sig_FFT_->setBufferFraction(0.2);
274  if (verbose_){
275  MELAout << "Value of signal m4l CB shape is " << sig_CB_->getVal() << endl;
276  MELAout << "Value of signal m4l BW shape is " << sig_BW_->getVal() << endl;
277  //sig_FFT_->Print("v");
278  //MELAout << "Value of signal m4l FFT shape is " << sig_FFT_->getVal() << endl;
279  }
280 
281  RooAbsReal* tmpint;
282 
283  tmpint = sig_CB_->createIntegral(RooArgSet(*m4l_rrv_), RooFit::Range("shape"));
284  norm_sig_CB_ =tmpint->getVal();
285  delete tmpint;
286  if (verbose_)MELAout << "Normalization of signal m4l CB shape is " << norm_sig_CB_ << endl;
287 
288  if (verbose_)MELAout << "\n---> Integrating Breit-Wigner:" << endl;
289  tmpint = sig_BW_->createIntegral(RooArgSet(*m4l_rrv_), RooFit::Range("shape"));
290  double norm_sig_BW_ =tmpint->getVal();
291  delete tmpint;
292  if (verbose_)MELAout << "Normalization of signal m4l BW shape is " << norm_sig_BW_ << endl;
293 
294  /*
295  if (verbose_)MELAout << "\n---> Integrating full signal:" << endl;
296  tmpint = sig_FFT_->createIntegral(RooArgSet(*m4l_rrv_), RooFit::Range("shape"));
297  norm_sig_FFT_=tmpint->getVal();
298  delete tmpint;
299  if (verbose_)MELAout << "Normalization of signal m4l shape is " << norm_sig_FFT_ << endl;
300  */
301 
302  if (verbose_)MELAout << "Reading background shape parameters" << endl;
303  vector<double> v_apars;
304  readBkgParsFromFile(v_apars);
305 
306  if (verbose_){
307  MELAout << "Size of vector with bkg shape pars is " << v_apars.size() << endl;
308  MELAout << "Param [0]=" << v_apars.at(0) << " [13]=" << v_apars.at(13) << endl;
309  }
310 
311  delete a0_qqZZ_; a0_qqZZ_=new RooRealVar("CMS_zz4l_a0_qqZZ", "CMS_zz4l_a0_qqZZ", v_apars.at(0), 0., 200.);
312  delete a1_qqZZ_; a1_qqZZ_=new RooRealVar("CMS_zz4l_a1_qqZZ", "CMS_zz4l_a1_qqZZ", v_apars.at(1), 0., 200.);
313  delete a2_qqZZ_; a2_qqZZ_=new RooRealVar("CMS_zz4l_a2_qqZZ", "CMS_zz4l_a2_qqZZ", v_apars.at(2), 0., 200.);
314  delete a3_qqZZ_; a3_qqZZ_=new RooRealVar("CMS_zz4l_a3_qqZZ", "CMS_zz4l_a3_qqZZ", v_apars.at(3), 0., 1.);
315  delete a4_qqZZ_; a4_qqZZ_=new RooRealVar("CMS_zz4l_a4_qqZZ", "CMS_zz4l_a4_qqZZ", v_apars.at(4), 0., 200.);
316  delete a5_qqZZ_; a5_qqZZ_=new RooRealVar("CMS_zz4l_a5_qqZZ", "CMS_zz4l_a5_qqZZ", v_apars.at(5), 0., 200.);
317  delete a6_qqZZ_; a6_qqZZ_=new RooRealVar("CMS_zz4l_a6_qqZZ", "CMS_zz4l_a6_qqZZ", v_apars.at(6), 0., 100.);
318  delete a7_qqZZ_; a7_qqZZ_=new RooRealVar("CMS_zz4l_a7_qqZZ", "CMS_zz4l_a7_qqZZ", v_apars.at(7), 0., 1.);
319  delete a8_qqZZ_; a8_qqZZ_=new RooRealVar("CMS_zz4l_a8_qqZZ", "CMS_zz4l_a8_qqZZ", v_apars.at(8), 0., 200.);
320  delete a9_qqZZ_; a9_qqZZ_=new RooRealVar("CMS_zz4l_a9_qqZZ", "CMS_zz4l_a9_qqZZ", v_apars.at(9), 0., 1.);
321  delete a10_qqZZ_; a10_qqZZ_=new RooRealVar("CMS_zz4l_a10_qqZZ", "CMS_zz4l_a10_qqZZ", v_apars.at(10), 0., 200.);
322  delete a11_qqZZ_; a11_qqZZ_=new RooRealVar("CMS_zz4l_a11_qqZZ", "CMS_zz4l_a11_qqZZ", v_apars.at(11), -100., 100.);
323  delete a12_qqZZ_; a12_qqZZ_=new RooRealVar("CMS_zz4l_a12_qqZZ", "CMS_zz4l_a12_qqZZ", v_apars.at(12), 0., 10000.);
324  delete a13_qqZZ_; a13_qqZZ_=new RooRealVar("CMS_zz4l_a13_qqZZ", "CMS_zz4l_a13_qqZZ", v_apars.at(13), 0., 1.);
325  a0_qqZZ_->setConstant(kTRUE);
326  a1_qqZZ_->setConstant(kTRUE);
327  a2_qqZZ_->setConstant(kTRUE);
328  a3_qqZZ_->setConstant(kTRUE);
329  a4_qqZZ_->setConstant(kTRUE);
330  a5_qqZZ_->setConstant(kTRUE);
331  a6_qqZZ_->setConstant(kTRUE);
332  a7_qqZZ_->setConstant(kTRUE);
333  a8_qqZZ_->setConstant(kTRUE);
334  a9_qqZZ_->setConstant(kTRUE);
335  a10_qqZZ_->setConstant(kTRUE);
336  a11_qqZZ_->setConstant(kTRUE);
337  a12_qqZZ_->setConstant(kTRUE);
338  a13_qqZZ_->setConstant(kTRUE);
339 
340 
341  delete qqZZ_pdf_; qqZZ_pdf_ = new MELAqqZZPdf_v2("bkg_qqzz", "bkg_qqzz", *m4l_rrv_, *a0_qqZZ_, *a1_qqZZ_, *a2_qqZZ_, *a3_qqZZ_,
342  *a4_qqZZ_, *a5_qqZZ_, *a6_qqZZ_, *a7_qqZZ_,
343  *a8_qqZZ_, *a9_qqZZ_, *a10_qqZZ_, *a11_qqZZ_,
344  *a12_qqZZ_, *a13_qqZZ_);
345 
346  tmpint = qqZZ_pdf_->createIntegral(RooArgSet(*m4l_rrv_), RooFit::Range("shape"));
347  norm_bkg_qqZZ_=tmpint->getVal();
348  delete tmpint;
349 }
350 
351 void SuperMELA::SetMH(double myMH){
352  mHVal_=myMH;
353  mH_rrv_->setVal(mHVal_);
354  if (verbose_) MELAout << "Setting MH to " << mHVal_ << endl;
355  init();
356 }
357 
358 void SuperMELA::SetPathToCards(string dirToCards){
359  pathToCards_=dirToCards;
360  if (verbose_) MELAout << "New path to cards is " << pathToCards_ << endl;
361 }
362 
363 void SuperMELA::splitLine(const string& rawoption, vector<string>& splitoptions, char delimiter){
364  string suboption=rawoption, result=rawoption;
365  string remnant;
366  while (result!=""){
367  size_t posEq = suboption.find(delimiter);
368  if (posEq!=string::npos){
369  result=suboption;
370  remnant=suboption.substr(posEq+1);
371  result.erase(result.begin()+posEq, result.end());
372  }
373  else{
374  result="";
375  remnant=suboption;
376  }
377  if (result!="") splitoptions.push_back(result);
378  suboption = remnant;
379  }
380  if (remnant!="") splitoptions.push_back(remnant);
381 }
382 
383 void SuperMELA::readSigSystFromFile(
384  double& str_mean_CB_err_e,
385  double& str_mean_CB_err_m,
386  double& str_sigma_CB_err_e,
387  double& str_sigma_CB_err_m
388  ){
389  bool mean_e_OK=false, sigma_e_OK=false;
390  bool mean_m_OK=false, sigma_m_OK=false;
391 
392  //open text file
393  string fCardName=pathToCards_+"inputs_"+strChan_+".txt";
394  if (verbose_) MELAout << "SuperMELA::readSigSystFromFile: Parsing signal shape systematics from input card " << fCardName.c_str() << endl;
395  ifstream card(fCardName.c_str(), ios::in);
396  if (!card.good()){
397  MELAerr << "SuperMELA::readSigSystFromFile: Input card " << fCardName << " is not good!" << endl;
398  assert(0);
399  }
400  else if (verbose_) MELAout << "SuperMELA::readSigSystFromFile: Input card " << fCardName << " can be read!" << endl;
401 
402  string line;
403  while (card.good()){
404  getline(card, line);
405  vector<string> fields;
406  splitLine(line, fields, ' ');
407 
408  if (fields.size()<2 || !(fields[0]=="systematic"&&fields[1]=="param")) continue;
409 
410  if (fields.size()!=4){
411  MELAerr << "Error in SuperMELA::readSigSystFromFile! Incorrect format of line " << line.c_str() << endl;
412  break;
413  }
414 
415  if (fields[2]=="CMS_zz4l_mean_m_sig"){
416  str_mean_CB_err_m = atof(fields.at(3).c_str());
417  mean_m_OK=true;
418  }
419  if (fields[2]=="CMS_zz4l_mean_e_sig"){
420  str_mean_CB_err_e = atof(fields.at(3).c_str());
421  mean_e_OK=true;
422  }
423  if (fields[2]=="CMS_zz4l_sigma_m_sig"){
424  str_sigma_CB_err_m = atof(fields.at(3).c_str());
425  sigma_m_OK=true;
426  }
427  if (fields[2]=="CMS_zz4l_sigma_e_sig"){
428  str_sigma_CB_err_e = atof(fields.at(3).c_str());
429  sigma_e_OK=true;
430  }
431 
432  if (mean_e_OK && sigma_e_OK && mean_m_OK && sigma_m_OK) break;
433  }
434 
435  try{
436  if ((!(mean_e_OK&&sigma_e_OK))&&(!(mean_m_OK&&sigma_m_OK))) throw 20;
437  }
438  catch (int e){
439  MELAout
440  << "Exception " << e << " in SuperMELA::readSigSystFromFile! Not all signal shape formulas were read " << mean_e_OK << " " << sigma_e_OK << " " << mean_m_OK << " " << sigma_m_OK
441  << endl;
442  }
443 
444  card.close();
445 }
446 
447 void SuperMELA::readSigParsFromFile(
448  string& str_mean_CB,
449  string& str_sigma_CB,
450  string& str_n_CB,
451  string& str_alpha_CB,
452  string& str_n2_CB,
453  string& str_alpha2_CB
454  ){
455  bool meanOK=false, sigmaOK=false, nOK=false, alphaOK=false, n2OK=false, alpha2OK=false;
456  //open text file
457  string fCardName=pathToCards_+"inputs_"+strChan_+".txt";
458  if (verbose_) MELAout << "Parsing input card " << fCardName.c_str() << endl;
459  ifstream card(fCardName.c_str(), ios::in);
460  string line;
461  while (card.good()){
462  getline(card, line);
463  vector<string> fields;
464  splitLine(line, fields, ' ');
465  if (fields.size()==0 || fields[0]!="signalShape")continue;
466  //ok, we found somethign interesting
467  if (fields.size()<3){
468  MELAout << "Error in SuperMELA::readSigParsFromFile! Incorrect format of line " << line.c_str() << " There should be at least three fields, I could find only " << fields.size() << " fields." << endl;
469  break;
470  }
471  if (fields.at(1)=="n_CB"){ str_n_CB=fields.at(2); nOK=true; }
472  if (fields.at(1)=="alpha_CB"){ str_alpha_CB=fields.at(2); alphaOK=true; }
473  if (fields.at(1)=="n2_CB"){ str_n2_CB=fields.at(2); n2OK=true; }
474  if (fields.at(1)=="alpha2_CB"){ str_alpha2_CB=fields.at(2); alpha2OK=true; }
475  if (fields.at(1)=="mean_CB"){ str_mean_CB=fields.at(2); meanOK=true; }
476  if (fields.at(1)=="sigma_CB"){ str_sigma_CB=fields.at(2); sigmaOK=true; }
477  if (verbose_) MELAout << fields.at(1) << " == " << fields.at(2) << endl;
478 
479  if (meanOK && sigmaOK && alphaOK && nOK && alpha2OK && n2OK)break;
480  }//end while loop on lines
481 
482  if (!(meanOK && sigmaOK && alphaOK && nOK && alpha2OK && n2OK)){
483  MELAout << "Exception in SuperMELA::readSigParsFromFile! Not all signal shape formulas were read " << meanOK << " " << sigmaOK << " " << alphaOK << " " << nOK << " " << alpha2OK << " " << n2OK << endl;
484  throw 20;
485  }
486 
487  card.close();
488 }
489 
490 void SuperMELA::readBkgParsFromFile(vector<double>& apars){
491  const int nPars=14;
492  int nFound=0;
493  apars.resize(14);
494  string fCardName=pathToCards_+"inputs_"+strChan_+".txt";
495  if (verbose_)MELAout << "Parsing input card " << fCardName.c_str() << endl;
496  ifstream card(fCardName.c_str(), ios::in);
497  string line;
498 
499  while (card.good()){
500  getline(card, line);
501  vector<string> fields;
502  splitLine(line, fields, ' ');
503  if (fields.size()==0 || fields[0]!="qqZZshape")continue;
504 
505  if (fields.size()<3){
506  MELAout << "Error in SuperMELA::readSigParsFromFile! Incorrect format of line \'" << line.c_str() << "\' . It contains " << fields.size() << "fields (it should be 3)" << endl;
507  break;
508  }
509 
510  stringstream ssip;
511  ssip << nFound;
512  if (fields[1]=="a"+ssip.str()+"_bkgd"){
513  apars[nFound]=atof(fields[2].c_str());
514  nFound++;
515  }
516  }
517 
518  try{
519  if (nFound!=nPars){
520  throw 30;
521  }
522  }
523  catch (int e){
524  if (e==30){
525  MELAerr << "Exception from void SuperMELA::readBkgParsFromFile(vector<double> apars ). Mismatched number of params of qqZZ shape read from card file " << fCardName.c_str() << " ---> " << nFound << " (it should be " << nPars << endl;
526  for (unsigned int j=0; j<apars.size(); j++){ MELAerr << apars[j] << " " << flush; }
527  MELAerr << endl;
528  }
529  }
530  card.close();
531 }
532 
533 void SuperMELA::calc_mZZ_range(const double mHVal, double& low_M, double& high_M){
534  //low_M=0.;
535  //high_M=sqrts_*1000.;
536 
537  const string MELAPKGPATH = TVar::GetMELAPath();
538  string path = MELAPKGPATH + "data/HiggsTotalWidth_YR3.txt";
539 
540  path.erase((find(path.rbegin(), path.rend(), '/').base()), path.end());
541  MELAHXSWidth myCSW(path.c_str());
542  double widthHVal = myCSW.HiggsWidth(mHVal);
543  double windowVal = max(widthHVal, 1.);
544  double lowside = 100.;
545  if (mHVal >= 275){ lowside = 180.; }
546  else { lowside = 100.; }
547  // Apply rounding
548  low_M = int(max((mHVal - 20.*windowVal), lowside)+0.5);
549  high_M = int(min((mHVal + 15.*windowVal), sqrts_*1000.)+0.5);
550 }
551 
552 std::pair<double, double> SuperMELA::M4lProb(double m4l){
553  double Psig =-1.;
554  double Pbkg =-1.;
555  if (m4l<lowMH_ || m4l>highMH_){
556  if (verbose_) MELAerr << "SuperMELA::M4lProb: m4l = " << m4l << " outside range [" << lowMH_ << ", " << highMH_ << "]. Setting SuperMELA to dummy values." << endl;
557  }
558  else{
559  m4l_rrv_->setVal(m4l);
560  Psig=sig_CB_->getVal() / norm_sig_CB_;
561  Pbkg=qqZZ_pdf_->getVal() / norm_bkg_qqZZ_;
562  if (verbose_){
563  MELAout << "SuperMELA::M4lProb: m4l = " << m4l
564  << ", m4lPsig = " << Psig
565  << ", m4lPbkg = " << Pbkg
566  << endl;
567  }
568  }
569  return make_pair(Psig, Pbkg);
570 }
571 
572 
573 std::pair<double, double> SuperMELA::M4lProb(std::pair<double, double> const& m4lPair){
574  double Psig =-1.;
575  double Pbkg =-1.;
576  if (
577  (m4lPair.first<lowMH_ || m4lPair.first>highMH_)
578  ||
579  (m4lPair.second<lowMH_ || m4lPair.second>highMH_)
580  ){
581  if (verbose_) MELAerr << "SuperMELA::M4lProb: m4l pairs " << m4lPair.first << " - " << m4lPair.second << " outside range [" << lowMH_ << ", " << highMH_ << "]. Setting SuperMELA to dummy values." << endl;
582  }
583  else{
584  m4l_rrv_->setVal(m4lPair.first);
585  Psig=sig_CB_->getVal() / norm_sig_CB_;
586  m4l_rrv_->setVal(m4lPair.second);
587  Pbkg=qqZZ_pdf_->getVal() / norm_bkg_qqZZ_;
588  if (verbose_){
589  MELAout << "SuperMELA::M4lProb: "
590  << "m4lPsig(" << m4lPair.first << ") = " << Psig
591  << ", m4lPbkg(" << m4lPair.second << ") = " << Pbkg
592  << endl;
593  }
594  }
595  return make_pair(Psig, Pbkg);
596 }
597 
598 bool SuperMELA::checkChannel(){
599  try{
600  if (strChan_!="4mu" && strChan_!="4e" && strChan_!="2e2mu" && strChan_!="2mu2e"){
601  throw 10;
602  }
603  }
604  catch (int e){
605  MELAerr << "Exception " << e << " from SuperMELA::SetDecayChannel(string myChan). Unrecognized string for decay channel: " << strChan_.c_str() << endl;
606  return false;
607  }
608 
609  if (strChan_=="4mu") ch_=0;
610  else if (strChan_=="4e") ch_=1;
611  else ch_=2;
612 
613  if (strChan_=="2mu2e") strChan_="2e2mu";
614  return true;
615 }
616 
SuperMELA::verbose_
bool verbose_
Definition: SuperMELA.h:60
SuperMELA::sigma_CB_err_
RooRealVar * sigma_CB_err_
Definition: SuperMELA.h:75
SuperMELA::sqrts_
double sqrts_
Definition: SuperMELA.h:56
SuperMELA::a13_qqZZ_
RooRealVar * a13_qqZZ_
Definition: SuperMELA.h:101
SuperMELA::SuperMELA
SuperMELA(double mH=125, std::string channel="4mu", double LHCsqrts=8.)
Definition: SuperMELA.cc:18
SuperMELA::a11_qqZZ_
RooRealVar * a11_qqZZ_
Definition: SuperMELA.h:99
SuperMELA::readBkgParsFromFile
void readBkgParsFromFile(std::vector< double > &apars)
Definition: SuperMELA.cc:490
SuperMELA::~SuperMELA
~SuperMELA()
Definition: SuperMELA.cc:74
MELARelBWUFParam
Definition: MELACombinePdfs.h:134
SuperMELA::meanTOT_CB_
RooFormulaVar * meanTOT_CB_
Definition: SuperMELA.h:72
SuperMELA::norm_bkg_qqZZ_
double norm_bkg_qqZZ_
Definition: SuperMELA.h:103
MELAHXSWidth
Definition: MELAHXSWidth.h:10
SuperMELA::a6_qqZZ_
RooRealVar * a6_qqZZ_
Definition: SuperMELA.h:94
SuperMELA::mean_CB_
RooFormulaVar * mean_CB_
Definition: SuperMELA.h:70
SuperMELA::SetMH
void SetMH(double myMH)
Definition: SuperMELA.cc:351
SuperMELA::readSigParsFromFile
void readSigParsFromFile(std::string &str_mean_CB, std::string &str_sigma_CB, std::string &str_n_CB, std::string &str_alpha_CB, std::string &str_n2_CB, std::string &str_alpha2_CB)
Definition: SuperMELA.cc:447
SuperMELA::a8_qqZZ_
RooRealVar * a8_qqZZ_
Definition: SuperMELA.h:96
SuperMELA::n_CB_
RooFormulaVar * n_CB_
Definition: SuperMELA.h:66
SuperMELA::highMH_
double highMH_
Definition: SuperMELA.h:57
SuperMELA::readSigSystFromFile
void readSigSystFromFile(double &str_mean_CB_err_e, double &str_mean_CB_err_m, double &str_sigma_CB_err_e, double &str_sigma_CB_err_m)
Definition: SuperMELA.cc:383
SuperMELA::sigma_CB_
RooFormulaVar * sigma_CB_
Definition: SuperMELA.h:71
SuperMELA::init
void init()
Definition: SuperMELA.cc:158
SuperMELA::qqZZ_pdf_
MELAqqZZPdf_v2 * qqZZ_pdf_
Definition: SuperMELA.h:102
SuperMELA::alpha2_CB_
RooFormulaVar * alpha2_CB_
Definition: SuperMELA.h:69
MELAOutputStreamer::close
void close()
Definition: MELAOutputStreamer.cc:145
SuperMELA::norm_sig_CB_
double norm_sig_CB_
Definition: SuperMELA.h:85
SuperMELA::a0_qqZZ_
RooRealVar * a0_qqZZ_
Definition: SuperMELA.h:88
MELAStreamHelpers::MELAout
MELAOutputStreamer MELAout
MELAHXSWidth.h
SuperMELA::SetPathToCards
void SetPathToCards(std::string dirToCards)
Definition: SuperMELA.cc:358
SuperMELA::a3_qqZZ_
RooRealVar * a3_qqZZ_
Definition: SuperMELA.h:91
SuperMELA::M4lProb
std::pair< double, double > M4lProb(double m4l)
Definition: SuperMELA.cc:552
SuperMELA::a4_qqZZ_
RooRealVar * a4_qqZZ_
Definition: SuperMELA.h:92
testME_all.int
int
Definition: testME_all.py:13
MELAHXSWidth::HiggsWidth
double HiggsWidth(double mH) const
Definition: MELAHXSWidth.cc:72
SuperMELA::splitLine
void splitLine(const std::string &rawoption, std::vector< std::string > &splitoptions, char delimiter)
Definition: SuperMELA.cc:363
SuperMELA::sig_FFT_
RooFFTConvPdf * sig_FFT_
Definition: SuperMELA.h:82
TVar::GetMELAPath
std::string GetMELAPath()
Definition: TVar.cc:121
SuperMELA::a5_qqZZ_
RooRealVar * a5_qqZZ_
Definition: SuperMELA.h:93
SuperMELA::a9_qqZZ_
RooRealVar * a9_qqZZ_
Definition: SuperMELA.h:97
SuperMELA::GetSigShapeSystematic
double GetSigShapeSystematic(std::string parName)
Definition: SuperMELA.cc:124
SuperMELA::mean_CB_err_
RooRealVar * mean_CB_err_
Definition: SuperMELA.h:74
SuperMELA::corr_mean_sig
RooRealVar * corr_mean_sig
Definition: SuperMELA.h:77
SuperMELA::a10_qqZZ_
RooRealVar * a10_qqZZ_
Definition: SuperMELA.h:98
SuperMELA::mH_rrv_
RooRealVar * mH_rrv_
Definition: SuperMELA.h:64
SuperMELA::GetSigShapeParameter
double GetSigShapeParameter(std::string parName)
Definition: SuperMELA.cc:140
SuperMELA::m4l_rrv_
RooRealVar * m4l_rrv_
Definition: SuperMELA.h:63
SuperMELA::mHVal_
double mHVal_
Definition: SuperMELA.h:55
SuperMELA::ch_
int ch_
Definition: SuperMELA.h:59
SuperMELA::calc_mZZ_range
void calc_mZZ_range(const double mHVal, double &low_M, double &high_M)
Definition: SuperMELA.cc:533
SuperMELA::lowMH_
double lowMH_
Definition: SuperMELA.h:57
SuperMELA::a7_qqZZ_
RooRealVar * a7_qqZZ_
Definition: SuperMELA.h:95
MELAStreamHelpers::MELAerr
MELAOutputStreamer MELAerr
SuperMELA::a1_qqZZ_
RooRealVar * a1_qqZZ_
Definition: SuperMELA.h:89
SuperMELA::sig_BW_
MELARelBWUFParam * sig_BW_
Definition: SuperMELA.h:81
SuperMELA::corr_sigma_sig
RooRealVar * corr_sigma_sig
Definition: SuperMELA.h:78
SuperMELA::n2_CB_
RooFormulaVar * n2_CB_
Definition: SuperMELA.h:68
SuperMELA::width_BW_
RooRealVar * width_BW_
Definition: SuperMELA.h:84
SuperMELA::strChan_
std::string strChan_
Definition: SuperMELA.h:58
TVar.hh
MELAqqZZPdf_v2
Definition: MELACombinePdfs.h:18
SuperMELA::sig_CB_
MELADoubleCB * sig_CB_
Definition: SuperMELA.h:80
SuperMELA::a12_qqZZ_
RooRealVar * a12_qqZZ_
Definition: SuperMELA.h:100
MELAStreamHelpers.hh
SuperMELA::pathToCards_
std::string pathToCards_
Definition: SuperMELA.h:61
SuperMELA::SetDecayChannel
void SetDecayChannel(std::string myChan)
Definition: SuperMELA.cc:113
MELADoubleCB
Definition: MELACombinePdfs.h:102
SuperMELA::mean_BW_
RooRealVar * mean_BW_
Definition: SuperMELA.h:83
SuperMELA::a2_qqZZ_
RooRealVar * a2_qqZZ_
Definition: SuperMELA.h:90
SuperMELA.h
SuperMELA::alpha_CB_
RooFormulaVar * alpha_CB_
Definition: SuperMELA.h:67
SuperMELA::checkChannel
bool checkChannel()
Definition: SuperMELA.cc:598
Generated by   doxygen 1.8.17