| 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> 48 #include <G4SystemOfUnits.hh>
49 #include "CexmcReimplementedGenbod.hh" 49 #include "CexmcReimplementedGenbod.hh"
50 #include "CexmcException.hh" 50 #include "CexmcException.hh"
51 51
52 52
53 namespace 53 namespace
54 { 54 {
55 G4int DoubleMax( const void * a, const v 55 G4int DoubleMax( const void * a, const void * b )
56 { 56 {
57 G4double aa( *( ( G4double * )a ) ); 57 G4double aa( *( ( G4double * )a ) );
58 G4double bb( *( ( G4double * )b ) ); 58 G4double bb( *( ( G4double * )b ) );
59 59
60 if ( aa > bb ) 60 if ( aa > bb )
61 return 1; 61 return 1;
62 62
63 if ( aa < bb ) 63 if ( aa < bb )
64 return -1; 64 return -1;
65 65
66 return 0; 66 return 0;
67 } 67 }
68 68
69 69
70 G4double PDK( G4double a, G4double b, G 70 G4double PDK( G4double a, G4double b, G4double c )
71 { 71 {
72 G4double x( ( a - b - c ) * ( a + b + 72 G4double x( ( a - b - c ) * ( a + b + c ) * ( a - b + c ) *
73 ( a + b - c ) ); 73 ( a + b - c ) );
74 x = std::sqrt( x ) / ( 2 * a ); 74 x = std::sqrt( x ) / ( 2 * a );
75 75
76 return x; 76 return x;
77 } 77 }
78 } 78 }
79 79
80 80
81 CexmcReimplementedGenbod::CexmcReimplementedGe 81 CexmcReimplementedGenbod::CexmcReimplementedGenbod() : maxWeight( 0. ),
82 nmbOfOutputParticles( 0 ) 82 nmbOfOutputParticles( 0 )
83 { 83 {
84 } 84 }
85 85
86 86
87 G4double CexmcReimplementedGenbod::Generate( 87 G4double CexmcReimplementedGenbod::Generate( void )
88 { 88 {
89 // Generate a random final state. 89 // Generate a random final state.
90 // The function returns the weigth of the 90 // The function returns the weigth of the current event.
91 // Note that Momentum, Energy units are Ge 91 // Note that Momentum, Energy units are Gev/C, GeV
92 92
93 G4double te_minus_tm( totalEnergy - total 93 G4double te_minus_tm( totalEnergy - totalMass );
94 G4double rno[ maxParticles ]; 94 G4double rno[ maxParticles ];
95 rno[ 0 ] = 0; 95 rno[ 0 ] = 0;
96 96
97 if ( nmbOfOutputParticles > 2 ) 97 if ( nmbOfOutputParticles > 2 )
98 { 98 {
99 for ( G4int i( 1 ); i < nmbOfOutputPa 99 for ( G4int i( 1 ); i < nmbOfOutputParticles - 1; ++i )
100 { 100 {
101 rno[ i ] = G4UniformRand(); 101 rno[ i ] = G4UniformRand();
102 } 102 }
103 qsort( rno + 1, nmbOfOutputParticles - 103 qsort( rno + 1, nmbOfOutputParticles - 2, sizeof( G4double ),
104 DoubleMax ); 104 DoubleMax );
105 } 105 }
106 rno[ nmbOfOutputParticles - 1 ] = 1; 106 rno[ nmbOfOutputParticles - 1 ] = 1;
107 107
108 G4double invMas[ maxParticles ]; 108 G4double invMas[ maxParticles ];
109 G4double sum( 0 ); 109 G4double sum( 0 );
110 110
111 for ( int i( 0 ); i < nmbOfOutputParticle 111 for ( int i( 0 ); i < nmbOfOutputParticles; ++i )
112 { 112 {
113 sum += outVec[ i ].mass / GeV; 113 sum += outVec[ i ].mass / GeV;
114 invMas[ i ] = rno[ i ] * te_minus_tm / 114 invMas[ i ] = rno[ i ] * te_minus_tm / GeV + sum;
115 } 115 }
116 116
117 // 117 //
118 //-----> compute the weight of the current 118 //-----> compute the weight of the current event
119 // 119 //
120 G4double wt( maxWeight ); 120 G4double wt( maxWeight );
121 G4double pd[ maxParticles ]; 121 G4double pd[ maxParticles ];
122 122
123 for ( int i( 0 ); i < nmbOfOutputParticle 123 for ( int i( 0 ); i < nmbOfOutputParticles - 1; ++i )
124 { 124 {
125 pd[ i ] = PDK( invMas[ i + 1 ], invMas 125 pd[ i ] = PDK( invMas[ i + 1 ], invMas[ i ],
126 outVec[ i + 1 ].mass / 126 outVec[ i + 1 ].mass / GeV );
127 wt *= pd[ i ]; 127 wt *= pd[ i ];
128 } 128 }
129 129
130 // 130 //
131 //-----> complete specification of event ( 131 //-----> complete specification of event (Raubold-Lynch method)
132 // 132 //
133 outVec[ 0 ].lVec->setPx( 0. ); 133 outVec[ 0 ].lVec->setPx( 0. );
134 outVec[ 0 ].lVec->setPy( pd[ 0 ] ); 134 outVec[ 0 ].lVec->setPy( pd[ 0 ] );
135 outVec[ 0 ].lVec->setPz( 0. ); 135 outVec[ 0 ].lVec->setPz( 0. );
136 outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] 136 outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] * pd[ 0 ] +
137 outVec[ 137 outVec[ 0 ].mass / GeV *
138 outVec[ 138 outVec[ 0 ].mass / GeV ) );
139 139
140 G4int i( 1 ); 140 G4int i( 1 );
141 141
142 while ( true ) 142 while ( true )
143 { 143 {
144 outVec[ i ].lVec->setPx( 0. ); 144 outVec[ i ].lVec->setPx( 0. );
145 outVec[ i ].lVec->setPy( -pd[ i - 1 ] 145 outVec[ i ].lVec->setPy( -pd[ i - 1 ] );
146 outVec[ i ].lVec->setPz( 0. ); 146 outVec[ i ].lVec->setPz( 0. );
147 outVec[ i ].lVec->setE( std::sqrt( pd[ 147 outVec[ i ].lVec->setE( std::sqrt( pd[ i - 1 ] * pd[ i - 1 ] +
148 out 148 outVec[ i ].mass / GeV *
149 out 149 outVec[ i ].mass / GeV ) );
150 150
151 G4double cZ( 2 * G4UniformRand() - 1 151 G4double cZ( 2 * G4UniformRand() - 1 );
152 G4double sZ( std::sqrt( 1 - cZ * cZ ) 152 G4double sZ( std::sqrt( 1 - cZ * cZ ) );
153 G4double angY( 2 * pi * G4UniformRand 153 G4double angY( 2 * pi * G4UniformRand() );
154 G4double cY( std::cos( angY ) ); 154 G4double cY( std::cos( angY ) );
155 G4double sY( std::sin( angY ) ); 155 G4double sY( std::sin( angY ) );
156 156
157 for ( int j( 0 ); j <= i; ++j ) 157 for ( int j( 0 ); j <= i; ++j )
158 { 158 {
159 G4LorentzVector * v( outVec[ j ]. 159 G4LorentzVector * v( outVec[ j ].lVec );
160 G4double x( v->px() ); 160 G4double x( v->px() );
161 G4double y( v->py() ); 161 G4double y( v->py() );
162 v->setPx( cZ * x - sZ * y ); 162 v->setPx( cZ * x - sZ * y );
163 v->setPy( sZ * x + cZ * y ); // 163 v->setPy( sZ * x + cZ * y ); // rotation around Z
164 x = v->px(); 164 x = v->px();
165 G4double z( v->pz() ); 165 G4double z( v->pz() );
166 v->setPx( cY * x - sY * z ); 166 v->setPx( cY * x - sY * z );
167 v->setPz( sY * x + cY * z ); // 167 v->setPz( sY * x + cY * z ); // rotation around Y
168 } 168 }
169 169
170 if ( i == nmbOfOutputParticles - 1 ) 170 if ( i == nmbOfOutputParticles - 1 )
171 break; 171 break;
172 172
173 G4double beta( pd[ i ] / std::sqrt( p 173 G4double beta( pd[ i ] / std::sqrt( pd[ i ] * pd[ i ] +
174 i 174 invMas[ i ] * invMas[ i ] ) );
175 for ( int j( 0 ); j <= i; ++j ) 175 for ( int j( 0 ); j <= i; ++j )
176 outVec[ j ].lVec->boost( 0, beta, 176 outVec[ j ].lVec->boost( 0, beta, 0 );
177 177
178 ++i; 178 ++i;
179 } 179 }
180 180
181 for ( int j( 0 ); j < nmbOfOutputParticle 181 for ( int j( 0 ); j < nmbOfOutputParticles; ++j )
182 *outVec[ j ].lVec *= GeV; 182 *outVec[ j ].lVec *= GeV;
183 183
184 // 184 //
185 //---> return the weigth of event 185 //---> return the weigth of event
186 // 186 //
187 return wt; 187 return wt;
188 } 188 }
189 189
190 190
191 void CexmcReimplementedGenbod::ParticleChange 191 void CexmcReimplementedGenbod::ParticleChangeHook( void )
192 { 192 {
193 nmbOfOutputParticles = outVec.size(); 193 nmbOfOutputParticles = outVec.size();
194 194
195 if ( nmbOfOutputParticles < 2 || nmbOfOutp 195 if ( nmbOfOutputParticles < 2 || nmbOfOutputParticles > maxParticles )
196 throw CexmcException( CexmcKinematicsE 196 throw CexmcException( CexmcKinematicsException );
197 197
198 SetMaxWeight(); 198 SetMaxWeight();
199 } 199 }
200 200
201 201
202 void CexmcReimplementedGenbod::FermiEnergyDep 202 void CexmcReimplementedGenbod::FermiEnergyDepStatusChangeHook( void )
203 { 203 {
204 SetMaxWeight(); 204 SetMaxWeight();
205 } 205 }
206 206
207 207
208 void CexmcReimplementedGenbod::SetMaxWeight( 208 void CexmcReimplementedGenbod::SetMaxWeight( void )
209 { 209 {
210 G4double te_minus_tm( totalEnergy - total 210 G4double te_minus_tm( totalEnergy - totalMass );
211 211
212 if ( fermiEnergyDepIsOn ) 212 if ( fermiEnergyDepIsOn )
213 { 213 {
214 // ffq[] = pi * (2*pi)**(N-2) / (N-2)! 214 // ffq[] = pi * (2*pi)**(N-2) / (N-2)!
215 G4double ffq[] = { 0 215 G4double ffq[] = { 0
216 ,3.141592, 19.73921, 62.0 216 ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
217 ,256.3704, 268.4705, 240. 217 ,256.3704, 268.4705, 240.9780, 189.2637
218 ,132.1308, 83.0202, 47. 218 ,132.1308, 83.0202, 47.4210, 24.8295
219 ,12.0006, 5.3858, 2.2 219 ,12.0006, 5.3858, 2.2560, 0.8859 };
220 maxWeight = 220 maxWeight =
221 std::pow( te_minus_tm / GeV, n 221 std::pow( te_minus_tm / GeV, nmbOfOutputParticles - 2 ) *
222 ffq[ nmbOfOutputParticles - 1 222 ffq[ nmbOfOutputParticles - 1 ] / ( totalEnergy / GeV );
223 } 223 }
224 else 224 else
225 { 225 {
226 G4double emmax( ( te_minus_tm + outVe 226 G4double emmax( ( te_minus_tm + outVec[ 0 ].mass ) / GeV );
227 G4double emmin( 0. ); 227 G4double emmin( 0. );
228 G4double wtmax( 1. ); 228 G4double wtmax( 1. );
229 229
230 for ( G4int i( 1 ); i < nmbOfOutputPa 230 for ( G4int i( 1 ); i < nmbOfOutputParticles; ++i )
231 { 231 {
232 emmin += outVec[ i - 1 ].mass / Ge 232 emmin += outVec[ i - 1 ].mass / GeV;
233 emmax += outVec[ i ].mass / GeV; 233 emmax += outVec[ i ].mass / GeV;
234 wtmax *= PDK( emmax, emmin, outVec 234 wtmax *= PDK( emmax, emmin, outVec[ i ].mass / GeV );
235 } 235 }
236 maxWeight = 1 / wtmax; 236 maxWeight = 1 / wtmax;
237 } 237 }
238 } 238 }
239 239
240 240