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

Last change on this file since 1201 was 1199, checked in by garnier, 15 years ago

nvx fichiers dans CVS
File size: 22.0 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.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.