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