LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
DecayGenFrame.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 <stdlib.h>
32 
33 namespace RestFrames {
34 
36  // DecayGenFrame class
38 
39  DecayGenFrame::DecayGenFrame(const std::string& sname,
40  const std::string& stitle)
41  : DecayFrame<GeneratorFrame>(sname,stitle)
42  {
43  m_CosDecayAngle = -2.;
44  m_DeltaPhiDecayPlane = -2.;
45  }
46 
47  DecayGenFrame::DecayGenFrame() : DecayFrame<GeneratorFrame>() {}
48 
49  DecayGenFrame::~DecayGenFrame() {}
50 
51  bool DecayGenFrame::IsSoundBody() const{
52  if(RFBase::IsSoundBody())
53  return true;
54  if(!RestFrame::IsSoundBody())
55  return false;
56  int Nchild = GetNChildren();
57  if(Nchild < 2 || GetParentFrame().IsEmpty())
58  return SetBody(false);
59  return SetBody(true);
60  }
61 
62  void DecayGenFrame::SetMass(double val){
63  SetMind(false);
64  if(val < 0.){
65  m_Log << LogWarning;
66  m_Log << "Unable to set mass to negative value ";
67  m_Log << val << ". Setting to zero." << LogEnd;
68  m_Mass = 0.;
69  } else {
70  m_Mass = val;
71  }
72 
73  }
74 
75  void DecayGenFrame::SetCosDecayAngle(double val){
76  if(val < 0.){
77  m_Log << LogWarning;
78  m_Log << "CosDecay angle must be in [-1, 1]: ";
79  m_Log << val << ". Setting to random." << LogEnd;
80  m_CosDecayAngle = -2.;
81  } else {
82  m_CosDecayAngle = val;
83  }
84  }
85 
86  void DecayGenFrame::SetDeltaPhiDecayPlane(double val){
87  while(val > acos(-1.)*2.) val -= acos(-1.)*2.;
88  while(val < 0.) val += acos(-1.)*2.;
89  m_DeltaPhiDecayPlane = val;
90  }
91 
92  void DecayGenFrame::ResetGenFrame(){
93  SetSpirit(false);
94  m_CosDecayAngle = -2.;
95  m_DeltaPhiDecayPlane = -2.;
96  }
97 
98  void DecayGenFrame::SetVariableMass(bool varymass) {
99  SetMind(false);
100  SetVariableMassMCMC(varymass);
101  }
102 
103  bool DecayGenFrame::InitializeGenAnalysis(){
104  if(!IsSoundBody())
105  return SetMind(false);
106 
107  double min_mass = GetMinimumMassMCMC();
108  if(m_Mass < min_mass && !IsVariableMassMCMC()){
109  m_Log << LogWarning;
110  m_Log << "Unable to initialize analysis: ";
111  m_Log << "decay frame mass (" << m_Mass << ") ";
112  m_Log << "is less than required child masses (";
113  m_Log << min_mass << ")" << LogEnd;
114  return SetMind(false);
115  }
116 
117  m_ChildIndexMCMC.clear();
118  m_ChildMassMCMC.clear();
119  m_ChildProbMCMC.clear();
120  double Mass = GetMass();
121  m_Log << LogInfo;
122  int N = GetNChildren();
123  for(int i = 0; i < N; i++){
124  double cmass = 0.;
125  double cprob = 1.;
126  GeneratorFrame& child = GetChildFrame(i);
127  if(child.IsVariableMassMCMC()){
128  child.GenerateMassMCMC(cmass, cprob, Mass);
129  SetMassMCMC(cmass, child);
130  m_ChildIndexMCMC.push_back(i);
131  m_ChildMassMCMC.push_back(cmass);
132  m_ChildProbMCMC.push_back(cprob);
133  } else {
134  cmass = child.GetMass();
135  }
136  Mass -= cmass;
137  }
138 
139  m_InterMassFracMCMC.clear();
140  m_InterMassFracMCMC.push_back(0.);
141  for(int i = 1; i < N-1; i++)
142  m_InterMassFracMCMC.push_back(GetRandom());
143  qsort((double*)(&m_InterMassFracMCMC[0])+1,N-2,sizeof(double),DoubleMax);
144  m_InterMassFracMCMC.push_back(1.);
145 
146  return SetMind(true);
147  }
148 
149  bool DecayGenFrame::IterateMCMC(){
150  int N = GetNChildren();
151 
152  std::vector<double> InterMassFrac;
153  InterMassFrac.push_back(0.);
154  for(int i = 1; i < N-1; i++)
155  InterMassFrac.push_back(GetRandom());
156  qsort((double*)(&InterMassFrac[0])+1,N-2,sizeof(double),DoubleMax);
157  InterMassFrac.push_back(1.);
158 
159  double probOld = GetProbMCMC();
160 
161  std::vector<double> InterMassFracOld = m_InterMassFracMCMC;
162  m_InterMassFracMCMC = InterMassFrac;
163 
164  double probNew = GetProbMCMC();
165 
166  if(probOld > 0.)
167  if(probNew/probOld < GetRandom())
168  m_InterMassFracMCMC = InterMassFracOld;
169 
170  int Nvar = m_ChildIndexMCMC.size();
171  for(int v = 0; v < Nvar; v++){
172  int index = m_ChildIndexMCMC[v];
173  GeneratorFrame& child = GetChildFrame(index);
174 
175  double massMax = GetMass();
176  for(int i = 0; i < N; i++)
177  if(i != index)
178  massMax -= GetChildFrame(i).GetMass();
179 
180  double ChildMass = 0.;
181  double ChildProb = 0.;
182  child.GenerateMassMCMC(ChildMass, ChildProb, massMax);
183  double probOld = child.GetProbMCMC(m_ChildMassMCMC[v]);
184  double probNew = child.GetProbMCMC(ChildMass);
185  probOld /= m_ChildProbMCMC[v];
186  probNew /= ChildProb;
187 
188  probOld *= GetProbMCMC();
189  SetMassMCMC(ChildMass, child);
190  probNew *= GetProbMCMC();
191 
192  if(probOld > 0){
193  if(probNew/probOld < GetRandom()){
194  SetMassMCMC(m_ChildMassMCMC[v], child);
195  } else {
196  m_ChildMassMCMC[v] = ChildMass;
197  m_ChildProbMCMC[v] = ChildProb;
198  }
199  }
200  }
201 
202  return SetMind(true);
203  }
204 
205  double DecayGenFrame::GetProbMCMC(double mass) const {
206  if(mass < 0.)
207  mass = GetMass();
208 
209  double SumChildMass = 0.;
210  int N = GetNChildren();
211  for(int i = 0; i < N; i++)
212  SumChildMass += GetChildFrame(i).GetMass();
213 
214  if(mass < SumChildMass)
215  return 0.;
216 
217  double ETOT = mass - SumChildMass;
218  std::vector<double> InterMass;
219  for(int i = 0; i < N; i++){
220  InterMass.push_back(m_InterMassFracMCMC[N-1-i]*ETOT + SumChildMass);
221  SumChildMass -= GetChildFrame(i).GetMass();
222  }
223 
224  double prob = 1.;
225  for(int i = 0; i < N-1; i++)
226  prob *= GetP(InterMass[i], InterMass[i+1], GetChildFrame(i).GetMass())/mass;
227 
228  prob /= mass*mass;
229 
230  return prob;
231  }
232 
233  void DecayGenFrame::GenerateMassMCMC(double& mass, double& prob,
234  double max) const {
235  int N = GetNChildren();
236  double SumMinChildMass = 0.;
237  double SumChildMass = 0.;
238  for(int i = 0; i < N; i++){
239  GeneratorFrame& child = GetChildFrame(i);
240  SumChildMass += child.GetMass();
241  if(!child.IsVariableMassMCMC())
242  SumMinChildMass += child.GetMass();
243  else
244  SumMinChildMass += child.GetMinimumMassMCMC();
245  }
246 
247  if(SumMinChildMass > max && max > 0){
248  mass = max;
249  prob = 0;
250  return;
251  }
252 
253  double T = SumChildMass;
254  double min = SumMinChildMass;
255  double SL = (T > 0 ? T/10. : 10.);
256  double SU = (T > 0 ? T : 100.);
257  SU = (max > 0 ? std::max(max/100.,T) : 100.);
258  double IL = SL*(1.-exp(-(T-min)/SL));
259  double IU = (max > 0 ? SU*(1.-exp(-(max-T)/SU)) : 1.);
260 
261  double R = GetRandom();
262  if(R > IL/(IL+IU)){
263  R = R*(IL+IU)/IU - IL/IU;
264  mass = T - log(1-R*IU/SU)*SU;
265  prob = exp(-(mass-T)/SU);
266  } else {
267  R = R*(IL+IU)/IL;
268  mass = T + log(1-R*IL/SL)*SL;
269  prob = exp((mass-T)/SL);
270  }
271 
272  }
273 
274  bool DecayGenFrame::GenerateFrame(){
275  if(!IsSoundMind()){
276  m_Log << LogWarning;
277  m_Log << "Unable to generate event for frame";
278  m_Log << LogEnd;
279  return SetSpirit(false);
280  }
281 
282  std::vector<double> ChildMasses;
283  double SumChildMass = 0.;
284  double cmass;
285  int N = GetNChildren();
286  for(int i = 0; i < N; i++){
287  cmass = GetChildFrame(i).GetMass();
288  ChildMasses.push_back(cmass);
289  SumChildMass += cmass;
290  }
291 
292  double ETOT = GetMass() - SumChildMass;
293  std::vector<double> InterMass;
294  for(int i = 0; i < N; i++){
295  InterMass.push_back(m_InterMassFracMCMC[N-1-i]*ETOT + SumChildMass);
296  SumChildMass -= ChildMasses[i];
297  }
298 
299 
300 
301  SetSpirit(true);
302 
303  TVector3 n_par = GetParentBoostVector();
304  TVector3 n_perp = GetParentFrame().GetDecayPlaneNormalVector(*this);
305 
306  if(n_par.Cross(n_perp).Mag() <= 0.){
307  n_par.SetXYZ(1.,0.,0.);
308  n_perp.SetXYZ(0.,1.,0.);
309  }
310 
311  std::vector<TLorentzVector> ChildVectors;
312  GenerateTwoBodyRecursive(InterMass, ChildMasses,
313  n_par, n_perp, ChildVectors);
314 
315  SetChildren(ChildVectors);
316 
317  return SetSpirit(true);
318  }
319 
320  void DecayGenFrame::GenerateTwoBodyRecursive(const std::vector<double>& M_p,
321  const std::vector<double>& M_c,
322  const TVector3& axis_par,
323  const TVector3& axis_perp,
324  std::vector<TLorentzVector>& P_c) {
325  TVector3 n_par = axis_par.Unit();
326  TVector3 n_perp = n_par.Cross(axis_perp.Cross(n_par)).Unit();
327 
328  int N_c = M_c.size();
329 
330  double m[2], Mp = M_p[0];
331  m[0] = M_c[0];
332  m[1] = M_p[1];
333  TVector3 V_c[2];
334 
335  double Pcm = GetP(Mp, m[0], m[1]);
336 
337  V_c[0] = Pcm*n_par;
338  V_c[1] = -Pcm*n_par;
339 
340  if(m_CosDecayAngle < -1.) m_CosDecayAngle = 1.-2.*GetRandom();
341  if(m_DeltaPhiDecayPlane < 0.) m_DeltaPhiDecayPlane = 2.*acos(-1.)*GetRandom();
342 
343  for(int i = 0; i < 2; i++) V_c[i].Rotate(-acos(m_CosDecayAngle),n_perp);
344  for(int i = 0; i < 2; i++) V_c[i].Rotate(-m_DeltaPhiDecayPlane,n_par);
345  m_CosDecayAngle = -2.;
346  m_DeltaPhiDecayPlane = -2.;
347 
348  TLorentzVector P_child[2];
349  for(int i = 0; i < 2; i++)
350  P_child[i].SetVectM(V_c[i], m[i]);
351  P_c.push_back(P_child[0]);
352 
353  if(N_c == 2){
354  P_c.push_back(P_child[1]);
355  return;
356  }
357 
358  // Recursively generate other two-body decays for N > 2
359  std::vector<double> M_pR;
360  std::vector<double> M_cR;
361  for(int i = 1; i < N_c; i++) M_pR.push_back(M_p[i]);
362  for(int i = 1; i < N_c; i++) M_cR.push_back(M_c[i]);
363  TVector3 boost = P_child[1].BoostVector();
364  std::vector<TLorentzVector> P_cR;
365  GenerateTwoBodyRecursive(M_pR, M_cR, boost, V_c[0].Cross(axis_par), P_cR);
366  for(int i = 0; i < N_c-1; i++) P_cR[i].Boost(boost);
367  for(int i = 0; i < N_c-1; i++) P_c.push_back(P_cR[i]);
368  }
369 
370  int DoubleMax(const void *a, const void *b){
371  Double_t aa = *((double*)a);
372  Double_t bb = *((double*)b);
373  if (aa > bb) return 1;
374  if (aa < bb) return -1;
375  return 0;
376  }
377 
378 }
virtual double GetMass() const
Get the mass of this frame.
virtual GeneratorFrame const & GetParentFrame() const
Returns the parent of this frame.
bool IsVariableMassMCMC() const
Frame is capable having a variable mass? (true/false)
virtual GeneratorFrame & GetChildFrame(int i=0) const
Get the frame of the i th child.