source: trunk/source/geometry/magneticfield/test/testPropagateMagField.cc@ 1325

Last change on this file since 1325 was 1231, checked in by garnier, 16 years ago

update from CVS
File size: 22.4 KB
Line 
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: testPropagateMagField.cc,v 1.33 2009/11/05 13:18:05 japost Exp $
28// GEANT4 tag $Name: $
29//
30//
31//
32// Started from testG4Navigator1.cc,v 1.7 1996/08/29 15:42 pkent
33// Locate & Step within simple boxlike geometry, both
34// with and without voxels. Parameterised volumes are included.
35
36#include <assert.h>
37// #include "ApproxEqual.hh"
38
39// Global defs
40#include "globals.hh"
41
42#include "G4Navigator.hh"
43
44#include "G4LogicalVolume.hh"
45#include "G4VPhysicalVolume.hh"
46#include "G4PVPlacement.hh"
47#include "G4PVParameterised.hh"
48#include "G4VPVParameterisation.hh"
49#include "G4Box.hh"
50
51#include "G4GeometryManager.hh"
52
53#include "G4RotationMatrix.hh"
54#include "G4ThreeVector.hh"
55
56#include "G4ios.hh"
57#include <iomanip>
58
59// Sample Parameterisation
60class G4LinScale : public G4VPVParameterisation
61{
62 virtual void ComputeTransformation(const G4int n,
63 G4VPhysicalVolume* pRep) const
64 {
65 pRep->SetTranslation(G4ThreeVector(0,(n-1)*15,0));
66 }
67
68 virtual void ComputeDimensions(G4Box &pBox,
69 const G4int n,
70 const G4VPhysicalVolume* ) const
71 {
72 pBox.SetXHalfLength(10);
73 pBox.SetYHalfLength(5+n);
74 pBox.SetZHalfLength(5+n);
75 }
76
77 virtual void ComputeDimensions(G4Tubs &,
78 const G4int ,
79 const G4VPhysicalVolume*) const {}
80 virtual void ComputeDimensions(G4Trd &,
81 const G4int,
82 const G4VPhysicalVolume*) const {}
83 virtual void ComputeDimensions(G4Cons &,
84 const G4int ,
85 const G4VPhysicalVolume*) const {}
86 virtual void ComputeDimensions(G4Trap &,
87 const G4int ,
88 const G4VPhysicalVolume*) const {}
89 virtual void ComputeDimensions(G4Hype &,
90 const G4int ,
91 const G4VPhysicalVolume*) const {}
92 virtual void ComputeDimensions(G4Orb &,
93 const G4int ,
94 const G4VPhysicalVolume*) const {}
95 virtual void ComputeDimensions(G4Sphere &,
96 const G4int ,
97 const G4VPhysicalVolume*) const {}
98 virtual void ComputeDimensions(G4Torus &,
99 const G4int ,
100 const G4VPhysicalVolume*) const {}
101 virtual void ComputeDimensions(G4Para &,
102 const G4int ,
103 const G4VPhysicalVolume*) const {}
104 virtual void ComputeDimensions(G4Polycone &,
105 const G4int ,
106 const G4VPhysicalVolume*) const {}
107 virtual void ComputeDimensions(G4Polyhedra &,
108 const G4int ,
109 const G4VPhysicalVolume*) const {}
110};
111G4LinScale myParam;
112
113// Build simple geometry:
114// 4 small cubes + 1 slab (all G4Boxes) are positioned inside a larger cuboid
115G4VPhysicalVolume* BuildGeometry()
116{
117
118 G4Box *myHugeBox= new G4Box("huge box",15*m,15*m,25*m);
119 G4Box *myBigBox= new G4Box("big cube",10*m,10*m,10*m);
120 G4Box *mySmallBox= new G4Box("smaller cube",2.5*m,2.5*m,2.5*m);
121 G4Box *myTinyBox= new G4Box("tiny cube",.25*m,.25*m,.25*m);
122
123 // G4Box *myVariableBox=
124 new G4Box("Variable Box",10,5,5);
125
126 // World Volume
127 //
128 G4LogicalVolume *worldLog=new G4LogicalVolume(myHugeBox,0,
129 "World",0,0,0);
130 // Logical with no material,field,
131 // sensitive detector or user limits
132
133 G4PVPlacement *worldPhys=new
134 G4PVPlacement(0,G4ThreeVector(0,0,0), "World",worldLog,
135 0,false,0);
136 // Note: no mother pointer set
137
138// Create the logical Volumes
139//
140// G4LogicalVolume(*pSolid, *pMaterial, Name, *pField, *pSDetector, *pULimits)
141//
142 G4LogicalVolume *BigBoxLog=new G4LogicalVolume(myBigBox,0,
143 "Crystal Box (large)",0,0,0);
144 G4LogicalVolume *smallBoxLog=new G4LogicalVolume(mySmallBox,0,
145 "Crystal Box (small)");
146 G4LogicalVolume *tinyBoxLog=new G4LogicalVolume(myTinyBox,0,
147 "Crystal Box (tiny)");
148
149
150// Place them.
151//
152// 1) Two big boxes in the world volume
153//
154 // G4PVPlacement *BigTg1Phys=
155 new G4PVPlacement(0,G4ThreeVector(0,0,-15*m),
156 "Big Target 1",BigBoxLog,
157 worldPhys,false,0);
158 // G4PVPlacement *BigTg2Phys=
159 new G4PVPlacement(0,G4ThreeVector(0,0, 15*m),
160 "Big Target 2",BigBoxLog,
161 worldPhys,false,0);
162
163// 2) Four (medium) boxes in X & Y near the origin of the world volume
164//
165 // G4PVPlacement *MedTg3a_Phys=
166 new G4PVPlacement(0,G4ThreeVector(0, 7.5*m,0),
167 "Target 3a",smallBoxLog,
168 worldPhys,false,0);
169 // G4PVPlacement *MedTg3b_Phys=
170 new G4PVPlacement(0,G4ThreeVector(0,-7.5*m,0),
171 "Target 3b",smallBoxLog,
172 worldPhys,false,0);
173 // G4PVPlacement *MedTg3c_Phys=
174 new G4PVPlacement(0,G4ThreeVector(-7.5*m,0,0),
175 "Target 3c",smallBoxLog,
176 worldPhys,false,0);
177 // G4PVPlacement *MedTg3d_Phys=
178 new G4PVPlacement(0,G4ThreeVector( 7.5*m,0,0),
179 "Target 3d",smallBoxLog,
180 worldPhys,false,0);
181
182
183// 3) Eight small boxes around the origin of the world volume
184// (in +-X, +-Y & +-Z)
185//
186 // G4PVPlacement *SmTg4a_Phys=
187 new G4PVPlacement
188 (0,G4ThreeVector( 0.3*m, 0.3*m,0.3*m), "Target 4a",tinyBoxLog,
189 worldPhys,false,0);
190 // G4PVPlacement *SmTg4b_Phys=
191 new G4PVPlacement
192 (0,G4ThreeVector( 0.3*m,-0.3*m,0.3*m), "Target 4b",tinyBoxLog,
193 worldPhys,false,0);
194 // G4PVPlacement *SmTg4c_Phys=
195 new G4PVPlacement
196 (0,G4ThreeVector(-0.3*m,-0.3*m,0.3*m), "Target 4c",tinyBoxLog,
197 worldPhys,false,0);
198 // G4PVPlacement *SmTg4d_Phys=
199 new G4PVPlacement
200 (0,G4ThreeVector(-0.3*m, 0.3*m,0.3*m), "Target 4d",tinyBoxLog,
201 worldPhys,false,0);
202
203 // G4PVPlacement *SmTg4e_Phys=
204 new G4PVPlacement
205 (0,G4ThreeVector( 0.3*m, 0.3*m,-0.3*m), "Target 4e",tinyBoxLog,
206 worldPhys,false,0);
207 // G4PVPlacement *SmTg4f_Phys=
208 new G4PVPlacement
209 (0,G4ThreeVector( 0.3*m,-0.3*m,-0.3*m), "Target 4f",tinyBoxLog,
210 worldPhys,false,0);
211 // G4PVPlacement *SmTg4g_Phys=
212 new G4PVPlacement
213 (0,G4ThreeVector(-0.3*m,-0.3*m,-0.3*m), "Target 4g",tinyBoxLog,
214 worldPhys,false,0);
215 // G4PVPlacement *SmTg4h_Phys=
216 new G4PVPlacement
217 (0,G4ThreeVector(-0.3*m, 0.3*m,-0.3*m), "Target 4h",tinyBoxLog,
218 worldPhys,false,0);
219
220 return worldPhys;
221}
222
223#include "G4UniformMagField.hh"
224#include "G4QuadrupoleMagField.hh"
225#include "G4CachedMagneticField.hh"
226
227#include "G4ChordFinder.hh"
228#include "G4PropagatorInField.hh"
229#include "G4MagneticField.hh"
230#include "G4FieldManager.hh"
231#include "G4TransportationManager.hh"
232#include "G4HelixExplicitEuler.hh"
233#include "G4HelixSimpleRunge.hh"
234#include "G4HelixImplicitEuler.hh"
235#include "G4ExactHelixStepper.hh"
236#include "G4ExplicitEuler.hh"
237#include "G4ImplicitEuler.hh"
238#include "G4SimpleRunge.hh"
239#include "G4SimpleHeum.hh"
240#include "G4ClassicalRK4.hh"
241#include "G4Mag_UsualEqRhs.hh"
242#include "G4CashKarpRKF45.hh"
243#include "G4RKG3_Stepper.hh"
244#include "G4ConstRK4.hh"
245#include "G4NystromRK4.hh"
246#include "G4HelixMixedStepper.hh"
247
248G4UniformMagField uniformMagField(10.*tesla, 0., 0.);
249G4QuadrupoleMagField quadrupoleMagField( 10.*tesla/(50.*cm) );
250G4CachedMagneticField myMagField( &quadrupoleMagField, 1.0 * cm);
251
252G4FieldManager* SetupField(G4int type)
253{
254 G4FieldManager *pFieldMgr;
255 G4ChordFinder *pChordFinder;
256 G4Mag_UsualEqRhs *fEquation = new G4Mag_UsualEqRhs(&myMagField);
257 G4MagIntegratorStepper *pStepper;
258
259 switch ( type )
260 {
261 case 0: pStepper = new G4ExplicitEuler( fEquation ); break;
262 case 1: pStepper = new G4ImplicitEuler( fEquation ); break;
263 case 2: pStepper = new G4SimpleRunge( fEquation ); break;
264 case 3: pStepper = new G4SimpleHeum( fEquation ); break;
265 case 4: pStepper = new G4ClassicalRK4( fEquation ); break;
266 case 5: pStepper = new G4HelixExplicitEuler( fEquation ); break;
267 case 6: pStepper = new G4HelixImplicitEuler( fEquation ); break;
268 case 7: pStepper = new G4HelixSimpleRunge( fEquation ); break;
269 case 8: pStepper = new G4CashKarpRKF45( fEquation ); break;
270 case 9: pStepper = new G4ExactHelixStepper( fEquation ); break;
271 case 10: pStepper = new G4RKG3_Stepper( fEquation ); break;
272 case 11: pStepper = new G4HelixMixedStepper( fEquation ); break;
273 case 12: pStepper = new G4ConstRK4( fEquation ); break;
274 case 13: pStepper = new G4NystromRK4( fEquation ); break;
275 default:
276 pStepper = 0;
277 G4cerr << " Stepper type provided is " << type << G4endl;
278 G4Exception(" Invalid value of stepper type");
279 break;
280 }
281
282 pFieldMgr= G4TransportationManager::GetTransportationManager()->
283 GetFieldManager();
284
285 pFieldMgr->SetDetectorField( &myMagField );
286
287 pChordFinder = new G4ChordFinder( &myMagField,
288 1.0e-2 * mm,
289 pStepper);
290 pChordFinder->SetVerbose(1); // ity();
291
292 pFieldMgr->SetChordFinder( pChordFinder );
293
294 return pFieldMgr;
295}
296
297G4PropagatorInField* SetupPropagator( G4int type)
298{
299 // G4FieldManager* fieldMgr=
300 SetupField( type) ;
301
302 // G4ChordFinder theChordFinder( &MagField, 0.05*mm ); // Default stepper
303
304 G4PropagatorInField *thePropagator =
305 G4TransportationManager::GetTransportationManager()->
306 GetPropagatorInField ();
307
308 // Let us test the new Minimum Epsilon Step functionality
309 // thePropagator -> SetMinimumEpsilonStep( 1.0e-3 ) ;
310 // thePropagator -> SetMaximumEpsilonStep( 1.0e-5 ) ;
311
312 return thePropagator;
313}
314
315// This is Done only for this test program ... the transportation does it.
316// The method is now obsolete -- as propagator in Field has this method,
317// in order to message the correct field manager's chord finder.
318//
319void ObsoleteSetChargeMomentumMass(G4double charge, G4double MomentumXc, G4double Mass)
320{
321 G4ChordFinder* pChordFinder;
322
323 pChordFinder= G4TransportationManager::GetTransportationManager()->
324 GetFieldManager()->GetChordFinder();
325
326 pChordFinder->SetChargeMomentumMass(
327 charge, // charge in e+ units
328 MomentumXc, // Momentum in Mev/c ?
329 Mass );
330}
331
332G4PropagatorInField *pMagFieldPropagator=0;
333//
334// Test Stepping
335//
336G4bool testG4PropagatorInField(G4VPhysicalVolume*, // *pTopNode,
337 G4int type)
338{
339 void report_endPV(G4ThreeVector Position,
340 G4ThreeVector UnitVelocity,
341 G4double step_len,
342 G4double physStep,
343 G4double safety,
344 G4ThreeVector EndPosition,
345 G4ThreeVector EndUnitVelocity,
346 G4int Step,
347 G4VPhysicalVolume* startVolume);
348
349 G4UniformMagField MagField(10.*tesla, 0., 0.);
350 G4Navigator *pNavig= G4TransportationManager::
351 GetTransportationManager()-> GetNavigatorForTracking();
352
353 pMagFieldPropagator= SetupPropagator(type);
354
355 pMagFieldPropagator->SetChargeMomentumMass( +1., // charge in e+ units
356 0.5 * proton_mass_c2, // Momentum in Mev/c
357 proton_mass_c2 );
358 // pNavig->SetWorldVolume(pTopNode);
359
360 G4VPhysicalVolume *located;
361 G4double step_len, physStep, safety;
362 G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
363 G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
364
365 // physStep=kInfinity;
366 G4ThreeVector Position(0.,0.,0.);
367 G4ThreeVector UnitMomentum(0.,0.6,0.8);
368 G4ThreeVector EndPosition, EndUnitMomentum;
369
370//
371// Test location & Step computation
372//
373 /* assert(located->GetName()=="World"); */
374 if( std::fabs(UnitMomentum.mag() - 1.0) > 1.e-8 )
375 {
376 G4cerr << "UnitMomentum.mag() - 1.0 = " << UnitMomentum.mag() - 1.0 <<
377 G4endl;
378 }
379
380 G4cout << G4endl;
381
382 for( int iparticle=0; iparticle < 2; iparticle++ )
383 {
384 physStep= 2.5 * mm ; // millimeters
385 Position = G4ThreeVector(0.,0.,0.)
386 + iparticle * G4ThreeVector(0.2, 0.3, 0.4);
387 UnitMomentum = (G4ThreeVector(0.,0.6,0.8)
388 + (float)iparticle * G4ThreeVector(0.1, 0.2, 0.3)).unit();
389
390 G4double momentum = (0.5+iparticle*10.0) * proton_mass_c2;
391
392 G4double kineticEnergy = momentum*momentum /
393 ( std::sqrt( momentum*momentum + proton_mass_c2 * proton_mass_c2 )
394 + proton_mass_c2 );
395 G4double velocity = momentum / ( proton_mass_c2 + kineticEnergy );
396 G4double labTof= 10.0*ns, properTof= 0.1*ns;
397 G4ThreeVector Spin(1.0, 0.0, 0.0);
398 // Momentum in Mev/c ?
399 pMagFieldPropagator->SetChargeMomentumMass(
400 +1, // charge in e+ units
401 momentum,
402 proton_mass_c2);
403 G4cout << G4endl;
404 G4cout << "Test PropagateMagField: ***********************" << G4endl
405 << " Starting New Particle with Position " << Position << G4endl
406 << " and UnitVelocity " << UnitMomentum << G4endl;
407 G4cout << " Momentum in GeV/c is " << momentum / GeV
408 << " = " << (0.5+iparticle*10.0)*proton_mass_c2 / MeV << " MeV"
409 << G4endl;
410
411
412 for( int istep=0; istep < 14; istep++ ){
413 // G4cerr << "UnitMomentum Magnitude is " << UnitMomentum.mag() << G4endl;
414 located= pNavig->LocateGlobalPointAndSetup(Position);
415 // G4cerr << "Starting Step " << istep << " in volume "
416 // << located->GetName() << G4endl;
417
418 G4FieldTrack initTrack( Position,
419 UnitMomentum,
420 0.0, // starting S curve len
421 kineticEnergy,
422 proton_mass_c2,
423 velocity,
424 labTof,
425 properTof,
426 0 // or &Spin
427 );
428
429 step_len=pMagFieldPropagator->ComputeStep( initTrack,
430 physStep,
431 safety,
432 located);
433 // --------------------
434 EndPosition= pMagFieldPropagator->EndPosition();
435 EndUnitMomentum= pMagFieldPropagator->EndMomentumDir();
436 // --------
437
438 if( std::fabs(EndUnitMomentum.mag2() - 1.0) > 1.e-8 )
439 G4cerr << "EndUnitMomentum.mag2() - 1.0 = " <<
440 EndUnitMomentum.mag2() - 1.0 << G4endl;
441
442 G4ThreeVector MoveVec = EndPosition - Position;
443 assert( MoveVec.mag() < physStep*(1.+1.e-9) );
444
445 // G4cout << " testPropagatorInField: After stepI " << istep << " : " << G4endl;
446 report_endPV(Position, UnitMomentum, step_len, physStep, safety,
447 EndPosition, EndUnitMomentum, istep, located );
448
449 assert(safety>=0);
450 pNavig->SetGeometricallyLimitedStep();
451 // pMagFieldPropagator->SetGeometricallyLimitedStep();
452
453 Position= EndPosition;
454 UnitMomentum= EndUnitMomentum;
455 physStep *= 2.;
456 } // ........................... end for ( istep )
457
458 myMagField.ReportStatistics();
459
460 } // .............................. end for ( iparticle )
461
462 return(1);
463}
464
465void report_endPV(G4ThreeVector Position,
466 G4ThreeVector InitialUnitVelocity,
467 G4double step_len,
468 G4double physStep,
469 G4double safety,
470 G4ThreeVector EndPosition,
471 G4ThreeVector EndUnitVelocity,
472 G4int Step,
473 G4VPhysicalVolume* startVolume)
474 // G4VPhysicalVolume* endVolume)
475{
476 const G4int verboseLevel=1;
477
478 if( Step == 0 && verboseLevel <= 3 )
479 {
480 G4cout.precision(6);
481 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
482 G4cout << std::setw( 5) << "Step#" << " "
483 << std::setw( 9) << "X(mm)" << " "
484 << std::setw( 9) << "Y(mm)" << " "
485 << std::setw( 9) << "Z(mm)" << " "
486 << std::setw( 9) << " N_x " << " "
487 << std::setw( 9) << " N_y " << " "
488 << std::setw( 9) << " N_z " << " "
489 << std::setw( 9) << " Delta|N|" << " "
490 << std::setw( 9) << " Delta(N_z) " << " "
491 // << std::setw( 9) << "KinE(MeV)" << " "
492 // << std::setw( 9) << "dE(MeV)" << " "
493 << std::setw( 9) << "StepLen" << " "
494 << std::setw( 9) << "PhsStep" << " "
495 << std::setw( 9) << "Safety" << " "
496 << std::setw(18) << "NextVolume" << " "
497 << G4endl;
498 }
499 //
500 //
501 if( verboseLevel > 3 )
502 {
503 G4cout << "End Position is " << EndPosition << G4endl
504 << " and UnitVelocity is " << EndUnitVelocity << G4endl;
505 G4cout << "Step taken was " << step_len
506 << " out of PhysicalStep= " << physStep << G4endl;
507 G4cout << "Final safety is: " << safety << G4endl;
508
509 G4cout << "Chord length = " << (EndPosition-Position).mag() << G4endl;
510 G4cout << G4endl;
511 }
512 else // if( verboseLevel > 0 )
513 {
514 G4cout.precision(6);
515 G4cout << std::setw( 5) << Step << " "
516 << std::setw( 9) << Position.x() << " "
517 << std::setw( 9) << Position.y() << " "
518 << std::setw( 9) << Position.z() << " "
519 << std::setw( 9) << EndUnitVelocity.x() << " "
520 << std::setw( 9) << EndUnitVelocity.y() << " "
521 << std::setw( 9) << EndUnitVelocity.z() << " ";
522 G4cout.precision(2);
523 G4cout
524 << std::setw( 9) << EndUnitVelocity.mag()-InitialUnitVelocity.mag() << " "
525 << std::setw( 9) << EndUnitVelocity.z() - InitialUnitVelocity.z() << " ";
526 // << std::setw( 9) << KineticEnergy << " "
527 // << std::setw( 9) << EnergyDifference << " "
528 G4cout.precision(6);
529 G4cout
530 << std::setw( 9) << step_len << " "
531 << std::setw( 9) << physStep << " "
532 << std::setw( 9) << safety << " ";
533 if( startVolume != 0) {
534 G4cout << std::setw(12) << startVolume->GetName() << " ";
535 } else {
536 G4cout << std::setw(12) << "OutOfWorld" << " ";
537 }
538#if 0
539 if( endVolume != 0)
540 G4cout << std::setw(12) << endVolume()->GetName() << " ";
541 else
542 G4cout << std::setw(12) << "OutOfWorld" << " ";
543#endif
544 G4cout << G4endl;
545 }
546}
547
548// Main program
549// -------------------------------
550int main(int argc, char **argv)
551{
552 G4VPhysicalVolume *myTopNode;
553 G4int type, optim, optimSaf;
554 G4bool optimiseVoxels=true;
555 G4bool optimisePiFwithSafety=true;
556
557 type = 8 ;
558 G4cout << " Arguments: stepper-no optimise-Voxels optimise-PiF-with-safety" << G4endl;
559
560 if( argc >= 2 ){
561 type = atoi(argv[1]);
562 }
563
564 if( argc >=3 ){
565 optim= atoi(argv[2]);
566 if( optim == 0 ) { optimiseVoxels = false; }
567 }
568
569 if( argc >=4 ){
570 optimSaf= atoi(argv[3]);
571 if( optimSaf == 0 ) { optimisePiFwithSafety= false; }
572 }
573
574 G4cout << " Testing with stepper number " << type << G4endl;
575 G4cout << " " ;
576 G4cout << " voxel optimisation " ;
577 // if (optimiseVoxels) G4cout << "On";
578 // else G4cout << "Off";
579 G4cout << (optimiseVoxels ? "On" : "Off") << G4endl;
580 G4cout << " " ;
581 G4cout << " Propagator safety optim " ;
582 // const char* OnOff= (optimisePiFwithSafety ? "on" : "off") ;
583 // G4cout << OnOff << G4endl;
584 G4cout << (optimisePiFwithSafety ? "On" : "Off") << G4endl;
585
586 // Create the geometry & field
587 myTopNode=BuildGeometry(); // Build the geometry
588
589 G4Navigator *pNavig= G4TransportationManager::
590 GetTransportationManager()-> GetNavigatorForTracking();
591 pNavig->SetWorldVolume(myTopNode);
592
593 G4GeometryManager::GetInstance()->CloseGeometry(false);
594
595 // Setup the propagator (will be overwritten by testG4Propagator ...)
596 pMagFieldPropagator= SetupPropagator(type);
597 G4cout << " Using default values for "
598 << " Min Eps = " << pMagFieldPropagator->GetMinimumEpsilonStep()
599 << " and "
600 << " MaxEps = " << pMagFieldPropagator->GetMaximumEpsilonStep()
601 << G4endl;
602
603 pMagFieldPropagator->SetUseSafetyForOptimization(optimisePiFwithSafety);
604
605// Do the tests without voxels
606 G4cout << " Test with no voxels" << G4endl;
607 testG4PropagatorInField(myTopNode, type);
608
609 pMagFieldPropagator->SetUseSafetyForOptimization(optimiseVoxels);
610 pMagFieldPropagator->SetVerboseLevel( 1 );
611
612// Repeat tests but with full voxels
613 G4cout << " Test with full voxels" << G4endl;
614
615 G4GeometryManager::GetInstance()->OpenGeometry();
616 G4GeometryManager::GetInstance()->CloseGeometry(true);
617
618 testG4PropagatorInField(myTopNode, type);
619
620 G4GeometryManager::GetInstance()->OpenGeometry();
621
622 G4cout << G4endl
623 << "----------------------------------------------------------"
624 << G4endl;
625
626// Repeat tests with full voxels and modified parameters
627 G4cout << "Test with more accurate parameters " << G4endl;
628
629 G4double maxEpsStep= 0.001;
630 G4double minEpsStep= 2.5e-8;
631 G4cout << " Setting values for Min Eps = " << minEpsStep
632 << " and MaxEps = " << maxEpsStep << G4endl;
633
634 pMagFieldPropagator->SetMaximumEpsilonStep(maxEpsStep);
635 pMagFieldPropagator->SetMinimumEpsilonStep(minEpsStep);
636
637 G4GeometryManager::GetInstance()->OpenGeometry();
638 G4GeometryManager::GetInstance()->CloseGeometry(true);
639
640 testG4PropagatorInField(myTopNode, type);
641
642 G4GeometryManager::GetInstance()->OpenGeometry();
643
644 optimiseVoxels = ! optimiseVoxels;
645// Repeat tests but with the opposite optimisation choice
646 G4cout << " Now test with optimisation " ;
647 if (optimiseVoxels) G4cout << "on";
648 else G4cout << "off";
649 G4cout << G4endl;
650
651 pMagFieldPropagator->SetUseSafetyForOptimization(optimiseVoxels);
652 testG4PropagatorInField(myTopNode, type);
653
654 G4GeometryManager::GetInstance()->OpenGeometry();
655
656
657
658
659 return 0;
660}
661
662
663
Note: See TracBrowser for help on using the repository browser.