MA5SandBox: atlas_susy_2019_08.cpp

File atlas_susy_2019_08.cpp, 24.7 KB (added by Benjamin Fuks, 6 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/atlas_susy_2019_08.h"
2using namespace MA5;
3using namespace std;
4
5
6
7double mT(MALorentzVector &v1)
8{
9 return sqrt(v1.M2() + v1.Perp2());
10}
11
12/*
13double mTcalc(MALorentzVector &v1, MALorentzVector &v2)
14{
15 // m_T^2(p_1, p_2) =& m_1^2 + m_2^2 + 2 m_T^{(1)} m_T^{(2)} - 2 |p_T^{(1)}| |p_T^{(2)}| \cos \phi_{12}
16
17
18
19 return sqrt(v1.M2()+v2.M2() + 2.0*mT(v1)*mT(v2) - 2.0*v1.Pt()*v2.Pt()*cos(v1.DeltaPhi(v2)));
20
21};
22*/
23
24// we only use mTcalc for leptonp and pTmiss so fudge it
25double mTcalc(MALorentzVector &v1, MALorentzVector &v2)
26{
27 // m_T^2(p_1, p_2) =& m_1^2 + m_2^2 + 2 m_T^{(1)} m_T^{(2)} - 2 |p_T^{(1)}| >
28 return sqrt(v1.M2() + 2.0*mT(v1)*v2.Pt() - 2.0*v1.Pt()*v2.Pt()*cos(v1.DeltaPhi(v2)));
29
30
31};
32
33
34
35
36
37double mCTcalc(MALorentzVector &v1, MALorentzVector &v2)
38{
39 // m_{CT}^2(p_1, p_2) =& m_1^2 + m_2^2 + 2 m_T^{(1)} m_T^{(2)} + 2 |p_T^{(1)}| |p_T^{(2)}| \cos \phi_{12}
40
41 return sqrt(v1.M2()+v2.M2() + 2.0*mT(v1)*mT(v2) + 2.0*v1.Pt()*v2.Pt()*cos(v1.DeltaPhi(v2)));
42
43};
44
45
46
47
48MAfloat64 sumET(const RecLeptonFormat* part,
49 const std::vector<RecTrackFormat>& towers,
50 const MAfloat64& DR)
51{
52 MAfloat64 sumET=0.;
53 MAuint32 counter=0;
54
55 // Loop over the tracks
56 for (MAuint32 i=0;i<towers.size();i++)
57 {
58 const RecTrackFormat& tower = towers[i];
59
60
61 // Cut on the DR
62 if (part->momentum().DeltaR(tower.momentum()) > DR) continue;
63
64 // Sum
65 sumET += tower.et();
66 counter++;
67 }
68
69 // return PT sum of towers in the cone
70 return sumET;
71};
72
73
74MAfloat64 sumET(const RecLeptonFormat* part,
75 const std::vector<RecTowerFormat>& towers,
76 const MAfloat64& DR)
77{
78 MAfloat64 sumET=0.;
79 MAuint32 counter=0;
80
81 // Loop over the tracks
82 for (MAuint32 i=0;i<towers.size();i++)
83 {
84 const RecTowerFormat& tower = towers[i];
85
86
87 // Cut on the DR
88 if (part->momentum().DeltaR(tower.momentum()) > DR) continue;
89
90 // Sum
91 sumET += tower.et();
92 counter++;
93 }
94
95 // return PT sum of towers in the cone
96 return sumET;
97};
98
99
100
101
102MAfloat64 sumET(const RecLeptonFormat* part,
103 const std::vector<RecParticleFormat>& towers,
104 const MAfloat64& DR)
105{
106 MAfloat64 sumET=0.;
107 MAuint32 counter=0;
108
109 // Loop over the tracks
110 for (MAuint32 i=0;i<towers.size();i++)
111 {
112 const RecParticleFormat& tower = towers[i];
113
114
115 // Cut on the DR
116 if (part->momentum().DeltaR(tower.momentum()) > DR) continue;
117
118 // Sum
119 sumET += tower.et();
120 counter++;
121 }
122
123 // return PT sum of towers in the cone
124 return sumET;
125};
126
127
128MAfloat64 sumETIsolation(const RecLeptonFormat* part, const RecEventFormat* event, const MAfloat64& DR)
129 {
130 if (part==0) return 0;
131 if (event==0) return 0;
132 MAfloat64 sum=0.;
133 sum += sumET(part,event->EFlowTracks(),DR);
134 sum += sumET(part,event->EFlowPhotons(),DR);
135 sum += sumET(part,event->EFlowNeutralHadrons(),DR);
136 sum -= part->et();
137 return sum;
138 }
139
140
141//PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(),iso_dR,0,IsolationEFlow::TRACK_COMPONENT);
142
143
144bool IsLooseTightLepton(const RecLeptonFormat* Lep, const EventFormat& event,
145 const double DeltaRp, const double DeltaRE, const double pratio, const double Eratio)
146{
147 return true;
148 // varconeXX for momentum, coneXX for transverse energy
149 double pt = Lep->pt();
150 if(pt < 1.0) // just to be sure ... cut is much stricter than this on pt anyway
151 {
152 return false;
153 };
154 double DRp = std::min(10.0/pt,DeltaRp);
155 //double imini = fabs((PHYSICS->Isol->eflow->sumIsolation(Lep,event.rec(),DRp,0.0,IsolationEFlow::ALL_COMPONENTS))/pt);
156 double imini = (PHYSICS->Isol->eflow->sumIsolation(Lep,event.rec(),DRp,1.0,IsolationEFlow::ALL_COMPONENTS))/pt;
157 //double imini = (PHYSICS->Isol->eflow->sumIsolation(Lep,event.rec(),DRp,0.0,IsolationEFlow::ALL_COMPONENTS))/pt;
158 if (imini > pratio)
159 {
160 return false;
161
162 };
163
164 double jmini = sumETIsolation(Lep,event.rec(),DeltaRE)/pt;
165
166 if (jmini > Eratio)
167 {
168 return false;
169
170 };
171
172
173 return true;
174}
175
176
177bool IsHighPTCaloOnly(const RecLeptonFormat* Lep, const EventFormat& event,
178 const double DeltaRE)
179{
180 return true;
181 // varconeXX for momentum, coneXX for transverse energy
182 double pt = Lep->pt();
183 if(pt < 1.0) // We only apply this for pt> 200, but here we're just being safe
184 {
185 return false;
186 };
187 //double DRp = std::min(10.0/pt,DeltaRp);
188 double ETratio=std::max(0.015*pt,3.5);
189
190 double jmini = sumETIsolation(Lep,event.rec(),DeltaRE);
191
192 if (jmini > ETratio)
193 {
194 return false;
195
196 };
197
198
199 return true;
200}
201
202
203
204
205
206// Overlap Removal
207
208
209template<typename T1, typename T2> std::vector<const T1*>
210 Removal(std::vector<const T1*> &v1, std::vector<const T2*> &v2,
211 const double &drmin)
212{
213 // Determining with objects should be removed
214
215 std::vector<bool> mask(v1.size(),false);
216
217 for (unsigned int j=0;j<v1.size();j++)
218 for (unsigned int i=0;i<v2.size();i++)
219 if (v2[i]->dr(v1[j]) < drmin)
220 {
221 mask[j]=true;
222 break;
223 }
224 // Building the cleaned container
225 std::vector<const T1*> cleaned_v1;
226 for (unsigned int i=0;i<v1.size();i++)
227 if (!mask[i]) cleaned_v1.push_back(v1[i]);
228
229 return cleaned_v1;
230}
231
232
233template<typename T1> std::vector<const T1*> SelfRemoval(std::vector<const T1*> &v1, const double &drmin)
234{
235 // Determining with objects should be removed
236
237 std::vector<bool> mask(v1.size(),false);
238
239 double tdr;
240 for (unsigned int j=0;j<v1.size();j++){
241 for (unsigned int i=j+1;i<v1.size();i++) {
242
243 tdr=v1[i]->dr(v1[j]);
244 if (tdr < drmin)
245 {
246 if(v1[i]->pt() > v1[j]->pt())
247 {
248 mask[j]=true;
249 }
250else
251{
252 mask[i]=true;
253}
254// break;
255 }
256
257 }
258 }
259
260
261 // Building the cleaned container
262 std::vector<const T1*> cleaned_v1;
263 for (unsigned int i=0;i<v1.size();i++)
264 if (!mask[i]) cleaned_v1.push_back(v1[i]);
265
266 return cleaned_v1;
267}
268
269
270//Remove "Lepton Jets" via a crude criterion
271template<typename T1, typename T2> std::vector<const T1*>
272 RemoveFakeJets(std::vector<const T1*> &v1, std::vector<const T2*> &v2)
273{
274 // somewhat like the jet cleaning except we choose a matching criterion based on the energies too
275 // Determining which objects should be removed
276 double drmin=0.05;
277 std::vector<bool> mask(v1.size(),false);
278 // cout << "Cleaning " << v1.size() << " jets" << endl;
279 for (unsigned int i=0;i<v2.size();i++) // loop over leptons
280 {
281 int jet_index=-1; // match each lepton to only one lepton jet
282 double dr=0.0;
283 double dr_closest=0.0;
284
285 for (unsigned int j=0;j<v1.size();j++) // loop over jets
286 {
287 if(!v1[j]->btag())
288 {
289 dr = v2[i]->dr(v1[j]);
290 /*
291 if(dr < 0.01){
292 cout << "Small dr! " << dr << endl;
293 }
294 */
295 if (dr < drmin)
296 {
297 double leptonE=v2[i]->e();
298 double jetE=v1[j]->e();
299
300 if(((leptonE > 0.4*jetE) && (leptonE < 1.2 * jetE)) || (dr < 0.05))
301 {
302 if (jet_index ==-1 || dr < dr_closest)
303 {
304 jet_index = j;
305 dr_closest=dr;
306 }
307
308 }
309 }
310 }
311 }
312
313 if(jet_index !=-1)
314 {
315 mask[jet_index] = true;
316 }
317
318
319 }
320 // Building the cleaned container
321 std::vector<const T1*> cleaned_v1;
322 for (unsigned int i=0;i<v1.size();i++)
323 if (!mask[i]) cleaned_v1.push_back(v1[i]);
324
325 return cleaned_v1;
326}
327
328// Remove leptons from around b-jets since they're probably fake, unless they are very energetic
329
330void CleanBJets(std::vector<const RecJetFormat*> &jets, std::vector<const RecLeptonFormat*> &leptons)
331{
332 std::vector<bool> maskLepton(leptons.size(),false);
333 std::vector<const RecLeptonFormat*> cleaned_Lepton;
334 double drmin = 0.2;
335 for(auto jet: jets)
336 {
337
338 if(!jet->btag()) continue;
339 double jetE=jet->e();
340 for(unsigned int i=0; i<leptons.size(); i++)
341 {
342
343 if(jet->dr(leptons[i]) < drmin)
344 {
345 double lepE=leptons[i]->e();
346 if(lepE < 1.2*jetE)
347 {
348 maskLepton[i] = true;
349 }
350 }
351
352 };
353
354 }
355
356 for(unsigned int i=0;i<leptons.size();i++)
357 {
358 if(!maskLepton[i]) {
359 cleaned_Lepton.push_back(leptons[i]);
360 }
361 }
362
363
364 leptons=cleaned_Lepton;
365}
366
367
368
369
370void JetOrLepton(std::vector<const RecJetFormat*> &v1, std::vector<const RecLeptonFormat*> &v2)
371{
372 // Determining with objects should be removed
373 double drmin=0.4;
374 std::vector<bool> maskJet(v1.size(),false);
375 std::vector<bool> maskLepton(v2.size(),false);
376 for (unsigned int j=0;j<v1.size();j++)
377 {
378
379 for (unsigned int i=0;i<v2.size();i++)
380 {
381 if (v2[i]->dr(v1[j]) < drmin)
382 {
383 double leptonE=v2[i]->e();
384 double jetE=v1[j]->e();
385 if(leptonE > 1.1*jetE) continue;
386 if((leptonE > 0.5*jetE) && ( leptonE < 1.1*jetE))
387 {
388 if(!v1[j]->btag())
389 {
390 maskJet[j]=true;
391 }
392
393 }
394 else
395 {
396
397
398 maskLepton[i]=true;
399
400 }
401 }
402 }
403 }
404 // Building the cleaned container
405 std::vector<const RecJetFormat*> cleaned_Jet;
406 std::vector<const RecLeptonFormat*> cleaned_Lepton;
407 for (unsigned int i=0;i<v1.size();i++)
408 if (!maskJet[i]) cleaned_Jet.push_back(v1[i]);
409
410 for (unsigned int i=0;i<v2.size();i++)
411 if (!maskLepton[i]) cleaned_Lepton.push_back(v2[i]);
412
413 // v1=cleaned_Jet;
414 v2=cleaned_Lepton;
415}
416
417
418
419
420// Overlap Removal for electrons with Delta R < min(0.4, 0.04 + 10 GeV/pT) around Jet
421template<typename T1, typename T2> std::vector<const T1*>
422 RemovalJE(std::vector<const T1*> &v1, std::vector<const T2*> &v2)
423{
424 // Determining with objects should be removed
425 double EPT,drmin;
426 std::vector<bool> mask(v1.size(),false);
427 for (unsigned int j=0;j<v1.size();j++)
428 {
429 EPT=v1[j]->pt();
430 drmin=0.04+10/EPT;
431 if(drmin > 0.4) drmin=0.4;
432
433 for (unsigned int i=0;i<v2.size();i++)
434 if (v2[i]->dr(v1[j]) < drmin)
435 {
436 mask[j]=true;
437 break;
438 }
439 }
440 // Building the cleaned container
441 std::vector<const T1*> cleaned_v1;
442 for (unsigned int i=0;i<v1.size();i++)
443 if (!mask[i]) cleaned_v1.push_back(v1[i]);
444
445 return cleaned_v1;
446}
447
448
449
450// -----------------------------------------------------------------------------
451// Initialize
452// function called one time at the beginning of the analysis
453// -----------------------------------------------------------------------------
454bool atlas_susy_2019_08::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
455{
456 // initialize variables, histos
457 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
458 INFO << " <> Analysis: ATLAS-SUSY-2019-08 <>" << endmsg;
459 INFO << " <> arXiv:1909.09226 <>" << endmsg;
460 INFO << " <> chargino + neutralino->W->lv+H->bb <>" << endmsg;
461 INFO << " <> Recasters: M. Goodsell <>" << endmsg;
462 INFO << " <> Contact: goodsell@lpthe.jussieu.fr <>" << endmsg;
463 INFO << " <> Based on MadAnalysis 5 v1.8 <>" << endmsg;
464 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
465
466 //engine(m_randomdevice());
467 this->engine = std::mt19937(time(NULL));
468 // Declaration of the signal regions (SR)
469
470
471 Manager()->AddRegionSelection("LMdisc");
472 Manager()->AddRegionSelection("MMdisc");
473 Manager()->AddRegionSelection("HMdisc");
474
475 Manager()->AddRegionSelection("LMlowCT");
476 Manager()->AddRegionSelection("LMmedCT");
477 Manager()->AddRegionSelection("LMhighCT");
478
479 Manager()->AddRegionSelection("MMlowCT");
480 Manager()->AddRegionSelection("MMmedCT");
481 Manager()->AddRegionSelection("MMhighCT");
482
483 Manager()->AddRegionSelection("HMlowCT");
484 Manager()->AddRegionSelection("HMmedCT");
485 Manager()->AddRegionSelection("HMhighCT");
486
487
488
489 //Preselection Cuts
490 //Manager()->AddCut("Preselection cuts");
491
492
493
494
495 std::string alldiscSRs[]={"LMdisc","MMdisc","HMdisc"};
496 std::string allSRs[]={"LMdisc","MMdisc","HMdisc","LMlowCT","LMmedCT","LMhighCT","MMlowCT","MMmedCT","MMhighCT","HMlowCT","HMmedCT","HMhighCT"};
497
498
499 std::string LMexcl[]={"LMlowCT","LMmedCT","LMhighCT"};
500 std::string allLM[]={"LMdisc","LMlowCT","LMmedCT","LMhighCT"};
501
502 std::string MMexcl[]={"MMlowCT","MMmedCT","MMhighCT"};
503 std::string allMM[]={"MMdisc","MMlowCT","MMmedCT","MMhighCT"};
504
505 std::string HMexcl[]={"HMlowCT","HMmedCT","HMhighCT"};
506 std::string allHM[]={"HMdisc","HMlowCT","HMmedCT","HMhighCT"};
507
508
509 std::string lowCT[]={"LMlowCT","MMlowCT","HMlowCT"};
510 std::string medCT[]={"LMmedCT","MMmedCT","HMmedCT"};
511 std::string highCT[]={"LMhighCT","MMhighCT","HMhighCT"};
512
513
514
515
516
517
518
519 // for cutflows, preselection cuts are split out:
520
521 //Manager()->AddCut("Njets25 >=2",allSRs2);
522
523
524 //Manager()->AddCut("1 signal lepton",allSRs);
525
526 Manager()->AddCut("Njets25 >=2",allSRs);
527 Manager()->AddCut("1 signal lepton",allSRs);
528 Manager()->AddCut("Second baseline lepton veto",allSRs);
529 Manager()->AddCut("mT>50",allSRs);
530 Manager()->AddCut("ET>180",allSRs);
531 Manager()->AddCut("Njets<=3",allSRs);
532 Manager()->AddCut("2 bjets",allSRs);
533 Manager()->AddCut("mbb>50",allSRs);
534 Manager()->AddCut("ET>240",allSRs);
535 Manager()->AddCut("mbb",allSRs);
536
537 Manager()->AddCut("mlb > 120",allHM);
538
539 Manager()->AddCut("mT in [100,160]",LMexcl);
540 Manager()->AddCut("mT in [160,240]",MMexcl);
541 Manager()->AddCut("mT > 100","LMdisc");
542 Manager()->AddCut("mT > 160","MMdisc");
543 Manager()->AddCut("mT > 240",allHM);
544
545
546
547 Manager()->AddCut("minCT",alldiscSRs);
548
549 Manager()->AddCut("mCT in [180,230]",lowCT);
550 Manager()->AddCut("mCT in [230,280]",medCT);
551 Manager()->AddCut("mCT > 280",highCT);
552
553
554
555
556 // Event counter
557 //Nevents = 0;
558 return true;
559}
560
561// -----------------------------------------------------------------------------
562// Finalize
563// function called one time at the end of the analysis
564// -----------------------------------------------------------------------------
565void atlas_susy_2019_08::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
566{
567 cout << "BEGIN Finalization" << endl;
568 // saving histos
569 cout << "END Finalization" << endl;
570}
571
572// -----------------------------------------------------------------------------
573// Execute
574// function called each time one event is read
575// -----------------------------------------------------------------------------
576bool atlas_susy_2019_08::Execute(SampleFormat& sample, const EventFormat& event)
577{
578
579 // Nevents += 1;
580
581 double myWeight=1.;
582 if (!Configuration().IsNoEventWeight()) myWeight=event.mc()->weight();
583 else if(event.mc()->weight()!=0.) myWeight=event.mc()->weight();
584 else
585 {
586 //////WARNING << "Found one event with a zero weight. Skipping...\n";
587 return false;
588 }
589
590
591 Manager()->InitializeForNewEvent(myWeight);
592
593 // The event loop starts here
594 if(event.rec()!=0)
595 {
596
597 // Containers
598 std::vector<const RecLeptonFormat*> InitialElectrons, InitialMuons, Electrons, Muons, BaselineLeptons, BaselineElectrons, BaselineMuons, SoftElectrons, SoftMuons;
599 std::vector<const RecJetFormat*> SoftJets, PSoftJets,SignalJets, bJets, Jets25;
600
601 std::vector<const RecParticleFormat*> InitialPhotons;
602 SoftElectrons.clear();
603 Electrons.clear();
604 SoftMuons.clear();
605 Muons.clear();
606 SoftJets.clear();
607 SignalJets.clear();
608 InitialPhotons.clear();
609
610 MALorentzVector pTmiss = event.rec()->MET().momentum();
611
612 MALorentzVector precleanpt;
613 for(unsigned int j=0; j<event.rec()->jets().size(); j++)
614 {
615 const RecJetFormat *CurrentJet = &(event.rec()->jets()[j]);
616 if(CurrentJet->pt()>20. && fabs(CurrentJet->eta())<4.5) // these are soft Jets
617 {
618
619 SoftJets.push_back(CurrentJet);
620
621 }
622 else if(fabs(CurrentJet->eta()) > 2.5)
623 {
624 pTmiss+=CurrentJet->momentum(); // not picked up in the "soft term" -> becomes missing pt
625 }
626
627
628 }
629
630 Electrons.clear();
631 Muons.clear();
632 InitialElectrons.clear();
633 InitialMuons.clear();
634
635 for(unsigned int e=0; e<event.rec()->electrons().size(); e++)
636 {
637 const RecLeptonFormat *CurrentElectron = &(event.rec()->electrons()[e]);
638 if(CurrentElectron->pt() > 7.0)
639 {
640 InitialElectrons.push_back(CurrentElectron);
641 }
642
643 };
644 for(unsigned int e=0; e<event.rec()->muons().size(); e++)
645 {
646 const RecLeptonFormat *CurrentMuon = &(event.rec()->muons()[e]);
647 if(CurrentMuon->pt() > 6.0)
648 {
649 InitialMuons.push_back(CurrentMuon);
650 }
651 }
652
653 for(unsigned int e=0; e<event.rec()->photons().size(); e++)
654 {
655 const RecParticleFormat *CurrentPhoton = &(event.rec()->photons()[e]);
656 if(CurrentPhoton->pt() > 15.0)
657 {
658 InitialPhotons.push_back(CurrentPhoton);
659 }
660 }
661
662 // Remove leptons/photons masquerading as jets
663
664 SoftJets=RemoveFakeJets(SoftJets,InitialElectrons);
665 SoftJets=RemoveFakeJets(SoftJets,InitialMuons);
666 SoftJets=RemoveFakeJets(SoftJets,InitialPhotons);
667
668 // now emulate removal of jets through JVT and pileup
669 std::uniform_real_distribution<double> rd(0.0,1.0);
670 for (auto tJet : SoftJets)
671 {
672 double teff=1.0;
673 if(fabs(tJet->eta()) < 2.8)
674 {
675 double jetpt=tJet->pt();
676 if(fabs(tJet->eta()) > 2.5) teff=0.8;
677 else if(jetpt < 120.0)
678 {
679 if(tJet->btag()) teff=0.9; // have already identified the jet through vertex info, just apply the official JVT efficiency here
680 else teff = 0.8;
681
682
683 }
684
685 }
686 if(rd(engine) > teff)
687 {
688 // don't include in MET calculation
689 pTmiss+=tJet->momentum();
690 continue;
691 };
692 PSoftJets.push_back(tJet);
693 precleanpt+=tJet->momentum();
694 }
695 SoftJets.clear();
696 SoftJets=PSoftJets;
697
698 const std::vector<RecTrackFormat> tracks=(event.rec()->tracks());
699 std::map<int,double> lepton_to_vert;
700
701 // It's difficult to distinguish hadronic leptons ("background") from signal ones using delphes, because
702 // we lose jet constituents. ATLAS certainly uses displacement of the track as a criterion,
703 // but not the only one. I just have the displacement of the vertex. This is very good
704 // for most cases, but misses a few coming from neutral pions, omegas, etc (not many). You could also argue that it
705 // is cheating because I set the displacement as 0.1 mm. But otherwise we get too many baseline leptons.
706 for(auto track: tracks)
707 {
708 int tracktag=track.delphesTags()[0];
709
710 const MCParticleFormat* mymc = track.mc();
711
712 double pos=mymc->decay_vertex().P();
713
714 lepton_to_vert[tracktag]=pos;
715
716 }
717
718
719 for(auto CurrentElectron : InitialElectrons)
720 {
721
722 int etag=CurrentElectron->delphesTags()[0];
723
724 double pos=lepton_to_vert[etag];
725 double sint=CurrentElectron->pt()/CurrentElectron->p();
726 double dzsint=fabs(CurrentElectron->dz()*sint);
727 //double d0=fabs(CurrentElectron->d0());
728
729 if((pos < 0.1) && (dzsint < 0.5) && (CurrentElectron->pt()>7.) && (fabs(CurrentElectron->eta())<2.47))
730 {
731 // Loose -> Baseline
732 // then check tight, tight + HighPTCaloOnly for pt > 200
733 // dRp, dRE, ratiop, ratio
734 // loose is 0.2, 0.2, 0.15, 0.2
735 // tight is 0.2, 0.2, 0.06, 0.06
736
737 if(IsLooseTightLepton(CurrentElectron, event,0.2, 0.2, 0.15, 0.2))
738 {
739 // at least it's loose
740 SoftElectrons.push_back(CurrentElectron);
741 precleanpt+=CurrentElectron->momentum();
742 }
743
744 }
745 }
746
747
748 for(auto CurrentMuon : InitialMuons)
749 {
750 int mtag=CurrentMuon->delphesTags()[0];
751
752 double pos=lepton_to_vert[mtag];
753
754 double sint=CurrentMuon->pt()/CurrentMuon->p();
755 double dzsint=fabs(CurrentMuon->dz()*sint);
756 //double d0=fabs(CurrentMuon->d0());
757
758 if((pos < 0.1) && (dzsint < 0.5) && (CurrentMuon->pt()>6.) && (fabs(CurrentMuon->eta())<2.7))
759 {
760
761 SoftMuons.push_back(CurrentMuon);
762 precleanpt+=CurrentMuon->momentum();
763 }
764
765 }
766
767 //////////////////////////
768 // cleaning as described in the paper
769 BaselineElectrons=SelfRemoval(SoftElectrons,0.01);
770 BaselineElectrons=Removal(BaselineElectrons,SoftMuons,0.01);
771 BaselineElectrons=RemovalJE(BaselineElectrons,SoftJets);
772 SoftJets=Removal(SoftJets,SoftMuons,0.2);
773 BaselineMuons=RemovalJE(SoftMuons,SoftJets);
774
775 //////////////////////////
776 // now we see how much momentum has been subtracted away by the cleaning process, while sorting into Baseline and Signal leptons
777 BaselineLeptons.clear();
778 for (auto tElec : BaselineElectrons)
779 {
780 precleanpt-=tElec->momentum();
781 // apply tight criterion for the vertex
782 //double d0sigma=0.001*(11.0 + 94.0/tElec->pt());
783 double d0sigma=0.01+0.04/tElec->pt();
784
785 bool d0ratio = (tElec->d0()/d0sigma < 5.0);
786
787 if(d0ratio && IsLooseTightLepton(tElec, event, 0.2, 0.2, 0.06, 0.06))
788 {
789
790 if((tElec->pt()<200.0) || IsHighPTCaloOnly(tElec, event,0.2))
791 {
792 Electrons.push_back(tElec);
793 }
794 else
795 {
796 BaselineLeptons.push_back(tElec);
797 }
798 }
799 else
800 {
801 BaselineLeptons.push_back(tElec);
802 }
803
804 }
805 for (auto tMuon : BaselineMuons)
806 {
807 precleanpt-=tMuon->momentum();
808 //double d0sigma=0.001*(11.0 + 94.0/tMuon->pt());
809 double d0sigma=0.01+0.04/tMuon->pt();
810 bool d0ratio = (tMuon->d0()/d0sigma < 3.0);
811
812 if((d0ratio) && (tMuon->pt()>6.) && (fabs(tMuon->eta())<2.5) && IsLooseTightLepton(tMuon, event, 0.3, 0.2, 0.15, 0.3))
813 {
814 Muons.push_back(tMuon);
815 }
816 else
817 {
818 BaselineLeptons.push_back(tMuon);
819 }
820 }
821
822
823 MALorentzVector bJetMomentum;
824
825 int nj25=0;
826
827 for (auto tJet : SoftJets)
828 {
829 precleanpt-=tJet->momentum();
830 if(tJet->pt() > 25.0)
831 {
832
833 nj25++;
834 }
835
836
837 if(tJet->pt()>30. && fabs(tJet->eta())<2.8)
838 {
839 SignalJets.push_back(tJet);
840
841
842 if((tJet->btag()) && (fabs(tJet->eta())< 2.5))
843 {
844 bJetMomentum+=tJet->momentum();
845 bJets.push_back(tJet);
846
847 }
848
849
850 }
851
852 }
853
854 SORTER->sort(SignalJets);
855 SORTER->sort(bJets);
856 pTmiss+=precleanpt; // this is the difference between pre and post cleaning procedures -> becomes missing energy
857 double MET = pTmiss.Pt();
858
859 unsigned int NJets = SignalJets.size();
860
861
862 //Preselection cuts
863
864
865
866 bool IsMET180= (MET>180.);
867 bool IsNJets2 = (nj25>1);
868 bool IsMET = (MET>240.);
869 int nleptons=Electrons.size()+Muons.size();
870
871
872 bool IsLeptons = (nleptons == 1);
873 bool IsNJets = (NJets < 4);
874 bool isBJets = (bJets.size() == 2);
875
876 unsigned int NbaselineL= BaselineLeptons.size();
877
878 bool NoBaselineLeptons= (NbaselineL == 0);
879
880 bool isFromHiggs = false;
881 bool mbb50 = false;
882 if(isBJets)
883 {
884 double mbb = bJetMomentum.M();
885 mbb50=(mbb>50.0);
886 isFromHiggs = (mbb > 100.) && (mbb < 140.);
887 }
888
889
890 if(!Manager()->ApplyCut(IsNJets2,"Njets25 >=2")) return true;
891 if(!Manager()->ApplyCut(IsLeptons,"1 signal lepton")) return true;
892 if(!Manager()->ApplyCut(NoBaselineLeptons,"Second baseline lepton veto")) return true;
893 MALorentzVector leptonp;
894 if (Electrons.size() == 1)
895 {
896 // only one, so lepton pt is electron pt
897 leptonp=Electrons[0]->momentum();
898 }
899 else
900 {
901 if (Muons.size() > 0)
902 {
903 leptonp=Muons[0]->momentum();
904 }
905 else
906 {
907 leptonp.clear();
908 }
909
910 }
911
912 // mT defined in the paper like this, only true for massless leptons ...
913 //double mT = sqrt(2.0*leptonp.Pt()*MET*(1.0-cos(leptonp.DeltaPhi(pTmiss))));
914
915 double mT = mTcalc(leptonp,pTmiss);
916 bool IsMT50 = (mT > 50.);
917 if(!Manager()->ApplyCut(IsMT50,"mT>50")) return true;
918 if(!Manager()->ApplyCut(IsMET180,"ET>180")) return true;
919 if(!Manager()->ApplyCut(IsNJets,"Njets<=3")) return true;
920 if(!Manager()->ApplyCut(isBJets,"2 bjets")) return true;
921 if(!Manager()->ApplyCut(mbb50,"mbb>50")) return true;
922 if(!Manager()->ApplyCut(IsMET,"ET>240")) return true;
923 if(!Manager()->ApplyCut(isFromHiggs,"mbb")) return true;
924
925
926 MALorentzVector bjetmom1,bjetmom2;
927
928 if(bJets.size() == 2)
929 {
930 bjetmom1=bJets[0]->momentum();
931 bjetmom2=bJets[1]->momentum();
932 }
933
934 // mCT defined in the paper like this, only true for massless jets. Pseudocode uses some pre-defined function. Who to believe?
935 //double mCT = sqrt(2.0*(bjetmom1.Pt())*(bjetmom2.Pt())*(1.0+cos(bjetmom1.DeltaPhi(bjetmom2))));
936
937 double mCT = mCTcalc(bjetmom1,bjetmom2);
938 double mlb = (bjetmom1+leptonp).M();
939
940 // cuts that depend on the SR
941
942 if(!Manager()->ApplyCut(mlb> 120.0,"mlb > 120")) return true;
943
944 if(!Manager()->ApplyCut((mT > 100.0) && (mT <= 160.0),"mT in [100,160]")) return true;
945 if(!Manager()->ApplyCut((mT > 160.0) && (mT <= 240.0),"mT in [160,240]")) return true;
946 if(!Manager()->ApplyCut(mT > 100.0,"mT > 100")) return true;
947 if(!Manager()->ApplyCut(mT > 160.0,"mT > 160")) return true;
948 if(!Manager()->ApplyCut(mT > 240.0,"mT > 240")) return true;
949
950
951 if(!Manager()->ApplyCut(mCT> 180.0,"minCT")) return true;
952 if(!Manager()->ApplyCut((mCT > 180.0) && (mCT <= 230.0),"mCT in [180,230]")) return true;
953 if(!Manager()->ApplyCut((mCT > 230.0) && (mCT <= 280.0),"mCT in [230,280]")) return true;
954 if(!Manager()->ApplyCut(mCT> 280.0,"mCT > 280")) return true;
955
956 }
957 return true;
958}
959