14 namespace ColliderBit {
19 class Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb :
public Analysis {
23 static constexpr
const char* detector =
"ATLAS";
26 #define CUTFLOWMAP(X) \ 28 X(electron_eq_4_jets) \ 29 X(electron_met_gt_40) \ 30 X(electron_eq_2_bjets) \ 44 const std::vector<std::string> cutflowNames = {
CUTFLOWMAP(
g)};
57 const std::vector<std::string> varNames = {
VARMAP(
g)};
63 std::map<std::string, std::vector<double>> varResults;
65 std::map<std::string, int> cutflows;
71 std::vector<double> calcNuPz(
double Mw, P4 metMom, P4 lepMom)
73 double mu =
sqr(Mw) / 2.0 + metMom.px() * lepMom.px() + metMom.py() * lepMom.py();
74 double denom = lepMom.E2() - lepMom.pz2();
75 double a = mu * lepMom.pz() / denom;
77 double b = (lepMom.E2() * metMom.E2() -
sqr(mu)) / denom;
85 return {a + std::sqrt(delta), a - std::sqrt(delta)};
89 P4 getBestHadronicTop(
90 std::vector<const Jet *> bJets,
91 std::vector<const Jet *> lightJets,
99 auto prob = [&width, &mean](P4 particle) {
100 return 1 - std::erf(1.0 * std::abs(particle.m() - mean) / (std::sqrt(2.0) * width));
105 std::vector<double> nuPzChoices = calcNuPz(80.0, metMom, leptonMom);
107 for (
double nuPz : nuPzChoices)
109 double nuE = std::sqrt(
sqr(metMom.px()) +
sqr(metMom.py()) +
sqr(nuPz));
110 nu.setPE(metMom.px(), metMom.py(), nuPz, nuE);
111 P4 WLep = leptonMom + nu;
113 for (
const Jet* firstBJet : bJets)
115 for (
const Jet* secondBJet : bJets)
117 if (firstBJet == secondBJet)
121 P4 topLep = *firstBJet + WLep;
123 for (
const Jet* firstLightJet : lightJets)
125 for (
const Jet* secondLightJet : lightJets)
128 if (firstLightJet == secondLightJet)
132 P4 WHad = *firstLightJet + *secondLightJet;
133 P4 topHad = *secondBJet + WHad;
135 double newPTotal = prob(topHad) * prob(WHad) * prob(topLep) * prob(WLep);
136 if (newPTotal > pTotal)
140 bestHadronicTop = topHad;
147 return bestHadronicTop;
150 double calcMt(P4 metVec, P4 lepVec)
152 double Met = metVec.pT();
153 double pTLep = lepVec.pT();
158 double calcSqrtSSubMin(P4 visibleSubsystem, P4 invisbleSubsystem)
160 double visiblePart = std::sqrt(
sqr(visibleSubsystem.m()) +
sqr(visibleSubsystem.pT()));
161 double invisiblePart = invisbleSubsystem.pT();
162 double twoDimensionalVecSum =
163 sqr(visibleSubsystem.px() + invisbleSubsystem.px())
164 +
sqr(visibleSubsystem.py() + invisbleSubsystem.py());
165 return std::sqrt(
sqr(visiblePart + invisiblePart) - twoDimensionalVecSum);
169 std::vector<const Jet*>& jets,
170 std::vector<const Jet*>* bJets,
171 std::vector<const Jet*>* lightJets)
174 const std::vector<double> a = {0, 10.};
175 const std::vector<double>
b = {0, 10000.};
176 const std::vector<double> c = {0.60};
177 BinnedFn2D<double> _eff2d(a,b,c);
178 for (
const Jet* jet : jets)
180 bool hasTag =
has_tag(_eff2d, jet->eta(), jet->pT());
181 if(jet->btag() && hasTag && jet->abseta() < 2.5)
183 bJets->push_back(jet);
187 lightJets->push_back(jet);
196 Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb()
198 set_analysis_name(
"ATLAS_7TeV_1OR2LEPStop_4_7invfb");
206 void run(
const HEPUtils::Event* event)
211 incrementCut(Total_events);
212 std::vector<const Particle*> electrons =
event->electrons();
213 std::vector<const Particle*> muons =
event->muons();
214 std::vector<const Jet*> jets =
event->jets();
220 std::vector<const Jet*> bJets, lightJets;
221 getBJets(jets, &bJets, &lightJets);
239 nLeptons = leptons.size(),
241 nBJets = bJets.size(),
242 nLightJets = lightJets.size();
244 double Met =
event->met();
245 const P4& metVec =
event->missingmom();
254 if (electrons.size() == 1) a = electron_eq_4_jets, b = electron_met_gt_40, c = electron_eq_2_bjets;
255 if (muons.size() == 1) a = muon_eq_4_jets, b = muon_met_gt_40, c = muon_eq_2_bjets;
276 if (nLeptons == 1 && nBJets >= 2 && nLightJets >= 2 && Met > 40)
278 double mT = calcMt(metVec, leptons[0]->mom());
280 auto isValidTop = [](
double mean,
double width,
double mass) {
return mass < mean - 0.5 * width;};
284 double mean = 0.0, width = 0.0;
287 P4 visibleSubsystem = *leptons[0] + *lightJets[0] + *lightJets[1] + *bJets[0] + *bJets[1];
288 double sqrtSsubMin = calcSqrtSSubMin(visibleSubsystem, metVec);
289 bool isOneLep =
false;
291 if (electrons.size() == 1 && electrons[0]->pT() > 25)
293 mean = 168.4, width = 18.0;
294 hadronicTop = getBestHadronicTop(bJets, lightJets, *electrons[0], metVec, width, mean);
299 if (muons.size() == 1 && muons[0]->pT() > 20)
301 mean = 168.2, width = 18.6;
302 hadronicTop = getBestHadronicTop(bJets, lightJets, *muons[0], metVec, width, mean);
307 bool validTop = isValidTop(mean, width, hadronicTop.m());
311 varResults[varNames[mTopHad]].push_back(hadronicTop.m());
312 varResults[varNames[oneLepSqrtS]].push_back(sqrtSsubMin);
316 if (isOneLep && validTop && sqrtSsubMin < 250)
318 num1LSR +=
event->weight();
319 if (electrons.size() == 1) incrementCut(electron_sr);
320 if (muons.size() == 1) incrementCut(muon_sr);
327 P4 ll = *leptons[0] + *leptons[1];
329 incrementCut(twoLep_met_gt_40);
333 incrementCut(twoLep_gt_1_bjet);
336 incrementCut(mll_lt_81);
338 if (mll < 81 && mll > 30)
340 incrementCut(mll_gt_30_lt_81);
348 P4 ll = *leptons[0] + *leptons[1];
351 P4 visibleSubsystem = *leptons[0] + *leptons[1] + *jets[0] + *jets[1];
353 double sqrtSsubMin = calcSqrtSSubMin(visibleSubsystem, metVec);
355 double mlljj = visibleSubsystem.m();
356 varResults[varNames[mllKey]].push_back(mll);
357 if (mll > 30 && mll < 81)
359 bool isTwoLeptonEvent =
false;
362 if (electrons.size() == 2 && electrons[0]->pT() > 25)
364 isTwoLeptonEvent =
true;
368 if (muons.size() == 2 && muons[0]->pT() > 20)
370 isTwoLeptonEvent =
true;
374 if (electrons.size() == 1 && muons.size() == 1 && (electrons[0]->pT() > 25 || muons[0]->pT() > 20))
376 isTwoLeptonEvent =
true;
379 if (isTwoLeptonEvent)
381 varResults[varNames[twoLepSqrtS]].push_back(sqrtSsubMin);
383 if (sqrtSsubMin < 225)
385 num2LSR1 +=
event->weight();
386 incrementCut(num_2lsr1);
388 if (sqrtSsubMin < 235 && mlljj < 140)
390 num2LSR2 +=
event->weight();
391 incrementCut(num_2lsr2);
403 void combine(
const Analysis* other)
406 const Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb* specificOther =
dynamic_cast<const Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb*
>(other);
407 num1LSR += specificOther->num1LSR;
408 num2LSR1 += specificOther->num2LSR1;
409 num2LSR2 += specificOther->num2LSR2;
413 void collect_results()
417 add_result(SignalRegionData(
"1LSR", 50, {num1LSR, 0.}, {38., 7.}));
418 add_result(SignalRegionData(
"2LSR1", 123, {num2LSR1, 0.}, {115., 15.}));
419 add_result(SignalRegionData(
"2LSR2", 47, {num2LSR2, 0.}, {46., 7.}));
436 void analysis_specific_reset()
441 for (std::string varName : varNames)
443 varResults[varName] = {};
463 void incrementCut(
int cutIndex)
465 cutflows[cutflowNames[cutIndex]]++;
470 double scale_by = 1.0;
471 cout <<
"SAVE_START:cuts" << endl;
472 cout <<
"CUT;RAW;SCALED;%" << endl;
473 double initialCut = cutflows[cutflowNames[Total_events]];
475 for (std::string name : cutflowNames) {
476 thisCut = cutflows[name];
477 cout << name.c_str() <<
";" 479 << thisCut * scale_by <<
";" 480 << 100. * thisCut / initialCut
483 cout <<
"SAVE_END" << endl;
static std::vector< const HEPUtils::Jet * > filterPtEta(std::vector< const HEPUtils::Jet *> jets, double pT, double absEta)
Utils for ColliderBit analyses.
bool has_tag(const HEPUtils::BinnedFn2D< double > &effmap, double eta, double pt)
Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map.
static std::vector< const HEPUtils::Jet * > jetLeptonOverlapRemoval(std::vector< const HEPUtils::Jet *> jets, std::vector< const HEPUtils::Particle *> leptons, double dR)
static bool sortParticlesByPt(const HEPUtils::Particle *particle1, const HEPUtils::Particle *particle2)
Funk delta(std::string arg, double pos, double width)
static bool muonFilter7TeV(const std::vector< const HEPUtils::Particle *> &muons)
void run(bool, bool, bool, int, double, double, int, int, int, int, int, double, char[], int, int[], bool, bool, bool, bool, double, int, double(*)(double *, int, int, void *), void(*)(int, int, int, double *, double *, double *, double, double, double, void *), void *)
static bool sortJetsByPt(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
EXPORT_SYMBOLS double sqr(double a)
returns square of double - saves tedious repetition
static std::vector< const HEPUtils::Jet * > filterMaxEta(const std::vector< const HEPUtils::Jet *> &jets, double maxAbsEta)
The Cutflow and Cutflows classes.
static std::vector< const HEPUtils::Particle * > getSortedLeptons(const std::vector< std::vector< const HEPUtils::Particle *>> allLeptons)
void applyTightIDElectronSelection(std::vector< const HEPUtils::Particle *> &electrons)
Efficiency function for Tight ID electrons.
static double dot2D(HEPUtils::P4 mom1, HEPUtils::P4 mom2)
static bool oppositeSign(const HEPUtils::Particle *a, const HEPUtils::Particle *b)
def combine(region_file, pip_file)
static std::vector< const HEPUtils::Particle * > leptonJetOverlapRemoval(std::vector< const HEPUtils::Particle *> leptons, std::vector< const HEPUtils::Jet *> jets, double dR)
DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb)
TODO: see if we can use this one:
Class for ColliderBit analyses.
Functions that do super fast ATLAS detector simulation based on four-vector smearing.