ChargeExchangeMC/include/CexmcChargeExchangeProductionModel.hh

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:  CexmcChargeExchangeProductionModel.hh
0030  *
0031  *    Description:  charge exchange physics itself
0032  *
0033  *        Version:  1.0
0034  *        Created:  01.11.2009 00:30:46
0035  *       Revision:  none
0036  *       Compiler:  gcc
0037  *
0038  *         Author:  Alexey Radkov (),
0039  *        Company:  PNPI
0040  *
0041  * =============================================================================
0042  */
0043
0044 #ifndef CEXMC_CHARGE_EXCHANGE_PRODUCTION_MODEL_HH
0045 #define CEXMC_CHARGE_EXCHANGE_PRODUCTION_MODEL_HH
0046
0047 #include <G4HadronicInteraction.hh>
0048 #include <G4HadFinalState.hh>
0049 #include <G4HadProjectile.hh>
0050 #include <G4Nucleus.hh>
0051 #include <G4Proton.hh>
0052 #include <G4Neutron.hh>
0053 #include "CexmcProductionModel.hh"
0054 #include "CexmcGenbod.hh"
0055 #include "CexmcReimplementedGenbod.hh"
0056 #include "CexmcException.hh"
0057
0058
0059 template  < typename  OutputParticle >
0060 class  CexmcChargeExchangeProductionModel : public G4HadronicInteraction,
0061                                             public CexmcProductionModel
0062 {
0063     public:
0064         CexmcChargeExchangeProductionModel();
0065
0066         ~CexmcChargeExchangeProductionModel();
0067
0068     public:
0069         G4HadFinalState *  ApplyYourself( const G4HadProjectile &  projectile,
0070                                           G4Nucleus &  targetNucleus );
0071
0072     private:
0073         G4double                    nucleusParticleMass;
0074
0075         CexmcPhaseSpaceGenerator *  phaseSpaceGenerator;
0076 };
0077
0078
0079 template  < typename  OutputParticle >
0080 CexmcChargeExchangeProductionModel< OutputParticle >::
0081                                         CexmcChargeExchangeProductionModel() :
0082     G4HadronicInteraction( CexmcChargeExchangeInteractionName ),
0083     CexmcProductionModel( CexmcChargeExchangeProductionModelName ),
0084     nucleusParticleMass( 0 ), phaseSpaceGenerator( NULL )
0085 {
0086     incidentParticle = G4PionMinus::Definition();
0087     nucleusParticle = G4Proton::Definition();
0088     outputParticle = OutputParticle::Definition();
0089     nucleusOutputParticle = G4Neutron::Definition();
0090
0091     nucleusParticleMass = nucleusParticle->GetPDGMass();
0092
0093     productionModelData.incidentParticle = incidentParticle;
0094     productionModelData.nucleusParticle = nucleusParticle;
0095     productionModelData.outputParticle = outputParticle;
0096     productionModelData.nucleusOutputParticle = nucleusOutputParticle;
0097
0098     CexmcPhaseSpaceInVector   inVec;
0099
0100     inVec.push_back( &productionModelData.incidentParticleSCM );
0101     inVec.push_back( &productionModelData.nucleusParticleSCM );
0102
0103     CexmcPhaseSpaceOutVector  outVec;
0104
0105     outVec.push_back( CexmcPhaseSpaceOutVectorElement(
0106                             &productionModelData.outputParticleSCM,
0107                             outputParticle->GetPDGMass() ) );
0108     outVec.push_back( CexmcPhaseSpaceOutVectorElement(
0109                             &productionModelData.nucleusOutputParticleSCM,
0110                             nucleusOutputParticle->GetPDGMass() ) );
0111
0112 #ifdef CEXMC_USE_GENBOD
0113     phaseSpaceGenerator = new CexmcGenbod;
0114 #else
0115     phaseSpaceGenerator = new CexmcReimplementedGenbod;
0116 #endif
0117
0118     phaseSpaceGenerator->SetParticles( inVec, outVec );
0119 }
0120
0121
0122 template  < typename  OutputParticle >
0123 CexmcChargeExchangeProductionModel< OutputParticle >::
0124                                         ~CexmcChargeExchangeProductionModel()
0125 {
0126     delete phaseSpaceGenerator;
0127 }
0128
0129
0130 template  < typename  OutputParticle >
0131 G4HadFinalState *  CexmcChargeExchangeProductionModel< OutputParticle >::
0132                             ApplyYourself( const G4HadProjectile &  projectile,
0133                                            G4Nucleus &  targetNucleus )
0134 {
0135     theParticleChange.Clear();
0136
0137     G4double           kinEnergy( projectile.GetKineticEnergy() );
0138     G4HadProjectile &  theProjectile( const_cast< G4HadProjectile & >(
0139                                                                 projectile ) );
0140     const G4LorentzRotation &  projToLab(
0141                                     const_cast< const G4LorentzRotation & >(
0142                                             theProjectile.GetTrafoToLab() ) );
0143     productionModelData.incidentParticleLAB = projectile.Get4Momentum();
0144     productionModelData.incidentParticleLAB.transform( projToLab );
0145     productionModelData.nucleusParticleLAB.setPx( 0 );
0146     productionModelData.nucleusParticleLAB.setPy( 0 );
0147     productionModelData.nucleusParticleLAB.setPz( 0 );
0148     productionModelData.nucleusParticleLAB.setE( nucleusParticleMass );
0149
0150     if ( fermiMotionIsOn )
0151     {
0152         G4ThreeVector  targetNucleusMomentum(
0153                                         targetNucleus.GetFermiMomentum() );
0154         G4double       targetNucleusEnergy(
0155                             std::sqrt( targetNucleusMomentum.mag2() +
0156                                 nucleusParticleMass * nucleusParticleMass ) );
0157         productionModelData.nucleusParticleLAB = G4LorentzVector(
0158                                 targetNucleusMomentum, targetNucleusEnergy );
0159     }
0160     productionModelData.nucleusParticleLAB.transform( projToLab );
0161     G4LorentzVector  lVecSum( productionModelData.incidentParticleLAB +
0162                               productionModelData.nucleusParticleLAB );
0163     G4ThreeVector    boostVec( lVecSum.boostVector() );
0164
0165     productionModelData.incidentParticleSCM =
0166             productionModelData.incidentParticleLAB;
0167     productionModelData.nucleusParticleSCM =
0168             productionModelData.nucleusParticleLAB;
0169
0170     productionModelData.incidentParticleSCM.boost( -boostVec );
0171     productionModelData.nucleusParticleSCM.boost( -boostVec );
0172
0173     triggeredAngularRanges.clear();
0174
0175     if ( ! phaseSpaceGenerator->CheckKinematics() )
0176     {
0177         theParticleChange.SetEnergyChange( kinEnergy );
0178         theParticleChange.SetMomentumChange(
0179                                     projectile.Get4Momentum().vect().unit());
0180         return &theParticleChange;
0181     }
0182
0183     do
0184     {
0185         phaseSpaceGenerator->Generate();
0186         G4double  cosTheta( productionModelData.outputParticleSCM.cosTheta() );
0187         for ( CexmcAngularRangeList::iterator  k( angularRanges.begin() );
0188                                                 k != angularRanges.end(); ++k )
0189         {
0190             if ( cosTheta <= k->top && cosTheta > k->bottom )
0191                 triggeredAngularRanges.push_back( CexmcAngularRange(
0192                                                 k->top, k->bottom, k->index ) );
0193         }
0194     } while ( triggeredAngularRanges.empty() );
0195
0196     productionModelData.outputParticleLAB =
0197             productionModelData.outputParticleSCM;
0198     productionModelData.nucleusOutputParticleLAB =
0199             productionModelData.nucleusOutputParticleSCM;
0200
0201     productionModelData.outputParticleLAB.boost( boostVec );
0202     productionModelData.nucleusOutputParticleLAB.boost( boostVec );
0203
0204     theParticleChange.SetStatusChange( stopAndKill );
0205     theParticleChange.SetEnergyChange( 0.0 );
0206
0207     G4DynamicParticle *  secOutParticle( new G4DynamicParticle(
0208                             outputParticle,
0209                             productionModelData.outputParticleLAB ) );
0210     theParticleChange.AddSecondary( secOutParticle );
0211     G4DynamicParticle *  secNeutron( new G4DynamicParticle(
0212                             nucleusOutputParticle,
0213                             productionModelData.nucleusOutputParticleLAB ) );
0214     theParticleChange.AddSecondary( secNeutron );
0215
0216     /* projectile->GetDefinition() shall always be identical to incidentParticle
0217      * as far as CexmcHadronicProcess::IsApplicable() will check that only
0218      * incidentParticle is allowed. Here is mostly unnecessary assignment. */
0219     productionModelData.incidentParticle = projectile.GetDefinition();
0220
0221     return &theParticleChange;
0222 }
0223
0224
0225 #endif
0226