LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
RestFrame.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 "RestFrames/RestFrame.hh"
32 
33 namespace RestFrames {
34 
35  int RestFrame::m_class_key = 0;
36 
37  RestFrame::RestFrame() : RFBase() {
38  m_Type = kVanillaFrame;
39  m_Log.SetSource("RestFrame "+GetName());
40  }
41 
42  RestFrame::~RestFrame(){
43  Clear();
44  }
45 
46  RestFrame::RestFrame(const std::string& sname, const std::string& stitle)
47  : RFBase(sname, stitle, RestFrame::m_class_key++)
48  {
49  m_Log.SetSource("RestFrame "+GetName());
50  m_Type = kVanillaFrame;
51  m_ParentFramePtr = nullptr;
52  m_ParentBoost.SetXYZ(0.,0.,0.);
53  }
54 
55  void RestFrame::Clear(){
56  SetParentFrame();
57  RemoveChildFrames();
58  }
59 
60  RestFrame& RestFrame::Empty(){
61  return ReconstructionFrame::Empty();
62  }
63 
64  ConstRestFrameList const& RestFrame::EmptyList(){
65  return m_EmptyList;
66  }
67 
68  TVector3 RestFrame::m_Axis = TVector3(0.,0.,1.);
69 
70  void RestFrame::SetAxis(const TVector3& axis){
71  RestFrame::m_Axis = axis;
72  }
73 
74  TVector3 const& RestFrame::GetAxis(){
75  return RestFrame::m_Axis;
76  }
77 
78  RestFrameList RestFrame::operator + (RestFrame& frame){
79  RestFrameList list;
80  list.Add(frame);
81  list.Add(*this);
82  return list;
83  }
84 
85  RestFrameList RestFrame::operator + (const RestFrameList& frames){
86  RestFrameList list = frames;
87  list.Add(*this);
88  return list;
89  }
90 
91  FrameType RestFrame::GetType() const {
92  return m_Type;
93  }
94 
95  bool RestFrame::IsVisibleFrame() const {
96  return m_Type == kVisibleFrame;
97  }
98 
99  bool RestFrame::IsInvisibleFrame() const {
100  return m_Type == kInvisibleFrame;
101  }
102 
103  bool RestFrame::IsDecayFrame() const {
104  return m_Type == kDecayFrame;
105  }
106 
107  bool RestFrame::IsLabFrame() const {
108  return m_Type == kLabFrame;
109  }
110 
111  bool RestFrame::IsRecoFrame() const {
112  return m_Ana == kRecoFrame;
113  }
114 
115  bool RestFrame::IsGenFrame() const {
116  return m_Ana == kGenFrame;
117  }
118 
119  std::string RestFrame::PrintString(LogType type) const {
120  std::string output = RFBase::PrintString(type);
121  if(IsLabFrame())
122  output += " Frame Type: Lab \n";
123  if(IsDecayFrame())
124  output += " Frame Type: Decay \n";
125  if(IsVisibleFrame())
126  output += " Frame Type: Visible \n";
127  if(IsInvisibleFrame())
128  output += " Frame Type: Invisible \n";
129  if(IsGenFrame())
130  output += " Ana Type: Generator \n";
131  if(IsRecoFrame())
132  output += " Ana Type: Reconstruction \n";
133  return output;
134  }
135 
136  bool RestFrame::IsSoundBody() const {
137  if(RFBase::IsSoundBody()) return true;
138  int Nchild = GetNChildren();
139  for(int i = 0; i < Nchild; i++)
140  if(!GetChildFrame(i)){
141  m_Log << LogWarning << "Empty child frame:";
142  m_Log << Log(GetChildFrame(i)) << LogEnd;
143  return SetBody(false);
144  }
145  return SetBody(true);
146  }
147 
148  bool RestFrame::InitializeTreeRecursive() {
149  if(!IsSoundBody()) return false;
150  int Nchild = GetNChildren();
151  for(int i = 0; i < Nchild; i++){
152  if(!GetChildFrame(i).InitializeTreeRecursive()){
153  m_Log << LogWarning;
154  m_Log << "Problem with recursive tree structure from frame: ";
155  m_Log << Log(GetChildFrame(i)) << LogEnd;
156  return SetBody(false);
157  }
158  }
159  return true;
160  }
161 
162  void RestFrame::RemoveChildFrame(RestFrame& frame){
163  SetBody(false);
164  bool contains = m_ChildFrames.Contains(frame);
165  m_ChildBoosts.erase(&frame);
166  m_ChildFrames.Remove(frame);
167  if(contains)
168  frame.SetParentFrame();
169  }
170 
171  void RestFrame::RemoveChildFrames(){
172  SetBody(false);
173  while(GetNChildren() > 0)
174  RemoveChildFrame(m_ChildFrames[0]);
175  m_ChildFrames.Clear();
176  m_ChildBoosts.clear();
177  }
178 
179  void RestFrame::SetParentFrame(RestFrame& frame){
180  if(IsEmpty()) return;
181 
182  SetBody(false);
183 
184  RestFrame* prevPtr = m_ParentFramePtr;
185  if(m_ParentFramePtr)
186  if(*m_ParentFramePtr != frame){
187  m_ParentFramePtr = nullptr;
188  prevPtr->RemoveChildFrame(*this);
189  }
190  if(!frame)
191  m_ParentFramePtr = nullptr;
192  else
193  m_ParentFramePtr = &frame;
194  }
195 
196  void RestFrame::AddChildFrame(RestFrame& frame){
197  if(IsEmpty()) return;
198 
199  SetBody(false);
200 
201  if(!frame){
202  m_Log << LogWarning;
203  m_Log << "Cannot add empty frame as child.";
204  m_Log << LogEnd;
205  return;
206  }
207  if(frame.IsLabFrame()){
208  m_Log << LogWarning;
209  m_Log << "Cannot add LabFrame frame as child:";
210  m_Log << Log(frame) << LogEnd;
211  return;
212  }
213  if(!m_ChildFrames.Add(frame)){
214  m_Log << LogWarning;
215  m_Log << "Unable to add child frame:";
216  m_Log << Log(frame) << LogEnd;
217  return;
218  }
219  frame.SetParentFrame(*this);
220  m_ChildBoosts[&frame] = m_Empty3Vector;
221  }
222 
223  void RestFrame::AddChildFrames(const RestFrameList& frames){
224  int N = frames.GetN();
225  for(int i = 0; i < N; i++)
226  AddChildFrame(frames[i]);
227  }
228 
229  int RestFrame::GetNChildren() const {
230  return m_ChildFrames.GetN();
231  }
232 
233  RestFrame& RestFrame::GetChildFrame(int i) const {
234  int Nchild = GetNChildren();
235  if(i >= Nchild || i < 0){
236  m_Log << LogWarning;
237  m_Log << "Cannot GetChildFrame(" << i << "). ";
238  m_Log << "No " << i << "th child" << LogEnd;
239  }
240  return m_ChildFrames[i];
241  }
242 
243  RestFrame const& RestFrame::GetParentFrame() const {
244  if(m_ParentFramePtr)
245  return *m_ParentFramePtr;
246  else
247  return RestFrame::Empty();
248  }
249 
250  RestFrameList const& RestFrame::GetChildFrames() const {
251  return m_ChildFrames;
252  }
253 
254  RestFrame const& RestFrame::GetLabFrame() const {
255  if(IsLabFrame())
256  return *this;
257 
258  if(!GetParentFrame()){
259  m_Log << LogWarning;
260  m_Log << "Unable to find LabFrame above this frame. ";
261  m_Log << "No parent frame set" << LogEnd;
262  return RestFrame::Empty();
263  }
264  return m_ParentFramePtr->GetLabFrame();
265  }
266 
267  RestFrame const& RestFrame::GetSiblingFrame() const {
268  if(!IsSoundBody())
269  return RestFrame::Empty();
270 
271  if(IsLabFrame())
272  return RestFrame::Empty();
273 
274  if(!GetParentFrame())
275  return RestFrame::Empty();
276 
277  int Nsib = GetParentFrame().GetNChildren();
278  for(int s = 0; s < Nsib; s++){
279  if(IsSame(m_ParentFramePtr->GetChildFrame(s)))
280  continue;
281  return m_ParentFramePtr->GetChildFrame(s);
282  }
283  return RestFrame::Empty();
284  }
285 
286  int RestFrame::GetNDescendants() const {
287  if(!IsSoundBody()) return 0.;
288 
289  int Nchild = GetNChildren();
290  if(Nchild == 0) return 1;
291  int Nd = 0;
292  for(int i = 0; i < Nchild; i++){
293  Nd += GetChildFrame(i).GetNDescendants();
294  }
295  return Nd;
296  }
297 
298  void RestFrame::SetChildBoostVector(RestFrame& frame, const TVector3& boost) {
299  if(!m_ChildFrames.Contains(frame)){
300  m_Log << LogWarning;
301  m_Log << "Unable to set child's boost vector. ";
302  m_Log << "Frame is not among children:";
303  m_Log << Log(frame) << LogEnd;
304  return;
305  }
306  m_ChildBoosts[&frame] = boost;
307  frame.SetParentBoostVector(-1.*boost);
308  }
309 
310  void RestFrame::SetParentBoostVector(const TVector3& boost) {
311  if(!GetParentFrame()){
312  m_Log << LogWarning;
313  m_Log << "Unable to set parent boost vector. ";
314  m_Log << "No parent frame set.";
315  m_Log << LogEnd;
316  return;
317  }
318  m_ParentBoost = boost;
319  }
320 
321  TVector3 const& RestFrame::GetChildBoostVector(RestFrame& frame) const {
322  if(!m_ChildFrames.Contains(frame)){
323  m_Log << LogWarning;
324  m_Log << "Unable to get child's boost vector. ";
325  m_Log << "Frame is not among children:";
326  m_Log << Log(frame) << LogEnd;
327  return m_Empty3Vector;
328  }
329  return m_ChildBoosts[&frame];
330  }
331 
332  TVector3 const& RestFrame::GetParentBoostVector() const {
333  return m_ParentBoost;
334  }
335 
336  RestFrameList RestFrame::GetListFrames(FrameType type) const {
337  RestFrameList frames;
338  FillListFramesRecursive(frames,type);
339  return frames;
340  }
341 
342  RestFrameList RestFrame::GetListVisibleFrames() const {
343  return GetListFrames(kVisibleFrame);
344  }
345 
346  RestFrameList RestFrame::GetListInvisibleFrames() const {
347  return GetListFrames(kInvisibleFrame);
348  }
349 
350  void RestFrame::FillListFramesRecursive(RestFrameList& frames, FrameType type) const {
351  if(frames.Contains(*this)) return;
352  if(type == GetType() || type == kLabFrame) frames.Add((RestFrame&)(*m_This));
353  int Nchild = GetNChildren();
354  for(int i = 0; i < Nchild; i++)
355  GetChildFrame(i).FillListFramesRecursive(frames, type);
356  }
357 
358  bool RestFrame::IsCircularTree(std::vector<RFKey>& keys) const {
359  int Nkey = keys.size();
360  for(int i = 0; i < Nkey; i++){
361  if(keys[i] == GetKey()){
362  m_Log << LogWarning;
363  m_Log << "This RestFrame appears more than once in the tree:";
364  m_Log << Log(*this) << LogEnd;
365  return true;
366  }
367  }
368  keys.push_back(GetKey());
369  int Nchild = GetNChildren();
370  for(int i = 0; i < Nchild; i++)
371  if(GetChildFrame(i).IsCircularTree(keys))
372  return true;
373 
374  return false;
375  }
376 
377  bool RestFrame::FindPathToFrame(const RestFrame& dest_frame,
378  const RestFrame& prev_frame,
379  std::vector<const TVector3*>& boosts) const {
380  if(IsSame(dest_frame)) return true;
381 
382  std::vector<const RestFrame*> try_frames;
383  std::vector<const TVector3*> try_boosts;
384 
385  if(!GetParentFrame().IsEmpty()){
386  try_frames.push_back(&GetParentFrame());
387  try_boosts.push_back(&GetParentBoostVector());
388  }
389  int Nchild = GetNChildren();
390  for(int i = 0; i < Nchild; i++){
391  try_frames.push_back(&GetChildFrame(i));
392  try_boosts.push_back(&GetChildBoostVector(GetChildFrame(i)));
393  }
394 
395  int Ntry = try_frames.size();
396  for(int i = 0; i < Ntry; i++){
397  const RestFrame* nextPtr = try_frames[i];
398  if(nextPtr->IsSame(prev_frame)) continue;
399  boosts.push_back(try_boosts[i]);
400  if(nextPtr->FindPathToFrame(dest_frame,*this,boosts))
401  return true;
402  boosts.pop_back();
403  }
404  return false;
405  }
406 
407  void RestFrame::SetFourVector(const TLorentzVector& V,
408  const RestFrame& frame){
409  if(!IsSoundBody()){
410  UnSoundBody(RF_FUNCTION);
411  return;
412  }
413  if(!frame)
414  m_ProdFramePtr = &GetLabFrame();
415  else
416  m_ProdFramePtr = &frame;
417 
418  m_P = V;
419  }
420 
422  // User Analysis functions
424 
425  RFCharge RestFrame::GetCharge() const {
426  RFCharge charge;
427 
428  if(!IsSoundBody()){
429  UnSoundBody(RF_FUNCTION);
430  return charge;
431  }
432 
433  int Nchild = GetNChildren();
434  if(Nchild == 0) return charge;
435  for(int i = 0; i < Nchild; i++){
436  charge += GetChildFrame(i).GetCharge();
437  }
438  return charge;
439  }
440 
441  double RestFrame::GetMass() const {
442  if(!IsSoundSpirit()){
443  UnSoundSpirit(RF_FUNCTION);
444  return 0.;
445  }
446  return m_P.M();
447  }
448 
449  TLorentzVector RestFrame::GetFourVector(const RestFrame& frame) const {
450  if(!IsSoundSpirit()){
451  UnSoundSpirit(RF_FUNCTION);
452  return m_Empty4Vector;
453  }
454 
455  if(!frame)
456  return GetFourVector(GetLabFrame());
457 
458  if(!GetProductionFrame()){
459  m_Log << LogWarning;
460  m_Log << "Unable to get four vector. ";
461  m_Log << "Production frame is not defined." << LogEnd;
462  return m_Empty4Vector;
463  }
464 
465  TLorentzVector V = m_P;
466  if(V.E() <= 0.)
467  return m_Empty4Vector;
468  if(frame == GetProductionFrame()) return V;
469 
470  std::vector<const TVector3*> boosts;
471  if(!GetProductionFrame().
472  FindPathToFrame(frame, RestFrame::Empty(), boosts)){
473  m_Log << LogWarning;
474  m_Log << "Unable to get four vector. ";
475  m_Log << "Cannot find a path to frame " << frame.GetName();
476  m_Log << " from frame " << GetProductionFrame().GetName() << LogEnd;
477  return m_Empty4Vector;
478  }
479 
480  int Nboost = boosts.size();
481  for(int i = 0; i < Nboost; i++)
482  V.Boost(-1.*(*boosts[i]));
483  return V;
484  }
485 
486  TLorentzVector RestFrame::GetVisibleFourVector(const RestFrame& frame) const {
487  if(!IsSoundSpirit()){
488  UnSoundSpirit(RF_FUNCTION);
489  return m_Empty4Vector;
490  }
491 
492  if(!frame){
493  if(!GetLabFrame())
494  return m_Empty4Vector;
495  else
496  return GetVisibleFourVector(GetLabFrame());
497  }
498 
499  if(IsVisibleFrame())
500  return GetFourVector(frame);
501 
502  TLorentzVector V(0.,0.,0.,0.);
503  int Nc = GetNChildren();
504  for(int c = 0; c < Nc; c++){
505  RestFrameList frames = GetChildFrame(c).GetListVisibleFrames();
506  int Nf = frames.GetN();
507  for(int f = 0; f < Nf; f++)
508  V += frames[f].GetFourVector(frame);
509  }
510  return V;
511  }
512 
513  TLorentzVector RestFrame::GetInvisibleFourVector(const RestFrame& frame) const {
514  if(!IsSoundSpirit()){
515  UnSoundSpirit(RF_FUNCTION);
516  return m_Empty4Vector;
517  }
518 
519  if(!frame){
520  if(!GetLabFrame())
521  return m_Empty4Vector;
522  else
523  return GetInvisibleFourVector(GetLabFrame());
524  }
525 
526  if(IsInvisibleFrame())
527  return GetFourVector(frame);
528 
529  TLorentzVector V(0.,0.,0.,0.);
530  int Nc = GetNChildren();
531  for(int c = 0; c < Nc; c++){
532  RestFrameList frames = GetChildFrame(c).GetListInvisibleFrames();
533  int Nf = frames.GetN();
534  for(int f = 0; f < Nf; f++)
535  V += frames[f].GetFourVector(frame);
536  }
537  return V;
538  }
539 
540  double RestFrame::GetEnergy(const RestFrame& frame) const {
541  if(!IsSoundSpirit()){
542  UnSoundSpirit(RF_FUNCTION);
543  return 0.;
544  }
545  return GetFourVector(frame).E();
546  }
547 
548  double RestFrame::GetMomentum(const RestFrame& frame) const {
549  if(!IsSoundSpirit()){
550  UnSoundSpirit(RF_FUNCTION);
551  return 0.;
552  }
553  return GetFourVector(frame).P();
554  }
555 
556  double RestFrame::GetTransverseMomentum(const RestFrame& frame,
557  const TVector3& axis,
558  const RestFrame& axis_frame) const {
559  if(!IsSoundSpirit()){
560  UnSoundSpirit(RF_FUNCTION);
561  return 0.;
562  }
563 
564  if(frame == axis_frame){
565  TVector3 V = GetFourVector(frame).Vect();
566  return (V-V.Dot(axis.Unit())*axis.Unit()).Mag();
567  }
568 
569  TLorentzVector Pthis = GetFourVector(axis_frame);
570 
571  TLorentzVector Pframe;
572  if(!frame || frame.IsLabFrame()){
573  Pframe = axis_frame.GetFourVector(frame);
574  Pframe.SetVectM(-Pframe.Vect(),Pframe.M());
575  } else {
576  Pframe = frame.GetFourVector(axis_frame);
577  }
578 
579  TVector3 boost_par = Pframe.BoostVector();
580 
581  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
582  Pthis.Boost(-boost_par);
583  Pframe.Boost(-boost_par);
584  TVector3 boost_perp = Pframe.BoostVector();
585  Pthis.Boost(-boost_perp);
586 
587  TVector3 V = Pthis.Vect();
588  return (V-V.Dot(axis.Unit())*axis.Unit()).Mag();
589  }
590 
591  TLorentzVector RestFrame::GetFourVector(const TLorentzVector& P,
592  const RestFrame& def_frame) const {
593  if(!IsSoundSpirit()){
594  UnSoundSpirit(RF_FUNCTION);
595  return m_Empty4Vector;
596  }
597 
598  if(IsSame(def_frame) || (!def_frame && IsLabFrame()))
599  return P;
600 
601  TLorentzVector Pret = P;
602 
603  std::vector<const TVector3*> boosts;
604  if(!def_frame){
605  if(!GetLabFrame().
606  FindPathToFrame(*this, RestFrame::Empty(), boosts)){
607  m_Log << LogWarning;
608  m_Log << "Unable to get four vector. ";
609  m_Log << "Cannot find a path to frame " << GetName();
610  m_Log << " from frame " << GetLabFrame().GetName() << LogEnd;
611  return m_Empty4Vector;
612  }
613  } else {
614  if(!def_frame.
615  FindPathToFrame(*this, RestFrame::Empty(), boosts)){
616  m_Log << LogWarning;
617  m_Log << "Unable to get four vector. ";
618  m_Log << "Cannot find a path to frame " << GetName();
619  m_Log << " from frame " << GetLabFrame().GetName() << LogEnd;
620  return m_Empty4Vector;
621  }
622  }
623  int Nboost = boosts.size();
624  for(int i = 0; i < Nboost; i++)
625  Pret.Boost(-1.*(*boosts[i]));
626  return Pret;
627  }
628 
629  double RestFrame::GetTransverseMomentum(const TLorentzVector& P,
630  const TVector3& axis,
631  const RestFrame& axis_frame) const {
632  if(!IsSoundSpirit()){
633  UnSoundSpirit(RF_FUNCTION);
634  return 0.;
635  }
636 
637  if(IsLabFrame() && (!axis_frame || axis_frame.IsLabFrame())){
638  TVector3 V = P.Vect();
639  return (V-V.Dot(axis.Unit())*axis.Unit()).Mag();
640  }
641 
642  TLorentzVector Pret = P;
643 
644  // move P to axis_frame
645  if(!axis_frame.IsLabFrame() && !axis_frame.IsEmpty()){
646  std::vector<const TVector3*> boosts;
647  if(!GetLabFrame().
648  FindPathToFrame(axis_frame, RestFrame::Empty(), boosts)){
649  m_Log << LogWarning;
650  m_Log << "Unable to get four vector. ";
651  m_Log << "Cannot find a path to frame " << axis_frame.GetName();
652  m_Log << " from frame " << GetLabFrame().GetName() << LogEnd;
653  return 0.;
654  }
655  int Nboost = boosts.size();
656  for(int i = 0; i < Nboost; i++)
657  Pret.Boost(-1.*(*boosts[i]));
658  }
659 
660  if(*this == axis_frame){
661  TVector3 V = Pret.Vect();
662  return (V-V.Dot(axis.Unit())*axis.Unit()).Mag();
663  }
664 
665  TLorentzVector Pthis;
666  if(IsLabFrame()){
667  Pthis = axis_frame.GetFourVector();
668  Pthis.SetVectM(-Pthis.Vect(),Pthis.M());
669  } else {
670  Pthis = GetFourVector(axis_frame);
671  }
672 
673  TVector3 boost_par = Pthis.BoostVector();
674 
675  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
676  Pret.Boost(-boost_par);
677  Pthis.Boost(-boost_par);
678  TVector3 boost_perp = Pthis.BoostVector();
679  Pret.Boost(-boost_perp);
680 
681  TVector3 V = Pret.Vect();
682  return (V-V.Dot(axis.Unit())*axis.Unit()).Mag();
683  }
684 
685 
686  int RestFrame::GetFrameDepth(const RestFrame& frame) const {
687  if(!IsSoundSpirit()){
688  UnSoundSpirit(RF_FUNCTION);
689  return 0.;
690  }
691 
692  if(!frame) return -1;
693  if(IsSame(frame)) return 0.;
694  int Nchild = GetNChildren();
695  for(int i = 0; i < Nchild; i++){
696  int depth = GetChildFrame(i).GetFrameDepth(frame);
697  if(depth >= 0) return depth+1;
698  }
699  return -1;
700  }
701 
702  RestFrame const& RestFrame::GetFrameAtDepth(int depth, const RestFrame& frame) const {
703  if(!IsSoundSpirit()){
704  UnSoundSpirit(RF_FUNCTION);
705  return RestFrame::Empty();
706  }
707 
708  if(!frame || depth < 1)
709  return RestFrame::Empty();
710 
711  int N = GetNChildren();
712  for(int i = 0; i < N; i++){
713  RestFrame& child = GetChildFrame(i);
714  if(child.GetListFrames().Contains(frame)){
715  if(depth == 1) return child;
716  else return child.GetFrameAtDepth(depth-1,frame);
717  }
718  }
719  return RestFrame::Empty();
720  }
721 
722  double RestFrame::GetVisibleShape() const {
723  if(!IsSoundSpirit()){
724  UnSoundSpirit(RF_FUNCTION);
725  return 0.;
726  }
727  double Psum = 0.;
728  double Xsum = 0.;
729  int N = GetNChildren();
730  for(int i = 0; i < N; i++)
731  Psum += GetChildFrame(i).GetVisibleFourVector(*this).P();
732  for(int i = 0; i < N-1; i++){
733  TVector3 P1 = GetChildFrame(i).GetVisibleFourVector(*this).Vect();
734  for(int j = i+1; j < N; j++){
735  TVector3 P2 = GetChildFrame(j).GetVisibleFourVector(*this).Vect();
736  Xsum += (P1.Mag()+P2.Mag())*(P1.Mag()+P2.Mag())-(P1-P2).Mag2();
737  }
738  }
739  if(Psum > 0.)
740  return sqrt(Xsum)/Psum;
741  else
742  return 0.;
743  }
744 
745  double RestFrame::GetSumVisibleMomentum() const {
746  if(!IsSoundSpirit()){
747  UnSoundSpirit(RF_FUNCTION);
748  return 0.;
749  }
750 
751  double ret = 0.;
752  int N = GetNChildren();
753  for(int i = 0; i < N; i++)
754  ret += GetChildFrame(i).GetVisibleFourVector(*this).P();
755 
756  return ret;
757  }
758 
759  double RestFrame::GetSumInvisibleMomentum() const {
760  if(!IsSoundSpirit()){
761  UnSoundSpirit(RF_FUNCTION);
762  return 0.;
763  }
764 
765  double ret = 0.;
766  int N = GetNChildren();
767  for(int i = 0; i < N; i++)
768  ret += GetChildFrame(i).GetVisibleFourVector(*this).P();
769 
770  return ret;
771  }
772 
773  RestFrame const& RestFrame::GetProductionFrame() const {
774  if(!IsSoundSpirit()){
775  UnSoundSpirit(RF_FUNCTION);
776  return RestFrame::Empty();
777  }
778 
779  if(m_ProdFramePtr)
780  return *m_ProdFramePtr;
781  else
782  return RestFrame::Empty();
783  }
784 
785  TVector3 RestFrame::GetBoostInParentFrame() const {
786  TVector3 V(0.,0.,0.);
787 
788  if(!IsSoundSpirit()){
789  UnSoundSpirit(RF_FUNCTION);
790  return V;
791  }
792 
793  if(!GetParentFrame()) return V;
794  return -1.*GetParentBoostVector();
795  }
796 
797  double RestFrame::GetGammaInParentFrame() const {
798  if(!IsSoundSpirit()){
799  UnSoundSpirit(RF_FUNCTION);
800  return 0.;
801  }
802 
803  TVector3 vbeta = GetBoostInParentFrame();
804  double beta = std::min(1.,vbeta.Mag());
805  return 1./sqrt(1.-beta*beta);
806  }
807 
808  double RestFrame::GetCosDecayAngle(const RestFrame& frame) const {
809  if(!IsSoundSpirit()){
810  UnSoundSpirit(RF_FUNCTION);
811  return 0.;
812  }
813 
814  TVector3 V1 = GetParentBoostVector().Unit();
815  if(IsLabFrame())
816  V1 = RestFrame::GetAxis();
817  TVector3 V2;
818  if(!frame.IsEmpty())
819  V2 = frame.GetFourVector(*this).Vect().Unit();
820  else
821  if(GetNChildren() < 1)
822  return 0.;
823  else
824  V2 = GetChildFrame(0).GetFourVector(*this).Vect().Unit();
825 
826  return V1.Dot(V2);
827  }
828 
829  TVector3 RestFrame::GetDecayPlaneNormalVector(const RestFrame& frame) const {
830  TVector3 n = RestFrame::GetAxis();
831 
832  if(!IsSoundSpirit()){
833  UnSoundSpirit(RF_FUNCTION);
834  return n;
835  }
836 
837  if(!frame)
838  if(GetNChildren() < 1)
839  return n;
840 
841  TVector3 V1, V2;
842  if(!IsLabFrame()){
843  if(!frame)
844  V1 = GetChildFrame(0).GetFourVector(GetParentFrame()).Vect();
845  else
846  V1 = frame.GetFourVector(GetParentFrame()).Vect();
847  V2 = GetFourVector(GetParentFrame()).Vect();
848  } else {
849  if(!frame)
850  V1 = GetChildFrame(0).GetFourVector(*this).Vect().Unit();
851  else
852  V1 = frame.GetFourVector(*this).Vect().Unit();
853  V2 = n;
854  }
855  TVector3 ret = V1.Cross(V2);
856  if(ret.Mag() > 0)
857  return ret.Unit();
858 
859  std::vector<TVector3> tries;
860  tries.push_back(TVector3(1.,0.,0.));
861  tries.push_back(TVector3(0.,1.,0.));
862  tries.push_back(TVector3(0.,0.,1.));
863  for(int i = 0; i < 3; i++){
864  V2 = tries[i];
865  ret = V1.Cross(V2);
866  if(ret.Mag() > 0)
867  return ret.Unit();
868  }
869 
870  return n;
871  }
872 
873  double RestFrame::GetDeltaPhiDecayPlanes(const RestFrame& frame) const {
874  if(!IsSoundSpirit()){
875  UnSoundSpirit(RF_FUNCTION);
876  return 0.;
877  }
878 
879  if(!frame) return 0.;
880  if(GetNChildren() < 1) return 0.;
881 
882  TVector3 vNorm_frame = frame.GetDecayPlaneNormalVector();
883  TVector3 vNorm_this = GetDecayPlaneNormalVector();
884  double dphi = vNorm_this.Angle(vNorm_frame);
885 
886  if(frame.GetFourVector(*this).Vect().Cross(vNorm_frame).Dot(vNorm_this) < 0.){
887  dphi = TMath::Pi()*2. - dphi;
888  }
889 
890  return dphi;
891  }
892 
893  double RestFrame::GetDeltaPhiDecayAngle(const TVector3& axis, const RestFrame& frame) const {
894  if(!IsSoundSpirit()){
895  UnSoundSpirit(RF_FUNCTION);
896  return 0.;
897  }
898 
899  if(IsLabFrame())
900  return 0.;
901 
902  TLorentzVector Pthis = GetFourVector(frame);
903  TLorentzVector Pchild = GetChildFrame(0).GetFourVector(frame);
904 
905  TVector3 boost_par = Pthis.BoostVector();
906  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
907  Pthis.Boost(-boost_par);
908  Pchild.Boost(-boost_par);
909  TVector3 boost_perp = Pthis.BoostVector();
910  Pchild.Boost(-boost_perp);
911 
912  TVector3 V1 = boost_perp;
913  TVector3 V2 = Pchild.Vect();
914  V1 = V1 - V1.Dot(axis.Unit())*axis.Unit();
915  V2 = V2 - V2.Dot(axis.Unit())*axis.Unit();
916  return V1.Angle(V2);
917  }
918 
919  // Get angle between 'this' boost and visible children in plane
920  // perpendicular to 3-vector 'axis', where axis is defined
921  // in 'framePtr' (default gives lab frame).
922  double RestFrame::GetDeltaPhiBoostVisible(const TVector3& axis, const RestFrame& frame) const {
923  if(!IsSoundSpirit()){
924  UnSoundSpirit(RF_FUNCTION);
925  return 0.;
926  }
927 
928  TLorentzVector Pvis = GetVisibleFourVector(frame);
929  TLorentzVector Pthis = GetFourVector(frame);
930 
931  TVector3 boost_par = Pthis.BoostVector();
932  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
933  Pthis.Boost(-boost_par);
934  Pvis.Boost(-boost_par);
935  TVector3 boost_perp = Pthis.BoostVector();
936  Pvis.Boost(-boost_perp);
937 
938  TVector3 V = Pvis.Vect();
939  V = V - V.Dot(axis.Unit())*axis.Unit();
940 
941  return V.Angle(boost_perp);
942  }
943 
944  // Get angle between 'this' decay axis (defined by first child)
945  // and visible children in plane
946  // perpendicular to 3-vector 'axis', where axis is defined
947  // in 'framePtr' (default gives lab frame).
948  double RestFrame::GetDeltaPhiDecayVisible(const TVector3& axis, const RestFrame& frame) const {
949  if(!IsSoundSpirit()){
950  UnSoundSpirit(RF_FUNCTION);
951  return 0.;
952  }
953 
954  if(GetNChildren() < 1) return 0.;
955 
956  TLorentzVector Pvis = GetVisibleFourVector(frame);
957  TLorentzVector Pchild = GetChildFrame(0).GetFourVector(frame);
958  TLorentzVector Pthis = GetFourVector(frame);
959 
960  TVector3 boost_par = Pthis.BoostVector();
961  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
962  Pthis.Boost(-boost_par);
963  Pvis.Boost(-boost_par);
964  Pchild.Boost(-boost_par);
965  TVector3 boost_perp = Pthis.BoostVector();
966  Pvis.Boost(-boost_perp);
967  Pchild.Boost(-boost_perp);
968 
969  TVector3 Vv = Pvis.Vect();
970  Vv = Vv - Vv.Dot(axis.Unit())*axis.Unit();
971  TVector3 Vc = Pchild.Vect();
972  Vc = Vc - Vc.Dot(axis.Unit())*axis.Unit();
973 
974  return Vv.Angle(Vc);
975  }
976 
977  // Get angle between the visible portions of children 1 and 2
978  // in the plane perpendicular to 3-vector 'axis', where
979  // axis is defined in 'framePtr' (default gives lab frame).
980  double RestFrame::GetDeltaPhiVisible(const TVector3& axis, const RestFrame& frame) const {
981  if(!IsSoundSpirit()){
982  UnSoundSpirit(RF_FUNCTION);
983  return 0.;
984  }
985 
986  if(GetNChildren() != 2) return 0.;
987 
988  TLorentzVector Pthis = GetFourVector(frame);
989  TLorentzVector P1 = GetChildFrame(0).GetVisibleFourVector(frame);
990  TLorentzVector P2 = GetChildFrame(1).GetVisibleFourVector(frame);
991 
992  TVector3 boost_par = Pthis.BoostVector();
993  boost_par = boost_par.Dot(axis.Unit())*axis.Unit();
994  Pthis.Boost(-boost_par);
995  P1.Boost(-boost_par);
996  P2.Boost(-boost_par);
997  TVector3 boost_perp = Pthis.BoostVector();
998  P1.Boost(-boost_perp);
999  P2.Boost(-boost_perp);
1000 
1001  TVector3 V1 = P1.Vect();
1002  TVector3 V2 = P2.Vect();
1003  V1 = V1 - V1.Dot(axis.Unit())*axis.Unit();
1004  V2 = V2 - V2.Dot(axis.Unit())*axis.Unit();
1005 
1006  return V1.Angle(V2);
1007  }
1008 
1009  const ConstRestFrameList RestFrame::m_EmptyList;
1010 
1011 }
virtual std::string PrintString(LogType type) const
String of information associated with object.
Definition: RFBase.cc:146