LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
GeneratorFrame.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 <math.h>
31 #include <TRandom3.h>
32 #include <TDatime.h>
35 
36 namespace RestFrames {
37 
39  // GeneratorFrame class methods
41  GeneratorFrame::GeneratorFrame(const std::string& sname,
42  const std::string& stitle)
43  : RestFrame(sname, stitle)
44  {
45  m_Ana = kGenFrame;
46  m_VarMassMCMC = false;
47  m_Mass = 0.;
48 
49  m_Ngen = 0;
50  m_Npass = 0;
51 
52  m_PCut = -1.;
53  m_PtCut = -1.;
54  m_EtaCut = -1.;
55  m_minMassCut = -1.;
56  m_maxMassCut = -1.;
57  m_doCuts = false;
58  m_doPCut = false;
59  m_doPtCut = false;
60  m_doEtaCut = false;
61  m_dominMassCut = false;
62  m_domaxMassCut = false;
63 
64  TDatime now;
65  int today = now.GetDate();
66  int clock = now.GetTime();
67  int key = GetKey().GetKey();
68  int seed = today+clock+key;
69  m_Random = new TRandom3(seed);
70  }
71 
72  GeneratorFrame::GeneratorFrame()
73  : RestFrame() { }
74 
75  GeneratorFrame::~GeneratorFrame(){
76  if(m_Random) delete m_Random;
77  }
78 
80  RestFrame::Clear();
81  }
82 
85  return VisibleGenFrame::Empty();
86  }
87 
88  double GeneratorFrame::GetMass() const {
89  return std::max(m_Mass, 0.);
90  }
91 
93  if(!frame) return;
94  if(!frame.IsGenFrame()) return;
95  RestFrame::AddChildFrame(frame);
96  }
97 
99  if(!frame) return;
100  if(!frame.IsGenFrame()) return;
101  RestFrame::SetParentFrame(frame);
102  }
103 
105  const RestFrame& frame = RestFrame::GetParentFrame();
106  if(!frame.IsEmpty())
107  return static_cast<const GeneratorFrame&>(frame);
108  else
109  return GeneratorFrame::Empty();
110  }
111 
113  RestFrame& frame = RestFrame::GetChildFrame(i);
114  if(!frame.IsEmpty())
115  return static_cast<GeneratorFrame&>(frame);
116  else
117  return GeneratorFrame::Empty();
118  }
119 
120  bool GeneratorFrame::ClearEventRecursive(){
121  ResetGenFrame();
122  if(!IsSoundMind())
123  return false;
124 
125  int Nf = GetNChildren();
126  for(int i = 0; i < Nf; i++)
127  if(!GetChildFrame(i).ClearEventRecursive())
128  return false;
129 
130  return true;
131  }
132 
133  bool GeneratorFrame::AnalyzeEventRecursive(){
134  if(!IsSoundMind()){
135  UnSoundMind(RF_FUNCTION);
136  return SetSpirit(false);
137  }
138  if(!GenerateFrame()){
139  m_Log << LogWarning;
140  m_Log << "Unable to generate event for this frame.";
141  m_Log << LogEnd;
142  return SetSpirit(false);
143  }
144 
145  int Nchild = GetNChildren();
146  for(int i = 0; i < Nchild; i++)
147  if(!GetChildFrame(i).AnalyzeEventRecursive()){
148  return SetSpirit(false);
149  }
150 
151  return SetSpirit(true);
152  }
153 
154  void GeneratorFrame::SetChildren(const std::vector<TLorentzVector>& P_children){
155  int N = P_children.size();
156  for(int i = 0; i < N; i++){
157  TLorentzVector P = P_children[i];
158  TVector3 B_child = P.BoostVector();
159 
160  SetChildBoostVector(GetChildFrame(i), B_child);
161  GetChildFrame(i).SetFourVector(P,*this);
162  }
163  }
164 
165  bool GeneratorFrame::InitializeGenAnalysis(){
166  return SetMind(true);
167  }
168 
169  bool GeneratorFrame::InitializeAnalysisRecursive(){
170  if(!IsSoundBody()){
171  UnSoundBody(RF_FUNCTION);
172  return SetMind(false);
173  }
174 
175  if(!InitializeGenAnalysis())
176  return SetMind(false);
177 
178  m_Ngen = 0;
179  m_Npass = 0;
180 
181  int Nchild = GetNChildren();
182  for(int i = 0; i < Nchild; i++)
183  if(!GetChildFrame(i).InitializeAnalysisRecursive()){
184  m_Log << LogWarning;
185  m_Log << "Unable to recursively initialize analysis for frame:";
186  m_Log << Log(GetChildFrame(i)) << LogEnd;
187  return SetMind(false);
188  }
189 
190  return SetMind(true);
191  }
192 
193  double GeneratorFrame::GetRandom() const {
194  return m_Random->Rndm();
195  }
196 
197  double GeneratorFrame::GetGaus(double mu, double sig) const {
198  return m_Random->Gaus(mu,sig);
199  }
200 
201  bool GeneratorFrame::IterateMCMC(){
202  return true;
203  }
204 
205  bool GeneratorFrame::IterateRecursiveMCMC(){
206  if(!IsSoundMind()){
207  UnSoundMind(RF_FUNCTION);
208  return SetMind(false);
209  }
210 
211  if(!IterateMCMC())
212  return SetMind(false);
213 
214  int N = GetNChildren();
215  for(int i = 0; i < N; i++)
216  if(!GetChildFrame(i).IterateRecursiveMCMC())
217  return SetMind(false);
218 
219  return SetMind(true);
220  }
221 
222  void GeneratorFrame::SetVariableMassMCMC(bool var){
223  m_VarMassMCMC = var;
224  }
225 
227  return m_VarMassMCMC;
228  }
229 
230  double GeneratorFrame::GetMinimumMassMCMC() const {
231  if(!IsSoundBody()){
232  UnSoundBody(RF_FUNCTION);
233  return SetBody(false);
234  }
235 
236  double mass = 0.;
237  int N = GetNChildren();
238  for(int i = 0; i < N; i++)
239  mass += GetChildFrame(i).GetMinimumMassMCMC();
240 
241  if(!IsVariableMassMCMC())
242  mass = std::max(GetMass(),mass);
243 
244  return mass;
245  }
246 
247  void GeneratorFrame::GenerateMassMCMC(double& mass, double& prob,
248  double max) const {
249  mass = 0.;
250  prob = 1.;
251  }
252 
253  void GeneratorFrame::SetMassMCMC(double val){
254  if(val < 0.){
255  m_Log << LogWarning;
256  m_Log << "Unable to set mass to negative value ";
257  m_Log << val << ". Setting to zero." << LogEnd;
258  m_Mass = 0.;
259  } else {
260  m_Mass = val;
261  }
262  }
263 
264  void GeneratorFrame::SetMassMCMC(double mass,
265  GeneratorFrame& frame) const {
266  frame.SetMassMCMC(mass);
267  }
268 
269  double GeneratorFrame::GetProbMCMC(double mass) const {
270  return 1.;
271  }
272 
274  if(IsLabFrame()){
275  m_Log << LogInfo << std::endl;
276  m_Log << "Total events generated: " << m_Ngen << std::endl;
277  m_Log << "Events in acceptance: " << m_Npass << std::endl;
278  m_Log << "Generator efficiency: ";
279  m_Log << 100.*double(m_Npass)/double(m_Ngen) << " %";
280  m_Log << std::endl << LogEnd;
281  }
282 
283  if(m_doCuts){
284  m_Log << LogInfo;
285  m_Log << "Acceptance cuts for frame:" << std::endl;
286  if(m_dominMassCut || m_domaxMassCut){
287  m_Log << " ";
288  if(m_dominMassCut)
289  m_Log << m_minMassCut << " < ";
290  m_Log << "mass";
291  if(m_domaxMassCut)
292  m_Log << " < " << m_maxMassCut;
293  m_Log << std::endl;
294  }
295  if(m_doPCut){
296  m_Log << " P > " << m_PCut << std::endl;
297  }
298  if(m_doPtCut){
299  m_Log << " Pt > " << m_PtCut << std::endl;
300  }
301  if(m_doEtaCut){
302  m_Log << " |Eta| < " << m_EtaCut << std::endl;
303  }
304  m_Log << "Acceptance efficiency = ";
305  m_Log << 100.*double(m_Npass)/double(m_Ngen) << " %";
306  m_Log << std::endl << LogEnd;
307  }
308 
309  int N = GetNChildren();
310  for(int i = 0; i < N; i++)
312  }
313 
314  bool GeneratorFrame::EventInAcceptance() const {
315  if(!IsSoundSpirit()){
316  UnSoundSpirit(RF_FUNCTION);
317  return SetSpirit(false);
318  }
319 
320  bool pass = true;
321  if(m_doCuts){
322  TLorentzVector P = GetFourVector();
323  if(m_doPCut)
324  if(P.P() < (1.-1e-10)*m_PCut)
325  pass = false;
326  if(m_doPtCut)
327  if(P.Pt() < (1.-1e-10)*m_PtCut)
328  pass = false;
329  if(m_doEtaCut)
330  if(fabs(P.Eta()) > (1.+1e-10)*m_EtaCut)
331  pass = false;
332  if(m_dominMassCut)
333  if(m_Mass < (1.-1e-10)*m_minMassCut)
334  pass = false;
335  if(m_domaxMassCut)
336  if(m_Mass > (1.+1e-10)*m_maxMassCut)
337  pass = false;
338  }
339 
340  bool evt_pass = pass;
341  bool pass_c;
342  int N = GetNChildren();
343  for(int i = 0; i < N; i++){
344  pass_c = GetChildFrame(i).EventInAcceptance();
345  evt_pass = evt_pass && pass_c;
346  }
347 
348  if(IsLabFrame())
349  pass = evt_pass;
350 
351  m_Ngen++;
352  if(pass)
353  m_Npass++;
354 
355  return evt_pass;
356  }
357 
358  void GeneratorFrame::SetPCut(double cut){
359  if(cut <= 0.) return;
360  m_PCut = cut;
361  m_doPCut = true;
362  m_doCuts = true;
363  }
364 
365  void GeneratorFrame::SetPtCut(double cut){
366  if(cut <= 0.) return;
367  m_PtCut = cut;
368  m_doPtCut = true;
369  m_doCuts = true;
370  }
371 
372  void GeneratorFrame::SetEtaCut(double cut){
373  m_EtaCut = fabs(cut);
374  m_doEtaCut = true;
375  m_doCuts = true;
376  }
377 
378  void GeneratorFrame::SetMassWindowCut(double min, double max){
379  if(min > 0 || max > 0)
380  m_doCuts = true;
381  else
382  return;
383 
384  if(min > 0){
385  m_dominMassCut = true;
386  m_minMassCut = min;
387  }
388 
389  if(max > 0){
390  m_domaxMassCut = true;
391  m_maxMassCut = max;
392  }
393  }
394 
395  void GeneratorFrame::RemovePCut(){
396  m_doPCut = false;
397  m_doCuts = m_doPCut || m_doPtCut || m_doEtaCut ||
398  m_dominMassCut || m_domaxMassCut;
399  }
400 
401  void GeneratorFrame::RemovePtCut(){
402  m_doPtCut = false;
403  m_doCuts = m_doPCut || m_doPtCut || m_doEtaCut ||
404  m_dominMassCut || m_domaxMassCut;
405  }
406 
407  void GeneratorFrame::RemoveEtaCut(){
408  m_doEtaCut = false;
409  m_doCuts = m_doPCut || m_doPtCut || m_doEtaCut ||
410  m_dominMassCut || m_domaxMassCut;
411  }
412 
413  void GeneratorFrame::RemoveMassWindowCut(){
414  m_dominMassCut = false;
415  m_domaxMassCut = false;
416  m_doCuts = m_doPCut || m_doPtCut || m_doEtaCut ||
417  m_dominMassCut || m_domaxMassCut;
418  }
419 
420 }
virtual double GetMass() const
Get the mass of this frame.
virtual GeneratorFrame const & GetParentFrame() const
Returns the parent of this frame.
bool IsVariableMassMCMC() const
Frame is capable having a variable mass? (true/false)
virtual void AddChildFrame(RestFrame &frame)
Add a child RestFrame to this frame.
virtual void SetParentFrame(RestFrame &frame=RestFrame::Empty())
Set the parent frame for this frame.
virtual void Clear()
Clears GeneratorFrame of all connections to other objects.
RFKey GetKey() const
gets object identification key
Definition: RFBase.cc:100
abstract base class for all Frame objects
static GeneratorFrame & Empty()
Returns empty instance of class.
virtual GeneratorFrame & GetChildFrame(int i=0) const
Get the frame of the i th child.
static VisibleGenFrame & Empty()
Returns empty instance of class.
void PrintGeneratorEfficiency() const
Print generator efficiency information.