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