Bug Summary

File:libraries/AMPTOOLS_MCGEN/NBodyPhaseSpaceFactory.cc
Location:line 91, column 47
Description:The left operand of '*' is a garbage value

Annotated Source Code

1/*
2 Generates N-body phase space decays using the Raubold Lynch
3 method (F.James CERN 68-15 (1968).
4 Borrows liberally from ROOT class TGenPhaseSpace as well as AcquRoot
5 class TMCGenerator (J.R.M.Annand -- http://nuclear.gla.ac.uk/~acqusys/doc).
6
7 -- C.Tarbert
8*/
9
10#include <vector>
11#include <cstdlib>
12#include <cassert>
13#include <iostream>
14#include <algorithm>
15
16#include "TLorentzVector.h"
17#include "TLorentzRotation.h"
18#include "TMath.h"
19
20#include "AMPTOOLS_MCGEN/NBodyPhaseSpaceFactory.h"
21
22const double NBodyPhaseSpaceFactory::kPi = 3.14159;
23
24NBodyPhaseSpaceFactory::NBodyPhaseSpaceFactory( double parentMass, const vector<double>& childMass ) :
25 m_parentMass( parentMass ),
26 m_childMass( childMass )
27{
28 m_Nd = (int)childMass.size();
29}
30
31vector<TLorentzVector>
32NBodyPhaseSpaceFactory::generateDecay(bool uniformWeights) {
33
34
35 vector<TLorentzVector> child( m_Nd );
36
37 double Tcm = m_parentMass;
38 int n, m;
39 for( n=0; n<m_Nd; n++ ){ Tcm -= m_childMass[n]; }
1
Loop condition is true. Entering loop body
→
2
←
Loop condition is false. Execution continues on line 40
→
40 assert( Tcm > 0. )((Tcm > 0.) ? static_cast<void> (0) : __assert_fail (
"Tcm > 0.", "libraries/AMPTOOLS_MCGEN/NBodyPhaseSpaceFactory.cc"
, 40, __PRETTY_FUNCTION__));
41
42 double emmax = Tcm + m_childMass[0];
43 double emmin = 0;
44 double wt = 1;
45 for (n=1; n<m_Nd; n++) {
3
←
Loop condition is false. Execution continues on line 50
→
46 emmin += m_childMass[n-1];
47 emmax += m_childMass[n];
48 wt *= pdk(emmax, emmin, m_childMass[n]);
49 }
50 double WtMax = 1/wt; // max weight for gating events
51
52 double rnd[m_Nd];
53 double invMas[m_Nd];
54 double pd[m_Nd];
55 int irnd[m_Nd];
56 rnd[0] = 0.0; rnd[m_Nd-1] = 1.0;
57 do{
12
←
Loop condition is false. Exiting loop
→
58 switch( m_Nd ){
4
←
Control jumps to the 'default' case at line 59
→
59 default:
60 for( n=1; n<m_Nd-1; n++){ rnd[n] = random( 0., 1. ); }
5
←
Loop condition is false. Execution continues on line 61
→
61 for( n=0; n<m_Nd; n++ ){ irnd[n] = n; }
6
←
Loop condition is true. Entering loop body
→
7
←
Loop condition is false. Execution continues on line 62
→
62 sort( irnd, irnd+m_Nd, CompareAsc<const double*>(rnd) ); // sort random numbers ascending
63 rnd[m_Nd-1] = 1.0;
64 break;
8
←
 Execution continues on line 73
→
65 case 3: // 3 decay particles
66 rnd[1] = random(0.0,1.0);
67 irnd[0] = 0; irnd[1] = 1; irnd[2] = 2;
68 break;
69 case 2:
70 irnd[0] = 0; irnd[1] = 1;
71 break;
72 }
73 wt = 0.0;
74 for (n=0; n<m_Nd; n++) {
9
←
Loop condition is true. Entering loop body
→
10
←
Loop condition is false. Execution continues on line 78
→
75 wt += m_childMass[n];
76 invMas[n] = rnd[irnd[n]]*Tcm + wt;
77 }
78 wt = WtMax;
79 for (n=0; n<m_Nd-1; n++) {
11
←
Loop condition is false. Execution continues on line 83
→
80 pd[n] = pdk(invMas[n+1],invMas[n],m_childMass[n+1]);
81 wt *= pd[n];
82 }
83 }while( uniformWeights && (random(0.0,WtMax) > wt) );
84
85 if(uniformWeights) m_lastWt = 1.0;
13
←
Taking false branch
→
86 else m_lastWt = wt;
87
88 //
89 // Specification of 4-momenta (Raubold-Lynch method)
90 //
91 child[0].SetPxPyPzE(0, pd[0], 0 , sqrt(pd[0]*pd[0]+m_childMass[0]*m_childMass[0]) );
14
←
The left operand of '*' is a garbage value
92 for(n=1;;){
93 child[n].SetPxPyPzE(0, -pd[n-1], 0 ,
94 sqrt(pd[n-1]*pd[n-1]+m_childMass[n]*m_childMass[n]) );
95
96 double cosZ = random(-1.,1.);
97 double angY = random(0.0, 2.*kPi);
98 for (m=0; m<=n; m++) {
99 child[m].RotateZ( acos(cosZ) );
100 child[m].RotateY( angY );
101 }
102 if( n == m_Nd-1 ) break;
103 double beta = pd[n] / sqrt(pd[n]*pd[n] + invMas[n]*invMas[n]);
104 for (m=0; m<=n; m++) child[m].Boost(0,beta,0);
105 n++;
106 }
107
108 return child;
109}
110
111double
112NBodyPhaseSpaceFactory::pdk( double a, double b, double c ) const {
113
114 double x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c);
115 x = sqrt(x)/(2*a);
116 return x;
117
118}
119
120double
121NBodyPhaseSpaceFactory::random( double low, double hi ) const {
122
123 return( ( hi - low ) * gRandom->Uniform() + low );
124}