| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file field/BlineTracer/src/G4BlineTracer.cc
27 /// \brief Implementation of the G4BlineTracer class
28 //
29 //
30 //
31 //
32 // --------------------------------------------------------------------
33 //
34 // G4BlineTracer implementation
35 //
36 // --------------------------------------------------------------------
37 // Author: Laurent Desorgher (desorgher@phim.unibe.ch)
38 // Created - 2003-10-06
39 // --------------------------------------------------------------------
40
41 #include "G4BlineTracer.hh"
42
43 #include "G4BlineEquation.hh"
44 #include "G4BlineEventAction.hh"
45 #include "G4BlinePrimaryGeneratorAction.hh"
46 #include "G4BlineSteppingAction.hh"
47 #include "G4BlineTracerMessenger.hh"
48 #include "G4CashKarpRKF45.hh"
49 #include "G4ChordFinder.hh"
50 #include "G4FieldManager.hh"
51 #include "G4LogicalVolumeStore.hh"
52 #include "G4MagIntegratorDriver.hh"
53 #include "G4PropagatorInField.hh"
54 #include "G4RunManager.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4TransportationManager.hh"
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60 G4BlineTracer::G4BlineTracer()
61 {
62 fMessenger = new G4BlineTracerMessenger(this);
63 fSteppingAction = new G4BlineSteppingAction(this);
64 fEventAction = new G4BlineEventAction(this);
65 fPrimaryGeneratorAction = new G4BlinePrimaryGeneratorAction();
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70 G4BlineTracer::~G4BlineTracer()
71 {
72 delete fMessenger;
73 delete fSteppingAction;
74 delete fEventAction;
75 delete fPrimaryGeneratorAction;
76 for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
77 if (fVecEquationOfMotion[i]) delete fVecEquationOfMotion[i];
78 if (fVecChordFinders[i]) delete fVecChordFinders[i];
79 }
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84 void G4BlineTracer::BeginOfRunAction(const G4Run*) {}
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 void G4BlineTracer::EndOfRunAction(const G4Run*) {}
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91
92 void G4BlineTracer::ComputeBlines(G4int n_of_lines)
93 {
94 // the first time ResetChordFinders should be called
95 //
96 if (!fWas_ResetChordFinders_already_called) {
97 ResetChordFinders();
98 fWas_ResetChordFinders_already_called = true;
99 }
100
101 // Replace the user action by the ad-hoc actions for Blines
102
103 G4RunManager* theRunManager = G4RunManager::GetRunManager();
104 auto user_run_action = (G4UserRunAction*)theRunManager->GetUserRunAction();
105 theRunManager->SetUserAction(this);
106
107 auto user_stepping_action = (G4UserSteppingAction*)theRunManager->GetUserSteppingAction();
108 theRunManager->SetUserAction(fSteppingAction);
109
110 auto userPrimaryAction =
111 (G4VUserPrimaryGeneratorAction*)theRunManager->GetUserPrimaryGeneratorAction();
112 if (userPrimaryAction) fPrimaryGeneratorAction->SetUserPrimaryAction(userPrimaryAction);
113 theRunManager->SetUserAction(fPrimaryGeneratorAction);
114
115 auto user_event_action = (G4UserEventAction*)theRunManager->GetUserEventAction();
116 theRunManager->SetUserAction(fEventAction);
117
118 auto user_tracking_action = (G4UserTrackingAction*)theRunManager->GetUserTrackingAction();
119 G4UserTrackingAction* aNullTrackingAction = nullptr;
120 theRunManager->SetUserAction(aNullTrackingAction);
121
122 auto user_stacking_action = (G4UserStackingAction*)theRunManager->GetUserStackingAction();
123 G4UserStackingAction* aNullStackingAction = nullptr;
124 theRunManager->SetUserAction(aNullStackingAction);
125
126 // replace the user defined chordfinder by the element of fVecChordFinders
127
128 std::vector<G4ChordFinder*> user_chord_finders;
129 std::vector<G4double> user_largest_acceptable_step;
130 for (size_t i = 0; i < fVecChordFinders.size(); i++) {
131 user_largest_acceptable_step.push_back(-1.);
132 if (fVecChordFinders[i]) {
133 user_chord_finders.push_back(fVecFieldManagers[i]->GetChordFinder());
134 fVecChordFinders[i]->SetDeltaChord(user_chord_finders[i]->GetDeltaChord());
135 fVecFieldManagers[i]->SetChordFinder(fVecChordFinders[i]);
136 }
137 else
138 user_chord_finders.push_back(nullptr);
139 }
140
141 // I have tried to use the smooth line filter ability but I could not obtain
142 // a smooth trajectory in the G4TrajectoryContainer after an event
143 // Another solution for obtaining a smooth trajectory is to limit
144 // the LargestAcceptableStep in the G4PropagatorInField object.
145 // This is the solution I used.
146
147 // Old solution:
148 // G4TransportationManager::GetTransportationManager()
149 // ->GetPropagatorInField()->SetTrajectoryFilter(fTrajectoryFilter);
150
151 // New solution:
152 // set the largest_acceptable_step to max_step:length
153
154 G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();
155 G4double previous_largest_acceptable_step =
156 tmanager->GetPropagatorInField()->GetLargestAcceptableStep();
157
158 tmanager->GetPropagatorInField()->SetLargestAcceptableStep(fMaxTrackingStep);
159
160 // Start the integration of n_of_lines different magnetic field lines
161
162 for (G4int il = 0; il < n_of_lines; il++) {
163 // for each magnetic field line we integrate once backward and once
164 // forward from the same starting point
165
166 // backward integration
167
168 for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
169 if (fVecEquationOfMotion[i]) fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(true);
170 }
171 theRunManager->BeamOn(1);
172
173 // forward integration
174
175 for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
176 if (fVecEquationOfMotion[i])
177 fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(false);
178 }
179 theRunManager->BeamOn(1);
180 }
181
182 // Remove trajectory filter to PropagatorInField
183 // It was for old solution when using smooth trajectory filter
184
185 // tmanager->GetPropagatorInField()->SetTrajectoryFilter(0);
186
187 // back to User defined actions and other parameters
188 // -------------------------------------------------
189
190 tmanager->GetPropagatorInField()->SetLargestAcceptableStep(previous_largest_acceptable_step);
191
192 // return to User actions
193
194 theRunManager->SetUserAction(user_run_action);
195 theRunManager->SetUserAction(user_event_action);
196 theRunManager->SetUserAction(userPrimaryAction);
197 theRunManager->SetUserAction(user_stepping_action);
198 theRunManager->SetUserAction(user_tracking_action);
199 theRunManager->SetUserAction(user_stacking_action);
200
201 // set user defined chord finders and largest acceptable step
202
203 for (size_t i = 0; i < fVecFieldManagers.size(); i++) {
204 if (user_chord_finders[i]) fVecFieldManagers[i]->SetChordFinder(user_chord_finders[i]);
205 }
206 }
207
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 ////////////////////////////////////////////////////////////////
210
211 /*
212 G4bool G4BlineTracer::CheckMagneticFields()
213 {
214 // Check FieldManagers
215
216 G4TransportationManager* tmanager =
217 G4TransportationManager::GetTransportationManager();
218
219 if (fVecFieldManagers[0] != tmanager->GetFieldManager())
220 return false;
221 if (fVecMagneticFields[0] != tmanager->GetFieldManager()->GetDetectorField())
222 return false;
223 G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
224
225 std::vector<G4FieldManagers*> LogicalVolumeFields;
226 size_t j=0;
227 for (size_t i=0; i<theVolumeStore.size();i++)
228 {
229 if (theVolumeStore[i]->GetFieldManager())
230 {
231 j++;
232 if (j >= fVecFieldManagers.size()) return false;
233 if (fVecFieldManagers[j] != theVolumeStore[i]->GetFieldManager())
234 return false;
235 if (fVecMagneticFields[j] !=
236 theVolumeStore[i]->GetFieldManager()->GetDetectorField())
237 return false;
238 }
239 }
240 if (j<fVecFieldManagers.size()) return false;
241
242 return true;
243 }
244 */
245
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247
248 void G4BlineTracer::ResetChordFinders()
249 {
250 for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
251 delete fVecEquationOfMotion[i];
252 delete fVecChordFinders[i];
253 }
254
255 fVecChordFinders.clear();
256 fVecFieldManagers.clear();
257 fVecMagneticFields.clear();
258 fVecEquationOfMotion.clear();
259
260 // global field
261
262 fVecChordFinders.push_back(nullptr);
263 fVecMagneticFields.push_back(nullptr);
264 fVecEquationOfMotion.push_back(nullptr);
265 fVecFieldManagers.push_back(
266 G4TransportationManager::GetTransportationManager()->GetFieldManager());
267 if (fVecFieldManagers[0]) {
268 fVecMagneticFields[0] = (G4MagneticField*)fVecFieldManagers[0]->GetDetectorField();
269 if (fVecMagneticFields[0]) {
270 fVecEquationOfMotion[0] = new G4BlineEquation(fVecMagneticFields[0]);
271 auto pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[0]);
272 auto pIntgrDriver =
273 new G4MagInt_Driver(0.01 * mm, pStepper, pStepper->GetNumberOfVariables());
274 fVecChordFinders[0] = new G4ChordFinder(pIntgrDriver);
275 }
276 }
277
278 // local fields
279
280 G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
281
282 size_t j = 0;
283 for (size_t i = 0; i < theVolumeStore->size(); i++) {
284 if ((*theVolumeStore)[i]->GetFieldManager()) {
285 j++;
286 fVecFieldManagers.push_back(((*theVolumeStore)[i])->GetFieldManager());
287 fVecMagneticFields.push_back((G4MagneticField*)fVecFieldManagers[j]->GetDetectorField());
288 fVecEquationOfMotion.push_back(nullptr);
289 fVecChordFinders.push_back(nullptr);
290 if (fVecMagneticFields[j]) {
291 fVecEquationOfMotion[j] = new G4BlineEquation(fVecMagneticFields[j]);
292 auto pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[j]);
293 auto pIntgrDriver =
294 new G4MagInt_Driver(.01 * mm, pStepper, pStepper->GetNumberOfVariables());
295 fVecChordFinders[j] = new G4ChordFinder(pIntgrDriver);
296 }
297 }
298 }
299 }
300