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

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

update from CVS
File size: 22.8 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: testProPerpSpin.cc,v 1.17 2006/06/29 18:25:02 gunter 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 "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
271 pFieldMgr->SetChordFinder( pChordFinder );
272
273 return pFieldMgr;
274}
275
276G4PropagatorInField* SetupPropagator( G4int type)
277{
278 // G4FieldManager* fieldMgr=
279 SetupField( type) ;
280
281 // G4ChordFinder theChordFinder( &MagField, 0.05*mm ); // Default stepper
282
283 G4PropagatorInField *thePropagator =
284 G4TransportationManager::GetTransportationManager()->
285 GetPropagatorInField ();
286
287 // Let us test the new Minimum Epsilon Step functionality
288 thePropagator -> SetMinimumEpsilonStep( 1.0e-8 ) ;
289 thePropagator -> SetMaximumEpsilonStep( 1.0e-8 ) ;
290
291 return thePropagator;
292}
293
294// This is Done only for this test program ... the transportation does it.
295// The method is now obsolete -- as propagator in Field has this method,
296// in order to message the correct field manager's chord finder.
297//
298void ObsoleteSetChargeMomentumMass(G4double charge, G4double MomentumXc, G4double Mass)
299{
300 G4ChordFinder* pChordFinder;
301
302 pChordFinder= G4TransportationManager::GetTransportationManager()->
303 GetFieldManager()->GetChordFinder();
304
305 pChordFinder->SetChargeMomentumMass(
306 charge, // charge in e+ units
307 MomentumXc, // Momentum in Mev/c ?
308 Mass );
309}
310
311//
312// Test Stepping
313//
314G4bool testG4PropagatorInField(G4VPhysicalVolume *pTopNode, G4int ) // type
315{
316 void report_endPV(G4ThreeVector Position,
317 G4ThreeVector UnitVelocity,
318 G4ThreeVector Spin,
319 G4double step_len,
320 G4double physStep,
321 G4double safety,
322 G4ThreeVector EndPosition,
323 G4ThreeVector EndUnitVelocity,
324 G4ThreeVector EndSpin,
325 G4int Step,
326 G4VPhysicalVolume* startVolume);
327
328 G4UniformMagField MagField(10.*tesla, 0., 0.); // Tesla Defined ?
329 G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
330 G4Navigator* pNavig= transportMgr-> GetNavigatorForTracking();
331 // G4Navigator *pNavig= G4TransportationManager::
332 // GetTransportationManager()-> GetNavigatorForTracking();
333
334 // G4PropagatorInField *pMagFieldPropagator= SetupPropagator(type);
335 G4PropagatorInField *pMagFieldPropagator=
336 transportMgr->GetPropagatorInField();
337
338 pMagFieldPropagator->SetChargeMomentumMass(
339 +1., // charge in e+ units
340 0.1*GeV, // Momentum in Mev/c ?
341 0.105658387*GeV );
342 pNavig->SetWorldVolume(pTopNode);
343
344 G4VPhysicalVolume *located;
345 G4double step_len, physStep, safety;
346 G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
347 G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
348
349 // physStep=kInfinity;
350 G4ThreeVector Position(0.,0.,0.);
351 G4ThreeVector UnitMomentum(0.,0.6,0.8);
352 G4ThreeVector EndPosition, EndUnitMomentum;
353
354//
355// Test location & Step computation
356//
357 /* assert(located->GetName()=="World"); */
358
359 const G4double threshold= 1.e-6;
360
361 if( std::fabs(UnitMomentum.mag() - 1.0) > threshold )
362 {
363 G4cout << "UnitMomentum.mag() - 1.0 = " << UnitMomentum.mag() - 1.0 <<
364 G4endl;
365 }
366
367 G4cout << G4endl;
368
369 for( int iparticle=0; iparticle < 2; iparticle++ )
370 {
371 physStep= 2.5 * mm ; // millimeters
372 Position = G4ThreeVector(0.,0.,0.)
373 + iparticle * G4ThreeVector(0.2, 0.3, 0.4);
374 UnitMomentum = (G4ThreeVector(0.,0.6,0.8)
375 + (float)iparticle * G4ThreeVector(0.1, 0.2, 0.3)).unit();
376
377 G4double momentum_val= (0.5+iparticle*1.0) * 0.1*GeV; // As energy/c
378 G4double rest_mass = 0.105658387*GeV; // A muon
379
380 G4double kineticEnergy = momentum_val*momentum_val /
381 ( std::sqrt( momentum_val*momentum_val + rest_mass * rest_mass )
382 + rest_mass );
383 G4double labTof= 10.0*ns, properTof= 0.1*ns;
384 pMagFieldPropagator->SetChargeMomentumMass(
385 +1, // charge in e+ units
386 momentum_val,
387 rest_mass);
388
389 G4double beta = momentum_val / std::sqrt( rest_mass*rest_mass + momentum_val*momentum_val );
390 G4double velocity_magnitude = beta * c_light;
391 // G4ThreeVector Velocity = UnitMomentum * velocity_magnitude;
392
393 G4cout << G4endl;
394 G4cout << "Test PropagateMagField: ***********************" << G4endl
395 << " Starting New Particle with Position " << Position << G4endl
396 << " and UnitVelocity " << UnitMomentum << G4endl;
397 G4cout << " Momentum in MeV/c is "<< (0.5+iparticle*1.0)*0.1*GeV/MeV;
398 G4cout << G4endl;
399
400 // G4ThreeVector initialSpin = UnitMomentum;
401 // G4ThreeVector initialSpin(UnitMomentum.y(),-UnitMomentum.x(),UnitMomentum.z());
402 G4ThreeVector initialSpin = (UnitMomentum.cross(G4ThreeVector(1.0,0.,0.))).unit();
403
404 G4cout << " Initial spin is set to = " << initialSpin << G4endl;
405 G4cout << " Initial dot product of spin and momentum = "
406 << initialSpin.dot(UnitMomentum) << G4endl;
407
408 for( int istep=0; istep < 14; istep++ ){
409 // // G4cout << "UnitMomentum Magnitude is " << UnitMomentum.mag() << G4endl;
410 located= pNavig->LocateGlobalPointAndSetup(Position);
411 // Is the following better ?? It would need "changes"
412 // located= pMagFieldPropagator->LocateGlobalPointAndSetup(Position);
413 // G4cout << "Starting Step " << istep << " in volume "
414 // << located->GetName() << G4endl;
415
416 G4ThreeVector spinValue( initialSpin );
417
418 G4FieldTrack initTrack( Position,
419 UnitMomentum,
420 0.0, // starting S curve len
421 kineticEnergy,
422 rest_mass,
423 velocity_magnitude,
424 labTof,
425 properTof,
426 &spinValue
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 G4FieldTrack EndFieldTrack= pMagFieldPropagator->GetEndState();
438 G4ThreeVector EndSpin= EndFieldTrack.GetSpin();
439 G4ThreeVector EndUnitMomentum = EndFieldTrack.GetMomentumDir();
440
441// G4cout << " EndPosition " << EndPosition << G4endl;
442// G4cout << " EndUnitMomentum " << EndUnitMomentum << G4endl;
443// G4cout << " initialSpin " << initialSpin.mag() << G4endl;
444// G4cout << " EndSpin " << EndSpin.mag() << G4endl;
445
446 if( std::fabs(EndUnitMomentum.mag2() - 1.0) > threshold )
447 G4cout << "EndUnitMomentum.mag2() - 1.0 = " <<
448 EndUnitMomentum.mag2() - 1.0 << G4endl;
449
450 // In this case spin should be parallel (equal) to momentum
451 G4double endDot= EndSpin.dot(EndUnitMomentum) ;
452 G4cout << " dot product of spin and momentum = " << endDot << G4endl;
453 if( std::fabs(endDot) > threshold ){
454 G4cout << " $$$$$$$$$$$$ Spin dot Momentum is above threshold of "
455 << threshold << G4endl ;
456 G4cout << " Spin dot UnitMomentum= " << endDot << " ";
457 G4cout << " Spin magnitude= " << EndSpin.mag() << " "
458 << "(-1=" << (EndSpin.mag()-1.0) << ") ";
459 G4cout << " UnitMom mag= " << EndUnitMomentum.mag() << " "
460 << "(-1=" << (EndUnitMomentum.mag()-1.0) << ") ";
461 G4cout << G4endl;
462 }
463 G4ThreeVector MoveVec = EndPosition - Position;
464 assert( MoveVec.mag() < physStep*(1.+1.e-9) );
465
466 // G4cout << " testPropagatorInField: After stepI " << istep << " : " << G4endl;
467 report_endPV(Position, UnitMomentum, initialSpin, step_len, physStep, safety,
468 EndPosition, EndUnitMomentum, EndSpin, istep, located );
469
470 assert(safety>=0);
471 pNavig->SetGeometricallyLimitedStep();
472 // pMagFieldPropagator->SetGeometricallyLimitedStep();
473
474 Position= EndPosition;
475 // Velocity= EndVelocity;
476 UnitMomentum= EndUnitMomentum;
477 initialSpin = EndSpin;
478
479 physStep *= 2.;
480 } // ........................... end for ( istep )
481 } // .............................. end for ( iparticle )
482
483 return(1);
484}
485
486// int main(int argc, char** argv)
487int main(int argc, char **argv)
488{
489 G4VPhysicalVolume *myTopNode;
490 G4int type;
491 myTopNode=BuildGeometry(); // Build the geometry
492 G4GeometryManager::GetInstance()->CloseGeometry(false);
493
494 type = 8 ;
495
496 if( argc == 2 )
497 type = atoi(argv[1]);
498
499 // Setup the Propagator and register it with the Transportation Manager
500 G4PropagatorInField *pMagFieldPropagator= SetupPropagator(type);
501
502 G4cout << " Using the following values for "
503 << " Min Eps = " << pMagFieldPropagator->GetMinimumEpsilonStep()
504 << " and "
505 << " Max Eps = " << pMagFieldPropagator->GetMaximumEpsilonStep()
506 << G4endl;
507
508// Do the tests without voxels
509 G4cout << " Test with no voxels" << G4endl;
510 testG4PropagatorInField(myTopNode, type);
511
512// Repeat tests but with full voxels
513 G4cout << " Test with full voxels" << G4endl;
514 G4GeometryManager::GetInstance()->OpenGeometry();
515 G4GeometryManager::GetInstance()->CloseGeometry(true);
516
517 testG4PropagatorInField(myTopNode, type);
518
519 G4GeometryManager::GetInstance()->OpenGeometry();
520 return 0;
521}
522
523
524void report_endPV(G4ThreeVector Position,
525 G4ThreeVector , // UnitVelocity,
526 G4ThreeVector , // Spin,
527 G4double step_len,
528 G4double physStep,
529 G4double safety,
530 G4ThreeVector EndPosition,
531 G4ThreeVector EndUnitVelocity,
532 G4ThreeVector EndSpin,
533 G4int Step,
534 G4VPhysicalVolume* startVolume)
535 // G4VPhysicalVolume* endVolume)
536{
537 const G4int verboseLevel=1;
538
539 if( Step == 0 && verboseLevel <= 3 )
540 {
541 G4cout.precision(3);
542 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
543 G4cout << std::setw( 5) << "Step#" << " "
544 << std::setw( 9) << "X(mm)" << " "
545 << std::setw( 9) << "Y(mm)" << " "
546 << std::setw( 9) << "Z(mm)" << " "
547 << std::setw( 7) << " N_x " << " "
548 << std::setw( 7) << " N_y " << " "
549 << std::setw( 7) << " N_z " << " "
550 << std::setw( 7) << " S_x " << " "
551 << std::setw( 7) << " S_y " << " "
552 << std::setw( 7) << " S_z " << " "
553 // << std::setw( 9) << "KinE(MeV)" << " "
554 // << std::setw( 9) << "dE(MeV)" << " "
555 << std::setw( 9) << "StepLen" << " "
556 << std::setw( 9) << "PhsStep" << " "
557 << std::setw( 9) << "Safety" << " "
558 << std::setw(18) << "NextVolume" << " "
559 << G4endl;
560 }
561 //
562 //
563 if( verboseLevel > 3 )
564 {
565 G4cout << "End Position is " << EndPosition << G4endl
566 << " and UnitVelocity is " << EndUnitVelocity << G4endl;
567 G4cout << "Step taken was " << step_len
568 << " out of PhysicalStep= " << physStep << G4endl;
569 G4cout << "Final safety is: " << safety << G4endl;
570
571 G4cout << "Chord length = " << (EndPosition-Position).mag() << G4endl;
572 G4cout << G4endl;
573 }
574 else // if( verboseLevel > 0 )
575 {
576 G4cout.precision(3);
577 G4cout << std::setw( 5) << Step << " "
578 << std::setw( 9) << Position.x() << " "
579 << std::setw( 9) << Position.y() << " "
580 << std::setw( 9) << Position.z() << " "
581 << std::setw( 7) << EndUnitVelocity.x() << " "
582 << std::setw( 7) << EndUnitVelocity.y() << " "
583 << std::setw( 7) << EndUnitVelocity.z() << " "
584 << std::setw( 7) << EndSpin.x() << " "
585 << std::setw( 7) << EndSpin.y() << " "
586 << std::setw( 7) << EndSpin.z() << " "
587 // << std::setw( 9) << KineticEnergy << " "
588 // << std::setw( 9) << EnergyDifference << " "
589 << std::setw( 9) << step_len << " "
590 << std::setw( 9) << physStep << " "
591 << std::setw( 9) << safety << " ";
592 if( startVolume != 0) {
593 G4cout << std::setw(12) << startVolume->GetName() << " ";
594 } else {
595 G4cout << std::setw(12) << "OutOfWorld" << " ";
596 }
597
598#if 0
599 if( endVolume != 0)
600 {
601 G4cout << std::setw(12) << endVolume()->GetName() << " ";
602 }
603 else
604 {
605 G4cout << std::setw(12) << "OutOfWorld" << " ";
606 }
607#endif
608 G4cout << G4endl;
609 }
610}
611
612int readin_particle( )
613{
614 static const
615 double pmass[5] = {
616 0.00051099906 , // electron
617 0.105658389 , // muon
618 0.13956995 , // pion
619 0.493677 , // kaon
620 0.93827231 // proton
621 } ;
622 const double cSpeed = 299792458.0 ; // light speed in m/s
623 const double pi = 3.141592653589793238 ;
624 int pCharge, i ;
625 double pMomentum, pTeta, pPhi, h ;
626 G4cout<<"Enter particle type: 0 - electron, 1 - muon, 2 - pion, \n"
627 <<"3 - kaon, 4 - proton "<< G4endl ;
628 G4cin>>i ;
629 double pMass = pmass[i] ;
630 G4cout<<"Enter particle charge in units of the positron charge "<< G4endl ;
631 G4cin>>pCharge ;
632 G4cout<<"Enter particle momentum in GeV/c"<<G4endl ;
633 G4cin>>pMomentum ;
634 G4cout<<"Enter particle teta & phi in degrees"<<G4endl ;
635 G4cin>>pTeta ;
636 G4cin>>pPhi ;
637 G4cout<<"Enter particle Step in centimeters"<<G4endl ;
638 G4cin>>h ;
639
640 h *= 10.; // G4 units are in millimeters.
641
642 double betaGamma = pMomentum/pMass ;
643 double pSpeed = betaGamma*cSpeed/std::sqrt(1 + betaGamma*betaGamma) ;
644 double pEnergy = pMomentum*cSpeed/pSpeed ;
645 pEnergy *= 1.60217733e-10 ; // energy in J (SI units)
646 pTeta *= pi/180 ;
647 pPhi *= pi/180 ;
648
649#if 0
650 for(i=0;i<3;i++) ystart[i] = 0 ; // initial coordinates
651 ystart[3] = pSpeed*std::sin(pTeta)*std::cos(pPhi) ; // and speeds
652 ystart[4] = pSpeed*std::sin(pTeta)*std::sin(pPhi) ;
653 ystart[5] = pSpeed*std::cos(pTeta) ;
654#endif
655
656 return 1;
657}
658
659
Note: See TracBrowser for help on using the repository browser.