MA5SandBox: atlas_susy_2018_04.cpp

File atlas_susy_2018_04.cpp, 12.5 KB (added by Benjamin Fuks, 5 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/atlas_susy_2018_04.h"
2#include <iostream>
3#include <ctime>
4using namespace MA5;
5
6// User defined function: overlap Removal
7template<typename T1, typename T2> std::vector<const T1*>
8 Removal(std::vector<const T1*> &, std::vector<const T2*> &, const double &);
9
10// -----------------------------------------------------------------------------
11// Initialize
12// function called one time at the beginning of the analysis
13// -----------------------------------------------------------------------------
14bool atlas_susy_2018_04::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
15{
16 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
17 INFO << " <> Analysis: ATLAS SUSY 2018 04 <>" << endmsg;
18 INFO << " <> SUSY, stau, hadronic decaying ditau + MET @ 13 TeV, 139 fb^-1 <>" << endmsg;
19 INFO << " <> arXiv:1911.06660 <>" << endmsg;
20 INFO << " <> Recasted by : Jongwon Lim, Chih-Ting Lu, <>" << endmsg;
21 INFO << " <> Jae-hyeon Park, Jiwon Park <>" << endmsg;
22 INFO << " <> Contact : jongwon.lim@cern.ch, timluyu@gmail.com, <>" << endmsg;
23 INFO << " <> jhpark@kias.re.kr, jiwon.park@cern.ch <>" << endmsg;
24 INFO << " <> Based on MadAnalysis 5 v1.8 and above <>" << endmsg;
25 INFO << " <> For more information, see <>" << endmsg;
26 INFO << " <> http://madanalysis.irmp.ucl.ac.be/wiki/PublicAnalysisDatabase <>" << endmsg;
27 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
28
29 // ========================= //
30 // ===== Signal region ===== //
31 // ========================= //
32
33 Manager()->AddRegionSelection("SRlow");
34 Manager()->AddRegionSelection("SRhigh");
35
36 // ====================== //
37 // ===== Selections ===== //
38 // ====================== //
39
40 // Baseline cut
41 Manager()->AddCut("Baseline cut");
42
43 // Trigger and offline cuts
44 Manager()->AddCut("asymmetric di-$\\tau$ trigger", "SRlow");
45 Manager()->AddCut("di-$\\tau +E^{miss}_{T}$ trigger", "SRhigh");
46
47 // Common selections
48 Manager()->AddCut("2 medium $\\tau$ (OS) and 3rd medium $\\tau$ veto");
49 Manager()->AddCut("b-jet veto");
50 Manager()->AddCut("light lepton veto");
51 Manager()->AddCut("Z/H veto");
52
53 // SR Low Mass selections
54 Manager()->AddCut("$75 < E^{miss}_{T} < 150$ GeV", "SRlow");
55 Manager()->AddCut("2 tight $\\tau$ (OS)", "SRlow");
56 Manager()->AddCut("$|\\Delta\\phi(\\tau_{1},\\tau_{2})|>0.8$ [rad],low", "SRlow");
57 Manager()->AddCut("$\\Delta R(\\tau_{1},\\tau_{2})<3.2$,low", "SRlow");
58 Manager()->AddCut("$m_{T2}>70$ GeV,low", "SRlow");
59
60 // SR High Mass selections
61 Manager()->AddCut("$\\geq 1$ tight $\\tau$", "SRhigh");
62 Manager()->AddCut("$|\\Delta\\phi(\\tau_{1},\\tau_{2})|>0.8$ [rad],high", "SRhigh");
63 Manager()->AddCut("$\\Delta R(\\tau_{1},\\tau_{2})<3.2$,high", "SRhigh");
64 Manager()->AddCut("$m_{T2}>70$ GeV,high", "SRhigh");
65
66 // ====================== //
67 // ===== Histograms ===== //
68 // ====================== //
69 Manager()->AddHisto("SRlow_mT2", 5,70.0,120., "SRlow");
70 Manager()->AddHisto("SRhigh_mT2", 5,70.0,220., "SRhigh");
71
72 return true;
73}
74
75// -----------------------------------------------------------------------------
76// Finalize
77// function called one time at the end of the analysis
78// -----------------------------------------------------------------------------
79void atlas_susy_2018_04::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files){}
80
81// -----------------------------------------------------------------------------
82// Execute
83// function called each time one event is read
84// -----------------------------------------------------------------------------
85bool atlas_susy_2018_04::Execute(SampleFormat& sample, const EventFormat& event)
86{
87 // Event weight
88 double myWeight;
89 if ( Configuration().IsNoEventWeight() ) myWeight=1.;
90 else if ( event.mc()->weight()!=0. ) myWeight=event.mc()->weight();
91 else return false;
92 myWeight=1.;
93 Manager()->InitializeForNewEvent(myWeight);
94
95 // Security for empty events
96 if ( event.rec()==0 ) return true;
97
98 // ================================ //
99 // ===== Event reconstruction ===== //
100 // ================================ //
101
102 //// Jets ////
103 std::vector <const RecJetFormat*> Jets;
104 unsigned int nb = 0;
105 for( unsigned int ij=0; ij<event.rec()->jets().size(); ij++ ){
106 const RecJetFormat *Jet = &(event.rec()->jets()[ij]);
107
108 if( Jet->pt() > 20. && Jet->abseta() < 2.8 ){
109 Jets.push_back(Jet);
110 if( Jet->btag() && Jet->abseta() < 2.5) nb++;
111 }
112 }
113
114 //// Electrons ////
115 std::vector<const RecLeptonFormat*> SignalElectrons;
116 for( unsigned int ie=0; ie<event.rec()->electrons().size(); ie++ ){
117 const RecLeptonFormat *Lep = &(event.rec()->electrons()[ie]);
118
119 // Kinematics
120 double pt = Lep->pt();
121
122 // Isolation
123 double iso_dR = std::min(10./pt,0.2);
124 double iso_tracks = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(), iso_dR, 0., IsolationEFlow::TRACK_COMPONENT);
125 double iso_all = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(), 0.2, 0., IsolationEFlow::ALL_COMPONENTS);
126 bool iso = (iso_tracks<0.15 && iso_all<0.20);
127 if( pt>200. ) iso = (iso_all<std::max(0.015, 3.5/pt));
128
129 // Signal electrons
130 if( Lep->abseta() < 2.47 && pt > 17. && iso) SignalElectrons.push_back(Lep);
131 }
132
133
134 //// Muons ////
135 std::vector<const RecLeptonFormat*> SignalMuons;
136 for( unsigned int im=0; im<event.rec()->muons().size(); im++ ){
137 const RecLeptonFormat *Lep = &(event.rec()->muons()[im]);
138
139 // Kinematics
140 double pt = Lep->pt();
141
142 // Isolation
143 double iso_dR = std::min(10./pt,0.3);
144 double iso_tracks = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(), iso_dR, 0., IsolationEFlow::TRACK_COMPONENT);
145 double iso_all = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(), 0.2, 0., IsolationEFlow::ALL_COMPONENTS);
146 bool iso = (iso_tracks<0.15 && iso_all<0.30);
147
148 // Signal leptons
149 if( Lep->abseta() < 2.7 && pt > 14. && iso ) SignalMuons.push_back(Lep);
150 }
151
152
153 //// Taus ////
154 std::vector<const RecTauFormat*> SignalTaus;
155 for( unsigned int it=0; it<event.rec()->taus().size(); it++ ){
156 const RecTauFormat * CurrentTau = &(event.rec()->taus()[it]);
157
158 // Kinematics
159 double pt = CurrentTau->pt();
160 double eta = CurrentTau->abseta();
161
162 if( pt > 20. && eta < 2.5 && (eta < 1.37 || eta > 1.52) ){
163 SignalTaus.push_back(CurrentTau);
164 }
165 }
166
167
168 //// MET ////
169 MALorentzVector pTmiss = event.rec()->MET().momentum();
170 pTmiss.SetPz(0.); // Remove eta component
171 double MET = pTmiss.Pt();
172
173
174 // Cleaning of the jet collection and overlap removal
175 SignalTaus = Removal(SignalTaus, SignalElectrons, 0.2);
176 SignalTaus = Removal(SignalTaus, SignalMuons, 0.2);
177 SignalElectrons = Removal(SignalElectrons, SignalMuons, 0.01);
178 Jets = PHYSICS->Isol->JetCleaning(Jets, SignalElectrons, 0.2);
179 Jets = PHYSICS->Isol->JetCleaning(Jets, SignalMuons, 0.2);
180 SignalElectrons = Removal(SignalElectrons, Jets, 0.4);
181 SignalMuons = Removal(SignalMuons, Jets, 0.4);
182 Jets = Removal(Jets, SignalTaus, 0.2);
183
184
185 // ============================ //
186 // ===== Event selection ===== //
187 // =========================== //
188
189 // Select '2018' events according to the lumi weight
190 // 139 /fb = 36.2 + 44.3 + 58.5
191 event_num++;
192 srand(event_num);
193 bool is2018 = (rand() % 100)/100. < 0.421;
194
195 //// Baseline cut ////
196 if( !Manager()->ApplyCut(SignalTaus.size() == 2 && SignalTaus[0]->pt() > 50. && SignalTaus[1]->pt() > 40.,"Baseline cut") ) return true;
197
198 //// SRlow cut-1 : asymmetric di-tau trigger. ////a
199 // HLT online eff for tau candidates identified by the offline medium tau identification ~ 0.9,
200 // assuming that trigger object matching is done. ref: ATLAS-CONF-2017-029 Figure 13
201 myWeight = myWeight * 0.8 * 0.9 * 0.9;
202 Manager()->SetCurrentEventWeight(myWeight);
203 if( is2018 ){
204 if( !Manager()->ApplyCut(SignalTaus[0]->pt() > 95. && SignalTaus[1]->pt() > 75.,"asymmetric di-$\\tau$ trigger") )
205 return true;
206 }
207 else{
208 if( !Manager()->ApplyCut(SignalTaus[0]->pt() > 95. && SignalTaus[1]->pt() > 60.,"asymmetric di-$\\tau$ trigger") )
209 return true;
210 }
211
212 //// SRhigh cut-1 : di-tau +mET trigger. ////
213 if( is2018 ){
214 if( !Manager()->ApplyCut(SignalTaus[0]->pt() > 75. && MET > 150., "di-$\\tau +E^{miss}_{T}$ trigger") )
215 return true;
216 }
217 else{
218 if( !Manager()->ApplyCut(MET > 150., "di-$\\tau +E^{miss}_{T}$ trigger") )
219 return true;
220 }
221
222 // Since we don't know the efficiency of '2 taus to be
223 // medium tagged when 2 taus passed ditau(+met) trigger,
224 // it is implemented via reweighting.
225 // 0.7 is the ratio between MA5 and ATLAS efficiency
226 // from Nraw of the 2 medium tau cut, with OS selection
227 myWeight = myWeight * 0.7;
228 Manager()->SetCurrentEventWeight(myWeight);
229 //// Common cut-1 : 2 medium taus (OS) and 3rd medium tau veto. ////
230 if( !Manager()->ApplyCut(SignalTaus[0]->charge()*SignalTaus[1]->charge() < 0., "2 medium $\\tau$ (OS) and 3rd medium $\\tau$ veto") )
231 return true;
232
233 //// Common cut-2 : b-jet veto. ////
234 if( !Manager()->ApplyCut(nb == 0,"b-jet veto") ) return true;
235
236 //// Common cut-3 : light lepton veto. ////
237 if( !Manager()->ApplyCut(SignalElectrons.size() == 0 && SignalMuons.size() == 0, "light lepton veto") ) return true;
238
239 //// Common cut-4 : Z/H veto. ////
240 double mtata=0;
241 mtata=(SignalTaus[1]->momentum()+SignalTaus[0]->momentum()).M();
242 if( !Manager()->ApplyCut(mtata > 120.0, "Z/H veto") ) return true;
243
244 //// SRlow cut-2 : 75 < ETmiss < 150 GeV. ////
245 if( !Manager()->ApplyCut((MET > 75. && MET < 150.), "$75 < E^{miss}_{T} < 150$ GeV") ) return true;
246
247 //// SRlow cut-3 : 2 tight taus (OS) - implementation via reweighting ////
248 // 0.7 is ratio of N_weighted before and after 2-tight-tau cut in SR-lowMass at
249 // https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2018-04/tabaux_01.pdf
250 double tight_low=0.7;
251 Manager()->SetCurrentEventWeight(myWeight * tight_low);
252 if( !Manager()->ApplyCut(true, "2 tight $\\tau$ (OS)") ) return true;
253
254 //// SRlow cut-4 : |dphi(ta1,ta2)|>0.8 [rad]. ////
255 double DeltaPhiTau = SignalTaus[0]->dphi_0_pi(SignalTaus[1]);
256 if( !Manager()->ApplyCut(DeltaPhiTau > 0.8, "$|\\Delta\\phi(\\tau_{1},\\tau_{2})|>0.8$ [rad],low") ) return true;
257
258 //// SRlow cut-5 : dR(ta1,ta2)<3.2. ////
259 if( !Manager()->ApplyCut(SignalTaus[0]->dr(SignalTaus[1])<3.2, "$\\Delta R(\\tau_{1},\\tau_{2})<3.2$,low"))
260 return true;
261
262 //// SRlow cut-6 : mT2>70 GeV. ////
263 double mt2_high = PHYSICS->Transverse->MT2(SignalTaus[0],SignalTaus[1],event.rec()->MET(),0.);
264 if( !Manager()->ApplyCut(mt2_high > 70, "$m_{T2}>70$ GeV,low") ) return true;
265
266 // Histograms
267 Manager()->FillHisto("SRlow_mT2",mt2_high);
268
269
270 //// SRhigh cut-3 : >= 1 tight tau. ////
271 //eff of 2 tight tau = 0.7 = p, from 120 GeV mass point
272 //eff of 1 tight tau = sqrt(p)
273 //eff of at least 1 tau = p^2 + 2*(1-p)*p = 0.91
274 Manager()->SetCurrentEventWeight(myWeight * 0.91);
275 if( !Manager()->ApplyCut(true, "$\\geq 1$ tight $\\tau$") ) return true;
276
277 //// SRhigh cut-4 : |dphi(ta1,ta2)|>0.8 [rad]. ////
278 DeltaPhiTau = SignalTaus[0]->dphi_0_pi(SignalTaus[1]);
279 if( !Manager()->ApplyCut(DeltaPhiTau > 0.8, "$|\\Delta\\phi(\\tau_{1},\\tau_{2})|>0.8$ [rad],high") ) return true;
280
281 //// SRhigh cut-5 : dR(ta1,ta2)<3.2. ////
282 if( !Manager()->ApplyCut(SignalTaus[0]->dr(SignalTaus[1])< 3.2, "$\\Delta R(\\tau_{1},\\tau_{2})<3.2$,high"))
283 return true;
284
285 //// SRhigh cut-6 : mT2>70 GeV. ////
286 if( !Manager()->ApplyCut(mt2_high > 70, "$m_{T2}>70$ GeV,high") ) return true;
287
288 // Histograms
289 Manager()->FillHisto("SRhigh_mT2",mt2_high);
290
291 return true;
292}
293
294// User-defined functions for overlap removal
295template<typename T1, typename T2> std::vector<const T1*>
296 Removal(std::vector<const T1*> &v1, std::vector<const T2*> &v2,
297 const double &drmin)
298{
299 // Determining which objects should be removed
300 std::vector<bool> mask(v1.size(),false);
301 for (unsigned int j=0;j<v1.size();j++)
302 for (unsigned int i=0;i<v2.size();i++)
303 if (v2[i]->dr(v1[j]) < drmin)
304 {
305 mask[j]=true;
306 break;
307 }
308
309 // Building the cleaned container
310 std::vector<const T1*> cleaned_v1;
311 for (unsigned int i=0;i<v1.size();i++)
312 if (!mask[i]) cleaned_v1.push_back(v1[i]);
313
314 return cleaned_v1;
315}