ATLAS Offline Software
TtresNeutrinoBuilder.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3  */
4 
5 #include <TMinuit.h>
6 #include "TopEventReconstructionTools/TtresNeutrinoBuilder.h"
7 
8 //_________________________________________________________________________________________________
9 void delta2_fcn(Int_t&, Double_t*, Double_t&, Double_t*, Int_t);
10 
11 //_________________________________________________________________________________________________
12 TtresNeutrinoBuilder::TtresNeutrinoBuilder(std::string units) {
13  m_debug = 0;
14  if (units == "MeV") m_Units = 1000.;
15  else if (units == "GeV") m_Units = 1.;
16  else std::cerr << "WARNING in NeutrinoBuilder: units different from GeV or MeV" << std::endl;
17 }
18 
19 //_________________________________________________________________________________________________
20 TtresNeutrinoBuilder::~TtresNeutrinoBuilder() {
21  if (m_debug > 0) std::cout << "in destructor " << std::endl;
22 }
23 
24 //_________________________________________________________________________________________________
25 TtresNeutrinoBuilder::TtresNeutrinoBuilder(const TtresNeutrinoBuilder& other) {
26  m_debug = other.m_debug;
27 }
28 
29 //_________________________________________________________________________________________________
30 TtresNeutrinoBuilder& TtresNeutrinoBuilder::operator = (const TtresNeutrinoBuilder& other) {
31  if (this != &other) {
32  m_debug = other.m_debug;
33  }
34  return *this;
35 }
36 
37 //_________________________________________________________________________________________________
38 // JDM - fix compiler warnings
39 void delta2_fcn(Int_t& /*npar*/, Double_t* /*grad*/, Double_t& f, Double_t* par, Int_t /*iflag*/) {
40  Double_t delta2 = 0;
41  Double_t alpha = par[0];
42  Double_t r = par[1];
43  Double_t dphi = par[2];
44  Double_t l_pt = par[3];
45  Double_t l_m = par[4];
46  Double_t n_px = par[5];
47  Double_t n_py = par[6];
48 
49  r /= sqrt(l_pt * l_pt + l_m * l_m) - l_pt * cos(dphi + alpha);
50  TLorentzVector* neut = new TLorentzVector(n_px, n_py, 0., 0.);
51  neut->SetE(neut->P());
52 
53  TLorentzVector* neut_new =
54  new TLorentzVector(r * neut->P() * cos(neut->Phi() + alpha), r * neut->P() * sin(neut->Phi() + alpha), 0., 0.);
55  neut_new->SetE(neut_new->P());
56 
57  delta2 = pow((neut_new->Px() - neut->Px()), 2) + pow((neut_new->Py() - neut->Py()), 2);
58  r *= sqrt(l_pt * l_pt + l_m * l_m) - l_pt * cos(dphi + alpha);
59  delete neut;
60  delete neut_new;
61  f = delta2;
62 }
63 
64 //_________________________________________________________________________________________________
65 double TtresNeutrinoBuilder::fitAlpha(const TLorentzVector* L, const Double_t met, const Double_t metphi) {
66  // initialize
67  double m_mWpdg = 80.4 * m_Units;
68  Double_t pxNu = met * cos(metphi);
69  Double_t pyNu = met * sin(metphi);
70  Double_t ptNu = met;
71 
72  TMinuit* fit = new TMinuit(7);
73 
74  fit->SetFCN(delta2_fcn);
75  int ierr = 0;
76  double arglist[1] = {
77  -1
78  };
79  fit->mnexcm("SET PRIN", arglist, 1, ierr);
80 
81  // Initialise the parameters
82  std::string par_name[7] = {
83  "alpha", "r", "dphi", "l_pt", "l_m", "n_px", "n_py"
84  };
85  Double_t par_ival[7] = {
86  0., (m_mWpdg * m_mWpdg - L->M() * L->M()) / (2 * ptNu), metphi - L->Phi(), L->Pt(), L->M(), pxNu, pyNu
87  };
88  Double_t par_step[7] = {
89  0.1, 0., 0., 0., 0., 0., 0.
90  };
91  Double_t par_min[7] = {
92  -3.15, 0., -3.15, 0., 0., -10000., -10000.
93  };
94  Double_t par_max[7] = {
95  3.15, 1., 3.15, 10000., 80., 10000., 10000.
96  };
97  for (Int_t i = 0; i < 7; i++) {
98  fit->DefineParameter(i, par_name[i].c_str(), par_ival[i], par_step[i], par_min[i], par_max[i]);
99  if (i != 0) {
100  fit->FixParameter(i);
101  }
102  }
103  fit->SetPrintLevel(-1);
104  fit->Migrad();
105  Double_t a, e_a;
106  Int_t ret = fit->GetParameter(0, a, e_a);
107  delete fit;
108  if (ret > 0) {
109  return a;
110  } else {
111  std::cout << "Error in fit of alpha for met correction" << std::endl;
112  return 0.;
113  }
114 }
115 
116 // In case of negative discriminant, decrease the MET
117 //_________________________________________________________________________________________________
118 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_Scaling(const TLorentzVector* L,
119  const TLorentzVector* MET) {
120  return this->candidatesFromWMass_Scaling(L, MET->Px(), MET->Py());
121 }
122 
123 // In case of negative discriminant, decrease the MET
124 //_________________________________________________________________________________________________
125 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_Scaling(const TLorentzVector* L,
126  const double met_etx,
127  const double met_ety) {
128  if (m_debug > 0) std::cout << "entering NeutrinoBuilder::candidatesFromWMass_Scaling()" << std::endl;
129  std::vector<TLorentzVector*> NC;
130  this->candidatesFromWMass_Scaling(L, met_etx, met_ety, NC);
131  return NC;
132 }
133 
134 // In case of negative discriminant, decrease the MET
135 //_________________________________________________________________________________________________
136 bool TtresNeutrinoBuilder::candidatesFromWMass_Scaling(const TLorentzVector* L, const double met_etx,
137  const double met_ety, std::vector<TLorentzVector*>& NC) {
138  if (m_debug > 0) std::cout << "entering NeutrinoBuilder::candidatesFromWMass_Scaling()" << std::endl;
139  // initialize
140  double delta;
141  double lambda;
142  double pxNu;
143  double pyNu;
144  double ptNu = 21;
145  double alpha;
146  double beta;
147  double gamma;
148  double m_mWpdg = 80.4 * m_Units;
149  int count = 0;
150 
151  pxNu = met_etx;
152  pyNu = met_ety;
153  do {
154  // solve the quadratic equation
155  ++count;
156  ptNu = sqrt(pxNu * pxNu + pyNu * pyNu);
157  alpha = pow(m_mWpdg, 2) + pow((pxNu + L->Px()), 2) + pow((pyNu + L->Py()), 2) - pow(L->E(), 2);
158 
159  beta = 0.5 * (alpha - pow(ptNu, 2) + pow(L->Pz(), 2));
160  gamma = -(beta * beta - (pow(L->E(), 2) * pow(ptNu, 2))) / (pow(L->E(), 2) - pow(L->Pz(), 2));
161  lambda = 2 * beta * L->Pz() / (pow(L->E(), 2) - pow(L->Pz(), 2));
162  delta = pow(lambda, 2) - 4 * gamma;
163  if (m_debug >
164  1) std::cout << "count = " << count << "\tptNu = " << ptNu << "\tdelta = " << delta << "\tm_mWpdg = " <<
165  m_mWpdg <<
166  std::endl;
167  // if no real solution, reduce Pt by 0.1 GeV and recompute pxNu & pyNu
168  // in consequence
169  if (delta < 0) {
170  if (ptNu > 100.) {
171  double Ptmiss = ptNu;
172  ptNu -= 100.;
173  pxNu *= ptNu / Ptmiss;
174  pyNu *= ptNu / Ptmiss;
175  }
176  } // end of case delta is negatively defined
177  } while ((delta < 0) && (count < 10000) && (ptNu > 20 * 1000));
178  if (delta < 0) {
179  if (m_debug > 0) std::cout << "no solution delta still<0" << std::endl;
180  return false;
181  }
182  delta = sqrt(delta);
183  // instantiate Neutrino
184  double pz = (lambda - delta) / 2.0;
185  double e = sqrt(pxNu * pxNu + pyNu * pyNu + pz * pz);
186  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pz, e);
187  pz = (lambda + delta) / 2.0;
188  e = sqrt(pxNu * pxNu + pyNu * pyNu + pz * pz);
189  TLorentzVector* nu2 = new TLorentzVector(pxNu, pyNu, pz, e);
190  NC.push_back(nu1);
191  NC.push_back(nu2);
192  if (m_debug > 0) std::cout << "quitting NeutrinoBuilder::candidatesFromWMass_Scaling()" << std::endl;
193  return true;
194 }
195 
196 // In case of negative discriminant, decrese the MET
197 //_________________________________________________________________________________________________
198 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_Rotation(const TLorentzVector* L,
199  const TLorentzVector* MET,
200  const bool useSmallestPz) {
201  return this->candidatesFromWMass_Rotation(L, MET->Pt(), MET->Phi(), useSmallestPz);
202 }
203 
204 // In case of negative discriminant, decrese the MET
205 //_________________________________________________________________________________________________
206 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_Rotation(const TLorentzVector* L,
207  const Double_t met,
208  const Double_t metphi,
209  const bool useSmallestPz) {
210  if (m_debug > 0) std::cout << "entering candidatesFromWMassRotation()" << std::endl;
211 
212  // initialize
213  Double_t m_mWpdg = 80.4 * m_Units;
214  Double_t pxNu = met * cos(metphi);
215  Double_t pyNu = met * sin(metphi);
216  Double_t pzNu = -1000000;
217  Double_t ptNu = met;
218  Double_t eNu;
219 
220  std::vector<TLorentzVector*> NC;
221 
222  Double_t c1 = m_mWpdg * m_mWpdg - L->M() * L->M() + 2 * (L->Px() * pxNu + L->Py() * pyNu);
223  Double_t b1 = 2 * L->Pz();
224 
225  Double_t A = 4 * pow(L->E(), 2) - b1 * b1;
226  Double_t B = -2 * c1 * b1;
227  Double_t C = 4 * pow(L->E(), 2) * ptNu * ptNu - c1 * c1;
228  Double_t discr = B * B - 4 * A * C;
229  Double_t r = 1;
230 
231  Double_t sol1, sol2;
232  if (discr > 0) {
233  sol1 = (-B + sqrt(discr)) / (2 * A);
234  sol2 = (-B - sqrt(discr)) / (2 * A);
235  } else {
236  Double_t alpha = fitAlpha(L, met, metphi);
237  Double_t dphi = metphi - L->Phi();
238  r =
239  (pow(m_mWpdg,
240  2) - pow(L->M(), 2)) / (2 * ptNu * (sqrt(pow(L->Pt(), 2) + pow(L->M(), 2)) - L->Pt() * cos(dphi + alpha)));
241 
242  Double_t old_p = ptNu;
243  Double_t old_phi = metphi;
244  pxNu = r * old_p * cos(old_phi + alpha);
245  pyNu = r * old_p * sin(old_phi + alpha);
246  ptNu = sqrt(pxNu * pxNu + pyNu * pyNu);
247 
248  c1 = m_mWpdg * m_mWpdg - pow(L->M(), 2) + 2 * (L->Px() * pxNu + L->Py() * pyNu);
249  B = -2 * c1 * b1;
250  C = 4 * pow(L->E(), 2) * ptNu * ptNu - c1 * c1;
251  discr = B * B - 4 * A * C;
252 
253  sol1 = -B / (2 * A);
254  sol2 = -B / (2 * A);
255  }
256 
257  if (useSmallestPz) {
258  pzNu = (fabs(sol1) > fabs(sol2)) ? sol2 : sol1;
259 
260  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
261  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
262  NC.push_back(nu1);
263  } else {
264  pzNu = sol1;
265  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
266  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
267  pzNu = sol2;
268  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
269  TLorentzVector* nu2 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
270  NC.push_back(nu1);
271  NC.push_back(nu2);
272  }
273 
274  if (m_debug > 0) std::cout << "quitting NeutrinoBuilder::candidatesFromWMassRotation() : " << NC.size() << std::endl;
275  return NC;
276 }
277 
278 // In case of negative discriminant, use the real part
279 //_________________________________________________________________________________________________
280 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_RealPart(const TLorentzVector* L,
281  const TLorentzVector* MET,
282  const bool useSmallestPz) {
283  return this->candidatesFromWMass_RealPart(L, MET->Pt(), MET->Phi(), useSmallestPz);
284 }
285 
286 // In case of negative discriminant, use the real part
287 //_________________________________________________________________________________________________
288 std::vector<TLorentzVector*> TtresNeutrinoBuilder::candidatesFromWMass_RealPart(const TLorentzVector* L, Double_t met,
289  Double_t metphi,
290  const bool useSmallestPz) {
291  if (m_debug > 0) std::cout << "entering candidatesFromWMass_RealPart()" << std::endl;
292 
293  // initialize
294  Double_t m_mWpdg = 80.4 * m_Units;
295  Double_t pxNu = met * cos(metphi);
296  Double_t pyNu = met * sin(metphi);
297  Double_t pzNu = -1000000;
298  Double_t ptNu = met;
299  Double_t eNu;
300 
301  std::vector<TLorentzVector*> NC;
302 
303  Double_t c1 = m_mWpdg * m_mWpdg - L->M() * L->M() + 2 * (L->Px() * pxNu + L->Py() * pyNu);
304  Double_t b1 = 2 * L->Pz();
305 
306  Double_t A = 4 * pow(L->E(), 2) - b1 * b1;
307  Double_t B = -2 * c1 * b1;
308  Double_t C = 4 * pow(L->E(), 2) * ptNu * ptNu - c1 * c1;
309  Double_t discr = B * B - 4 * A * C;
310 
311  Double_t sol1, sol2;
312  if (discr > 0) {
313  sol1 = (-B + sqrt(discr)) / (2 * A);
314  sol2 = (-B - sqrt(discr)) / (2 * A);
315  }
316  //if discr<0
317  else {
318  pzNu = -B / (2 * A);
319 
320  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
321  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
322  NC.push_back(nu1);
323  return NC;
324  }
325 
326  if (useSmallestPz) {
327  pzNu = (fabs(sol1) > fabs(sol2)) ? sol2 : sol1;
328 
329  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
330  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
331  NC.push_back(nu1);
332  } else {
333  pzNu = sol1;
334  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
335  TLorentzVector* nu1 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
336  pzNu = sol2;
337  eNu = sqrt(pxNu * pxNu + pyNu * pyNu + pzNu * pzNu);
338  TLorentzVector* nu2 = new TLorentzVector(pxNu, pyNu, pzNu, eNu);
339  NC.push_back(nu1);
340  NC.push_back(nu2);
341  }
342 
343  if (m_debug > 0) std::cout << "quitting NeutrinoBuilder::candidatesFromWMass_RealPart() : " << NC.size() << std::endl;
344 
345  return NC;
346 }
347 
348 Double_t TtresNeutrinoBuilder::getDiscriminant(const TLorentzVector* L, Double_t met, Double_t metphi) {
349  // initialize
350  Double_t m_mWpdg = 80.4 * m_Units;
351  Double_t pxNu = met * cos(metphi);
352  Double_t pyNu = met * sin(metphi);
353  Double_t ptNu = met;
354 
355  Double_t c1 = m_mWpdg * m_mWpdg - L->M() * L->M() + 2 * (L->Px() * pxNu + L->Py() * pyNu);
356  Double_t b1 = 2 * L->Pz();
357 
358  Double_t A = 4 * pow(L->E(), 2) - b1 * b1;
359  Double_t B = -2 * c1 * b1;
360  Double_t C = 4 * pow(L->E(), 2) * ptNu * ptNu - c1 * c1;
361  Double_t discr = B * B - 4 * A * C;
362 
363  return discr;
364 }
beamspotman.r
def r
Definition: beamspotman.py:676
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
TtresNeutrinoBuilder::~TtresNeutrinoBuilder
virtual ~TtresNeutrinoBuilder()
Definition: TtresNeutrinoBuilder.cxx:20
extractSporadic.c1
c1
Definition: extractSporadic.py:134
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
DMTest::C
C_v1 C
Definition: C.h:26
TtresNeutrinoBuilder::m_Units
double m_Units
Definition: TtresNeutrinoBuilder.h:41
TRT_PAI_gasdata::NC
const int NC
Number of levels for Carbon.
Definition: TRT_PAI_gasdata.h:237
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
TtresNeutrinoBuilder
Definition: TtresNeutrinoBuilder.h:16
TtresNeutrinoBuilder::TtresNeutrinoBuilder
TtresNeutrinoBuilder(std::string units)
Definition: TtresNeutrinoBuilder.cxx:12
met
Definition: IMETSignificance.h:24
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:26
lumiFormat.i
int i
Definition: lumiFormat.py:92
ret
T ret(T t)
Definition: rootspy.cxx:260
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
TCS::MET
@ MET
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/L1TopoCommon/Types.h:16
TtresNeutrinoBuilder::getDiscriminant
Double_t getDiscriminant(const TLorentzVector *, const Double_t, const Double_t)
Definition: TtresNeutrinoBuilder.cxx:348
MuonValidation_CreateResolutionProfiles.fit
def fit(h, emin, emax)
Definition: MuonValidation_CreateResolutionProfiles.py:69
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
dqt_zlumi_alleff_HIST.B
B
Definition: dqt_zlumi_alleff_HIST.py:110
perfmonmt-refit.units
string units
Definition: perfmonmt-refit.py:77
delta2_fcn
void delta2_fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t)
Definition: TtresNeutrinoBuilder.cxx:39
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
InDetDD::other
@ other
Definition: InDetDD_Defs.h:16
TtresNeutrinoBuilder::operator=
TtresNeutrinoBuilder & operator=(const TtresNeutrinoBuilder &)
Definition: TtresNeutrinoBuilder.cxx:30
TtresNeutrinoBuilder::m_debug
int m_debug
Definition: TtresNeutrinoBuilder.h:40
TtresNeutrinoBuilder::candidatesFromWMass_Scaling
bool candidatesFromWMass_Scaling(const TLorentzVector *, double, Double_t, std::vector< TLorentzVector * > &)
TtresNeutrinoBuilder::candidatesFromWMass_Rotation
std::vector< TLorentzVector * > candidatesFromWMass_Rotation(const TLorentzVector *, const Double_t, const Double_t, const bool)
Definition: TtresNeutrinoBuilder.cxx:206
TtresNeutrinoBuilder::candidatesFromWMass_RealPart
std::vector< TLorentzVector * > candidatesFromWMass_RealPart(const TLorentzVector *, const Double_t, const Double_t, const bool)
Definition: TtresNeutrinoBuilder.cxx:288
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
TtresNeutrinoBuilder::fitAlpha
double fitAlpha(const TLorentzVector *, const Double_t, const Double_t)
Definition: TtresNeutrinoBuilder.cxx:65
TtresNeutrinoBuilder.h
Generated on Tue Jul 2 2024 21:32:06 for ATLAS Offline Software by   doxygen 1.8.18