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