-
Notifications
You must be signed in to change notification settings - Fork 2
/
HbbHbb_LMRSelection_chi2_150.cc
495 lines (450 loc) · 22.2 KB
/
HbbHbb_LMRSelection_chi2_150.cc
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TSystem.h>
#include <TROOT.h>
#include <TLorentzVector.h>
#include <iostream>
#include <vector>
//#include "PDFs/BtagCalibrationStandalone.cc"
#include "HbbHbb_Component_SignalPurity.cc"
#include "HbbHbb_Component_KinFit.cc"
double jet_pT_cut1=30.;
double mean_H1_mass_=150;
double sigma_H1_mass_=20;
double mean_H2_mass_=150;
double sigma_H2_mass_=20;
float mass_mean( int mass )
{
if( mass == 260 ) return 243;
else if( mass == 300 ) return 276;
else if( mass == 350 ) return 318;
else if( mass == 400 ) return 363;
else if( mass == 450 ) return 413;
else if( mass == 500 ) return 463;
else if( mass == 550 ) return 530;
else return mass;
}
float mass_mean_reg( int mass )
{
if( mass == 260 ) return 257;
else if( mass == 300 ) return 291;
else if( mass == 350 ) return 336;
else if( mass == 400 ) return 387;
else if( mass == 450 ) return 433;
else if( mass == 500 ) return 481;
else if( mass == 550 ) return 531;
else return mass;
}
float mass_sigma( int mass )
{
if( mass == 260 ) return 20;
else if( mass == 300 ) return 25;
else if( mass == 350 ) return 31;
else if( mass == 400 ) return 36;
else if( mass == 450 ) return 32;
else if( mass == 500 ) return 35;
else if( mass == 550 ) return 35;
else return 30;
}
float mass_sigma_reg( int mass )
{
if( mass == 260 ) return 18;
else if( mass == 300 ) return 22;
else if( mass == 350 ) return 26;
else if( mass == 400 ) return 24;
else if( mass == 450 ) return 27;
else if( mass == 500 ) return 31;
else if( mass == 550 ) return 30;
else return 30;
}
TLorentzVector fillTLorentzVector(double pT, double eta, double phi, double M)
{
TLorentzVector jet_p4;
jet_p4.SetPtEtaPhiM(pT, eta, phi, M);
return jet_p4;
}
void HbbHbb_LMRSelection_chi2_150(std::string type, std::string sample, int signal_mass = 300, bool reg = false )
{
std::string inputfilename="../PreSelected_"+sample+".root";
TChain *tree=new TChain("tree");
tree->Add(inputfilename.c_str());
//std::cout<<"Opened input file "<<inputfilename<<std::endl;
// Book variables
int evt;
float eventWeight;
int nJets, nGenBQuarkFromH;
float jet_btagCSV[100];
float jet_pT[100], jet_eta[100], jet_phi[100], jet_mass[100];
float genBQuarkFromH_pT[100],genBQuarkFromH_eta[100],genBQuarkFromH_phi[100],genBQuarkFromH_mass[100];
float jet_regressed_pT[100], jet_flavor[100];
std::vector<unsigned int> *jetIndex_CentralpT40btag_deepCSVOrder=0;
// Retrieve variables
tree->SetBranchAddress("evt", &evt);
tree->SetBranchAddress("eventWeight", &(eventWeight));
tree->SetBranchAddress("nJet", &(nJets));
tree->SetBranchAddress("Jet_btagCSV", &(jet_btagCSV));
tree->SetBranchAddress("Jet_pt", &(jet_pT));
tree->SetBranchAddress("Jet_eta", &(jet_eta));
tree->SetBranchAddress("Jet_phi", &(jet_phi));
tree->SetBranchAddress("Jet_mass", &(jet_mass));
tree->SetBranchAddress("Jet_regressed_pt", &(jet_regressed_pT));
tree->SetBranchAddress("jetIndex_CentralpT40btag_deepCSVOrder", &(jetIndex_CentralpT40btag_deepCSVOrder));
if(type!="Data"){
tree->SetBranchAddress("nGenBQuarkFromH", &(nGenBQuarkFromH));
tree->SetBranchAddress("GenBQuarkFromH_pt", &(genBQuarkFromH_pT));
tree->SetBranchAddress("GenBQuarkFromH_eta", &(genBQuarkFromH_eta));
tree->SetBranchAddress("GenBQuarkFromH_phi", &(genBQuarkFromH_phi));
tree->SetBranchAddress("GenBQuarkFromH_mass", &(genBQuarkFromH_mass));
tree->SetBranchAddress("Jet_mcFlavour", &(jet_flavor));
}
// Book histograms
TH1F *h_H1_mass = new TH1F("h_H1_mass", "; m_{H1} (GeV)", 100, 50., 250.);
TH1F *h_H1_pT = new TH1F("h_H1_pT", "; H1 p_{T} (GeV/c)", 800, 0., 800.);
TH1F *h_H2_mass = new TH1F("h_H2_mass", "; m_{H2} (GeV)", 100, 50., 250.);
TH1F *h_H2_pT = new TH1F("h_H2_pT", "; H2 p_{T} (GeV/c)", 800, 0., 800.);
TH1F *h_HH_balance = new TH1F("h_HH_balance", "; (#vec{p}_{H1} + #vec{p}_{H2} - #vec{p}_{X}^{gen})_{T} GeV", 200, 0, 200.);
TH2F *h_mH1_mH2_asym = new TH2F("h_mH1_mH2_asym", "; m_{H1} (GeV); m_{H2} (GeV)", 300, 0., 300., 300, 0., 300.);
TH1F *h_H1_mass_biasCorrected = new TH1F("h_H1_mass_biasCorrected", "; Bias Corrected m_{H1} (GeV)", 300, 0., 300.);
TH1F *h_H1_pT_biasCorrected = new TH1F("h_H1_pT_biasCorrected", "; H1 p_{T} (GeV/c)", 800, 0., 800.);
TH1F *h_H2_mass_biasCorrected = new TH1F("h_H2_mass_biasCorrected", "; Bias Corrected m_{H2} (GeV)", 300, 0., 300.);
TH1F *h_H2_pT_biasCorrected = new TH1F("h_H2_pT_biasCorrected", "; H2 p_{T} (GeV/c)", 800, 0., 800.);
TH1F *h_HH_balance_biasCorrected = new TH1F("h_HH_balance_biasCorrected", "; (#vec{p}_{H1} + #vec{p}_{H2} - #vec{p}_{X}^{gen})_{T} GeV", 200, 0, 200.);
TH2F *h_mH1_mH2_asym_biasCorrected = new TH2F("h_mH1_mH2_asym_biasCorrected", "; m_{H1} (GeV); m_{H2} (GeV)", 300, 0., 300., 300, 0., 300.);
TH1F *h_GenX_pT = new TH1F("h_GenX_pT", "; (#vec{p}_{H1} + #vec{p}_{H2})_{T} GeV", 200, 0., 800.);
TH1F *h_kinFitchi2=new TH1F("h_kinFitchi2", "; Event 4 jet kinematic #chi^2", 200, 0., 10.);
TH1F *h_chi=new TH1F("h_chi", "; HH #chi", 100, 0, 100);
TH1F *h_chi_biasCorrected=new TH1F("h_chi_biasCorrected", "; HH #chi", 100, 0, 100);
TH1F *h_mX_SR = new TH1F("h_mX_SR", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR->Sumw2();
TH1F *h_mX_SR_purity0 = new TH1F("h_mX_SR_purity0", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity0->Sumw2();
TH1F *h_mX_SR_purity1 = new TH1F("h_mX_SR_purity1", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity1->Sumw2();
TH1F *h_mX_SR_purity2 = new TH1F("h_mX_SR_purity2", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity2->Sumw2();
TH1F *h_mX_SR_purity3 = new TH1F("h_mX_SR_purity3", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity3->Sumw2();
TH1F *h_mX_SR_purity4 = new TH1F("h_mX_SR_purity4", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity4->Sumw2();
TH1F *h_mX_SR_purity5 = new TH1F("h_mX_SR_purity5", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_purity5->Sumw2();
TH1F *h_mX_SR_biasCorrected = new TH1F("h_mX_SR_biasCorrected", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_biasCorrected->Sumw2();
TH1F *h_mX_SR_kinFit = new TH1F("h_mX_SR_kinFit", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit->Sumw2();
TH1F *h_mX_SR_kinFit_purity0 = new TH1F("h_mX_SR_kinFit_purity0", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity0->Sumw2();
TH1F *h_mX_SR_kinFit_purity1 = new TH1F("h_mX_SR_kinFit_purity1", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity1->Sumw2();
TH1F *h_mX_SR_kinFit_purity2 = new TH1F("h_mX_SR_kinFit_purity2", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity2->Sumw2();
TH1F *h_mX_SR_kinFit_purity3 = new TH1F("h_mX_SR_kinFit_purity3", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity3->Sumw2();
TH1F *h_mX_SR_kinFit_purity4 = new TH1F("h_mX_SR_kinFit_purity4", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity4->Sumw2();
TH1F *h_mX_SR_kinFit_purity5 = new TH1F("h_mX_SR_kinFit_purity5", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SR_kinFit_purity5->Sumw2();
TH1F *h_HH_balance_kinFit = new TH1F("h_HH_balance_kinFit", "; (#vec{p}_{H1} + #vec{p}_{H2} - #vec{p}_{X}^{gen})_{T} GeV", 200, 0, 200.);
TH1F *h_mX_SB = new TH1F("h_mX_SB", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SB->Sumw2();
TH1F *h_mX_SB_biasCorrected = new TH1F("h_mX_SB_biasCorrected", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SB_biasCorrected->Sumw2();
TH1F *h_mX_SB_kinFit = new TH1F("h_mX_SB_kinFit", "; m_{X} (GeV)", 3000, 0., 3000.); h_mX_SB_kinFit->Sumw2();
std::string Old_histfilename="../Histograms_Preselected_"+sample+".root";
std::string histfilename="Histograms_LMR_150"+sample+".root";
gSystem->Exec(("cp "+Old_histfilename+" "+histfilename).c_str());
TFile *tFile1=new TFile((Old_histfilename).c_str(), "READ");
TH1F h_Cuts=*((TH1F*)((TH1F*)tFile1->Get("h_Cuts"))->Clone("h_Cuts"));
tFile1->Close();
float SR_tot=0, SB_tot=0;
int masses[]={260,300,350,400,450,500,550};
std::map<int,float> SB_m, SR_m;
for( auto m : masses ){ SR_m[m]=0; SB_m[m]=0; }
// Event loop
int nEvents=tree->GetEntries();
double nCut4=0, nCut5=0, nCutGen=0;
for (int i=0; i<tree->GetEntries(); ++i)
{
tree->GetEvent(i);
bool foundHH=false;
double chi2_old=200.;
double m_diff_old=100.;
int H1jet1_i=-1, H1jet2_i=-1;
int H2jet1_i=-1, H2jet2_i=-1;
for (unsigned int j=0; j<jetIndex_CentralpT40btag_deepCSVOrder->size(); ++j)
{
unsigned int j_jetIndex=jetIndex_CentralpT40btag_deepCSVOrder->at(j);
TLorentzVector jet1_p4, jet2_p4, jet3_p4, jet4_p4;
jet1_p4=fillTLorentzVector(jet_regressed_pT[j_jetIndex], jet_eta[j_jetIndex], jet_phi[j_jetIndex], jet_mass[j_jetIndex]);
if (jet1_p4.Pt()>jet_pT_cut1)
{
for (unsigned int k=0; k<jetIndex_CentralpT40btag_deepCSVOrder->size(); ++k)
{
unsigned int k_jetIndex=jetIndex_CentralpT40btag_deepCSVOrder->at(k);
jet2_p4=fillTLorentzVector(jet_regressed_pT[k_jetIndex], jet_eta[k_jetIndex], jet_phi[k_jetIndex], jet_mass[k_jetIndex]);
if (k_jetIndex!=j_jetIndex && jet2_p4.Pt()>jet_pT_cut1)
{
for (unsigned int l=0; l<jetIndex_CentralpT40btag_deepCSVOrder->size(); ++l)
{
unsigned int l_jetIndex=jetIndex_CentralpT40btag_deepCSVOrder->at(l);
jet3_p4=fillTLorentzVector(jet_regressed_pT[l_jetIndex], jet_eta[l_jetIndex], jet_phi[l_jetIndex], jet_mass[l_jetIndex]);
if (l_jetIndex!=k_jetIndex && l_jetIndex!=j_jetIndex && jet3_p4.Pt()>jet_pT_cut1)
{
for (unsigned int m=0; m<jetIndex_CentralpT40btag_deepCSVOrder->size(); ++m)
{
unsigned int m_jetIndex=jetIndex_CentralpT40btag_deepCSVOrder->at(m);
jet4_p4=fillTLorentzVector(jet_regressed_pT[m_jetIndex], jet_eta[m_jetIndex], jet_phi[m_jetIndex], jet_mass[m_jetIndex]);
if (m_jetIndex!=j_jetIndex && m_jetIndex!=k_jetIndex && m_jetIndex!=l_jetIndex && jet4_p4.Pt()>jet_pT_cut1)
{
// swap if H pT is odd in second decimal place
if (int((jet1_p4+jet2_p4).Pt()*100.) % 2 == 1)
{
swap(j_jetIndex, l_jetIndex);
swap(k_jetIndex, m_jetIndex);
swap(jet1_p4, jet3_p4);
swap(jet2_p4, jet4_p4);
}
double deltaR1=jet1_p4.DeltaR(jet2_p4);
double deltaR2=jet3_p4.DeltaR(jet4_p4);
TLorentzVector diJet1_p4=jet1_p4+jet2_p4;
TLorentzVector diJet2_p4=jet3_p4+jet4_p4;
double mH1=diJet1_p4.M();
double mH2=diJet2_p4.M();
double m_diff=fabs(mH1-mH2);
double chi2=pow((mH1-mean_H1_mass_)/sigma_H1_mass_, 2)+pow((mH2-mean_H2_mass_)/sigma_H2_mass_, 2);
if (chi2< chi2_old && ((110.<mH1 && mH1<190.) && (110.<mH2 && mH2<190.)))
{
H1jet1_i=j_jetIndex;
H1jet2_i=k_jetIndex;
H2jet1_i=l_jetIndex;
H2jet2_i=m_jetIndex;
chi2_old=chi2;
m_diff_old=m_diff;
foundHH=true;
}
} // Conditions on 4th jet
} // Loop over 4th jet
} // Conditions on 3rd jet
} // Loop over 3rd jet
} // Conditions on 2nd jet
} // Loop over 2nd jet
} // Condition of 1st jet
} // Loop over 1st jet
if (foundHH)
{
nCut4+=eventWeight;
double chi=pow(chi2_old, 0.5);
h_chi->Fill(chi, eventWeight);
TLorentzVector jet1_p4=fillTLorentzVector(jet_regressed_pT[H1jet1_i], jet_eta[H1jet1_i], jet_phi[H1jet1_i], jet_mass[H1jet1_i]);
TLorentzVector jet2_p4=fillTLorentzVector(jet_regressed_pT[H1jet2_i], jet_eta[H1jet2_i], jet_phi[H1jet2_i], jet_mass[H1jet2_i]);
TLorentzVector jet3_p4=fillTLorentzVector(jet_regressed_pT[H2jet1_i], jet_eta[H2jet1_i], jet_phi[H2jet1_i], jet_mass[H2jet1_i]);
TLorentzVector jet4_p4=fillTLorentzVector(jet_regressed_pT[H2jet2_i], jet_eta[H2jet2_i], jet_phi[H2jet2_i], jet_mass[H2jet2_i]);
TLorentzVector jet1_p4_unregressed=fillTLorentzVector(jet_pT[H1jet1_i], jet_eta[H1jet1_i], jet_phi[H1jet1_i], jet_mass[H1jet1_i]);
TLorentzVector jet2_p4_unregressed=fillTLorentzVector(jet_pT[H1jet2_i], jet_eta[H1jet2_i], jet_phi[H1jet2_i], jet_mass[H1jet2_i]);
TLorentzVector jet3_p4_unregressed=fillTLorentzVector(jet_pT[H2jet1_i], jet_eta[H2jet1_i], jet_phi[H2jet1_i], jet_mass[H2jet1_i]);
TLorentzVector jet4_p4_unregressed=fillTLorentzVector(jet_pT[H2jet2_i], jet_eta[H2jet2_i], jet_phi[H2jet2_i], jet_mass[H2jet2_i]);
// Fill histograms before bias correction
TLorentzVector H1_p4=jet1_p4+jet2_p4;
TLorentzVector H2_p4=jet3_p4+jet4_p4;
TLorentzVector X_p4=H1_p4+H2_p4;
double pTH1=H1_p4.Pt();
double pTH2=H2_p4.Pt();
double mH1=H1_p4.M();
double mH2=H2_p4.M();
h_H1_mass->Fill(mH1, eventWeight);
h_H1_pT->Fill(pTH1, eventWeight);
h_H2_mass->Fill(mH2, eventWeight);
h_H2_pT->Fill(pTH2, eventWeight);
h_mH1_mH2_asym->Fill(mH1, mH2, eventWeight);
// Apply bias correction
TLorentzVector jet1_p4_biasCorrected=biasEt_signal(jet1_p4_unregressed);
TLorentzVector jet2_p4_biasCorrected=biasEt_signal(jet2_p4_unregressed);
TLorentzVector jet3_p4_biasCorrected=biasEt_signal(jet3_p4_unregressed);
TLorentzVector jet4_p4_biasCorrected=biasEt_signal(jet4_p4_unregressed);
// Fill histograms after bias correction
TLorentzVector H1_p4_biasCorrected=jet1_p4_biasCorrected+jet2_p4_biasCorrected;
TLorentzVector H2_p4_biasCorrected=jet3_p4_biasCorrected+jet4_p4_biasCorrected;
TLorentzVector X_p4_biasCorrected=H1_p4_biasCorrected+H2_p4_biasCorrected;
double pTH1_biasCorrected=H1_p4_biasCorrected.Pt();
double pTH2_biasCorrected=H2_p4_biasCorrected.Pt();
double mH1_biasCorrected=H1_p4_biasCorrected.M();
double mH2_biasCorrected=H2_p4_biasCorrected.M();
h_H1_mass_biasCorrected->Fill(mH1_biasCorrected, eventWeight);
h_H1_pT_biasCorrected->Fill(pTH1_biasCorrected, eventWeight);
h_H2_mass_biasCorrected->Fill(mH2_biasCorrected, eventWeight);
h_H2_pT_biasCorrected->Fill(pTH2_biasCorrected, eventWeight);
h_mH1_mH2_asym_biasCorrected->Fill(mH1_biasCorrected, mH2_biasCorrected, eventWeight);
// Check purity of jet selection here // FIX THIS to check against kin fit reco jets.
TLorentzVector b1_p4;
TLorentzVector b2_p4;
TLorentzVector b3_p4;
TLorentzVector b4_p4;
int purity=-3;
if (type=="Signal")
{
if (nGenBQuarkFromH==4)
{
b1_p4=fillTLorentzVector(genBQuarkFromH_pT[0], genBQuarkFromH_eta[0], genBQuarkFromH_phi[0], genBQuarkFromH_mass[0]);
b2_p4=fillTLorentzVector(genBQuarkFromH_pT[1], genBQuarkFromH_eta[1], genBQuarkFromH_phi[1], genBQuarkFromH_mass[1]);
b3_p4=fillTLorentzVector(genBQuarkFromH_pT[2], genBQuarkFromH_eta[2], genBQuarkFromH_phi[2], genBQuarkFromH_mass[2]);
b4_p4=fillTLorentzVector(genBQuarkFromH_pT[3], genBQuarkFromH_eta[3], genBQuarkFromH_phi[3], genBQuarkFromH_mass[3]);
TLorentzVector j[4]={jet1_p4, jet2_p4, jet3_p4, jet4_p4};
TLorentzVector b[4]={b1_p4, b2_p4, b3_p4, b4_p4};
int jMatchedbindex[4]={-1, -1, -1, -1};
purity=purityTest(j, b, jMatchedbindex);
}
else
{
std::cout<<"ERROR: This is a signal sample without 4 gen b from H."<<std::endl;
}
}
double oppSign=(mH1-mean_H1_mass_)*(mH2-mean_H2_mass_);
if (chi<=1) // Signal Region
{
nCut5+=eventWeight;
SR_tot++;
for( auto m : masses )
{
if( !reg && X_p4.M() > mass_mean(m)-mass_sigma(m)*3 && X_p4.M() < mass_mean(m)+mass_sigma(m)*3 ) SR_m[m]++;
else if( reg && X_p4.M() > mass_mean_reg(m)-mass_sigma_reg(m)*3 && X_p4.M() < mass_mean_reg(m)+mass_sigma_reg(m)*3 ) SR_m[m]++;
}
// Apply kinematic constraint
// jet1_p4, jet2_p4, jet3_p4, jet4_p4 will change values
double kinFitchi2=constrainHH_signalMeasurement(&jet1_p4, &jet2_p4, &jet3_p4, &jet4_p4);
// double kinFitchi2=constrainHH_afterRegression(&jet1_p4, &jet2_p4, &jet3_p4, &jet4_p4);
h_kinFitchi2->Fill(kinFitchi2, eventWeight);
TLorentzVector X_p4_kinFit=(jet1_p4+jet2_p4+jet3_p4+jet4_p4);
h_mX_SR->Fill(X_p4.M(), eventWeight);
h_mX_SR_biasCorrected->Fill(X_p4_biasCorrected.M(), eventWeight);
h_mX_SR_kinFit->Fill(X_p4_kinFit.M(), eventWeight);
if (purity==-1)
{
h_mX_SR_purity5->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity5->Fill(X_p4_kinFit.M(), eventWeight);
}
else if (purity==0)
{
h_mX_SR_purity0->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity0->Fill(X_p4_kinFit.M(), eventWeight);
}
else if (purity==1)
{
h_mX_SR_purity1->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity1->Fill(X_p4_kinFit.M(), eventWeight);
}
else if (purity==2)
{
h_mX_SR_purity2->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity2->Fill(X_p4_kinFit.M(), eventWeight);
}
else if (purity==3)
{
h_mX_SR_purity3->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity3->Fill(X_p4_kinFit.M(), eventWeight);
}
else if (purity==4)
{
h_mX_SR_purity4->Fill(X_p4.M(), eventWeight);
h_mX_SR_kinFit_purity4->Fill(X_p4_kinFit.M(), eventWeight);
}
// Fill HH pT balancing histograms
if (type=="Signal" && nGenBQuarkFromH==4)
{
TLorentzVector gen_X_p4 = b1_p4 + b2_p4 + b3_p4 + b4_p4;
h_HH_balance->Fill((X_p4 - gen_X_p4).Pt(), eventWeight);
h_HH_balance_biasCorrected->Fill((X_p4_biasCorrected - gen_X_p4).Pt(), eventWeight);
h_HH_balance_kinFit->Fill((X_p4_kinFit - gen_X_p4).Pt(), eventWeight);
}
}
else if (1<chi && chi<2. && oppSign<0) // Sideband Region
{
SB_tot++;
for( auto m : masses )
{
if( !reg && X_p4.M() > mass_mean(m)-mass_sigma(m)*3 && X_p4.M() < mass_mean(m)+mass_sigma(m)*3 ) SB_m[m]++;
else if( reg && X_p4.M() > mass_mean_reg(m)-mass_sigma_reg(m)*3 && X_p4.M() < mass_mean_reg(m)+mass_sigma_reg(m)*3 ) SB_m[m]++;
}
// Apply kinematic constraint
// jet1_p4, jet2_p4, jet3_p4, jet4_p4 will change values
double kinFitchi2=constrainHH_signalMeasurement(&jet1_p4, &jet2_p4, &jet3_p4, &jet4_p4);
TLorentzVector X_p4_kinFit=(jet1_p4+jet2_p4+jet3_p4+jet4_p4);
h_mX_SB->Fill(X_p4.M(), eventWeight);
h_mX_SB_biasCorrected->Fill(X_p4_biasCorrected.M(), eventWeight);
h_mX_SB_kinFit->Fill(X_p4_kinFit.M(), eventWeight);
}
}
if (i%(nEvents/10)==0) std::cout<<int(i*100./nEvents)+1<<"% of "<<nEvents<<" events have been processed."<<std::endl;
} // Event loop
if(type!="Data") std::cout << signal_mass << " " << SR_m[signal_mass] << std::endl;
else for( auto m : masses ) std::cout << m << " " << SB_m[m]*SR_tot/SB_tot << std::endl;
h_Cuts.Fill(9, nCut4); // HH Candidates
h_Cuts.Fill(11, nCut5); // SR
TFile *tFile2=new TFile(histfilename.c_str(), "UPDATE");
tFile2->Delete("h_Cuts;1");
h_H1_mass->Write();
h_H1_pT->Write();
h_H2_mass->Write();
h_H2_pT->Write();
h_HH_balance->Write();
h_mH1_mH2_asym->Write();
h_H1_mass_biasCorrected->Write();
h_H1_pT_biasCorrected->Write();
h_H2_mass_biasCorrected->Write();
h_H2_pT_biasCorrected->Write();
h_HH_balance_biasCorrected->Write();
h_mH1_mH2_asym_biasCorrected->Write();
h_GenX_pT->Write();
h_kinFitchi2->Write();
h_chi->Write();
h_chi_biasCorrected->Write();
if(type!="Data"){ h_mX_SR->Write();
h_mX_SR_biasCorrected->Write();
h_mX_SR_purity5->Write();
h_mX_SR_purity0->Write();
h_mX_SR_purity1->Write();
h_mX_SR_purity2->Write();
h_mX_SR_purity3->Write();
h_mX_SR_purity4->Write();
//h_mX_SR_kinFit->Write();
}
h_HH_balance_kinFit->Write();
if(type!="Data"){ h_mX_SR_kinFit_purity0->Write();
h_mX_SR_kinFit_purity1->Write();
h_mX_SR_kinFit_purity2->Write();
h_mX_SR_kinFit_purity3->Write();
h_mX_SR_kinFit_purity4->Write();
h_mX_SR_kinFit_purity5->Write(); }
h_mX_SB->Write();
h_mX_SB_biasCorrected->Write();
h_Cuts.Write();
h_mX_SB_kinFit->Write();
h_mX_SR_kinFit->Write();
tFile2->Write();
tFile2->Close();
std::cout<<"Wrote output file "<<histfilename<<std::endl;
std::cout<<"=== Cut Efficiencies === "<<std::endl;
std::cout<<"Number of events after finding HH candidate (btag && pT>40 GeV && |eta|<2.5) = "<<nCut4<<std::endl;
std::cout<<"Number of events in SR = "<<nCut5<<std::endl;
std::cout<<"========================"<<std::endl;
delete h_H1_mass;
delete h_H1_pT;
delete h_H2_mass;
delete h_H2_pT;
delete h_HH_balance;
delete h_mH1_mH2_asym;
delete h_H1_mass_biasCorrected;
delete h_H1_pT_biasCorrected;
delete h_H2_mass_biasCorrected;
delete h_H2_pT_biasCorrected;
delete h_HH_balance_biasCorrected;
delete h_mH1_mH2_asym_biasCorrected;
delete h_GenX_pT;
delete h_kinFitchi2;
delete h_chi;
delete h_chi_biasCorrected;
delete h_mX_SR;
delete h_mX_SR_biasCorrected;
delete h_mX_SR_purity5;
delete h_mX_SR_purity0;
delete h_mX_SR_purity1;
delete h_mX_SR_purity2;
delete h_mX_SR_purity3;
delete h_mX_SR_purity4;
delete h_mX_SR_kinFit;
delete h_HH_balance_kinFit;
delete h_mX_SR_kinFit_purity0;
delete h_mX_SR_kinFit_purity1;
delete h_mX_SR_kinFit_purity2;
delete h_mX_SR_kinFit_purity3;
delete h_mX_SR_kinFit_purity4;
delete h_mX_SR_kinFit_purity5;
delete h_mX_SB;
delete h_mX_SB_biasCorrected;
delete h_mX_SB_kinFit;
}