ChargeExchangeMC/src/CexmcHadronicProcess.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:  CexmcHadronicProcess.cc
0030  *
0031  *    Description:  hadronic process with production model
0032  *
0033  *        Version:  1.0
0034  *        Created:  31.10.2009 23:54:38
0035  *       Revision:  none
0036  *       Compiler:  gcc
0037  *
0038  *         Author:  Alexey Radkov (),
0039  *        Company:  PNPI
0040  *
0041  * ============================================================================
0042  */
0043
0044 #include <G4ParticleChange.hh>
0045 #include <G4ParticleDefinition.hh>
0046 #include <G4HadronicInteraction.hh>
0047 #include <G4Track.hh>
0048 #include <G4Step.hh>
0049 #include <G4Element.hh>
0050 #include <G4StableIsotopes.hh>
0051 #include <G4TrackStatus.hh>
0052 #include "CexmcHadronicProcess.hh"
0053 #include "CexmcProductionModel.hh"
0054 #include "CexmcIncidentParticleTrackInfo.hh"
0055 #include "CexmcException.hh"
0056
0057
0058 CexmcHadronicProcess::CexmcHadronicProcess( const G4String &  name ) :
0059     G4HadronicProcess( name ), productionModel( NULL ), interaction( NULL ),
0060     theTotalResult( NULL ), isInitialized( false )
0061 {
0062     theTotalResult = new G4ParticleChange();
0063     SetIntegral(false);
0064 }
0065
0066
0067 CexmcHadronicProcess::~CexmcHadronicProcess()
0068 {
0069     delete theTotalResult;
0070 }
0071
0072
0073 void  CexmcHadronicProcess::RegisterProductionModel(
0074                                                 CexmcProductionModel *  model )
0075 {
0076     productionModel = model;
0077
0078     interaction = dynamic_cast< G4HadronicInteraction * >( productionModel );
0079
0080     if ( ! interaction )
0081         throw CexmcException( CexmcIncompatibleProductionModel );
0082
0083     G4HadronicProcess::RegisterMe( interaction );
0084 }
0085
0086
0087 void  CexmcHadronicProcess::CalculateTargetNucleus(
0088                                                 const G4Material *  material )
0089 {
0090     G4int  numberOfElements( material->GetNumberOfElements() );
0091     if ( numberOfElements > 1 )
0092     {
0093         G4cout << CEXMC_LINE_START "WARNING: Number of elements in target "
0094                   "material is more than 1.\n              Only the first "
0095                   "element will be chosen for target nucleus" << G4endl;
0096     }
0097
0098     const G4Element *  element( material->GetElement( 0 ) );
0099     G4double           ZZ( element->GetZ() );
0100     G4int              Z( G4int( ZZ + 0.5 ) );
0101
0102     G4StableIsotopes  stableIsotopes;
0103     G4int             index( stableIsotopes.GetFirstIsotope( Z ) );
0104     G4double          AA( stableIsotopes.GetIsotopeNucleonCount( index ) );
0105
0106     targetNucleus.SetParameters( AA, ZZ );
0107 }
0108
0109
0110 void  CexmcHadronicProcess::FillTotalResult( G4HadFinalState *  hadFinalState,
0111                                              const G4Track &  track )
0112 {
0113     G4int  numberOfSecondaries( hadFinalState->GetNumberOfSecondaries() );
0114
0115     theTotalResult->Clear();
0116     theTotalResult->Initialize( track );
0117     theTotalResult->SetSecondaryWeightByProcess( true );
0118     theTotalResult->ProposeLocalEnergyDeposit(
0119                                     hadFinalState->GetLocalEnergyDeposit() );
0120     theTotalResult->SetNumberOfSecondaries( numberOfSecondaries );
0121     theTotalResult->ProposeEnergy( hadFinalState->GetEnergyChange() );
0122     theTotalResult->ProposeTrackStatus( fAlive );
0123     if ( hadFinalState->GetStatusChange() == stopAndKill )
0124         theTotalResult->ProposeTrackStatus( fStopAndKill );
0125
0126     for ( G4int  i( 0 ); i < numberOfSecondaries; ++i )
0127     {
0128         G4double   time( hadFinalState->GetSecondary( i )->GetTime() );
0129         if ( time < 0 )
0130             time = track.GetGlobalTime();
0131
0132         G4Track *  newTrack( new G4Track(
0133                              hadFinalState->GetSecondary( i )->GetParticle(),
0134                              time, track.GetPosition() ) );
0135
0136         G4double   newWeight( track.GetWeight() *
0137                               hadFinalState->GetSecondary( i )->GetWeight() );
0138         newTrack->SetWeight( newWeight );
0139         newTrack->SetTouchableHandle( track.GetTouchableHandle() );
0140         theTotalResult->AddSecondary( newTrack );
0141     }
0142
0143     hadFinalState->Clear();
0144 }
0145
0146
0147 G4VParticleChange *  CexmcHadronicProcess::PostStepDoIt( const G4Track &  track,
0148                                                          const G4Step & )
0149 {
0150     G4TrackStatus  trackStatus( track.GetTrackStatus() );
0151
0152     if ( trackStatus != fAlive && trackStatus != fSuspend )
0153     {
0154         theTotalResult->Clear();
0155         theTotalResult->Initialize( track );
0156
0157         return theTotalResult;
0158     }
0159
0160     /* NB: the target nucleus is chosen only once, it means that it will always
0161      * have same Z and A, practically the first stable isotope of the first
0162      * element in elements vector will be chosen. This simplification prompts
0163      * the user to choose simple single-element material for the target, for
0164      * example liquid hydrogen. On the other hand target nucleus is supposedly
0165      * only needed if user decides to turn Fermi motion on, so this
0166      * simplification should not be very harmful */
0167     if ( ! isInitialized )
0168     {
0169         CalculateTargetNucleus( track.GetMaterial() );
0170         isInitialized = true;
0171     }
0172
0173     G4HadProjectile    projectile( track );
0174     G4HadFinalState *  result( interaction->ApplyYourself( projectile,
0175                                                            targetNucleus ) );
0176     FillTotalResult( result, track );
0177
0178     if ( theTotalResult->GetTrackStatus() != fStopAndKill )
0179     {
0180         CexmcTrackInfo *  trackInfo( static_cast< CexmcTrackInfo * >(
0181                                                 track.GetUserInformation() ) );
0182
0183         if ( trackInfo &&
0184              trackInfo->GetTypeInfo() == CexmcIncidentParticleTrackType )
0185         {
0186             CexmcIncidentParticleTrackInfo *  theTrackInfo(
0187                 static_cast< CexmcIncidentParticleTrackInfo * >( trackInfo ) );
0188             theTrackInfo->SetNeedsTrackLengthResampling();
0189         }
0190     }
0191
0192     return theTotalResult;
0193 }
0194
0195
0196 G4bool  CexmcHadronicProcess::IsApplicable(
0197                                         const G4ParticleDefinition &  particle )
0198 {
0199     if ( ! productionModel )
0200         return false;
0201
0202     G4ParticleDefinition *  incidentParticle(
0203                                     productionModel->GetIncidentParticle() );
0204
0205     if ( ! incidentParticle )
0206         return false;
0207
0208     return particle == *incidentParticle;
0209 }
0210