file analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp
[No description available]
Namespaces
| Name |
|---|
| Gambit TODO: see if we can use this one: |
| Gambit::ColliderBit |
Classes
| Name | |
|---|---|
| class | Gambit::ColliderBit::Analysis_ATLAS_13TeV_3b_36invfb |
| class | Gambit::ColliderBit::Analysis_ATLAS_13TeV_3b_discoverySR_36invfb |
Source code
///
/// \author Are Raklev
/// \date 2018 June
///
/// Based on the search presented in 1806.04030.
/// Only the high mass analysis is implemented here.
/// This analysis has overlapping exclusion and discovery signal regions,
/// the discovery regions are separated into a derived class.
/// *********************************************
#include <vector>
#include <cmath>
#include <memory>
#include <iomanip>
#include "gambit/ColliderBit/analyses/Analysis.hpp"
#include "gambit/ColliderBit/ATLASEfficiencies.hpp"
using namespace std;
// TODO: See if adding muons to jets gives some improvement.
namespace Gambit {
namespace ColliderBit {
class Analysis_ATLAS_13TeV_3b_36invfb : public Analysis {
protected:
// Signal region map
std::map<string, EventCounter> _counters = {
// Exclusion regions, disjoint
{"SR-3b-meff1-A", EventCounter("SR-3b-meff1-A")},
{"SR-3b-meff2-A", EventCounter("SR-3b-meff2-A")},
{"SR-3b-meff3-A", EventCounter("SR-3b-meff3-A")},
{"SR-4b-meff1-A", EventCounter("SR-4b-meff1-A")},
{"SR-4b-meff1-B", EventCounter("SR-4b-meff1-B")},
{"SR-4b-meff2-A", EventCounter("SR-4b-meff2-A")},
{"SR-4b-meff2-B", EventCounter("SR-4b-meff2-B")},
// Discovery regions, SR-4b-meff1-A and SR-4b-meff2-A are subsets
{"SR-4b-meff1-A-disc", EventCounter("SR-4b-meff1-A-disc")},
};
private:
// Cut-flows
size_t NCUTS;
vector<int> cutFlowVector;
vector<string> cutFlowVector_str;
vector<double> cutFlowVectorATLAS;
public:
// Required detector sim
static constexpr const char* detector = "ATLAS";
static bool sortByPT(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
Analysis_ATLAS_13TeV_3b_36invfb() {
set_analysis_name("ATLAS_13TeV_3b_36invfb");
set_luminosity(36.1);
NCUTS=14;
for(size_t i=0;i<NCUTS;i++){
cutFlowVector.push_back(0);
cutFlowVectorATLAS.push_back(0);
cutFlowVector_str.push_back("");
}
}
// The following section copied from Analysis_ATLAS_1LEPStop_20invfb.cpp
void JetLeptonOverlapRemoval(vector<const HEPUtils::Jet*> &jetvec, vector<const HEPUtils::Particle*> &lepvec, double DeltaRMax) {
//Routine to do jet-lepton check
//Discards jets if they are within DeltaRMax of a lepton
vector<const HEPUtils::Jet*> Survivors;
for(unsigned int itjet = 0; itjet < jetvec.size(); itjet++) {
bool overlap = false;
HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
double dR;
dR=jetmom.deltaR_eta(lepmom);
if(fabs(dR) <= DeltaRMax) overlap=true;
}
if(overlap) continue;
Survivors.push_back(jetvec.at(itjet));
}
jetvec=Survivors;
return;
}
void LeptonJetOverlapRemoval(vector<const HEPUtils::Particle*> &lepvec, vector<const HEPUtils::Jet*> &jetvec) {
//Routine to do lepton-jet check
//Discards leptons if they are within dR of a jet as defined in analysis paper
vector<const HEPUtils::Particle*> Survivors;
for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
bool overlap = false;
HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
for(unsigned int itjet= 0; itjet < jetvec.size(); itjet++) {
HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
double dR;
double DeltaRMax = std::min(0.4, 0.04 + 10 / lepmom.pT());
dR=jetmom.deltaR_eta(lepmom);
if(fabs(dR) <= DeltaRMax) overlap=true;
}
if(overlap) continue;
Survivors.push_back(lepvec.at(itlep));
}
lepvec=Survivors;
return;
}
// Calculate transverse mass
double mTrans(HEPUtils::P4 pmiss, HEPUtils::P4 jet) {
double mT = sqrt( pow(pmiss.pT()+jet.pT(),2) - pow(pmiss.px()+jet.px(),2) - pow(pmiss.py()+jet.py(),2) );
//cout << "pTmiss " << pmiss.pT() << " jetpT " << jet.pT() << endl;
//cout << "pxmiss " << pmiss.px() << "pxjet " << jet.px() << " pymiss " << pmiss.py() << " pyjet " << jet.py() << endl;
return mT;
}
void run(const HEPUtils::Event* event) {
// Get the missing energy in the event
double met = event->met();
HEPUtils::P4 metVec = event->missingmom();
// Now define vectors of baseline objects, including:
// - retrieval of electron, muon and jets from the event
// - application of basic pT and eta cuts
// Electrons
vector<const HEPUtils::Particle*> electrons;
for (const HEPUtils::Particle* electron : event->electrons()) {
if (electron->pT() > 5.
&& fabs(electron->eta()) < 2.47)
electrons.push_back(electron);
}
// Apply electron efficiency
ATLAS::applyElectronEff(electrons);
// Muons
vector<const HEPUtils::Particle*> muons;
for (const HEPUtils::Particle* muon : event->muons()) {
if (muon->pT() > 5.
&& fabs(muon->eta()) < 2.5)
muons.push_back(muon);
}
// Apply muon efficiency
ATLAS::applyMuonEff(muons);
vector<const HEPUtils::Jet*> candJets;
for (const HEPUtils::Jet* jet : event->jets()) {
if (jet->pT() > 20. && fabs(jet->eta()) < 2.8)
candJets.push_back(jet);
}
// Jets
vector<const HEPUtils::Jet*> bJets;
vector<const HEPUtils::Jet*> nonbJets;
// Find b-jets
double btag = 0.77; double cmisstag = 1/6.; double misstag = 1./134.;
for (const HEPUtils::Jet* jet : candJets) {
// Tag
if( jet->btag() && random_bool(btag) ) bJets.push_back(jet);
// Misstag c-jet
else if( jet->ctag() && random_bool(cmisstag) ) bJets.push_back(jet);
// Misstag light jet
else if( random_bool(misstag) ) bJets.push_back(jet);
// Non b-jet
else nonbJets.push_back(jet);
}
// Overlap removal
JetLeptonOverlapRemoval(nonbJets,electrons,0.2);
LeptonJetOverlapRemoval(electrons,nonbJets);
JetLeptonOverlapRemoval(nonbJets,muons,0.2);
LeptonJetOverlapRemoval(muons,nonbJets);
// Find veto leptons with pT > 20 GeV
vector<const HEPUtils::Particle*> vetoElectrons;
for (const HEPUtils::Particle* electron : electrons) {
if (electron->pT() > 20.) vetoElectrons.push_back(electron);
}
vector<const HEPUtils::Particle*> vetoMuons;
for (const HEPUtils::Particle* muon : muons) {
if (muon->pT() > 20.) vetoMuons.push_back(muon);
}
// Restrict jets to pT > 25 GeV after overlap removal
vector<const HEPUtils::Jet*> bJets_survivors;
for (const HEPUtils::Jet* jet : bJets) {
if(jet->pT() > 25.) bJets_survivors.push_back(jet);
}
vector<const HEPUtils::Jet*> nonbJets_survivors;
for (const HEPUtils::Jet* jet : nonbJets) {
if(jet->pT() > 25.) nonbJets_survivors.push_back(jet);
}
vector<const HEPUtils::Jet*> jet_survivors;
jet_survivors = nonbJets_survivors;
for (const HEPUtils::Jet* jet : bJets) {
jet_survivors.push_back(jet);
}
std::sort(jet_survivors.begin(), jet_survivors.end(), sortByPT);
// Number of objects
size_t nbJets = bJets_survivors.size();
size_t nnonbJets = nonbJets_survivors.size();
size_t nJets = nbJets + nnonbJets;
//size_t nJets = jet_survivors.size();
size_t nMuons=vetoMuons.size();
size_t nElectrons=vetoElectrons.size();
size_t nLeptons = nElectrons+nMuons;
// Loop over jets to find angle wrt to missing momentum
double phi4min = 7;
for(int i = 0; i < min(4,(int)nJets); i++){
double phi = jet_survivors.at(i)->mom().deltaPhi(metVec);
if(phi < phi4min) phi4min = phi;
}
// Collect the four signal jets.
vector<const HEPUtils::Jet*> signalJets;
for(const HEPUtils::Jet* jet : bJets_survivors){
if(signalJets.size() < 4) signalJets.push_back(jet);
}
for(const HEPUtils::Jet* jet : nonbJets_survivors){
if(signalJets.size() < 4) signalJets.push_back(jet);
}
// Effective mass (using the four jets used in Higgses)
double meff = met;
for(const HEPUtils::Jet* jet : signalJets){
meff += jet->pT();
}
// Find Higgs candidates
double mlead = 0; double msubl = 0;
double m1 = 0; double m2 = 0;
double Rbbmax = 10;
if(signalJets.size() == 4){
double R11 = signalJets.at(0)->mom().deltaR_eta(signalJets.at(1)->mom());
double R12 = signalJets.at(2)->mom().deltaR_eta(signalJets.at(3)->mom());
double DR1 = max(R11,R12);
//cout << DR1 << " " << R11 << " " << R12 << endl;
double R21 = signalJets.at(0)->mom().deltaR_eta(signalJets.at(2)->mom());
double R22 = signalJets.at(1)->mom().deltaR_eta(signalJets.at(3)->mom());
double DR2 = max(R21,R22);
//cout << DR2 << " " << R21 << " " << R22 << endl;
double R31 = signalJets.at(0)->mom().deltaR_eta(signalJets.at(3)->mom());
double R32 = signalJets.at(1)->mom().deltaR_eta(signalJets.at(2)->mom());
double DR3 = max(R31,R32);
//cout << DR3 << " " << R31 << " " << R32 << endl;
//cout << endl;
if( DR1 < DR2 && DR1 < DR3 ){
m1 = (signalJets.at(0)->mom()+signalJets.at(1)->mom()).m();
m2 = (signalJets.at(2)->mom()+signalJets.at(3)->mom()).m();
Rbbmax = DR1;
}
else if( DR2 < DR1 && DR2 < DR3 ){
m1 = (signalJets.at(0)->mom()+signalJets.at(2)->mom()).m();
m2 = (signalJets.at(1)->mom()+signalJets.at(3)->mom()).m();
Rbbmax = DR2;
}
else{
m1 = (signalJets.at(0)->mom()+signalJets.at(3)->mom()).m();
m2 = (signalJets.at(1)->mom()+signalJets.at(2)->mom()).m();
Rbbmax = DR3;
}
mlead = max(m1,m2); msubl = min(m1,m2);
//cout << mlead << " " << msubl << endl;
}
// Transverse mass for leading b-jets
double mTmin = 10E6;
for(int i = 0; i < min(3,(int)nbJets); i++){
double mT = mTrans(metVec,bJets_survivors.at(i)->mom());
if(mT < mTmin) mTmin = mT;
}
//cout << "mTmin " << mTmin << endl;
// Increment cutFlowVector elements
// Cut flow strings
// cutFlowVector_str[0] = "No cuts ";
// cutFlowVector_str[1] = "Trigger, $E_T^{miss} > 200$ GeV";
// cutFlowVector_str[2] = "$\\Delta\\phi_{min}^{4j} > 0.4$";
// cutFlowVector_str[3] = "$N_{lep} = 0$";
// cutFlowVector_str[4] = "$N_{jet} \\ge 4$, $N_{jet} \\le 5$";
// cutFlowVector_str[5] = "$110 < m(h_1)< 150$ GeV";
// cutFlowVector_str[6] = "$90 < m(h_2)< 140$ GeV$";
// cutFlowVector_str[7] = "$m_{T,min}^{b-jets}> 130$ GeV";
// cutFlowVector_str[8] = "$m_{eff} > 1100$ GeV";
// cutFlowVector_str[9] = "$N_{b-jets} \\ge 3$";
// cutFlowVector_str[10] = "$0.4 \\le \\Delta R_{max}^{bb} \\le 1.4$";
// cutFlowVector_str[11] = "m_{eff} > 600 GeV";
// cutFlowVector_str[12] = "$N_{b-jet} \\ge 4$";
// cutFlowVector_str[13] = "$0.4 \\le \\Delta R_{max}^{bb} \\le 1.4$";
// Cut flow from paper
// Higgsino 300 GeV
// cutFlowVectorATLAS[0] = 10276.0;
// cutFlowVectorATLAS[1] = 1959.1;
// cutFlowVectorATLAS[2] = 1533.0;
// cutFlowVectorATLAS[3] = 1319.3;
// cutFlowVectorATLAS[4] = 664.9;
// cutFlowVectorATLAS[5] = 249.3;
// cutFlowVectorATLAS[6] = 123.0;
// cutFlowVectorATLAS[7] = 74.3;
// cutFlowVectorATLAS[8] = 4.0;
// cutFlowVectorATLAS[9] = 1.5;
// cutFlowVectorATLAS[10] = 1.4;
// cutFlowVectorATLAS[11] = 90.2;
// cutFlowVectorATLAS[12] = 15.6;
// cutFlowVectorATLAS[13] = 6.8;
// Higgsino 500 GeV
// cutFlowVectorATLAS[0] = 1220.7;
// cutFlowVectorATLAS[1] = 739.0;
// cutFlowVectorATLAS[2] = 647.1;
// cutFlowVectorATLAS[3] = 548.2;
// cutFlowVectorATLAS[4] = 291.9;
// cutFlowVectorATLAS[5] = 133.5;
// cutFlowVectorATLAS[6] = 78.0;
// cutFlowVectorATLAS[7] = 64.1;
// cutFlowVectorATLAS[8] = 12.0;
// cutFlowVectorATLAS[9] = 5.7;
// cutFlowVectorATLAS[10] = 4.8;
// cutFlowVectorATLAS[11] = 74.3;
// cutFlowVectorATLAS[12] = 15.0;
// cutFlowVectorATLAS[13] = 9.7;
// Higgsino 800 GeV
// cutFlowVectorATLAS[0] = 124.9;
// cutFlowVectorATLAS[1] = 101.9;
// cutFlowVectorATLAS[2] = 89.5;
// cutFlowVectorATLAS[3] = 73.7;
// cutFlowVectorATLAS[4] = 39.4;
// cutFlowVectorATLAS[5] = 19.0;
// cutFlowVectorATLAS[6] = 13.4;
// cutFlowVectorATLAS[7] = 11.8;
// cutFlowVectorATLAS[8] = 8.1;
// cutFlowVectorATLAS[9] = 3.8;
// cutFlowVectorATLAS[10] = 3.6;
// cutFlowVectorATLAS[11] = 13.3;
// cutFlowVectorATLAS[12] = 2.3;
// cutFlowVectorATLAS[13] = 2.0;
// Apply cutflow
// for(size_t j=0;j<NCUTS;j++){
// if(
// (j==0) ||
//
// (j==1 && met > 200.) ||
//
// (j==2 && met > 200 && phi4min > 0.4) ||
//
// (j==3 && met > 200 && phi4min > 0.4 && nLeptons == 0) ||
//
// (j==4 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5)) ||
//
// (j==5 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150.) ||
//
// (j==6 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140.) ||
//
// (j==7 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && mTmin > 130.) ||
//
// (j==8 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && mTmin > 130. && meff > 1100.) ||
//
// (j==9 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && mTmin > 130. && meff > 1100. && nbJets >= 3) ||
//
// (j==10 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && mTmin > 130. && meff > 1100. && nbJets >= 3 && Rbbmax > 0.4 && Rbbmax < 1.4) ||
//
// (j==11 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && meff > 600.) ||
//
// (j==12 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && meff > 600. && nbJets >= 4) ||
//
// (j==13 && met > 200 && phi4min > 0.4 && nLeptons == 0 && (nJets == 4 || nJets == 5) && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && meff > 600. && nbJets >= 4 && Rbbmax > 0.4 && Rbbmax < 1.4)
//
// ) cutFlowVector[j]++;
// }
// Now increment signal region variables
// First exclusion regions
if(nbJets == 3 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && mTmin > 150. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff > 600. && meff < 850.) _counters.at("SR-3b-meff1-A").add_event(event);
if(nbJets == 3 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && mTmin > 150. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff > 850. && meff < 1100.) _counters.at("SR-3b-meff2-A").add_event(event);
if(nbJets >= 3 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && mTmin > 130. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff > 1100.) _counters.at("SR-3b-meff3-A").add_event(event);
if(nbJets >= 4 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && meff > 600. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff < 850.) _counters.at("SR-4b-meff1-A").add_event(event);
if(nbJets >= 4 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && meff > 600. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 1.4 && Rbbmax < 2.4 && meff < 850.) _counters.at("SR-4b-meff1-B").add_event(event);
if(nbJets >= 4 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 6 && meff > 850. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff < 1100.) _counters.at("SR-4b-meff2-A").add_event(event);
if(nbJets >= 4 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 6 && meff > 850. && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 1.4 && Rbbmax < 2.4 && meff < 1100.) _counters.at("SR-4b-meff2-B").add_event(event);
// Discovery regions
if(nbJets >= 4 && met > 200 && nLeptons == 0 && phi4min > 0.4 && nJets >= 4 && nJets <= 5 && mlead > 110. && mlead < 150. && msubl > 90. && msubl < 140. && Rbbmax > 0.4 && Rbbmax < 1.4 && meff > 600.) _counters.at("SR-4b-meff1-A-disc").add_event(event);
return;
} // End of analyze
/// Combine the variables of another copy of this analysis (typically on another thread) into this one.
void combine(const Analysis* other)
{
const Analysis_ATLAS_13TeV_3b_36invfb* specificOther
= dynamic_cast<const Analysis_ATLAS_13TeV_3b_36invfb*>(other);
for (auto& pair : _counters) { pair.second += specificOther->_counters.at(pair.first); }
if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
for (size_t j=0; j<NCUTS; j++) {
cutFlowVector[j] += specificOther->cutFlowVector[j];
cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
}
}
virtual void collect_results() {
// // DEBUG
// double L = 36.1;
//// double xsec = 284.65; // 300 GeV
//// double xsec = 33.81; // 500 GeV
// double xsec = 3.460; // 800 GeV
// cout << "DEBUG:" << endl;
// for (size_t i=0; i<NCUTS; i++)
// {
// double ATLAS_abs = cutFlowVectorATLAS[i];
//
// double eff = (double)cutFlowVector[i] / (double)cutFlowVector[0];
//
// double GAMBIT_scaled = eff * xsec * L;
//
// double ratio = GAMBIT_scaled/ATLAS_abs;
// cout << "DEBUG 1: i: " << i << ": " << setprecision(4) << ATLAS_abs << "\t" << GAMBIT_scaled << "\t" << "\t" << ratio << "\t\t" << cutFlowVector_str[i] << endl;
// }
// cout << "DEBUG:" << endl;
// Now fill a results object with the results for each SR
// Only exclusion regions here
add_result(SignalRegionData(_counters.at("SR-3b-meff1-A"), 4., {2.5, 1.0}));
add_result(SignalRegionData(_counters.at("SR-3b-meff2-A"), 3., {2.0, 0.5}));
add_result(SignalRegionData(_counters.at("SR-3b-meff3-A"), 0., {0.8, 0.5}));
add_result(SignalRegionData(_counters.at("SR-4b-meff1-A"), 1., {0.43, 0.31}));
add_result(SignalRegionData(_counters.at("SR-4b-meff1-B"), 2., {2.6, 0.9}));
add_result(SignalRegionData(_counters.at("SR-4b-meff2-A"), 1., {0.43, 0.27}));
add_result(SignalRegionData(_counters.at("SR-4b-meff2-B"), 0., {1.3, 0.6}));
return;
}
void analysis_specific_reset() {
// Clear signal regions
for (auto& pair : _counters) { pair.second.reset(); }
// Clear cut flow vector
std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
}
};
DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_3b_36invfb)
//
// Class for collecting results for discovery regions as a derived class
//
class Analysis_ATLAS_13TeV_3b_discoverySR_36invfb : public Analysis_ATLAS_13TeV_3b_36invfb {
public:
Analysis_ATLAS_13TeV_3b_discoverySR_36invfb() {
set_analysis_name("ATLAS_13TeV_3b_discoverySR_36invfb");
}
virtual void collect_results() {
add_result(SignalRegionData(_counters.at("SR-4b-meff1-A-disc"), 2., {0.7, 0.5}));
}
};
// Factory fn
DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_3b_discoverySR_36invfb)
}
}
Updated on 2023-06-26 at 21:36:56 +0000