include/Geant4/G4DNAOneStepThermalizationModel.hpp

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 // Author: Mathieu Karamitros
0028 //
0029 // WARNING : This class is released as a prototype.
0030 // It might strongly evolve or even disappear in the next releases.
0031 //
0032 // History:
0033 // -----------
0034 // 13 Nov 2016 M.Karamitros created
0035 //
0036 // -------------------------------------------------------------------
0037
0038 #include "G4DNAChemistryManager.hh"
0039 #include "G4DNAMolecularMaterial.hh"
0040 #include "G4DNAWaterExcitationStructure.hh"
0041 #include "G4ITNavigator.hh"
0042 #include "G4Navigator.hh"
0043 #include "G4NistManager.hh"
0044 #include "G4ParticleChangeForGamma.hh"
0045 #include "G4PhysicalConstants.hh"
0046 #include "G4SystemOfUnits.hh"
0047 #include "G4TransportationManager.hh"
0048
0049 #include <memory>
0050
0051 //#define MODEL_VERBOSE
0052
0053 //------------------------------------------------------------------------------
0054
0055 template<typename MODEL>
0056 G4TDNAOneStepThermalizationModel<MODEL>::
0057 G4TDNAOneStepThermalizationModel(const G4ParticleDefinition*,
0058                                  const G4String& nam) :
0059 G4VEmModel(nam)
0060 {
0061   fVerboseLevel = 0;
0062   SetLowEnergyLimit(0.);
0063   G4DNAWaterExcitationStructure exStructure;
0064   SetHighEnergyLimit(exStructure.ExcitationEnergy(0));
0065   fpParticleChangeForGamma = nullptr;
0066   fpWaterDensity = nullptr;
0067 }
0068
0069 //------------------------------------------------------------------------------
0070
0071 template<typename MODEL>
0072 G4TDNAOneStepThermalizationModel<MODEL>::~G4TDNAOneStepThermalizationModel()
0073 = default;
0074
0075 //------------------------------------------------------------------------------
0076 template<typename MODEL>
0077 void G4TDNAOneStepThermalizationModel<MODEL>::
0078 Initialise(const G4ParticleDefinition* particleDefinition,
0079            const G4DataVector&)
0080 {
0081 #ifdef MODEL_VERBOSE
0082   if(fVerboseLevel)
0083     G4cout << "Calling G4DNAOneStepThermalizationModel::Initialise()"
0084            << G4endl;
0085 #endif
0086   if (particleDefinition->GetParticleName() != "e-")
0087   {
0088     G4ExceptionDescription errMsg;
0089     errMsg << "G4DNAOneStepThermalizationModel can only be applied "
0090     "to electrons";
0091     G4Exception("G4DNAOneStepThermalizationModel::CrossSectionPerVolume",
0092                 "G4DNAOneStepThermalizationModel001",
0093                 FatalErrorInArgument,errMsg);
0094     return;
0095   }
0096
0097   if(!fIsInitialised)
0098   {
0099     fIsInitialised = true;
0100     fpParticleChangeForGamma = GetParticleChangeForGamma();
0101   }
0102
0103   G4Navigator* navigator =
0104   G4TransportationManager::GetTransportationManager()->
0105   GetNavigatorForTracking();
0106
0107   fpNavigator = std::make_unique<G4Navigator>();
0108
0109   if(navigator != nullptr){ // add these checks for testing mode
0110     auto world=navigator->GetWorldVolume();
0111     if(world != nullptr){
0112       fpNavigator->SetWorldVolume(world);
0113       //fNavigator->NewNavigatorState();
0114     }
0115   }
0116
0117   fpWaterDensity =
0118   G4DNAMolecularMaterial::Instance()->
0119   GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
0120 }
0121
0122 //------------------------------------------------------------------------------
0123 template<typename MODEL>
0124 G4double G4TDNAOneStepThermalizationModel<MODEL>::
0125 CrossSectionPerVolume(const G4Material* material,
0126                       const G4ParticleDefinition*,
0127                       G4double ekin,
0128                       G4double,
0129                       G4double)
0130 {
0131 #ifdef MODEL_VERBOSE
0132   if(fVerboseLevel > 1)
0133     G4cout << "Calling CrossSectionPerVolume() of G4DNAOneStepThermalizationModel"
0134     << G4endl;
0135 #endif
0136
0137   if(ekin > HighEnergyLimit()){
0138     return 0.0;
0139   }
0140
0141   G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
0142
0143   if(waterDensity!= 0.0){
0144     return DBL_MAX;
0145   }
0146   return 0.;
0147 }
0148
0149 //------------------------------------------------------------------------------
0150 template<typename MODEL>
0151 double G4TDNAOneStepThermalizationModel<MODEL>::GetRmean(double k){
0152   return MODEL::GetRmean(k);
0153 }
0154
0155
0156 //------------------------------------------------------------------------------
0157
0158 template<typename MODEL>
0159 void G4TDNAOneStepThermalizationModel<MODEL>::
0160 GetPenetration(G4double k, G4ThreeVector& displacement)
0161 {
0162   return MODEL::GetPenetration(k, displacement);
0163 }
0164
0165 //------------------------------------------------------------------------------
0166 template<typename MODEL>
0167 void G4TDNAOneStepThermalizationModel<MODEL>::
0168 SampleSecondaries(std::vector<G4DynamicParticle*>*,
0169                   const G4MaterialCutsCouple*,
0170                   const G4DynamicParticle* particle,
0171                   G4double,
0172                   G4double)
0173 {
0174 #ifdef MODEL_VERBOSE
0175   if(fVerboseLevel)
0176     G4cout << "Calling SampleSecondaries() of G4DNAOneStepThermalizationModel"
0177     << G4endl;
0178 #endif
0179
0180   G4double k = particle->GetKineticEnergy();
0181
0182   if (k <= HighEnergyLimit())
0183   {
0184     fpParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
0185     fpParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
0186
0187     if(G4DNAChemistryManager::IsActivated())
0188     {
0189       G4ThreeVector displacement(0,0,0);
0190       GetPenetration(k, displacement);
0191
0192       //______________________________________________________________
0193       const G4Track * theIncomingTrack =
0194       fpParticleChangeForGamma->GetCurrentTrack();
0195       G4ThreeVector finalPosition(theIncomingTrack->GetPosition()+displacement);
0196
0197       fpNavigator->SetWorldVolume(theIncomingTrack->GetTouchable()->
0198                                  GetVolume(theIncomingTrack->GetTouchable()->
0199                                            GetHistoryDepth()));
0200
0201       double displacementMag = displacement.mag();
0202       double safety = DBL_MAX;
0203       G4ThreeVector direction = displacement/displacementMag;
0204
0205       //--
0206       // 6/09/16 - recupere de molecular dissocation
0207       double mag_displacement = displacement.mag();
0208       G4ThreeVector displacement_direction = displacement/mag_displacement;
0209
0210       //     double step = DBL_MAX;
0211       //     step = fNavigator->CheckNextStep(theIncomingTrack->GetPosition(),
0212       //                                     displacement_direction,
0213       //                                     mag_displacement,
0214       //                                     safety);
0215       //
0216       //
0217       //     if(safety < mag_displacement)
0218       //     {
0219       ////       mag_displacement = prNewSafety;
0220       //       finalPosition = theIncomingTrack->GetPosition()
0221       //       + (displacement/displacementMag)*safety*0.80;
0222       //     }
0223       //--
0224
0225       fpNavigator->ResetHierarchyAndLocate(theIncomingTrack->GetPosition(),
0226                                           direction,
0227                                           *((G4TouchableHistory*)
0228                                             theIncomingTrack->GetTouchable()));
0229
0230       fpNavigator->ComputeStep(theIncomingTrack->GetPosition(),
0231                               displacement/displacementMag,
0232                               displacementMag,
0233                               safety);
0234
0235       if(safety <= displacementMag)
0236       {
0237         finalPosition = theIncomingTrack->GetPosition()
0238         + (displacement/displacementMag)*safety*0.80;
0239       }
0240
0241       G4DNAChemistryManager::Instance()->CreateSolvatedElectron(theIncomingTrack,
0242                                                                 &finalPosition);
0243
0244       fpParticleChangeForGamma->SetProposedKineticEnergy(25.e-3*eV);
0245     }
0246   }
0247 }