Korea2017-04-CMS-SUS-16-041: Recasting_SUS16041-1.cpp

File Recasting_SUS16041-1.cpp, 24.2 KB (added by Chalons, 8 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/Recasting_SUS16041.h"
2using namespace MA5;
3using namespace std;
4
5// -- user-defined functions -- //
6double Calc_MiniIso( const RecLeptonFormat* lepton, const EventFormat& event );
7double Calc_PtRatio( const RecLeptonFormat* lepton, const EventFormat& event );
8double Calc_PtRel( const RecLeptonFormat* lepton, const EventFormat& event );
9std::vector<const RecJetFormat*> Find_MatchedJets( const RecLeptonFormat* lepton, const EventFormat& event );
10
11bool Find_DiLeptonPair( std::vector<const RecLeptonFormat*> &baseLeptons, vector<double> &vec_M, vector<double> &vec_MT, const EventFormat& event );
12double Calc_MT( const RecLeptonFormat* lepton, const EventFormat& event );
13
14double Check_FlavorComposition( vector<const RecLeptonFormat*> baseLeptons );
15
16
17// -----------------------------------------------------------------------------
18// Initialize
19// function called one time at the beginning of the analysis
20// -----------------------------------------------------------------------------
21bool Recasting_SUS16041::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
22{
23 cout << "BEGIN Initialization" << endl;
24
25 nEvent = 0;
26 SumWeight = 0;
27 // Flag_Verbose = true;
28 Flag_Verbose = false;
29 // -- define signal regions -- //
30 Manager()->AddRegionSelection("SR15b_onZ_HT600toInf_MET150to300_MT120toInf");
31 Manager()->AddRegionSelection("SR15b_offZ_HT600toInf_MET150to300_MT120toInf");
32
33 Manager()->AddRegionSelection("SR16a_onZ_MET300toInf_MT0to120");
34 Manager()->AddRegionSelection("SR16a_offZ_MET300toInf_MT0to120");
35
36 Manager()->AddRegionSelection("SR16b_onZ_MET300toInf_MT120toInf");
37 Manager()->AddRegionSelection("SR16b_offZ_MET300toInf_MT120toInf");
38
39 Manager()->AddRegionSelection("SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf");
40 Manager()->AddRegionSelection("SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf");
41
42 Manager()->AddRegionSelection("SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf");
43 Manager()->AddRegionSelection("SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf");
44
45 std::string SR_All[] =
46 {
47 "SR15b_onZ_HT600toInf_MET150to300_MT120toInf",
48 "SR15b_offZ_HT600toInf_MET150to300_MT120toInf",
49 "SR16a_onZ_MET300toInf_MT0to120",
50 "SR16a_offZ_MET300toInf_MT0to120",
51 "SR16b_onZ_MET300toInf_MT120toInf",
52 "SR16b_offZ_MET300toInf_MT120toInf",
53 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
54 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
55 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
56 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
57 };
58
59 std::string SR_onZ[] =
60 {
61 "SR15b_onZ_HT600toInf_MET150to300_MT120toInf",
62 "SR16a_onZ_MET300toInf_MT0to120",
63 "SR16b_onZ_MET300toInf_MT120toInf",
64 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
65 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
66 };
67
68 std::string SR_offZ[] =
69 {
70 "SR15b_offZ_HT600toInf_MET150to300_MT120toInf",
71 "SR16a_offZ_MET300toInf_MT0to120",
72 "SR16b_offZ_MET300toInf_MT120toInf",
73 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
74 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
75 };
76
77 std::string SR_nbjet0to2[] =
78 {
79 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
80 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf"
81 };
82
83 std::string SR_nbjet3toInf[] =
84 {
85 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
86 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
87 };
88
89 std::string SR_MET150to300[] =
90 {
91 "SR15b_onZ_HT600toInf_MET150to300_MT120toInf",
92 "SR15b_offZ_HT600toInf_MET150to300_MT120toInf"
93 };
94
95 std::string SR_MET300toInf[] =
96 {
97 "SR16a_onZ_MET300toInf_MT0to120",
98 "SR16a_offZ_MET300toInf_MT0to120",
99 "SR16b_onZ_MET300toInf_MT120toInf",
100 "SR16b_offZ_MET300toInf_MT120toInf"
101 };
102
103 std::string SR_MET250toInf[] =
104 {
105 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
106 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf"
107 };
108
109 std::string SR_MET50toInf[] =
110 {
111 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
112 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
113 };
114
115 std::string SR_HT600toInf[] =
116 {
117 "SR15b_onZ_HT600toInf_MET150to300_MT120toInf",
118 "SR15b_offZ_HT600toInf_MET150to300_MT120toInf"
119 };
120
121 std::string SR_HT200toInf[] =
122 {
123 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
124 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf"
125 };
126
127 std::string SR_HT60toInf[] =
128 {
129 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
130 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf"
131 };
132
133 std::string SR_MT0to120[] =
134 {
135 "SR16a_onZ_MET300toInf_MT0to120",
136 "SR16a_offZ_MET300toInf_MT0to120",
137 };
138
139 std::string SR_MT120toInf[] =
140 {
141 "SR15b_onZ_HT600toInf_MET150to300_MT120toInf",
142 "SR15b_offZ_HT600toInf_MET150to300_MT120toInf",
143 "SR16b_onZ_MET300toInf_MT120toInf",
144 "SR16b_offZ_MET300toInf_MT120toInf",
145 "SSR1_offZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
146 "SSR2_offZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
147 "SSR3_onZ_nbjet0to2_HT200toInf_MET_250toInf_MT120toInf",
148 "SSR4_onZ_nbjet3toInf_HT60toInf_MET_50toInf_MT120toInf",
149 };
150
151 Manager()->AddCut( "nGoodLepton >= 3", SR_All );
152 Manager()->AddCut( "nJet >= 2", SR_All );
153 Manager()->AddCut( "MET > 50 (70) GeV", SR_All );
154
155 Manager()->AddCut("onZ", SR_onZ);
156 Manager()->AddCut("offZ", SR_offZ);
157
158 Manager()->AddCut("nbjet <= 2", SR_nbjet0to2);
159 Manager()->AddCut("nbjet >= 3", SR_nbjet3toInf);
160
161 Manager()->AddCut("150 < MET < 300 GeV", SR_MET150to300 );
162 Manager()->AddCut("MET > 300 GeV", SR_MET300toInf );
163 Manager()->AddCut("MET > 250 GeV", SR_MET250toInf );
164 Manager()->AddCut("MET > 50 GeV", SR_MET50toInf );
165
166 Manager()->AddCut("HT > 600 GeV", SR_HT600toInf );
167 Manager()->AddCut("HT > 200 GeV", SR_HT200toInf );
168 Manager()->AddCut("HT > 60 GeV", SR_HT60toInf );
169
170 Manager()->AddCut("MT < 120 GeV", SR_MT0to120 );
171 Manager()->AddCut("MT > 120 GeV", SR_MT120toInf );
172
173
174 // -- histogram: for validation -- //
175 Manager()->AddHisto("MiniIso_Elec", 1000, 0, 100 );
176 Manager()->AddHisto("PtRatio_Elec", 1000, 0, 100 );
177 Manager()->AddHisto("PtRel_Elec", 1000, 0, 100 );
178 Manager()->AddHisto("MiniIso_Mu", 1000, 0, 100 );
179 Manager()->AddHisto("PtRatio_Mu", 1000, 0, 100 );
180 Manager()->AddHisto("PtRel_Mu", 1000, 0, 100 );
181
182 // -- histograms: after baseline selection -- //
183 Manager()->AddHisto("# jets (on-Z)", 6, 2, 8 );
184 Manager()->AddHisto("# b-jets (on-Z)", 5, 0, 5 );
185 Manager()->AddHisto("HT (on-Z)", 9, 60, 660 );
186 Manager()->AddHisto("MT (on-Z)", 10, 0, 200 );
187 Manager()->AddHisto("MET (on-Z)", 18, 50, 500 );
188 Manager()->AddHisto("Leading lepton pT (on-Z)", 14, 20, 300 );
189 Manager()->AddHisto("Sub-leading lepton pT (on-Z)", 9, 10, 100 );
190 Manager()->AddHisto("Trailing lepton pT (on-Z)", 9, 10, 100 );
191 Manager()->AddHisto("Flavor composition (on-Z)", 5, 0, 5 );
192
193 Manager()->AddHisto("# jets (off-Z)", 6, 2, 8 );
194 Manager()->AddHisto("# b-jets (off-Z)", 5, 0, 5 );
195 Manager()->AddHisto("HT (off-Z)", 9, 60, 660 );
196 Manager()->AddHisto("MT (off-Z)", 10, 0, 200 );
197 Manager()->AddHisto("MET (off-Z)", 18, 50, 500 );
198 Manager()->AddHisto("Leading lepton pT (off-Z)", 14, 20, 300 );
199 Manager()->AddHisto("Sub-leading lepton pT (off-Z)", 9, 10, 100 );
200 Manager()->AddHisto("Trailing lepton pT (off-Z)", 9, 10, 100 );
201 Manager()->AddHisto("Flavor composition (off-Z)", 5, 0, 5 );
202
203 cout << "END Initialization" << endl;
204
205 return true;
206}
207
208// -----------------------------------------------------------------------------
209// Finalize
210// function called one time at the end of the analysis
211// -----------------------------------------------------------------------------
212void Recasting_SUS16041::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
213{
214 cout << "BEGIN Finalization" << endl;
215 // saving histos
216 cout << "# total events processed: " << nEvent << endl;
217 cout << "\tSum of weights: " << SumWeight << endl;
218
219 cout << "END Finalization" << endl;
220}
221
222// -----------------------------------------------------------------------------
223// Execute
224// function called each time one event is read
225// -----------------------------------------------------------------------------
226bool Recasting_SUS16041::Execute(SampleFormat& sample, const EventFormat& event)
227{
228 if(event.rec()!=0)
229 {
230 double myEventWeight;
231 if(Configuration().IsNoEventWeight()) myEventWeight=1.;
232 else if(event.mc()->weight()!=0.) myEventWeight = event.mc()->weight();
233 else
234 {
235 WARNING << "Found one event with a zero weight. Skipping..." << endmsg;
236 return false;
237 }
238
239 myEventWeight = 1.0; // -- do not use weight! ... we will use only LO events (no negative weights) -- //
240 Manager()->InitializeForNewEvent(myEventWeight);
241
242 // -- count # events & sum of weights -- //
243 nEvent++;
244 SumWeight += myEventWeight;
245
246 if( Flag_Verbose )
247 {
248 cout << "##### " << nEvent << "th event #####" << endl;
249 cout << "Sum of weights: " << SumWeight << endl;
250 cout << "\tmyEventWeight: " << myEventWeight << endl;
251 }
252
253 std::vector<const RecLeptonFormat*> baseElectrons, baseMuons, baseLeptons;
254 std::vector<const RecJetFormat*> baseJets;
255
256 // -- baseline jets (anti-kt, r=0.4 simulated @ Delphes) -- //
257 double HT = 0; // -- sum of jet pT for jets with pT > 30 GeV -- //
258 int nbjet = 0;
259 for(unsigned int ii=0; ii<event.rec()->jets().size(); ii++)
260 {
261 const RecJetFormat *CurrentJet = &(event.rec()->jets()[ii]);
262 double pt = CurrentJet->pt();
263 double abseta = fabs(CurrentJet->eta());
264
265 if(pt > 30.0 && abseta < 2.4)
266 baseJets.push_back(CurrentJet);
267
268 // -- HT -- //
269 if( pt > 30.0 ) HT = HT + pt;
270
271 // -- # b-tagged jets -- //
272 if( CurrentJet->btag() && pt > 25.0 && abseta < 2.4 ) nbjet++;
273 }
274
275 // -- base electron selection -- //
276 for(unsigned int ii=0; ii<event.rec()->electrons().size(); ii++)
277 {
278 const RecLeptonFormat *CurrentElectron = &(event.rec()->electrons()[ii]);
279 // double pt = CurrentElectron->pt();
280 double abseta = fabs(CurrentElectron->eta());
281 double MiniIso = Calc_MiniIso(CurrentElectron, event);
282 double PtRatio = Calc_PtRatio(CurrentElectron, event);
283 double PtRel = Calc_PtRel(CurrentElectron, event);
284
285 Manager()->FillHisto( "MiniIso_Elec", MiniIso );
286 Manager()->FillHisto( "PtRatio_Elec", PtRatio );
287 Manager()->FillHisto( "PtRel_Elec", PtRel );
288
289 if( Flag_Verbose )
290 cout << "[" << ii << "th electron] MiniIso: " << MiniIso << ", PtRatio: " << PtRatio << ", PtRel: " << PtRel << endl;
291
292 if( abseta < 2.5 && MiniIso < 0.12 && (PtRatio > 0.76 || PtRel > 7.2) )
293 {
294 baseElectrons.push_back(CurrentElectron);
295 baseLeptons.push_back(CurrentElectron);
296 }
297 }
298
299 // -- base muon selection -- //
300 for(unsigned int ii=0; ii<event.rec()->muons().size(); ii++)
301 {
302 const RecLeptonFormat *CurrentMuon = &(event.rec()->muons()[ii]);
303 // double pt = CurrentMuon->pt();
304 double abseta = fabs(CurrentMuon->eta());
305 double MiniIso = Calc_MiniIso(CurrentMuon, event);
306 double PtRatio = Calc_PtRatio(CurrentMuon, event);
307 double PtRel = Calc_PtRel(CurrentMuon, event);
308
309 Manager()->FillHisto( "MiniIso_Mu", MiniIso );
310 Manager()->FillHisto( "PtRatio_Mu", PtRatio );
311 Manager()->FillHisto( "PtRel_Mu", PtRel );
312
313 if( Flag_Verbose )
314 cout << "[" << ii << "th muon] MiniIso: " << MiniIso << ", PtRatio: " << PtRatio << ", PtRel: " << PtRel << endl;
315
316 if( abseta < 2.4 && MiniIso < 0.16 && (PtRatio > 0.69 || PtRel > 6.0) )
317 {
318 baseMuons.push_back(CurrentMuon);
319 baseLeptons.push_back(CurrentMuon);
320 }
321 }
322
323 // -- test whether this event contains at least 3 leptons & passing kinematic cuts & dilepton mass > 12 GeV -- //
324 bool Flag_PassAcc = false;
325 int nBaseLepton = (int)baseLeptons.size();
326 if( baseLeptons.size() >= 3 )
327 {
328 // -- test pt cuts -- //
329 bool Flag_PassPtCut = false;
330 SORTER->sort( baseLeptons );
331
332 double PtCut_1st = 0;
333 double PtCut_2nd = 0;
334 double PtCut_3rd = 0;
335 if( HT < 300 )
336 {
337 PtCut_1st = 25;
338
339 // -- sub-leading cut: depend on the lepton flavor -- //
340 if( baseLeptons[1]->isMuon() ) PtCut_2nd = 10;
341 if( baseLeptons[1]->isElectron() ) PtCut_2nd = 15;
342
343 PtCut_3rd = 10;
344 }
345 else if( HT > 300 )
346 {
347 //-- leading cut: depend on the lepton flavor -- //
348 if( baseLeptons[0]->isMuon() ) PtCut_1st = 10;
349 if( baseLeptons[0]->isElectron() ) PtCut_1st = 15;
350
351 //-- leading cut: depend on the lepton flavor -- //
352 if( baseLeptons[1]->isMuon() ) PtCut_2nd = 10;
353 if( baseLeptons[1]->isElectron() ) PtCut_2nd = 15;
354
355 PtCut_3rd = 10;
356 }
357
358 if( baseLeptons[0]->pt() > PtCut_1st &&
359 baseLeptons[1]->pt() > PtCut_2nd &&
360 baseLeptons[2]->pt() > PtCut_3rd ) Flag_PassPtCut = true;
361
362 // -- test the dilepton(OS, SF) mass of all possible pairs: M(ll) > 12 GeV -- //
363 bool Flag_OSMass = true; // -- it will be false if one pair below M < 12 GeV is found -- //
364
365 vector<double> vec_M;
366 vector<double> vec_MT_temp;
367 bool Flag_PairFound = Find_DiLeptonPair( baseLeptons, vec_M, vec_MT_temp, event );
368 if( Flag_PairFound )
369 {
370 for( const auto& M : vec_M )
371 {
372 if( M <= 12.0 ) // -- reject the event with dilepton pair mass less than 12 GeV -- //
373 {
374 Flag_OSMass = false;
375 break;
376 }
377 }
378 }
379
380 if( Flag_PassPtCut && Flag_OSMass ) Flag_PassAcc = true;
381 }
382
383 // -- test on-Z or off-Z -- //
384 bool Flag_onZ = false;
385 double MT = 0; // -- will be decided in on-Z test (this value is related to on-Z condition) -- //
386 vector<double> vec_M;
387 vector<double> vec_MT;
388 bool Flag_PairFound = Find_DiLeptonPair( baseLeptons, vec_M, vec_MT, event );
389
390 double MT_NotInv_Smallest = 999; // -- only used when on-Z case -- //
391 if( Flag_PairFound )
392 {
393 int nPair = (int)vec_M.size();
394 for(int i_pair=0; i_pair<nPair; i_pair++)
395 {
396 double M_pair = vec_M[i_pair];
397 double MT_NotInv_temp = vec_MT[i_pair];
398
399 // cout << i_pair << "th pair: M = " << M_pair << ", MT = " << MT_NotInv_temp << ", HT = " << HT << endl;
400
401 if( M_pair > 76 && M_pair < 106 )
402 {
403 Flag_onZ = true;
404 if( MT_NotInv_temp < MT_NotInv_Smallest )
405 MT_NotInv_Smallest = MT_NotInv_temp;
406 }
407 }
408 }
409
410 // -- assign MT of this event -- //
411 if( Flag_onZ )
412 MT = MT_NotInv_Smallest;
413 else // -- off-Z: take the smallest MT value -- //
414 {
415 double MT_Smallest = 999;
416 for(int i_lep=0; i_lep<nBaseLepton; i_lep++)
417 {
418 double MT_temp = Calc_MT(baseLeptons[i_lep], event);
419 if( MT_temp < MT_Smallest )
420 MT_Smallest = MT_temp;
421 }
422
423 MT = MT_Smallest;
424 }
425
426 if( Flag_Verbose )
427 cout << "MT: "<< MT << endl;
428
429 // -- apply cuts following cut flow order-- //
430 double MET = event.rec()->MET().pt();
431 unsigned int nJet = baseJets.size();
432
433 if( Flag_Verbose )
434 cout << "MET: "<< MET << endl;
435
436 if(!Manager()->ApplyCut(Flag_PassAcc, "nGoodLepton >= 3"))
437 return true;
438
439 if(!Manager()->ApplyCut(nJet >= 2, "nJet >= 2"))
440 return true;
441
442 double METCut = Flag_onZ ? 70 : 50;
443 if(!Manager()->ApplyCut(MET > METCut, "MET > 50 (70) GeV"))
444 return true;
445
446 if( Flag_onZ )
447 {
448 Manager()->FillHisto("# jets (on-Z)", nJet );
449 Manager()->FillHisto("# b-jets (on-Z)", nbjet );
450 Manager()->FillHisto("HT (on-Z)", HT );
451 Manager()->FillHisto("MT (on-Z)", MT );
452 Manager()->FillHisto("MET (on-Z)", MET );
453 Manager()->FillHisto("Leading lepton pT (on-Z)", baseLeptons[0]->pt() );
454 Manager()->FillHisto("Sub-leading lepton pT (on-Z)", baseLeptons[1]->pt() );
455 Manager()->FillHisto("Trailing lepton pT (on-Z)", baseLeptons[2]->pt() );
456 double Flavor = Check_FlavorComposition( baseLeptons );
457 Manager()->FillHisto("Flavor composition (on-Z)", 0.1 ); // -- 1st bin: total -- //
458 Manager()->FillHisto("Flavor composition (on-Z)", Flavor );
459 }
460 else // -- off-Z -- //
461 {
462 Manager()->FillHisto("# jets (off-Z)", nJet );
463 Manager()->FillHisto("# b-jets (off-Z)", nbjet );
464 Manager()->FillHisto("HT (off-Z)", HT );
465 Manager()->FillHisto("MT (off-Z)", MT );
466 Manager()->FillHisto("MET (off-Z)", MET );
467 Manager()->FillHisto("Leading lepton pT (off-Z)", baseLeptons[0]->pt() );
468 Manager()->FillHisto("Sub-leading lepton pT (off-Z)", baseLeptons[1]->pt() );
469 Manager()->FillHisto("Trailing lepton pT (off-Z)", baseLeptons[2]->pt() );
470 double Flavor = Check_FlavorComposition( baseLeptons );
471 Manager()->FillHisto("Flavor composition (off-Z)", 0.1 ); // -- 1st bin: total -- //
472 Manager()->FillHisto("Flavor composition (off-Z)", Flavor );
473 }
474
475 // -- signal region -- //
476 if(!Manager()->ApplyCut(Flag_onZ, "onZ"))
477 return true;
478
479 if(!Manager()->ApplyCut(!Flag_onZ, "offZ"))
480 return true;
481
482 if(!Manager()->ApplyCut(nbjet <= 2, "nbjet <= 2"))
483 return true;
484
485 if(!Manager()->ApplyCut(nbjet >= 3, "nbjet >= 3"))
486 return true;
487
488 if(!Manager()->ApplyCut(MET > 150.0 && MET < 300.0, "150 < MET < 300 GeV"))
489 return true;
490
491 if(!Manager()->ApplyCut(MET > 300.0, "MET > 300 GeV"))
492 return true;
493
494 if(!Manager()->ApplyCut(MET > 250.0, "MET > 250 GeV"))
495 return true;
496
497 if(!Manager()->ApplyCut(MET > 50.0, "MET > 50 GeV"))
498 return true;
499
500 if(!Manager()->ApplyCut(HT > 600.0, "HT > 600 GeV"))
501 return true;
502
503 if(!Manager()->ApplyCut(HT > 200.0, "HT > 200 GeV"))
504 return true;
505
506 if(!Manager()->ApplyCut(HT > 60.0, "HT > 60 GeV"))
507 return true;
508
509 if(!Manager()->ApplyCut(MT < 120.0, "MT < 120 GeV"))
510 return true;
511
512 if(!Manager()->ApplyCut(MT > 120.0, "MT > 120 GeV"))
513 return true;
514
515 // SORTER->sort( baseJets );
516 // SORTER->sort(sigElectrons);
517 // SORTER->sort(sigMuons);
518 // SORTER->sort(sigLeptons);
519 if( Flag_Verbose )
520 cout << "\n" << endl;
521 }
522
523 return true;
524}
525
526double Calc_MiniIso( const RecLeptonFormat* lepton, const EventFormat& event )
527{
528 double RelIso_Combined = 0;
529 double pt = lepton->pt();
530 double dR = 10.0 / std::min( std::max(pt, 50.0), 200.0 );
531
532 double ALLIso = PHYSICS->Isol->eflow->sumIsolation(lepton, event.rec(), dR, 0., IsolationEFlow::ALL_COMPONENTS);
533
534 // RelIso_Combined = (ChargedIso + NeutralIso + PhotonIso ) / pt;
535 RelIso_Combined = (ALLIso) / pt;
536
537 // if( RelIso_Combined > 1 )
538 // {
539 // double ChargedIso = PHYSICS->Isol->eflow->sumIsolation(lepton, event.rec(), dR, 0., IsolationEFlow::TRACK_COMPONENT);
540 // double NeutralIso = PHYSICS->Isol->eflow->sumIsolation(lepton, event.rec(), dR, 0., IsolationEFlow::NEUTRAL_COMPONENT);
541 // double PhotonIso = PHYSICS->Isol->eflow->sumIsolation(lepton, event.rec(), dR, 0., IsolationEFlow::PHOTON_COMPONENT);
542
543 // cout << "Lepton pT: " << pt << endl;
544 // cout << "dR: " << dR << endl;
545 // cout << "ChargedIso: " << ChargedIso << endl;
546 // cout << "NeutralIso: " << NeutralIso << endl;
547 // cout << "PhotonIso: " << PhotonIso << endl;
548 // cout << "ALLIso: "<< ALLIso << endl;
549 // cout << "RelIso_Combined: " << RelIso_Combined << endl;
550 // }
551
552
553
554 return RelIso_Combined;
555}
556
557double Calc_PtRatio( const RecLeptonFormat* lepton, const EventFormat& event )
558{
559 double PtRatio = 0;
560 std::vector<const RecJetFormat*> MatchedJets = Find_MatchedJets( lepton, event );
561
562 unsigned int nJets = MatchedJets.size();
563 if( nJets == 0 ) // -- there's no matched jets -- //
564 PtRatio = 1;
565
566 else if( nJets == 1 )
567 PtRatio = lepton->pt() / MatchedJets[0]->pt();
568
569 else // -- more than 1 matched jets -- //
570 {
571 SORTER->sort( MatchedJets );
572 PtRatio = lepton->pt() / MatchedJets[0]->pt(); // -- take the largest-pT jets -- //
573 }
574
575 return PtRatio;
576}
577
578double Calc_PtRel( const RecLeptonFormat* lepton, const EventFormat& event )
579{
580 double PtRel = 0;
581 std::vector<const RecJetFormat*> MatchedJets = Find_MatchedJets( lepton, event );
582 unsigned int nJets = MatchedJets.size();
583 if( nJets == 0 ) // -- there's no matched jets -- //
584 PtRel = 0;
585
586 else if( nJets == 1 )
587 {
588 MALorentzVector P_lep = lepton->momentum();
589 MALorentzVector P_jet = MatchedJets[0]->momentum();
590 MALorentzVector P_axis = P_jet - P_lep;
591 PtRel = P_axis.Dot(P_lep) / P_axis.Mag();
592 }
593
594 else
595 {
596 SORTER->sort( MatchedJets ); // -- take the jet with largest pT -- //
597 MALorentzVector P_lep = lepton->momentum();
598 MALorentzVector P_jet = MatchedJets[0]->momentum();
599 MALorentzVector P_axis = P_jet - P_lep;
600 PtRel = P_axis.Dot(P_lep) / P_axis.Mag();
601 }
602
603 return PtRel;
604}
605
606std::vector<const RecJetFormat*> Find_MatchedJets( const RecLeptonFormat* lepton, const EventFormat& event )
607{
608 std::vector<const RecJetFormat*> MatchedJets;
609
610 for(unsigned int ii=0; ii<event.rec()->jets().size(); ii++)
611 {
612 const RecJetFormat *CurrentJet = &(event.rec()->jets()[ii]);
613
614 double dR_lj = lepton->dr(CurrentJet);
615 if( dR_lj < 0.4 && CurrentJet->pt() > 5.0 )
616 MatchedJets.push_back( CurrentJet );
617 }
618
619 return MatchedJets;
620}
621
622bool Find_DiLeptonPair( std::vector<const RecLeptonFormat*> &baseLeptons, vector<double> &vec_M, vector<double> &vec_MT, const EventFormat& event )
623{
624 bool Flag_PairFound = false;
625 int nBaseLepton = (int)baseLeptons.size();
626 for(int i_lep=0; i_lep<nBaseLepton; i_lep++)
627 {
628 const RecLeptonFormat* lepton_ith = baseLeptons[i_lep];
629 for(int j_lep=i_lep+1; j_lep<nBaseLepton; j_lep++)
630 {
631 const RecLeptonFormat* lepton_jth = baseLeptons[j_lep];
632
633 bool Flag_SF = (lepton_ith->isMuon() && lepton_jth->isMuon()) || (lepton_ith->isElectron() && lepton_jth->isElectron());
634 bool Flag_OS = lepton_ith->charge() != lepton_jth->charge();
635
636 // cout << "Flag_SF: " << Flag_SF << ", Flag_OS: " << Flag_OS << endl;
637
638 if( Flag_SF && Flag_OS )
639 {
640 Flag_PairFound = true;
641 double M = ( lepton_ith->momentum() + lepton_jth->momentum() ).M();
642 vec_M.push_back( M );
643
644 double MT_Smallest = 999; // -- take the smallest one -- //
645 for(int k_lep=0; k_lep<nBaseLepton; k_lep++) // -- iteration over the other leptons not involved with Z mass -- //
646 {
647 if( (k_lep != i_lep) && (k_lep != j_lep) )
648 {
649 double MT_temp = Calc_MT(baseLeptons[k_lep], event);
650 if( MT_temp < MT_Smallest ) MT_Smallest = MT_temp;
651 }
652 }
653 vec_MT.push_back( MT_Smallest );
654 }
655 } // -- end of j-th lepton iteration -- //
656 } // -- end of i-th lepton iteration -- //
657
658 return Flag_PairFound;
659}
660
661double Calc_MT( const RecLeptonFormat* lepton, const EventFormat& event )
662{
663 double MT = 0;
664 double lep_pt = lepton->pt();
665 double lep_phi = lepton->phi();
666 double MET_pt = event.rec()->MET().pt();
667 double MET_phi = event.rec()->MET().phi();
668
669 MT = sqrt( 2*lep_pt*MET_pt*(1 - std::cos(lep_phi-MET_phi)) );
670 return MT;
671}
672
673double Check_FlavorComposition( vector<const RecLeptonFormat*> baseLeptons )
674{
675 if( baseLeptons.size() < 3 )
676 {
677 cout << "[Check_FlavorComposition] # leptons (" << baseLeptons.size() << ") should be larger or equal to 3!" <<endl;
678 return 999;
679 }
680
681 int nMuon = 0;
682 int nElectron = 0;
683 for(int i=0; i<3; i++)
684 {
685 if( baseLeptons[i]->isElectron() ) nElectron++;
686 if( baseLeptons[i]->isMuon() ) nMuon++;
687 }
688
689 if( nMuon == 3 ) return 1 + 0.1; // -- 2nd bin -- //
690 if( nMuon == 2 && nElectron == 1 ) return 2 + 0.1; // -- 3rd bin -- //
691 if( nMuon == 1 && nElectron == 2 ) return 3 + 0.1; // -- 4th bin -- //
692 if( nElectron == 3 ) return 4 + 0.1; // -- 5th bin -- //
693
694 return 999;
695}