| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4FieldTrack implementation <<
27 // 26 //
28 // Author: John Apostolakis, CERN - First vers << 27 // $Id: G4FieldTrack.cc,v 1.14 2007/10/03 15:34:42 japost Exp $
>> 28 // GEANT4 tag $Name: geant4-09-02 $
>> 29 //
29 // ------------------------------------------- 30 // -------------------------------------------------------------------
30 31
31 #include "G4FieldTrack.hh" 32 #include "G4FieldTrack.hh"
32 33
33 std::ostream& operator<<( std::ostream& os, co 34 std::ostream& operator<<( std::ostream& os, const G4FieldTrack& SixVec)
34 { 35 {
35 const G4double* SixV = SixVec.SixVector; << 36 const G4double *SixV = SixVec.SixVector;
36 const G4int precPos= 9; // For position <<
37 const G4int precEp= 9; // For Energy / <<
38 const G4int precLen= 12; // For Length a <<
39 const G4int precSpin= 9; // For polarisa <<
40 const G4int precTime= 6; // For time of <<
41 const G4long oldpr= os.precision(precPos) <<
42 os << " ( "; 37 os << " ( ";
43 os << " X= " << SixV[0] << " " << SixV[1] << 38 os << " X= " << SixV[0] << " " << SixV[1] << " " << SixV[2] << " "; // Position
44 << SixV[2] << " "; // Posit << 39 os << " V= " << SixV[3] << " " << SixV[4] << " " << SixV[5] << " "; // Momentum
45 os.precision(precEp); << 40 os << " v2= " << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude
46 os << " P= " << SixV[3] << " " << SixV[4] << 41 os << " mdm= " << SixVec.fMomentumDir.mag();
47 << SixV[5] << " "; // Momen <<
48 os << " Pmag= " <<
49 << G4ThreeVector(SixV[3], SixV[4], Six <<
50 os << " Ekin= " << SixVec.fKineticEnergy <<
51 os.precision(precLen); <<
52 os << " l= " << SixVec.GetCurveLength(); 42 os << " l= " << SixVec.GetCurveLength();
53 os.precision(6); <<
54 os << " m0= " << SixVec.fRestMass_c2; <<
55 os << " (Pdir-1)= " << SixVec.fMomentumD <<
56 if( SixVec.fLabTimeOfFlight > 0.0 ) <<
57 { <<
58 os.precision(precTime); <<
59 } <<
60 else <<
61 { <<
62 os.precision(3); <<
63 } <<
64 os << " t_lab= " << SixVec.fLabTimeOfF <<
65 os << " t_proper= " << SixVec.fProperTime <<
66 G4ThreeVector pol= SixVec.GetPolarization <<
67 if( pol.mag2() > 0.0 ) <<
68 { <<
69 os.precision(precSpin); <<
70 os << " PolV= " << pol; // SixVec.GetP <<
71 } <<
72 else <<
73 { <<
74 os << " PolV= (0,0,0) "; <<
75 } <<
76 os << " ) "; 43 os << " ) ";
77 os.precision(oldpr); <<
78 return os; 44 return os;
79 } 45 }
80 46
81 G4FieldTrack::G4FieldTrack( const G4ThreeVecto 47 G4FieldTrack::G4FieldTrack( const G4ThreeVector& pPosition,
82 G4double LaboratoryTimeO 48 G4double LaboratoryTimeOfFlight,
83 const G4ThreeVector& pMomentumDirect 49 const G4ThreeVector& pMomentumDirection,
84 G4double kineticEnergy, 50 G4double kineticEnergy,
85 G4double restMass_c2, 51 G4double restMass_c2,
86 G4double charge, 52 G4double charge,
87 const G4ThreeVector& vecPolarization << 53 const G4ThreeVector& Spin,
88 G4double magnetic_dipole 54 G4double magnetic_dipole_moment,
89 G4double << 55 G4double curve_length )
90 G4double << 56 : fKineticEnergy(kineticEnergy),
91 : fDistanceAlongCurve(curve_length), <<
92 fKineticEnergy(kineticEnergy), <<
93 fRestMass_c2(restMass_c2), 57 fRestMass_c2(restMass_c2),
94 fLabTimeOfFlight(LaboratoryTimeOfFlight), 58 fLabTimeOfFlight(LaboratoryTimeOfFlight),
95 fProperTimeOfFlight(0.), << 59 // fProperTimeOfFlight(0.0),
96 // fMomentumDir(pMomentumDirection), 60 // fMomentumDir(pMomentumDirection),
97 fChargeState( charge, magnetic_dipole_mome << 61 fChargeState( charge, magnetic_dipole_moment )
98 // fChargeState( charge, magnetic_dipole_m <<
99 // fPDGSpin( pdgSpin ) <<
100 { 62 {
101 UpdateFourMomentum( kineticEnergy, pMomentum << 63 G4double momentum = std::sqrt(kineticEnergy*kineticEnergy
102 // Sets momentum direction as well. << 64 +2.0*restMass_c2*kineticEnergy);
>> 65
>> 66 G4ThreeVector pMomentum= momentum * pMomentumDirection;
>> 67 SetCurvePnt( pPosition, pMomentum, curve_length );
>> 68 // Sets momentum direction as well.
>> 69
>> 70 // Set the momentum direction again - keeping value from argument exactly
>> 71 fMomentumDir=pMomentumDirection;
103 72
104 SetPosition( pPosition ); << 73 InitialiseSpin( Spin );
105 SetPolarization( vecPolarization ); << 74
>> 75 // fpChargeState = new G4ChargeState( charge, magnetic_dipole_moment );
106 } 76 }
107 77
108 G4FieldTrack::G4FieldTrack( const G4ThreeVecto 78 G4FieldTrack::G4FieldTrack( const G4ThreeVector& pPosition,
109 const G4ThreeVecto 79 const G4ThreeVector& pMomentumDirection,
110 G4double 80 G4double curve_length,
111 G4double 81 G4double kineticEnergy,
112 const G4double 82 const G4double restMass_c2,
113 G4double, 83 G4double, // velocity
114 G4double 84 G4double pLaboratoryTimeOfFlight,
115 G4double 85 G4double pProperTimeOfFlight,
116 const G4ThreeVecto << 86 const G4ThreeVector* pSpin)
117 G4double << 87 : fKineticEnergy(kineticEnergy),
118 : fDistanceAlongCurve(curve_length), <<
119 fKineticEnergy(kineticEnergy), <<
120 fRestMass_c2(restMass_c2), 88 fRestMass_c2(restMass_c2),
121 fLabTimeOfFlight(pLaboratoryTimeOfFlight), 89 fLabTimeOfFlight(pLaboratoryTimeOfFlight),
122 fProperTimeOfFlight(pProperTimeOfFlight), 90 fProperTimeOfFlight(pProperTimeOfFlight),
123 fChargeState( DBL_MAX, DBL_MAX, -1.0 ) // << 91 // fMomentumDir(pMomentumDirection),
>> 92 fChargeState( DBL_MAX ) // charge not set
124 { 93 {
125 UpdateFourMomentum( kineticEnergy, pMomentum << 94 G4double momentum = std::sqrt(kineticEnergy*kineticEnergy
126 // Sets momentum direction as well. << 95 +2.0*restMass_c2*kineticEnergy);
127 << 96 G4ThreeVector pMomentum= momentum * pMomentumDirection;
128 SetPosition( pPosition ); << 97
129 fChargeState.SetPDGSpin( pdgSpin ); << 98 SetCurvePnt( pPosition, pMomentum, curve_length );
130 << 99 // Sets momentum direction as well.
131 G4ThreeVector PolarVec(0.0, 0.0, 0.0); << 100
132 if( pPolarization != nullptr ) { PolarVec= << 101 // Set the momentum direction again
133 SetPolarization( PolarVec ); << 102 // -- to avoid numerical issues from multiplying by momentum and dividing again
>> 103 fMomentumDir=pMomentumDirection;
>> 104
>> 105 G4ThreeVector Spin(0.0, 0.0, 0.0);
>> 106 if( !pSpin ) Spin= G4ThreeVector(0.,0.,0.);
>> 107 else Spin= *pSpin;
>> 108 InitialiseSpin( Spin );
>> 109
>> 110 // fpChargeState = new G4ChargeState( DBL_MAX ); // charge not yet set !!
134 } 111 }
135 112
136 G4FieldTrack::G4FieldTrack( char ) 113 G4FieldTrack::G4FieldTrack( char ) // Nothing is set !!
137 : fKineticEnergy(0.), fRestMass_c2(0.), fLab << 114 : fRestMass_c2(0.0), fLabTimeOfFlight(0.0),
138 fProperTimeOfFlight(0.), fChargeState( DBL << 115 fChargeState( DBL_MAX ) // charge not set
139 { 116 {
140 G4ThreeVector Zero(0.0, 0.0, 0.0); 117 G4ThreeVector Zero(0.0, 0.0, 0.0);
141 SetCurvePnt( Zero, Zero, 0.0 ); 118 SetCurvePnt( Zero, Zero, 0.0 );
142 SetPolarization( Zero ); << 119 InitialiseSpin( Zero );
143 // fInitialMomentumMag = 0.00; // Invalid << 120 // fpChargeState = new G4ChargeState( DBL_MAX );
144 // fLastMomentumMag = 0.0; <<
145 } 121 }
146 122
147 void G4FieldTrack:: 123 void G4FieldTrack::
148 SetChargeAndMoments(G4double charge, 124 SetChargeAndMoments(G4double charge,
149 G4double magnetic_dipole_moment, // def << 125 G4double magnetic_dipole_moment, // default= DBL_MAX - do not change
150 G4double electric_dipole_moment, // dit << 126 G4double electric_dipole_moment, // ditto
151 G4double magnetic_charge ) // dit << 127 G4double magnetic_charge ) // ditto
152 { 128 {
153 fChargeState.SetChargesAndMoments( charge, << 129 fChargeState.SetChargeAndMoments( charge, magnetic_dipole_moment,
154 magnetic_ << 130 electric_dipole_moment, magnetic_charge );
155 electric_ <<
156 magnetic_ <<
157 131
158 // NOTE: Leaves Spin unchanged ! << 132 // fpChargeState->SetChargeAndMoments( charge, magnetic_dipole_moment,
159 // << 133 // electric_dipole_moment, magnetic_charge );
160 // G4double pdgSpin= fChargeState.GetSpin(); <<
161 // New Property of ChargeState (not well doc <<
162 134
163 // IDEA: Improve the implementation using ha << 135 // TO-DO: Improve the implementation using handles
164 // -- and handle to the old one (which can 136 // -- and handle to the old one (which can be shared by other copies) and
165 // must not be left to hang loose 137 // must not be left to hang loose
166 // 138 //
167 // fpChargeState= new G4ChargeState( charge 139 // fpChargeState= new G4ChargeState( charge, magnetic_dipole_moment,
168 // electric_dipole_moment, magneti 140 // electric_dipole_moment, magnetic_charge );
169 } <<
170 <<
171 // Load values from array <<
172 // <<
173 // Note that momentum direction must-be/is nor <<
174 // <<
175 void G4FieldTrack::LoadFromArray(const G4doubl <<
176 G4int n <<
177 { <<
178 // Fill the variables not integrated with ze <<
179 // <<
180 G4double valArr[ncompSVEC]; <<
181 for(G4int i=0; i<noVarsIntegrated; ++i) <<
182 { <<
183 valArr[i] = valArrIn[i]; <<
184 } <<
185 for(G4int i=noVarsIntegrated; i<ncompSVEC; + <<
186 { <<
187 valArr[i] = 0.0; <<
188 } <<
189 <<
190 SixVector[0] = valArr[0]; <<
191 SixVector[1] = valArr[1]; <<
192 SixVector[2] = valArr[2]; <<
193 SixVector[3] = valArr[3]; <<
194 SixVector[4] = valArr[4]; <<
195 SixVector[5] = valArr[5]; <<
196 <<
197 G4ThreeVector Momentum(valArr[3],valArr[4],v <<
198 <<
199 G4double momentum_square= Momentum.mag2(); <<
200 fMomentumDir= Momentum.unit(); <<
201 <<
202 fKineticEnergy = momentum_square <<
203 / (std::sqrt(momentum_square+ <<
204 + fRestMass_c2 ); <<
205 // The above equation is stable for small <<
206 <<
207 // The following components may or may not b <<
208 // integrated over -- integration is optiona <<
209 // fKineticEnergy = valArr[6]; <<
210 <<
211 fLabTimeOfFlight = valArr[7]; <<
212 fProperTimeOfFlight = valArr[8]; <<
213 G4ThreeVector vecPolarization= G4ThreeVector <<
214 SetPolarization( vecPolarization ); <<
215 <<
216 // fMomentumDir=G4ThreeVector(valArr[13],val <<
217 // fDistanceAlongCurve= valArr[]; <<
218 } 141 }
219 142