ChargeExchangeMC/src/CexmcReimplementedGenbod.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /*
0027  * =============================================================================
0028  *
0029  *       Filename:  CexmcReimplementedGenbod.cc
0030  *
0031  *    Description:  reimplemented GENBOD
0032  *                  (mostly adopted from ROOT TGenPhaseSpace)
0033  *
0034  *        Version:  1.0
0035  *        Created:  08.09.2010 18:52:39
0036  *       Revision:  none
0037  *       Compiler:  gcc
0038  *
0039  *         Author:  Alexey Radkov (),
0040  *        Company:  PNPI
0041  *
0042  * =============================================================================
0043  */
0044
0045 #include <cmath>
0046 #include <Randomize.hh>
0047 #include <G4PhysicalConstants.hh>
0048 #include <G4SystemOfUnits.hh>
0049 #include "CexmcReimplementedGenbod.hh"
0050 #include "CexmcException.hh"
0051
0052
0053 namespace
0054 {
0055     G4int  DoubleMax( const void *  a, const void *  b )
0056     {
0057        G4double  aa( *( ( G4double * )a ) );
0058        G4double  bb( *( ( G4double * )b ) );
0059
0060        if ( aa > bb )
0061            return  1;
0062
0063        if ( aa < bb )
0064            return -1;
0065
0066        return 0;
0067     }
0068
0069
0070     G4double  PDK( G4double  a, G4double  b, G4double  c )
0071     {
0072         G4double  x( ( a - b - c ) * ( a + b + c ) * ( a - b + c ) *
0073                      ( a + b - c ) );
0074         x = std::sqrt( x ) / ( 2 * a );
0075
0076         return x;
0077     }
0078 }
0079
0080
0081 CexmcReimplementedGenbod::CexmcReimplementedGenbod() : maxWeight( 0. ),
0082     nmbOfOutputParticles( 0 )
0083 {
0084 }
0085
0086
0087 G4double  CexmcReimplementedGenbod::Generate( void )
0088 {
0089     // Generate a random final state.
0090     // The function returns the weigth of the current event.
0091     // Note that Momentum, Energy units are Gev/C, GeV
0092
0093     G4double  te_minus_tm( totalEnergy - totalMass );
0094     G4double  rno[ maxParticles ];
0095     rno[ 0 ] = 0;
0096
0097     if ( nmbOfOutputParticles > 2 )
0098     {
0099         for ( G4int  i( 1 ); i < nmbOfOutputParticles - 1; ++i )
0100         {
0101             rno[ i ] = G4UniformRand();
0102         }
0103         qsort( rno + 1, nmbOfOutputParticles - 2, sizeof( G4double ),
0104                DoubleMax );
0105     }
0106     rno[ nmbOfOutputParticles - 1 ] = 1;
0107
0108     G4double  invMas[ maxParticles ];
0109     G4double  sum( 0 );
0110
0111     for ( int  i( 0 ); i < nmbOfOutputParticles; ++i )
0112     {
0113         sum += outVec[ i ].mass / GeV;
0114         invMas[ i ] = rno[ i ] * te_minus_tm / GeV + sum;
0115     }
0116
0117     //
0118     //-----> compute the weight of the current event
0119     //
0120     G4double  wt( maxWeight );
0121     G4double  pd[ maxParticles ];
0122
0123     for ( int  i( 0 ); i < nmbOfOutputParticles - 1; ++i )
0124     {
0125         pd[ i ] = PDK( invMas[ i + 1 ], invMas[ i ],
0126                        outVec[ i + 1 ].mass / GeV );
0127         wt *= pd[ i ];
0128     }
0129
0130     //
0131     //-----> complete specification of event (Raubold-Lynch method)
0132     //
0133     outVec[ 0 ].lVec->setPx( 0. );
0134     outVec[ 0 ].lVec->setPy( pd[ 0 ] );
0135     outVec[ 0 ].lVec->setPz( 0. );
0136     outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] * pd[ 0 ] +
0137                                        outVec[ 0 ].mass / GeV *
0138                                        outVec[ 0 ].mass / GeV ) );
0139
0140     G4int  i( 1 );
0141
0142     while ( true )
0143     {
0144         outVec[ i ].lVec->setPx( 0. );
0145         outVec[ i ].lVec->setPy( -pd[ i - 1 ] );
0146         outVec[ i ].lVec->setPz( 0. );
0147         outVec[ i ].lVec->setE( std::sqrt( pd[ i - 1 ] * pd[ i - 1 ] +
0148                                            outVec[ i ].mass / GeV *
0149                                            outVec[ i ].mass / GeV ) );
0150
0151         G4double  cZ( 2 * G4UniformRand() - 1 );
0152         G4double  sZ( std::sqrt( 1 - cZ * cZ ) );
0153         G4double  angY( 2 * pi * G4UniformRand() );
0154         G4double  cY( std::cos( angY ) );
0155         G4double  sY( std::sin( angY ) );
0156
0157         for ( int  j( 0 ); j <= i; ++j )
0158         {
0159             G4LorentzVector *  v( outVec[ j ].lVec );
0160             G4double           x( v->px() );
0161             G4double           y( v->py() );
0162             v->setPx( cZ * x - sZ * y );
0163             v->setPy( sZ * x + cZ * y );   // rotation around Z
0164             x = v->px();
0165             G4double           z( v->pz() );
0166             v->setPx( cY * x - sY * z );
0167             v->setPz( sY * x + cY * z );   // rotation around Y
0168         }
0169
0170         if ( i == nmbOfOutputParticles - 1 )
0171             break;
0172
0173         G4double  beta( pd[ i ] / std::sqrt( pd[ i ] * pd[ i ] +
0174                                              invMas[ i ] * invMas[ i ] ) );
0175         for ( int  j( 0 ); j <= i; ++j )
0176             outVec[ j ].lVec->boost( 0, beta, 0 );
0177
0178         ++i;
0179     }
0180
0181     for ( int  j( 0 ); j < nmbOfOutputParticles; ++j )
0182         *outVec[ j ].lVec *= GeV;
0183
0184     //
0185     //---> return the weigth of event
0186     //
0187     return wt;
0188 }
0189
0190
0191 void  CexmcReimplementedGenbod::ParticleChangeHook( void )
0192 {
0193     nmbOfOutputParticles = outVec.size();
0194
0195     if ( nmbOfOutputParticles < 2 || nmbOfOutputParticles > maxParticles )
0196         throw CexmcException( CexmcKinematicsException );
0197
0198     SetMaxWeight();
0199 }
0200
0201
0202 void  CexmcReimplementedGenbod::FermiEnergyDepStatusChangeHook( void )
0203 {
0204     SetMaxWeight();
0205 }
0206
0207
0208 void  CexmcReimplementedGenbod::SetMaxWeight( void )
0209 {
0210     G4double  te_minus_tm( totalEnergy - totalMass );
0211
0212     if ( fermiEnergyDepIsOn )
0213     {
0214         // ffq[] = pi * (2*pi)**(N-2) / (N-2)!
0215         G4double  ffq[] = { 0
0216                      ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
0217                      ,256.3704, 268.4705, 240.9780, 189.2637
0218                      ,132.1308,  83.0202,  47.4210,  24.8295
0219                      ,12.0006,   5.3858,   2.2560,   0.8859 };
0220         maxWeight =
0221                 std::pow( te_minus_tm / GeV, nmbOfOutputParticles - 2 ) *
0222                 ffq[ nmbOfOutputParticles - 1 ] / ( totalEnergy / GeV );
0223     }
0224     else
0225     {
0226         G4double  emmax( ( te_minus_tm + outVec[ 0 ].mass ) / GeV );
0227         G4double  emmin( 0. );
0228         G4double  wtmax( 1. );
0229
0230         for ( G4int  i( 1 ); i < nmbOfOutputParticles; ++i )
0231         {
0232             emmin += outVec[ i - 1 ].mass / GeV;
0233             emmax += outVec[ i ].mass / GeV;
0234             wtmax *= PDK( emmax, emmin, outVec[ i ].mass / GeV );
0235         }
0236         maxWeight = 1 / wtmax;
0237     }
0238 }
0239