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

Last change on this file since 1320 was 1231, checked in by garnier, 14 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.