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