33 namespace RestFrames {
40 const std::string& stitle)
43 m_CosDecayAngle = -2.;
44 m_DeltaPhiDecayPlane = -2.;
49 DecayGenFrame::~DecayGenFrame() {}
51 bool DecayGenFrame::IsSoundBody()
const{
52 if(RFBase::IsSoundBody())
54 if(!RestFrame::IsSoundBody())
58 return SetBody(
false);
66 m_Log <<
"Unable to set mass to negative value ";
67 m_Log << val <<
". Setting to zero." << LogEnd;
78 m_Log <<
"CosDecay angle must be in [-1, 1]: ";
79 m_Log << val <<
". Setting to random." << LogEnd;
80 m_CosDecayAngle = -2.;
82 m_CosDecayAngle = val;
87 while(val > acos(-1.)*2.) val -= acos(-1.)*2.;
88 while(val < 0.) val += acos(-1.)*2.;
89 m_DeltaPhiDecayPlane = val;
92 void DecayGenFrame::ResetGenFrame(){
94 m_CosDecayAngle = -2.;
95 m_DeltaPhiDecayPlane = -2.;
100 SetVariableMassMCMC(varymass);
103 bool DecayGenFrame::InitializeGenAnalysis(){
105 return SetMind(
false);
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);
117 m_ChildIndexMCMC.clear();
118 m_ChildMassMCMC.clear();
119 m_ChildProbMCMC.clear();
123 for(
int i = 0; i < N; 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);
134 cmass = child.GetMass();
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.);
146 return SetMind(
true);
149 bool DecayGenFrame::IterateMCMC(){
152 std::vector<double> InterMassFrac;
153 InterMassFrac.push_back(0.);
154 for(
int i = 1; i < N-1; i++)
156 qsort((
double*)(&InterMassFrac[0])+1,N-2,
sizeof(
double),DoubleMax);
157 InterMassFrac.push_back(1.);
161 std::vector<double> InterMassFracOld = m_InterMassFracMCMC;
162 m_InterMassFracMCMC = InterMassFrac;
168 m_InterMassFracMCMC = InterMassFracOld;
170 int Nvar = m_ChildIndexMCMC.size();
171 for(
int v = 0; v < Nvar; v++){
172 int index = m_ChildIndexMCMC[v];
176 for(
int i = 0; i < N; i++)
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;
189 SetMassMCMC(ChildMass, child);
194 SetMassMCMC(m_ChildMassMCMC[v], child);
196 m_ChildMassMCMC[v] = ChildMass;
197 m_ChildProbMCMC[v] = ChildProb;
202 return SetMind(
true);
209 double SumChildMass = 0.;
211 for(
int i = 0; i < N; i++)
214 if(mass < SumChildMass)
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);
225 for(
int i = 0; i < N-1; i++)
236 double SumMinChildMass = 0.;
237 double SumChildMass = 0.;
238 for(
int i = 0; i < N; i++){
240 SumChildMass += child.
GetMass();
242 SumMinChildMass += child.
GetMass();
247 if(SumMinChildMass > max && max > 0){
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.);
263 R = R*(IL+IU)/IU - IL/IU;
264 mass = T - log(1-R*IU/SU)*SU;
265 prob = exp(-(mass-T)/SU);
268 mass = T + log(1-R*IL/SL)*SL;
269 prob = exp((mass-T)/SL);
274 bool DecayGenFrame::GenerateFrame(){
277 m_Log <<
"Unable to generate event for frame";
279 return SetSpirit(
false);
282 std::vector<double> ChildMasses;
283 double SumChildMass = 0.;
286 for(
int i = 0; i < N; i++){
288 ChildMasses.push_back(cmass);
289 SumChildMass += cmass;
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];
303 TVector3 n_par = GetParentBoostVector();
306 if(n_par.Cross(n_perp).Mag() <= 0.){
307 n_par.SetXYZ(1.,0.,0.);
308 n_perp.SetXYZ(0.,1.,0.);
311 std::vector<TLorentzVector> ChildVectors;
312 GenerateTwoBodyRecursive(InterMass, ChildMasses,
313 n_par, n_perp, ChildVectors);
315 SetChildren(ChildVectors);
317 return SetSpirit(
true);
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();
328 int N_c = M_c.size();
330 double m[2], Mp = M_p[0];
335 double Pcm = GetP(Mp, m[0], m[1]);
340 if(m_CosDecayAngle < -1.) m_CosDecayAngle = 1.-2.*
GetRandom();
341 if(m_DeltaPhiDecayPlane < 0.) m_DeltaPhiDecayPlane = 2.*acos(-1.)*
GetRandom();
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.;
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]);
354 P_c.push_back(P_child[1]);
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]);
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;