source: trunk/examples/extended/field/BlineTracer/src/G4BlineTracer.cc@ 1036

Last change on this file since 1036 was 807, checked in by garnier, 17 years ago

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