LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
CombinatoricJigsaw.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 
32 
33 namespace RestFrames {
34 
35  CombinatoricJigsaw::CombinatoricJigsaw(const std::string& sname,
36  const std::string& stitle,
37  int Ncomb, int Nobject)
38  : Jigsaw(sname, stitle, Ncomb, Nobject),
39  m_Ncomb(Ncomb), m_Nobj(Nobject)
40  {
41  m_Type = kCombinatoricJigsaw;
42  for(int i = 0; i < m_Ncomb; i++){
43  m_ChildStates += GetNewChildState();
44  }
45  }
46 
47  CombinatoricJigsaw::CombinatoricJigsaw()
48  : Jigsaw(), m_Ncomb(0), m_Nobj(0) {}
49 
50  CombinatoricJigsaw::~CombinatoricJigsaw() {}
51 
53  Jigsaw::Clear();
54  }
55 
56  CombinatoricJigsaw& CombinatoricJigsaw::Empty(){
57  return MinMassesCombJigsaw::Empty();
58  }
59 
60  void CombinatoricJigsaw::SetGroup(Group& group){
61  if(!group) return;
62  if(!group.IsCombinatoricGroup()) return;
63  Jigsaw::SetGroup(group);
64  }
65 
66  CombinatoricGroup& CombinatoricJigsaw::GetGroup() const {
67  if(!Jigsaw::GetGroup())
68  return CombinatoricGroup::Empty();
69  else
70  return static_cast<CombinatoricGroup&>(Jigsaw::GetGroup());
71  }
72 
73  void CombinatoricJigsaw::SetParentState(State& state){
74  if(!state) return;
75  if(!state.IsCombinatoricState()) return;
76  Jigsaw::SetParentState(state);
77  }
78 
79  CombinatoricState const& CombinatoricJigsaw::GetParentState() const {
80  if(!Jigsaw::GetParentState())
81  return CombinatoricState::Empty();
82  else
83  return static_cast<const CombinatoricState&>(Jigsaw::GetParentState());
84  }
85 
86  CombinatoricState& CombinatoricJigsaw::GetChildState(int i) const {
87  if(!Jigsaw::GetChildState(i))
88  return CombinatoricState::Empty();
89  else
90  return static_cast<CombinatoricState&>(Jigsaw::GetChildState(i));
91  }
92 
93  void CombinatoricJigsaw::AddCombFrame(const RestFrame& frame, int i){
94  if(!frame) return;
95 
96  ConstRestFrameList frames =
97  frame.GetListVisibleFrames();
98  int N = frames.GetN();
99  for(int f = 0; f < N; f++)
100  AddChildFrame(frames[f], i);
101  }
102 
103  void CombinatoricJigsaw::AddCombFrames(const ConstRestFrameList& frames, int i){
104  int N = frames.GetN();
105  for(int f = 0; f < N; f++)
106  AddCombFrame(frames[f], i);
107  }
108 
109  void CombinatoricJigsaw::AddObjectFrame(const RestFrame& frame, int i){
110  if(!frame) return;
111 
112  ConstRestFrameList frames =
113  frame.GetListVisibleFrames()+
114  frame.GetListInvisibleFrames();
115  int N = frames.GetN();
116  for(int f = 0; f < N; f++)
117  AddDependancyFrame(frames[f], i);
118  }
119 
120  void CombinatoricJigsaw::AddObjectFrames(const ConstRestFrameList& frames, int i){
121  int N = frames.GetN();
122  for(int f = 0; f < N; f++)
123  AddObjectFrame(frames[f], i);
124  }
125 
126  void CombinatoricJigsaw::SetCombCharge(const RFCharge& charge, int i){
127  if(i < 0 || i >= m_Ncomb)
128  return;
129 
130  m_ChargeForChild[i] = charge;
131  }
132 
133  void CombinatoricJigsaw::SetCombCharge(int charge, int i){
134  SetCombCharge(RFCharge(charge), i);
135  }
136 
137  void CombinatoricJigsaw::SetCombCharge(int charge_num, int charge_den, int i){
138  SetCombCharge(RFCharge(charge_num, charge_den), i);
139  }
140 
141  void CombinatoricJigsaw::UnsetCombCharge(int i){
142  if(m_ChargeForChild.count(i) > 0)
143  m_ChargeForChild.erase(i);
144  }
145 
146  void CombinatoricJigsaw::SetObjectCharge(const RFCharge& charge, int i){
147  if(i < 0 || i >= m_Nobj)
148  return;
149 
150  m_ChargeForObject[i] = charge;
151  }
152 
153  void CombinatoricJigsaw::SetObjectCharge(int charge, int i){
154  SetObjectCharge(RFCharge(charge), i);
155  }
156 
157  void CombinatoricJigsaw::SetObjectCharge(int charge_num, int charge_den, int i){
158  SetObjectCharge(RFCharge(charge_num, charge_den), i);
159  }
160 
161  void CombinatoricJigsaw::UnsetObjectCharge(int i){
162  if(m_ChargeForObject.count(i) > 0)
163  m_ChargeForObject.erase(i);
164  }
165 
166  bool CombinatoricJigsaw::InitializeJigsawExecutionList(JigsawList& exec_jigsaws){
167  if(!IsSoundMind()) return false;
168  if(exec_jigsaws.Contains(*this)) return true;
169 
170  m_ExecuteJigsaws.Clear();
171 
172  // Add group dependancy jigsaws first
173  JigsawList group_jigsaws;
174  FillGroupJigsawDependancies(group_jigsaws);
175  group_jigsaws -= *this;
176 
177  int Ngroup = group_jigsaws.GetN();
178  for(int i = Ngroup-1; i >= 0; i--){
179  Jigsaw& jigsaw = group_jigsaws[i];
180  m_DependancyJigsaws -= jigsaw;
181  if(!exec_jigsaws.Contains(jigsaw))
182  if(!jigsaw.InitializeJigsawExecutionList(exec_jigsaws))
183  return SetMind(false);
184  }
185  // Satisfy dependancy jigsaws
186  while(m_DependancyJigsaws.GetN() > 0){
187  Jigsaw& jigsaw = m_DependancyJigsaws[m_DependancyJigsaws.GetN()-1];
188  if(exec_jigsaws.Contains(jigsaw)){
189  m_DependancyJigsaws -= jigsaw;
190  continue;
191  }
192  if(!jigsaw.DependsOnJigsaw(*this)){
193  if(!jigsaw.InitializeJigsawExecutionList(exec_jigsaws)){
194  return SetMind(false);
195  }
196  m_DependancyJigsaws -= jigsaw;
197  continue;
198  }
199  JigsawList temp_exec_jigsaws = exec_jigsaws;
200  temp_exec_jigsaws += m_ExecuteJigsaws;
201  temp_exec_jigsaws += *this;
202  if(!jigsaw.InitializeJigsawExecutionList(temp_exec_jigsaws))
203  return SetMind(false);
204  temp_exec_jigsaws -= *this;
205  temp_exec_jigsaws -= exec_jigsaws;
206 
207  m_DependancyJigsaws -= temp_exec_jigsaws;
208  m_ExecuteJigsaws += temp_exec_jigsaws;
209  }
210 
211  exec_jigsaws += *this;
212  exec_jigsaws += m_ExecuteJigsaws;
213  return true;
214  }
215 
216  CombinatoricState& CombinatoricJigsaw::GetNewChildState(){
217  char strn[10];
218  sprintf(strn,"%d",GetNChildren());
219  std::string name = GetName()+"_"+std::string(strn);
220  CombinatoricState* statePtr = new CombinatoricState(name, name);
221  AddDependent(statePtr);
222  return *statePtr;
223  }
224 
225  bool CombinatoricJigsaw::ExecuteDependancyJigsaws(){
226  int N = m_ExecuteJigsaws.GetN();
227  for(int i = 0; i < N; i++)
228  if(!m_ExecuteJigsaws[i].AnalyzeEvent())
229  return false;
230  return true;
231  }
232 
233  bool CombinatoricJigsaw::AnalyzeEvent(){
234  if(!IsSoundMind() || !GetGroup())
235  return SetSpirit(false);
236 
237  if(!InitializeCombinatoric()){
238  m_Log << LogWarning;
239  m_Log << "Problem initializing event info" << LogEnd;
240  return SetSpirit(false);
241  }
242 
243  if(!LoopCombinatoric()){
244  m_Log << LogWarning;
245  m_Log << "Problem looping over combinatorics" << LogEnd;
246  return SetSpirit(false);
247  }
248 
249  return SetSpirit(true);
250  }
251 
252  bool CombinatoricJigsaw::InitializeAnalysis(){
253  if(!Jigsaw::InitializeAnalysis())
254  return SetMind(false);
255 
256  CombinatoricGroup& group = GetGroup();
257 
258  m_NForChild.clear();
259  m_NExclusive.clear();
260 
261  m_NinputTOT = 0;
262  m_NExclusiveTOT = true;
263 
264  for(int i = 0; i < m_Ncomb; i++){
265  RestFrameList const& frames = GetChildState(i).GetListFrames();
266  int Nf = frames.GetN();
267  int NTOT = 0;
268  bool exclTOT = true;
269  for(int f = 0; f < Nf; f++){
270  int N = -1;
271  bool excl = false;
272  group.GetNElementsForFrame(frames[f], N, excl);
273  if(N < 0) return SetMind(false);
274  NTOT += N;
275  exclTOT = exclTOT && excl;
276  }
277  m_NForChild.push_back(NTOT);
278  m_NExclusive.push_back(exclTOT);
279  m_NinputTOT += NTOT;
280  m_NExclusiveTOT = m_NExclusiveTOT && exclTOT;
281  }
282 
283  return SetMind(true);
284  }
285 
286  bool CombinatoricJigsaw::InitializeCombinatoric(){
287  if(!IsSoundMind())
288  return SetSpirit(false);
289 
290  if(!GetGroup())
291  return SetSpirit(false);
292 
293  if(!GetParentState())
294  return SetSpirit(false);
295 
296  m_InputStates.Clear();
297  m_InputStates = GetParentState().GetElements();
298 
299  if(m_InputStates.GetN() < m_NinputTOT){
300  m_Log << LogWarning;
301  m_Log << "Unable to execute Jigsaw. ";
302  m_Log << "Insufficienct number of inputs: ";
303  m_Log << m_NinputTOT << " (required) != ";
304  m_Log << m_InputStates.GetN() << " (provided)";
305  m_Log << LogEnd;
306  return SetSpirit(false);
307  }
308 
309  if(m_NExclusiveTOT &&
310  (m_InputStates.GetN() != m_NinputTOT)){
311  m_Log << LogWarning;
312  m_Log << "Unable to execute Jigsaw. ";
313  m_Log << "Incorrect number of exclusive inputs: ";
314  m_Log << m_NinputTOT << " (required) != ";
315  m_Log << m_InputStates.GetN() << " (provided)";
316  m_Log << LogEnd;
317  return SetSpirit(false);
318  }
319 
320  return SetSpirit(true);
321  }
322 
323  int CombinatoricJigsaw::GetNInputStates() const {
324  return m_InputStates.GetN();
325  }
326 
327  VisibleState& CombinatoricJigsaw::GetInputState(int i) const {
328  return m_InputStates[i];
329  }
330 
331  int CombinatoricJigsaw::GetNinputForChild(int i) const {
332  if(i < 0 || i >= m_Ncomb)
333  return 0;
334  return m_NForChild[i];
335  }
336 
337  bool CombinatoricJigsaw::IsNinputExclForChild(int i) const {
338  if(i < 0 || i >= m_Ncomb)
339  return false;
340  return m_NExclusive[i];
341  }
342 
343  bool CombinatoricJigsaw::IsChargeSetForChild(int i) const {
344  return !(m_ChargeForChild.count(i) == 0);
345  }
346 
347  RFCharge CombinatoricJigsaw::GetChargeForChild(int i) const {
348  if(IsChargeSetForChild(i))
349  return m_ChargeForChild[i];
350  else
351  return RFCharge();
352  }
353 
354  bool CombinatoricJigsaw::IsChargeSetForObject(int i) const {
355  return !(m_ChargeForObject.count(i) == 0);
356  }
357 
358  RFCharge CombinatoricJigsaw::GetChargeForObject(int i) const {
359  if(IsChargeSetForObject(i))
360  return m_ChargeForObject[i];
361  else
362  return RFCharge();
363  }
364 
365  bool CombinatoricJigsaw::LoopCombinatoric(){
366  int Ninput = m_InputStates.GetN();
367 
368  int N_comb = 1;
369  for(int i = 0; i < Ninput; i++) N_comb *= m_Ncomb;
370 
371  int c_min = -1;
372  double metric_min = -1;
373 
374  for(int c = 0; c < N_comb; c++){
375  int key = c;
376  for(int i = 0; i < m_Ncomb; i++)
377  GetChildState(i).ClearElements();
378 
379  // set output states for combinatoric;
380  for(int i = 0; i < Ninput; i++){
381  int ihem = key%m_Ncomb;
382  key /= m_Ncomb;
383  GetChildState(ihem).AddElement(m_InputStates[i]);
384  }
385 
386  // check validity of combinatoric
387  bool valid = true;
388  for(int i = 0; i < m_Ncomb; i++){
389  if(IsNinputExclForChild(i)){
390  if(GetChildState(i).GetNElements() != GetNinputForChild(i)){
391  valid = false;
392  break;
393  }
394  } else {
395  if(GetChildState(i).GetNElements() < GetNinputForChild(i)){
396  valid = false;
397  break;
398  }
399  }
400  if(IsChargeSetForChild(i)){
401  if(GetChildState(i).GetCharge() != GetChargeForChild(i)){
402  valid = false;
403  break;
404  }
405  }
406  }
407  if(!valid)
408  continue;
409 
410  // Execute depedancy Jigsaws for this combintoric
411  if(!ExecuteDependancyJigsaws())
412  continue;
413 
414  // check validity of objects
415  for(int i = 0; i < m_Nobj; i++){
416  if(IsChargeSetForObject(i)){
417  if(GetDependancyStates(i).GetCharge() != GetChargeForObject(i)){
418  valid = false;
419  break;
420  }
421  }
422  }
423  if(!valid)
424  continue;
425 
426  // Evaluate metric for this combinatoric
427  double metric;
428  if(!EvaluateMetric(metric))
429  continue;
430 
431  if(metric < metric_min || c_min < 0){
432  metric_min = metric;
433  c_min = c;
434  }
435 
436  if((metric < metric_min && metric >= 0.) || c_min < 0){
437  metric_min = metric;
438  c_min = c;
439  }
440  }
441 
442  if(c_min < 0){
443  m_Log << LogWarning;
444  m_Log << "Unable to find combinatoric satisfying ";
445  m_Log << "requested conditions." << LogEnd;
446  return SetSpirit(false);
447  }
448 
449  // Set outputs to best combinatoric
450  for(int i = 0; i < m_Ncomb; i++)
451  GetChildState(i).ClearElements();
452  int key = c_min;
453  for(int i = 0; i < Ninput; i++){
454  int ihem = key%m_Ncomb;
455  key /= m_Ncomb;
456  GetChildState(ihem).AddElement(m_InputStates[i]);
457  }
458 
459  // Execute depedancy Jigsaws
460  if(!ExecuteDependancyJigsaws())
461  return SetSpirit(false);
462 
463  return SetSpirit(true);
464  }
465 
466  bool CombinatoricJigsaw::IsSoundBody() const {
467  if(RFBase::IsSoundBody())
468  return true;
469 
470  if(!Jigsaw::IsSoundBody())
471  return SetBody(false);
472 
473  for(int i = 0; i < m_Nobj; i++){
474  if(GetDependancyFrames(i).GetN() <= 0){
475  m_Log << LogWarning;
476  m_Log << "Empty collection of object frames: " << i;
477  m_Log << LogEnd;
478  return SetBody(false);
479  }
480  }
481 
482  return SetBody(true);
483  }
484 
485 }
std::string GetName() const
Returns object name.
Definition: RFBase.cc:104
virtual void Clear()
Clears Jigsaw of all connections to other objects.
void AddDependent(RFBase *dep)
pointer to RFBase object owned by this one
Definition: RFBase.cc:88
virtual void Clear()
Clears Jigsaw of all connections to other objects.
Definition: Jigsaw.cc:74