LOGO

RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
ResonanceGenFrame.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 <stdlib.h>
32 
33 namespace RestFrames {
34 
36  // ResonanceGenFrame class
38 
39  ResonanceGenFrame::ResonanceGenFrame(const std::string& sname,
40  const std::string& stitle)
41  : DecayGenFrame(sname,stitle)
42  {
43  m_PoleMass = 0.;
44  m_Width = 0.;
45  }
46 
47  ResonanceGenFrame::ResonanceGenFrame() : DecayGenFrame()
48  {
49  m_PoleMass = 0.;
50  m_Width = 0.;
51  }
52 
53  ResonanceGenFrame::~ResonanceGenFrame() {}
54 
55  ResonanceGenFrame& ResonanceGenFrame::Empty(){
56  return ResonanceGenFrame::m_Empty;
57  }
58 
59  void ResonanceGenFrame::SetMass(double val){
60  SetMind(false);
61 
62  if(val < 0.){
63  m_Log << LogWarning;
64  m_Log << "Unable to set mass to negative value ";
65  m_Log << val << ". Setting to zero." << LogEnd;
66  m_PoleMass = 0.;
67  m_Mass = 0.;
68  } else {
69  m_PoleMass = val;
70  m_Mass = val;
71  }
72  }
73 
74  void ResonanceGenFrame::SetWidth(double val){
75  SetMind(false);
76 
77  if(val < 0.){
78  m_Log << LogWarning;
79  m_Log << "Unable to set width to negative value ";
80  m_Log << val << ". Setting to zero." << LogEnd;
81  m_Width = 0.;
82  SetVariableMassMCMC(false);
83  } else {
84  m_Width = val;
85  SetVariableMassMCMC(true);
86  }
87  }
88 
89  void ResonanceGenFrame::SetVariableMass(bool varymass) {
90  SetMind(false);
91 
92  if(varymass){
93  if(m_Width > 0.){
94  SetVariableMassMCMC(true);
95  }
96  else {
97  m_Log << LogWarning;
98  m_Log << "Unable to set variable mass. ";
99  m_Log << "Resonance width is set to zero. " << LogEnd;
100  }
101  } else {
102  m_Mass = m_PoleMass;
103  SetVariableMassMCMC(false);
104  }
105  }
106 
107 
108  double ResonanceGenFrame::GetPoleMass() const {
109  return m_PoleMass;
110  }
111 
112  double ResonanceGenFrame::GetWidth() const {
113  return m_Width;
114  }
115 
116  double ResonanceGenFrame::GetProbMCMC(double mass) const {
117  if(mass < 0)
118  mass = GetMass();
119 
120  double den = mass*mass-m_PoleMass*m_PoleMass;
121  den *= den;
122  den += m_PoleMass*m_PoleMass*m_Width*m_Width;
123 
124  if(den > 0)
125  return (DecayGenFrame::GetProbMCMC(mass)*mass*mass)/den;
126  else
127  return 0.;
128  }
129 
130  void ResonanceGenFrame::GenerateMassMCMC(double& mass, double&
131  prob, double max) const {
132  double min = 0.;
133  int N = GetNChildren();
134  for(int i = 0; i < N; i++)
135  min += GetChildFrame(i).GetMass();
136 
137  if(((max < min) && (max > 0)) || m_Width <= 0.){
138  mass = 0.;
139  prob = 1.;
140  return;
141  }
142 
143  if(min <= 0)
144  min = 0.;
145 
146  double M2 = m_PoleMass*m_PoleMass;
147  double MW = m_PoleMass*m_Width;
148  double Imin = atan((min*min-M2)/MW);
149  double Imax;
150  if(max <= 0)
151  Imax = TMath::Pi()/2.;
152  else
153  Imax = atan((max*max-M2)/MW);
154 
155  mass = sqrt(M2 + MW*tan(Imin+GetRandom()*(Imax-Imin)));
156 
157  double den = mass*mass-m_PoleMass*m_PoleMass;
158  den *= den;
159  den += m_PoleMass*m_PoleMass*m_Width*m_Width;
160  if(den > 0)
161  prob = 1./den;
162  else
163  prob = 0.;
164  }
165 
166  ResonanceGenFrame ResonanceGenFrame::m_Empty;
167 }
virtual double GetMass() const
Get the mass of this frame.
virtual GeneratorFrame & GetChildFrame(int i=0) const
Get the frame of the i th child.