SFS: cms_exo_19_010.cpp

File cms_exo_19_010.cpp, 45.4 KB (added by Benjamin Fuks, 5 years ago)
Line 
1/*
2
3 CMS EXO-19-010 disappearing track search 139 fb^-1
4 -- arXiv:2004.05153
5 Recast By Mark Goodsell (goodsell@lpthe.jussieu.fr)
6 -- arXiv:2106.08815
7
8*/
9
10
11#include "SampleAnalyzer/User/Analyzer/cms_exo_19_010.h"
12using namespace MA5;
13using namespace std;
14
15
16
17void get_which_hits(const std::vector<double> &radii, std::vector<bool> &outhits, double maxR, double Rdec, double Rprod)
18{
19 //std::vector<bool> outhits;
20 for(auto discR : radii)
21 {
22 if(discR > maxR) return;
23
24 if(Rprod > discR)
25 {
26 outhits.push_back(false);
27 continue;
28 }
29
30 if(Rdec < discR)
31 {
32 outhits.push_back(false);
33 continue;
34 }
35
36 outhits.push_back(true);
37
38 }
39
40
41}
42
43std::vector<bool> get_truth_hits(const RecTrackFormat* p)
44{
45 // This function to determine how many tracker layers are hit
46
47 // strategy: Determine truth level and then apply per-hit efficiencies
48 // Start with pixel detector and work outwards: TIB, TID, TEC, TOB
49 // Pixel detector: https://lss.fnal.gov/archive/design/fermilab-design-2012-02.pdf figure 1.2, 2.1, table 2.1
50 // Tracker: 1405.6569
51
52 std::vector<bool> outhits;
53 outhits.clear();
54
55
56 //MAfloat64 abseta=fabs(p->etaCalo());
57 MAfloat64 abseta=fabs(p->mc()->momentum().Eta());
58
59 // don't bother with large eta
60 if(abseta > 2.5)
61 {
62 //outhits.push_back(false);
63 return outhits;
64 }
65
66 MALorentzVector vstart = p->ProductionVertex();
67
68 double ppz = fabs(p->momentum().Pz());
69 double ppT=p->momentum().Pt();
70
71 double ootantheta=0.0;
72 if(ppz > 0.0)
73 {
74 ootantheta = ppz/ppT; // nb for abseta > 2.5 we don't worry about zeroes here
75 }
76
77 double vdecz, vdecR;
78 if(p->mc()->ctau() > 1.0e-3)
79 {
80 MALorentzVector vdec = p->mc()->decay_vertex();
81 vdecz=fabs(vdec.Pz());
82 vdecR=vdec.Pt();
83 }
84 else
85 {
86 vdecR=1.0e6;
87 vdecz=vdecR*ootantheta;
88 }
89
90 // for now assume tracks start at z=0
91 double vstartz=fabs(vstart.Pz());
92 double vstartR=vstart.Pt();
93
94
95
96
97 // Pixel
98 // https://arxiv.org/pdf/0911.5434.pdf for
99 // Pixel barrel layers at e = 3, 6.8, 10.2, 16 cm
100 // pixel z extends to |z| = 53.3/2 cm
101 // This means |eta| < 1.3 is entirely inside the pixel barrel
102 // Pre-upgrade: The endcap disks, extending from 6 to 15 cm in radius, are placed at z = ±35.5 cm and z = ±48.5 cm.
103 // Post-upgrade: three endcap disks, see figure 2.1 in https://lss.fnal.gov/archive/design/fermilab-design-2012-02.pdf
104 // let us guess that these are at 30, 40, 50 cm
105
106 //1710.03842: Theradii of the four barrel layers are 2.9 cm, 6.8 cm, 10.9 cm, and 16.0 cm.
107 //The three forward disks arepositioned alongzat 3.2 cm, 3.9 cm, and 4.8 cm.
108 // Each disk is composed of two rings of moduleswith average radii of 12.8 cm and 7.8 cm.
109
110 //static vector<double> pixelbarrelR= {30.0,68.0,102.0,160.0};
111 static vector<double> pixelbarrelR= {29.0,68.0,109.0,160.0};
112
113
114 //static vector<double> pixelForwardZ = {300.0,400.0,500.0};
115 static vector<double> pixelForwardZ = {320.0,390.0,480.0};
116
117
118 // cannot escape inner barrel of pixel altogether
119 double maxpixelR;
120 if(ootantheta > 1)
121 {
122 maxpixelR=266.5/ootantheta;
123 }
124 else
125 {
126 maxpixelR=200.0;
127 }
128
129 get_which_hits(pixelbarrelR, outhits, maxpixelR, vdecR, vstartR);
130
131 // check if we need pixel endcaps
132 if(abseta > 1.3)
133 {
134
135 double maxpixelZ=150.0*ootantheta;
136
137 get_which_hits(pixelForwardZ, outhits, maxpixelZ, vdecz, vstartz);
138
139 }
140
141 if(outhits.size() < 4)
142 {
143 //std::cout << "Track shorter than 4 pixel hits! " << abseta << std::endl;
144 // Have I done something wrong for this case? Here I will pad the tracks with "false"
145 // because we first check if there are any missing pixel hits
146 for(int i=0; i < 4-outhits.size(); i++)
147 {
148 outhits.push_back(false);
149 }
150
151 }
152
153
154
155
156 // Now for tracker barrel, see CMS 2008 report for some of these
157 /*
158 The Tracker Inner Barrel (TIB)
159and Disks (TID) cover r < 55 cm and | z | < 118 cm, and are composed of four barrel layers,
160supplemented by three disks at each end. The Tracker Outer Barrel (TOB) covers r > 55 cm and
161| z | < 118 cm and consists of six barrel layers. The Tracker EndCaps (TEC) cover the region 124 <
162| z | < 282 cm. Each TEC is composed of nine disks, each containing up to seven concentric
163rings of silicon strip modules, yielding a range of resolutions similar to that of the TOB.
164*/
165
166
167 //static vector<double> TIBR= {230.0, 300.0, 400.0,500.0};
168 // See CMS 2008 report page 64, these exetend -700 to 700 mm
169 static vector<double> TIBR= {255.0, 339.0, 418.5,498.0};
170
171 // The TID± are assemblies of three disks placed in z between ±800 mm and ±900 mm. The
172 // disks are identical and each one consists of three rings which span the radius from roughly 200 mm
173 // to 500 mm.
174 // Here I use values inferred from Figure 3.34 in the CMS 2008 report
175 static vector<double> TIDZ= {775.0,900.0,1025.0};
176
177 // Close!!
178 //static vector<double> TOBR = {600.0,700.0,800.0,890.0,980.0,1060.0};
179 // See CMS 2008 report page 68
180 static vector<double> TOBR = {608.0,692.0,780.0,868.0,965.0,1080.0};
181
182 // Inferred from Figure 3.34 in the CMS 2008 report
183 static vector<double> TECZ= {1250.0,1400.0,1550.0,1700.0,1950.0,2000.0,2225.0,2450.0,2700.0};
184
185 // inner barrel
186 double maxIBR;
187
188 // inner barrel up to |z| = 580mm, r < 550
189 //if(ootantheta < 1.05)
190 // inner barrel up to |z| = 700mm, r < 550 -> 1/tanbeta < 700/550 = 1.27
191 if(ootantheta < 1.27)
192 {
193 maxIBR=550.0;
194 get_which_hits(TIBR, outhits, maxIBR, vdecR, vstartR);
195 }
196 else
197 {
198 maxIBR=550.0/ootantheta;
199 get_which_hits(TIBR, outhits, maxIBR, vdecR, vstartR);
200 //outer inner endcap, only possible if ootantheta > 1.05
201
202 //double maxIBZ=580.0*ootantheta;
203 double maxIBZ=700.0*ootantheta;
204 get_which_hits(TIDZ, outhits, maxIBZ, vdecz, vstartz);
205 }
206
207
208 // now outer barrel, r > 55 cm and | z | < 118 cm
209 if(ootantheta < 2.145)
210 {
211 double maxOBR=1200.0;
212 get_which_hits(TOBR, outhits, maxOBR, vdecR, vstartR);
213
214 }
215 else
216 {
217 double maxOBR=1180.0/ootantheta;
218 get_which_hits(TOBR, outhits, maxOBR, vdecR, vstartR);
219
220 // Now check the endcaps
221
222 double maxTEZ=1180.0*ootantheta;
223 get_which_hits(TECZ, outhits, maxTEZ, vdecz, vstartz);
224 }
225
226
227
228 return outhits;
229}
230
231class charged_track {
232 private:
233
234 public:
235 const RecTrackFormat* p; // in future change this to a vector or something else
236 std::vector<bool> _hits;
237
238 charged_track() { }
239
240 charged_track(const RecTrackFormat* q)
241 {
242 p = q;
243 _hits = get_truth_hits(p);
244
245 };
246
247 charged_track(const RecTrackFormat &q, std::mt19937 &engine, std::uniform_real_distribution<double> &rd)
248 {
249 p = &q;
250 std::vector<bool> newhits = get_truth_hits(p);
251
252 // Now model efficiencies: first a per-hit efficiency
253 // 94.5 percent is really for the earlier data; the later data claims 99% although I don't know whether to believe it
254 int nhits=0;
255 for(auto hit : newhits)
256 {
257 if(hit)
258 {
259 double prob = rd(engine);
260 //cout << prob << ", " ;
261 if(prob < 0.945)
262 {
263 _hits.push_back(true);
264 nhits++;
265 //std::cout << "1,";
266 continue;
267 }
268
269 }
270 _hits.push_back(false);
271 //std::cout << "0,";
272 }
273 //cout << std::endl;
274
275 // Next: we know that tracks with a smaller number of total hits
276 // see 1405.6569 page 21, "As a consequence, weaker selection criteria can be applied for tracks having many hit layers"
277 // Apply 70% chance of reconstructing if nhits=4, 80% if nhits=5, 90% if 6 100% if > 6
278
279 if(nhits < 7)
280 {
281 double recoprob=1.0;
282 if(nhits == 4) recoprob=0.7;
283 else if (nhits == 5) recoprob = 0.8;
284 else if (nhits == 6) recoprob = 0.9;
285
286 if(rd(engine) > recoprob)
287 {
288 _hits.clear();
289
290
291 }
292
293 }
294
295
296
297 }
298
299 // for some reason etaCalo is giving 0
300 //double abseta() { return fabs(p->etaCalo()); }
301 double eta() {return p->mc()->momentum().Eta(); }
302 double abseta() {return fabs(this->eta());}
303 //double eta() {return p->etaCalo(); }
304 double Pt() { return p->mc()->momentum().Pt();}
305
306 double E() { return p->mc()->momentum().E(); }
307 // This probably not initialised either
308 //double phi() { return p->phiCalo(); }
309 double phi() { return p->mc()->momentum().Phi(); }
310 int abspid() {return abs(p->mc()->pdgid());}
311 int pid() {return p->mc()->pdgid();}
312 const MALorentzVector& momentum() const { return p->momentum(); }
313
314};
315
316
317
318
319
320// -----------------------------------------------------------------------------
321// Initialize
322// function called one time at the beginning of the analysis
323// -----------------------------------------------------------------------------
324bool cms_exo_19_010::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
325{
326 cout << "BEGIN Initialization" << endl;
327
328 PHYSICS->mcConfig().Reset();
329
330
331 // definition of the multiparticle "hadronic"
332 PHYSICS->mcConfig().AddHadronicId(-20543);
333 PHYSICS->mcConfig().AddHadronicId(-20533);
334 PHYSICS->mcConfig().AddHadronicId(-20523);
335 PHYSICS->mcConfig().AddHadronicId(-20513);
336 PHYSICS->mcConfig().AddHadronicId(-20433);
337 PHYSICS->mcConfig().AddHadronicId(-20423);
338 PHYSICS->mcConfig().AddHadronicId(-20413);
339 PHYSICS->mcConfig().AddHadronicId(-20323);
340 PHYSICS->mcConfig().AddHadronicId(-20313);
341 PHYSICS->mcConfig().AddHadronicId(-20213);
342 PHYSICS->mcConfig().AddHadronicId(-10543);
343 PHYSICS->mcConfig().AddHadronicId(-10541);
344 PHYSICS->mcConfig().AddHadronicId(-10533);
345 PHYSICS->mcConfig().AddHadronicId(-10531);
346 PHYSICS->mcConfig().AddHadronicId(-10523);
347 PHYSICS->mcConfig().AddHadronicId(-10521);
348 PHYSICS->mcConfig().AddHadronicId(-10513);
349 PHYSICS->mcConfig().AddHadronicId(-10511);
350 PHYSICS->mcConfig().AddHadronicId(-10433);
351 PHYSICS->mcConfig().AddHadronicId(-10431);
352 PHYSICS->mcConfig().AddHadronicId(-10423);
353 PHYSICS->mcConfig().AddHadronicId(-10421);
354 PHYSICS->mcConfig().AddHadronicId(-10413);
355 PHYSICS->mcConfig().AddHadronicId(-10411);
356 PHYSICS->mcConfig().AddHadronicId(-10323);
357 PHYSICS->mcConfig().AddHadronicId(-10321);
358 PHYSICS->mcConfig().AddHadronicId(-10313);
359 PHYSICS->mcConfig().AddHadronicId(-10311);
360 PHYSICS->mcConfig().AddHadronicId(-10213);
361 PHYSICS->mcConfig().AddHadronicId(-10211);
362 PHYSICS->mcConfig().AddHadronicId(-5554);
363 PHYSICS->mcConfig().AddHadronicId(-5544);
364 PHYSICS->mcConfig().AddHadronicId(-5542);
365 PHYSICS->mcConfig().AddHadronicId(-5534);
366 PHYSICS->mcConfig().AddHadronicId(-5532);
367 PHYSICS->mcConfig().AddHadronicId(-5524);
368 PHYSICS->mcConfig().AddHadronicId(-5522);
369 PHYSICS->mcConfig().AddHadronicId(-5514);
370 PHYSICS->mcConfig().AddHadronicId(-5512);
371 PHYSICS->mcConfig().AddHadronicId(-5503);
372 PHYSICS->mcConfig().AddHadronicId(-5444);
373 PHYSICS->mcConfig().AddHadronicId(-5442);
374 PHYSICS->mcConfig().AddHadronicId(-5434);
375 PHYSICS->mcConfig().AddHadronicId(-5432);
376 PHYSICS->mcConfig().AddHadronicId(-5424);
377 PHYSICS->mcConfig().AddHadronicId(-5422);
378 PHYSICS->mcConfig().AddHadronicId(-5414);
379 PHYSICS->mcConfig().AddHadronicId(-5412);
380 PHYSICS->mcConfig().AddHadronicId(-5403);
381 PHYSICS->mcConfig().AddHadronicId(-5401);
382 PHYSICS->mcConfig().AddHadronicId(-5342);
383 PHYSICS->mcConfig().AddHadronicId(-5334);
384 PHYSICS->mcConfig().AddHadronicId(-5332);
385 PHYSICS->mcConfig().AddHadronicId(-5324);
386 PHYSICS->mcConfig().AddHadronicId(-5322);
387 PHYSICS->mcConfig().AddHadronicId(-5314);
388 PHYSICS->mcConfig().AddHadronicId(-5312);
389 PHYSICS->mcConfig().AddHadronicId(-5303);
390 PHYSICS->mcConfig().AddHadronicId(-5301);
391 PHYSICS->mcConfig().AddHadronicId(-5242);
392 PHYSICS->mcConfig().AddHadronicId(-5232);
393 PHYSICS->mcConfig().AddHadronicId(-5224);
394 PHYSICS->mcConfig().AddHadronicId(-5222);
395 PHYSICS->mcConfig().AddHadronicId(-5214);
396 PHYSICS->mcConfig().AddHadronicId(-5212);
397 PHYSICS->mcConfig().AddHadronicId(-5203);
398 PHYSICS->mcConfig().AddHadronicId(-5201);
399 PHYSICS->mcConfig().AddHadronicId(-5142);
400 PHYSICS->mcConfig().AddHadronicId(-5132);
401 PHYSICS->mcConfig().AddHadronicId(-5122);
402 PHYSICS->mcConfig().AddHadronicId(-5114);
403 PHYSICS->mcConfig().AddHadronicId(-5112);
404 PHYSICS->mcConfig().AddHadronicId(-5103);
405 PHYSICS->mcConfig().AddHadronicId(-5101);
406 PHYSICS->mcConfig().AddHadronicId(-4444);
407 PHYSICS->mcConfig().AddHadronicId(-4434);
408 PHYSICS->mcConfig().AddHadronicId(-4432);
409 PHYSICS->mcConfig().AddHadronicId(-4424);
410 PHYSICS->mcConfig().AddHadronicId(-4422);
411 PHYSICS->mcConfig().AddHadronicId(-4414);
412 PHYSICS->mcConfig().AddHadronicId(-4412);
413 PHYSICS->mcConfig().AddHadronicId(-4403);
414 PHYSICS->mcConfig().AddHadronicId(-4334);
415 PHYSICS->mcConfig().AddHadronicId(-4332);
416 PHYSICS->mcConfig().AddHadronicId(-4324);
417 PHYSICS->mcConfig().AddHadronicId(-4322);
418 PHYSICS->mcConfig().AddHadronicId(-4314);
419 PHYSICS->mcConfig().AddHadronicId(-4312);
420 PHYSICS->mcConfig().AddHadronicId(-4303);
421 PHYSICS->mcConfig().AddHadronicId(-4301);
422 PHYSICS->mcConfig().AddHadronicId(-4232);
423 PHYSICS->mcConfig().AddHadronicId(-4224);
424 PHYSICS->mcConfig().AddHadronicId(-4222);
425 PHYSICS->mcConfig().AddHadronicId(-4214);
426 PHYSICS->mcConfig().AddHadronicId(-4212);
427 PHYSICS->mcConfig().AddHadronicId(-4203);
428 PHYSICS->mcConfig().AddHadronicId(-4201);
429 PHYSICS->mcConfig().AddHadronicId(-4132);
430 PHYSICS->mcConfig().AddHadronicId(-4122);
431 PHYSICS->mcConfig().AddHadronicId(-4114);
432 PHYSICS->mcConfig().AddHadronicId(-4112);
433 PHYSICS->mcConfig().AddHadronicId(-4103);
434 PHYSICS->mcConfig().AddHadronicId(-4101);
435 PHYSICS->mcConfig().AddHadronicId(-3334);
436 PHYSICS->mcConfig().AddHadronicId(-3324);
437 PHYSICS->mcConfig().AddHadronicId(-3322);
438 PHYSICS->mcConfig().AddHadronicId(-3314);
439 PHYSICS->mcConfig().AddHadronicId(-3312);
440 PHYSICS->mcConfig().AddHadronicId(-3303);
441 PHYSICS->mcConfig().AddHadronicId(-3224);
442 PHYSICS->mcConfig().AddHadronicId(-3222);
443 PHYSICS->mcConfig().AddHadronicId(-3214);
444 PHYSICS->mcConfig().AddHadronicId(-3212);
445 PHYSICS->mcConfig().AddHadronicId(-3203);
446 PHYSICS->mcConfig().AddHadronicId(-3201);
447 PHYSICS->mcConfig().AddHadronicId(-3122);
448 PHYSICS->mcConfig().AddHadronicId(-3114);
449 PHYSICS->mcConfig().AddHadronicId(-3112);
450 PHYSICS->mcConfig().AddHadronicId(-3103);
451 PHYSICS->mcConfig().AddHadronicId(-3101);
452 PHYSICS->mcConfig().AddHadronicId(-2224);
453 PHYSICS->mcConfig().AddHadronicId(-2214);
454 PHYSICS->mcConfig().AddHadronicId(-2212);
455 PHYSICS->mcConfig().AddHadronicId(-2203);
456 PHYSICS->mcConfig().AddHadronicId(-2114);
457 PHYSICS->mcConfig().AddHadronicId(-2112);
458 PHYSICS->mcConfig().AddHadronicId(-2103);
459 PHYSICS->mcConfig().AddHadronicId(-2101);
460 PHYSICS->mcConfig().AddHadronicId(-1114);
461 PHYSICS->mcConfig().AddHadronicId(-1103);
462 PHYSICS->mcConfig().AddHadronicId(-545);
463 PHYSICS->mcConfig().AddHadronicId(-543);
464 PHYSICS->mcConfig().AddHadronicId(-541);
465 PHYSICS->mcConfig().AddHadronicId(-535);
466 PHYSICS->mcConfig().AddHadronicId(-533);
467 PHYSICS->mcConfig().AddHadronicId(-531);
468 PHYSICS->mcConfig().AddHadronicId(-525);
469 PHYSICS->mcConfig().AddHadronicId(-523);
470 PHYSICS->mcConfig().AddHadronicId(-521);
471 PHYSICS->mcConfig().AddHadronicId(-515);
472 PHYSICS->mcConfig().AddHadronicId(-513);
473 PHYSICS->mcConfig().AddHadronicId(-511);
474 PHYSICS->mcConfig().AddHadronicId(-435);
475 PHYSICS->mcConfig().AddHadronicId(-433);
476 PHYSICS->mcConfig().AddHadronicId(-431);
477 PHYSICS->mcConfig().AddHadronicId(-425);
478 PHYSICS->mcConfig().AddHadronicId(-423);
479 PHYSICS->mcConfig().AddHadronicId(-421);
480 PHYSICS->mcConfig().AddHadronicId(-415);
481 PHYSICS->mcConfig().AddHadronicId(-413);
482 PHYSICS->mcConfig().AddHadronicId(-411);
483 PHYSICS->mcConfig().AddHadronicId(-325);
484 PHYSICS->mcConfig().AddHadronicId(-323);
485 PHYSICS->mcConfig().AddHadronicId(-321);
486 PHYSICS->mcConfig().AddHadronicId(-315);
487 PHYSICS->mcConfig().AddHadronicId(-313);
488 PHYSICS->mcConfig().AddHadronicId(-311);
489 PHYSICS->mcConfig().AddHadronicId(-215);
490 PHYSICS->mcConfig().AddHadronicId(-213);
491 PHYSICS->mcConfig().AddHadronicId(-211);
492 PHYSICS->mcConfig().AddHadronicId(111);
493 PHYSICS->mcConfig().AddHadronicId(113);
494 PHYSICS->mcConfig().AddHadronicId(115);
495 PHYSICS->mcConfig().AddHadronicId(130);
496 PHYSICS->mcConfig().AddHadronicId(211);
497 PHYSICS->mcConfig().AddHadronicId(213);
498 PHYSICS->mcConfig().AddHadronicId(215);
499 PHYSICS->mcConfig().AddHadronicId(221);
500 PHYSICS->mcConfig().AddHadronicId(223);
501 PHYSICS->mcConfig().AddHadronicId(225);
502 PHYSICS->mcConfig().AddHadronicId(310);
503 PHYSICS->mcConfig().AddHadronicId(311);
504 PHYSICS->mcConfig().AddHadronicId(313);
505 PHYSICS->mcConfig().AddHadronicId(315);
506 PHYSICS->mcConfig().AddHadronicId(321);
507 PHYSICS->mcConfig().AddHadronicId(323);
508 PHYSICS->mcConfig().AddHadronicId(325);
509 PHYSICS->mcConfig().AddHadronicId(331);
510 PHYSICS->mcConfig().AddHadronicId(333);
511 PHYSICS->mcConfig().AddHadronicId(335);
512 PHYSICS->mcConfig().AddHadronicId(411);
513 PHYSICS->mcConfig().AddHadronicId(413);
514 PHYSICS->mcConfig().AddHadronicId(415);
515 PHYSICS->mcConfig().AddHadronicId(421);
516 PHYSICS->mcConfig().AddHadronicId(423);
517 PHYSICS->mcConfig().AddHadronicId(425);
518 PHYSICS->mcConfig().AddHadronicId(431);
519 PHYSICS->mcConfig().AddHadronicId(433);
520 PHYSICS->mcConfig().AddHadronicId(435);
521 PHYSICS->mcConfig().AddHadronicId(441);
522 PHYSICS->mcConfig().AddHadronicId(443);
523 PHYSICS->mcConfig().AddHadronicId(445);
524 PHYSICS->mcConfig().AddHadronicId(511);
525 PHYSICS->mcConfig().AddHadronicId(513);
526 PHYSICS->mcConfig().AddHadronicId(515);
527 PHYSICS->mcConfig().AddHadronicId(521);
528 PHYSICS->mcConfig().AddHadronicId(523);
529 PHYSICS->mcConfig().AddHadronicId(525);
530 PHYSICS->mcConfig().AddHadronicId(531);
531 PHYSICS->mcConfig().AddHadronicId(533);
532 PHYSICS->mcConfig().AddHadronicId(535);
533 PHYSICS->mcConfig().AddHadronicId(541);
534 PHYSICS->mcConfig().AddHadronicId(543);
535 PHYSICS->mcConfig().AddHadronicId(545);
536 PHYSICS->mcConfig().AddHadronicId(551);
537 PHYSICS->mcConfig().AddHadronicId(553);
538 PHYSICS->mcConfig().AddHadronicId(555);
539 PHYSICS->mcConfig().AddHadronicId(1103);
540 PHYSICS->mcConfig().AddHadronicId(1114);
541 PHYSICS->mcConfig().AddHadronicId(2101);
542 PHYSICS->mcConfig().AddHadronicId(2103);
543 PHYSICS->mcConfig().AddHadronicId(2112);
544 PHYSICS->mcConfig().AddHadronicId(2114);
545 PHYSICS->mcConfig().AddHadronicId(2203);
546 PHYSICS->mcConfig().AddHadronicId(2212);
547 PHYSICS->mcConfig().AddHadronicId(2214);
548 PHYSICS->mcConfig().AddHadronicId(2224);
549 PHYSICS->mcConfig().AddHadronicId(3101);
550 PHYSICS->mcConfig().AddHadronicId(3103);
551 PHYSICS->mcConfig().AddHadronicId(3112);
552 PHYSICS->mcConfig().AddHadronicId(3114);
553 PHYSICS->mcConfig().AddHadronicId(3122);
554 PHYSICS->mcConfig().AddHadronicId(3201);
555 PHYSICS->mcConfig().AddHadronicId(3203);
556 PHYSICS->mcConfig().AddHadronicId(3212);
557 PHYSICS->mcConfig().AddHadronicId(3214);
558 PHYSICS->mcConfig().AddHadronicId(3222);
559 PHYSICS->mcConfig().AddHadronicId(3224);
560 PHYSICS->mcConfig().AddHadronicId(3303);
561 PHYSICS->mcConfig().AddHadronicId(3312);
562 PHYSICS->mcConfig().AddHadronicId(3314);
563 PHYSICS->mcConfig().AddHadronicId(3322);
564 PHYSICS->mcConfig().AddHadronicId(3324);
565 PHYSICS->mcConfig().AddHadronicId(3334);
566 PHYSICS->mcConfig().AddHadronicId(4101);
567 PHYSICS->mcConfig().AddHadronicId(4103);
568 PHYSICS->mcConfig().AddHadronicId(4112);
569 PHYSICS->mcConfig().AddHadronicId(4114);
570 PHYSICS->mcConfig().AddHadronicId(4122);
571 PHYSICS->mcConfig().AddHadronicId(4132);
572 PHYSICS->mcConfig().AddHadronicId(4201);
573 PHYSICS->mcConfig().AddHadronicId(4203);
574 PHYSICS->mcConfig().AddHadronicId(4212);
575 PHYSICS->mcConfig().AddHadronicId(4214);
576 PHYSICS->mcConfig().AddHadronicId(4222);
577 PHYSICS->mcConfig().AddHadronicId(4224);
578 PHYSICS->mcConfig().AddHadronicId(4232);
579 PHYSICS->mcConfig().AddHadronicId(4301);
580 PHYSICS->mcConfig().AddHadronicId(4303);
581 PHYSICS->mcConfig().AddHadronicId(4312);
582 PHYSICS->mcConfig().AddHadronicId(4314);
583 PHYSICS->mcConfig().AddHadronicId(4322);
584 PHYSICS->mcConfig().AddHadronicId(4324);
585 PHYSICS->mcConfig().AddHadronicId(4332);
586 PHYSICS->mcConfig().AddHadronicId(4334);
587 PHYSICS->mcConfig().AddHadronicId(4403);
588 PHYSICS->mcConfig().AddHadronicId(4412);
589 PHYSICS->mcConfig().AddHadronicId(4414);
590 PHYSICS->mcConfig().AddHadronicId(4422);
591 PHYSICS->mcConfig().AddHadronicId(4424);
592 PHYSICS->mcConfig().AddHadronicId(4432);
593 PHYSICS->mcConfig().AddHadronicId(4434);
594 PHYSICS->mcConfig().AddHadronicId(4444);
595 PHYSICS->mcConfig().AddHadronicId(5101);
596 PHYSICS->mcConfig().AddHadronicId(5103);
597 PHYSICS->mcConfig().AddHadronicId(5112);
598 PHYSICS->mcConfig().AddHadronicId(5114);
599 PHYSICS->mcConfig().AddHadronicId(5122);
600 PHYSICS->mcConfig().AddHadronicId(5132);
601 PHYSICS->mcConfig().AddHadronicId(5142);
602 PHYSICS->mcConfig().AddHadronicId(5201);
603 PHYSICS->mcConfig().AddHadronicId(5203);
604 PHYSICS->mcConfig().AddHadronicId(5212);
605 PHYSICS->mcConfig().AddHadronicId(5214);
606 PHYSICS->mcConfig().AddHadronicId(5222);
607 PHYSICS->mcConfig().AddHadronicId(5224);
608 PHYSICS->mcConfig().AddHadronicId(5232);
609 PHYSICS->mcConfig().AddHadronicId(5242);
610 PHYSICS->mcConfig().AddHadronicId(5301);
611 PHYSICS->mcConfig().AddHadronicId(5303);
612 PHYSICS->mcConfig().AddHadronicId(5312);
613 PHYSICS->mcConfig().AddHadronicId(5314);
614 PHYSICS->mcConfig().AddHadronicId(5322);
615 PHYSICS->mcConfig().AddHadronicId(5324);
616 PHYSICS->mcConfig().AddHadronicId(5332);
617 PHYSICS->mcConfig().AddHadronicId(5334);
618 PHYSICS->mcConfig().AddHadronicId(5342);
619 PHYSICS->mcConfig().AddHadronicId(5401);
620 PHYSICS->mcConfig().AddHadronicId(5403);
621 PHYSICS->mcConfig().AddHadronicId(5412);
622 PHYSICS->mcConfig().AddHadronicId(5414);
623 PHYSICS->mcConfig().AddHadronicId(5422);
624 PHYSICS->mcConfig().AddHadronicId(5424);
625 PHYSICS->mcConfig().AddHadronicId(5432);
626 PHYSICS->mcConfig().AddHadronicId(5434);
627 PHYSICS->mcConfig().AddHadronicId(5442);
628 PHYSICS->mcConfig().AddHadronicId(5444);
629 PHYSICS->mcConfig().AddHadronicId(5503);
630 PHYSICS->mcConfig().AddHadronicId(5512);
631 PHYSICS->mcConfig().AddHadronicId(5514);
632 PHYSICS->mcConfig().AddHadronicId(5522);
633 PHYSICS->mcConfig().AddHadronicId(5524);
634 PHYSICS->mcConfig().AddHadronicId(5532);
635 PHYSICS->mcConfig().AddHadronicId(5534);
636 PHYSICS->mcConfig().AddHadronicId(5542);
637 PHYSICS->mcConfig().AddHadronicId(5544);
638 PHYSICS->mcConfig().AddHadronicId(5554);
639 PHYSICS->mcConfig().AddHadronicId(10111);
640 PHYSICS->mcConfig().AddHadronicId(10113);
641 PHYSICS->mcConfig().AddHadronicId(10211);
642 PHYSICS->mcConfig().AddHadronicId(10213);
643 PHYSICS->mcConfig().AddHadronicId(10221);
644 PHYSICS->mcConfig().AddHadronicId(10223);
645 PHYSICS->mcConfig().AddHadronicId(10311);
646 PHYSICS->mcConfig().AddHadronicId(10313);
647 PHYSICS->mcConfig().AddHadronicId(10321);
648 PHYSICS->mcConfig().AddHadronicId(10323);
649 PHYSICS->mcConfig().AddHadronicId(10331);
650 PHYSICS->mcConfig().AddHadronicId(10333);
651 PHYSICS->mcConfig().AddHadronicId(10411);
652 PHYSICS->mcConfig().AddHadronicId(10413);
653 PHYSICS->mcConfig().AddHadronicId(10421);
654 PHYSICS->mcConfig().AddHadronicId(10423);
655 PHYSICS->mcConfig().AddHadronicId(10431);
656 PHYSICS->mcConfig().AddHadronicId(10433);
657 PHYSICS->mcConfig().AddHadronicId(10441);
658 PHYSICS->mcConfig().AddHadronicId(10443);
659 PHYSICS->mcConfig().AddHadronicId(10511);
660 PHYSICS->mcConfig().AddHadronicId(10513);
661 PHYSICS->mcConfig().AddHadronicId(10521);
662 PHYSICS->mcConfig().AddHadronicId(10523);
663 PHYSICS->mcConfig().AddHadronicId(10531);
664 PHYSICS->mcConfig().AddHadronicId(10533);
665 PHYSICS->mcConfig().AddHadronicId(10541);
666 PHYSICS->mcConfig().AddHadronicId(10543);
667 PHYSICS->mcConfig().AddHadronicId(10551);
668 PHYSICS->mcConfig().AddHadronicId(10553);
669 PHYSICS->mcConfig().AddHadronicId(20113);
670 PHYSICS->mcConfig().AddHadronicId(20213);
671 PHYSICS->mcConfig().AddHadronicId(20223);
672 PHYSICS->mcConfig().AddHadronicId(20313);
673 PHYSICS->mcConfig().AddHadronicId(20323);
674 PHYSICS->mcConfig().AddHadronicId(20333);
675 PHYSICS->mcConfig().AddHadronicId(20413);
676 PHYSICS->mcConfig().AddHadronicId(20423);
677 PHYSICS->mcConfig().AddHadronicId(20433);
678 PHYSICS->mcConfig().AddHadronicId(20443);
679 PHYSICS->mcConfig().AddHadronicId(20513);
680 PHYSICS->mcConfig().AddHadronicId(20523);
681 PHYSICS->mcConfig().AddHadronicId(20533);
682 PHYSICS->mcConfig().AddHadronicId(20543);
683 PHYSICS->mcConfig().AddHadronicId(20553);
684 PHYSICS->mcConfig().AddHadronicId(100443);
685 PHYSICS->mcConfig().AddHadronicId(100553);
686 PHYSICS->mcConfig().AddHadronicId(9900440);
687 PHYSICS->mcConfig().AddHadronicId(9900441);
688 PHYSICS->mcConfig().AddHadronicId(9900443);
689 PHYSICS->mcConfig().AddHadronicId(9900551);
690 PHYSICS->mcConfig().AddHadronicId(9900553);
691 PHYSICS->mcConfig().AddHadronicId(9910441);
692 PHYSICS->mcConfig().AddHadronicId(9910551);
693
694 // definition of the multiparticle "invisible"
695 PHYSICS->mcConfig().AddInvisibleId(-16);
696 PHYSICS->mcConfig().AddInvisibleId(-14);
697 PHYSICS->mcConfig().AddInvisibleId(-12);
698 PHYSICS->mcConfig().AddInvisibleId(12);
699 PHYSICS->mcConfig().AddInvisibleId(14);
700 PHYSICS->mcConfig().AddInvisibleId(16);
701 PHYSICS->mcConfig().AddInvisibleId(1000022);
702 PHYSICS->mcConfig().AddInvisibleId(1000039);
703
704 // Initializing PhysicsService for RECO
705 PHYSICS->recConfig().Reset();
706
707 // initialize variables, histos
708
709 cout << "--------------------------------------------------" <<std::endl;
710 cout << "-- CMS disappearing track search 139 fb^-1 --" <<std::endl;
711 cout << "-- arXiv:2004.05153 --" << std::endl;
712 cout << "-- By Mark Goodsell (goodsell@lpthe.jussieu.fr) --" << std::endl;
713 cout << "--------------------------------------------------" <<std::endl;
714 this->engine = std::mt19937(time(NULL));
715 // std::string preliminary[] ={"Njets25 >=2","1 signal lepton","Second baseline lepton veto","mT>50","ET>180","Njets <=3","2 bjets","mbb>50","ET>240","mbb"};
716 //const std::vector<std::string> allSRs={"SR1_2017","SR1_2018A","SR1_2018B","SR2_2017","SR2_2018A","SR2_2018B","SR3_2017","SR3_2018A","SR3_2018B"};
717
718
719 const std::vector<std::string> vecallSRs={"SR3_2015","SR3_2016A","SR3_2016B","SR1_2017","SR1_2018A","SR1_2018B","SR2_2017","SR2_2018A","SR2_2018B","SR3_2017","SR3_2018A","SR3_2018B"};
720 std::string allSRs[]={"SR3_2015","SR3_2016A","SR3_2016B","SR1_2017","SR1_2018A","SR1_2018B","SR2_2017","SR2_2018A","SR2_2018B","SR3_2017","SR3_2018A","SR3_2018B"};
721
722 std::string allSRs2[]={"SR1_2017","SR1_2018A","SR1_2018B","SR2_2017","SR2_2018A","SR2_2018B","SR3_2017","SR3_2018A","SR3_2018B"};
723
724 /*
725 std::string allSR1[]={"SR1_2017","SR1_2018A","SR1_2018B"};
726 std::string allSR2[]={"SR2_2017","SR2_2018A","SR2_2018B"};
727 std::string allSR3[]={"SR3_2017","SR3_2018A","SR3_2018B"};
728 */
729
730 std::string all2016[]={"SR3_2016A","SR3_2016B"};
731 std::string all2017[]={"SR1_2017","SR2_2017","SR3_2017"};
732 std::string all2018A[]={"SR1_2018A","SR2_2018A","SR3_2018A"};
733 std::string all2018B[]={"SR1_2018B","SR2_2018B","SR3_2018B"};
734
735 for(std::string region : vecallSRs)
736 {
737 Manager()->AddRegionSelection(region);
738 }
739
740
741
742
743 Manager()->AddCut("MET2015","SR3_2015");
744 Manager()->AddCut("MET",allSRs2);
745
746 Manager()->AddCut("OneGoodJet",allSRs);
747
748 Manager()->AddCut("JetAngleCuts",allSRs);
749 Manager()->AddCut("|Delta phi(leading jet, pTmiss)| > 0.5",allSRs);
750
751 Manager()->AddCut(">=1 track with |eta| < 2.1",allSRs);
752 Manager()->AddCut(">=1 track with p_T > 55",allSRs);
753 Manager()->AddCut(">=1 track passing fiducial selections",allSRs);
754 Manager()->AddCut(">=1 track with > 3 or 4 pixel hits",allSRs);
755 Manager()->AddCut(">=1 track with no missing inner/middle hits",allSRs);
756 Manager()->AddCut(">=1 track with relative isolation",allSRs);
757 Manager()->AddCut(">=1 track with d0 < 0.2 mm",allSRs);
758 Manager()->AddCut(">=1 track with dz < 5.0 mm",allSRs);
759 Manager()->AddCut(">=1 track with DR(track, jet) > 0.5",allSRs);
760 Manager()->AddCut(">=1 track with DR(track, electron) > 0.15",allSRs);
761 Manager()->AddCut(">=1 track with DR(track, muon) > 0.15",allSRs);
762 Manager()->AddCut(">=1 track with DR(track, tau) > 0.15",allSRs);
763 Manager()->AddCut(">=1 track with Ecalo < 10",allSRs);
764 Manager()->AddCut(">=1 track with >=3 missing outer hits",allSRs);
765
766 Manager()->AddCut("HEMVeto",all2018B);
767
768
769
770 Manager()->AddCut(">5 Layers 2015","SR3_2015");
771 // Manager()->AddCut(">5 Layers 2016","SR3_2016");
772 Manager()->AddCut(">5 Layers 2016",all2016);
773
774 Manager()->AddCut("4 Layers 2017","SR1_2017");
775 Manager()->AddCut("4 Layers 2018A","SR1_2018A");
776 Manager()->AddCut("4 Layers 2018B","SR1_2018B");
777
778 Manager()->AddCut("5 Layers 2017","SR2_2017");
779 Manager()->AddCut("5 Layers 2018A","SR2_2018A");
780 Manager()->AddCut("5 Layers 2018B","SR2_2018B");
781
782 Manager()->AddCut(">5 Layers 2017","SR3_2017");
783 Manager()->AddCut(">5 Layers 2018A","SR3_2018A");
784 Manager()->AddCut(">5 Layers 2018B","SR3_2018B");
785
786 Manager()->AddCut("Reweight 2015","SR3_2015");
787 Manager()->AddCut("Reweight 2016A","SR3_2016A");
788 Manager()->AddCut("Reweight 2016B","SR3_2016B");
789 Manager()->AddCut("Reweight 2017",all2017);
790 Manager()->AddCut("Reweight 2018A",all2018A);
791 Manager()->AddCut("Reweight 2018B",all2018B);
792
793 double totallumi=2.7+8.3+27.4+42.0+21.0+39.0;
794
795 lumiratios["2015"] = 2.7/totallumi;
796 lumiratios["2016A"] = 8.3/totallumi;
797 lumiratios["2016B"] = 27.4/totallumi;
798 lumiratios["2017"] = 42.0/totallumi;
799 lumiratios["2018A"] = 21.0/totallumi;
800 lumiratios["2018B"] = 39.0/totallumi;
801
802 cout << "END Initialization" << endl;
803 return true;
804}
805
806// -----------------------------------------------------------------------------
807// Finalize
808// function called one time at the end of the analysis
809// -----------------------------------------------------------------------------
810void cms_exo_19_010::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
811{
812 cout << "BEGIN Finalization" << endl;
813 // saving histos
814 cout << "END Finalization" << endl;
815}
816
817// -----------------------------------------------------------------------------
818// Execute
819// function called each time one event is read
820// -----------------------------------------------------------------------------
821bool cms_exo_19_010::Execute(SampleFormat& sample, const EventFormat& event)
822{
823
824
825 double myWeight=1.;
826 if (!Configuration().IsNoEventWeight()) myWeight=event.mc()->weight();
827 else if(event.mc()->weight()!=0.) myWeight=event.mc()->weight();
828 else
829 {
830 //////WARNING << "Found one event with a zero weight. Skipping...\n";
831 return false;
832 }
833
834
835 Manager()->InitializeForNewEvent(myWeight);
836
837 // The event loop starts here
838 if(event.rec() ==0) return false;
839
840 MALorentzVector pTmiss,pTmissNoMu;
841
842 std::vector<const RecTrackFormat*> charginos;
843
844 // Count signal electrons, muons, jets etc
845
846
847static std::uniform_real_distribution<double> rd(0.0,1.0);
848bool chargino50=false;
849
850std::vector<charged_track*> tracks;
851
852MALorentzVector muonp;
853MALorentzVector charginop;
854 muonp.clear();
855 charginop.clear();
856
857 for(auto muon : event.rec()->muons())
858 {
859 muonp+=muon.momentum();
860
861 }
862
863double maxtrackpt=0.0;
864double chpt=0.0;
865
866// Now need to fill up the candidate charged tracks ...
867
868 for(auto &track : event.rec()->tracks())
869 {
870
871 const MCParticleFormat* chargino = track.mc();
872
873 // NOT for tracks that are final states: MA5 will have already included them in the ET calculation
874 //if ( (PHYSICS->Id->IsFinalState(chargino)) || (chargino->decay_vertex().Pt() > 7.0e3) || (fabs(chargino->decay_vertex().Pz()) > 1.1e4))
875 if ( (chargino->decay_vertex().Pt() > 7.0e3) || (fabs(chargino->decay_vertex().Pz()) > 1.1e4))
876 {
877 if(!(PHYSICS->Id->IsHadronic(chargino)))
878 {
879 // gets reconstructed as a muon
880 charginop+=chargino->momentum();
881 //new_fake_muons.push_back(chargino);
882 }
883 }
884
885 chpt=chargino->momentum().Pt();
886
887 if(chpt > maxtrackpt) maxtrackpt=chpt;
888
889 //if(fabs(track.etaCalo()) < 2.5)
890 if(fabs(chargino->momentum().Eta()) < 2.5)
891 {
892 charged_track* newtrack = new charged_track(track,engine,rd);
893 tracks.push_back(newtrack);
894
895 if((!chargino50) && (chargino->momentum().Pt() >50.0)) chargino50=true;
896
897
898 }
899
900
901
902 }
903
904 //////// NOW DO MET CALC!
905
906 //pTmiss = event.rec()->MET().momentum();
907 // pTmissNoMu=pTmiss+muonp;
908
909 pTmissNoMu = event.rec()->MET().momentum();
910
911 pTmiss=pTmissNoMu-charginop;
912
913 pTmissNoMu = pTmissNoMu +muonp;
914 double MET=pTmiss.Pt();
915 double METnomu=pTmissNoMu.Pt();
916
917
918
919 bool passMET=false;
920 bool passMET2015=false;
921
922 if((chargino50) && (MET > 105.0)) passMET=true;
923 if(MET > 120.0) passMET = true;
924 if(METnomu < 120.0)
925 {
926 passMET = false;
927 }
928 else
929 {
930 passMET = true;
931 }
932
933 if((chargino50) && (MET > 75.0) ) passMET2015=true;
934 if((MET > 90.0) || (METnomu > 90.0)) passMET2015 = true;
935 if(MET < 100.0) passMET2015 = false;
936
937
938 double eff=1.0;
939 // 1903.06078 figure 5, since the first cut on MET is the trigger. It uses HLT with 90 GeV or 120 GeV
940 if(passMET2015)
941 {
942
943 double tMET=MET;
944 if (METnomu > MET) tMET=METnomu;
945 if(tMET < 90.0)
946 {
947 eff=0.0;
948 }
949 else if (tMET < 250.0)
950 {
951 eff = (tMET-50.0)/200.0;
952 }
953 else
954 {
955 eff=1.0;
956 }
957
958 if(rd(engine) > eff) passMET2015=false;
959 }
960
961 if(passMET)
962 {
963 // Use METnomu?
964
965 double tMET=MET;
966 if (METnomu > MET) tMET=METnomu;
967
968 if(tMET < 120.0)
969 {
970 eff=0.0;
971 }
972 else if (tMET < 280.0)
973 {
974 eff = (tMET-120.0)/160.0;
975 }
976 else
977 {
978 eff=1.0;
979 }
980
981 if(rd(engine) > eff) passMET=false;
982 }
983
984 if(!(Manager()->ApplyCut(passMET,"MET"))) return true;
985 if(!(Manager()->ApplyCut(passMET2015,"MET2015"))) return true;
986
987 if((!passMET) && (!passMET2015)) {for(auto track: tracks) {delete track;}; return true;} ;
988
989
990
991 bool onegoodjet = false;
992 bool JetMetDeltaphi=true;
993 bool passJetAngleCuts = true;
994
995 if(event.rec()->jets().size() ==0 ) {for(auto track: tracks) {delete track;}; return true;} ;
996
997 std::vector<const RecJetFormat*> Jets;
998 for (int i=0; i < event.rec()->jets().size(); i++)
999 {
1000 const RecJetFormat *CurrentJet = &(event.rec()->jets()[i]);
1001 Jets.push_back(CurrentJet);
1002 }
1003
1004 SORTER->sort(Jets);
1005 std::vector<const RecJetFormat*> GoodJets;
1006
1007 if(fabs(Jets[0]->momentum().DeltaPhi(pTmiss)) < 0.5) JetMetDeltaphi = false;
1008
1009 filterPhaseSpace(Jets,30.0,2.4);
1010
1011 for(int m = 0; m< Jets.size() ; m++)
1012 {
1013
1014 const RecJetFormat* tjet=Jets[m];
1015
1016
1017 if((!onegoodjet) && (fabs(tjet->eta())< 2.4) && (tjet->pt() > 110.0))
1018 {
1019 onegoodjet=true;
1020 }
1021 if(tjet->pt() > 30.0)
1022 {
1023 GoodJets.push_back(tjet);
1024
1025 }
1026 for(int n=m+1; n < Jets.size(); n++)
1027 {
1028
1029 if(fabs(tjet->momentum().DeltaPhi(Jets[n]->momentum())) > 2.5)
1030 {
1031 passJetAngleCuts = false;
1032 break;
1033 }
1034
1035 }
1036
1037 }
1038
1039
1040 if(!(Manager()->ApplyCut(onegoodjet,"OneGoodJet"))) return true;
1041 if(!onegoodjet) {for(auto track: tracks) {delete track;}; return true;} ;
1042
1043
1044 if(!(Manager()->ApplyCut(passJetAngleCuts,"JetAngleCuts"))) return true;
1045 if(!(Manager()->ApplyCut(JetMetDeltaphi,"|Delta phi(leading jet, pTmiss)| > 0.5"))) return true;
1046 if((!passJetAngleCuts) || (!JetMetDeltaphi)) {for(auto track: tracks) {delete track;}; return true;} ;
1047
1048
1049
1050 if(tracks.size() < 1) return true;
1051 bool good2015_3=false;
1052 bool good2016_3=false;
1053 bool good2017_1=false;
1054 bool good2017_2=false;
1055 bool good2017_3=false;
1056 bool good2018A_1=false;
1057 bool good2018A_2=false;
1058 bool good2018A_3=false;
1059 bool good2018B_1=false;
1060 bool good2018B_2=false;
1061 bool good2018B_3=false;
1062
1063 bool tracketa=false;
1064 bool trackpT=false;
1065 bool trackfiducial=false;
1066 bool trackpixel=false;
1067 bool trackinnermiddle=false;
1068 bool trackisolation=false;
1069 bool trackd0=false;
1070 bool trackdz=false;
1071 bool trackDRjet=false;
1072
1073
1074 std::vector<charged_track*> tracks2;
1075
1076 for(auto chargino: tracks)
1077 {
1078 if((chargino->abseta() > 2.1) ) {delete chargino; continue;};
1079
1080 tracketa=true;
1081
1082
1083 if(chargino->Pt() < 55.0) {delete chargino; continue;};
1084 trackpT=true;
1085
1086
1087
1088 double cheta=chargino->abseta();
1089
1090
1091 // apply conditions on the track to avoid low efficiency parts of muon chamber/ECAL
1092 // muon chamber
1093 //if((decayrad > 4020.0) && (decayrad < 7380.0))
1094 //{
1095 if((cheta >0.15) && (cheta < 0.35) ) {delete chargino; continue;};
1096 if((cheta >1.55) && (cheta < 1.85) ) {delete chargino; continue;};
1097 //}
1098 // ECAL
1099 //if((decayrad > 1290.0) && (decayrad < 1790.0))
1100 //{
1101 if((cheta >1.42) && (cheta < 1.65) ) {delete chargino; continue;};
1102
1103 trackfiducial=true;
1104
1105
1106 // Now we have >=1 track with > 3 or 4 pixel hits, actually this means we have to have all 4 closest hits
1107
1108 if(chargino->_hits.size() < 4) {delete chargino; continue;};
1109
1110 bool charginotrackpixel=true;
1111
1112 for(int i=0; i<4; i++)
1113 {
1114 if(!chargino->_hits[i])
1115 {
1116 charginotrackpixel=false;
1117 break;
1118 }
1119 }
1120 if(!charginotrackpixel) {delete chargino; continue;};
1121 trackpixel=true;
1122 bool charginoinnermiddle=true;
1123 // now check for missing inner/middle hits
1124 if(chargino->_hits.size() > 4)
1125 {
1126 bool foundoutside=false;
1127 for(int i=chargino->_hits.size() -1; i>3; i--)
1128 {
1129 if(foundoutside)
1130 {
1131 if(!chargino->_hits[i])
1132 {
1133 charginoinnermiddle=false;
1134 break;
1135 }
1136
1137 }
1138 else
1139 {
1140 if(chargino->_hits[i])
1141 {
1142 foundoutside=true;
1143 }
1144 }
1145 }
1146
1147 }
1148
1149 if(!charginoinnermiddle) {delete chargino; continue;};
1150
1151 trackinnermiddle=true;
1152
1153 double isopt=fabs(chargino->p->isolCones()[0].sumPT());
1154 //double isopt=chargino->p->isolCones()[0].sumPT()-chargino->Pt(); // MA5 isolation includes its own ...
1155 //double isopt = sumpTisolation(chargino->p->momentum(),Event,0.3);
1156
1157 if( (isopt/chargino->Pt()) > 0.05) {delete chargino; continue;};
1158
1159 trackisolation=true;
1160
1161 double absd0=fabs(chargino->p->d0());
1162 double absdz=fabs(chargino->p->dz());
1163
1164 if(absd0 > 0.2) {delete chargino; continue;};
1165 trackd0=true;
1166 if(absdz > 5.0) {delete chargino; continue;};
1167 trackdz=true;
1168
1169 tracks2.push_back(chargino);
1170 }
1171
1172 if(!(Manager()->ApplyCut(tracketa,">=1 track with |eta| < 2.1"))) return true;
1173 if(!(Manager()->ApplyCut(trackpT,">=1 track with p_T > 55"))) return true;
1174 if(!(Manager()->ApplyCut(trackfiducial,">=1 track passing fiducial selections"))) return true;
1175 if(!(Manager()->ApplyCut(trackpixel,">=1 track with > 3 or 4 pixel hits"))) return true;
1176 if(!(Manager()->ApplyCut(trackinnermiddle,">=1 track with no missing inner/middle hits"))) return true;
1177 if(!(Manager()->ApplyCut(trackisolation,">=1 track with relative isolation"))) return true;
1178 if(!(Manager()->ApplyCut(trackd0,">=1 track with d0 < 0.2 mm"))) return true;
1179 if(!(Manager()->ApplyCut(trackdz,">=1 track with dz < 5.0 mm"))) return true;
1180
1181
1182
1183if(tracks2.size() < 1) {for(auto track: tracks2) {delete track;}; return true;} ;
1184
1185
1186
1187tracks=FullRemoval(tracks2,GoodJets,0.5);
1188 if(!(Manager()->ApplyCut((tracks.size() > 0),">=1 track with DR(track, jet) > 0.5"))) return true;
1189if(tracks.size() ==0) return true;
1190 tracks=FullRemoval(tracks,event.rec()->electrons(),0.15);
1191 if(!(Manager()->ApplyCut((tracks.size() > 0),">=1 track with DR(track, electron) > 0.15"))) return true;
1192if(tracks.size() ==0) return true;
1193
1194// This cut should remove long-lived charginos which get reconstructed as a muon
1195tracks=FullRemoval(tracks,event.rec()->muons(),0.15);
1196
1197std::vector<charged_track*> tracks3;
1198std::vector<charged_track*> fake_muons;
1199 for(auto chargino: tracks)
1200 {
1201 MALorentzVector vdec = chargino->p->mc()->decay_vertex();
1202 //if((vdec.pT() > 4000.0) || (fabs(vdec.pz()) > 6000.0 )) {delete chargino; continue;}
1203 //if((vdec.pT() > 7000.0) || (fabs(vdec.pz()) > 11000.0 )) {delete chargino; continue;}
1204 // if it decays outside muon chamber and is not hadronic (i.e. caught in the calo)
1205 if((!(PHYSICS->Id->IsHadronic(chargino->p->mc()))) &&( (vdec.Pt() > 7000.0) || (fabs(vdec.Pz()) > 11000.0 )))
1206 {
1207 fake_muons.push_back(chargino);
1208 }
1209 else
1210 {
1211 tracks3.push_back(chargino);
1212 }
1213 }
1214 //tracks=tracks3;
1215 tracks=FullRemoval(tracks3,fake_muons,0.15);
1216 for(auto fmu : fake_muons)
1217 {
1218 delete fmu;
1219 }
1220
1221 if(!(Manager()->ApplyCut((tracks.size() > 0),">=1 track with DR(track, muon) > 0.15"))) return true;
1222if(tracks.size() ==0) return true;
1223 tracks=FullRemoval(tracks,event.rec()->taus(),0.15);
1224 if(!(Manager()->ApplyCut((tracks.size() > 0),">=1 track with DR(track, tau) > 0.15"))) return true;
1225
1226if(tracks.size() ==0) return true;
1227
1228bool trackEcalo=false;
1229bool trackmissingouter=false;
1230
1231
1232
1233 for(auto chargino: tracks)
1234 {
1235
1236 trackEcalo=true;
1237
1238 /*
1239 Now we do the Ecalo, unlike in hackanalysis. The Ecalo is supposed to be total calorimeter
1240 energy, rather than sumET, but MA5 v1.9 doesn't have that option, so we do what we can.
1241 I will also include a cut on whether the track is a hadron or electron/tau, under the assumption that
1242 in those cases there is a large energy deposit (we should already have thrown it out
1243 under those cases but ok ...
1244
1245 */
1246
1247 if (PHYSICS->Id->IsHadronic(chargino->p->mc()))
1248 {
1249 trackEcalo=false;
1250 continue;
1251 };
1252
1253
1254 if(chargino->abspid() < 25)
1255 {
1256 trackEcalo=false;
1257 continue;
1258 }
1259 double isoET=fabs(chargino->p->isolCones()[1].sumET());
1260
1261 if(isoET > 10.0)
1262 {
1263 trackEcalo=false;
1264 continue;
1265 }
1266
1267
1268
1269 double chphi=chargino->phi();
1270 double cheta=chargino->eta();
1271 double acheta=fabs(cheta);
1272
1273 bool chtrackmissingouter=true;
1274 int nhits=(chargino->_hits.size())-1;
1275 for(int i=0; i<3; i++)
1276 {
1277 if(chargino->_hits[nhits-i]) // 3 outermost hits must be missing
1278 {
1279 chtrackmissingouter=false;
1280 break;
1281 }
1282
1283 }
1284 if(!chtrackmissingouter) continue;
1285 trackmissingouter=true;
1286 bool good2017 = true;
1287 bool good2018 = true;
1288
1289
1290
1291
1292 if((cheta >0.0) && (cheta < 1.42) && (chphi < 3.142) && (chphi > 2.7) ) good2017=false;
1293 if((cheta >0.0) && (cheta < 1.42) && (chphi < 0.8) && (chphi > 0.4) ) good2018=false;
1294
1295 // Now we've already selected tracks that have no missing hits from start to end, so just count layers
1296
1297 int nlayers=0;
1298
1299 for(auto hit : chargino->_hits)
1300 {
1301 if(hit) nlayers++;
1302 }
1303
1304 if(nlayers == 4)
1305 {
1306 if(good2017) good2017_1=true;
1307 if(good2018)
1308 {
1309 good2018A_1=true;
1310 good2018B_1=true;
1311 }
1312 }
1313 else if(nlayers == 5)
1314 {
1315 if(good2017) good2017_2=true;
1316 if(good2018)
1317 {
1318 good2018A_2=true;
1319 good2018B_2=true;
1320 }
1321
1322 }
1323 else // nlayers > 5
1324 {
1325 if(nlayers > 7)
1326 {
1327 good2015_3=true;
1328 good2016_3=true;
1329 }
1330 if(good2017) good2017_3=true;
1331
1332 if(good2018)
1333 {
1334 good2018A_3=true;
1335 good2018B_3=true;
1336 }
1337 }
1338
1339 }
1340
1341 if(!(Manager()->ApplyCut(trackEcalo,">=1 track with Ecalo < 10"))) return true;
1342 if(!(Manager()->ApplyCut(trackmissingouter,">=1 track with >=3 missing outer hits"))) return true;
1343
1344
1345double phipTmiss = pTmiss.Phi();
1346bool GoodpTPhi=true;
1347if((phipTmiss > -1.6) && (phipTmiss <-0.6)) GoodpTPhi=false;
1348
1349 if(!(Manager()->ApplyCut(GoodpTPhi,"HEMVeto"))) return true;
1350//if(!GoodpTPhi) return;
1351
1352
1353
1354 if(!(Manager()->ApplyCut(good2017_1,"4 Layers 2017"))) return true;
1355 if(!(Manager()->ApplyCut(good2018A_1,"4 Layers 2018A"))) return true;
1356 if(!(Manager()->ApplyCut(good2018B_1,"4 Layers 2018B"))) return true;
1357
1358 if(!(Manager()->ApplyCut(good2017_2,"5 Layers 2017"))) return true;
1359 if(!(Manager()->ApplyCut(good2018A_2,"5 Layers 2018A"))) return true;
1360 if(!(Manager()->ApplyCut(good2018B_2,"5 Layers 2018B"))) return true;
1361
1362 if(!(Manager()->ApplyCut(good2015_3,">5 Layers 2015"))) return true;
1363 if(!(Manager()->ApplyCut(good2016_3,">5 Layers 2016"))) return true;
1364 if(!(Manager()->ApplyCut(good2017_3,">5 Layers 2017"))) return true;
1365 if(!(Manager()->ApplyCut(good2018A_3,">5 Layers 2018A"))) return true;
1366 if(!(Manager()->ApplyCut(good2018B_3,">5 Layers 2018B"))) return true;
1367for(auto track: tracks) {delete track;};
1368
1369
1370// Now apply reweightings
1371
1372//Manager()->SetCurrentEventWeight(eventWeight);
1373
1374 const std::string rwghtname="Reweight ";
1375 for(auto period : alldataperiods)
1376 {
1377 Manager()->SetCurrentEventWeight(myWeight*lumiratios[period]);
1378 if(!(Manager()->ApplyCut(true,rwghtname+period))) return true;
1379 }
1380
1381 /*
1382 Manager()->SetCurrentEventWeight(myWeight*lumiratios["2015"]);
1383 if(!(Manager()->ApplyCut(true,"Reweight 2015"))) return true;
1384
1385 Manager()->SetCurrentEventWeight(myWeight*lumiratios["2016A"]);
1386 if(!(Manager()->ApplyCut(true,"Reweight 2016A"))) return true;
1387
1388 Manager()->SetCurrentEventWeight(myWeight*lumiratios["2016B"]);
1389 if(!(Manager()->ApplyCut(true,"Reweight 2016B"))) return true;
1390
1391 Manager()->SetCurrentEventWeight(myWeight*lumiratios["2017"]);
1392 if(!(Manager()->ApplyCut(true,"Reweight 2016B"))) return true;
1393 */
1394
1395 return true;
1396
1397}
1398