LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
ContraBoostInvJigsaw.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 
31 
32 namespace RestFrames {
33 
34  ContraBoostInvJigsaw::ContraBoostInvJigsaw(const std::string& sname,
35  const std::string& stitle) :
36  InvisibleJigsaw(sname, stitle, 2, 2)
37  {
38  m_InvMassDependancy = true;
39  }
40 
41  ContraBoostInvJigsaw::ContraBoostInvJigsaw() : InvisibleJigsaw() {}
42 
43  ContraBoostInvJigsaw::~ContraBoostInvJigsaw() {}
44 
45  ContraBoostInvJigsaw& ContraBoostInvJigsaw::Empty(){
46  return ContraBoostInvJigsaw::m_Empty;
47  }
48 
49  double ContraBoostInvJigsaw::GetMinimumMass() const {
50  if(!IsSoundMind())
51  return 0.;
52 
53  double Minv1 = GetChildState(0).GetMinimumMass();
54  double Minv2 = GetChildState(1).GetMinimumMass();
55  TLorentzVector Pvis1 = GetDependancyStates(0).GetFourVector();
56  TLorentzVector Pvis2 = GetDependancyStates(1).GetFourVector();
57  double Mvis1 = std::max(Pvis1.M(), 0.);
58  double Mvis2 = std::max(Pvis2.M(), 0.);
59  double M12 = (Pvis1+Pvis2).M();
60 
61  if(Minv1 < -0.5 && Minv2 < -0.5) // children can go tachyonic
62  return 2.*GetP(M12,Mvis1,Mvis2);
63 
64  Minv1 = std::max(Minv1,0.);
65  Minv2 = std::max(Minv2,0.);
66 
67  double Minvmax = std::max(0.,std::max(Minv1,Minv2));
68 
69  double Mvismin = std::min(Mvis1,Mvis2);
70  double Mvismax = std::max(Mvis1,Mvis2);
71 
72  if(Minv1 < Mvis2 && Minv2 < Mvis1){
73  if(Minvmax <= Mvismin)
74  return sqrt( M12*M12 + 4.*(Minvmax-Mvismin)*(Minvmax+Mvismax) );
75  return M12;
76  }
77 
78  if(Mvismin <= 0.0 && Minvmax > 0.)
79  return M12;
80 
81  return M12*(1.+sqrt(std::max(Minv1*Minv1-Mvis2*Mvis2,
82  Minv2*Minv2-Mvis1*Mvis1))/Mvismin);
83  }
84 
85  bool ContraBoostInvJigsaw::AnalyzeEvent(){
86  if(!IsSoundMind())
87  return SetSpirit(false);
88 
89  TLorentzVector Pvis1 = GetDependancyStates(0).GetFourVector();
90  TLorentzVector Pvis2 = GetDependancyStates(1).GetFourVector();
91  TLorentzVector INV = GetParentState().GetFourVector();
92 
93  // go to the rest frame of (Pvis1+Pvis2+INV system)
94  TVector3 Boost = (Pvis1+Pvis2+INV).BoostVector();
95  Pvis1.Boost(-Boost);
96  Pvis2.Boost(-Boost);
97  INV.Boost(-Boost);
98 
99  double E1 = Pvis1.E();
100  double E2 = Pvis2.E();
101  double m1 = std::max(0.,Pvis1.M());
102  double m2 = std::max(0.,Pvis2.M());
103  TVector3 P1 = Pvis1.Vect();
104  TVector3 P2 = Pvis2.Vect();
105 
106  double Minv1 = GetChildState(0).GetMinimumMass();
107  double Minv2 = GetChildState(1).GetMinimumMass();
108  double Minv = std::max(0.,std::max(Minv1,Minv2));
109  double Mvis = std::min(m1,m2);
110 
111  double c1 = 1.;
112  double c2 = 1.;
113  if(Minv < Mvis){
114  double MC2 = 2.*(E1*E2 + P1.Dot(P2));
115  double k1 = (m1+m2)*(m1-m2)*(1.-Minv/Mvis) + MC2-2*m1*m2 + (m1+m2)*fabs(m1-m2)*Minv/Mvis;
116  double k2 = -(m1+m2)*(m1-m2)*(1.-Minv/Mvis) + MC2-2*m1*m2 + (m1+m2)*fabs(m1-m2)*Minv/Mvis;
117  double Xbar = sqrt( (k1+k2)*(k1+k2)*(MC2*MC2-4*m1*m1*m2*m2) +
118  16.*Minv*Minv*(k1*k1*m1*m1 + k2*k2*m2*m2 + k1*k2*MC2) );
119  double K = ( fabs(k1*m1*m1-k2*m2*m2) - 0.5*fabs(k2-k1)*MC2 + 0.5*Xbar )/
120  (k1*k1*m1*m1 + k2*k2*m2*m2 + k1*k2*MC2);
121  c1 = 0.5*(1.+K*k1);
122  c2 = 0.5*(1.+K*k2);
123  }
124 
125  double sumE = E1+E2;
126  double sumcE = c1*E1+c2*E2;
127 
128  double N = (sqrt(sumE*sumE-(Pvis1+Pvis2).M2()+INV.M2())+sumE)/sumcE/2.;
129 
130  c1 *= N;
131  c2 *= N;
132 
133  TLorentzVector INV1,INV2;
134  double Einv1 = (c1-1.)*E1 + c2*E2;
135  double Einv2 = c1*E1 + (c2-1.)*E2;
136  TVector3 Pinv1 = (c1-1.)*P1 - c2*P2;
137  TVector3 Pinv2 = (c2-1.)*P2 - c1*P1;
138 
139  INV1.SetPxPyPzE(Pinv1.X(),Pinv1.Y(),Pinv1.Z(),Einv1);
140  INV2.SetPxPyPzE(Pinv2.X(),Pinv2.Y(),Pinv2.Z(),Einv2);
141 
142  if(Minv1 >= 0. && INV1.M() < Minv1)
143  INV1.SetVectM(Pinv1,Minv1);
144  if(Minv2 >= 0. && INV2.M() < Minv2)
145  INV2.SetVectM(Pinv2,Minv2);
146 
147  INV1.Boost(Boost);
148  INV2.Boost(Boost);
149 
150  GetChildState(0).SetFourVector(INV1);
151  GetChildState(1).SetFourVector(INV2);
152 
153  return SetSpirit(true);
154  }
155 
156  ContraBoostInvJigsaw ContraBoostInvJigsaw::m_Empty;
157 
158 }