MA5SandBox: cms_exo_17_030.cpp

File cms_exo_17_030.cpp, 23.7 KB (added by Benjamin Fuks, 5 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/cms_exo_17_030.h"
2#define cms_exo_17_030_PID 1000021
3using namespace MA5;
4using namespace std;
5
6// -----------------------------------------------------------------------------
7// Initialize
8// function called one time at the beginning of the analysis
9// -----------------------------------------------------------------------------
10bool cms_exo_17_030::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
11{
12 // initialize variables, histos
13 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
14 INFO << " <> Analysis: cms_exo_17_030 <>" << endmsg;
15 INFO << " <> arXiv:1810.10092 <>" << endmsg;
16 INFO << " <> (multijet) <>" << endmsg;
17 INFO << " <> Recasters: KANG Yechan,Kim Jihun, <>" << endmsg;
18 INFO << " <> Choi Jin, Yun SooHyun <>" << endmsg;
19 INFO << " <> Contact: choij@cern.ch <>" << endmsg;
20 INFO << " <> soohyun.yun@cern.ch <>" << endmsg;
21 INFO << " <> jihun.k@cern.ch <>" << endmsg;
22 INFO << " <> Based on MadAnalysis 5 v1.8 <>" << endmsg;
23 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
24
25 // Initial setting for the mother particle
26 WARNING << "Triplets will be matched to the particle with PID " << cms_exo_17_030_PID << endmsg;
27 WARNING << "Change cms_exo_17_030_PID for your triplets" << endmsg;
28
29 // Declaration of the signal regions (SR)
30 Manager()->AddRegionSelection("Mg_200to400");
31 Manager()->AddRegionSelection("Mg_400to700");
32 Manager()->AddRegionSelection("Mg_700to1200");
33 Manager()->AddRegionSelection("Mg_1200to2000");
34
35 // initialize the bins
36 for (unsigned int i = 0; i < nBins_SR1; i++) {
37 unsigned int this_bin_left = (firstBin_SR1 + i*binsize_SR1);
38 unsigned int this_bin_right = this_bin_left + binsize_SR1;
39 string this_bin = "SR1: " + to_string(this_bin_left) + "to" + to_string(this_bin_right);
40 string this_cut = "SR1: " + to_string(this_bin_left) + " < M(jjj) < " + to_string(this_bin_right);
41 bins_SR1.push_back(this_bin);
42 cuts_SR1.push_back(this_cut);
43 Manager()->AddRegionSelection(bins_SR1.at(i));
44
45 }
46
47 for (unsigned int i = 0; i < nBins_SR2; i++) {
48 unsigned int this_bin_left = (firstBin_SR2 + i*binsize_SR2);
49 unsigned int this_bin_right = this_bin_left + binsize_SR2;
50 string this_bin = "SR2: " + to_string(this_bin_left) + "to" + to_string(this_bin_right);
51 string this_cut = "SR2: " + to_string(this_bin_left) + " < M(jjj) < " + to_string(this_bin_right);
52 bins_SR2.push_back(this_bin);
53 cuts_SR2.push_back(this_cut);
54 Manager()->AddRegionSelection(bins_SR2.at(i));
55 }
56
57 for (unsigned int i = 0; i < nBins_SR3; i++) {
58 unsigned int this_bin_left = (firstBin_SR3 + i*binsize_SR3);
59 unsigned int this_bin_right = this_bin_left + binsize_SR3;
60 string this_bin = "SR3: " + to_string(this_bin_left) + "to" + to_string(this_bin_right);
61 string this_cut = "SR3: " + to_string(this_bin_left) + " < M(jjj) < " + to_string(this_bin_right);
62 bins_SR3.push_back(this_bin);
63 cuts_SR3.push_back(this_cut);
64 Manager()->AddRegionSelection(bins_SR3.at(i));
65 }
66
67 for (unsigned int i = 0; i < nBins_SR4; i++) {
68 unsigned int this_bin_left = (firstBin_SR4 + i*binsize_SR4);
69 unsigned int this_bin_right = this_bin_left + binsize_SR4;
70 string this_bin = "SR4: " + to_string(this_bin_left) + "to" + to_string(this_bin_right);
71 string this_cut = "SR4: " + to_string(this_bin_left) + " < M(jjj) < " + to_string(this_bin_right);
72 bins_SR4.push_back(this_bin);
73 cuts_SR4.push_back(this_cut);
74 Manager()->AddRegionSelection(bins_SR4.at(i));
75 }
76
77 // Signal Region Partitions
78 // All Regions
79 std::string All_Region[103];
80 All_Region[0] = "Mg_200to400";
81 All_Region[1] = "Mg_400to700";
82 All_Region[2] = "Mg_700to1200";
83 All_Region[3] = "Mg_1200to2000";
84
85 for (unsigned int i = 0; i < nBins_SR1; i++) {
86 All_Region[4+i] = bins_SR1.at(i);
87 }
88 for (unsigned int i = 0; i < nBins_SR2; i++) {
89 All_Region[4+nBins_SR1+i] = bins_SR2.at(i);
90 }
91 for (unsigned int i = 0; i < nBins_SR3; i++) {
92 All_Region[4+nBins_SR1+nBins_SR2+i] = bins_SR3.at(i);
93 }
94 for (unsigned int i = 0; i < nBins_SR4; i++) {
95 All_Region[4+nBins_SR1+nBins_SR2+nBins_SR3+i] = bins_SR4.at(i);
96 }
97
98 // Low Mass Regions
99 std::string Low_Mass_Region[64];
100 Low_Mass_Region[0] = "Mg_200to400";
101 Low_Mass_Region[1] = "Mg_400to700";
102 for (unsigned int i = 0; i < nBins_SR1; i++)
103 Low_Mass_Region[2+i] = bins_SR1.at(i);
104 for (unsigned int i = 0; i < nBins_SR2; i++)
105 Low_Mass_Region[2+nBins_SR1+i] = bins_SR2.at(i);
106
107 // High Mass Regions
108 std::string High_Mass_Region[39];
109 High_Mass_Region[0] = "Mg_700to1200";
110 High_Mass_Region[1] = "Mg_1200to2000";
111 for (unsigned int i = 0; i < nBins_SR3; i++)
112 High_Mass_Region[2+i] = bins_SR3.at(i);
113 for (unsigned int i = 0; i < nBins_SR4; i++)
114 High_Mass_Region[2+nBins_SR3+i] = bins_SR4.at(i);
115
116 // SR1
117 std::string SR1[31];
118 SR1[0] = "Mg_200to400";
119 for (unsigned int i = 0; i < nBins_SR1; i++)
120 SR1[1+i] = bins_SR1.at(i);
121
122 // SR2
123 std::string SR2[33];
124 SR2[0] = "Mg_400to700";
125 for (unsigned int i = 0; i < nBins_SR2; i++)
126 SR2[1+i] = bins_SR2.at(i);
127
128 // SR3
129 std::string SR3[18];
130 SR3[0] = "Mg_700to1200";
131 for (unsigned int i = 0; i < nBins_SR3; i++)
132 SR3[1+i] = bins_SR3.at(i);
133
134 // SR4
135 std::string SR4[21];
136 SR4[0] = "Mg_1200to2000";
137 for (unsigned int i = 0; i < nBins_SR4; i++)
138 SR4[1+i] = bins_SR4.at(i);
139
140 // Declaration of the preselection and HT cuts
141 Manager()->AddCut("ALL: Njets>=6", All_Region);
142
143 Manager()->AddCut("LOW: preselection", Low_Mass_Region);
144 Manager()->AddCut("HIGH: preselection", High_Mass_Region);
145
146 Manager()->AddCut("LOW: HT > 650GeV", Low_Mass_Region);
147 Manager()->AddCut("HIGH: HT > 900GeV", High_Mass_Region);
148
149 // Declaration of the sixth jet-pt cut
150 Manager()->AddCut("SR1: pt(j6) > 40GeV", SR1);
151 Manager()->AddCut("SR2: pt(j6) > 50GeV", SR2);
152 Manager()->AddCut("SR3: pt(j6) > 125GeV", SR3);
153 Manager()->AddCut("SR4: pt(j6) > 175GeV", SR4);
154
155 // Declaration of the D^2[6,3+3,2] cut
156 Manager()->AddCut("SR1: D^2[6,3+3,2] < 1.25", SR1);
157 Manager()->AddCut("SR2: D^2[6,3+3,2] < 1.0", SR2);
158 Manager()->AddCut("SR3: D^2[6,3+3,2] < 0.9", SR3);
159 Manager()->AddCut("SR4: D^2[6,3+3,2] < 0.75", SR4);
160
161 // Declaration of the Am cut - TripletPair cut
162 Manager()->AddCut("SR1: Am < 0.25", SR1);
163 Manager()->AddCut("SR2: Am < 0.175", SR2);
164 Manager()->AddCut("SR3: Am < 0.15", SR3);
165 Manager()->AddCut("SR4: Am < 0.15", SR4);
166
167 // Declaration of the Delta cut - Triplet cut
168 Manager()->AddCut("SR1: Delta > 250GeV", SR1);
169 Manager()->AddCut("SR2: Delta > 180GeV", SR2);
170 Manager()->AddCut("SR3: Delta > 20GeV", SR3);
171 Manager()->AddCut("SR4: Delta > -120GeV", SR4);
172
173 // Declaration of the D^2[3,2] cut - Triplet cut
174 Manager()->AddCut("SR1: D^2[3,2] < 0.05", SR1);
175 Manager()->AddCut("SR2: D^2[3,2] < 0.175", SR2);
176 Manager()->AddCut("SR3: D^2[3,2] < 0.2", SR3);
177 Manager()->AddCut("SR4: D^2[3,2] < 0.25", SR4);
178
179 // mass partition
180 Manager()->AddHisto("SR1: M(jjj)", 200, 200., 400.);
181 Manager()->AddHisto("SR2: M(jjj)", 300, 400., 700.);
182 Manager()->AddHisto("SR3: M(jjj)", 500, 700., 1200.);
183 Manager()->AddHisto("SR4: M(jjj)", 800, 1200., 2000.);
184
185 // Declaration of mass mins
186 for (unsigned int i = 0; i < nBins_SR1; i++)
187 Manager()->AddCut(cuts_SR1[i], bins_SR1[i]);
188
189 for (unsigned int i = 0; i < nBins_SR2; i++)
190 Manager()->AddCut(cuts_SR2[i], bins_SR2[i]);
191
192 for (unsigned int i = 0; i < nBins_SR3; i++)
193 Manager()->AddCut(cuts_SR3[i], bins_SR3[i]);
194
195 for (unsigned int i = 0; i < nBins_SR4; i++)
196 Manager()->AddCut(cuts_SR4[i], bins_SR4[i]);
197
198 return true;
199}
200
201// -----------------------------------------------------------------------------
202// Finalize
203// function called one time at the end of the analysis
204// -----------------------------------------------------------------------------
205void cms_exo_17_030::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files) {}
206
207// -----------------------------------------------------------------------------
208// Execute
209// function called each time one event is read
210// -----------------------------------------------------------------------------
211bool cms_exo_17_030::Execute(SampleFormat& sample, const EventFormat& event)
212{
213 // Event weight
214 double weight;
215 if (Configuration().IsNoEventWeight()) weight = 1.;
216 else if (event.mc()->weight() != 0.) weight = event.mc()->weight();
217 else return false;
218
219 const double nTrips = 20.;
220 const int nRegions = 4;
221 Manager()->InitializeForNewEvent(weight*nTrips);
222
223 const double ptCut[] = {30, 30, 50, 50};
224 const double HTcut[] = {650, 650, 900, 900};
225 const double j6ptCut[] = {40, 50, 125, 175};
226 const double mds6332Cut[] = {1.25, 1., 0.9, 0.75};
227 const double asymmCut[] = {0.25, 0.175, 0.15, 0.15};
228 const double deltaCut[] = {250, 180, 20, -120};
229 const double mds32Cut[] = {0.05, 0.175, 0.2, 0.25};
230
231 // The event loop start here
232 if(event.rec() == 0) return true;
233
234 // Select jets satisfying pt&eta cut
235 const double etaCut = 2.4;
236 const double basePtCut = 20.;
237 JetCollection jets = jetSelection(event, basePtCut, etaCut);
238 SORTER->sort(jets);
239
240 // Nj cut for low and high mass regions
241 if(!Manager()->ApplyCut(jets.size() >= 6, "ALL: Njets>=6")) return true;
242
243 // confirm jet pt cut
244 if(!Manager()->ApplyCut(jets.at(5)->pt() > ptCut[0], "LOW: preselection")) return true;
245 if(!Manager()->ApplyCut(jets.at(5)->pt() > ptCut[2], "HIGH: preselection")) return true;
246
247 // HT, j6pt cut
248 double H_T = HT(jets);
249 double j6pt = jets.at(5)->pt();
250
251 // HT cut for low and high mass regions
252 if(!Manager()->ApplyCut(H_T > HTcut[0], "LOW: HT > 650GeV")) return true;
253 if(!Manager()->ApplyCut(H_T > HTcut[2], "HIGH: HT > 900GeV")) return true;
254
255 // 6th Jet Pt cut for each regions
256 if(!Manager()->ApplyCut(j6pt > j6ptCut[0], "SR1: pt(j6) > 40GeV") ) return true;
257 if(!Manager()->ApplyCut(j6pt > j6ptCut[1], "SR2: pt(j6) > 50GeV") ) return true;
258 if(!Manager()->ApplyCut(j6pt > j6ptCut[2], "SR3: pt(j6) > 125GeV")) return true;
259 if(!Manager()->ApplyCut(j6pt > j6ptCut[3], "SR4: pt(j6) > 175GeV")) return true;
260
261 // Mds6332 cut
262 double evtMds6332 = mds6332(jets);
263
264 if(!Manager()->ApplyCut(evtMds6332 < mds6332Cut[0], "SR1: D^2[6,3+3,2] < 1.25")) return true;
265 if(!Manager()->ApplyCut(evtMds6332 < mds6332Cut[1], "SR2: D^2[6,3+3,2] < 1.0") ) return true;
266 if(!Manager()->ApplyCut(evtMds6332 < mds6332Cut[2], "SR3: D^2[6,3+3,2] < 0.9") ) return true;
267 if(!Manager()->ApplyCut(evtMds6332 < mds6332Cut[3], "SR4: D^2[6,3+3,2] < 0.75")) return true;
268
269
270 // Mass asymmetry cut
271 PairCollection tripPairs = makePairCollection(jets);
272 TripletCollection trips[nRegions];
273 for (int i = 0; i < nRegions; i++) {
274 trips[i] = pairSelection(tripPairs, asymmCut[i]);
275 trips[i] = GenMatchedTriplets(event, trips[i]);
276 }
277
278 double nTrips_passAsymm[nRegions];
279 for (int i = 0; i < nRegions; i++)
280 nTrips_passAsymm[i] = trips[i].size();
281
282 // starting to implement triplet size as weight
283 Manager()->SetCurrentEventWeight(weight*nTrips_passAsymm[0]);
284 if(!Manager()->ApplyCut(trips[0].size() != 0, "SR1: Am < 0.25")) return true;
285 Manager()->SetCurrentEventWeight(weight*nTrips_passAsymm[1]);
286 if(!Manager()->ApplyCut(trips[1].size() != 0, "SR2: Am < 0.175")) return true;
287 Manager()->SetCurrentEventWeight(weight*nTrips_passAsymm[2]);
288 if(!Manager()->ApplyCut(trips[2].size() != 0, "SR3: Am < 0.15")) return true;
289 Manager()->SetCurrentEventWeight(weight*nTrips_passAsymm[3]);
290 Manager()->SetCurrentEventWeight(weight*trips[3].size());
291 if(!Manager()->ApplyCut(trips[3].size() != 0, "SR4: Am < 0.15")) return true;
292
293 // Delta cut
294 double nTrips_passDelta[nRegions];
295 for (int i = 0; i < nRegions; i++) {
296 // select triplets
297 trips[i] = deltaSelection(trips[i], deltaCut[i]);
298
299 // update nTrips
300 nTrips_passDelta[i] = trips[i].size();
301 }
302
303 Manager()->SetCurrentEventWeight(weight*nTrips_passDelta[0]);
304 if(!Manager()->ApplyCut(trips[0].size() != 0, "SR1: Delta > 250GeV")) return true;
305 Manager()->SetCurrentEventWeight(weight*nTrips_passDelta[1]);
306 if(!Manager()->ApplyCut(trips[1].size() != 0, "SR2: Delta > 180GeV")) return true;
307 Manager()->SetCurrentEventWeight(weight*nTrips_passDelta[2]);
308 if(!Manager()->ApplyCut(trips[2].size() != 0, "SR3: Delta > 20GeV")) return true;
309 Manager()->SetCurrentEventWeight(weight*nTrips_passDelta[3]);
310 if(!Manager()->ApplyCut(trips[3].size() != 0, "SR4: Delta > -120GeV")) return true;
311
312 // mds32 cut
313 double nTrips_passMDS32[nRegions];
314 for (int i = 0; i < nRegions; i++) {
315 // select triplets
316 trips[i] = mds32Selection(trips[i], mds32Cut[i]);
317
318 // update nTrips
319 nTrips_passMDS32[i] = trips[i].size();
320 }
321
322 Manager()->SetCurrentEventWeight(weight*nTrips_passMDS32[0]);
323 if(!Manager()->ApplyCut(trips[0].size() != 0, "SR1: D^2[3,2] < 0.05")) return true;
324 Manager()->SetCurrentEventWeight(weight*nTrips_passMDS32[1]);
325 if(!Manager()->ApplyCut(trips[1].size() != 0, "SR2: D^2[3,2] < 0.175")) return true;
326 Manager()->SetCurrentEventWeight(weight*nTrips_passMDS32[2]);
327 if(!Manager()->ApplyCut(trips[2].size() != 0, "SR3: D^2[3,2] < 0.2")) return true;
328 Manager()->SetCurrentEventWeight(weight*nTrips_passMDS32[3]);
329 if(!Manager()->ApplyCut(trips[3].size() != 0, "SR4: D^2[3,2] < 0.25")) return true;
330
331 // Filling Histo
332 Manager()->SetCurrentEventWeight(weight);
333 unsigned int nTrips_SR1[nBins_SR1] = {0}, nTrips_SR2[nBins_SR2] = {0}, nTrips_SR3[nBins_SR3] = {0}, nTrips_SR4[nBins_SR4] = {0};
334
335 for (unsigned int i = 0; i < trips[0].size(); i++) {
336 Manager()->FillHisto("SR1: M(jjj)", mass(trips[0].at(i)));
337
338 // prepare for region selector
339 double this_mass = mass(trips[0].at(i));
340 unsigned int this_bin = (((int)this_mass) - firstBin_SR1) / binsize_SR1;
341 if (! (this_bin&&this_bin<nBins_SR1))
342 continue;
343 nTrips_SR1[this_bin]++;
344 }
345
346 for (unsigned int i = 0; i < trips[1].size(); i++) {
347 Manager()->FillHisto("SR2: M(jjj)", mass(trips[1].at(i)));
348
349 // prepare for region selector
350 double this_mass = mass(trips[1].at(i));
351 unsigned int this_bin = (((int)this_mass)-firstBin_SR2) / binsize_SR2;
352 if (! (this_bin&&this_bin<nBins_SR2))
353 continue;
354 nTrips_SR2[this_bin]++;
355 }
356
357 for (unsigned int i = 0; i < trips[2].size(); i++) {
358 Manager()->FillHisto("SR3: M(jjj)", mass(trips[2].at(i)));
359
360 // prepare for region selector
361 double this_mass = mass(trips[2].at(i));
362 unsigned int this_bin = (((int)this_mass)-firstBin_SR3) / binsize_SR3;
363 if (! (this_bin&&this_bin<nBins_SR3))
364 continue;
365 nTrips_SR3[this_bin]++;
366 }
367
368 for (unsigned int i = 0; i < trips[3].size(); i++) {
369 Manager()->FillHisto("SR4: M(jjj)", mass(trips[3].at(i)));
370
371 // prepare for region selector
372 double this_mass = mass(trips[3].at(i));
373 unsigned int this_bin = (((int)this_mass)-firstBin_SR4) / binsize_SR4;
374 if (! (this_bin&&this_bin<nBins_SR4))
375 continue;
376 nTrips_SR4[this_bin]++;
377 }
378
379 // Final mass bin
380 for (unsigned int i = 0; i < nBins_SR1; i++) {
381 Manager()->SetCurrentEventWeight(weight*(double)nTrips_SR1[i]);
382 if(!Manager()->ApplyCut(nTrips_SR1[i] != 0, cuts_SR1[i])) return true;
383 }
384
385 for (unsigned int i = 0; i < nBins_SR2; i++) {
386 Manager()->SetCurrentEventWeight(weight*(double)nTrips_SR2[i]);
387 if(!Manager()->ApplyCut(nTrips_SR2[i] != 0, cuts_SR2[i])) return true;
388 }
389
390 for (unsigned int i = 0; i < nBins_SR3; i++) {
391 Manager()->SetCurrentEventWeight(weight*(double)nTrips_SR3[i]);
392 if(!Manager()->ApplyCut(nTrips_SR3[i] != 0, cuts_SR3[i])) return true;
393 }
394
395 for (unsigned int i = 0; i < nBins_SR4; i++) {
396 Manager()->SetCurrentEventWeight(weight*(double)nTrips_SR4[i]);
397 if (!Manager()->ApplyCut(nTrips_SR4[i] != 0, cuts_SR4[i])) return true;
398 }
399
400 return true;
401}
402
403//==== Functions to use
404//==== 1. kinematic variables
405MALorentzVector cms_exo_17_030::momentum(const Triplet &trip) {
406 return trip[0]->momentum() + trip[1]->momentum() + trip[2]->momentum();
407}
408
409double cms_exo_17_030::mass(const Triplet &trip) {
410 MALorentzVector mom = momentum(trip);
411 return mom.M();
412}
413
414double cms_exo_17_030::HT(const JetCollection &jetcoll) {
415 double this_HT = 0;
416 for (unsigned int i = 0; i < jetcoll.size(); i++)
417 this_HT += jetcoll.at(i)->pt();
418
419 return this_HT;
420}
421
422double cms_exo_17_030::dalitz32(const Triplet &trip, const int &idx1, const int &idx2) {
423 double M = pow((trip[idx1]->momentum()+trip[idx2]->momentum()).M(), 2);
424 double M_123 = pow((trip[0]->momentum()+trip[1]->momentum()+trip[2]->momentum()).M(), 2);
425 double M_1 = pow((trip[0]->momentum()).M(), 2);
426 double M_2 = pow((trip[1]->momentum()).M(), 2);
427 double M_3 = pow((trip[2]->momentum()).M(), 2);
428 return M/(M_123+M_1+M_2+M_3);
429}
430
431double cms_exo_17_030::mds32(const Triplet &trip) {
432 double c = 1./sqrt(3.);
433 double out = 0.;
434 for (int i = 0; i < 3; i++) {
435 for (int j = i+1; j < 3; j++) {
436 out += pow( sqrt(dalitz32(trip, i, j))-c, 2);
437 }
438 }
439 return out;
440}
441
442double cms_exo_17_030::dalitz6332(const JetCollection &jetcoll, const int &idx1, const int &idx2, const int &idx3) {
443 double M = pow((jetcoll.at(idx1)->momentum() + jetcoll.at(idx2)->momentum() + jetcoll.at(idx3)->momentum() ).M(), 2);
444 MALorentzVector PTv;
445 double massSum = 0;
446 for (int i = 0; i < 6; i++) {
447 PTv += jetcoll.at(i)->momentum();
448 massSum += pow((jetcoll.at(i)->momentum()).M(), 2);
449 }
450 double M_inv = pow(PTv.M(), 2);
451
452 return M / (4*M_inv + 6*massSum);
453}
454
455double cms_exo_17_030::mds6332(const JetCollection &jetcoll) {
456 double out = 0; double c = 1./sqrt(20.);
457
458 for (int i = 0; i < 6; i++) {
459 for (int j = i+1; j < 6; j++) {
460 for (int k = j+1; k < 6; k++) {
461 Triplet this_trip = {jetcoll.at(i), jetcoll.at(j), jetcoll.at(k)};
462 double this_sum = dalitz6332(jetcoll, i, j, k) + mds32(this_trip);
463 this_sum = sqrt(this_sum) - c;
464 out += pow(this_sum, 2);
465 }
466 }
467 }
468 return out;
469}
470
471double cms_exo_17_030::massAsymm(const TripletPair &pair) {
472 double mass1 = mass(pair.first);
473 double mass2 = mass(pair.second);
474
475 return fabs(mass1 - mass2) / (mass1 + mass2);
476}
477
478double cms_exo_17_030::delta(const Triplet &trip) {
479 double this_HT = 0.;
480 for (unsigned int i = 0; i < trip.size(); i++)
481 this_HT += trip.at(i)->pt();
482
483 double this_mass = mass(trip);
484
485 return this_HT - this_mass;
486}
487
488//==== 2. selections
489JetCollection cms_exo_17_030::jetSelection(const EventFormat &event, const double &ptCut, const double &etaCut) {
490 JetCollection out;
491
492 for (unsigned int i = 0; i < (event.rec()->jets()).size(); i++) {
493 const RecJetFormat& this_jet = (event.rec()->jets()).at(i);
494 if ((this_jet.pt() > ptCut) && (fabs(this_jet.eta()) < etaCut))
495 out.push_back(&this_jet);
496 }
497
498 return out;
499}
500
501TripletCollection cms_exo_17_030::pairSelection(const PairCollection &pairs, const double &asymmCut) {
502 TripletCollection out;
503
504 for (unsigned int i = 0; i < pairs.size(); i++) {
505 const TripletPair& this_pair = pairs.at(i);
506 if (massAsymm(this_pair) < asymmCut) {
507 out.push_back(this_pair.first);
508 out.push_back(this_pair.second);
509 }
510 }
511
512 return out;
513}
514
515TripletCollection cms_exo_17_030::deltaSelection(const TripletCollection &trips, const double &deltaCut) {
516 TripletCollection out;
517
518 for (unsigned int i = 0; i < trips.size(); i++) {
519 const Triplet& this_trip = trips.at(i);
520 if (delta(this_trip) > deltaCut)
521 out.push_back(this_trip);
522 }
523
524 return out;
525}
526
527TripletCollection cms_exo_17_030::mds32Selection(const TripletCollection &trips, const double &mdsCut) {
528 TripletCollection out;
529
530 for (unsigned int i = 0; i < trips.size(); i++) {
531 const Triplet& this_trip = trips.at(i);
532 if (mds32(this_trip) < mdsCut)
533 out.push_back(this_trip);
534 }
535
536 return out;
537}
538
539//==== 3. tools
540PairCollection cms_exo_17_030::makePairCollection(const JetCollection &jetcoll) {
541 PairCollection out;
542
543 TripletCollection trips;
544 for (int i = 0; i < 6; i++) {
545 for (int j = i+1; j < 6; j++) {
546 for (int k = j+1; k < 6; k++) {
547 trips.push_back({jetcoll.at(i), jetcoll.at(j), jetcoll.at(k)});
548 }
549 }
550 }
551
552 for (int i = 0; i < 10; i++) out.push_back(make_pair(trips.at(i), trips.at(19-i)));
553 return out;
554}
555
556TripletCollection cms_exo_17_030::GenMatchedTriplets(const EventFormat &event, const TripletCollection &trips) {
557 JetCollection matched_jets1, matched_jets2, matched_jets3, matched_jets4;
558 TripletCollection matched_trips;
559 const vector<MCParticleFormat> &gens = event.mc()->particles();
560 vector<const MCParticleFormat*> p_gluinos;
561
562 for (unsigned int i = 0; i < gens.size(); i++) {
563 const MCParticleFormat& gen = gens.at(i);
564 if(gen.mothers().size()==0) continue;
565 if (1<=gen.pdgid()&&gen.pdgid()<=4&&(gen.mothers()).at(0)->pdgid()==cms_exo_17_030_PID)
566 p_gluinos.push_back((gen.mothers()).at(0));
567 if (-4<=gen.pdgid()&&gen.pdgid()<=-1&&(gen.mothers()).at(0)->pdgid()==cms_exo_17_030_PID)
568 p_gluinos.push_back((gen.mothers()).at(0));
569 }
570
571 sort(p_gluinos.begin(), p_gluinos.end());
572 vector<const MCParticleFormat*>::iterator last;
573 last = unique(p_gluinos.begin(), p_gluinos.end());
574 p_gluinos.erase(last, p_gluinos.end());
575
576 // start matching
577 // triplets should be matched to the "same" mother gluino -> use the pointer comparision
578 bool matched;
579 for (unsigned int i = 0; i < trips.size(); i++) {
580 for (unsigned int j = 0; j < gens.size(); j++) {
581 const Triplet& trip = trips.at(i);
582 const MCParticleFormat& gen = gens.at(j);
583 if (gen.mothers().size()==0 || p_gluinos.size()==0) continue;
584 else if (1<=gen.pdgid()&&gen.pdgid()<=4&&((gen.mothers()).at(0)==p_gluinos.at(0))) {
585 for (unsigned int k = 0; k < trip.size(); k++) {
586 const RecJetFormat* jet = trip.at(k);
587 matched = (jet->momentum()).DeltaR(gen.momentum()) < 0.3;
588 if (matched) matched_jets1.push_back(jet);
589 }
590 }
591 else if (1<=gen.pdgid()&&gen.pdgid()<=4&&((gen.mothers()).at(0)==p_gluinos.at(1))) {
592 for (unsigned int k = 0; k < trip.size(); k++) {
593 const RecJetFormat* jet = trip.at(k);
594 matched = (jet->momentum()).DeltaR(gen.momentum()) < 0.3;
595 if (matched) matched_jets2.push_back(jet);
596 }
597 }
598 else if (-4<=gen.pdgid()&&gen.pdgid()<=-1&&(gen.mothers()).at(0)==p_gluinos.at(0)) {
599 for (unsigned int k = 0; k < trip.size(); k++) {
600 const RecJetFormat* jet = trip.at(k);
601 matched = (jet->momentum()).DeltaR(gen.momentum()) < 0.3;
602 if (matched) matched_jets3.push_back(jet);
603 }
604 }
605 else if (-4<=gen.pdgid()&&gen.pdgid()<=-1&&(gen.mothers()).at(0)==p_gluinos.at(1)) {
606 for (unsigned int k = 0; k < trip.size(); k++) {
607 const RecJetFormat* jet = trip.at(k);
608 matched = (jet->momentum()).DeltaR(gen.momentum()) < 0.3;
609 if (matched) matched_jets4.push_back(jet);
610 }
611 }
612
613 }
614
615 // remove double matching
616 JetCollection::iterator last1, last2, last3, last4;
617 sort(matched_jets1.begin(), matched_jets1.end());
618 last1 = unique(matched_jets1.begin(), matched_jets1.end());
619 matched_jets1.erase(last1, matched_jets1.end());
620
621 sort(matched_jets2.begin(), matched_jets2.end());
622 last2 = unique(matched_jets2.begin(), matched_jets2.end());
623 matched_jets2.erase(last2, matched_jets2.end());
624
625 sort(matched_jets3.begin(), matched_jets3.end());
626 last3 = unique(matched_jets3.begin(), matched_jets3.end());
627 matched_jets3.erase(last3, matched_jets3.end());
628
629 sort(matched_jets4.begin(), matched_jets4.end());
630 last4 = unique(matched_jets4.begin(), matched_jets4.end());
631 matched_jets4.erase(last4, matched_jets4.end());
632
633 if (matched_jets1.size() == 3) matched_trips.push_back({matched_jets1.at(0), matched_jets1.at(1), matched_jets1.at(2)});
634 if (matched_jets2.size() == 3) matched_trips.push_back({matched_jets2.at(0), matched_jets2.at(1), matched_jets2.at(2)});
635 if (matched_jets3.size() == 3) matched_trips.push_back({matched_jets3.at(0), matched_jets3.at(1), matched_jets3.at(2)});
636 if (matched_jets4.size() == 3) matched_trips.push_back({matched_jets4.at(0), matched_jets4.at(1), matched_jets4.at(2)});
637
638 // prepare for next iteration
639 matched_jets1.clear();
640 matched_jets2.clear();
641 matched_jets3.clear();
642 matched_jets4.clear();
643 }
644 return matched_trips;
645}
646