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