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