gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp
Go to the documentation of this file.
1 #include <iomanip>
2 
6 
7 using namespace std;
8 
9 // The CMS 1 lepton DM + top pair analysis (20fb^-1)
10 
11 // based on: https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsB2G14004
12 
13 // Code by Martin White, Guy Pitman
14 // Known issues:
15 // a) Impossible to test results against CMS due to the impossibility of reproducing their model information (even after contacting CMS). Note that the variables used have been debugged in other contexts however.
16 // b) Overlap removal is not applied (CMS do not use it, but we don't exactly use their particle flow technique either)
17 // c) Jets here need kT radius of 0.5 not 0.4
18 
19 namespace Gambit {
20  namespace ColliderBit {
21 
22 
23 
24  //Puts dphi in the range -pi to pi
25  double _Phi_mpi_pi(double x){
26  while (x >= M_PI) x -= 2*M_PI;
27  while (x < -M_PI) x += 2*M_PI;
28  return x;
29  }
30 
31 
32  class Analysis_CMS_8TeV_1LEPDMTOP_20invfb : public Analysis {
33  private:
34 
35  // Numbers passing cuts
36  double _numSR;
37 
38  vector<int> cutFlowVector;
39  vector<string> cutFlowVector_str;
40  int NCUTS; //=24;
41 
42  public:
43 
44  // Required detector sim
45  static constexpr const char* detector = "CMS";
46 
47  Analysis_CMS_8TeV_1LEPDMTOP_20invfb()
48  : _numSR(0),
49  NCUTS(6)
50  {
51  set_analysis_name("CMS_8TeV_1LEPDMTOP_20invfb");
52  set_luminosity(19.7);
53 
54  for (int i=0; i<NCUTS; i++) {
55  cutFlowVector.push_back(0);
56  cutFlowVector_str.push_back("");
57  }
58  }
59 
60  double SmallestdPhi(std::vector<const HEPUtils::Jet*> jets,double phi_met)
61  {
62  if (jets.size()<2) return(999);
63  double dphi1 = std::acos(std::cos(jets.at(0)->phi()-phi_met));
64  double dphi2 = std::acos(std::cos(jets.at(1)->phi()-phi_met));
65  // double dphi3 = 999;
66  //if (jets.size() > 2 && jets[2]->pT() > 40.)
67  // dphi3 = std::acos(std::cos(jets[2]->phi() - phi_met));
68  double min1 = std::min(dphi1, dphi2);
69 
70  return min1;
71 
72  }
73 
74  void run(const HEPUtils::Event* event) {
75 
76  // Missing energy
77  HEPUtils::P4 ptot = event->missingmom();
78  double met = event->met();
79 
80  // Baseline electrons
81  vector<const HEPUtils::Particle*> baselineElectrons;
82  for (const HEPUtils::Particle* electron : event->electrons()) {
83  if (electron->pT() > 30. && fabs(electron->eta()) < 2.5) {
84  baselineElectrons.push_back(electron);
85  }
86  }
87 
88  // Apply electron efficiency
89  CMS::applyElectronEff(baselineElectrons);
90 
91  // Baseline muons
92  vector<const HEPUtils::Particle*> baselineMuons;
93  for (const HEPUtils::Particle* muon : event->muons()) {
94  if (muon->pT() > 30. && fabs(muon->eta()) < 2.1) {
95  baselineMuons.push_back(muon);
96  }
97  }
98 
99  // Apply muon efficiency
100  CMS::applyMuonEff(baselineMuons);
101 
102  // All baseline leptons
103  vector<const HEPUtils::Particle*> baselineLeptons = baselineElectrons;
104  baselineLeptons.insert(baselineLeptons.end(), baselineMuons.begin(), baselineMuons.end() );
105 
106  vector<const HEPUtils::Jet*> baselineJets;
107  //vector<LorentzVector> jets;
108  vector<HEPUtils::P4> jets;
109  vector<const HEPUtils::Jet*> bJets;
110  vector<bool> btag;
111 
112  const std::vector<double> a = {0,10.};
113  const std::vector<double> b = {0,10000.};
114  const std::vector<double> c = {0.60};
115  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
116 
117  for (const HEPUtils::Jet* jet : event->jets()) {
118  if (jet->pT() > 30. && fabs(jet->eta()) < 4.0) {
119  baselineJets.push_back(jet);
120  //LorentzVector j1 (jet->mom().px(),jet->mom().py(),jet->mom().pz(),jet->mom().E()) ;
121  //jets.push_back(j1);
122  jets.push_back(jet->mom());
123  bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
124  bool isB=false;
125 
126  if(jet->btag() && hasTag && fabs(jet->eta()) < 2.4 && jet->pT() > 30.) {
127  isB=true;
128  bJets.push_back(jet);
129  }
130  btag.push_back(isB);
131  }
132  }
133 
134  // Calculate common variables and cuts first
135  //applyTightIDElectronSelection(signalElectrons);
136 
137  //int nElectrons = signalElectrons.size();
138  //int nMuons = signalMuons.size();
139  int nJets = baselineJets.size();
140  int nLeptons = baselineLeptons.size();
141  int nBJets = bJets.size();
142 
143  //Preselection cuts
144  bool passPresel=false;
145  if(nLeptons==1 &&
146  nJets>=3 &&
147  nBJets>=1 &&
148  met > 160.)passPresel=true;
149 
150  //Calculate mT
151  HEPUtils::P4 lepVec;
152  double mT=0;
153  if(nLeptons==1){
154  lepVec=baselineLeptons[0]->mom();
155  mT=sqrt(2.*lepVec.pT()*met*(1. - cos(_Phi_mpi_pi(lepVec.phi()-ptot.phi()))));
156  }
157 
158  //Calculate MT2W
159  double MT2W=0;
160  // double MT2W_HU=0;
161  if (nJets > 1 && nLeptons==1) {
162  HEPUtils::P4 lepVec;
163  lepVec=baselineLeptons[0]->mom();
164  //LorentzVector lep (lepVec.px(),lepVec.py(),lepVec.pz(),lepVec.E());
165  float phi=float (ptot.phi());
166  //MT2W=calculateMT2w(jets, btag, lep, met, phi);
167  MT2W=calculateMT2wHepUtils(jets,btag,lepVec,met,phi);
168  }
169 
170  //Calculate dPhi variable
171  float phi=float (ptot.phi());
172  double dPhiMin12=SmallestdPhi(baselineJets,phi);
173 
174  //Cuts
175  //MET > 320
176  //MT > 160
177  //MT2W > 300
178  //dPhiMin12 > 1.2
179 
180  cutFlowVector_str[0] = "No cuts ";
181  cutFlowVector_str[1] = "Presel ";
182  cutFlowVector_str[2] = "MET > 320 GeV ";
183  cutFlowVector_str[3] = "MT > 160 GeV ";
184  cutFlowVector_str[4] = "MT2W > 300 GeV ";
185  cutFlowVector_str[5] = "dPhiMin12 > 1.2 ";
186 
187  for(int j=0;j<NCUTS;j++){
188  if(
189  (j==0) ||
190 
191  (j==1 && passPresel) ||
192 
193  (j==2 && passPresel && met > 320.) ||
194 
195  (j==3 && passPresel && met > 320. && mT > 160.) ||
196 
197  (j==4 && passPresel && met > 320. && mT > 160. && MT2W > 300.) ||
198 
199  (j==5 && passPresel && met > 320. && mT > 160. && MT2W > 300. && dPhiMin12 > 1.2))
200 
201  cutFlowVector[j]++;
202  }
203 
204  //We're now ready to apply the cuts for each signal region
205  //_numSR1, _numSR2, _numSR3;
206 
207  if(passPresel && met > 320. && mT > 160. && MT2W > 300. && dPhiMin12 > 1.2) _numSR += event->weight();
208 
209  return;
210  }
211 
213  void combine(const Analysis* other)
214  {
215  const Analysis_CMS_8TeV_1LEPDMTOP_20invfb* specificOther
216  = dynamic_cast<const Analysis_CMS_8TeV_1LEPDMTOP_20invfb*>(other);
217  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
218  for (int j=0; j<NCUTS; j++) {
219  cutFlowVector[j] += specificOther->cutFlowVector[j];
220  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
221  }
222  _numSR += specificOther->_numSR;
223  }
224 
225 
226  double loglikelihood() {
228  return 0;
229  }
230 
231  void collect_results() {
232 
233  // add_result(SignalRegionData("SR label", n_obs, {n_sig_MC, n_sig_MC_sys}, {n_bkg, n_bkg_err}));
234  add_result(SignalRegionData("SR", 18., {_numSR, 0.}, { 16.4, 3.48}));
235 
236  return;
237  }
238 
239 
240  protected:
241  void analysis_specific_reset() {
242  _numSR = 0;
243  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
244  }
245 
246  };
247 
248 
249  DEFINE_ANALYSIS_FACTORY(CMS_8TeV_1LEPDMTOP_20invfb)
250 
251 
252  }
253 }
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.
Definition: Utils.hpp:207
#define DEFINE_ANALYSIS_FACTORY(ANAME)
For analysis factory function definition.
Definition: Analysis.hpp:124
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 *)
STL namespace.
START_MODEL b
Definition: demo.hpp:270
Functions that do super fast CMS detector simulation based on four-vector smearing.
double calculateMT2wHepUtils(vector< HEPUtils::P4 > &jets, vector< bool > &btag, HEPUtils::P4 &lep, float met, float metphi)
Definition: mt2w.cc:3
A class for collider analyses within ColliderBit.
Definition: Analysis.hpp:41
A simple container for the result of one signal region from one analysis.
void applyMuonEff(std::vector< const HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon efficiency.
def combine(region_file, pip_file)
Definition: colouring.py:169
void applyElectronEff(std::vector< const HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Class for ColliderBit analyses.