| 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 /// \file field/BlineTracer/src/G4BlineTracer. 26 /// \file field/BlineTracer/src/G4BlineTracer.cc
27 /// \brief Implementation of the G4BlineTracer 27 /// \brief Implementation of the G4BlineTracer class
28 // 28 //
29 // 29 //
30 // 30 //
31 // 31 //
32 // ------------------------------------------- 32 // --------------------------------------------------------------------
33 // 33 //
34 // G4BlineTracer implementation 34 // G4BlineTracer implementation
35 // 35 //
36 // ------------------------------------------- 36 // --------------------------------------------------------------------
37 // Author: Laurent Desorgher (desorgher@phim.u 37 // Author: Laurent Desorgher (desorgher@phim.unibe.ch)
38 // Created - 2003-10-06 38 // Created - 2003-10-06
39 // ------------------------------------------- 39 // --------------------------------------------------------------------
40 40
41 #include "G4BlineTracer.hh" 41 #include "G4BlineTracer.hh"
42 << 42 #include "G4BlineTracerMessenger.hh"
43 #include "G4BlineEquation.hh" <<
44 #include "G4BlineEventAction.hh" <<
45 #include "G4BlinePrimaryGeneratorAction.hh" 43 #include "G4BlinePrimaryGeneratorAction.hh"
>> 44 #include "G4BlineEventAction.hh"
46 #include "G4BlineSteppingAction.hh" 45 #include "G4BlineSteppingAction.hh"
47 #include "G4BlineTracerMessenger.hh" << 46 #include "G4BlineEquation.hh"
48 #include "G4CashKarpRKF45.hh" << 47
49 #include "G4ChordFinder.hh" << 48 #include "G4RunManager.hh"
>> 49 #include "G4TransportationManager.hh"
50 #include "G4FieldManager.hh" 50 #include "G4FieldManager.hh"
51 #include "G4LogicalVolumeStore.hh" <<
52 #include "G4MagIntegratorDriver.hh" <<
53 #include "G4PropagatorInField.hh" 51 #include "G4PropagatorInField.hh"
54 #include "G4RunManager.hh" << 52 #include "G4MagIntegratorDriver.hh"
>> 53 #include "G4CashKarpRKF45.hh"
>> 54 #include "G4LogicalVolumeStore.hh"
>> 55 #include "G4ChordFinder.hh"
55 #include "G4SystemOfUnits.hh" 56 #include "G4SystemOfUnits.hh"
56 #include "G4TransportationManager.hh" <<
57 57
58 //....oooOO0OOooo........oooOO0OOooo........oo << 58 //////////////////////////////////////////////////////////////////
59 59
60 G4BlineTracer::G4BlineTracer() 60 G4BlineTracer::G4BlineTracer()
61 { 61 {
62 fMessenger = new G4BlineTracerMessenger(this 62 fMessenger = new G4BlineTracerMessenger(this);
63 fSteppingAction = new G4BlineSteppingAction( << 63 fSteppingAction = new G4BlineSteppingAction(this) ;
64 fEventAction = new G4BlineEventAction(this); 64 fEventAction = new G4BlineEventAction(this);
65 fPrimaryGeneratorAction = new G4BlinePrimary 65 fPrimaryGeneratorAction = new G4BlinePrimaryGeneratorAction();
>> 66 fMaxTrackingStep =1000.*m;
>> 67 fWas_ResetChordFinders_already_called=false;
66 } 68 }
67 69
68 //....oooOO0OOooo........oooOO0OOooo........oo << 70 ///////////////////////////////////////////////////////////////////////
69 71
70 G4BlineTracer::~G4BlineTracer() 72 G4BlineTracer::~G4BlineTracer()
71 { 73 {
72 delete fMessenger; 74 delete fMessenger;
73 delete fSteppingAction; 75 delete fSteppingAction;
74 delete fEventAction; << 76 delete fEventAction;
75 delete fPrimaryGeneratorAction; 77 delete fPrimaryGeneratorAction;
76 for (size_t i = 0; i < fVecEquationOfMotion. << 78 for (size_t i=0; i< fVecEquationOfMotion.size();i++)
>> 79 {
77 if (fVecEquationOfMotion[i]) delete fVecEq 80 if (fVecEquationOfMotion[i]) delete fVecEquationOfMotion[i];
78 if (fVecChordFinders[i]) delete fVecChordF 81 if (fVecChordFinders[i]) delete fVecChordFinders[i];
79 } 82 }
80 } << 83 }
81 84
82 //....oooOO0OOooo........oooOO0OOooo........oo << 85 ////////////////////////////////////////////////////////////////////
83 86
84 void G4BlineTracer::BeginOfRunAction(const G4R << 87 void G4BlineTracer::BeginOfRunAction(const G4Run*)
>> 88 {
>> 89 }
85 90
86 //....oooOO0OOooo........oooOO0OOooo........oo << 91 ///////////////////////////////////////////////////////////////////////
87 92
88 void G4BlineTracer::EndOfRunAction(const G4Run << 93 void G4BlineTracer::EndOfRunAction(const G4Run*)
>> 94 {
>> 95 }
89 96
90 //....oooOO0OOooo........oooOO0OOooo........oo << 97 ////////////////////////////////////////////////////////////////
91 98
92 void G4BlineTracer::ComputeBlines(G4int n_of_l 99 void G4BlineTracer::ComputeBlines(G4int n_of_lines)
93 { 100 {
94 // the first time ResetChordFinders should b << 101 //the first time ResetChordFinders should be called
95 // 102 //
96 if (!fWas_ResetChordFinders_already_called) << 103 if (!fWas_ResetChordFinders_already_called)
>> 104 {
97 ResetChordFinders(); 105 ResetChordFinders();
98 fWas_ResetChordFinders_already_called = tr << 106 fWas_ResetChordFinders_already_called=true;
99 } 107 }
100 108
101 // Replace the user action by the ad-hoc act 109 // Replace the user action by the ad-hoc actions for Blines
102 << 110
103 G4RunManager* theRunManager = G4RunManager:: << 111 G4RunManager* theRunManager = G4RunManager::GetRunManager();
104 auto user_run_action = (G4UserRunAction*)the << 112 G4UserRunAction* user_run_action =
>> 113 (G4UserRunAction*)theRunManager->GetUserRunAction();
105 theRunManager->SetUserAction(this); 114 theRunManager->SetUserAction(this);
106 115
107 auto user_stepping_action = (G4UserSteppingA << 116 G4UserSteppingAction* user_stepping_action =
>> 117 (G4UserSteppingAction*)theRunManager->GetUserSteppingAction();
108 theRunManager->SetUserAction(fSteppingAction 118 theRunManager->SetUserAction(fSteppingAction);
109 << 119
110 auto userPrimaryAction = << 120 G4VUserPrimaryGeneratorAction* userPrimaryAction =
111 (G4VUserPrimaryGeneratorAction*)theRunMana 121 (G4VUserPrimaryGeneratorAction*)theRunManager->GetUserPrimaryGeneratorAction();
112 if (userPrimaryAction) fPrimaryGeneratorActi << 122 if (userPrimaryAction)
>> 123 fPrimaryGeneratorAction->SetUserPrimaryAction(userPrimaryAction);
113 theRunManager->SetUserAction(fPrimaryGenerat 124 theRunManager->SetUserAction(fPrimaryGeneratorAction);
114 125
115 auto user_event_action = (G4UserEventAction* << 126 G4UserEventAction* user_event_action =
>> 127 (G4UserEventAction*)theRunManager->GetUserEventAction();
116 theRunManager->SetUserAction(fEventAction); 128 theRunManager->SetUserAction(fEventAction);
117 << 129
118 auto user_tracking_action = (G4UserTrackingA << 130 G4UserTrackingAction* user_tracking_action =
119 G4UserTrackingAction* aNullTrackingAction = << 131 (G4UserTrackingAction*)theRunManager->GetUserTrackingAction();
>> 132 G4UserTrackingAction* aNullTrackingAction = 0;
120 theRunManager->SetUserAction(aNullTrackingAc 133 theRunManager->SetUserAction(aNullTrackingAction);
121 134
122 auto user_stacking_action = (G4UserStackingA << 135 G4UserStackingAction* user_stacking_action =
123 G4UserStackingAction* aNullStackingAction = << 136 (G4UserStackingAction*)theRunManager->GetUserStackingAction();
>> 137 G4UserStackingAction* aNullStackingAction = 0;
124 theRunManager->SetUserAction(aNullStackingAc 138 theRunManager->SetUserAction(aNullStackingAction);
125 139
126 // replace the user defined chordfinder by t << 140 // replace the user defined chordfinder by the element of fVecChordFinders
127 << 141
128 std::vector<G4ChordFinder*> user_chord_finde 142 std::vector<G4ChordFinder*> user_chord_finders;
129 std::vector<G4double> user_largest_acceptabl 143 std::vector<G4double> user_largest_acceptable_step;
130 for (size_t i = 0; i < fVecChordFinders.size << 144 for (size_t i=0;i<fVecChordFinders.size();i++)
>> 145 {
131 user_largest_acceptable_step.push_back(-1. 146 user_largest_acceptable_step.push_back(-1.);
132 if (fVecChordFinders[i]) { << 147 if (fVecChordFinders[i])
>> 148 {
133 user_chord_finders.push_back(fVecFieldMa 149 user_chord_finders.push_back(fVecFieldManagers[i]->GetChordFinder());
134 fVecChordFinders[i]->SetDeltaChord(user_ 150 fVecChordFinders[i]->SetDeltaChord(user_chord_finders[i]->GetDeltaChord());
135 fVecFieldManagers[i]->SetChordFinder(fVe 151 fVecFieldManagers[i]->SetChordFinder(fVecChordFinders[i]);
136 } 152 }
137 else << 153 else user_chord_finders.push_back(0);
138 user_chord_finders.push_back(nullptr); <<
139 } 154 }
140 155
141 // I have tried to use the smooth line filte << 156 // I have tried to use the smooth line filter ability but I could not obtain
142 // a smooth trajectory in the G4TrajectoryCo 157 // a smooth trajectory in the G4TrajectoryContainer after an event
143 // Another solution for obtaining a smooth t 158 // Another solution for obtaining a smooth trajectory is to limit
144 // the LargestAcceptableStep in the G4Propag 159 // the LargestAcceptableStep in the G4PropagatorInField object.
145 // This is the solution I used. << 160 // This is the solution I used.
146 161
147 // Old solution: 162 // Old solution:
148 // G4TransportationManager::GetTransportatio 163 // G4TransportationManager::GetTransportationManager()
149 // ->GetPropagatorInField()->SetTrajecto 164 // ->GetPropagatorInField()->SetTrajectoryFilter(fTrajectoryFilter);
150 165
151 // New solution: 166 // New solution:
152 // set the largest_acceptable_step to max_st 167 // set the largest_acceptable_step to max_step:length
153 168
154 G4TransportationManager* tmanager = G4Transp << 169 G4TransportationManager* tmanager =
>> 170 G4TransportationManager::GetTransportationManager();
155 G4double previous_largest_acceptable_step = 171 G4double previous_largest_acceptable_step =
156 tmanager->GetPropagatorInField()->GetLarge 172 tmanager->GetPropagatorInField()->GetLargestAcceptableStep();
157 173
158 tmanager->GetPropagatorInField()->SetLargest << 174 tmanager->GetPropagatorInField()
>> 175 ->SetLargestAcceptableStep(fMaxTrackingStep);
159 176
160 // Start the integration of n_of_lines diffe 177 // Start the integration of n_of_lines different magnetic field lines
161 178
162 for (G4int il = 0; il < n_of_lines; il++) { << 179 for (G4int il=0; il<n_of_lines;il++)
>> 180 {
163 // for each magnetic field line we integra 181 // for each magnetic field line we integrate once backward and once
164 // forward from the same starting point 182 // forward from the same starting point
165 183
166 // backward integration 184 // backward integration
167 185
168 for (size_t i = 0; i < fVecEquationOfMotio << 186 for (size_t i=0; i< fVecEquationOfMotion.size();i++)
169 if (fVecEquationOfMotion[i]) fVecEquatio << 187 {
>> 188 if (fVecEquationOfMotion[i])
>> 189 fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(true);
170 } 190 }
171 theRunManager->BeamOn(1); 191 theRunManager->BeamOn(1);
172 192
173 // forward integration 193 // forward integration
174 194
175 for (size_t i = 0; i < fVecEquationOfMotio << 195 for (size_t i=0; i < fVecEquationOfMotion.size();i++)
176 if (fVecEquationOfMotion[i]) << 196 {
>> 197 if (fVecEquationOfMotion[i])
177 fVecEquationOfMotion[i]->SetBackwardDi 198 fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(false);
178 } 199 }
179 theRunManager->BeamOn(1); 200 theRunManager->BeamOn(1);
180 } 201 }
181 202
182 // Remove trajectory filter to PropagatorInF 203 // Remove trajectory filter to PropagatorInField
183 // It was for old solution when using smooth 204 // It was for old solution when using smooth trajectory filter
184 205
185 // tmanager->GetPropagatorInField()->SetTraj 206 // tmanager->GetPropagatorInField()->SetTrajectoryFilter(0);
186 207
187 // back to User defined actions and other pa 208 // back to User defined actions and other parameters
188 // ----------------------------------------- 209 // -------------------------------------------------
189 210
190 tmanager->GetPropagatorInField()->SetLargest << 211 tmanager->GetPropagatorInField()
>> 212 ->SetLargestAcceptableStep(previous_largest_acceptable_step);
191 213
192 // return to User actions 214 // return to User actions
193 215
194 theRunManager->SetUserAction(user_run_action 216 theRunManager->SetUserAction(user_run_action);
195 theRunManager->SetUserAction(user_event_acti 217 theRunManager->SetUserAction(user_event_action);
196 theRunManager->SetUserAction(userPrimaryActi 218 theRunManager->SetUserAction(userPrimaryAction);
197 theRunManager->SetUserAction(user_stepping_a 219 theRunManager->SetUserAction(user_stepping_action);
198 theRunManager->SetUserAction(user_tracking_a 220 theRunManager->SetUserAction(user_tracking_action);
199 theRunManager->SetUserAction(user_stacking_a 221 theRunManager->SetUserAction(user_stacking_action);
200 222
201 // set user defined chord finders and larges 223 // set user defined chord finders and largest acceptable step
202 224
203 for (size_t i = 0; i < fVecFieldManagers.siz << 225 for (size_t i=0;i<fVecFieldManagers.size();i++)
204 if (user_chord_finders[i]) fVecFieldManage << 226 {
205 } << 227 if (user_chord_finders[i])
>> 228 fVecFieldManagers[i]->SetChordFinder(user_chord_finders[i]);
>> 229 }
206 } 230 }
207 231
208 //....oooOO0OOooo........oooOO0OOooo........oo <<
209 ////////////////////////////////////////////// 232 ////////////////////////////////////////////////////////////////
210 233
211 /* 234 /*
212 G4bool G4BlineTracer::CheckMagneticFields() 235 G4bool G4BlineTracer::CheckMagneticFields()
213 { 236 {
214 // Check FieldManagers 237 // Check FieldManagers
215 238
216 G4TransportationManager* tmanager = 239 G4TransportationManager* tmanager =
217 G4TransportationManager::GetTransportation 240 G4TransportationManager::GetTransportationManager();
218 241
219 if (fVecFieldManagers[0] != tmanager->GetFie 242 if (fVecFieldManagers[0] != tmanager->GetFieldManager())
220 return false; 243 return false;
221 if (fVecMagneticFields[0] != tmanager->GetFi 244 if (fVecMagneticFields[0] != tmanager->GetFieldManager()->GetDetectorField())
222 return false; 245 return false;
223 G4LogicalVolumeStore* theVolumeStore = G4Log << 246 G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
224 << 247
225 std::vector<G4FieldManagers*> LogicalVolumeF 248 std::vector<G4FieldManagers*> LogicalVolumeFields;
226 size_t j=0; 249 size_t j=0;
227 for (size_t i=0; i<theVolumeStore.size();i++ 250 for (size_t i=0; i<theVolumeStore.size();i++)
228 { 251 {
229 if (theVolumeStore[i]->GetFieldManager()) 252 if (theVolumeStore[i]->GetFieldManager())
230 { 253 {
231 j++; 254 j++;
232 if (j >= fVecFieldManagers.size()) retur 255 if (j >= fVecFieldManagers.size()) return false;
233 if (fVecFieldManagers[j] != theVolumeSto 256 if (fVecFieldManagers[j] != theVolumeStore[i]->GetFieldManager())
234 return false; 257 return false;
235 if (fVecMagneticFields[j] != 258 if (fVecMagneticFields[j] !=
236 theVolumeStore[i]->GetFieldManager() 259 theVolumeStore[i]->GetFieldManager()->GetDetectorField())
237 return false; 260 return false;
238 } 261 }
239 } 262 }
240 if (j<fVecFieldManagers.size()) return false 263 if (j<fVecFieldManagers.size()) return false;
241 264
242 return true; 265 return true;
243 } 266 }
244 */ 267 */
245 268
246 //....oooOO0OOooo........oooOO0OOooo........oo << 269 ////////////////////////////////////////////////////////////////
247 270
248 void G4BlineTracer::ResetChordFinders() 271 void G4BlineTracer::ResetChordFinders()
249 { 272 {
250 for (size_t i = 0; i < fVecEquationOfMotion. << 273 for (size_t i=0; i<fVecEquationOfMotion.size();i++)
>> 274 {
251 delete fVecEquationOfMotion[i]; 275 delete fVecEquationOfMotion[i];
252 delete fVecChordFinders[i]; 276 delete fVecChordFinders[i];
253 } << 277 }
254 278
255 fVecChordFinders.clear(); 279 fVecChordFinders.clear();
256 fVecFieldManagers.clear(); 280 fVecFieldManagers.clear();
257 fVecMagneticFields.clear(); 281 fVecMagneticFields.clear();
258 fVecEquationOfMotion.clear(); 282 fVecEquationOfMotion.clear();
259 283
260 // global field 284 // global field
261 285
262 fVecChordFinders.push_back(nullptr); << 286 fVecChordFinders.push_back(0);
263 fVecMagneticFields.push_back(nullptr); << 287 fVecMagneticFields.push_back(0);
264 fVecEquationOfMotion.push_back(nullptr); << 288 fVecEquationOfMotion.push_back(0);
265 fVecFieldManagers.push_back( << 289 fVecFieldManagers.push_back(G4TransportationManager::GetTransportationManager()
266 G4TransportationManager::GetTransportation << 290 ->GetFieldManager());
267 if (fVecFieldManagers[0]) { << 291 if (fVecFieldManagers[0])
268 fVecMagneticFields[0] = (G4MagneticField*) << 292 {
269 if (fVecMagneticFields[0]) { << 293 fVecMagneticFields[0] =
>> 294 (G4MagneticField*) fVecFieldManagers[0]->GetDetectorField();
>> 295 if (fVecMagneticFields[0])
>> 296 {
270 fVecEquationOfMotion[0] = new G4BlineEqu 297 fVecEquationOfMotion[0] = new G4BlineEquation(fVecMagneticFields[0]);
271 auto pStepper = new G4CashKarpRKF45(fVec << 298 G4CashKarpRKF45* pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[0]);
272 auto pIntgrDriver = << 299 G4MagInt_Driver* pIntgrDriver =
273 new G4MagInt_Driver(0.01 * mm, pSteppe << 300 new G4MagInt_Driver(0.01*mm,pStepper,pStepper->GetNumberOfVariables());
274 fVecChordFinders[0] = new G4ChordFinder( 301 fVecChordFinders[0] = new G4ChordFinder(pIntgrDriver);
275 } 302 }
276 } << 303 }
277 304
278 // local fields << 305 // local fields
279 306
280 G4LogicalVolumeStore* theVolumeStore = G4Log 307 G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
281 308
282 size_t j = 0; << 309 size_t j=0;
283 for (size_t i = 0; i < theVolumeStore->size( << 310 for (size_t i=0; i<theVolumeStore->size();i++)
284 if ((*theVolumeStore)[i]->GetFieldManager( << 311 {
>> 312 if ((*theVolumeStore)[i]->GetFieldManager())
>> 313 {
285 j++; 314 j++;
286 fVecFieldManagers.push_back(((*theVolume 315 fVecFieldManagers.push_back(((*theVolumeStore)[i])->GetFieldManager());
287 fVecMagneticFields.push_back((G4Magnetic << 316 fVecMagneticFields.push_back((G4MagneticField*)
288 fVecEquationOfMotion.push_back(nullptr); << 317 fVecFieldManagers[j]->GetDetectorField());
289 fVecChordFinders.push_back(nullptr); << 318 fVecEquationOfMotion.push_back(0);
290 if (fVecMagneticFields[j]) { << 319 fVecChordFinders.push_back(0);
291 fVecEquationOfMotion[j] = new G4BlineE << 320 if (fVecMagneticFields[j])
292 auto pStepper = new G4CashKarpRKF45(fV << 321 {
293 auto pIntgrDriver = << 322 fVecEquationOfMotion[j]= new G4BlineEquation(fVecMagneticFields[j]);
294 new G4MagInt_Driver(.01 * mm, pStepp << 323 G4CashKarpRKF45* pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[j]);
>> 324 G4MagInt_Driver* pIntgrDriver =
>> 325 new G4MagInt_Driver(.01*mm,pStepper,pStepper->GetNumberOfVariables());
295 fVecChordFinders[j] = new G4ChordFinde 326 fVecChordFinders[j] = new G4ChordFinder(pIntgrDriver);
296 } << 327 }
297 } 328 }
298 } << 329 }
299 } 330 }
300 331