LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
MinMassesSqInvJigsaw.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 <TMatrixD.h>
31 #include <TDecompSVD.h>
32 
33 #include "RestFrames/MinMassesSqInvJigsaw.hh"
35 
36 namespace RestFrames {
37 
38  MinMassesSqInvJigsaw::MinMassesSqInvJigsaw(const std::string& sname,
39  const std::string& stitle,
40  int N_vis_inv_pair) :
41  InvisibleJigsaw(sname, stitle, N_vis_inv_pair, N_vis_inv_pair),
42  m_Npair(N_vis_inv_pair) {}
43 
44  MinMassesSqInvJigsaw::MinMassesSqInvJigsaw()
45  : InvisibleJigsaw(), m_Npair(0) {}
46 
47  MinMassesSqInvJigsaw::~MinMassesSqInvJigsaw() {}
48 
49  MinMassesSqInvJigsaw& MinMassesSqInvJigsaw::Empty(){
50  return MinMassesSqInvJigsaw::m_Empty;
51  }
52 
53  double MinMassesSqInvJigsaw::GetMinimumMass() const {
54  if(!IsSoundMind())
55  return 0.;
56 
57  if(m_Npair <= 0)
58  return 0.;
59  if(m_Npair == 1)
60  return std::max(0.,GetChildState(0).GetMinimumMass());
61 
62  TLorentzVector PvisTOT(0.,0.,0.,0.);
63  m_Pvis.clear();
64  for(int i = 0; i < m_Npair; i++){
65  m_Pvis.push_back(GetDependancyStates(i+m_Npair).GetFourVector());
66  PvisTOT += m_Pvis[i];
67  }
68 
69  double ECM = 0.;
70  double Minv;
71  TVector3 Boost = PvisTOT.BoostVector();
72  for(int i = 0; i < m_Npair; i++){
73  m_Pvis[i].Boost(-Boost);
74  Minv = std::max(0.,GetChildState(i).GetMinimumMass());
75  ECM += sqrt(m_Pvis[i].P()*m_Pvis[i].P() + Minv*Minv);
76  }
77 
78  return ECM;
79  }
80 
81  bool MinMassesSqInvJigsaw::AnalyzeEvent(){
82  if(!IsSoundMind())
83  return SetSpirit(false);
84 
85  if(m_Npair <= 1)
86  return false;
87 
88  TLorentzVector INV = GetParentState().GetFourVector();
89  double Minv = INV.M();
90 
91  TLorentzVector VIS(0.,0.,0.,0.);
92  m_Pvis.clear();
93  for(int i = 0; i < m_Npair; i++){
94  m_Pvis.push_back(GetDependancyStates(i).GetFourVector());
95  VIS += m_Pvis[i];
96  }
97 
98  // go to INV+VIS CM frame
99  TVector3 BoostCM = (VIS+INV).BoostVector();
100  INV.Boost(-BoostCM);
101  VIS.Boost(-BoostCM);
102  for(int i = 0; i < m_Npair; i++)
103  m_Pvis[i].Boost(-BoostCM);
104 
105  // INV states defined in the VIS frame
106  TVector3 BoostVIS = VIS.BoostVector();
107  m_Pinv.clear();
108  m_Minv.clear();
109  for(int i = 0; i < m_Npair; i++){
110  m_Minv.push_back(std::max(0.,GetChildState(i).GetMinimumMass()));
111  m_Pinv.push_back(m_Pvis[i]);
112  m_Pinv[i].Boost(-BoostVIS);
113  m_Pinv[i].SetVectM(m_Pinv[i].Vect(), m_Minv[i]);
114  }
115 
116  // VIS states in INV frame
117  TVector3 BoostINV = INV.BoostVector();
118  for(int i = 0; i < m_Npair; i++)
119  m_Pvis[i].Boost(-BoostINV);
120 
121  if(m_Npair == 2){
122  TVector3 Vdiff = (m_Pvis[0].Vect()-m_Pvis[1].Vect()).Unit();
123  double pinv = GetP(Minv, m_Minv[0], m_Minv[1]);
124  m_Pinv[0].SetVectM( pinv*Vdiff, m_Minv[0]);
125  m_Pinv[1].SetVectM(-pinv*Vdiff, m_Minv[1]);
126  } else {
127  ApplyOptimalRotation();
128 
129  double k = GetPScale(Minv);
130  for(int i = 0; i < m_Npair; i++)
131  m_Pinv[i].SetVectM(k*m_Pinv[i].Vect(), m_Minv[i]);
132  }
133 
134  // return to original frame
135  for(int i = 0; i < m_Npair; i++){
136  m_Pinv[i].Boost(BoostINV);
137  m_Pinv[i].Boost(BoostCM);
138  GetChildState(i).SetFourVector(m_Pinv[i]);
139  }
140 
141  return SetSpirit(true);
142  }
143 
144  // Based on Newton-Raphson root finding
145  double MinMassesSqInvJigsaw::GetPScale(double Minv){
146  std::vector<double> Pinv2;
147  double Ek = 0.;
148  double fk = 0.;
149  double dfk = 0.;
150  for(int i = 0; i < m_Npair; i++){
151  Pinv2.push_back(m_Pinv[i].P()*m_Pinv[i].P());
152  Ek += sqrt(m_Minv[i]*m_Minv[i]+Pinv2[i]);
153  fk += m_Minv[i];
154  }
155  if(fk > Minv || Ek <= 0.) return 0.;
156 
157  double k2 = Minv/Ek;
158  k2 *= k2;
159  double dk2 = k2;
160  int count = 0;
161  while(fabs(dk2) >= 1e-10 && count < 100){
162  fk = -Minv;
163  dfk = 0.;
164  for(int i = 0; i < m_Npair; i++){
165  Ek = sqrt(m_Minv[i]*m_Minv[i]+k2*Pinv2[i]);
166  fk += Ek;
167  dfk += Ek > 0 ? Pinv2[i]/Ek : 0.;
168  }
169  dk2 = 2.*fk/dfk;
170  k2 -= dk2;
171  count++;
172  }
173 
174  return sqrt(std::max(0.,k2));
175  }
176 
177  // Optimal rotation found using Singular Value Decomposition
178  void MinMassesSqInvJigsaw::ApplyOptimalRotation(){
179  // first, check dimensionality of points
180  TVector3 Z(0.,0.,0.);
181  int index = 0;
182 
183  while(Z.Mag() <= 0. && index < m_Npair){
184  Z = m_Pinv[index].Vect().Cross(m_Pvis[index].Vect());
185  index++;
186  }
187  if(Z.Mag() <= 0.) return; // already aligned
188  Z = Z.Unit();
189 
190  bool b_2D = true;
191  for(int i = 0; i < m_Npair; i++){
192  if(Z.Dot(m_Pvis[i].Vect().Unit()) > 1e-8){
193  b_2D = false;
194  break;
195  }
196  if(Z.Dot(m_Pinv[i].Vect().Unit()) > 1e-8){
197  b_2D = false;
198  break;
199  }
200  }
201 
202  // two dimensional problem, one rotation
203  if(b_2D){
204  TVector3 X(0.,0.,0.);
205  for(int i = 0; i < m_Npair; i++)
206  X += m_Pvis[i].Vect() - m_Pinv[i].Vect();
207  if(X.Mag() <= 0) return; // can't improve
208 
209  X = X.Unit();
210  TVector3 Y = Z.Cross(X).Unit();
211 
212  double num = 0.;
213  double den = 0.;
214  for(int i = 0; i < m_Npair; i++){
215  num += m_Pvis[i].Vect().Dot(Y)*
216  (m_Pvis[i].Vect().Dot(X)-m_Pinv[i].Vect().Dot(X));
217  den += m_Pvis[i].Vect().Dot(m_Pinv[i].Vect());
218  }
219 
220  double theta = atan2(num,den);
221  for(int i = 0; i < m_Npair; i++)
222  m_Pinv[i].Rotate(theta, -Z);
223 
224  return;
225  }
226 
227  // three dimensional problem - R from SVD
228  TMatrixD H(3,3);
229  double val;
230  for(int i = 0; i < 3; i++)
231  for(int j = 0; j < 3; j++)
232  H(i,j) = 0.;
233 
234  for(int p = 0; p < m_Npair; p++){
235  H(0,0) += m_Pinv[p].Px()*m_Pvis[p].Px();
236  H(1,0) += m_Pinv[p].Px()*m_Pvis[p].Py();
237  H(2,0) += m_Pinv[p].Px()*m_Pvis[p].Pz();
238  H(0,1) += m_Pinv[p].Py()*m_Pvis[p].Px();
239  H(1,1) += m_Pinv[p].Py()*m_Pvis[p].Py();
240  H(2,1) += m_Pinv[p].Py()*m_Pvis[p].Pz();
241  H(0,2) += m_Pinv[p].Pz()*m_Pvis[p].Px();
242  H(1,2) += m_Pinv[p].Pz()*m_Pvis[p].Py();
243  H(2,2) += m_Pinv[p].Pz()*m_Pvis[p].Pz();
244  }
245 
246  TDecompSVD SVD(H);
247  SVD.Decompose();
248  TMatrixD UT(TMatrixD::kTransposed,SVD.GetU());
249  TMatrixD R(SVD.GetV(),TMatrixD::kMult,UT);
250 
251  TVector3 V;
252  for(int p = 0; p < m_Npair; p++){
253  V = m_Pinv[p].Vect();
254  V.SetXYZ(R(0,0)*V.Px()+R(1,0)*V.Py()+R(2,0)*V.Pz(),
255  R(0,1)*V.Px()+R(1,1)*V.Py()+R(2,1)*V.Pz(),
256  R(0,2)*V.Px()+R(1,2)*V.Py()+R(2,2)*V.Pz());
257  m_Pinv[p].SetVectM(V, m_Pinv[p].M());
258  }
259 
260  return;
261  }
262 
263  MinMassesSqInvJigsaw MinMassesSqInvJigsaw::m_Empty;
264 
265 }