LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
CombinedCBInvJigsaw.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/CombinedCBInvJigsaw.hh"
31 #include "RestFrames/ContraBoostInvJigsaw.hh"
32 
33 namespace RestFrames {
34 
35  CombinedCBInvJigsaw::CombinedCBInvJigsaw(const std::string& sname,
36  const std::string& stitle,
37  int N_CBjigsaw)
38  : InvisibleJigsaw(sname, stitle, 2*N_CBjigsaw, 2*N_CBjigsaw),
39  m_NCB(N_CBjigsaw)
40  {
41  m_InvMassDependancy = true;
42  }
43 
44  CombinedCBInvJigsaw::~CombinedCBInvJigsaw() {}
45 
46  void CombinedCBInvJigsaw::AddJigsaw(const ContraBoostInvJigsaw& jigsaw, int ijigsaw){
47  if(!jigsaw) return;
48  if(ijigsaw < 0 || ijigsaw >= m_NCB) return;
49 
50  AddInvisibleFrames(jigsaw.GetChildFrames(0), 2*ijigsaw+0);
51  AddInvisibleFrames(jigsaw.GetChildFrames(1), 2*ijigsaw+1);
52  AddVisibleFrames(jigsaw.GetDependancyFrames(0), 2*ijigsaw+0);
53  AddVisibleFrames(jigsaw.GetDependancyFrames(1), 2*ijigsaw+1);
54  }
55 
56  double CombinedCBInvJigsaw::GetMinimumMass() const {
57  if(m_NCB < 1) return 0.;
58  if(m_NCB == 1) return GetCBMinimumMass(0);
59 
60  double Mmin2 = 0.;
61 
62  std::vector<double> Minv2;
63  for(int i = 0; i < m_NCB; i++){
64  Minv2.push_back(GetCBMinimumMass(i));
65  Minv2[i] *= Minv2[i];
66  Mmin2 += Minv2[i];
67  }
68 
69  std::vector<TLorentzVector> Pvis;
70  std::vector<double> Mvis2;
71  for(int i = 0; i < m_NCB; i++){
72  Pvis.push_back(GetDependancyStates(2*i+0).GetFourVector()+
73  GetDependancyStates(2*i+1).GetFourVector());
74  Mvis2.push_back(Pvis[i].M2());
75  }
76 
77  for(int i = 0; i < m_NCB-1; i++){
78  for(int j = i+1; j < m_NCB; j++){
79  if(Minv2[i] > Mvis2[i] || Minv2[j] > Mvis2[j]){
80  Mmin2 += std::max(0.,(Pvis[i]+Pvis[j]).M2()-Mvis2[i]-Mvis2[j])*
81  std::max(((Minv2[i] > 0. && Mvis2[i] <= 0.) ? 1. : Minv2[i]/Mvis2[i]),
82  ((Minv2[j] > 0. && Mvis2[j] <= 0.) ? 1. : Minv2[j]/Mvis2[j]));
83 
84  } else {
85  Mmin2 += std::max(0.,(Pvis[i]+Pvis[j]).M2()-Mvis2[i]-Mvis2[j]);
86  Mmin2 += std::max(Minv2[i] - Mvis2[i], Minv2[j] - Mvis2[j])*2.;
87  }
88  }
89  }
90 
91  return sqrt(std::max(0., Mmin2));
92  }
93 
94  double CombinedCBInvJigsaw::GetCBMinimumMass(int i) const {
95  if(!IsSoundMind())
96  return 0.;
97 
98  double Minv1 = GetChildState(2*i+0).GetMinimumMass();
99  double Minv2 = GetChildState(2*i+1).GetMinimumMass();
100  TLorentzVector Pvis1 = GetDependancyStates(2*i+0).GetFourVector();
101  TLorentzVector Pvis2 = GetDependancyStates(2*i+1).GetFourVector();
102  double Mvis1 = std::max(Pvis1.M(), 0.);
103  double Mvis2 = std::max(Pvis2.M(), 0.);
104  double M12 = (Pvis1+Pvis2).M();
105 
106  if(Minv1 < -0.5 && Minv2 < -0.5) // children can go tachyonic
107  return 2.*GetP(M12,Mvis1,Mvis2);
108 
109  Minv1 = std::max(Minv1,0.);
110  Minv2 = std::max(Minv2,0.);
111 
112  double Minvmax = std::max(0.,std::max(Minv1,Minv2));
113 
114  double Mvismin = std::min(Mvis1,Mvis2);
115  double Mvismax = std::max(Mvis1,Mvis2);
116 
117  if(Minv1 < Mvis2 && Minv2 < Mvis1){
118  if(Minvmax <= Mvismin)
119  return sqrt( M12*M12 + 4.*(Minvmax-Mvismin)*(Minvmax+Mvismax) );
120  return M12;
121  }
122 
123  if(Mvismin <= 0.0 && Minvmax > 0.)
124  return M12;
125 
126  return M12*(1.+sqrt(std::max(Minv1*Minv1-Mvis2*Mvis2,
127  Minv2*Minv2-Mvis1*Mvis1))/Mvismin);
128  }
129 
130  bool CombinedCBInvJigsaw::AnalyzeEvent(){
131  if(!IsSoundMind())
132  return SetSpirit(false);
133 
134  if(m_NCB < 1)
135  return SetSpirit(false);
136 
137  TLorentzVector INV = GetParentState().GetFourVector();
138  TLorentzVector VIS(0.,0.,0.,0.);
139 
140  std::vector<TLorentzVector> Pvis[2];
141  for(int i = 0; i < m_NCB; i++){
142  Pvis[0].push_back( GetDependancyStates(2*i+0).GetFourVector() );
143  Pvis[1].push_back( GetDependancyStates(2*i+1).GetFourVector() );
144  VIS += Pvis[0][i]+Pvis[1][i];
145  }
146 
147  // go to the rest frame of (VIS+INV system)
148  TVector3 Boost = (VIS+INV).BoostVector();
149  for(int i = 0; i < m_NCB; i++){
150  Pvis[0][i].Boost(-Boost);
151  Pvis[1][i].Boost(-Boost);
152  }
153 
154  std::vector<double> c[2];
155  std::vector<double> E[2];
156  std::vector<double> mvis[2];
157  std::vector<double> minv[2];
158  std::vector<double> Minv2;
159  double mvismin, minvmax;
160  double m1, m2, k1, k2;
161  double MC2, Xbar, K;
162  double N, sumE, sumcE;
163  for(int i = 0; i < m_NCB; i++){
164  E[0].push_back(Pvis[0][i].E());
165  E[1].push_back(Pvis[1][i].E());
166  mvis[0].push_back(std::max(0.,Pvis[0][i].M()));
167  mvis[1].push_back(std::max(0.,Pvis[1][i].M()));
168  minv[0].push_back( GetChildState(2*i+0).GetMinimumMass() );
169  minv[1].push_back( GetChildState(2*i+1).GetMinimumMass() );
170  Minv2.push_back( GetCBMinimumMass(i) );
171  Minv2[i] *= Minv2[i];
172 
173  minvmax = std::max(0.,std::max(minv[0][i],minv[1][i]));
174  mvismin = std::min(mvis[0][i],mvis[1][i]);
175 
176  c[0].push_back(1.);
177  c[1].push_back(1.);
178  if(minvmax < mvismin){
179  m1 = mvis[0][i];
180  m2 = mvis[1][i];
181  MC2 = 2.*(E[0][i]*E[1][i] + Pvis[0][i].Vect().Dot(Pvis[1][i].Vect()));
182  k1 = (m1+m2)*(m1-m2)*(1.-minvmax/mvismin) + MC2-2*m1*m2 +
183  (m1+m2)*fabs(m1-m2)*minvmax/mvismin;
184  k2 = -(m1+m2)*(m1-m2)*(1.-minvmax/mvismin) + MC2-2*m1*m2 +
185  (m1+m2)*fabs(m1-m2)*minvmax/mvismin;
186  Xbar = sqrt( (k1+k2)*(k1+k2)*(MC2*MC2-4*m1*m1*m2*m2) +
187  16.*minvmax*minvmax*(k1*k1*m1*m1 + k2*k2*m2*m2 + k1*k2*MC2) );
188  K = ( fabs(k1*m1*m1-k2*m2*m2) - 0.5*fabs(k2-k1)*MC2 + 0.5*Xbar )/
189  (k1*k1*m1*m1 + k2*k2*m2*m2 + k1*k2*MC2);
190  c[0][i] = 0.5*(1.+K*k1);
191  c[1][i] = 0.5*(1.+K*k2);
192  }
193 
194  sumE = E[0][i] + E[1][i];
195  sumcE = c[0][i]*E[0][i] + c[1][i]*E[1][i];
196 
197  N = (sqrt(sumE*sumE+Minv2[i]-(Pvis[0][i]+Pvis[1][i]).M2())+sumE)/sumcE/2.;
198  //N = sumkE > 0. ? sqrt(std::max(0.,Minv2[i]-(Pvis[0][i]+Pvis[1][i]).M2()+sumE*sumE))/sumkE : 0.;
199 
200  c[0][i] *= N;
201  c[1][i] *= N;
202  }
203 
204  sumE = 0.;
205  sumcE = 0.;
206  for(int i = 0; i < m_NCB; i++){
207  sumE += E[0][i] + E[1][i];
208  sumcE += c[0][i]*E[0][i] + c[1][i]*E[1][i];
209  }
210 
211  N = (sqrt(sumE*sumE + INV.M2() - VIS.M2())+sumE)/sumcE/2.;
212  //N = sumkE > 0. ? sqrt(std::max(0.,INV.M2()-VIS.M2()+sumE*sumE))/sumkE : 0.;
213 
214  TLorentzVector INV1, INV2;
215  double Einv1, Einv2;
216  TVector3 Pinv1, Pinv2;
217  for(int i = 0; i < m_NCB; i++){
218  c[0][i] *= N;
219  c[1][i] *= N;
220 
221  Einv1 = (c[0][i]-1.)*E[0][i] + c[1][i]*E[1][i];
222  Einv2 = (c[1][i]-1.)*E[1][i] + c[0][i]*E[0][i];
223  Pinv1 = (c[0][i]-1.)*Pvis[0][i].Vect() - c[1][i]*Pvis[1][i].Vect();
224  Pinv2 = (c[1][i]-1.)*Pvis[1][i].Vect() - c[0][i]*Pvis[0][i].Vect();
225 
226  INV1.SetPxPyPzE(Pinv1.X(),Pinv1.Y(),Pinv1.Z(),Einv1);
227  INV2.SetPxPyPzE(Pinv2.X(),Pinv2.Y(),Pinv2.Z(),Einv2);
228 
229  if(minv[0][i] >= 0. && INV1.M() < minv[0][i])
230  INV1.SetVectM(Pinv1, minv[0][i]);
231  if(minv[1][i] >= 0. && INV2.M() < minv[1][i])
232  INV2.SetVectM(Pinv2, minv[1][i]);
233 
234  INV1.Boost(Boost);
235  INV2.Boost(Boost);
236 
237  GetChildState(2*i+0).SetFourVector(INV1);
238  GetChildState(2*i+1).SetFourVector(INV2);
239  }
240 
241  return SetSpirit(true);
242  }
243 
244 }
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1