LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
ppLabGenFrame.cc
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2015, 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 
31 
32 namespace RestFrames {
33 
35  // ppLabGenFrame class
37 
38  ppLabGenFrame::ppLabGenFrame(const std::string& sname,
39  const std::string& stitle) :
40  LabGenFrame(sname, stitle)
41  {
42  m_deltaLogX = 0;
43 
44  // default is 13 TeV LHC
45  m_Ep1 = 6500.;
46  m_Ep2 = 6500.;
47 
48  // default parton initial state is q-qbar
49  m_PDFqqbar = true;
50  m_PDFgg = false;
51  m_PDFgq = false;
52  m_PDFqq = false;
53  }
54 
55  ppLabGenFrame::~ppLabGenFrame() {}
56 
58  m_deltaLogX = 0;
60  }
61 
62  void ppLabGenFrame::SetPDFqqbar(){
63  SetMind(false);
64  m_PDFqqbar = true;
65  m_PDFgg = false;
66  m_PDFgq = false;
67  m_PDFqq = false;
68  }
69 
70  void ppLabGenFrame::SetPDFgg(){
71  SetMind(false);
72  m_PDFqqbar = false;
73  m_PDFgg = true;
74  m_PDFgq = false;
75  m_PDFqq = false;
76  }
77 
78  void ppLabGenFrame::SetPDFgq(){
79  SetMind(false);
80  m_PDFqqbar = false;
81  m_PDFgg = false;
82  m_PDFgq = true;
83  m_PDFqq = false;
84  }
85 
86  void ppLabGenFrame::SetPDFqq(){
87  SetMind(false);
88  m_PDFqqbar = false;
89  m_PDFgg = false;
90  m_PDFgq = false;
91  m_PDFqq = true;
92  }
93 
94  void ppLabGenFrame::SetEnergyP1(double E){
95  SetMind(false);
96  if(E > 0)
97  m_Ep1 = E;
98  }
99 
100  void ppLabGenFrame::SetEnergyP2(double E){
101  SetMind(false);
102  if(E > 0)
103  m_Ep2 = E;
104  }
105 
106  double ppLabGenFrame::GetEnergyP1() const {
107  return m_Ep1;
108  }
109 
110  double ppLabGenFrame::GetEnergyP2() const {
111  return m_Ep2;
112  }
113 
114  bool ppLabGenFrame::InitializeGenAnalysis(){
115  if(!IsSoundBody()){
116  UnSoundBody(RF_FUNCTION);
117  return SetMind(false);
118  }
119 
120  double sqrtS = 2.*sqrt(m_Ep1*m_Ep2);
121  m_MaxM = sqrtS;
122 
123  double Mmin = GetMinimumMassMCMC();
124  if(Mmin > sqrtS){
125  m_Log << LogWarning;
126  m_Log << "Unable to initialize event generation. ";
127  m_Log << "sqrt(S) of " << sqrtS << " ";
128  m_Log << "is less than minimum child mass of " << Mmin;
129  m_Log << LogEnd;
130  return SetMind(false);
131  }
132 
133  m_deltaLogX = 0.;
134 
135  if(!LabGenFrame::InitializeGenAnalysis())
136  return SetMind(false);
137 
138  return SetMind(true);
139  }
140 
141  bool ppLabGenFrame::IterateMCMC(){
142  double deltaLogX = GetRandom()*2.-1.;
143  double deltaLogXOld = m_deltaLogX;
144 
145  double probOld = GetProbMCMC();
146  m_deltaLogX = deltaLogX;
147  double probNew = GetProbMCMC();
148 
149  if(probOld > 0)
150  if(probNew/probOld < GetRandom())
151  m_deltaLogX = deltaLogXOld;
152 
153  LabGenFrame::IterateMCMC();
154 
155  double C = m_ChildMassMCMC*m_ChildMassMCMC/4./m_Ep1/m_Ep2;
156  if(C > 0.){
157  double expo = -m_deltaLogX*log(C);
158  double Xp1 = sqrt(C*exp(expo));
159  double Xp2 = sqrt(C*exp(-expo));
160  SetLongitudinalMomentum(m_Ep1*Xp1 - m_Ep2*Xp2);
161  } else
162  SetLongitudinalMomentum(0.);
163 
164  return SetMind(true);
165  }
166 
167  double ppLabGenFrame::GetProbMCMC(double mass) const {
168  if(mass < 0)
169  mass = GetChildFrame().GetMass();
170 
171  double C = mass*mass/(4.*m_Ep1*m_Ep2);
172  double expo = -m_deltaLogX*log(C);
173 
174  double Xp1 = sqrt(C*exp(expo));
175  double Xp2 = sqrt(C*exp(-expo));
176 
177  double prob = 1.;
178  if(m_PDFqqbar)
179  prob *= pPDF_q(Xp1)*pPDF_qbar(Xp2)+pPDF_q(Xp2)*pPDF_qbar(Xp1);
180  if(m_PDFgg)
181  prob *= pPDF_g(Xp1)*pPDF_g(Xp2);
182  if(m_PDFgq)
183  prob *= pPDF_q(Xp1)*pPDF_g(Xp2)+pPDF_q(Xp2)*pPDF_g(Xp1);
184  if(m_PDFqq)
185  prob *= pPDF_q(Xp1)*pPDF_q(Xp2);
186 
187  return prob;
188  }
189 
190  double ppLabGenFrame::pPDF_q(double x) const {
191  if(x <= 0. || x >= 1.)
192  return 0.;
193  return fabs((1/x)*pow(x,m_PDF_eta_1)*pow(1-x,m_PDF_eta_2)*
194  (1.+m_PDF_eps_u*sqrt(x)+m_PDF_g_u*x));
195  }
196 
197  double ppLabGenFrame::pPDF_qbar(double x) const {
198  if(x <= 0. || x >= 1.)
199  return 0.;
200  return fabs((1/x)*pow(x,m_PDF_del_S)*pow(1-x,m_PDF_eta_S)*
201  (1.+m_PDF_eps_S*sqrt(x)+m_PDF_g_S*x));
202  }
203 
204  double ppLabGenFrame::pPDF_g(double x) const {
205  if(x <= 0. || x >= 1.)
206  return 0.;
207  return fabs((1/x)*(m_PDF_A_g*pow(x,m_PDF_del_g)*pow(1-x,m_PDF_eta_g)*
208  (1.+m_PDF_eps_g*sqrt(x)+m_PDF_g_g*x)+
209  m_PDF_A_g1*pow(x,m_PDF_del_g1)*pow(1-x,m_PDF_eta_g1)));
210  }
211 
212  double ppLabGenFrame::m_PDF_eta_1 = 0.27871;
213  double ppLabGenFrame::m_PDF_eta_2 = 3.3627;
214  double ppLabGenFrame::m_PDF_eps_u = 4.4343;
215  double ppLabGenFrame::m_PDF_g_u = 38.599;
216  double ppLabGenFrame::m_PDF_del_S = -0.11912;
217  double ppLabGenFrame::m_PDF_eta_S = 9.4189;
218  double ppLabGenFrame::m_PDF_eps_S = -2.6287;
219  double ppLabGenFrame::m_PDF_g_S = 18.065;
220  double ppLabGenFrame::m_PDF_A_g = 3.4055;
221  double ppLabGenFrame::m_PDF_del_g = -0.12178;
222  double ppLabGenFrame::m_PDF_eta_g = 2.9278;
223  double ppLabGenFrame::m_PDF_eps_g = -2.3210;
224  double ppLabGenFrame::m_PDF_g_g = 1.9233;
225  double ppLabGenFrame::m_PDF_A_g1 = -1.6189;
226  double ppLabGenFrame::m_PDF_del_g1 = -0.23999;
227  double ppLabGenFrame::m_PDF_eta_g1 = 24.792;
228 
229 }
virtual void Clear()
Clears GeneratorFrame of all connections to other objects.
virtual void Clear()
Clears GeneratorFrame of all connections to other objects.
Definition: LabGenFrame.cc:56
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.