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.
test_JVBF.c
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <cmath>
4 #include <string>
5 #include <vector>
6 #include <fstream>
7 #include <cstdlib>
8 #include <iomanip>
9 #include "TMath.h"
10 #include "TLorentzVector.h"
11 #include "TLorentzRotation.h"
12 #include "TFile.h"
13 #include "TTree.h"
14 #include "TChain.h"
15 #include "TString.h"
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TH3F.h"
19 #include "TH1D.h"
20 #include "TH2D.h"
21 #include "TH3D.h"
22 #include "TProfile.h"
23 #include "TF1.h"
24 #include "TGraphErrors.h"
25 #include <JHUGenMELA/MELA/interface/Mela.h>
26 #include <JHUGenMELA/MELA/src/computeAngles.h>
27 
28 using namespace RooFit;
29 
30 using namespace std;
31 
32 namespace {
33  const float Zmass = 91.1876;
34  const float gamZ = 2.5;
35  const float M_muon = 0.105658389;
36  const float M_electron = 0.00051099907;
37  const float M_tau = 1.777;
38  const int PDG_electron=11,PDG_muon=13,PDG_tau=15;
39 }
40 
41 
42 
43 float getMCFMMELAWeight(Mela& myMela, int lepId[4], float angularOrdered[8], double ggcoupl[2],int useConstant=0){
44  float myprob=1.0;
45  int myflavor=-1;
46  if(abs(lepId[0])==abs(lepId[1]) &&
47  abs(lepId[0])==abs(lepId[2]) &&
48  abs(lepId[0])==abs(lepId[3])){
49  if(abs(lepId[0])==11) myflavor=1;
50  else myflavor=2;
51  }
52  else myflavor=3;
53 
54  if(myflavor>=0) myMela.computeP(angularOrdered[0],angularOrdered[1],angularOrdered[2],angularOrdered[3],
55  angularOrdered[4],angularOrdered[5],angularOrdered[6],angularOrdered[7],
56  myflavor,
57  myprob,
58  useConstant
59  );
60  return myprob;
61 };
62 float getJHUGenMELAWeight(Mela& myMela, int lepId[4], float angularOrdered[8], double selfDHvvcoupl[SIZE_HVV][2]){
63  float myprob=1.0;
64  int myflavor=-1;
65  if(abs(lepId[0])==abs(lepId[1]) &&
66  abs(lepId[0])==abs(lepId[2]) &&
67  abs(lepId[0])==abs(lepId[3])){
68  if(abs(lepId[0])==11) myflavor=1;
69  else myflavor=2;
70  }
71  else myflavor=3;
72 
73  if(myflavor>=0) myMela.computeP(angularOrdered[0],angularOrdered[1],angularOrdered[2],angularOrdered[3],
74  angularOrdered[4],angularOrdered[5],angularOrdered[6],angularOrdered[7],
75  myflavor,
76  selfDHvvcoupl,
77  myprob
78  );
79  return myprob;
80 };
81 float getSuperMELA(Mela& myMela, int lepId[4], float mZZ, TVar::SuperMelaSyst syst){
82  float myprob=1.0;
83  TVar::LeptonFlavor myflavor=TVar::Flavor_Dummy;
84  if(abs(lepId[0])==abs(lepId[1]) &&
85  abs(lepId[0])==abs(lepId[2]) &&
86  abs(lepId[0])==abs(lepId[3])){
87  if(abs(lepId[0])==11) myflavor=TVar::Flavor_4e;
88  else myflavor=TVar::Flavor_4mu;
89  }
90  else myflavor=TVar::Flavor_2e2mu;
91 
92  if(myflavor!=TVar::Flavor_Dummy) myMela.computePM4l(mZZ,
93  myflavor,
94  syst,
95  myprob
96  );
97  return myprob;
98 };
99 
100 double eta_to_costheta ( double eta ) { return cos(atan(exp(-eta))*2.); }
101 
102 void test_JVBF(int erg_tev=8, float mSample=0, bool isggH=false){
103  float mPOLE=mSample;
104  if (mPOLE<=0) mPOLE=125.6;
105  float wPOLE=4.15e-3;
106  char TREE_NAME[] = "SelectedTree";
107 
108  // TVar::VerbosityLevel verbosity = TVar::INFO;
109 
110  Mela mela(erg_tev, mPOLE);
111 
112  TFile* foutput;
113  if (!isggH){
114  if (mSample>0) foutput = new TFile(Form("HZZ4lTree_jvbfMELA_H%.0f_%iTeV.root", mSample, erg_tev), "recreate");
115  else foutput = new TFile(Form("HZZ4lTree_jvbfMELA_HAll_%iTeV.root", erg_tev), "recreate");
116  }
117  else{
118  if (mSample>0) foutput = new TFile(Form("HZZ4lTree_jvbfMELA_ggH%.0f_%iTeV.root", mSample, erg_tev), "recreate");
119  else foutput = new TFile(Form("HZZ4lTree_jvbfMELA_ggHAll_%iTeV.root", erg_tev), "recreate");
120  }
121  TLorentzVector nullFourVector(0, 0, 0, 0);
122 
123  float MC_weight_noxsec;
124 
125  float pjvbf_VAJHU;
126  float pjvbf_VAJHU_first;
127  float pjvbf_VAJHU_second;
128  float phj_VAJHU_first;
129  float phj_VAJHU_second;
130  float pAux_vbf;
131  float pAux_vbf_first;
132  float pAux_vbf_second;
133 
134  float jet1Pt, jet2Pt;
135  float jet1Eta, jet2Eta;
136  float jet1Phi, jet2Phi;
137  float jet1E, jet2E;
138  float jet1Pt_Fake, jet2Pt_Fake;
139  float jet1Eta_Fake, jet2Eta_Fake;
140  float jet1Phi_Fake, jet2Phi_Fake;
141  float jet1E_Fake, jet2E_Fake;
142  float jet1px, jet1py, jet1pz;
143  float jet2px, jet2py, jet2pz;
144  float ZZPx, ZZPy, ZZPz, ZZE, dR;
145  short NJets30;
146  std::vector<double> * JetPt=0;
147  std::vector<double> * JetEta=0;
148  std::vector<double> * JetPhi=0;
149  std::vector<double> * JetMass=0;
150  std::vector<double> myJetPt;
151  std::vector<double> myJetCosTheta;
152  std::vector<double> myJetEta;
153  std::vector<double> myJetPhi;
154  std::vector<double> myJetMass;
155  TBranch* bJetPt=0;
156  TBranch* bJetEta=0;
157  TBranch* bJetPhi=0;
158  TBranch* bJetMass=0;
159 
160  float ZZMass, ZZPt, ZZPhi, ZZEta;
161 
162  int GenLep1Id, GenLep2Id, GenLep3Id, GenLep4Id;
163 
164  TChain* tree = new TChain(TREE_NAME);
165  char* user_folder[3]={ "4mu", "4e", "2mu2e" };
166  TString cinput_main = "/scratch0/hep/ianderso/CJLST/140519/PRODFSR";
167  if (erg_tev==8) cinput_main.Append("_8TeV");
168 
169 // TString cinput_main = "/afs/cern.ch/work/u/usarica/HZZ4l-125p6-FullAnalysis/LHC_";
170 // cinput_main.Append(Form("%iTeV", erg_tev));
171 
172  for (int ff=0; ff<3; ff++){
173  if (!isggH){
174  if (mSample>0) tree->Add(Form("%s/%s/HZZ4lTree_VBFH%.0f.root", cinput_main.Data(), user_folder[ff], mSample));
175  else tree->Add(Form("%s/%s/HZZ4lTree_VBFH*.root", cinput_main.Data(), user_folder[ff]));
176  }
177  else{
178  if (mSample>0) tree->Add(Form("%s/%s/HZZ4lTree_minloH%.0f.root", cinput_main.Data(), user_folder[ff], mSample));
179  else tree->Add(Form("%s/%s/HZZ4lTree_minloH*.root", cinput_main.Data(), user_folder[ff]));
180  }
181  }
182  tree->SetBranchAddress("MC_weight_noxsec", &MC_weight_noxsec);
183  tree->SetBranchAddress("NJets30", &NJets30);
184  tree->SetBranchAddress("JetPt", &JetPt, &bJetPt);
185  tree->SetBranchAddress("JetEta", &JetEta, &bJetEta);
186  tree->SetBranchAddress("JetPhi", &JetPhi, &bJetPhi);
187  tree->SetBranchAddress("JetMass", &JetMass, &bJetMass);
188  tree->SetBranchAddress("ZZMass", &ZZMass);
189  tree->SetBranchAddress("ZZPt", &ZZPt);
190  tree->SetBranchAddress("ZZEta", &ZZEta);
191  tree->SetBranchAddress("ZZPhi", &ZZPhi);
192 
193  TTree* newtree = new TTree("TestTree", "");
194  newtree->Branch("MC_weight_noxsec", &MC_weight_noxsec);
195  newtree->Branch("ZZMass", &ZZMass);
196  newtree->Branch("pAux_vbf", &pAux_vbf);
197  newtree->Branch("pAux_vbf_first", &pAux_vbf_first);
198  newtree->Branch("pAux_vbf_second", &pAux_vbf_second);
199  newtree->Branch("pjvbf_VAJHU", &pjvbf_VAJHU);
200  newtree->Branch("pjvbf_VAJHU_first", &pjvbf_VAJHU_first);
201  newtree->Branch("pjvbf_VAJHU_second", &pjvbf_VAJHU_second);
202  newtree->Branch("phj_VAJHU_first", &phj_VAJHU_first);
203  newtree->Branch("phj_VAJHU_second", &phj_VAJHU_second);
204  newtree->Branch("NJets30", &NJets30);
205  newtree->Branch("jet1Pt", &jet1Pt);
206  newtree->Branch("jet1Eta", &jet1Eta);
207  newtree->Branch("jet1Phi", &jet1Phi);
208  newtree->Branch("jet1E", &jet1E);
209  newtree->Branch("jet2Pt", &jet2Pt);
210  newtree->Branch("jet2Eta", &jet2Eta);
211  newtree->Branch("jet2Phi", &jet2Phi);
212  newtree->Branch("jet2E", &jet2E);
213  newtree->Branch("jet1_Fake_Pt", &jet1Pt_Fake);
214  newtree->Branch("jet1_Fake_Eta", &jet1Eta_Fake);
215  newtree->Branch("jet1_Fake_Phi", &jet1Phi_Fake);
216  newtree->Branch("jet1_Fake_E", &jet1E_Fake);
217  newtree->Branch("jet2_Fake_Pt", &jet2Pt_Fake);
218  newtree->Branch("jet2_Fake_Eta", &jet2Eta_Fake);
219  newtree->Branch("jet2_Fake_Phi", &jet2Phi_Fake);
220  newtree->Branch("jet2_Fake_E", &jet2E_Fake);
221  newtree->Branch("JetPt", &myJetPt);
222 // newtree->Branch("JetCosTheta", &myJetCosTheta);
223  newtree->Branch("JetPhi", &myJetPhi);
224  newtree->Branch("JetMass", &myJetMass);
225 
226 
227  int nEntries = tree->GetEntries();
228  double selfDHggcoupl[SIZE_HGG][2] ={ { 0 } };
229  double selfDHvvcoupl[SIZE_HVV_VBF][2] ={ { 0 } };
230  double selfDHwwcoupl[SIZE_HWW_VBF][2] ={ { 0 } };
231  double ggvvcoupl[2]={ 0, 0 };
232  mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, TVar::ZZGG);
233  int recorded=0;
234  for (int ev = 0; ev < nEntries; ev++){
235  pjvbf_VAJHU=-1;
236  pjvbf_VAJHU_first=-1;
237  pjvbf_VAJHU_second=-1;
238  jet1Pt=-1;
239  jet2Pt=-1;
240  jet1Eta=0;
241  jet2Eta=0;
242 
243  tree->GetEntry(ev);
244  GenLep1Id=11;
245  GenLep2Id=-11;
246  GenLep3Id=11;
247  GenLep4Id=-11;
248 
249  myJetPt.clear();
250  myJetEta.clear();
251  myJetPhi.clear();
252  myJetMass.clear();
253  myJetCosTheta.clear();
254 
255  TLorentzVector jet1(0, 0, 1e-3, 1e-3), jet2(0, 0, 1e-3, 1e-3), higgs(0, 0, 0, 0);
256  TLorentzVector p4[3], jets[2];
257  higgs.SetPtEtaPhiM(ZZPt, ZZEta, ZZPhi, ZZMass);
258  for (int i = 0; i < NJets30; i++){
259  myJetPt.push_back(JetPt->at(i));
260  myJetEta.push_back(JetEta->at(i));
261  myJetPhi.push_back(JetPhi->at(i));
262  myJetMass.push_back(JetMass->at(i));
263  myJetCosTheta.push_back(eta_to_costheta(JetEta->at(i)));
264  }
265 
266  int filled = 0;
267  if (myJetPt.size()>=1){
268 
269  jets[0].SetPxPyPzE(0, 0, 0, 1);
270  jets[1].SetPxPyPzE(0, 0, 0, 1);
271  for (int i = 0; i < myJetPt.size(); i++){
272  jets[filled].SetPtEtaPhiM(myJetPt[i], myJetEta[i], myJetPhi[i], myJetMass[i]);
273 
274  if (filled==0){
275  double jetE = jets[filled].Energy();
276  double jetP = jets[filled].P();
277  double ratio = (jetP>0 ? jetE/jetP : 1);
278  ratio = 1.;
279  jet1.SetPxPyPzE(jets[filled].Px()*ratio, jets[filled].Py()*ratio, jets[filled].Pz()*ratio, jetE);
280  filled++;
281  jet1Pt = jet1.Pt();
282  jet1Eta = jet1.Eta();
283  jet1Phi = jet1.Phi();
284  jet1E = jet1.E();
285  jet2.SetXYZT(0, 0, 0, 0);
286  }
287  else if(filled==1){
288  double jetE = jets[filled].Energy();
289  double jetP = jets[filled].P();
290  double ratio = (jetP>0 ? jetE/jetP : 1);
291  ratio = 1.;
292  jet2.SetXYZT(jets[filled].Px()*ratio, jets[filled].Py()*ratio, jets[filled].Pz()*ratio, jetE);
293  filled++;
294  jet2Pt = jet2.Pt();
295  jet2Eta = jet2.Eta();
296  jet2Phi = jet2.Phi();
297  jet2E = jet2.E();
298  }
299  else continue;
300 
301 /*
302  if (filled == 0){
303  if (jets[filled].Pt()>jet1.Pt()){
304  jet1.SetPxPyPzE(jets[filled].Px(), jets[filled].Py(), jets[filled].Pz(), jetE);
305  jet1Pt = jet1.Pt();
306  }
307  if (i == myJetPt.size() - 1){
308  filled++;
309  i = 0;
310  }
311  }
312  else{
313  if (jets[filled].Pt()<jet1.Pt() && jets[filled].Pt()>jet2.Pt()){
314  jet2.SetPxPyPzE(jets[filled].Px(), jets[filled].Py(), jets[filled].Pz(), jetE);
315  jet2Pt = jet2.Pt();
316  }
317  if (i == myJetPt.size() - 1){
318  filled++;
319  }
320  }
321  if (filled == 2) break;
322 */
323  }
324 //cos(atan(exp(-jet2Eta))*2)
325 
326  if (filled == 2){
327  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
328  mela.computeProdP(jet1, 2, jet2, 2, higgs, 25, nullFourVector, 0, pjvbf_VAJHU);
329  mela.get_PAux(pAux_vbf);
330 
331  TLorentzVector pTotal;
332  pTotal.SetXYZT(0, 0, 0, 0);
333  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
334  mela.computeProdP(jet1, 2, pTotal, 2, higgs, 25, nullFourVector, 0, pjvbf_VAJHU_first);
335  mela.get_PAux(pAux_vbf_first);
336  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JH);
337  mela.computeProdP(jet1, 2, pTotal, 2, higgs, 25, nullFourVector, 0, phj_VAJHU_first);
338  mela::computeFakeJet(jet1, higgs, pTotal);
339  jet2Pt_Fake = pTotal.Pt();
340  jet2Eta_Fake = pTotal.Eta();
341  jet2Phi_Fake = pTotal.Phi();
342  jet2E_Fake = pTotal.E();
343 
344  pTotal.SetXYZT(0, 0, 0, 0);
345  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
346  mela.computeProdP(pTotal, 2, jet2, 2, higgs, 25, nullFourVector, 0, pjvbf_VAJHU_second);
347  mela.get_PAux(pAux_vbf_second);
348  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JH);
349  mela.computeProdP(pTotal, 2, jet2, 2, higgs, 25, nullFourVector, 0, phj_VAJHU_second);
350  mela::computeFakeJet(jet2, higgs, pTotal);
351  jet1Pt_Fake = pTotal.Pt();
352  jet1Eta_Fake = pTotal.Eta();
353  jet1Phi_Fake = pTotal.Phi();
354  jet1E_Fake = pTotal.E();
355 
356  newtree->Fill();
357  }
358  else if (filled == 1){
359 /*
360  TLorentzVector pTotal = higgs+jet1;
361  pTotal.SetVect(-pTotal.Vect());
362  pTotal.SetE(pTotal.P());
363  jet2 = pTotal;
364 
365  jet2Pt = jet2.Pt();
366  jet2Eta = jet2.Eta();
367  jet2Phi = jet2.Phi();
368  jet2E = jet2.E();
369 
370 
371  jet1Pt_Fake = jet1.Pt();
372  jet1Eta_Fake = jet1.Eta();
373  jet1Phi_Fake = jet1.Phi();
374  jet1E_Fake = jet1.E();
375  jet2Pt_Fake = jet2.Pt();
376  jet2Eta_Fake = jet2.Eta();
377  jet2Phi_Fake = jet2.Phi();
378  jet2E_Fake = jet2.E();
379 
380  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
381  mela.computeProdP(jet1, 2, jet2, 2, higgs, 25, nullFourVector, 0, pjvbf_VAJHU);
382 */
383 //
384  jet2.SetXYZT(0,0,0,0);
385  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
386  mela.computeProdP(jet1, 2, jet2, 2, higgs, 25, nullFourVector, 0, pjvbf_VAJHU);
387  mela.get_PAux(pAux_vbf);
388  mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, TVar::JH);
389  mela.computeProdP(jet1, 2, jet2, 2, higgs, 25, nullFourVector, 0, phj_VAJHU_first);
390  mela::computeFakeJet(jet1, higgs, jet2);
391 /*
392  cout << "TEST:"
393  << " Higgs Pz: " << higgs.Pz()
394  << " Higgs P: " << higgs.P()
395  << " Higgs E: " << higgs.T()
396  << " Jet 1 Pz: " << jet1.Pz()
397  << " Jet 1 P: " << jet1.P()
398  << " Jet 1 E: " << jet1.T()
399  << " Jet 2 Pz: " << jet2.Pz()
400  << " Jet 2 P: " << jet2.P()
401  << " Jet 2 E: " << jet2.T() << '\n' << endl;
402 */
403 
404  jet2Pt = jet2.Pt();
405  jet2Eta = jet2.Eta();
406  jet2Phi = jet2.Phi();
407  jet2E = jet2.E();
408 
409  jet1Pt_Fake = jet1.Pt();
410  jet1Eta_Fake = jet1.Eta();
411  jet1Phi_Fake = jet1.Phi();
412  jet1E_Fake = jet1.E();
413  jet2Pt_Fake = jet2.Pt();
414  jet2Eta_Fake = jet2.Eta();
415  jet2Phi_Fake = jet2.Phi();
416  jet2E_Fake = jet2.E();
417 //
418  pjvbf_VAJHU_first = pjvbf_VAJHU;
419  pjvbf_VAJHU_second = pjvbf_VAJHU;
420  pAux_vbf_first = pAux_vbf;
421  pAux_vbf_second = pAux_vbf;
422  phj_VAJHU_second = phj_VAJHU_first;
423 
424  newtree->Fill();
425  }
426  }
427  }
428 
429 
430  foutput->WriteTObject(newtree);
431  delete newtree;
432  foutput->Close();
433  delete tree;
434 }
435 
436 void fit_JVBF(int erg_tev=8, float mSample=0){
437  float mPOLE=mSample;
438  if (mPOLE<=0) mPOLE=125.6;
439  float wPOLE=4.15e-3;
440  char TREE_NAME[] = "TestTree";
441 
442  TString cinput;
443  if (mSample>0) cinput = Form("HZZ4lTree_jvbfMELA_H%.0f_%iTeV.root", mSample, erg_tev);
444  else cinput = Form("HZZ4lTree_jvbfMELA_HAll_%iTeV.root", erg_tev);
445 
446  TFile* foutput;
447  if (mSample>0) foutput = new TFile(Form("jvbfMELA_Fits_H%.0f_%iTeV.root", mSample, erg_tev), "recreate");
448  else foutput = new TFile(Form("jvbfMELA_fits_wide_%iTeV.root", erg_tev), "recreate");
449 
450  TLorentzVector nullFourVector(0, 0, 0, 0);
451 
452  float MC_weight_noxsec;
453 
454  float pjvbf_VAJHU;
455  float pjvbf_VAJHU_first;
456  float pjvbf_VAJHU_second;
457 
458  float jet1Pt, jet2Pt;
459  float jet1Eta, jet2Eta;
460  float jet1Phi, jet2Phi;
461  float jet1E, jet2E;
462  float jet1Pt_Fake, jet2Pt_Fake;
463  float jet1Eta_Fake, jet2Eta_Fake;
464  float jet1Phi_Fake, jet2Phi_Fake;
465  float jet1E_Fake, jet2E_Fake;
466  float jet1px, jet1py, jet1pz;
467  float jet2px, jet2py, jet2pz;
468  float ZZPx, ZZPy, ZZPz, ZZE, dR;
469  short NJets30;
470 
471  float ZZMass;
472 
473  TChain* tree = new TChain(TREE_NAME);
474  tree->Add(cinput);
475  tree->SetBranchAddress("MC_weight_noxsec", &MC_weight_noxsec);
476  tree->SetBranchAddress("ZZMass", &ZZMass);
477  tree->SetBranchAddress("pjvbf_VAJHU", &pjvbf_VAJHU);
478  tree->SetBranchAddress("pjvbf_VAJHU_first", &pjvbf_VAJHU_first);
479  tree->SetBranchAddress("pjvbf_VAJHU_second", &pjvbf_VAJHU_second);
480  tree->SetBranchAddress("NJets30", &NJets30);
481  tree->SetBranchAddress("jet1Pt", &jet1Pt);
482  tree->SetBranchAddress("jet1Eta", &jet1Eta);
483  tree->SetBranchAddress("jet1Phi", &jet1Phi);
484  tree->SetBranchAddress("jet1E", &jet1E);
485  tree->SetBranchAddress("jet2Pt", &jet2Pt);
486  tree->SetBranchAddress("jet2Eta", &jet2Eta);
487  tree->SetBranchAddress("jet2Phi", &jet2Phi);
488  tree->SetBranchAddress("jet2E", &jet2E);
489  tree->SetBranchAddress("jet1_Fake_Pt", &jet1Pt_Fake);
490  tree->SetBranchAddress("jet1_Fake_Eta", &jet1Eta_Fake);
491  tree->SetBranchAddress("jet1_Fake_Phi", &jet1Phi_Fake);
492  tree->SetBranchAddress("jet1_Fake_E", &jet1E_Fake);
493  tree->SetBranchAddress("jet2_Fake_Pt", &jet2Pt_Fake);
494  tree->SetBranchAddress("jet2_Fake_Eta", &jet2Eta_Fake);
495  tree->SetBranchAddress("jet2_Fake_Phi", &jet2Phi_Fake);
496  tree->SetBranchAddress("jet2_Fake_E", &jet2E_Fake);
497 
498 
499  const int nbinsx=6;
500  double bins_mzz[nbinsx+1]={70, 140, 190, 240, 340, 520, 20000};
501  double sum_w[nbinsx]={ 0 };
502  double sum_mzz[nbinsx][2]={ { 0 } };
503  const int nbins_eta = 8;
504  double bins_eta[nbins_eta+1]={ 0, 0.75, 1.5, 2.25, 3, 4, 5.5, 10.5, 12 };
505 
506  TProfile* hProb[nbinsx];
507  TH1D* hProb_all;
508  TProfile* hEta[nbinsx+1];
509  TH1D* hDist[nbinsx+1];
510  for (int bin=0; bin<nbinsx+1; bin++){
511  if (bin<nbinsx) hProb[bin] = new TProfile(Form("njets1_pjvbf_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
512  else hProb_all = new TH1D("njets1_pjvbf_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
513  if (bin<nbinsx) hEta[bin] = new TProfile(Form("njets1_eta_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
514  else hEta[bin] = new TProfile("njets1_eta_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
515  if (bin<nbinsx) hDist[bin] = new TH1D(Form("njets1_dist_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
516  else hDist[bin] = new TH1D("njets1_dist_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
517  if (bin<nbinsx) hProb[bin]->Sumw2();
518  else hProb_all->Sumw2();
519  hEta[bin]->Sumw2();
520  hDist[bin]->Sumw2();
521  }
522 
523  int nEntries = tree->GetEntries();
524  cout << nEntries << endl;
525  for (int ev = 0; ev < nEntries; ev++){
526  tree->GetEntry(ev);
527 
528  int massbin=-1;
529  for (int bin=0; bin<nbinsx; bin++){
530  if (ZZMass>=bins_mzz[bin] && ZZMass<bins_mzz[bin+1]){
531  massbin=bin;
532  }
533  }
534  if (NJets30==1){
535  hProb[massbin]->Fill(fabs(jet2Eta_Fake), pjvbf_VAJHU_second, MC_weight_noxsec);
536  hEta[massbin]->Fill(fabs(jet2Eta_Fake), fabs(jet2Eta_Fake), MC_weight_noxsec);
537  hDist[massbin]->Fill(fabs(jet2Eta_Fake), MC_weight_noxsec);
538 
539  hEta[nbinsx]->Fill(fabs(jet2Eta_Fake), fabs(jet2Eta_Fake), MC_weight_noxsec);
540  hDist[nbinsx]->Fill(fabs(jet2Eta_Fake), MC_weight_noxsec);
541  }
542  }
543 
544  for (int bin=0; bin<nbinsx; bin++){
545  hProb[bin]->Scale(1./hProb[bin]->Integral("width"));
546  }
547  for (int bin=0; bin<nbinsx+1; bin++){
548  for (int bineta=1; bineta<=nbins_eta; bineta++){
549  double bincontent = hDist[bin]->GetBinContent(bineta);
550  double binerror = hDist[bin]->GetBinError(bineta);
551  double binwidth = hDist[bin]->GetXaxis()->GetBinWidth(bineta);
552 
553  bincontent /= binwidth;
554  binerror /= binwidth;
555 
556  hDist[bin]->SetBinContent(bineta, bincontent);
557  hDist[bin]->SetBinError(bineta, binerror);
558  }
559  hDist[bin]->Scale(1./hDist[bin]->Integral("width"));
560  }
561  for (int bineta=1; bineta<=nbins_eta; bineta++){
562  double sum_invwsq=0;
563  double sum_prob_invwsq=0;
564  for (int bin=0; bin<nbinsx; bin++){
565  double bincontent = hProb[bin]->GetBinContent(bineta);
566  double binerror = hProb[bin]->GetBinError(bineta);
567 
568  double invwsq = 0, prob_invwsq = 0;
569  if (binerror!=0){
570  invwsq = 1./pow(binerror, 2);
571  prob_invwsq = bincontent*invwsq;
572  }
573  sum_invwsq += invwsq;
574  sum_prob_invwsq += prob_invwsq;
575  }
576 
577  if (sum_invwsq>0){
578  sum_prob_invwsq /= sum_invwsq;
579  sum_invwsq = 1./sqrt(sum_invwsq);
580  }
581  hProb_all->SetBinContent(bineta, sum_prob_invwsq);
582  hProb_all->SetBinError(bineta, sum_invwsq);
583  }
584  hProb_all->Scale(1./hProb_all->Integral("width"));
585 
586 
587  double eta_center_all[nbins_eta];
588  double prob_center_all[nbins_eta];
589  double eta_center_all_err[nbins_eta];
590  double prob_center_all_err[nbins_eta];
591  int nNonEmpty_all = 0;
592  for (int bineta=1; bineta<=nbins_eta; bineta++){
593  if (hProb_all->GetBinError(bineta)>0){
594  eta_center_all[nNonEmpty_all] = hEta[nbinsx]->GetBinContent(bineta);
595  eta_center_all_err[nNonEmpty_all] = hEta[nbinsx]->GetBinError(bineta);
596  prob_center_all[nNonEmpty_all] = hProb_all->GetBinContent(bineta);
597  prob_center_all_err[nNonEmpty_all] = hProb_all->GetBinError(bineta);
598  cout << nNonEmpty_all << '\t' << prob_center_all[nNonEmpty_all] << '\t' << prob_center_all_err[nNonEmpty_all] << endl;
599  nNonEmpty_all++;
600  }
601  }
602 
603  TGraphErrors* tg_all = new TGraphErrors(nNonEmpty_all, eta_center_all, prob_center_all, eta_center_all_err, prob_center_all_err);
604  tg_all->SetNameTitle("tg_njets1_pjvbf_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]));
605  tg_all->GetXaxis()->SetTitle("#eta");
606  tg_all->GetYaxis()->SetTitle("<P_{VBF}>");
607 
608  TF1* fit_tg_all = new TF1("fit_tg_njets1_pjvbf_all", "gaus*(1+gaus(3))", bins_eta[0], bins_eta[nbins_eta]);
609  const int npars = 6;
610  double mypars[npars]={ 0.155, 5.193, 1.46, 0.362, 3.524, 0.715 };
611  fit_tg_all->SetParameters(mypars);
612  fit_tg_all->SetParLimits(0, 0, 1);
613  fit_tg_all->SetParLimits(3, 0, 1);
614  fit_tg_all->SetParLimits(1, 0, 20);
615  fit_tg_all->SetParLimits(4, 0, 20);
616  fit_tg_all->SetParLimits(2, 0, 20);
617  fit_tg_all->SetParLimits(5, 0, 20);
618  tg_all->Fit(fit_tg_all);
619 
620  double* par_postfit = fit_tg_all->GetParameters();
621  double* parerr_postfit = fit_tg_all->GetParErrors();
622  double parIndex[npars];
623  double parIndexErr[npars]={ 0 };
624  for (int ip=0; ip<npars; ip++) parIndex[ip] = (double)ip;
625 
626  double integral = 0.5*par_postfit[0]*(1.+TMath::Erf(par_postfit[1]/par_postfit[2]/TMath::Sqrt(2)));
627  double integral2 = 1./(2.*TMath::Sqrt(2.*TMath::Pi()*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))));
628  integral2 *= par_postfit[0]*par_postfit[3]*TMath::Exp(-0.5*pow(par_postfit[1]-par_postfit[4], 2)/(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2)))*
629  (1.+TMath::Erf((TMath::Power(par_postfit[2], 2)*par_postfit[4]+TMath::Power(par_postfit[5], 2)*par_postfit[1]) / (par_postfit[2]*par_postfit[5]*TMath::Sqrt(2.*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))))));
630  integral += integral2;
631 
632  cout << integral << endl;
633  par_postfit[0] /= integral;
634  parerr_postfit[0] /= integral;
635 
636  TGraphErrors* tg_pars = new TGraphErrors(npars, parIndex, par_postfit, parIndexErr, parerr_postfit);
637  tg_pars->SetNameTitle("njets1_pjvbf_pars", "gaus*(1+gaus(3))");
638  tg_pars->GetXaxis()->SetTitle("Parameter index");
639  tg_pars->GetYaxis()->SetTitle("Parameter");
640  foutput->WriteTObject(tg_pars);
641  delete tg_pars;
642 
643  foutput->WriteTObject(tg_all);
644  delete tg_all;
645 
646  for (int bin=0; bin<nbinsx+1; bin++){
647  if (bin<nbinsx) foutput->WriteTObject(hProb[bin]);
648  else foutput->WriteTObject(hProb_all);
649  foutput->WriteTObject(hEta[bin]);
650  foutput->WriteTObject(hDist[bin]);
651  delete hDist[bin];
652  delete hEta[bin];
653  if (bin<nbinsx) delete hProb[bin];
654  else delete hProb_all;
655  }
656  delete tree;
657 
658  foutput->Close();
659 }
660 
661 /*
662 gaus*(1+gaus(3))
663 ->
664 Integrate[u*g[a,b,x]*(1 + v*h[c,d.x]), {x, 0, Infinity},Assumptions -> b > 0 && d > 0 && a > 0 && c > 0]
665 ==
666 1/2 u (1 + Erf[a/(Sqrt[2] b)]) + (
667 E^(-((a - c)^2/(2 (b^2 + d^2))))
668 u v (1 + Erf[(b^2 c + a d^2)/(Sqrt[2] b d Sqrt[b^2 + d^2])]))/(
669 2 Sqrt[b^2 + d^2] Sqrt[2 \[Pi]])
670 */
anonymous_namespace{test_JVBF.c}::gamZ
const float gamZ
Definition: test_JVBF.c:34
TVar::JJVBF
@ JJVBF
Definition: TVar.hh:72
anonymous_namespace{test_JVBF.c}::M_tau
const float M_tau
Definition: test_JVBF.c:37
anonymous_namespace{test_JVBF.c}::M_electron
const float M_electron
Definition: test_JVBF.c:36
anonymous_namespace{test_JVBF.c}::PDG_muon
const int PDG_muon
Definition: test_JVBF.c:38
getJHUGenMELAWeight
float getJHUGenMELAWeight(Mela &myMela, int lepId[4], float angularOrdered[8], double selfDHvvcoupl[SIZE_HVV][2])
Definition: test_JVBF.c:62
anonymous_namespace{test_JVBF.c}::M_muon
const float M_muon
Definition: test_JVBF.c:35
mela
Definition: mela.py:1
test_JVBF
void test_JVBF(int erg_tev=8, float mSample=0, bool isggH=false)
Definition: test_JVBF.c:102
TVar::ZZGG
@ ZZGG
Definition: TVar.hh:61
anonymous_namespace{TCouplingsBase.hh}::SIZE_HGG
@ SIZE_HGG
Definition: TCouplingsBase.hh:40
modparameters::ev
real(8), parameter, public ev
Definition: mod_Parameters.F90:97
testME_all.ratio
ratio
Definition: testME_all.py:136
Mela::computeP
void computeP(float &prob, bool useConstant=true)
Computes the probability for the probabilities on the decay side of things using the constituent daug...
Definition: Mela.cc:1183
TVar::HSMHiggs
@ HSMHiggs
Definition: TVar.hh:126
Mela::computePM4l
void computePM4l(TVar::SuperMelaSyst syst, float &prob)
Definition: Mela.cc:1864
selfDHggcoupl
double selfDHggcoupl[nSupportedHiggses][SIZE_HGG][2]
Definition: raw_names.txt:1
TVar::JHUGen
@ JHUGen
Definition: TVar.hh:57
dd_global::cout
integer cout
Definition: DD_global.F90:21
Mela
Definition: Mela.h:48
TUtil::computeFakeJet
void computeFakeJet(TLorentzVector const &realJet, TLorentzVector const &others, TLorentzVector &fakeJet)
Definition: TUtil.cc:156
getSuperMELA
float getSuperMELA(Mela &myMela, int lepId[4], float mZZ, TVar::SuperMelaSyst syst)
Definition: test_JVBF.c:81
fit_JVBF
void fit_JVBF(int erg_tev=8, float mSample=0)
Definition: test_JVBF.c:436
TVar::SuperMelaSyst
SuperMelaSyst
Definition: TVar.hh:186
TVar::SelfDefine_spin0
@ SelfDefine_spin0
Definition: TVar.hh:180
anonymous_namespace{test_JVBF.c}::PDG_electron
const int PDG_electron
Definition: test_JVBF.c:38
getMCFMMELAWeight
float getMCFMMELAWeight(Mela &myMela, int lepId[4], float angularOrdered[8], double ggcoupl[2], int useConstant=0)
Definition: test_JVBF.c:43
anonymous_namespace{test_JVBF.c}::PDG_tau
const int PDG_tau
Definition: test_JVBF.c:38
SIZE_HVV
@ SIZE_HVV
Definition: raw_couplings.txt:57
eta_to_costheta
double eta_to_costheta(double eta)
Definition: test_JVBF.c:100
anonymous_namespace{test_JVBF.c}::Zmass
const float Zmass
Definition: test_JVBF.c:33
selfDHwwcoupl
double selfDHwwcoupl[nSupportedHiggses][SIZE_HVV][2]
Definition: raw_names.txt:9
Generated by   doxygen 1.8.17