Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file 29 // 29 // 30 // CERN, Geneva, Switzerland 30 // CERN, Geneva, Switzerland 31 // 31 // 32 // File name: G4RKPropagation.cc 32 // File name: G4RKPropagation.cc 33 // 33 // 34 // Author: Alessandro Brunengo (Al 34 // Author: Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it) 35 // 35 // 36 // Creation date: 6 June 2000 36 // Creation date: 6 June 2000 37 // ------------------------------------------- 37 // ------------------------------------------------------------------- 38 #include "G4RKPropagation.hh" 38 #include "G4RKPropagation.hh" 39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh" 40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 41 // nuclear fields 41 // nuclear fields 42 #include "G4VNuclearField.hh" 42 #include "G4VNuclearField.hh" 43 #include "G4ProtonField.hh" 43 #include "G4ProtonField.hh" 44 #include "G4NeutronField.hh" 44 #include "G4NeutronField.hh" 45 #include "G4AntiProtonField.hh" 45 #include "G4AntiProtonField.hh" 46 #include "G4KaonPlusField.hh" 46 #include "G4KaonPlusField.hh" 47 #include "G4KaonMinusField.hh" 47 #include "G4KaonMinusField.hh" 48 #include "G4KaonZeroField.hh" 48 #include "G4KaonZeroField.hh" 49 #include "G4PionPlusField.hh" 49 #include "G4PionPlusField.hh" 50 #include "G4PionMinusField.hh" 50 #include "G4PionMinusField.hh" 51 #include "G4PionZeroField.hh" 51 #include "G4PionZeroField.hh" 52 #include "G4SigmaPlusField.hh" 52 #include "G4SigmaPlusField.hh" 53 #include "G4SigmaMinusField.hh" 53 #include "G4SigmaMinusField.hh" 54 #include "G4SigmaZeroField.hh" 54 #include "G4SigmaZeroField.hh" 55 // particles properties 55 // particles properties 56 #include "G4Proton.hh" 56 #include "G4Proton.hh" 57 #include "G4Neutron.hh" 57 #include "G4Neutron.hh" 58 #include "G4AntiProton.hh" 58 #include "G4AntiProton.hh" 59 #include "G4KaonPlus.hh" 59 #include "G4KaonPlus.hh" 60 #include "G4KaonMinus.hh" 60 #include "G4KaonMinus.hh" 61 #include "G4KaonZero.hh" 61 #include "G4KaonZero.hh" 62 #include "G4PionPlus.hh" 62 #include "G4PionPlus.hh" 63 #include "G4PionMinus.hh" 63 #include "G4PionMinus.hh" 64 #include "G4PionZero.hh" 64 #include "G4PionZero.hh" 65 #include "G4SigmaPlus.hh" 65 #include "G4SigmaPlus.hh" 66 #include "G4SigmaMinus.hh" 66 #include "G4SigmaMinus.hh" 67 #include "G4SigmaZero.hh" 67 #include "G4SigmaZero.hh" 68 68 69 #include "globals.hh" 69 #include "globals.hh" 70 70 71 #include "G4KM_OpticalEqRhs.hh" 71 #include "G4KM_OpticalEqRhs.hh" 72 #include "G4KM_NucleonEqRhs.hh" 72 #include "G4KM_NucleonEqRhs.hh" 73 #include "G4ClassicalRK4.hh" 73 #include "G4ClassicalRK4.hh" 74 #include "G4MagIntegratorDriver.hh" 74 #include "G4MagIntegratorDriver.hh" 75 75 76 #include "G4LorentzRotation.hh" 76 #include "G4LorentzRotation.hh" 77 77 78 // unsigned EncodingHashFun(const G4int& aEnco 78 // unsigned EncodingHashFun(const G4int& aEncoding); 79 79 80 G4RKPropagation::G4RKPropagation() : 80 G4RKPropagation::G4RKPropagation() : 81 theOuterRadius(0), theNucleus(0), 81 theOuterRadius(0), theNucleus(0), 82 theFieldMap(0), theEquationMap(0), 82 theFieldMap(0), theEquationMap(0), 83 theField(0) 83 theField(0) 84 { } 84 { } 85 85 86 86 87 G4RKPropagation::~G4RKPropagation() 87 G4RKPropagation::~G4RKPropagation() 88 { 88 { 89 // free theFieldMap memory 89 // free theFieldMap memory 90 if(theFieldMap) delete_FieldsAndMap(theFiel 90 if(theFieldMap) delete_FieldsAndMap(theFieldMap); 91 91 92 // free theEquationMap memory 92 // free theEquationMap memory 93 if(theEquationMap) delete_EquationsAndMap(t 93 if(theEquationMap) delete_EquationsAndMap(theEquationMap); 94 94 95 if (theField) delete theField; 95 if (theField) delete theField; 96 } 96 } 97 97 98 //-------------------------------------------- 98 //---------------------------------------------------------------------------- 99 void G4RKPropagation::Init(G4V3DNucleus * nucl 99 void G4RKPropagation::Init(G4V3DNucleus * nucleus) 100 //-------------------------------------------- 100 //---------------------------------------------------------------------------- 101 { 101 { 102 102 103 // free theFieldMap memory 103 // free theFieldMap memory 104 if(theFieldMap) delete_FieldsAndMap(theFiel 104 if(theFieldMap) delete_FieldsAndMap(theFieldMap); 105 105 106 // free theEquationMap memory 106 // free theEquationMap memory 107 if(theEquationMap) delete_EquationsAndMap(t 107 if(theEquationMap) delete_EquationsAndMap(theEquationMap); 108 108 109 if (theField) delete theField; 109 if (theField) delete theField; 110 110 111 // Initialize the nuclear field map. 111 // Initialize the nuclear field map. 112 theNucleus = nucleus; 112 theNucleus = nucleus; 113 theOuterRadius = theNucleus->GetOuterRadius 113 theOuterRadius = theNucleus->GetOuterRadius(); 114 114 115 theFieldMap = new std::map <G4int, G4VNucle 115 theFieldMap = new std::map <G4int, G4VNuclearField*, std::less<G4int> >; 116 116 117 (*theFieldMap)[G4Proton::Proton()->GetPDGEn 117 (*theFieldMap)[G4Proton::Proton()->GetPDGEncoding()] = new G4ProtonField(theNucleus); 118 (*theFieldMap)[G4Neutron::Neutron()->GetPDG 118 (*theFieldMap)[G4Neutron::Neutron()->GetPDGEncoding()] = new G4NeutronField(theNucleus); 119 (*theFieldMap)[G4AntiProton::AntiProton()-> 119 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = new G4AntiProtonField(theNucleus); 120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetP 120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = new G4KaonPlusField(theNucleus); 121 (*theFieldMap)[G4KaonMinus::KaonMinus()->Ge 121 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = new G4KaonMinusField(theNucleus); 122 (*theFieldMap)[G4KaonZero::KaonZero()->GetP 122 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = new G4KaonZeroField(theNucleus); 123 (*theFieldMap)[G4PionPlus::PionPlus()->GetP 123 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = new G4PionPlusField(theNucleus); 124 (*theFieldMap)[G4PionMinus::PionMinus()->Ge 124 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = new G4PionMinusField(theNucleus); 125 (*theFieldMap)[G4PionZero::PionZero()->GetP 125 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()] = new G4PionZeroField(theNucleus); 126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->Ge 126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = new G4SigmaPlusField(theNucleus); 127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()-> 127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = new G4SigmaMinusField(theNucleus); 128 (*theFieldMap)[G4SigmaZero::SigmaZero()->Ge 128 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = new G4SigmaZeroField(theNucleus); 129 129 130 theEquationMap = new std::map <G4int, G4Mag 130 theEquationMap = new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >; 131 131 132 // theField needed by the design of G4Mag_e 132 // theField needed by the design of G4Mag_eqRhs 133 theField = new G4KM_DummyField; //Field 133 theField = new G4KM_DummyField; //Field not needed for integration 134 G4KM_OpticalEqRhs * opticalEq; 134 G4KM_OpticalEqRhs * opticalEq; 135 G4KM_NucleonEqRhs * nucleonEq; 135 G4KM_NucleonEqRhs * nucleonEq; 136 G4double mass; 136 G4double mass; 137 G4double opticalCoeff; 137 G4double opticalCoeff; 138 138 139 nucleonEq = new G4KM_NucleonEqRhs(theField, 139 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus); 140 mass = G4Proton::Proton()->GetPDGMass(); 140 mass = G4Proton::Proton()->GetPDGMass(); 141 nucleonEq->SetMass(mass); 141 nucleonEq->SetMass(mass); 142 (*theEquationMap)[G4Proton::Proton()->GetPD 142 (*theEquationMap)[G4Proton::Proton()->GetPDGEncoding()] = nucleonEq; 143 143 144 nucleonEq = new G4KM_NucleonEqRhs(theField, 144 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus); 145 mass = G4Neutron::Neutron()->GetPDGMass(); 145 mass = G4Neutron::Neutron()->GetPDGMass(); 146 nucleonEq->SetMass(mass); 146 nucleonEq->SetMass(mass); 147 (*theEquationMap)[G4Neutron::Neutron()->Get 147 (*theEquationMap)[G4Neutron::Neutron()->GetPDGEncoding()] = nucleonEq; 148 148 149 opticalEq = new G4KM_OpticalEqRhs(theField, 149 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 150 mass = G4AntiProton::AntiProton()->GetPDGMa 150 mass = G4AntiProton::AntiProton()->GetPDGMass(); 151 opticalCoeff = 151 opticalCoeff = 152 (*theFieldMap)[G4AntiProton::AntiProt 152 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()]->GetCoeff(); 153 opticalEq->SetFactor(mass,opticalCoeff); 153 opticalEq->SetFactor(mass,opticalCoeff); 154 (*theEquationMap)[G4AntiProton::AntiProton( 154 (*theEquationMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = opticalEq; 155 155 156 opticalEq = new G4KM_OpticalEqRhs(theField, 156 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 157 mass = G4KaonPlus::KaonPlus()->GetPDGMass() 157 mass = G4KaonPlus::KaonPlus()->GetPDGMass(); 158 opticalCoeff = 158 opticalCoeff = 159 (*theFieldMap)[G4KaonPlus::KaonPlus() 159 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()]->GetCoeff(); 160 opticalEq->SetFactor(mass,opticalCoeff); 160 opticalEq->SetFactor(mass,opticalCoeff); 161 (*theEquationMap)[G4KaonPlus::KaonPlus()->G 161 (*theEquationMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = opticalEq; 162 162 163 opticalEq = new G4KM_OpticalEqRhs(theField, 163 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 164 mass = G4KaonMinus::KaonMinus()->GetPDGMass 164 mass = G4KaonMinus::KaonMinus()->GetPDGMass(); 165 opticalCoeff = 165 opticalCoeff = 166 (*theFieldMap)[G4KaonMinus::KaonMinus 166 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()]->GetCoeff(); 167 opticalEq->SetFactor(mass,opticalCoeff); 167 opticalEq->SetFactor(mass,opticalCoeff); 168 (*theEquationMap)[G4KaonMinus::KaonMinus()- 168 (*theEquationMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = opticalEq; 169 169 170 opticalEq = new G4KM_OpticalEqRhs(theField, 170 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 171 mass = G4KaonZero::KaonZero()->GetPDGMass() 171 mass = G4KaonZero::KaonZero()->GetPDGMass(); 172 opticalCoeff = 172 opticalCoeff = 173 (*theFieldMap)[G4KaonZero::KaonZero() 173 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()]->GetCoeff(); 174 opticalEq->SetFactor(mass,opticalCoeff); 174 opticalEq->SetFactor(mass,opticalCoeff); 175 (*theEquationMap)[G4KaonZero::KaonZero()->G 175 (*theEquationMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = opticalEq; 176 176 177 opticalEq = new G4KM_OpticalEqRhs(theField, 177 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 178 mass = G4PionPlus::PionPlus()->GetPDGMass() 178 mass = G4PionPlus::PionPlus()->GetPDGMass(); 179 opticalCoeff = 179 opticalCoeff = 180 (*theFieldMap)[G4PionPlus::PionPlus() 180 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()]->GetCoeff(); 181 opticalEq->SetFactor(mass,opticalCoeff); 181 opticalEq->SetFactor(mass,opticalCoeff); 182 (*theEquationMap)[G4PionPlus::PionPlus()->G 182 (*theEquationMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = opticalEq; 183 183 184 opticalEq = new G4KM_OpticalEqRhs(theField, 184 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 185 mass = G4PionMinus::PionMinus()->GetPDGMass 185 mass = G4PionMinus::PionMinus()->GetPDGMass(); 186 opticalCoeff = 186 opticalCoeff = 187 (*theFieldMap)[G4PionMinus::PionMinus 187 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()]->GetCoeff(); 188 opticalEq->SetFactor(mass,opticalCoeff); 188 opticalEq->SetFactor(mass,opticalCoeff); 189 (*theEquationMap)[G4PionMinus::PionMinus()- 189 (*theEquationMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = opticalEq; 190 190 191 opticalEq = new G4KM_OpticalEqRhs(theField, 191 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 192 mass = G4PionZero::PionZero()->GetPDGMass() 192 mass = G4PionZero::PionZero()->GetPDGMass(); 193 opticalCoeff = 193 opticalCoeff = 194 (*theFieldMap)[G4PionZero::PionZero() 194 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()]->GetCoeff(); 195 opticalEq->SetFactor(mass,opticalCoeff); 195 opticalEq->SetFactor(mass,opticalCoeff); 196 (*theEquationMap)[G4PionZero::PionZero()->G 196 (*theEquationMap)[G4PionZero::PionZero()->GetPDGEncoding()] = opticalEq; 197 197 198 opticalEq = new G4KM_OpticalEqRhs(theField, 198 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass 199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass(); 200 opticalCoeff = 200 opticalCoeff = 201 (*theFieldMap)[G4SigmaPlus::SigmaPlus 201 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()]->GetCoeff(); 202 opticalEq->SetFactor(mass,opticalCoeff); 202 opticalEq->SetFactor(mass,opticalCoeff); 203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()- 203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = opticalEq; 204 204 205 opticalEq = new G4KM_OpticalEqRhs(theField, 205 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMa 206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMass(); 207 opticalCoeff = 207 opticalCoeff = 208 (*theFieldMap)[G4SigmaMinus::SigmaMin 208 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()]->GetCoeff(); 209 opticalEq->SetFactor(mass,opticalCoeff); 209 opticalEq->SetFactor(mass,opticalCoeff); 210 (*theEquationMap)[G4SigmaMinus::SigmaMinus( 210 (*theEquationMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = opticalEq; 211 211 212 opticalEq = new G4KM_OpticalEqRhs(theField, 212 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus); 213 mass = G4SigmaZero::SigmaZero()->GetPDGMass 213 mass = G4SigmaZero::SigmaZero()->GetPDGMass(); 214 opticalCoeff = 214 opticalCoeff = 215 (*theFieldMap)[G4SigmaZero::SigmaZero 215 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()]->GetCoeff(); 216 opticalEq->SetFactor(mass,opticalCoeff); 216 opticalEq->SetFactor(mass,opticalCoeff); 217 (*theEquationMap)[G4SigmaZero::SigmaZero()- 217 (*theEquationMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = opticalEq; 218 } 218 } 219 219 220 220 221 //#define debug_1_RKPropagation 1 221 //#define debug_1_RKPropagation 1 222 //-------------------------------------------- 222 //---------------------------------------------------------------------------- 223 void G4RKPropagation::Transport(G4KineticTrack 223 void G4RKPropagation::Transport(G4KineticTrackVector & active, 224 //-------------------------------------- 224 //---------------------------------------------------------------------------- 225 const G4KineticTrackVector &, 225 const G4KineticTrackVector &, 226 G4double timeStep) 226 G4double timeStep) 227 { 227 { 228 // reset momentum transfer to field 228 // reset momentum transfer to field 229 theMomentumTranfer=G4ThreeVector(0,0,0); 229 theMomentumTranfer=G4ThreeVector(0,0,0); 230 230 231 // Loop over tracks 231 // Loop over tracks 232 232 233 std::vector<G4KineticTrack *>::iterator i; 233 std::vector<G4KineticTrack *>::iterator i; 234 for(i = active.begin(); i != active.end(); 234 for(i = active.begin(); i != active.end(); ++i) 235 { 235 { 236 G4double currTimeStep = timeStep; 236 G4double currTimeStep = timeStep; 237 G4KineticTrack * kt = *i; 237 G4KineticTrack * kt = *i; 238 G4int encoding = kt->GetDefinition()->Ge 238 G4int encoding = kt->GetDefinition()->GetPDGEncoding(); 239 239 240 std::map <G4int, G4VNuclearField*, std:: 240 std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding); 241 241 242 G4VNuclearField* currentField=0; 242 G4VNuclearField* currentField=0; 243 if ( fieldIter != theFieldMap->end() ) c 243 if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second; 244 244 245 // debug 245 // debug 246 // if ( timeStep > 1e30 ) { 246 // if ( timeStep > 1e30 ) { 247 // G4cout << " Name :" << kt->GetDefini 247 // G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl; 248 // } 248 // } 249 249 250 // Get the time of intersections with th 250 // Get the time of intersections with the nucleus surface. 251 G4double t_enter, t_leave; 251 G4double t_enter, t_leave; 252 // if the particle does not intersecate 252 // if the particle does not intersecate with the nucleus go to next particle 253 if(!GetSphereIntersectionTimes(kt, t_ent 253 if(!GetSphereIntersectionTimes(kt, t_enter, t_leave)) 254 { 254 { 255 kt->SetState(G4KineticTrack::miss_nuc 255 kt->SetState(G4KineticTrack::miss_nucleus); 256 continue; 256 continue; 257 } 257 } 258 258 259 259 260 #ifdef debug_1_RKPropagation 260 #ifdef debug_1_RKPropagation 261 G4cout <<" kt,timeStep, Intersection tim 261 G4cout <<" kt,timeStep, Intersection times tenter, tleave " 262 <<kt<< " / state= " <<kt->GetState 262 <<kt<< " / state= " <<kt->GetState() <<" / " <<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl; 263 #endif 263 #endif 264 264 265 // if the particle is already outside nu 265 // if the particle is already outside nucleus go to next @@GF should never happen? check! 266 // does happen for particles added as l 266 // does happen for particles added as late.... 267 // if(t_leave < 0 ) 267 // if(t_leave < 0 ) 268 // { 268 // { 269 // throw G4HadronicException(__FI 269 // throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a nucleus"); 270 // continue; 270 // continue; 271 // } 271 // } 272 272 273 // Apply a straight line propagation for 273 // Apply a straight line propagation for particle types 274 // not included in the model 274 // not included in the model 275 if( ! currentField ) 275 if( ! currentField ) 276 { 276 { 277 if(currTimeStep == DBL_MAX)currTimeSt 277 if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05; 278 FreeTransport(kt, currTimeStep); 278 FreeTransport(kt, currTimeStep); 279 if ( currTimeStep >= t_leave ) 279 if ( currTimeStep >= t_leave ) 280 { 280 { 281 if ( kt->GetState() == G4KineticTr 281 if ( kt->GetState() == G4KineticTrack::inside ) 282 { kt->SetState(G4KineticTrack::gon 282 { kt->SetState(G4KineticTrack::gone_out); } 283 else 283 else 284 { kt->SetState(G4KineticTrack::mis 284 { kt->SetState(G4KineticTrack::miss_nucleus);} 285 } else if (kt->GetState() == G4Kineti 285 } else if (kt->GetState() == G4KineticTrack::outside && currTimeStep >= t_enter ){ 286 kt->SetState(G4KineticTrack::insid 286 kt->SetState(G4KineticTrack::inside); 287 } 287 } 288 288 289 continue; 289 continue; 290 } 290 } 291 291 292 if(t_enter > 0) // the particle is out. 292 if(t_enter > 0) // the particle is out. Transport free to the surface 293 { 293 { 294 if(t_enter > currTimeStep) // the pa 294 if(t_enter > currTimeStep) // the particle won't enter the nucleus 295 { 295 { 296 FreeTransport(kt, currTimeStep); 296 FreeTransport(kt, currTimeStep); 297 continue; 297 continue; 298 } 298 } 299 else 299 else 300 { 300 { 301 FreeTransport(kt, t_enter); // g 301 FreeTransport(kt, t_enter); // go to surface 302 currTimeStep -= t_enter; 302 currTimeStep -= t_enter; 303 t_leave -= t_enter; // time 303 t_leave -= t_enter; // time left to leave nucleus 304 // on the surface the particle loo 304 // on the surface the particle loose the barrier energy 305 // G4double newE = mom.e()-(*theF 305 // G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier(); 306 // GetField = Barrier + FermiP 306 // GetField = Barrier + FermiPotential 307 G4double newE = kt->GetTrackingMom 307 G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition()); 308 308 309 if(newE <= kt->GetActualMass()) / 309 if(newE <= kt->GetActualMass()) // the particle cannot enter the nucleus 310 { 310 { 311 // FixMe: should be "pushed bac 311 // FixMe: should be "pushed back?" 312 // for the moment take it 312 // for the moment take it past the nucleus, so we'll not worry next time.. 313 FreeTransport(kt, 1.1*t_leave); 313 FreeTransport(kt, 1.1*t_leave); // take past nucleus 314 kt->SetState(G4KineticTrack::mi 314 kt->SetState(G4KineticTrack::miss_nucleus); 315 // G4cout << "G4RKPropagatio 315 // G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl; 316 // G4cout << " enter nucleus 316 // G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl; 317 // G4cout << " the Field "<< 317 // G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl; 318 // G4cout << " the particl 318 // G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl; 319 continue; 319 continue; 320 } 320 } 321 // 321 // 322 G4double newP = std::sqrt(newE*new 322 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 323 G4LorentzVector new4Mom(newP*kt->G 323 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 324 G4ThreeVector transfer(kt->GetTrac 324 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 325 G4ThreeVector boost= transfer / st 325 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 326 new4Mom*=G4LorentzRotation(boost); 326 new4Mom*=G4LorentzRotation(boost); 327 kt->SetTrackingMomentum(new4Mom); 327 kt->SetTrackingMomentum(new4Mom); 328 kt->SetState(G4KineticTrack::insid 328 kt->SetState(G4KineticTrack::inside); 329 329 330 /* 330 /* 331 G4cout <<" Enter Nucleus - E/Field/Sum: " 331 G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / " 332 << (*theFieldMap)[encoding]->GetField 332 << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / " 333 << kt->GetTrackingMomentum().e()-currentF 333 << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition()) 334 << G4endl 334 << G4endl 335 << " Barrier / field just inside nucleus 335 << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())" 336 << (*theFieldMap)[encoding]->GetBarrier() 336 << (*theFieldMap)[encoding]->GetBarrier() << " / " 337 << (*theFieldMap)[encoding]->GetField(0.9 337 << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition()) 338 << G4endl; 338 << G4endl; 339 */ 339 */ 340 } 340 } 341 } 341 } 342 342 343 // FixMe: should I add a control on theC 343 // FixMe: should I add a control on theCutOnP here? 344 // Transport the particle into the nucle 344 // Transport the particle into the nucleus 345 // G4cerr << "RKPropagation t_leav 345 // G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl; 346 G4bool is_exiting=false; 346 G4bool is_exiting=false; 347 if(currTimeStep > t_leave) // particle 347 if(currTimeStep > t_leave) // particle will exit from the nucleus 348 { 348 { 349 currTimeStep = t_leave; 349 currTimeStep = t_leave; 350 is_exiting=true; 350 is_exiting=true; 351 } 351 } 352 352 353 #ifdef debug_1_RKPropagation 353 #ifdef debug_1_RKPropagation 354 G4cerr << "RKPropagation is_exiting?, t_ 354 G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl; 355 G4cout << "RKPropagation Ekin, field, pr 355 G4cout << "RKPropagation Ekin, field, projectile potential, p " 356 << kt->GetTrackingMomentum().e() - 356 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " " 357 << kt->GetPosition()<<" " 357 << kt->GetPosition()<<" " 358 << G4endl << currentField->GetFiel 358 << G4endl << currentField->GetField(kt->GetPosition()) << " " 359 << kt->GetProjectilePotential()<< 359 << kt->GetProjectilePotential()<< G4endl 360 << kt->GetTrackingMomentum() 360 << kt->GetTrackingMomentum() 361 << G4endl; 361 << G4endl; 362 #endif 362 #endif 363 363 364 G4LorentzVector momold=kt->GetTrackingMo 364 G4LorentzVector momold=kt->GetTrackingMomentum(); 365 G4ThreeVector posold=kt->GetPosition(); 365 G4ThreeVector posold=kt->GetPosition(); 366 366 367 // if (currentField->GetField(kt->Get 367 // if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() || 368 if (currTimeStep > 0 && 368 if (currTimeStep > 0 && 369 ! FieldTransport(kt, currTimeStep) 369 ! FieldTransport(kt, currTimeStep)) { 370 FreeTransport(kt,currTimeStep); 370 FreeTransport(kt,currTimeStep); 371 } 371 } 372 372 373 #ifdef debug_1_RKPropagation 373 #ifdef debug_1_RKPropagation 374 G4cout << "RKPropagation Ekin, field, p 374 G4cout << "RKPropagation Ekin, field, p " 375 << kt->GetTrackingMomentum().e() - 375 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " " 376 << G4endl << currentField->GetFiel 376 << G4endl << currentField->GetField(kt->GetPosition())<< G4endl 377 << kt->GetTrackingMomentum() 377 << kt->GetTrackingMomentum() 378 << G4endl 378 << G4endl 379 << "delta p " << momold-kt->GetTra 379 << "delta p " << momold-kt->GetTrackingMomentum() << G4endl 380 << "del pos " << posold-kt->GetPos 380 << "del pos " << posold-kt->GetPosition() 381 << G4endl; 381 << G4endl; 382 #endif 382 #endif 383 383 384 // complete the transport 384 // complete the transport 385 // FixMe: in some cases there could be a 385 // FixMe: in some cases there could be a significant 386 // part to do still in the nucleu 386 // part to do still in the nucleus, or we stepped to far... depending on 387 // slope of potential 387 // slope of potential 388 G4double t_in=-1, t_out=0; // set onto 388 G4double t_in=-1, t_out=0; // set onto boundary. 389 389 390 // should go out, or are already out by 390 // should go out, or are already out by a too long step.. 391 if(is_exiting || 391 if(is_exiting || 392 (GetSphereIntersectionTimes(kt, t_ 392 (GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 )) // particle is exiting 393 { 393 { 394 if(t_in < 0 && t_out >= 0) //still 394 if(t_in < 0 && t_out >= 0) //still inside, transport safely out. 395 { 395 { 396 // transport free to a position th 396 // transport free to a position that is surely out of the nucleus, to avoid 397 // a new transportation and a new 397 // a new transportation and a new adding the barrier next loop. 398 G4ThreeVector savePos = kt->GetPos 398 G4ThreeVector savePos = kt->GetPosition(); 399 FreeTransport(kt, t_out); 399 FreeTransport(kt, t_out); 400 // and evaluate the right the ener 400 // and evaluate the right the energy 401 G4double newE=kt->GetTrackingMomen 401 G4double newE=kt->GetTrackingMomentum().e(); 402 402 403 // G4cout << " V pos/savePos << " 403 // G4cout << " V pos/savePos << " 404 // << (*theFieldMap)[encoding]- 404 // << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / " 405 // << (*theFieldMap)[encoding]- 405 // << (*theFieldMap)[encoding]->GetField(savePos) 406 // << G4endl; 406 // << G4endl; 407 407 408 if ( std::abs(currentField->GetFie 408 if ( std::abs(currentField->GetField(savePos)) > 0. && 409 std::abs(currentField->GetFi 409 std::abs(currentField->GetField(kt->GetPosition())) > 0.) 410 { // FixMe GF: savePos/pos may be 410 { // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0 411 // This wrongly adds 411 // This wrongly adds or subtracts the Barrier here while 412 // this is done later 412 // this is done later. 413 newE += currentField->GetField( 413 newE += currentField->GetField(savePos) 414 - currentField 414 - currentField->GetField(kt->GetPosition()); 415 } 415 } 416 416 417 // G4cout << " go border nuc 417 // G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 418 418 419 if(newE < kt->GetActualMass()) 419 if(newE < kt->GetActualMass()) 420 { 420 { 421 #ifdef debug_1_RKPropagation 421 #ifdef debug_1_RKPropagation 422 G4cout << "RKPropagation-Transp 422 G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl; 423 G4cout << " cannot leave nucleu 423 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 424 #endif 424 #endif 425 if (kt->GetDefinition() == G4Pr 425 if (kt->GetDefinition() == G4Proton::Proton() || 426 kt->GetDefinition() == G4 426 kt->GetDefinition() == G4Neutron::Neutron() ) { 427 kt->SetState(G4KineticTrack: 427 kt->SetState(G4KineticTrack::captured); 428 } else { 428 } else { 429 kt->SetState(G4KineticTrack: 429 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix 430 } 430 } 431 continue; // the particle canno 431 continue; // the particle cannot exit the nucleus 432 } 432 } 433 G4double newP = std::sqrt(newE*new 433 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 434 G4LorentzVector new4Mom(newP*kt->G 434 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 435 G4ThreeVector transfer(kt->GetTrac 435 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 436 G4ThreeVector boost= transfer / st 436 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 437 new4Mom*=G4LorentzRotation(boost); 437 new4Mom*=G4LorentzRotation(boost); 438 kt->SetTrackingMomentum(new4Mom); 438 kt->SetTrackingMomentum(new4Mom); 439 } 439 } 440 // add the potential barrier 440 // add the potential barrier 441 // FixMe the Coulomb field is not par 441 // FixMe the Coulomb field is not parallel to mom, this is simple approximation 442 G4double newE = kt->GetTrackingMoment 442 G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition()); 443 if(newE < kt->GetActualMass()) 443 if(newE < kt->GetActualMass()) 444 { // the particle cannot exit the nu 444 { // the particle cannot exit the nucleus @@@ GF check. 445 #ifdef debug_1_RKPropagation 445 #ifdef debug_1_RKPropagation 446 G4cout << " cannot leave nucleus, 446 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl; 447 #endif 447 #endif 448 if (kt->GetDefinition() == G4Proto 448 if (kt->GetDefinition() == G4Proton::Proton() || 449 kt->GetDefinition() == G4Neu 449 kt->GetDefinition() == G4Neutron::Neutron() ) { 450 kt->SetState(G4KineticTrack::ca 450 kt->SetState(G4KineticTrack::captured); 451 } else { 451 } else { 452 kt->SetState(G4KineticTrack::go 452 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix 453 } 453 } 454 continue; 454 continue; 455 } 455 } 456 G4double newP = std::sqrt(newE*newE- 456 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass())); 457 G4LorentzVector new4Mom(newP*kt->GetT 457 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE); 458 G4ThreeVector transfer(kt->GetTrackin 458 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect()); 459 G4ThreeVector boost= transfer / std:: 459 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass())); 460 new4Mom*=G4LorentzRotation(boost); 460 new4Mom*=G4LorentzRotation(boost); 461 kt->SetTrackingMomentum(new4Mom); 461 kt->SetTrackingMomentum(new4Mom); 462 kt->SetState(G4KineticTrack::gone_out 462 kt->SetState(G4KineticTrack::gone_out); 463 } 463 } 464 464 465 } 465 } 466 466 467 } 467 } 468 468 469 469 470 //-------------------------------------------- 470 //---------------------------------------------------------------------------- 471 G4ThreeVector G4RKPropagation::GetMomentumTran 471 G4ThreeVector G4RKPropagation::GetMomentumTransfer() const 472 //-------------------------------------------- 472 //---------------------------------------------------------------------------- 473 { 473 { 474 return theMomentumTranfer; 474 return theMomentumTranfer; 475 } 475 } 476 476 477 477 478 //-------------------------------------------- 478 //---------------------------------------------------------------------------- 479 G4bool G4RKPropagation::FieldTransport(G4Kinet 479 G4bool G4RKPropagation::FieldTransport(G4KineticTrack * kt, const G4double timeStep) 480 //-------------------------------------------- 480 //---------------------------------------------------------------------------- 481 { 481 { >> 482 theMomentumTranfer=G4ThreeVector(0,0,0); 482 // G4cout <<"Stepper input"<<kt->GetTrac 483 // G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl; 483 // create the integrator stepper 484 // create the integrator stepper 484 // G4Mag_EqRhs * equation = mapIter->sec 485 // G4Mag_EqRhs * equation = mapIter->second; 485 G4Mag_EqRhs * equation = (*theEquationMap)[ 486 G4Mag_EqRhs * equation = (*theEquationMap)[kt->GetDefinition()->GetPDGEncoding()]; 486 G4MagIntegratorStepper * stepper = new G4Cl 487 G4MagIntegratorStepper * stepper = new G4ClassicalRK4(equation); 487 488 488 // create the integrator driver 489 // create the integrator driver 489 G4double hMin = 1.0e-25*second; // arbitr 490 G4double hMin = 1.0e-25*second; // arbitrary choice. Means 0.03 fm at c 490 G4MagInt_Driver * driver = new G4MagInt_Dri 491 G4MagInt_Driver * driver = new G4MagInt_Driver(hMin, stepper); 491 492 492 // Temporary: use driver->AccurateAdvance() 493 // Temporary: use driver->AccurateAdvance() 493 // create the G4FieldTrack needed by Accura 494 // create the G4FieldTrack needed by AccurateAdvance 494 G4double curveLength = 0; 495 G4double curveLength = 0; 495 G4FieldTrack track(kt->GetPosition(), 496 G4FieldTrack track(kt->GetPosition(), 496 kt->GetTrackingMomentum().vect().unit 497 kt->GetTrackingMomentum().vect().unit(), // momentum direction 497 curveLength, // curvelength 498 curveLength, // curvelength 498 kt->GetTrackingMomentum().e()-kt->Get 499 kt->GetTrackingMomentum().e()-kt->GetActualMass(), // kinetic energy 499 kt->GetActualMass(), // restmass 500 kt->GetActualMass(), // restmass 500 kt->GetTrackingMomentum().beta()*c_li 501 kt->GetTrackingMomentum().beta()*c_light); // velocity 501 // integrate 502 // integrate 502 G4double eps = 0.01; 503 G4double eps = 0.01; 503 // G4cout << "currTimeStep = " << currTi 504 // G4cout << "currTimeStep = " << currTimeStep << G4endl; 504 if(!driver->AccurateAdvance(track, timeStep 505 if(!driver->AccurateAdvance(track, timeStep, eps)) 505 { // cannot track this particle 506 { // cannot track this particle 506 #ifdef debug_1_RKPropagation 507 #ifdef debug_1_RKPropagation 507 std::cerr << "G4RKPropagation::FieldTran 508 std::cerr << "G4RKPropagation::FieldTransport() warning: integration error." 508 << G4endl << "position " << kt->Ge 509 << G4endl << "position " << kt->GetPosition() << " 4mom " <<kt->GetTrackingMomentum() 509 <<G4endl << " timestep " <<timeSte 510 <<G4endl << " timestep " <<timeStep 510 << G4endl; 511 << G4endl; 511 #endif 512 #endif 512 delete driver; 513 delete driver; 513 delete stepper; 514 delete stepper; 514 return false; 515 return false; 515 } 516 } 516 /* 517 /* 517 G4cout <<" E/Field/Sum be4 : " <<mom.e( 518 G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / " 518 << (*theFieldMap)[encoding]->GetFie 519 << (*theFieldMap)[encoding]->GetField(pos) << " / " 519 << mom.e()+(*theFieldMap)[encoding]->Ge 520 << mom.e()+(*theFieldMap)[encoding]->GetField(pos) 520 << G4endl; 521 << G4endl; 521 */ 522 */ 522 523 523 // Correct for momentum ( thus energy) tran << 524 // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nuclues frame. 524 G4ThreeVector MomentumTranfer = kt->GetTrac 525 G4ThreeVector MomentumTranfer = kt->GetTrackingMomentum().vect() - track.GetMomentum(); 525 G4ThreeVector boost= MomentumTranfer / std: 526 G4ThreeVector boost= MomentumTranfer / std::sqrt (MomentumTranfer.mag2() +sqr(theNucleus->GetMass())); 526 527 527 // update the kt 528 // update the kt 528 kt->SetPosition(track.GetPosition()); 529 kt->SetPosition(track.GetPosition()); 529 G4LorentzVector mom(track.GetMomentum(),std 530 G4LorentzVector mom(track.GetMomentum(),std::sqrt(track.GetMomentum().mag2() + sqr(kt->GetActualMass()))); 530 mom *= G4LorentzRotation( boost ); 531 mom *= G4LorentzRotation( boost ); 531 theMomentumTranfer += ( kt->GetTrackingMome 532 theMomentumTranfer += ( kt->GetTrackingMomentum() - mom ).vect(); 532 kt->SetTrackingMomentum(mom); 533 kt->SetTrackingMomentum(mom); 533 534 534 // G4cout <<"Stepper output"<<kt<<" "<<k 535 // G4cout <<"Stepper output"<<kt<<" "<<kt->GetTrackingMomentum()<<" "<<kt->GetPosition()<<G4endl; 535 /* 536 /* 536 * G4ThreeVector MomentumTranfer2=kt->Get 537 * G4ThreeVector MomentumTranfer2=kt->GetTrackingMomentum().vect() - mom.vect(); 537 * G4cout << " MomentumTransfer/corrected" 538 * G4cout << " MomentumTransfer/corrected" << MomentumTranfer << " " << MomentumTranfer.mag() 538 * << " " << MomentumTranfer2 << 539 * << " " << MomentumTranfer2 << " " << MomentumTranfer2.mag() << " " 539 * << MomentumTranfer-MomentumTranfer2 540 * << MomentumTranfer-MomentumTranfer2 << " "<< 540 * MomentumTranfer-MomentumTranfer2.mag 541 * MomentumTranfer-MomentumTranfer2.mag() << " " << G4endl; 541 * G4cout <<" E/Field/Sum aft : " <<mo 542 * G4cout <<" E/Field/Sum aft : " <<mom.e() << " / " 542 * << " / " << (*theFieldMap)[en 543 * << " / " << (*theFieldMap)[encoding]->GetField(pos)<< " / " 543 * << mom.e()+(*theFieldMap)[encoding] 544 * << mom.e()+(*theFieldMap)[encoding]->GetField(pos) 544 * << G4endl; 545 * << G4endl; 545 */ 546 */ 546 547 547 delete driver; 548 delete driver; 548 delete stepper; 549 delete stepper; 549 return true; 550 return true; 550 } 551 } 551 552 552 //-------------------------------------------- 553 //---------------------------------------------------------------------------- 553 G4bool G4RKPropagation::FreeTransport(G4Kineti 554 G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep) 554 //-------------------------------------------- 555 //---------------------------------------------------------------------------- 555 { 556 { 556 G4ThreeVector newpos = kt->GetPosition() + 557 G4ThreeVector newpos = kt->GetPosition() + 557 timeStep*c_light/kt->GetTrackingMomen 558 timeStep*c_light/kt->GetTrackingMomentum().e() * kt->GetTrackingMomentum().vect(); 558 kt->SetPosition(newpos); 559 kt->SetPosition(newpos); 559 return true; 560 return true; 560 } 561 } 561 562 562 /* 563 /* 563 G4bool G4RKPropagation::WillBeCaptured(const G 564 G4bool G4RKPropagation::WillBeCaptured(const G4KineticTrack * kt) 564 { 565 { 565 G4double radius = theOuterRadius; 566 G4double radius = theOuterRadius; 566 567 567 // evaluate the final energy. Il will be captu 568 // evaluate the final energy. Il will be captured if newE or newP < 0 568 G4ParticleDefinition * definition = kt->GetD 569 G4ParticleDefinition * definition = kt->GetDefinition(); 569 G4double mass = definition->GetPDGMass(); 570 G4double mass = definition->GetPDGMass(); 570 G4ThreeVector pos = kt->GetPosition(); 571 G4ThreeVector pos = kt->GetPosition(); 571 G4LorentzVector mom = kt->GetTrackingMomentu 572 G4LorentzVector mom = kt->GetTrackingMomentum(); 572 G4VNuclearField * field = (*theFieldMap)[def 573 G4VNuclearField * field = (*theFieldMap)[definition->GetPDGEncoding()]; 573 G4ThreeVector newPos(0, 0, radius); // to ge 574 G4ThreeVector newPos(0, 0, radius); // to get the field on the surface 574 575 575 G4double newE = mom.e()+field->GetField(pos) 576 G4double newE = mom.e()+field->GetField(pos)-field->GetField(newPos); 576 577 577 return ((newE < mass) ? false : true); 578 return ((newE < mass) ? false : true); 578 } 579 } 579 */ 580 */ 580 581 581 582 582 583 583 //-------------------------------------------- 584 //---------------------------------------------------------------------------- 584 G4bool G4RKPropagation::GetSphereIntersectionT 585 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4double radius, 585 //-------------------------------------- 586 //---------------------------------------------------------------------------- 586 const G4ThreeVector & currentPos, 587 const G4ThreeVector & currentPos, 587 const G4LorentzVector & momentum, 588 const G4LorentzVector & momentum, 588 G4double & t1, G4double & t2) 589 G4double & t1, G4double & t2) 589 { 590 { 590 G4ThreeVector speed = momentum.vect()/momen 591 G4ThreeVector speed = momentum.vect()/momentum.e(); // boost vector 591 G4double scalarProd = currentPos.dot(speed) 592 G4double scalarProd = currentPos.dot(speed); 592 G4double speedMag2 = speed.mag2(); 593 G4double speedMag2 = speed.mag2(); 593 G4double sqrtArg = scalarProd*scalarProd - 594 G4double sqrtArg = scalarProd*scalarProd - 594 speedMag2*(currentPos.mag2()-radius*r 595 speedMag2*(currentPos.mag2()-radius*radius); 595 if(sqrtArg <= 0.) // particle will not inte 596 if(sqrtArg <= 0.) // particle will not intersect the sphere 596 { 597 { 597 // G4cout << " GetSphereIntersection 598 // G4cout << " GetSphereIntersectionTimes sqrtArg negative: " << sqrtArg << G4endl; 598 return false; 599 return false; 599 } 600 } 600 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe 601 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light; 601 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe 602 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light; 602 return true; 603 return true; 603 } 604 } 604 605 605 //-------------------------------------------- 606 //---------------------------------------------------------------------------- 606 G4bool G4RKPropagation::GetSphereIntersectionT 607 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4KineticTrack * kt, 607 G4double & t1, G4double & t2) 608 G4double & t1, G4double & t2) 608 { 609 { 609 G4double radius = theOuterRadius + 3*fermi; 610 G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi 610 G4ThreeVector speed = kt->GetTrackingMoment 611 G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector 611 G4double scalarProd = kt->GetPosition().dot 612 G4double scalarProd = kt->GetPosition().dot(speed); 612 G4double speedMag2 = speed.mag2(); 613 G4double speedMag2 = speed.mag2(); 613 G4double sqrtArg = scalarProd*scalarProd - 614 G4double sqrtArg = scalarProd*scalarProd - 614 speedMag2*(kt->GetPosition().mag2()-r 615 speedMag2*(kt->GetPosition().mag2()-radius*radius); 615 if(sqrtArg <= 0.) // particle will not inte 616 if(sqrtArg <= 0.) // particle will not intersect the sphere 616 { 617 { 617 return false; 618 return false; 618 } 619 } 619 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe 620 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light; 620 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe 621 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light; 621 return true; 622 return true; 622 } 623 } 623 624 624 // Implementation methods 625 // Implementation methods 625 626 626 //-------------------------------------------- 627 //---------------------------------------------------------------------------- 627 void G4RKPropagation::delete_FieldsAndMap( 628 void G4RKPropagation::delete_FieldsAndMap( 628 //-------------------------------------- 629 //---------------------------------------------------------------------------- 629 std::map <G4int, G4VNuclearField *, std: 630 std::map <G4int, G4VNuclearField *, std::less<G4int> > * aMap) 630 { 631 { 631 if(aMap) 632 if(aMap) 632 { 633 { 633 std::map <G4int, G4VNuclearField *, std: 634 std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur; 634 for(cur = aMap->begin(); cur != aMap->en 635 for(cur = aMap->begin(); cur != aMap->end(); ++cur) 635 delete (*cur).second; 636 delete (*cur).second; 636 637 637 aMap->clear(); 638 aMap->clear(); 638 delete aMap; 639 delete aMap; 639 } 640 } 640 641 641 } 642 } 642 643 643 //-------------------------------------------- 644 //---------------------------------------------------------------------------- 644 void G4RKPropagation::delete_EquationsAndMap( 645 void G4RKPropagation::delete_EquationsAndMap( 645 //-------------------------------------- 646 //---------------------------------------------------------------------------- 646 std::map <G4int, G4Mag_EqRhs *, std::les 647 std::map <G4int, G4Mag_EqRhs *, std::less<G4int> > * aMap) 647 { 648 { 648 if(aMap) 649 if(aMap) 649 { 650 { 650 std::map <G4int, G4Mag_EqRhs *, std::les 651 std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur; 651 for(cur = aMap->begin(); cur != aMap->en 652 for(cur = aMap->begin(); cur != aMap->end(); ++cur) 652 delete (*cur).second; 653 delete (*cur).second; 653 654 654 aMap->clear(); 655 aMap->clear(); 655 delete aMap; 656 delete aMap; 656 } 657 } 657 } 658 } 658 659