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

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

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