Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/ChargeExchangeMC/src/CexmcReimplementedGenbod.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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