LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
MaxProbBreitWignerCombJigsaw.cc
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2016, Christopher Rogan
14 // This file is part of RestFrames.
15 //
16 // RestFrames is free software; you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation; either version 2 of the License, or
19 // (at your option) any later version.
20 //
21 // RestFrames is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with RestFrames. If not, see <http://www.gnu.org/licenses/>.
29 
30 #include "RestFrames/RestFrame.hh"
32 
33 namespace RestFrames {
34 
35  MaxProbBreitWignerCombJigsaw::MaxProbBreitWignerCombJigsaw(const std::string& sname,
36  const std::string& stitle,
37  int N_comb, int N_object)
38  : CombinatoricJigsaw(sname, stitle, N_comb, N_object),
39  m_Ncomb(N_comb), m_Nobj(N_object)
40  {
41  for(int i = 0; i < m_Nobj; i++){
42  m_Mass.push_back(0.);
43  m_Width.push_back(-1.);
44  }
45  }
46 
47  MaxProbBreitWignerCombJigsaw::~MaxProbBreitWignerCombJigsaw() {}
48 
49  void MaxProbBreitWignerCombJigsaw::SetPoleMass(double mass, int i){
50  if(i < 0 || i >= m_Nobj)
51  return;
52 
53  m_Mass[i] = std::max(mass, 0.);
54  }
55 
56  void MaxProbBreitWignerCombJigsaw::SetWidth(double width, int i){
57  if(i < 0 || i >= m_Nobj)
58  return;
59  if(width <= 0.)
60  m_Width[i] = -1.;
61  else
62  m_Width[i] = width;
63  }
64 
65  bool MaxProbBreitWignerCombJigsaw::EvaluateMetric(double& metric) const {
66  std::vector<TLorentzVector> P;
67  for(int i = 0; i < m_Nobj; i++)
68  P.push_back(GetDependancyStates(i).GetFourVector());
69 
70  double prob = 1.;
71  TLorentzVector SUM(0.,0.,0.,0.);
72  for(int i = 0; i < m_Nobj; i++)
73  SUM += P[i];
74  double M = SUM.M();
75  for(int i = 0; i < m_Nobj-1; i++){
76  SUM -= P[i];
77  prob *= GetP((P[i]+SUM).M(), P[i].M(), SUM.M())/M;
78  }
79 
80  double den;
81  for(int i = 0; i < m_Nobj; i++){
82  if(m_Width[i] > 0.){
83  den = P[i].M2()-m_Mass[i]*m_Mass[i];
84  den *= den;
85  den += m_Mass[i]*m_Mass[i]*m_Width[i]*m_Width[i];
86  if(den > 0.)
87  prob /= den;
88  }
89  }
90 
91  if(prob <= 0.)
92  return false;
93 
94  return 1./prob;
95  return true;
96  }
97 
98 }
99