source: trunk/source/geometry/magneticfield/test/testPropagateSpin.cc@ 1323

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

update from CVS
File size: 23.3 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: testPropagateSpin.cc,v 1.17 2006/06/29 18:25:06 gunter Exp $
[1231]28// GEANT4 tag $Name: $
[1199]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 // World Volume
126 //
127 G4LogicalVolume *worldLog=new G4LogicalVolume(myHugeBox,0,
128 "World",0,0,0);
129 // Logical with no material,field,
130 // sensitive detector or user limits
131
132 G4PVPlacement *worldPhys=new
133 G4PVPlacement(0,G4ThreeVector(0,0,0), "World",worldLog,
134 0,false,0);
135 // Note: no mother pointer set
136
137// Create the logical Volumes
138//
139// G4LogicalVolume(*pSolid, *pMaterial, Name, *pField, *pSDetector, *pULimits)
140//
141 G4LogicalVolume *BigBoxLog=new G4LogicalVolume(myBigBox,0,
142 "Crystal Box (large)",0,0,0);
143 G4LogicalVolume *smallBoxLog=new G4LogicalVolume(mySmallBox,0,
144 "Crystal Box (small)");
145 G4LogicalVolume *tinyBoxLog=new G4LogicalVolume(myTinyBox,0,
146 "Crystal Box (tiny)");
147
148
149// Place them.
150//
151// 1) Two big boxes in the world volume
152//
153 // G4PVPlacement *BigTg1Phys=
154 new G4PVPlacement(0,G4ThreeVector(0,0,-15*m),
155 "Big Target 1",BigBoxLog,
156 worldPhys,false,0);
157 // G4PVPlacement *BigTg2Phys=
158 new G4PVPlacement(0,G4ThreeVector(0,0, 15*m),
159 "Big Target 2",BigBoxLog,
160 worldPhys,false,0);
161
162// 2) Four (medium) boxes in X & Y near the origin of the world volume
163//
164 // G4PVPlacement *MedTg3a_Phys=
165 new G4PVPlacement(0,G4ThreeVector(0, 7.5*m,0),
166 "Target 3a",smallBoxLog,
167 worldPhys,false,0);
168 // G4PVPlacement *MedTg3b_Phys=
169 new G4PVPlacement(0,G4ThreeVector(0,-7.5*m,0),
170 "Target 3b",smallBoxLog,
171 worldPhys,false,0);
172 // G4PVPlacement *MedTg3c_Phys=
173 new G4PVPlacement(0,G4ThreeVector(-7.5*m,0,0),
174 "Target 3c",smallBoxLog,
175 worldPhys,false,0);
176 // G4PVPlacement *MedTg3d_Phys=
177 new G4PVPlacement(0,G4ThreeVector( 7.5*m,0,0),
178 "Target 3d",smallBoxLog,
179 worldPhys,false,0);
180
181
182// 3) Eight small boxes around the origin of the world volume
183// (in +-X, +-Y & +-Z)
184//
185 // G4PVPlacement *SmTg4a_Phys=
186 new G4PVPlacement
187 (0,G4ThreeVector( 0.3*m, 0.3*m,0.3*m), "Target 4a",tinyBoxLog,
188 worldPhys,false,0);
189 // G4PVPlacement *SmTg4b_Phys=
190 new G4PVPlacement
191 (0,G4ThreeVector( 0.3*m,-0.3*m,0.3*m), "Target 4b",tinyBoxLog,
192 worldPhys,false,0);
193 // G4PVPlacement *SmTg4c_Phys=
194 new G4PVPlacement
195 (0,G4ThreeVector(-0.3*m,-0.3*m,0.3*m), "Target 4c",tinyBoxLog,
196 worldPhys,false,0);
197 // G4PVPlacement *SmTg4d_Phys=
198 new G4PVPlacement
199 (0,G4ThreeVector(-0.3*m, 0.3*m,0.3*m), "Target 4d",tinyBoxLog,
200 worldPhys,false,0);
201
202 // G4PVPlacement *SmTg4e_Phys=
203 new G4PVPlacement
204 (0,G4ThreeVector( 0.3*m, 0.3*m,-0.3*m), "Target 4e",tinyBoxLog,
205 worldPhys,false,0);
206 // G4PVPlacement *SmTg4f_Phys=
207 new G4PVPlacement
208 (0,G4ThreeVector( 0.3*m,-0.3*m,-0.3*m), "Target 4f",tinyBoxLog,
209 worldPhys,false,0);
210 // G4PVPlacement *SmTg4g_Phys=
211 new G4PVPlacement
212 (0,G4ThreeVector(-0.3*m,-0.3*m,-0.3*m), "Target 4g",tinyBoxLog,
213 worldPhys,false,0);
214 // G4PVPlacement *SmTg4h_Phys=
215 new G4PVPlacement
216 (0,G4ThreeVector(-0.3*m, 0.3*m,-0.3*m), "Target 4h",tinyBoxLog,
217 worldPhys,false,0);
218
219 return worldPhys;
220}
221
222#include "G4ChordFinder.hh"
223#include "G4PropagatorInField.hh"
224#include "G4MagneticField.hh"
225#include "G4FieldManager.hh"
226#include "G4TransportationManager.hh"
227#include "G4HelixExplicitEuler.hh"
228#include "G4HelixSimpleRunge.hh"
229#include "G4HelixImplicitEuler.hh"
230#include "G4ExplicitEuler.hh"
231#include "G4ImplicitEuler.hh"
232#include "G4SimpleRunge.hh"
233#include "G4SimpleHeum.hh"
234#include "G4ClassicalRK4.hh"
235#include "G4Mag_SpinEqRhs.hh"
236#include "G4CashKarpRKF45.hh"
237#include "G4RKG3_Stepper.hh"
238
239G4UniformMagField myMagField(10.*tesla, 0., 0.);
240
241
242G4FieldManager* SetupField(G4int type)
243{
244 G4FieldManager *pFieldMgr;
245 G4ChordFinder *pChordFinder;
246 G4Mag_SpinEqRhs *fEquation = new G4Mag_SpinEqRhs(&myMagField);
247 G4MagIntegratorStepper *pStepper;
248
249 const int ncompspin=12;
250
251 switch ( type )
252 {
253 case 0: pStepper = new G4ExplicitEuler( fEquation, ncompspin ); break;
254 case 1: pStepper = new G4ImplicitEuler( fEquation, ncompspin ); break;
255 case 2: pStepper = new G4SimpleRunge( fEquation, ncompspin ); break;
256 case 3: pStepper = new G4SimpleHeum( fEquation, ncompspin ); break;
257 case 4: pStepper = new G4ClassicalRK4( fEquation, ncompspin ); break;
258 case 8: pStepper = new G4CashKarpRKF45( fEquation, ncompspin ); break;
259 default: pStepper = new G4ClassicalRK4( fEquation, ncompspin ); break;
260 }
261
262 pFieldMgr= G4TransportationManager::GetTransportationManager()->
263 GetFieldManager();
264
265 pFieldMgr->SetDetectorField( &myMagField );
266
267 pChordFinder = new G4ChordFinder( &myMagField,
268 1.0e-2 * mm,
269 pStepper);
270 pChordFinder->SetVerbose(1); // ity();
271
272 pFieldMgr->SetChordFinder( pChordFinder );
273
274 return pFieldMgr;
275}
276
277G4PropagatorInField* SetupPropagator( G4int type)
278{
279 // G4FieldManager* fieldMgr=
280 SetupField( type) ;
281
282 // G4ChordFinder theChordFinder( &MagField, 0.05*mm ); // Default stepper
283
284 G4PropagatorInField *thePropagator =
285 G4TransportationManager::GetTransportationManager()->
286 GetPropagatorInField ();
287
288 // Let us test the new Minimum Epsilon Step functionality
289 thePropagator -> SetMinimumEpsilonStep( 1.0e-7 ) ;
290 thePropagator -> SetMaximumEpsilonStep( 1.0e-7 ) ;
291
292 return thePropagator;
293}
294
295// This is Done only for this test program ... the transportation does it.
296// The method is now obsolete -- as propagator in Field has this method,
297// in order to message the correct field manager's chord finder.
298//
299void ObsoleteSetChargeMomentumMass(G4double charge, G4double MomentumXc, G4double Mass)
300{
301 G4ChordFinder* pChordFinder;
302
303 pChordFinder= G4TransportationManager::GetTransportationManager()->
304 GetFieldManager()->GetChordFinder();
305
306 pChordFinder->SetChargeMomentumMass(
307 charge, // charge in e+ units
308 MomentumXc, // Momentum in Mev/c ?
309 Mass );
310}
311
312G4PropagatorInField *pMagFieldPropagator;
313//
314// Test Stepping
315//
316G4bool testG4PropagatorInField(G4VPhysicalVolume*, // *pTopNode,
317 G4int ) // type)
318{
319 void report_endPV(G4ThreeVector Position,
320 G4ThreeVector UnitVelocity,
321 G4ThreeVector Spin,
322 G4double step_len,
323 G4double physStep,
324 G4double safety,
325 G4ThreeVector EndPosition,
326 G4ThreeVector EndUnitVelocity,
327 G4ThreeVector EndSpin,
328 G4int Step,
329 G4VPhysicalVolume* startVolume);
330
331 G4UniformMagField MagField(10.*tesla, 0., 0.); // Tesla Defined ?
332 G4TransportationManager* transpMgr = G4TransportationManager::
333 GetTransportationManager();
334 G4Navigator* pNavig= transpMgr-> GetNavigatorForTracking();
335
336 // pMagFieldPropagator= SetupPropagator(type);
337
338 G4cout << "Test G4PropInFld with "
339 << "optimise = "
340 << ( pMagFieldPropagator->GetUseSafetyForOptimization() ? "on" : "off" )
341 // << " Eps min= " << pMagFieldPropagator->GetMinimumEpsilonStep()
342 // << " & max= " << pMagFieldPropagator->GetMaximumEpsilonStep()
343 << G4endl;
344
345 const G4FieldManager* pFieldMgr= transpMgr->GetFieldManager();
346 G4cout << " The global field manager has the following parameters "
347 << G4endl;
348 G4cout << " Eps min= " << pFieldMgr->GetMinimumEpsilonStep()
349 << " & max= " << pFieldMgr->GetMaximumEpsilonStep()
350 << G4endl;
351 G4cout << " Delta Intersection= " << pFieldMgr->GetDeltaIntersection()
352 << " Delta One step = " << pFieldMgr->GetDeltaOneStep()
353 << G4endl;
354
355
356 pMagFieldPropagator->SetChargeMomentumMass(
357 +1., // charge in e+ units
358 0.1*GeV, // Momentum in Mev/c ?
359 0.105658387*GeV );
360 // pNavig->SetWorldVolume(pTopNode);
361
362
363 G4VPhysicalVolume *located;
364 G4double step_len, physStep, safety;
365 G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
366 G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
367
368 // physStep=kInfinity;
369 G4ThreeVector Position(0.,0.,0.);
370 G4ThreeVector UnitMomentum(0.,0.6,0.8);
371 G4ThreeVector EndPosition, EndUnitMomentum;
372
373//
374// Test location & Step computation
375//
376 /* assert(located->GetName()=="World"); */
377
378 const G4double threshold= 1.e-6;
379
380 if( std::fabs(UnitMomentum.mag() - 1.0) > threshold )
381 {
382 G4cout << "UnitMomentum.mag() - 1.0 = " << UnitMomentum.mag() - 1.0 <<
383 G4endl;
384 }
385
386 G4cout << G4endl;
387
388 for( int iparticle=0; iparticle < 2; iparticle++ )
389 {
390 physStep= 2.5 * mm ; // millimeters
391 Position = G4ThreeVector(0.,0.,0.)
392 + iparticle * G4ThreeVector(0.2, 0.3, 0.4);
393 // ->GetChordFinder().SetChargeAndMomentum(
394
395 G4double momentum_val= (0.5+iparticle*1.0) * 0.1*GeV; // As energy/c
396 G4double rest_mass = 0.105658387*GeV; // A muon
397
398 G4double momentum_sq = momentum_val * momentum_val;
399 G4double kineticEnergy = momentum_sq /
400 ( std::sqrt( momentum_sq + rest_mass * rest_mass )
401 + rest_mass );
402 G4double labTof= 10.0*ns, properTof= 0.1*ns;
403 pMagFieldPropagator->SetChargeMomentumMass(
404 +1, // charge in e+ units
405 momentum_val,
406 rest_mass);
407
408 UnitMomentum = (G4ThreeVector(0.,0.6,0.8)
409 + (float)iparticle * G4ThreeVector(0.1, 0.2, 0.3)).unit();
410 G4double beta = momentum_val / std::sqrt( rest_mass*rest_mass + momentum_val*momentum_val );
411 G4double VelocityMag = beta * c_light;
412 G4ThreeVector Velocity = VelocityMag * UnitMomentum ;
413
414 G4cout << G4endl;
415 G4cout << "Test PropagateMagField: ***********************" << G4endl
416 << " Starting New Particle with Position " << Position << G4endl
417 << " and UnitVelocity " << UnitMomentum << G4endl;
418 G4cout << " Momentum in MeV/c is "<< (0.5+iparticle*1.0)*0.1*GeV/MeV;
419 G4cout << G4endl;
420
421 G4ThreeVector initialSpin = UnitMomentum;
422
423 for( int istep=0; istep < 14; istep++ ){
424 // // G4cout << "UnitMomentum Magnitude is " << UnitMomentum.mag() << G4endl;
425 located= pNavig->LocateGlobalPointAndSetup(Position);
426 // Is the following better ?? It would need "changes"
427 // located= pMagFieldPropagator->LocateGlobalPointAndSetup(Position);
428 // G4cout << "Starting Step " << istep << " in volume "
429 // << located->GetName() << G4endl;
430
431 // G4FieldTrack stateVec( Position, Velocity, 0.0, 0.0,
432 // 0.0, 0.0, &initialSpin );
433
434 G4FieldTrack stateVec( Position,
435 UnitMomentum,
436 0.0, // starting S curve len
437 kineticEnergy,
438 rest_mass,
439 VelocityMag,
440 labTof,
441 properTof,
442 &initialSpin
443 );
444
445 step_len=pMagFieldPropagator->ComputeStep( stateVec,
446 physStep, safety
447 ,located);
448 // --------------------
449 EndPosition= pMagFieldPropagator->EndPosition();
450 EndUnitMomentum= pMagFieldPropagator->EndMomentumDir();
451 // --------
452 G4FieldTrack EndFieldTrack= pMagFieldPropagator->GetEndState();
453 G4ThreeVector EndSpin= EndFieldTrack.GetSpin();
454 G4ThreeVector EndUnitMomentum = EndFieldTrack.GetMomentumDir();
455
456// G4cout << " EndPosition " << EndPosition << G4endl;
457// G4cout << " EndUnitMomentum " << EndUnitMomentum << G4endl;
458// G4cout << " initialSpin " << initialSpin.mag() << G4endl;
459// G4cout << " EndSpin " << EndSpin.mag() << G4endl;
460
461 if( std::fabs(EndUnitMomentum.mag2() - 1.0) > threshold )
462 G4cout << "EndUnitMomentum.mag2() - 1.0 = " <<
463 EndUnitMomentum.mag2() - 1.0 << G4endl;
464
465 // In this case spin should be parallel (equal) to momentum
466 G4double magdiff= (EndUnitMomentum - EndSpin).mag();
467 if( magdiff > 1.e-8 ){
468 G4cout.precision(4);
469 G4cout << " Spin is not equal to Momentum "
470 << " Diff = " << magdiff << G4endl;
471 }
472
473 G4ThreeVector MoveVec = EndPosition - Position;
474 assert( MoveVec.mag() < physStep*(1.+1.e-9) );
475
476 // G4cout << " testPropagatorInField: After stepI " << istep << " : " << G4endl;
477 report_endPV(Position, UnitMomentum, initialSpin, step_len, physStep, safety,
478 EndPosition, EndUnitMomentum, EndSpin, istep, located );
479
480 assert(safety>=0);
481 pNavig->SetGeometricallyLimitedStep();
482 // pMagFieldPropagator->SetGeometricallyLimitedStep();
483
484 Position= EndPosition;
485 UnitMomentum= EndUnitMomentum;
486 initialSpin = EndSpin;
487
488 physStep *= 2.;
489 } // ........................... end for ( istep )
490 } // .............................. end for ( iparticle )
491
492 return(1);
493}
494
495void report_endPV(G4ThreeVector Position,
496 G4ThreeVector, // UnitVelocity,
497 G4ThreeVector, // Spin,
498 G4double step_len,
499 G4double physStep,
500 G4double safety,
501 G4ThreeVector EndPosition,
502 G4ThreeVector EndUnitVelocity,
503 G4ThreeVector EndSpin,
504 G4int Step,
505 G4VPhysicalVolume* startVolume)
506 // G4VPhysicalVolume* endVolume)
507{
508 const G4int verboseLevel=1;
509 G4int oldPrec= G4cout.precision(4);
510
511 if( Step == 0 && verboseLevel <= 3 )
512 {
513 // G4cout.precision(6);
514 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
515 G4cout << std::setw( 3) << "Stp#" << " "
516 << std::setw( 7) << "X(mm)" << " "
517 << std::setw( 7) << "Y(mm)" << " "
518 << std::setw( 7) << "Z(mm)" << " "
519 << std::setw( 7) << " N_x " << " "
520 << std::setw( 7) << " N_y " << " "
521 << std::setw( 7) << " N_z " << " "
522 << std::setw( 7) << " S_x " << " "
523 << std::setw( 7) << " S_y " << " "
524 << std::setw( 7) << " S_z " << " "
525 << std::setw( 9) << " |S-N|" << " "
526 << std::setw( 9) << " (S_z-N_z) " << " "
527 // << std::setw( 9) << "KinE(MeV)" << " "
528 // << std::setw( 9) << "dE(MeV)" << " "
529 << std::setw( 9) << "StepLen" << " "
530 << std::setw( 9) << "PhsStep" << " "
531 << std::setw( 9) << "Safety" << " "
532 << std::setw(18) << "NextVolume" << " "
533 << G4endl;
534 }
535 //
536 //
537 if( verboseLevel > 3 )
538 {
539 G4cout << "End Position is " << EndPosition << G4endl
540 << " and UnitVelocity is " << EndUnitVelocity << G4endl;
541 G4cout << "Step taken was " << step_len
542 << " out of PhysicalStep= " << physStep << G4endl;
543 G4cout << "Final safety is: " << safety << G4endl;
544
545 G4cout << "Chord length = " << (EndPosition-Position).mag() << G4endl;
546 G4cout << G4endl;
547 }
548 else // if( verboseLevel > 0 )
549 {
550 G4cout.precision(3); // 4 ?
551 G4cout << std::setw( 3) << Step << " "
552 << std::setw( 7) << Position.x() << " "
553 << std::setw( 7) << Position.y() << " "
554 << std::setw( 7) << Position.z() << " "
555 << std::setw( 7) << EndUnitVelocity.x() << " "
556 << std::setw( 7) << EndUnitVelocity.y() << " "
557 << std::setw( 7) << EndUnitVelocity.z() << " "
558 << std::setw( 7) << EndSpin.x() << " "
559 << std::setw( 7) << EndSpin.y() << " "
560 << std::setw( 7) << EndSpin.z() << " ";
561 G4cout.precision(2);
562 G4cout << std::setw( 8) << (EndSpin-EndUnitVelocity).mag() << " "
563 << std::setw( 8) << EndSpin.z() - EndUnitVelocity.z() << " ";
564 // << std::setw( 9) << KineticEnergy << " "
565 // << std::setw( 9) << EnergyDifference << " "
566 G4cout.precision(6);
567 G4cout << std::setw( 9) << step_len << " "
568 << std::setw( 9) << physStep << " ";
569 G4cout.precision(3); // could be 4 ?
570 G4cout << std::setw( 9) << safety << " ";
571 if( startVolume != 0) {
572 G4cout << std::setw(12) << startVolume->GetName() << " ";
573 } else {
574 G4cout << std::setw(12) << "OutOfWorld" << " ";
575 }
576
577#if 0
578 if( endVolume != 0) {
579 G4cout << std::setw(12) << endVolume()->GetName() << " ";
580 } else {
581 G4cout << std::setw(12) << "OutOfWorld" << " ";
582 }
583#endif
584 G4cout << G4endl;
585 }
586
587 G4cout.precision(oldPrec);
588}
589
590// Main program
591// -------------------------------
592int main(int argc, char **argv)
593{
594 G4VPhysicalVolume *myTopNode;
595 G4int type, optim;
596 G4bool optimise=true;
597
598 type = 8 ;
599
600 if( argc >= 2 )
601 type = atoi(argv[1]);
602
603 if( argc >=3 ){
604 optim= atoi(argv[2]);
605 if( optim == 0 ) { optimise = false; }
606 }
607
608 G4cout << " Testing with stepper number " << type;
609 G4cout << " and PiF safety optimisation " ;
610 if (optimise) G4cout << "on";
611 else G4cout << "off";
612 G4cout << G4endl;
613
614 // Create the geometry & field
615 myTopNode=BuildGeometry(); // Build the geometry
616
617 G4Navigator *pNavig= G4TransportationManager::
618 GetTransportationManager()-> GetNavigatorForTracking();
619 pNavig->SetWorldVolume(myTopNode);
620
621 G4GeometryManager::GetInstance()->CloseGeometry(false);
622
623 // Setup the propagator (will be overwritten by testG4Propagator ...)
624 pMagFieldPropagator= SetupPropagator(type);
625 G4cout << " Using default values for "
626 << " Min Eps = " << pMagFieldPropagator->GetMinimumEpsilonStep()
627 << " and "
628 << " MaxEps = " << pMagFieldPropagator->GetMaximumEpsilonStep()
629 << G4endl;
630
631 pMagFieldPropagator->SetUseSafetyForOptimization(optimise);
632
633// Do the tests without voxels
634 G4cout << " Test with no voxels" << G4endl;
635 testG4PropagatorInField(myTopNode, type);
636
637// Repeat tests but with full voxels
638 G4cout << " Test with full voxels" << G4endl;
639
640 G4GeometryManager::GetInstance()->OpenGeometry();
641 G4GeometryManager::GetInstance()->CloseGeometry(true);
642
643 testG4PropagatorInField(myTopNode, type);
644
645 G4GeometryManager::GetInstance()->OpenGeometry();
646
647 G4cout << G4endl
648 << "----------------------------------------------------------"
649 << G4endl;
650
651// Repeat tests with full voxels and modified parameters
652 G4cout << "Test with more accurate parameters " << G4endl;
653
654 G4double maxEpsStep= 0.001;
655 G4double minEpsStep= 2.5e-8;
656 G4cout << " Setting values for Min Eps = " << minEpsStep
657 << " and MaxEps = " << maxEpsStep << G4endl;
658
659 pMagFieldPropagator->SetMaximumEpsilonStep(maxEpsStep);
660 pMagFieldPropagator->SetMinimumEpsilonStep(minEpsStep);
661
662 G4GeometryManager::GetInstance()->OpenGeometry();
663 G4GeometryManager::GetInstance()->CloseGeometry(true);
664
665 testG4PropagatorInField(myTopNode, type);
666
667 G4GeometryManager::GetInstance()->OpenGeometry();
668
669 optimise = ! optimise;
670// Repeat tests but with the opposite optimisation choice
671 G4cout << " Now test with safety optimisation " ;
672 if (optimise) G4cout << "on";
673 else G4cout << "off";
674 G4cout << G4endl;
675
676 pMagFieldPropagator->SetUseSafetyForOptimization(optimise);
677 testG4PropagatorInField(myTopNode, type);
678
679 G4GeometryManager::GetInstance()->OpenGeometry();
680
681 return 0;
682}
683
Note: See TracBrowser for help on using the repository browser.