LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
LabGenFrame.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 
31 
32 namespace RestFrames {
33 
35  // LabGenFrame class
37  LabGenFrame::LabGenFrame(const std::string& sname,
38  const std::string& stitle) :
39  LabFrame<GeneratorFrame>(sname, stitle)
40  {
41  m_PT = 0.;
42  m_PToM = -1.;
43  m_PL = 0.;
44  m_Phi = -1.;
45 
46  m_MaxM = -1.;
47 
48  m_NBurnInMCMC = 1000;
49  m_NDiscardMCMC = 5;
50 
51  m_FailTolerance = 1000;
52  }
53 
54  LabGenFrame::~LabGenFrame() {}
55 
58  }
59 
60  void LabGenFrame::SetThreeVector(const TVector3& P){
61  m_PT = P.Pt();
62  m_PL = P.Z();
63  SetPhi(P.Phi());
64  }
65 
66  void LabGenFrame::SetPToverM(double val){
67  if(val < 0.){
68  m_Log << LogWarning;
69  m_Log << "Unable to set transverse momentum ";
70  m_Log << "to negative value: " << val << "*mass";
71  m_Log << LogEnd;
72  } else {
73  m_PToM = val;
74  m_PT = 0.;
75  }
76  }
77 
78  void LabGenFrame::SetTransverseMomentum(double val){
79  if(val < 0.){
80  m_Log << LogWarning;
81  m_Log << "Unable to set transverse momentum ";
82  m_Log << "to negative value: " << val;
83  m_Log << LogEnd;
84  } else {
85  m_PT = val;
86  m_PToM = -1.;
87  }
88  }
89 
90  void LabGenFrame::SetLongitudinalMomentum(double val){
91  m_PL = val;
92  }
93 
94  void LabGenFrame::SetPhi(double val){
95  while(val > acos(-1.)*2.) val -= acos(-1.)*2.;
96  while(val < 0.) val += acos(-1.)*2.;
97  m_Phi = val;
98  }
99 
100  void LabGenFrame::ResetGenFrame(){
101  SetSpirit(false);
102  m_Phi = -1.;
103  }
104 
105  bool LabGenFrame::InitializeAnalysis(){
106  m_Log << LogVerbose << "Initializing this tree for analysis..." << LogEnd;
107 
108  if(!IsSoundBody()){
109  UnSoundBody(RF_FUNCTION);
110  return SetMind(false);
111  }
112 
113  if(!InitializeAnalysisRecursive()){
114  m_Log << LogWarning << "...Unable to recursively initialize analysis" << LogEnd;
115  return SetMind(false);
116  }
117 
118  for(int i = 0; i < m_NBurnInMCMC; i++)
119  if(!IterateRecursiveMCMC()){
120  m_Log << LogWarning << "...Unable to recursively initialize analysis" << LogEnd;
121  return SetMind(false);
122  }
123 
124  m_Log << LogVerbose << "...Done" << LogEnd;
125  return SetMind(true);
126  }
127 
128  void LabGenFrame::SetN_MCMCBurnIn(int N){
129  SetMind(false);
130  m_NBurnInMCMC = std::max(0,N);
131  }
132 
133  void LabGenFrame::SetN_MCMCDiscard(int N){
134  SetMind(false);
135  m_NDiscardMCMC = std::max(1,N);
136  }
137 
138  void LabGenFrame::SetFailTolerance(int Nfail){
139  m_FailTolerance = Nfail;
140  }
141 
142  bool LabGenFrame::InitializeGenAnalysis(){
143  if(!IsSoundBody()){
144  UnSoundBody(RF_FUNCTION);
145  return SetMind(false);
146  }
147 
148  GeneratorFrame& child = GetChildFrame();
149  if(child.IsVariableMassMCMC()){
150  double ChildMass, ChildProb;
151  child.GenerateMassMCMC(ChildMass, ChildProb, m_MaxM);
152  m_ChildMassMCMC = ChildMass;
153  m_ChildProbMCMC = ChildProb;
154  SetMassMCMC(ChildMass, child);
155  } else {
156  m_ChildMassMCMC = child.GetMass();
157  m_ChildProbMCMC = 1.;
158  }
159 
160  return SetMind(true);
161  }
162 
163  bool LabGenFrame::IterateMCMC(){
164  GeneratorFrame& child = GetChildFrame();
165  if(child.IsVariableMassMCMC()){
166  double ChildMass, ChildProb = 0.;
167  child.GenerateMassMCMC(ChildMass, ChildProb, m_MaxM);
168 
169  double probOld = GetProbMCMC(m_ChildMassMCMC)*
170  child.GetProbMCMC(m_ChildMassMCMC)/m_ChildProbMCMC;
171 
172  double probNew = GetProbMCMC(ChildMass)*
173  child.GetProbMCMC(ChildMass)/ChildProb;
174 
175  if(probNew/probOld > GetRandom()){
176  m_ChildMassMCMC = ChildMass;
177  m_ChildProbMCMC = ChildProb;
178  SetMassMCMC(ChildMass, child);
179  } else {
180  SetMassMCMC(m_ChildMassMCMC, child);
181  }
182  }
183 
184  return SetMind(true);
185  }
186 
187  bool LabGenFrame::GenerateFrame(){
188  if(!IsSoundBody())
189  return false;
190 
191  TLorentzVector P;
192  double M = GetChildFrame().GetMass();
193  if(m_Phi < 0.) m_Phi = 2.*acos(-1.)*GetRandom();
194 
195  if(m_PToM > 0.)
196  P.SetPxPyPzE(m_PToM*M*cos(m_Phi), m_PToM*M*sin(m_Phi), m_PL,
197  sqrt(m_PL*m_PL + M*M*(1. + m_PToM*m_PToM)));
198  else
199  P.SetPxPyPzE(m_PT*cos(m_Phi), m_PT*sin(m_Phi), m_PL,
200  sqrt(m_PT*m_PT + m_PL*m_PL + M*M));
201  m_Phi = -1.;
202 
203  std::vector<TLorentzVector> ChildVector;
204  ChildVector.push_back(P);
205  SetChildren(ChildVector);
206 
207  return SetSpirit(true);
208  }
209 
210  bool LabGenFrame::ClearEvent(){
211  SetSpirit(false);
212  if(!IsSoundBody())
213  return false;
214  return ClearEventRecursive();
215  }
216 
217  bool LabGenFrame::AnalyzeEvent(){
218  bool pass = false;
219  int tries = 0;
220 
221  while(!pass){
222  for(int i = 0; i < m_NDiscardMCMC; i++)
223  if(!IterateRecursiveMCMC())
224  return SetSpirit(false);
225 
226  if(!AnalyzeEventRecursive()){
227  return SetSpirit(false);
228  }
229 
230  pass = EventInAcceptance();
231 
232  if(!pass){
233  tries++;
234  if(tries > m_FailTolerance &&
235  m_FailTolerance > 0){
236  m_Log << LogWarning;
237  m_Log << "Failed to generate event in ";
238  m_Log << "acceptance in " << tries;
239  m_Log << " tries. Giving up." << LogEnd;
240  return SetSpirit(false);
241  }
242  }
243  }
244 
245  return SetSpirit(true);
246  }
247 
248 }
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 void Clear()
Clears GeneratorFrame of all connections to other objects.
virtual GeneratorFrame & GetChildFrame(int i=0) const
Get the frame of the i th child.