source: trunk/source/geometry/magneticfield/test/OtherFields/testDelphiField.cc@ 1199

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

nvx fichiers dans CVS
File size: 19.6 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: testDelphiField.cc,v 1.7 2006/06/29 18:26:40 gunter Exp $
28// GEANT4 tag $Name: HEAD $
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#include "G4DELPHIMagField.hh"
58#include "G4QuadrupoleMagField.hh"
59#include "G4HarmonicPolMagField.hh"
60// #include "G4LineCurrentMagField.hh"
61
62#include "G4ios.hh"
63#include <iomanip>
64
65// Sample Parameterisation
66class G4LinScale : public G4VPVParameterisation
67{
68 virtual void ComputeTransformation(const G4int n,
69 G4VPhysicalVolume* pRep) const
70 {
71 pRep->SetTranslation(G4ThreeVector(0,(n-1)*15,0));
72 }
73
74 virtual void ComputeDimensions(G4Box &pBox,
75 const G4int n,
76 const G4VPhysicalVolume* pRep) const
77 {
78 pBox.SetXHalfLength(10);
79 pBox.SetYHalfLength(5+n);
80 pBox.SetZHalfLength(5+n);
81 }
82
83 virtual void ComputeDimensions(G4Tubs &,
84 const G4int ,
85 const G4VPhysicalVolume*) const {}
86 virtual void ComputeDimensions(G4Trd &,
87 const G4int,
88 const G4VPhysicalVolume*) const {}
89 virtual void ComputeDimensions(G4Cons &,
90 const G4int ,
91 const G4VPhysicalVolume*) const {}
92 virtual void ComputeDimensions(G4Trap &,
93 const G4int ,
94 const G4VPhysicalVolume*) const {}
95 virtual void ComputeDimensions(G4Hype &,
96 const G4int ,
97 const G4VPhysicalVolume*) const {}
98 virtual void ComputeDimensions(G4Sphere &,
99 const G4int ,
100 const G4VPhysicalVolume*) const {}
101 virtual void ComputeDimensions(G4Torus &,
102 const G4int ,
103 const G4VPhysicalVolume*) const {}
104 virtual void ComputeDimensions(G4Para &,
105 const G4int ,
106 const G4VPhysicalVolume*) const {}
107};
108G4LinScale myParam;
109
110// Build simple geometry:
111// 4 small cubes + 1 slab (all G4Boxes) are positioned inside a larger cuboid
112G4VPhysicalVolume* BuildGeometry()
113{
114
115 G4Box *myHugeBox= new G4Box("huge box",15*m,15*m,25*m);
116 G4Box *myBigBox= new G4Box("big cube",10*m,10*m,10*m);
117 G4Box *mySmallBox= new G4Box("smaller cube",2.5*m,2.5*m,2.5*m);
118 G4Box *myTinyBox= new G4Box("tiny cube",.25*m,.25*m,.25*m);
119
120 // G4Box *myVariableBox=
121 new G4Box("Variable Box",10,5,5);
122
123 // World Volume
124 //
125 G4LogicalVolume *worldLog=new G4LogicalVolume(myHugeBox,0,
126 "World",0,0,0);
127 // Logical with no material,field,
128 // sensitive detector or user limits
129
130 G4PVPlacement *worldPhys=new
131 G4PVPlacement(0,G4ThreeVector(0,0,0), "World",worldLog,
132 0,false,0);
133 // Note: no mother pointer set
134
135// Create the logical Volumes
136//
137// G4LogicalVolume(*pSolid, *pMaterial, Name, *pField, *pSDetector, *pULimits)
138//
139 G4LogicalVolume *BigBoxLog=new G4LogicalVolume(myBigBox,0,
140 "Crystal Box (large)",0,0,0);
141 G4LogicalVolume *smallBoxLog=new G4LogicalVolume(mySmallBox,0,
142 "Crystal Box (small)");
143 G4LogicalVolume *tinyBoxLog=new G4LogicalVolume(myTinyBox,0,
144 "Crystal Box (tiny)");
145
146
147// Place them.
148//
149// 1) Two big boxes in the world volume
150//
151 // G4PVPlacement *BigTg1Phys=
152 new G4PVPlacement(0,G4ThreeVector(0,0,-15*m),
153 "Big Target 1",BigBoxLog,
154 worldPhys,false,0);
155 // G4PVPlacement *BigTg2Phys=
156 new G4PVPlacement(0,G4ThreeVector(0,0, 15*m),
157 "Big Target 2",BigBoxLog,
158 worldPhys,false,0);
159
160// 2) Four (medium) boxes in X & Y near the origin of the world volume
161//
162 // G4PVPlacement *MedTg3a_Phys=
163 new G4PVPlacement(0,G4ThreeVector(0, 7.5*m,0),
164 "Target 3a",smallBoxLog,
165 worldPhys,false,0);
166 // G4PVPlacement *MedTg3b_Phys=
167 new G4PVPlacement(0,G4ThreeVector(0,-7.5*m,0),
168 "Target 3b",smallBoxLog,
169 worldPhys,false,0);
170 // G4PVPlacement *MedTg3c_Phys=
171 new G4PVPlacement(0,G4ThreeVector(-7.5*m,0,0),
172 "Target 3c",smallBoxLog,
173 worldPhys,false,0);
174 // G4PVPlacement *MedTg3d_Phys=
175 new G4PVPlacement(0,G4ThreeVector( 7.5*m,0,0),
176 "Target 3d",smallBoxLog,
177 worldPhys,false,0);
178
179
180// 3) Eight small boxes around the origin of the world volume
181// (in +-X, +-Y & +-Z)
182//
183 // G4PVPlacement *SmTg4a_Phys=
184 new G4PVPlacement
185 (0,G4ThreeVector( 0.3*m, 0.3*m,0.3*m), "Target 4a",tinyBoxLog,
186 worldPhys,false,0);
187 // G4PVPlacement *SmTg4b_Phys=
188 new G4PVPlacement
189 (0,G4ThreeVector( 0.3*m,-0.3*m,0.3*m), "Target 4b",tinyBoxLog,
190 worldPhys,false,0);
191 // G4PVPlacement *SmTg4c_Phys=
192 new G4PVPlacement
193 (0,G4ThreeVector(-0.3*m,-0.3*m,0.3*m), "Target 4c",tinyBoxLog,
194 worldPhys,false,0);
195 // G4PVPlacement *SmTg4d_Phys=
196 new G4PVPlacement
197 (0,G4ThreeVector(-0.3*m, 0.3*m,0.3*m), "Target 4d",tinyBoxLog,
198 worldPhys,false,0);
199
200 // G4PVPlacement *SmTg4e_Phys=
201 new G4PVPlacement
202 (0,G4ThreeVector( 0.3*m, 0.3*m,-0.3*m), "Target 4e",tinyBoxLog,
203 worldPhys,false,0);
204 // G4PVPlacement *SmTg4f_Phys=
205 new G4PVPlacement
206 (0,G4ThreeVector( 0.3*m,-0.3*m,-0.3*m), "Target 4f",tinyBoxLog,
207 worldPhys,false,0);
208 // G4PVPlacement *SmTg4g_Phys=
209 new G4PVPlacement
210 (0,G4ThreeVector(-0.3*m,-0.3*m,-0.3*m), "Target 4g",tinyBoxLog,
211 worldPhys,false,0);
212 // G4PVPlacement *SmTg4h_Phys=
213 new G4PVPlacement
214 (0,G4ThreeVector(-0.3*m, 0.3*m,-0.3*m), "Target 4h",tinyBoxLog,
215 worldPhys,false,0);
216
217 return worldPhys;
218}
219
220#include "G4ChordFinder.hh"
221#include "G4PropagatorInField.hh"
222#include "G4MagneticField.hh"
223#include "G4FieldManager.hh"
224#include "G4TransportationManager.hh"
225#include "G4HelixExplicitEuler.hh"
226#include "G4HelixSimpleRunge.hh"
227#include "G4HelixImplicitEuler.hh"
228#include "G4ExplicitEuler.hh"
229#include "G4ImplicitEuler.hh"
230#include "G4SimpleRunge.hh"
231#include "G4SimpleHeum.hh"
232#include "G4ClassicalRK4.hh"
233#include "G4Mag_UsualEqRhs.hh"
234#include "G4CashKarpRKF45.hh"
235#include "G4RKG3_Stepper.hh"
236
237G4MagneticField* pMagField;
238
239G4FieldManager* SetupField(G4int type)
240{
241 G4FieldManager *pFieldMgr;
242 G4ChordFinder *pChordFinder;
243
244 pMagField = new G4DELPHIMagField();
245 // pMagField = new G4HarmonicPolMagField();
246
247 G4Mag_UsualEqRhs *fEquation = new G4Mag_UsualEqRhs(pMagField);
248 G4MagIntegratorStepper *pStepper;
249
250 switch ( type )
251 {
252 case 0: pStepper = new G4ExplicitEuler( fEquation ); break;
253 case 1: pStepper = new G4ImplicitEuler( fEquation ); break;
254 case 2: pStepper = new G4SimpleRunge( fEquation ); break;
255 case 3: pStepper = new G4SimpleHeum( fEquation ); break;
256 case 4: pStepper = new G4ClassicalRK4( fEquation ); break;
257 case 5: pStepper = new G4HelixExplicitEuler( fEquation ); break;
258 case 6: pStepper = new G4HelixImplicitEuler( fEquation ); break;
259 case 7: pStepper = new G4HelixSimpleRunge( fEquation ); break;
260 case 8: pStepper = new G4CashKarpRKF45( fEquation ); break;
261 case 9: pStepper = new G4RKG3_Stepper( fEquation ); break;
262 default: pStepper = 0;
263 }
264
265 pFieldMgr= G4TransportationManager::GetTransportationManager()->
266 GetFieldManager();
267
268 pFieldMgr->SetDetectorField( pMagField );
269
270 pChordFinder = new G4ChordFinder( pMagField,
271 1.0e-2 * mm,
272 pStepper);
273
274 pFieldMgr->SetChordFinder( pChordFinder );
275
276 return pFieldMgr;
277}
278
279G4PropagatorInField* SetupPropagator( G4int type)
280{
281 // G4FieldManager* fieldMgr=
282 SetupField( type) ;
283
284 // G4ChordFinder theChordFinder( &MagField, 0.05*mm ); // Default stepper
285
286 G4PropagatorInField *thePropagator =
287 G4TransportationManager::GetTransportationManager()->
288 GetPropagatorInField ();
289
290 // Let us test the new Minimum Epsilon Step functionality
291 thePropagator -> SetMinimumEpsilonStep( 1.0e-4 ) ;
292
293 return thePropagator;
294}
295
296// This is Done only for this test program ... the transportation does it.
297//
298void SetChargeMomentumMass(G4double charge, G4double MomentumXc, G4double Mass)
299{
300 G4ChordFinder* pChordFinder;
301
302 pChordFinder= G4TransportationManager::GetTransportationManager()->
303 GetFieldManager()->GetChordFinder();
304
305 // pMagFieldPropagator->set_magnetic_field();
306 pChordFinder->SetChargeMomentumMass(
307 charge, // charge in e+ units
308 MomentumXc, // Momentum in Mev/c ?
309 Mass );
310}
311
312//
313// Test Stepping
314//
315G4bool testG4PropagatorInField(G4VPhysicalVolume *pTopNode, G4int type)
316{
317 void report_endPV(G4ThreeVector Position,
318 G4ThreeVector UnitVelocity,
319 G4double step_len,
320 G4double physStep,
321 G4double safety,
322 G4ThreeVector EndPosition,
323 G4ThreeVector EndUnitVelocity,
324 G4int Step,
325 G4VPhysicalVolume* startVolume);
326
327 G4UniformMagField MagField(10.*tesla, 0., 0.); // Tesla Defined ?
328 G4Navigator *pNavig= G4TransportationManager::
329 GetTransportationManager()-> GetNavigatorForTracking();
330 G4PropagatorInField *pMagFieldPropagator= SetupPropagator(type);
331
332 SetChargeMomentumMass( +1., // charge in e+ units
333 0.5 * proton_mass_c2, // Momentum in Mev/c ?
334 proton_mass_c2 );
335 pNavig->SetWorldVolume(pTopNode);
336
337
338 G4VPhysicalVolume *located;
339 G4double step_len, physStep, safety;
340 G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
341 G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
342
343 // physStep=kInfinity;
344 G4ThreeVector Position(0.,0.,0.);
345 G4ThreeVector UnitMomentum(0.,0.6,0.8);
346 G4ThreeVector EndPosition, EndUnitMomentum;
347
348//
349// Test location & Step computation
350//
351 /* assert(located->GetName()=="World"); */
352 if( std::fabs(UnitMomentum.mag() - 1.0) > 1.e-8 )
353 {
354 G4cerr << "UnitMomentum.mag() - 1.0 = " << UnitMomentum.mag() - 1.0 <<
355 G4endl;
356 }
357
358 G4cout << G4endl;
359
360 for( int iparticle=0; iparticle < 2; iparticle++ )
361 {
362 physStep= 2.5 * mm ; // millimeters
363 Position = G4ThreeVector(0.,0.,0.)
364 + iparticle * G4ThreeVector(0.2, 0.3, 0.4);
365 UnitMomentum = (G4ThreeVector(0.,0.6,0.8)
366 + (float)iparticle * G4ThreeVector(0.1, 0.2, 0.3)).unit();
367
368 G4double momentum = (0.5+iparticle*10.0) * proton_mass_c2;
369
370 G4double kineticEnergy = momentum*momentum /
371 ( std::sqrt( momentum*momentum + proton_mass_c2 * proton_mass_c2 )
372 + proton_mass_c2 );
373 G4double velocity = momentum / ( proton_mass_c2 + kineticEnergy );
374 G4double labTof= 10.0*ns, properTof= 0.1*ns;
375 G4ThreeVector Spin(1.0, 0.0, 0.0);
376 // Momentum in Mev/c ?
377 SetChargeMomentumMass(
378 +1, // charge in e+ units
379 momentum,
380 proton_mass_c2);
381 G4cout << G4endl;
382 G4cout << "Test PropagateMagField: ***********************" << G4endl
383 << " Starting New Particle with Position " << Position << G4endl
384 << " and UnitVelocity " << UnitMomentum << G4endl;
385 G4cout << " Momentum in GeV/c is "<< (0.5+iparticle*10.0)*proton_mass_c2;
386 G4cout << G4endl;
387
388
389 for( int istep=0; istep < 14; istep++ ){
390 // // G4cerr << "UnitMomentum Magnitude is " << UnitMomentum.mag() << G4endl;
391 located= pNavig->LocateGlobalPointAndSetup(Position);
392 // Is the following better ?? It would need "changes"
393 // located= pMagFieldPropagator->LocateGlobalPointAndSetup(Position);
394 // G4cerr << "Starting Step " << istep << " in volume "
395 // << located->GetName() << G4endl;
396
397 G4FieldTrack initTrack( Position,
398 UnitMomentum,
399 0.0, // starting S curve len
400 kineticEnergy,
401 proton_mass_c2,
402 velocity,
403 labTof,
404 properTof,
405 0 // or &Spin
406 );
407
408 step_len=pMagFieldPropagator->ComputeStep( initTrack,
409 physStep,
410 safety
411#ifdef G4MAG_CHECK_VOLUME
412 ,located);
413#else
414 );
415#endif
416 // --------------------
417 EndPosition= pMagFieldPropagator->EndPosition();
418 EndUnitMomentum= pMagFieldPropagator->EndMomentumDir();
419 // --------
420
421 if( std::fabs(EndUnitMomentum.mag2() - 1.0) > 1.e-8 )
422 G4cerr << "EndUnitMomentum.mag2() - 1.0 = " <<
423 EndUnitMomentum.mag2() - 1.0 << G4endl;
424
425 G4ThreeVector MoveVec = EndPosition - Position;
426 assert( MoveVec.mag() < physStep*(1.+1.e-9) );
427
428 // G4cout << " testPropagatorInField: After stepI " << istep << " : " << G4endl;
429 report_endPV(Position, UnitMomentum, step_len, physStep, safety,
430 EndPosition, EndUnitMomentum, istep, located );
431
432 assert(safety>=0);
433 pNavig->SetGeometricallyLimitedStep();
434 // pMagFieldPropagator->SetGeometricallyLimitedStep();
435
436 Position= EndPosition;
437 UnitMomentum= EndUnitMomentum;
438 physStep *= 2.;
439 } // ........................... end for ( istep )
440 } // .............................. end for ( iparticle )
441
442 return(1);
443}
444
445// int main(int argc, char** argv)
446int main(int argc, char **argv)
447{
448 G4VPhysicalVolume *myTopNode;
449 G4int type;
450 myTopNode=BuildGeometry(); // Build the geometry
451 G4GeometryManager::GetInstance()->CloseGeometry(false);
452
453 type = 8 ;
454
455 if( argc == 2 )
456 type = atoi(argv[1]);
457
458 testG4PropagatorInField(myTopNode, type);
459
460// Repeat tests but with full voxels
461 G4GeometryManager::GetInstance()->OpenGeometry();
462 G4GeometryManager::GetInstance()->CloseGeometry(true);
463
464 testG4PropagatorInField(myTopNode, type);
465
466 G4GeometryManager::GetInstance()->OpenGeometry();
467 return 0;
468}
469
470
471void report_endPV(G4ThreeVector Position,
472 G4ThreeVector UnitVelocity,
473 G4double step_len,
474 G4double physStep,
475 G4double safety,
476 G4ThreeVector EndPosition,
477 G4ThreeVector EndUnitVelocity,
478 G4int Step,
479 G4VPhysicalVolume* startVolume)
480 // G4VPhysicalVolume* endVolume)
481{
482 const G4int verboseLevel=1;
483
484 if( Step == 0 && verboseLevel <= 3 )
485 {
486 G4cout.precision(3);
487 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
488 G4cout << std::setw( 5) << "Step#" << " "
489 << std::setw( 9) << "X(mm)" << " "
490 << std::setw( 9) << "Y(mm)" << " "
491 << std::setw( 9) << "Z(mm)" << " "
492 << std::setw( 7) << " N_x " << " "
493 << std::setw( 7) << " N_y " << " "
494 << std::setw( 7) << " N_z " << " "
495 // << std::setw( 9) << "KinE(MeV)" << " "
496 // << std::setw( 9) << "dE(MeV)" << " "
497 << std::setw( 9) << "StepLen" << " "
498 << std::setw( 9) << "PhsStep" << " "
499 << std::setw( 9) << "Safety" << " "
500 << std::setw(18) << "NextVolume" << " "
501 << G4endl;
502 }
503 //
504 //
505 if( verboseLevel > 3 )
506 {
507 G4cout << "End Position is " << EndPosition << G4endl
508 << " and UnitVelocity is " << EndUnitVelocity << G4endl;
509 G4cout << "Step taken was " << step_len
510 << " out of PhysicalStep= " << physStep << G4endl;
511 G4cout << "Final safety is: " << safety << G4endl;
512
513 G4cout << "Chord length = " << (EndPosition-Position).mag() << G4endl;
514 G4cout << G4endl;
515 }
516 else // if( verboseLevel > 0 )
517 {
518 G4cout.precision(3);
519 G4cout << std::setw( 5) << Step << " "
520 << std::setw( 9) << Position.x() << " "
521 << std::setw( 9) << Position.y() << " "
522 << std::setw( 9) << Position.z() << " "
523 << std::setw( 7) << EndUnitVelocity.x() << " "
524 << std::setw( 7) << EndUnitVelocity.y() << " "
525 << std::setw( 7) << EndUnitVelocity.z() << " "
526 // << std::setw( 9) << KineticEnergy << " "
527 // << std::setw( 9) << EnergyDifference << " "
528 << std::setw( 9) << step_len << " "
529 << std::setw( 9) << physStep << " "
530 << std::setw( 9) << safety << " ";
531 if( startVolume != 0) {
532 G4cout << std::setw(12) << startVolume->GetName() << " ";
533 } else {
534 G4cout << std::setw(12) << "OutOfWorld" << " ";
535 }
536
537#if 0
538 if( endVolume != 0)
539 {
540 G4cout << std::setw(12) << endVolume()->GetName() << " ";
541 }
542 else
543 {
544 G4cout << std::setw(12) << "OutOfWorld" << " ";
545 }
546#endif
547 G4cout << G4endl;
548 }
549}
550
551int readin_particle( )
552{
553 static const
554 double pmass[5] = {
555 0.00051099906 , // electron
556 0.105658389 , // muon
557 0.13956995 , // pion
558 0.493677 , // kaon
559 0.93827231 // proton
560 } ;
561 const double cSpeed = 299792458.0 ; // light speed in m/s
562 const double pi = 3.141592653589793238 ;
563 int pCharge, i ;
564 double pMomentum, pTeta, pPhi, h ;
565 G4cout<<"Enter particle type: 0 - electron, 1 - muon, 2 - pion, \n"
566 <<"3 - kaon, 4 - proton "<< G4endl ;
567 G4cin>>i ;
568 double pMass = pmass[i] ;
569 G4cout<<"Enter particle charge in units of the positron charge "<< G4endl ;
570 G4cin>>pCharge ;
571 G4cout<<"Enter particle momentum in GeV/c"<<G4endl ;
572 G4cin>>pMomentum ;
573 G4cout<<"Enter particle teta & phi in degrees"<<G4endl ;
574 G4cin>>pTeta ;
575 G4cin>>pPhi ;
576 G4cout<<"Enter particle Step in centimeters"<<G4endl ;
577 G4cin>>h ;
578
579 h *= 10.; // G4 units are in millimeters.
580
581 double betaGamma = pMomentum/pMass ;
582 double pSpeed = betaGamma*cSpeed/std::sqrt(1 + betaGamma*betaGamma) ;
583 double pEnergy = pMomentum*cSpeed/pSpeed ;
584 pEnergy *= 1.60217733e-10 ; // energy in J (SI units)
585 pTeta *= pi/180 ;
586 pPhi *= pi/180 ;
587
588#if 0
589 for(i=0;i<3;i++) ystart[i] = 0 ; // initial coordinates
590 ystart[3] = pSpeed*std::sin(pTeta)*std::cos(pPhi) ; // and speeds
591 ystart[4] = pSpeed*std::sin(pTeta)*std::sin(pPhi) ;
592 ystart[5] = pSpeed*std::cos(pTeta) ;
593#endif
594
595 return 1;
596}
597
598
Note: See TracBrowser for help on using the repository browser.