LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
MaxProbBreitWignerInvJigsaw.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/MaxProbBreitWignerInvJigsaw.hh"
33 
34 namespace RestFrames {
35 
36  MaxProbBreitWignerInvJigsaw::MaxProbBreitWignerInvJigsaw(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, &MaxProbBreitWignerInvJigsaw::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  for(int i = 0; i < m_Npair; i++){
56  m_Mass.push_back(0.);
57  m_Width.push_back(-1.);
58  }
59  }
60 
61  MaxProbBreitWignerInvJigsaw::~MaxProbBreitWignerInvJigsaw(){
62  delete m_minimizer;
63  delete m_functor;
64  }
65 
66  double MaxProbBreitWignerInvJigsaw::GetMinimumMass() const {
67  if(!IsSoundMind())
68  return 0.;
69 
70  if(m_Npair <= 0)
71  return 0.;
72  if(m_Npair == 1)
73  return std::max(0.,GetChildState(0).GetMinimumMass());
74 
75  TLorentzVector PvisTOT(0.,0.,0.,0.);
76  m_Pvis.clear();
77  for(int i = 0; i < m_Npair; i++){
78  m_Pvis.push_back(GetDependancyStates(i+m_Npair).GetFourVector());
79  PvisTOT += m_Pvis[i];
80  }
81 
82  double ECM = 0.;
83  double Minv;
84  TVector3 Boost = PvisTOT.BoostVector();
85  for(int i = 0; i < m_Npair; i++){
86  m_Pvis[i].Boost(-Boost);
87  Minv = std::max(0.,GetChildState(i).GetMinimumMass());
88  ECM += sqrt(m_Pvis[i].P()*m_Pvis[i].P() + Minv*Minv);
89  }
90 
91  return ECM;
92  }
93 
94  bool MaxProbBreitWignerInvJigsaw::AnalyzeEvent(){
95  if(!IsSoundMind())
96  return SetSpirit(false);
97 
98  if(m_Npair <= 1)
99  return false;
100 
101  TLorentzVector INV = GetParentState().GetFourVector();
102  double Minv = INV.M();
103 
104  TLorentzVector VIS(0.,0.,0.,0.);
105  m_Pvis.clear();
106  for(int i = 0; i < m_Npair; i++){
107  m_Pvis.push_back(GetDependancyStates(i).GetFourVector());
108  VIS += m_Pvis[i];
109  }
110 
111  // go to INV+VIS CM frame
112  TVector3 BoostCM = (VIS+INV).BoostVector();
113  INV.Boost(-BoostCM);
114  VIS.Boost(-BoostCM);
115  for(int i = 0; i < m_Npair; i++)
116  m_Pvis[i].Boost(-BoostCM);
117 
118  // INV states defined in the VIS frame
119  TVector3 BoostVIS = VIS.BoostVector();
120  m_Pinv.clear();
121  m_Minv.clear();
122  for(int i = 0; i < m_Npair; i++){
123  m_Minv.push_back(std::max(0.,GetChildState(i).GetMinimumMass()));
124  m_Pinv.push_back(m_Pvis[i]);
125  m_Pinv[i].Boost(-BoostVIS);
126  m_Pinv[i].SetVectM(m_Pinv[i].Vect(), m_Minv[i]);
127  }
128 
129  // VIS states in INV frame
130  TVector3 BoostINV = INV.BoostVector();
131  for(int i = 0; i < m_Npair; i++)
132  m_Pvis[i].Boost(-BoostINV);
133 
134  if(m_Npair == 2){
135  TVector3 Vdiff = (m_Pvis[0].Vect()-m_Pvis[1].Vect()).Unit();
136  double pinv = GetP(Minv, m_Minv[0], m_Minv[1]);
137  m_Pinv[0].SetVectM( pinv*Vdiff, m_Minv[0]);
138  m_Pinv[1].SetVectM(-pinv*Vdiff, m_Minv[1]);
139  } else {
140  double k = GetPScale(Minv);
141  for(int i = 0; i < m_Npair; i++)
142  m_Pinv[i].SetVectM(k*m_Pinv[i].Vect(), m_Minv[i]);
143  }
144 
145  ApplyOptimalRotation();
146 
147  // return to original frame
148  for(int i = 0; i < m_Npair; i++){
149  m_Pinv[i].Boost(BoostINV);
150  m_Pinv[i].Boost(BoostCM);
151  GetChildState(i).SetFourVector(m_Pinv[i]);
152  }
153 
154  return SetSpirit(true);
155  }
156 
157  // Based on Newton-Raphson root finding
158  double MaxProbBreitWignerInvJigsaw::GetPScale(double Minv){
159  std::vector<double> Pinv2;
160  double Ek = 0.;
161  double fk = 0.;
162  double dfk = 0.;
163  for(int i = 0; i < m_Npair; i++){
164  Pinv2.push_back(m_Pinv[i].P()*m_Pinv[i].P());
165  Ek += sqrt(m_Minv[i]*m_Minv[i]+Pinv2[i]);
166  fk += m_Minv[i];
167  }
168  if(fk > Minv || Ek <= 0.) return 0.;
169 
170  double k2 = Minv/Ek;
171  k2 *= k2;
172  double dk2 = k2;
173  int count = 0;
174  while(fabs(dk2) >= 1e-10 && count < 100){
175  fk = -Minv;
176  dfk = 0.;
177  for(int i = 0; i < m_Npair; i++){
178  Ek = sqrt(m_Minv[i]*m_Minv[i]+k2*Pinv2[i]);
179  fk += Ek;
180  dfk += Ek > 0 ? Pinv2[i]/Ek : 0.;
181  }
182  dk2 = 2.*fk/dfk;
183  k2 -= dk2;
184  count++;
185  }
186 
187  return sqrt(std::max(0.,k2));
188  }
189 
190  // Optimal rotation found using Minuit2 ROOT interface
191  void MaxProbBreitWignerInvJigsaw::ApplyOptimalRotation(){
192  // first, check dimensionality of points
193  m_Z.SetXYZ(0.,0.,0.);
194  int index = 0;
195 
196  while(m_Z.Mag() <= 0. && index < m_Npair){
197  m_Z = m_Pinv[index].Vect().Cross(m_Pvis[index].Vect());
198  index++;
199  }
200  m_Z = m_Z.Unit();
201 
202  m_2D = true;
203  if(m_Npair > 2){
204  for(int i = 0; i < m_Npair; i++){
205  if(m_Z.Dot(m_Pvis[i].Vect().Unit()) > 1e-8){
206  m_2D = false;
207  break;
208  }
209  if(m_Z.Dot(m_Pinv[i].Vect().Unit()) > 1e-8){
210  m_2D = false;
211  break;
212  }
213  }
214  }
215 
216  m_minimizer->SetVariableValue(0, 0.);
217  m_minimizer->SetVariableValue(1, 0.);
218  m_minimizer->SetVariableValue(2, 0.);
219 
220  if(m_2D){
221  m_minimizer->FixVariable(1);
222  m_minimizer->FixVariable(2);
223  } else {
224  m_minimizer->ReleaseVariable(1);
225  m_minimizer->ReleaseVariable(2);
226  m_X.SetXYZ(0.,0.,0.);
227  for(int i = 0; i < m_Npair; i++)
228  m_X += m_Pvis[i].Vect() - m_Pinv[i].Vect();
229 
230  if(m_X.Mag() <= 0)
231  m_X = m_Pvis[0].Vect().Cross(m_Z);
232 
233  m_X = m_X.Unit();
234  m_Y = m_Z.Cross(m_X).Unit();
235  }
236 
237  m_minimizer->Minimize();
238  if(m_minimizer->Status() > 0 && !m_2D){
239  TVector3 temp = m_Z;
240  m_Z = m_X;
241  m_X = m_Y;
242  m_Y = temp;
243  m_minimizer->SetVariableValue(0, 0.);
244  m_minimizer->SetVariableValue(1, 0.);
245  m_minimizer->SetVariableValue(2, 0.);
246  m_minimizer->Minimize();
247  }
248 
249  const double* PHIs = m_minimizer->X();
250 
251  for(int i = 0; i < m_Npair; i++){
252  m_Pinv[i].Rotate(PHIs[0], m_Z);
253  if(!m_2D){
254  m_Pinv[i].Rotate(PHIs[1], m_Y);
255  m_Pinv[i].Rotate(PHIs[2], m_X);
256  }
257  }
258 
259  return;
260  }
261 
262  double MaxProbBreitWignerInvJigsaw::EvaluateMetric(const double* param){
263  std::vector<TLorentzVector> Pinv;
264  for(int i = 0; i < m_Npair; i++){
265  Pinv.push_back(m_Pinv[i]);
266  Pinv[i].Rotate(param[0], m_Z);
267  if(!m_2D){
268  Pinv[i].Rotate(param[1], m_Y);
269  Pinv[i].Rotate(param[2], m_X);
270  }
271  Pinv[i] += m_Pvis[i];
272  }
273 
274  double prob = 1.;
275  TLorentzVector SUM(0.,0.,0.,0.);
276  for(int i = 0; i < m_Npair; i++)
277  SUM += Pinv[i];
278  double M = SUM.M();
279  for(int i = 0; i < m_Npair-1; i++){
280  SUM -= Pinv[i];
281  prob *= GetP((Pinv[i]+SUM).M(), Pinv[i].M(), SUM.M())/M;
282  }
283 
284  double den;
285  for(int i = 0; i < m_Npair; i++){
286  if(m_Width[i] > 0.){
287  den = Pinv[i].M2()-m_Mass[i]*m_Mass[i];
288  den *= den;
289  den += m_Mass[i]*m_Mass[i]*m_Width[i]*m_Width[i];
290  if(den > 0.)
291  prob /= den;
292  }
293  }
294 
295  if(prob <= 0.)
296  return 1e10;
297 
298  return 1./prob;
299  }
300 
301  void MaxProbBreitWignerInvJigsaw::SetPoleMass(double mass, int i){
302  if(i < 0 || i >= m_Npair)
303  return;
304 
305  m_Mass[i] = std::max(mass, 0.);
306  }
307 
308  void MaxProbBreitWignerInvJigsaw::SetWidth(double width, int i){
309  if(i < 0 || i >= m_Npair)
310  return;
311  if(width <= 0.)
312  m_Width[i] = -1.;
313  else
314  m_Width[i] = width;
315  }
316 
317 }
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1