MA5SandBox: atlas_susy_2018_031.cpp

File atlas_susy_2018_031.cpp, 16.0 KB (added by Benjamin Fuks, 5 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/atlas_susy_2018_031.h"
2using namespace MA5;
3using namespace std;
4
5template<typename T1, typename T2> std::vector<const T1*>
6 Removal(std::vector<const T1*>&, std::vector<const T2*>&, const double &);
7// -----------------------------------------------------------------------------
8// Initialize
9// function called one time at the beginning of the analysis
10// -----------------------------------------------------------------------------
11bool atlas_susy_2018_031::Initialize(const MA5::Configuration& cfg,
12 const std::map<std::string,std::string>& parameters)
13{
14 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
15 INFO << " <> Analysis: ATLAS SUSY 2018 031 <>" << endmsg;
16 INFO << " <> Higgs, sbottom, Bjets + MET @ 13 TeV, 139 fb^-1 <>" << endmsg;
17 INFO << " <> arXiv:1908.03122 <>" << endmsg;
18 INFO << " <> Recasted by : J.Y. Araz <>" << endmsg;
19 INFO << " <> Contact : jack.araz@concordia.ca <>" << endmsg;
20 INFO << " <> Based on MadAnalysis 5 v1.8 and above <>" << endmsg;
21 INFO << " <> For more information, see <>" << endmsg;
22 INFO << " <> http://madanalysis.irmp.ucl.ac.be/wiki/PublicAnalysisDatabase <>" << endmsg;
23 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
24
25 // ========================= //
26 // ===== Signal region ===== //
27 // ========================= //
28
29 Manager()->AddRegionSelection("SRA");
30 Manager()->AddRegionSelection("SRA_L");
31 Manager()->AddRegionSelection("SRA_M");
32 Manager()->AddRegionSelection("SRA_H");
33
34 Manager()->AddRegionSelection("SRB");
35
36 Manager()->AddRegionSelection("SRC");
37 Manager()->AddRegionSelection("SRC_22");
38 Manager()->AddRegionSelection("SRC_24");
39 Manager()->AddRegionSelection("SRC_26");
40 Manager()->AddRegionSelection("SRC_28");
41
42 std::string SRA[] = {"SRA","SRA_L","SRA_M","SRA_H"};
43 std::string SRC[] = {"SRC","SRC_22","SRC_24","SRC_26","SRC_28"};
44 std::string SRABC[] = {"SRA","SRA_L","SRA_M","SRA_H","SRB",
45 "SRC","SRC_22","SRC_24","SRC_26","SRC_28"};
46 std::string SRAB[] = {"SRA","SRA_L","SRA_M","SRA_H","SRB"};
47 std::string SRAB1[] = {"SRA","SRB"};
48
49 // ====================== //
50 // ===== Selections ===== //
51 // ====================== //
52
53 // Common selections
54 Manager()->AddCut("$\\slashed{E}_T$ trigger", SRABC);
55 Manager()->AddCut("$N_{lep} = 0$", SRABC);
56
57 Manager()->AddCut("$N_{j} \\geq 6$", SRA);
58 Manager()->AddCut("$N_{j} \\geq 5$", "SRB");
59 Manager()->AddCut("$N_{j} \\geq 4$", SRC);
60
61 Manager()->AddCut("$N_{b} \\geq 4$", SRAB);
62 Manager()->AddCut("$N_{b} \\geq 3$", SRC);
63
64 Manager()->AddCut("$\\slashed{E}_T > 350$", SRAB);
65 Manager()->AddCut("$\\slashed{E}_T > 250$", SRC);
66
67 Manager()->AddCut("$min(\\Delta\\phi(j_{1-4},\\slashed{E}_T))>0.4$ [rad]",SRABC);
68
69 Manager()->AddCut("$\\tau^h$ veto", SRAB);
70
71 // SRA selections
72 Manager()->AddCut("$p^{b_1}_T > 200$ [GeV]", SRA);
73 Manager()->AddCut("$\\Delta R_{max}(b,b)>2.5$", SRA);
74 Manager()->AddCut("$\\Delta R_{max-min}(b,b)<2.5$", SRA);
75 Manager()->AddCut("$m(h_{cand})>80$ [GeV]", SRA);
76
77 // SRB selections
78 Manager()->AddCut("$m(h_{cand1},h_{cand2})_{avg}\\in[75,175]$ [GeV]", "SRB");
79 Manager()->AddCut("Leading jet non-b-tagged", "SRB");
80 Manager()->AddCut("$p^{j_1}_T > 350$ [GeV]", "SRB");
81 Manager()->AddCut("$|\\Delta\\phi(j_{1},\\slashed{E}_T)|>2.8$ [rad]", "SRB");
82
83 Manager()->AddCut("$m_{eff} > 1$ [TeV]", SRAB1);
84
85 Manager()->AddCut("$m_{eff} \\in [1,1.5]$ [TeV]", "SRA_L");
86 Manager()->AddCut("$m_{eff} \\in [1.5,2]$ [TeV]", "SRA_M");
87 Manager()->AddCut("$m_{eff} > 2$ [TeV]", "SRA_H");
88
89 // SRC selections
90 Manager()->AddCut("$\\slashed{E}_T$ Sig. $>22$", "SRC");
91 Manager()->AddCut("$\\slashed{E}_T$ Sig. $\\in [22,24]$", "SRC_22");
92 Manager()->AddCut("$\\slashed{E}_T$ Sig. $\\in [24,26]$", "SRC_24");
93 Manager()->AddCut("$\\slashed{E}_T$ Sig. $\\in [26,28]$", "SRC_26");
94 Manager()->AddCut("$\\slashed{E}_T$ Sig. $>28$", "SRC_28");
95
96 // ====================== //
97 // ===== Histograms ===== //
98 // ====================== //
99
100 Manager()->AddHisto("SRA_Meff", 11,800.0,3000., "SRA");
101 Manager()->AddHisto("SRA_Mh", 12,0.0,480., "SRA");
102
103 Manager()->AddHisto("SRB_PTj1", 9,50.0,950., "SRB");
104 Manager()->AddHisto("SRB_MhAvg", 16,50.0,450., "SRB");
105
106 Manager()->AddHisto("SRC_MET", 13,200.0,1500., "SRC");
107 Manager()->AddHisto("SRC_Sig", 19,17.0,36., "SRC");
108
109 return true;
110}
111
112// -----------------------------------------------------------------------------
113// Finalize
114// function called one time at the end of the analysis
115// -----------------------------------------------------------------------------
116void atlas_susy_2018_031::Finalize(const SampleFormat& summary,
117 const std::vector<SampleFormat>& files){}
118
119// -----------------------------------------------------------------------------
120// Execute
121// function called each time one event is read
122// -----------------------------------------------------------------------------
123bool atlas_susy_2018_031::Execute(SampleFormat& sample, const EventFormat& event)
124{
125 // Event weight
126 MAdouble64 EvWeight;
127 if(Configuration().IsNoEventWeight()) EvWeight=1.;
128 else if(event.mc()->weight()!=0.) EvWeight=event.mc()->weight();
129 else { return false;}
130 Manager()->InitializeForNewEvent(EvWeight);
131 // Empty event
132 if (event.rec()==0) {return true;}
133 event_num++;
134 std::vector<const RecJetFormat*> BaseJets, SignalBJets, SignalnonBJets, SignalJets;
135 std::vector<const RecLeptonFormat*> BaseMuons, BaseElectrons, SignalMuons, SignalElectrons;
136 std::vector<const RecTauFormat*> SignalTaus;
137 DEBUG << "============== Event " << event_num << endmsg;
138 // Jets
139 for(MAuint32 ij=0; ij<event.rec()->jets().size(); ij++)
140 {
141 const RecJetFormat *CurrentJet = &(event.rec()->jets()[ij]);
142 if ( CurrentJet->pt() > 20.0 && CurrentJet->abseta() < 4.8)
143 {
144 BaseJets.push_back(CurrentJet);
145 }
146 }
147
148 // Taus
149 // Candidates (tau) are identified as jets which have |eta| < 2.5 and
150 // less than five inner detector tracks of pT > 500 MeV.
151 for(MAuint32 ij=0; ij<event.rec()->taus().size(); ij++)
152 {
153 const RecTauFormat *CurrentTau = &(event.rec()->taus()[ij]);
154 if ( CurrentTau->pt() > 0.5 && CurrentTau->abseta()<2.5)
155 {
156 SignalTaus.push_back(CurrentTau);
157 }
158 }
159
160 // Electron with pT > 4.5 GeV & eta < 2.47
161 for(MAuint32 ie=0; ie<event.rec()->electrons().size(); ie++)
162 {
163 const RecLeptonFormat *CurrentElectron = &(event.rec()->electrons()[ie]);
164 if( CurrentElectron->pt()>4.5 && CurrentElectron->abseta() < 2.47)
165 BaseElectrons.push_back(CurrentElectron);
166 }
167
168 // Muons with pT > 4 GeV & eta < 2.5
169 for(MAuint32 im=0; im<event.rec()->muons().size(); im++)
170 {
171 const RecLeptonFormat *CurrentMuon = &(event.rec()->muons()[im]);
172 if( CurrentMuon->pt()>4. && CurrentMuon->abseta() < 2.5)
173 BaseMuons.push_back(CurrentMuon);
174 }
175
176 // Lepton-jet overlap removal
177 BaseJets = Removal(BaseJets, BaseElectrons, 0.2);
178 BaseElectrons = Removal(BaseElectrons, BaseJets, 0.4);
179 SignalMuons = Removal(BaseMuons, BaseJets, 0.4);
180
181 // Select signal jets and bjets from basejets
182 MAdouble64 HT = 0.0;
183 for (MAuint32 i=0; i<BaseJets.size(); i++)
184 {
185 if ( BaseJets[i]->pt() > 30.0 && BaseJets[i]->abseta() < 2.8)
186 {
187 SignalJets.push_back(BaseJets[i]);
188 HT += BaseJets[i]->pt();
189 if(!BaseJets[i]->btag())
190 {
191 SignalnonBJets.push_back(BaseJets[i]);
192 } else if ( BaseJets[i]->abseta() < 2.5 && BaseJets[i]->btag())
193 {
194 SignalBJets.push_back(BaseJets[i]);
195 }
196 }
197 }
198
199 // Select signal electrons and muons from base electrons and muons
200 MAdouble64 LeadLepPt = 0.0;
201 for (MAuint32 i=0; i<BaseElectrons.size(); i++)
202 {
203 if ( BaseElectrons[i]->pt() > 20.0 && BaseElectrons[i]->abseta()<2.47)
204 {
205 SignalElectrons.push_back(BaseElectrons[i]);
206 if (BaseElectrons[i]->pt() > LeadLepPt)
207 LeadLepPt = BaseElectrons[i]->pt();
208 }
209 }
210
211 for (MAuint32 i=0; i<BaseMuons.size(); i++)
212 {
213 if ( BaseMuons[i]->pt() > 20.0 && BaseMuons[i]->abseta()<2.5)
214 {
215 SignalMuons.push_back(BaseMuons[i]);
216 if (BaseMuons[i]->pt() > LeadLepPt)
217 LeadLepPt = BaseMuons[i]->pt();
218 }
219 }
220
221//TrkMiss
222 MALorentzVector TrkMiss;
223 for(MAuint32 i=0; i<event.rec()->tracks().size(); i++)
224 {
225 const RecTrackFormat *CurrentTrack = &(event.rec()->tracks()[i]);
226 TrkMiss -= CurrentTrack->momentum();
227 }
228
229
230 DEBUG << " * Event reconstructed properly..." << endmsg;
231 // The MET & Meff
232 MALorentzVector pTmiss = event.rec()->MET().momentum();
233 MAdouble64 MET = pTmiss.Pt();
234 MAdouble64 Meff = HT + MET;
235
236 // Ordering everything
237 SORTER->sort(SignalJets, PTordering);
238 SORTER->sort(SignalBJets, PTordering);
239 SORTER->sort(SignalnonBJets, PTordering);
240 SORTER->sort(SignalMuons, PTordering);
241 SORTER->sort(SignalElectrons, PTordering);
242
243 MAint32 nBaselineLeptons = BaseMuons.size() + BaseElectrons.size();
244 MAuint32 nJets = SignalJets.size();
245 MAint32 nbJets = SignalBJets.size();
246// MAint32 nLeptons = SignalMuons.size() + SignalElectrons.size();
247 MAbool nonBLead = false;
248
249 if(SignalnonBJets.size() > 0 && nbJets > 0)
250 {
251 if (SignalnonBJets[0]->pt() > SignalBJets[0]->pt())
252 {
253 nonBLead = true;
254 }
255 else nonBLead = false;
256 }
257
258 //minDPhis
259 MAdouble64 minDPhiJMET_4 = 99.;
260 MAdouble64 DPhiJMET_1 = 99.;
261 for (MAuint32 i=0; i<4; i++)
262 {
263 if (nJets >= i+1)
264 {
265 MAdouble64 DPhi = SignalJets[i]->dphi_0_pi(pTmiss);
266 if (DPhi < minDPhiJMET_4) minDPhiJMET_4 = DPhi;
267 }
268 }
269 if (nJets > 0) DPhiJMET_1 = SignalJets[0]->dphi_0_pi(pTmiss);
270
271 MAdouble64 minDPhiTauMet = 99.;
272 for (MAuint32 i=0; i<SignalTaus.size(); i++)
273 {
274 MAdouble64 DPhi = SignalTaus[i]->dphi_0_pi(pTmiss);
275 if (DPhi < minDPhiTauMet) minDPhiTauMet = DPhi;
276 }
277 DEBUG << " * Starting cut-flow..." << endmsg;
278 // Need to be figured out:
279 if (!Manager()->ApplyCut(TrkMiss.Pt()>0.5 && minDPhiTauMet>1.0471976, "$\\slashed{E}_T$ trigger")) return true;
280
281 if (!Manager()->ApplyCut(nBaselineLeptons == 0, "$N_{lep} = 0$")) return true;
282
283 if (!Manager()->ApplyCut(nJets >= 6, "$N_{j} \\geq 6$")) return true;
284 if (!Manager()->ApplyCut(nJets >= 5, "$N_{j} \\geq 5$")) return true;
285 if (!Manager()->ApplyCut(nJets >= 4, "$N_{j} \\geq 4$")) return true;
286
287 if (!Manager()->ApplyCut(nbJets >= 4, "$N_{b} \\geq 4$")) return true;
288 if (!Manager()->ApplyCut(nbJets >= 3, "$N_{b} \\geq 3$")) return true;
289
290 Manager()->FillHisto("SRC_MET", MET);
291 if (!Manager()->ApplyCut(MET > 350.0, "$\\slashed{E}_T > 350$")) return true;
292 if (!Manager()->ApplyCut(MET > 250.0, "$\\slashed{E}_T > 250$")) return true;
293
294 if (!Manager()->ApplyCut(minDPhiJMET_4 > 0.4, "$min(\\Delta\\phi(j_{1-4},\\slashed{E}_T))>0.4$ [rad]")) return true;
295
296 MAdouble64 METtoHT = 0.;
297 if (MET > 0. && HT > 0.) METtoHT = MET/sqrt(HT);
298
299 Manager()->FillHisto("SRC_Sig", METtoHT);
300 if (!Manager()->ApplyCut(METtoHT>22.0, "$\\slashed{E}_T$ Sig. $>22$")) return true;
301 if (!Manager()->ApplyCut(METtoHT>=22.0 && METtoHT<=24.0, "$\\slashed{E}_T$ Sig. $\\in [22,24]$")) return true;
302 if (!Manager()->ApplyCut(METtoHT>=24.0 && METtoHT<=26.0, "$\\slashed{E}_T$ Sig. $\\in [24,26]$")) return true;
303 if (!Manager()->ApplyCut(METtoHT>=26.0 && METtoHT<=28.0, "$\\slashed{E}_T$ Sig. $\\in [26,28]$")) return true;
304 if (!Manager()->ApplyCut(METtoHT > 28.0, "$\\slashed{E}_T$ Sig. $>28$")) return true;
305
306
307 // Tau veto is wrong!!!
308 if (!Manager()->ApplyCut(SignalTaus.size()==0, "$\\tau^h$ veto")) return true;
309
310 if (!Manager()->ApplyCut(SignalBJets[0]->pt() > 200.0, "$p^{b_1}_T > 200$ [GeV]")) return true;
311
312 // =============================== //
313 // ===== Max - Min Algorithm ===== //
314 // =============================== //
315
316 MAdouble64 maxDR = 0.;
317 MAdouble64 maxmin_M = -99.;
318 MAdouble64 max_M = -99.;
319 MAdouble64 max_avg_M = -99.;
320 MAuint32 imax = -1;
321 MAuint32 jmax = -1;
322 MAdouble64 maxminDR = 99.;
323 MAdouble64 maxmaxDR = -1.;
324 MAuint32 imin = -1;
325 MAuint32 jmin = -1;
326 MAuint32 imaxmax = -1;
327 MAuint32 jmaxmax = -1;
328 DEBUG << " * Starting Min-Max Algorithm..." << endmsg;
329 std::vector< std::vector<double> > DR_matrix;
330 for(MAuint32 i=0; i<SignalBJets.size(); i++)
331 {
332 std::vector<MAdouble64> temp;
333 for(MAuint32 j=0; j<=i; j++)
334 {
335 MAdouble64 DR = SignalBJets[i]->dr(SignalBJets[j]);
336 if (i !=j)
337 {
338 if (DR>maxDR)
339 {
340 maxDR = DR;
341 imax = i;
342 jmax = j;
343 }
344 }
345 temp.push_back(DR);
346 }
347 DR_matrix.push_back(temp);
348 }
349
350 max_M = (SignalBJets[imax]->momentum()+SignalBJets[jmax]->momentum()).M();
351
352 for(MAuint32 i=0; i<SignalBJets.size(); i++)
353 {
354 for(MAuint32 j=0; j<=i; j++)
355 {
356 if (i != j)
357 {
358 MAdouble64 DR = DR_matrix[i][j];
359 if (i != imax && j != jmax && DR > maxmaxDR)
360 {
361 imaxmax = i;
362 jmaxmax = j;
363 }
364 if (i != imax && j != jmax && DR < maxminDR)
365 {
366 maxminDR = DR;
367 imin = i;
368 jmin = j;
369 }
370 }
371 }
372 }
373
374 MAdouble64 maxmax_M = 0.;
375
376 if (SignalBJets.size() >= imaxmax && SignalBJets.size() >= jmaxmax && imaxmax >= 0 && jmaxmax >= 0)
377 maxmax_M = (SignalBJets[imaxmax]->momentum() + SignalBJets[jmaxmax]->momentum()).M();
378 if (max_M > 0. && maxmax_M > 0.)
379 max_avg_M = (max_M+maxmax_M)/2.0;
380 if (SignalBJets.size() >= imin && SignalBJets.size() >= jmin && imin >= 0 && jmin >= 0)
381 maxmin_M = (SignalBJets[imin]->momentum() + SignalBJets[jmin]->momentum()).M();
382 if (!Manager()->ApplyCut(maxDR > 2.5, "$\\Delta R_{max}(b,b)>2.5$")) return true;
383 if (!Manager()->ApplyCut(maxminDR < 2.5, "$\\Delta R_{max-min}(b,b)<2.5$")) return true;
384
385 Manager()->FillHisto("SRA_Mh", maxmin_M);
386 if (!Manager()->ApplyCut(maxmin_M > 80.0, "$m(h_{cand})>80$ [GeV]")) return true;
387
388 Manager()->FillHisto("SRB_MhAvg", max_avg_M);
389 if (!Manager()->ApplyCut(max_avg_M > 175.0 && max_avg_M < 75.0,
390 "$m(h_{cand1},h_{cand2})_{avg}\\in[75,175]$ [GeV]")) return true;
391 if (!Manager()->ApplyCut(nonBLead, "Leading jet non-b-tagged")) return true;
392
393 Manager()->FillHisto("SRB_PTj1", SignalJets[0]->pt());
394 if (!Manager()->ApplyCut(SignalJets[0]->pt()>350.0, "$p^{j_1}_T > 350$ [GeV]")) return true;
395 if (!Manager()->ApplyCut(DPhiJMET_1>2.8, "$|\\Delta\\phi(j_{1},\\slashed{E}_T)|>2.8$ [rad]")) return true;
396
397 Manager()->FillHisto("SRA_Meff", Meff);
398 if (!Manager()->ApplyCut(Meff>1000.0, "$m_{eff} > 1$ [TeV]")) return true;
399 if (!Manager()->ApplyCut(Meff>=1000.0 && Meff<=1500.0, "$m_{eff} \\in [1,1.5]$ [TeV]")) return true;
400 if (!Manager()->ApplyCut(Meff>=1500.0 && Meff<=2000.0, "$m_{eff} \\in [1.5,2]$ [TeV]")) return true;
401 if (!Manager()->ApplyCut(Meff>2000.0, "$m_{eff} > 2$ [TeV]")) return true;
402 DEBUG << " * All cuts applied properly..." << endmsg;
403 return true;
404}
405
406// Overlap Removal
407template<typename T1, typename T2> std::vector<const T1*>
408 Removal(std::vector<const T1*> &v1, std::vector<const T2*> &v2,
409 const MAdouble64 &drmin)
410{
411 // Determining with objects should be removed
412 std::vector<bool> mask(v1.size(),false);
413 for (MAuint32 j=0;j<v1.size();j++)
414 for (MAuint32 i=0;i<v2.size();i++)
415 if (v2[i]->dr(v1[j]) < drmin)
416 {
417 mask[j]=true;
418 break;
419 }
420
421 // Building the cleaned container
422 std::vector<const T1*> cleaned_v1;
423 for (MAuint32 i=0;i<v1.size();i++)
424 if (!mask[i]) cleaned_v1.push_back(v1[i]);
425
426 return cleaned_v1;
427}