1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
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 | |
22 | const double NBodyPhaseSpaceFactory::kPi = 3.14159; |
23 | |
24 | NBodyPhaseSpaceFactory::NBodyPhaseSpaceFactory( double parentMass, const vector<double>& childMass ) : |
25 | m_parentMass( parentMass ), |
26 | m_childMass( childMass ) |
27 | { |
28 | m_Nd = (int)childMass.size(); |
29 | } |
30 | |
31 | vector<TLorentzVector> |
32 | NBodyPhaseSpaceFactory::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; |
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) ); |
63 | rnd[m_Nd-1] = 1.0; |
64 | break; |
| 8 | | Execution continues on line 73 | |
|
65 | case 3: |
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; |
| |
86 | else m_lastWt = wt; |
87 | |
88 | |
89 | |
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 | |
111 | double |
112 | NBodyPhaseSpaceFactory::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 | |
120 | double |
121 | NBodyPhaseSpaceFactory::random( double low, double hi ) const { |
122 | |
123 | return( ( hi - low ) * gRandom->Uniform() + low ); |
124 | } |