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

Last change on this file since 1212 was 1199, checked in by garnier, 16 years ago

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