LOGO
RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
MinMassDiffInvJigsaw.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/Factory.h"
31 
32 #include "RestFrames/MinMassDiffInvJigsaw.hh"
33 
34 namespace RestFrames {
35 
36  MinMassDiffInvJigsaw::MinMassDiffInvJigsaw(const std::string& sname,
37  const std::string& stitle,
38  int N_vis_inv_pair) :
39  InvisibleJigsaw(sname, stitle, N_vis_inv_pair, N_vis_inv_pair),
40  m_Npair(N_vis_inv_pair)
41  {
42  m_minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
43  m_minimizer->SetMaxFunctionCalls(10000000);
44  m_minimizer->SetMaxIterations(100000);
45  m_minimizer->SetTolerance(0.001);
46  m_minimizer->SetPrintLevel(0);
47 
48  m_functor = new ROOT::Math::Functor(this, &MinMassDiffInvJigsaw::EvaluateMetric, 3);
49  m_minimizer->SetFunction(*m_functor);
50 
51  m_minimizer->SetVariable(0, "phi1", 0., 0.001);
52  m_minimizer->SetVariable(1, "phi2", 0., 0.001);
53  m_minimizer->SetVariable(2, "phi3", 0., 0.001);
54  }
55 
56  MinMassDiffInvJigsaw::~MinMassDiffInvJigsaw(){
57  delete m_minimizer;
58  delete m_functor;
59  }
60 
61  double MinMassDiffInvJigsaw::GetMinimumMass() const {
62  if(!IsSoundMind())
63  return 0.;
64 
65  if(m_Npair <= 0)
66  return 0.;
67  if(m_Npair == 1)
68  return std::max(0.,GetChildState(0).GetMinimumMass());
69 
70  TLorentzVector PvisTOT(0.,0.,0.,0.);
71  m_Pvis.clear();
72  for(int i = 0; i < m_Npair; i++){
73  m_Pvis.push_back(GetDependancyStates(i+m_Npair).GetFourVector());
74  PvisTOT += m_Pvis[i];
75  }
76 
77  double ECM = 0.;
78  double Minv;
79  TVector3 Boost = PvisTOT.BoostVector();
80  for(int i = 0; i < m_Npair; i++){
81  m_Pvis[i].Boost(-Boost);
82  Minv = std::max(0.,GetChildState(i).GetMinimumMass());
83  ECM += sqrt(m_Pvis[i].P()*m_Pvis[i].P() + Minv*Minv);
84  }
85 
86  return ECM;
87  }
88 
89  bool MinMassDiffInvJigsaw::AnalyzeEvent(){
90  if(!IsSoundMind())
91  return SetSpirit(false);
92 
93  if(m_Npair <= 1)
94  return false;
95 
96  TLorentzVector INV = GetParentState().GetFourVector();
97  double Minv = INV.M();
98 
99  TLorentzVector VIS(0.,0.,0.,0.);
100  m_Pvis.clear();
101  for(int i = 0; i < m_Npair; i++){
102  m_Pvis.push_back(GetDependancyStates(i).GetFourVector());
103  VIS += m_Pvis[i];
104  }
105 
106  // go to INV+VIS CM frame
107  TVector3 BoostCM = (VIS+INV).BoostVector();
108  INV.Boost(-BoostCM);
109  VIS.Boost(-BoostCM);
110  for(int i = 0; i < m_Npair; i++)
111  m_Pvis[i].Boost(-BoostCM);
112 
113  // INV states defined in the VIS frame
114  TVector3 BoostVIS = VIS.BoostVector();
115  m_Pinv.clear();
116  m_Minv.clear();
117  for(int i = 0; i < m_Npair; i++){
118  m_Minv.push_back(std::max(0.,GetChildState(i).GetMinimumMass()));
119  m_Pinv.push_back(m_Pvis[i]);
120  m_Pinv[i].Boost(-BoostVIS);
121  m_Pinv[i].SetVectM(m_Pinv[i].Vect(), m_Minv[i]);
122  }
123 
124  // VIS states in INV frame
125  TVector3 BoostINV = INV.BoostVector();
126  for(int i = 0; i < m_Npair; i++)
127  m_Pvis[i].Boost(-BoostINV);
128 
129  if(m_Npair == 2){
130  TVector3 Vdiff = (m_Pvis[0].Vect()-m_Pvis[1].Vect()).Unit();
131  double pinv = GetP(Minv, m_Minv[0], m_Minv[1]);
132  m_Pinv[0].SetVectM( pinv*Vdiff, m_Minv[0]);
133  m_Pinv[1].SetVectM(-pinv*Vdiff, m_Minv[1]);
134  } else {
135  double k = GetPScale(Minv);
136  for(int i = 0; i < m_Npair; i++)
137  m_Pinv[i].SetVectM(k*m_Pinv[i].Vect(), m_Minv[i]);
138  }
139 
140  ApplyOptimalRotation();
141 
142  // return to original frame
143  for(int i = 0; i < m_Npair; i++){
144  m_Pinv[i].Boost(BoostINV);
145  m_Pinv[i].Boost(BoostCM);
146  GetChildState(i).SetFourVector(m_Pinv[i]);
147  }
148 
149  return SetSpirit(true);
150  }
151 
152  // Based on Newton-Raphson root finding
153  double MinMassDiffInvJigsaw::GetPScale(double Minv){
154  std::vector<double> Pinv2;
155  double Ek = 0.;
156  double fk = 0.;
157  double dfk = 0.;
158  for(int i = 0; i < m_Npair; i++){
159  Pinv2.push_back(m_Pinv[i].P()*m_Pinv[i].P());
160  Ek += sqrt(m_Minv[i]*m_Minv[i]+Pinv2[i]);
161  fk += m_Minv[i];
162  }
163  if(fk > Minv || Ek <= 0.) return 0.;
164 
165  double k2 = Minv/Ek;
166  k2 *= k2;
167  double dk2 = k2;
168  int count = 0;
169  while(fabs(dk2) >= 1e-10 && count < 100){
170  fk = -Minv;
171  dfk = 0.;
172  for(int i = 0; i < m_Npair; i++){
173  Ek = sqrt(m_Minv[i]*m_Minv[i]+k2*Pinv2[i]);
174  fk += Ek;
175  dfk += Ek > 0 ? Pinv2[i]/Ek : 0.;
176  }
177  dk2 = 2.*fk/dfk;
178  k2 -= dk2;
179  count++;
180  }
181 
182  return sqrt(std::max(0.,k2));
183  }
184 
185  // Optimal rotation found using Minuit2 ROOT interface
186  void MinMassDiffInvJigsaw::ApplyOptimalRotation(){
187  // first, check dimensionality of points
188  m_Z.SetXYZ(0.,0.,0.);
189  int index = 0;
190 
191  while(m_Z.Mag() <= 0. && index < m_Npair){
192  m_Z = m_Pinv[index].Vect().Cross(m_Pvis[index].Vect());
193  index++;
194  }
195  m_Z = m_Z.Unit();
196 
197  m_2D = true;
198  if(m_Npair > 2){
199  for(int i = 0; i < m_Npair; i++){
200  if(m_Z.Dot(m_Pvis[i].Vect().Unit()) > 1e-8){
201  m_2D = false;
202  break;
203  }
204  if(m_Z.Dot(m_Pinv[i].Vect().Unit()) > 1e-8){
205  m_2D = false;
206  break;
207  }
208  }
209  }
210 
211  m_minimizer->SetVariableValue(0, 0.);
212  m_minimizer->SetVariableValue(1, 0.);
213  m_minimizer->SetVariableValue(2, 0.);
214 
215  if(m_2D){
216  m_minimizer->FixVariable(1);
217  m_minimizer->FixVariable(2);
218  } else {
219  m_minimizer->ReleaseVariable(1);
220  m_minimizer->ReleaseVariable(2);
221  m_X.SetXYZ(0.,0.,0.);
222  for(int i = 0; i < m_Npair; i++)
223  m_X += m_Pvis[i].Vect() - m_Pinv[i].Vect();
224 
225  if(m_X.Mag() <= 0)
226  m_X = m_Pvis[0].Vect().Cross(m_Z);
227 
228  m_X = m_X.Unit();
229  m_Y = m_Z.Cross(m_X).Unit();
230  }
231 
232  m_minimizer->Minimize();
233  if(m_minimizer->Status() > 0 && !m_2D){
234  TVector3 temp = m_Z;
235  m_Z = m_X;
236  m_X = m_Y;
237  m_Y = temp;
238  m_minimizer->SetVariableValue(0, 0.);
239  m_minimizer->SetVariableValue(1, 0.);
240  m_minimizer->SetVariableValue(2, 0.);
241  m_minimizer->Minimize();
242  }
243 
244  const double* PHIs = m_minimizer->X();
245 
246  for(int i = 0; i < m_Npair; i++){
247  m_Pinv[i].Rotate(PHIs[0], m_Z);
248  if(!m_2D){
249  m_Pinv[i].Rotate(PHIs[1], m_Y);
250  m_Pinv[i].Rotate(PHIs[2], m_X);
251  }
252  }
253 
254  return;
255  }
256 
257  double MinMassDiffInvJigsaw::EvaluateMetric(const double* param){
258  std::vector<TLorentzVector> Pinv;
259  for(int i = 0; i < m_Npair; i++){
260  Pinv.push_back(m_Pinv[i]);
261  Pinv[i].Rotate(param[0], m_Z);
262  if(!m_2D){
263  Pinv[i].Rotate(param[1], m_Y);
264  Pinv[i].Rotate(param[2], m_X);
265  }
266  Pinv[i] += m_Pvis[i];
267  }
268  double diff = 0.;
269  for(int i = 0; i < m_Npair-1; i++){
270  for(int j = i+1; j < m_Npair; j++){
271  diff += (Pinv[i].M()-Pinv[j].M())*(Pinv[i].M()-Pinv[j].M());
272  }
273  }
274 
275  return diff;
276  }
277 
278 }
RestFrames::State::SetFourVector
virtual void SetFourVector(const TLorentzVector &V)
Sets four-vector of this frame.
Definition: State.cc:161
RestFrames::MinMassDiffInvJigsaw::MinMassDiffInvJigsaw
MinMassDiffInvJigsaw(const std::string &sname, const std::string &stitle, int N_vis_inv_pair)
Standard constructor.
Definition: MinMassDiffInvJigsaw.cc:36
RestFrames::State::GetFourVector
virtual TLorentzVector GetFourVector() const
Returns four vector of this frame.
Definition: State.cc:165
RestFrames::MinMassDiffInvJigsaw::GetMinimumMass
virtual double GetMinimumMass() const
Returns minimum Lorentz invariant mass.
Definition: MinMassDiffInvJigsaw.cc:61
RestFrames::MinMassDiffInvJigsaw::AnalyzeEvent
virtual bool AnalyzeEvent()
Analyzes event for this jigsaw.
Definition: MinMassDiffInvJigsaw.cc:89
RestFrames::InvisibleJigsaw
Definition: InvisibleJigsaw.hh:39
MinMassDiffInvJigsaw.hh
Generated on Thu Dec 17 2020 14:04:48 for RestFrames by   doxygen 1.8.17