| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 * =========================================== 27 * =============================================================================
28 * 28 *
29 * Filename: CexmcReimplementedGenbod.c 29 * Filename: CexmcReimplementedGenbod.cc
30 * 30 *
31 * Description: reimplemented GENBOD 31 * Description: reimplemented GENBOD
32 * (mostly adopted from ROOT 32 * (mostly adopted from ROOT TGenPhaseSpace)
33 * 33 *
34 * Version: 1.0 34 * Version: 1.0
35 * Created: 08.09.2010 18:52:39 35 * Created: 08.09.2010 18:52:39
36 * Revision: none 36 * Revision: none
37 * Compiler: gcc 37 * Compiler: gcc
38 * 38 *
39 * Author: Alexey Radkov (), 39 * Author: Alexey Radkov (),
40 * Company: PNPI 40 * Company: PNPI
41 * 41 *
42 * =========================================== 42 * =============================================================================
43 */ 43 */
44 44
45 #include <cmath> 45 #include <cmath>
46 #include <Randomize.hh> 46 #include <Randomize.hh>
47 #include <G4PhysicalConstants.hh> 47 #include <G4PhysicalConstants.hh>
48 #include <G4SystemOfUnits.hh> <<
49 #include "CexmcReimplementedGenbod.hh" 48 #include "CexmcReimplementedGenbod.hh"
50 #include "CexmcException.hh" 49 #include "CexmcException.hh"
51 50
52 51
53 namespace 52 namespace
54 { 53 {
55 G4int DoubleMax( const void * a, const v 54 G4int DoubleMax( const void * a, const void * b )
56 { 55 {
57 G4double aa( *( ( G4double * )a ) ); 56 G4double aa( *( ( G4double * )a ) );
58 G4double bb( *( ( G4double * )b ) ); 57 G4double bb( *( ( G4double * )b ) );
59 58
60 if ( aa > bb ) 59 if ( aa > bb )
61 return 1; 60 return 1;
62 61
63 if ( aa < bb ) 62 if ( aa < bb )
64 return -1; 63 return -1;
65 64
66 return 0; 65 return 0;
67 } 66 }
68 67
69 68
70 G4double PDK( G4double a, G4double b, G 69 G4double PDK( G4double a, G4double b, G4double c )
71 { 70 {
72 G4double x( ( a - b - c ) * ( a + b + 71 G4double x( ( a - b - c ) * ( a + b + c ) * ( a - b + c ) *
73 ( a + b - c ) ); 72 ( a + b - c ) );
74 x = std::sqrt( x ) / ( 2 * a ); 73 x = std::sqrt( x ) / ( 2 * a );
75 74
76 return x; 75 return x;
77 } 76 }
78 } 77 }
79 78
80 79
81 CexmcReimplementedGenbod::CexmcReimplementedGe 80 CexmcReimplementedGenbod::CexmcReimplementedGenbod() : maxWeight( 0. ),
82 nmbOfOutputParticles( 0 ) 81 nmbOfOutputParticles( 0 )
83 { 82 {
84 } 83 }
85 84
86 85
87 G4double CexmcReimplementedGenbod::Generate( 86 G4double CexmcReimplementedGenbod::Generate( void )
88 { 87 {
89 // Generate a random final state. 88 // Generate a random final state.
90 // The function returns the weigth of the 89 // The function returns the weigth of the current event.
91 // Note that Momentum, Energy units are Ge 90 // Note that Momentum, Energy units are Gev/C, GeV
92 91
93 G4double te_minus_tm( totalEnergy - total 92 G4double te_minus_tm( totalEnergy - totalMass );
94 G4double rno[ maxParticles ]; 93 G4double rno[ maxParticles ];
95 rno[ 0 ] = 0; 94 rno[ 0 ] = 0;
96 95
97 if ( nmbOfOutputParticles > 2 ) 96 if ( nmbOfOutputParticles > 2 )
98 { 97 {
99 for ( G4int i( 1 ); i < nmbOfOutputPa 98 for ( G4int i( 1 ); i < nmbOfOutputParticles - 1; ++i )
100 { 99 {
101 rno[ i ] = G4UniformRand(); 100 rno[ i ] = G4UniformRand();
102 } 101 }
103 qsort( rno + 1, nmbOfOutputParticles - 102 qsort( rno + 1, nmbOfOutputParticles - 2, sizeof( G4double ),
104 DoubleMax ); 103 DoubleMax );
105 } 104 }
106 rno[ nmbOfOutputParticles - 1 ] = 1; 105 rno[ nmbOfOutputParticles - 1 ] = 1;
107 106
108 G4double invMas[ maxParticles ]; 107 G4double invMas[ maxParticles ];
109 G4double sum( 0 ); 108 G4double sum( 0 );
110 109
111 for ( int i( 0 ); i < nmbOfOutputParticle 110 for ( int i( 0 ); i < nmbOfOutputParticles; ++i )
112 { 111 {
113 sum += outVec[ i ].mass / GeV; 112 sum += outVec[ i ].mass / GeV;
114 invMas[ i ] = rno[ i ] * te_minus_tm / 113 invMas[ i ] = rno[ i ] * te_minus_tm / GeV + sum;
115 } 114 }
116 115
117 // 116 //
118 //-----> compute the weight of the current 117 //-----> compute the weight of the current event
119 // 118 //
120 G4double wt( maxWeight ); 119 G4double wt( maxWeight );
121 G4double pd[ maxParticles ]; 120 G4double pd[ maxParticles ];
122 121
123 for ( int i( 0 ); i < nmbOfOutputParticle 122 for ( int i( 0 ); i < nmbOfOutputParticles - 1; ++i )
124 { 123 {
125 pd[ i ] = PDK( invMas[ i + 1 ], invMas 124 pd[ i ] = PDK( invMas[ i + 1 ], invMas[ i ],
126 outVec[ i + 1 ].mass / 125 outVec[ i + 1 ].mass / GeV );
127 wt *= pd[ i ]; 126 wt *= pd[ i ];
128 } 127 }
129 128
130 // 129 //
131 //-----> complete specification of event ( 130 //-----> complete specification of event (Raubold-Lynch method)
132 // 131 //
133 outVec[ 0 ].lVec->setPx( 0. ); 132 outVec[ 0 ].lVec->setPx( 0. );
134 outVec[ 0 ].lVec->setPy( pd[ 0 ] ); 133 outVec[ 0 ].lVec->setPy( pd[ 0 ] );
135 outVec[ 0 ].lVec->setPz( 0. ); 134 outVec[ 0 ].lVec->setPz( 0. );
136 outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] 135 outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] * pd[ 0 ] +
137 outVec[ 136 outVec[ 0 ].mass / GeV *
138 outVec[ 137 outVec[ 0 ].mass / GeV ) );
139 138
140 G4int i( 1 ); 139 G4int i( 1 );
141 140
142 while ( true ) 141 while ( true )
143 { 142 {
144 outVec[ i ].lVec->setPx( 0. ); 143 outVec[ i ].lVec->setPx( 0. );
145 outVec[ i ].lVec->setPy( -pd[ i - 1 ] 144 outVec[ i ].lVec->setPy( -pd[ i - 1 ] );
146 outVec[ i ].lVec->setPz( 0. ); 145 outVec[ i ].lVec->setPz( 0. );
147 outVec[ i ].lVec->setE( std::sqrt( pd[ 146 outVec[ i ].lVec->setE( std::sqrt( pd[ i - 1 ] * pd[ i - 1 ] +
148 out 147 outVec[ i ].mass / GeV *
149 out 148 outVec[ i ].mass / GeV ) );
150 149
151 G4double cZ( 2 * G4UniformRand() - 1 150 G4double cZ( 2 * G4UniformRand() - 1 );
152 G4double sZ( std::sqrt( 1 - cZ * cZ ) 151 G4double sZ( std::sqrt( 1 - cZ * cZ ) );
153 G4double angY( 2 * pi * G4UniformRand 152 G4double angY( 2 * pi * G4UniformRand() );
154 G4double cY( std::cos( angY ) ); 153 G4double cY( std::cos( angY ) );
155 G4double sY( std::sin( angY ) ); 154 G4double sY( std::sin( angY ) );
156 155
157 for ( int j( 0 ); j <= i; ++j ) 156 for ( int j( 0 ); j <= i; ++j )
158 { 157 {
159 G4LorentzVector * v( outVec[ j ]. 158 G4LorentzVector * v( outVec[ j ].lVec );
160 G4double x( v->px() ); 159 G4double x( v->px() );
161 G4double y( v->py() ); 160 G4double y( v->py() );
162 v->setPx( cZ * x - sZ * y ); 161 v->setPx( cZ * x - sZ * y );
163 v->setPy( sZ * x + cZ * y ); // 162 v->setPy( sZ * x + cZ * y ); // rotation around Z
164 x = v->px(); 163 x = v->px();
165 G4double z( v->pz() ); 164 G4double z( v->pz() );
166 v->setPx( cY * x - sY * z ); 165 v->setPx( cY * x - sY * z );
167 v->setPz( sY * x + cY * z ); // 166 v->setPz( sY * x + cY * z ); // rotation around Y
168 } 167 }
169 168
170 if ( i == nmbOfOutputParticles - 1 ) 169 if ( i == nmbOfOutputParticles - 1 )
171 break; 170 break;
172 171
173 G4double beta( pd[ i ] / std::sqrt( p 172 G4double beta( pd[ i ] / std::sqrt( pd[ i ] * pd[ i ] +
174 i 173 invMas[ i ] * invMas[ i ] ) );
175 for ( int j( 0 ); j <= i; ++j ) 174 for ( int j( 0 ); j <= i; ++j )
176 outVec[ j ].lVec->boost( 0, beta, 175 outVec[ j ].lVec->boost( 0, beta, 0 );
177 176
178 ++i; 177 ++i;
179 } 178 }
180 179
181 for ( int j( 0 ); j < nmbOfOutputParticle << 180 for ( int i( 0 ); i < nmbOfOutputParticles; ++i )
182 *outVec[ j ].lVec *= GeV; << 181 *outVec[ i ].lVec *= GeV;
183 182
184 // 183 //
185 //---> return the weigth of event 184 //---> return the weigth of event
186 // 185 //
187 return wt; 186 return wt;
188 } 187 }
189 188
190 189
191 void CexmcReimplementedGenbod::ParticleChange 190 void CexmcReimplementedGenbod::ParticleChangeHook( void )
192 { 191 {
193 nmbOfOutputParticles = outVec.size(); 192 nmbOfOutputParticles = outVec.size();
194 193
195 if ( nmbOfOutputParticles < 2 || nmbOfOutp 194 if ( nmbOfOutputParticles < 2 || nmbOfOutputParticles > maxParticles )
196 throw CexmcException( CexmcKinematicsE 195 throw CexmcException( CexmcKinematicsException );
197 196
198 SetMaxWeight(); 197 SetMaxWeight();
199 } 198 }
200 199
201 200
202 void CexmcReimplementedGenbod::FermiEnergyDep 201 void CexmcReimplementedGenbod::FermiEnergyDepStatusChangeHook( void )
203 { 202 {
204 SetMaxWeight(); 203 SetMaxWeight();
205 } 204 }
206 205
207 206
208 void CexmcReimplementedGenbod::SetMaxWeight( 207 void CexmcReimplementedGenbod::SetMaxWeight( void )
209 { 208 {
210 G4double te_minus_tm( totalEnergy - total 209 G4double te_minus_tm( totalEnergy - totalMass );
211 210
212 if ( fermiEnergyDepIsOn ) 211 if ( fermiEnergyDepIsOn )
213 { 212 {
214 // ffq[] = pi * (2*pi)**(N-2) / (N-2)! 213 // ffq[] = pi * (2*pi)**(N-2) / (N-2)!
215 G4double ffq[] = { 0 214 G4double ffq[] = { 0
216 ,3.141592, 19.73921, 62.0 215 ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
217 ,256.3704, 268.4705, 240. 216 ,256.3704, 268.4705, 240.9780, 189.2637
218 ,132.1308, 83.0202, 47. 217 ,132.1308, 83.0202, 47.4210, 24.8295
219 ,12.0006, 5.3858, 2.2 218 ,12.0006, 5.3858, 2.2560, 0.8859 };
220 maxWeight = 219 maxWeight =
221 std::pow( te_minus_tm / GeV, n 220 std::pow( te_minus_tm / GeV, nmbOfOutputParticles - 2 ) *
222 ffq[ nmbOfOutputParticles - 1 221 ffq[ nmbOfOutputParticles - 1 ] / ( totalEnergy / GeV );
223 } 222 }
224 else 223 else
225 { 224 {
226 G4double emmax( ( te_minus_tm + outVe 225 G4double emmax( ( te_minus_tm + outVec[ 0 ].mass ) / GeV );
227 G4double emmin( 0. ); 226 G4double emmin( 0. );
228 G4double wtmax( 1. ); 227 G4double wtmax( 1. );
229 228
230 for ( G4int i( 1 ); i < nmbOfOutputPa 229 for ( G4int i( 1 ); i < nmbOfOutputParticles; ++i )
231 { 230 {
232 emmin += outVec[ i - 1 ].mass / Ge 231 emmin += outVec[ i - 1 ].mass / GeV;
233 emmax += outVec[ i ].mass / GeV; 232 emmax += outVec[ i ].mass / GeV;
234 wtmax *= PDK( emmax, emmin, outVec 233 wtmax *= PDK( emmax, emmin, outVec[ i ].mass / GeV );
235 } 234 }
236 maxWeight = 1 / wtmax; 235 maxWeight = 1 / wtmax;
237 } 236 }
238 } 237 }
239 238
240 239