JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
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
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:122
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