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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 23.3 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: testPropagateSpin.cc,v 1.17 2006/06/29 18:25:06 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
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.