LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
SelfAssemblingRecoFrame.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 
35 
36 namespace RestFrames {
37 
39  // SelfAssemblingRecoFrame class
41  SelfAssemblingRecoFrame::SelfAssemblingRecoFrame(const std::string& sname,
42  const std::string& stitle)
43  : DecayRecoFrame(sname,stitle)
44  {
45  m_RType = RDSelfAssembling;
46  m_IsAssembled = false;
47  m_Nvisible = 0;
48  m_Ndecay = 0;
49  m_NewEvent = true;
50  }
51 
52  SelfAssemblingRecoFrame::~SelfAssemblingRecoFrame() {}
53 
55  if(m_IsAssembled) Disassemble();
56  m_VisibleFrames.Clear();
57  m_DecayFrames.Clear();
59  }
60 
61  bool SelfAssemblingRecoFrame::ResetRecoFrame(){
62  if(!IsSoundMind()){
63  UnSoundMind(RF_FUNCTION);
64  return SetSpirit(false);
65  }
66  if(m_IsAssembled) Disassemble();
67  m_NewEvent = true;
68  return SetMind(true);
69  }
70 
72  m_ChildFrames_UnAssembled.Remove(frame);
74  }
75 
76  void SelfAssemblingRecoFrame::Disassemble(){
77  m_Nvisible = 0;
78  m_Ndecay = 0;
79 
80  // replace frames with unassembled ones
81  const LabRecoFrame& lab_frame = static_cast<const LabRecoFrame&>(GetLabFrame());
82  lab_frame.RemoveTreeStates(m_VisibleStates);
83  RestFrameList ChildFrames = m_ChildFrames_UnAssembled;
85  m_ChildFrames_UnAssembled = ChildFrames;
86  ClearNewFrames();
87  AddChildFrames(m_ChildFrames_UnAssembled);
88 
89  if(!InitializeTreeRecursive()){
90  m_Log << LogWarning;
91  m_Log << "Problem with recursive tree after disassembly";
92  m_Log << LogEnd;
93  SetBody(false);
94  return;
95  } else
96  SetBody(true);
97 
98  if(!InitializeAnalysisRecursive()){
99  m_Log << LogWarning;
100  m_Log << "Problem connecting states after disassembly";
101  m_Log << LogEnd;
102  SetMind(false);
103  return;
104  } else
105  SetMind(true);
106 
107  SetSpirit(false);
108  m_IsAssembled = false;
109  }
110 
111  void SelfAssemblingRecoFrame::Assemble(){
112  if(m_IsAssembled) Disassemble();
113  if(!IsSoundMind()){
114  UnSoundMind(RF_FUNCTION);
115  return;
116  }
117 
118  // new Visible States
119  m_VisibleStates.Clear();
120  // new Frames associated with States
121  std::vector<RestFrame*> frames;
122  // States' four-vector
123  std::vector<TLorentzVector> Ps;
124 
125  // clear unassembled lists
126  m_ChildFrames_UnAssembled.Clear();
127  m_ChildFrames_UnAssembled += GetChildFrames();
128 
129  int N = GetNChildren();
130  for(int i = 0; i < N; i++){
131  ReconstructionFrame& frame = GetChildFrame(i);
132  bool expand_frame = false;
133  if(GetChildStates(frame).GetN() == 1 && frame.IsVisibleFrame())
134  if(GetChildStates(frame)[0].IsCombinatoricState()){
135  expand_frame = true;
136  VisibleStateList const& elements =
137  static_cast<CombinatoricState&>(GetChildStates(frame)[0]).GetElements();
138  int Nelement = elements.GetN();
139  for(int e = 0; e < Nelement; e++){
140  VisibleState& element = elements[e];
141  VisibleRecoFrame& new_frame =
142  GetNewVisibleFrame(frame.GetName(),frame.GetTitle());
143  new_frame.SetCharge(element.GetCharge());
144  element.AddFrame(new_frame);
145  frames.push_back(&new_frame);
146  TLorentzVector V = element.GetFourVector();
147  if(V.M() < 0.) V.SetVectM(V.Vect(),0.);
148  Ps.push_back(V);
149  m_VisibleStates += element;
150  }
151  if(Nelement < 1){
152  expand_frame = false;
153  }
154  }
155  if(!expand_frame){
156  TLorentzVector V = GetChildStates(frame).GetFourVector();
157  if(V.M() < 0.) V.SetVectM(V.Vect(),0.);
158  Ps.push_back(V);
159  frames.push_back(&frame);
160  }
161  }
162 
163  RestFrameList ChildFrames = m_ChildFrames_UnAssembled;
165  m_ChildFrames_UnAssembled = ChildFrames;
166  AssembleRecursive(*this, frames, Ps);
167  if(!InitializeTreeRecursive()){
168  m_Log << LogWarning;
169  m_Log << "Problem with recursive tree after assembly";
170  m_Log << LogEnd;
171  SetBody(false);
172  return;
173  }
174 
175  SetMind(true);
176  const LabRecoFrame& lab_frame = static_cast<const LabRecoFrame&>(GetLabFrame());
177  lab_frame.AddTreeStates(m_VisibleStates);
178  if(!InitializeAnalysisRecursive()){
179  m_Log << LogWarning;
180  m_Log << "Problem connecting states after assembly";
181  m_Log << LogEnd;
182  SetMind(false);
183  return;
184  }
185 
186  m_IsAssembled = true;
187  }
188 
189  void SelfAssemblingRecoFrame::AssembleRecursive(RestFrame& frame,
190  std::vector<RestFrame*>& frames,
191  std::vector<TLorentzVector>& Ps){
192  int Ninput = frames.size();
193  if(Ninput <= 2){
194  for(int i = 0; i < Ninput; i++) frame.AddChildFrame(*frames[i]);
195  return;
196  }
197 
198  TLorentzVector TOT(0.,0.,0.,0.);
199  for(int i = 0; i < Ninput; i++) TOT += Ps[i];
200  TVector3 boost = TOT.BoostVector();
201  if(boost.Mag() > 1.)
202  boost.SetMagThetaPhi(0.,boost.Theta(),boost.Phi());
203  for(int i = 0; i < Ninput; i++){
204  if(Ps[i].M() < 0.)
205  Ps[i].SetVectM(Ps[i].Vect(),Ps[i].M());
206  Ps[i].Boost(-boost);
207  }
208 
209  int ip_max[2];
210  int jp_max[2];
211  for(int i = 0; i < 2; i++) ip_max[i] = -1;
212  for(int i = 0; i < 2; i++) jp_max[i] = -1;
213  double val_max = -1.;
214  // Loop over all 2-jet seed probes
215  int ip[2], jp[2];
216  for(ip[0] = 0; ip[0] < Ninput-1; ip[0]++){
217  for(ip[1] = ip[0]+1; ip[1] < Ninput; ip[1]++){
218  TVector3 nRef = Ps[ip[0]].Vect().Cross(Ps[ip[1]].Vect());
219  int Nhem[2];
220  TLorentzVector hem[2];
221  for(int i = 0; i < 2; i++){
222  Nhem[i] = 0;
223  hem[i].SetPxPyPzE(0.,0.,0.,0.);
224  }
225  // Loop over all jets
226  for(int i = 0; i < Ninput; i++){
227  if((i == ip[0]) || (i == ip[1])) continue;
228  int ihem = int(Ps[i].Vect().Dot(nRef) > 0.);
229  Nhem[ihem]++;
230  hem[ihem] += Ps[i];
231  }
232  // assign 2 probes
233  for(jp[0] = 0; jp[0] < 2; jp[0]++){
234  for(jp[1] = 0; jp[1] < 2; jp[1]++){
235  if(jp[0] == jp[1] && Nhem[(jp[0]+1)%2] == 0) continue;
236  TLorentzVector hem_probes[2];
237  for(int i = 0; i < 2; i++) hem_probes[i] = hem[i];
238  for(int i = 0; i < 2; i++) hem_probes[jp[i]] += Ps[ip[i]];
239  double val = hem_probes[0].P() + hem_probes[1].P();
240  if(val > val_max){
241  val_max = val;
242  for(int i = 0; i < 2; i++) ip_max[i] = ip[i];
243  for(int i = 0; i < 2; i++) jp_max[i] = jp[i];
244  }
245  }
246  }
247  }
248  }
249 
250  std::vector<RestFrame*> child_frames[2];
251  std::vector<TLorentzVector> child_Ps[2];
252  TLorentzVector hem[2];
253  for(int i = 0; i < 2; i++){
254  hem[i].SetPxPyPzE(0.,0.,0.,0.);
255  }
256 
257  for(int i = 0; i < 2; i++){
258  child_frames[jp_max[i]].push_back(frames[ip_max[i]]);
259  child_Ps[jp_max[i]].push_back(Ps[ip_max[i]]);
260  hem[jp_max[i]] += Ps[ip_max[i]];
261  }
262 
263  TVector3 nRef = Ps[ip_max[0]].Vect().Cross(Ps[ip_max[1]].Vect());
264  for(int i = 0; i < Ninput; i++){
265  if((i == ip_max[0]) || (i == ip_max[1])) continue;
266  int ihem = int(Ps[i].Vect().Dot(nRef) > 0.);
267  child_frames[ihem].push_back(frames[i]);
268  child_Ps[ihem].push_back(Ps[i]);
269  hem[ihem] += Ps[i];
270  }
271 
272  int flip = int(hem[1].M() > hem[0].M());
273  for(int i = 0; i < 2; i++){
274  int j = (i+flip)%2;
275  if(child_frames[j].size() == 1){
276  frame.AddChildFrame(*child_frames[j][0]);
277  } else {
278  RestFrame& new_frame = GetNewDecayFrame(GetName(),GetTitle());
279  frame.AddChildFrame(new_frame);
280  AssembleRecursive(new_frame, child_frames[j], child_Ps[j]);
281  }
282  }
283  }
284 
285  bool SelfAssemblingRecoFrame::ReconstructFrame(){
286  if(m_NewEvent){
287  m_NewEvent = false;
288  if(m_IsAssembled) Disassemble();
289  if(!AnalyzeEventRecursive()){
290  m_Log << LogWarning;
291  m_Log << "Unable to recursively analyze event with ";
292  m_Log << "disassembled SelfAssemblingRecoFrame" << LogEnd;
293  return SetSpirit(false);
294  }
295  Assemble();
296  }
297 
298  return ReconstructionFrame::ReconstructFrame();
299  }
300 
301  void SelfAssemblingRecoFrame::ClearNewFrames(){
302  int N = m_DecayFrames.GetN();
303  for(int i = 0; i < N; i++) m_DecayFrames[i].Clear();
304  N = m_VisibleFrames.GetN();
305  for(int i = 0; i < N; i++) m_VisibleFrames[i].Clear();
306  }
307 
308  DecayRecoFrame& SelfAssemblingRecoFrame::GetNewDecayFrame(const std::string& sname,
309  const std::string& stitle){
310  if(m_Ndecay < m_DecayFrames.GetN()){
311  m_Ndecay++;
312  return m_DecayFrames.Get(m_Ndecay-1);
313  }
314  char strn[10];
315  sprintf(strn,"%d",m_Ndecay+1);
316  std::string name = sname+"_"+std::string(strn);
317  std::string title = "#left("+stitle+"#right)_{"+std::string(strn)+"}";
318  DecayRecoFrame* framePtr = new DecayRecoFrame(name, title);
319 
320  m_DecayFrames.Add(*framePtr);
321  AddDependent(framePtr);
322  m_Ndecay++;
323  return *framePtr;
324  }
325 
326  VisibleRecoFrame& SelfAssemblingRecoFrame::GetNewVisibleFrame(const std::string& sname,
327  const std::string& stitle){
328  if(m_Nvisible < m_VisibleFrames.GetN()){
329  m_Nvisible++;
330  return m_VisibleFrames.Get(m_Nvisible-1);
331  }
332  char strn[10];
333  sprintf(strn,"%d",m_Nvisible+1);
334  std::string name = sname+"_"+std::string(strn);
335  std::string title = "#left("+stitle+"#right)_{"+std::string(strn)+"}";
336  VisibleRecoFrame* framePtr = new VisibleRecoFrame(name, title);
337 
338  m_VisibleFrames.Add(*framePtr);
339  AddDependent(framePtr);
340  m_Nvisible++;
341  return *framePtr;
342  }
343 
344  RestFrame const& SelfAssemblingRecoFrame::GetFrame(const RFKey& key) const {
345  if(!m_IsAssembled)
346  return RestFrame::Empty();
347 
348  int N = GetNChildren();
349  for(int i = 0; i < N; i++){
350  if(GetChildStates(i).Contains(key))
351  return GetChildStates(i).Get(key).GetListFrames()[0];
352  }
353 
354  return RestFrame::Empty();
355  }
356 
357 }
virtual void Clear()
Clears ReconstructionFrame of all connections to other objects.
virtual void Clear()
Clears ReconstructionFrame of all connections to other objects.
virtual void RemoveChildFrame(RestFrame &frame)
Remove a child of this frame.
abstract base class for all Frame objects
void RemoveChildFrame(RestFrame &frame)
Remove a child of this frame.
virtual ReconstructionFrame & GetChildFrame(int i=0) const
Get the frame of the i th child.
void RemoveChildFrames()
Remove all the children of this frame.