| 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 26 // G4FieldTrack implementation
27 // 27 //
28 // Author: John Apostolakis, CERN - First vers 28 // Author: John Apostolakis, CERN - First version, 14.10.1996
29 // ------------------------------------------- 29 // -------------------------------------------------------------------
30 30
31 #include "G4FieldTrack.hh" 31 #include "G4FieldTrack.hh"
32 32
33 std::ostream& operator<<( std::ostream& os, co 33 std::ostream& operator<<( std::ostream& os, const G4FieldTrack& SixVec)
34 { 34 {
35 const G4double* SixV = SixVec.SixVector; 35 const G4double* SixV = SixVec.SixVector;
36 const G4int precPos= 9; // For position 36 const G4int precPos= 9; // For position
37 const G4int precEp= 9; // For Energy / 37 const G4int precEp= 9; // For Energy / momentum
38 const G4int precLen= 12; // For Length a 38 const G4int precLen= 12; // For Length along track
39 const G4int precSpin= 9; // For polarisa 39 const G4int precSpin= 9; // For polarisation
40 const G4int precTime= 6; // For time of 40 const G4int precTime= 6; // For time of flight
41 const G4long oldpr= os.precision(precPos) << 41 const G4int oldpr= os.precision(precPos);
42 os << " ( "; 42 os << " ( ";
43 os << " X= " << SixV[0] << " " << SixV[1] 43 os << " X= " << SixV[0] << " " << SixV[1] << " "
44 << SixV[2] << " "; // Posit 44 << SixV[2] << " "; // Position
45 os.precision(precEp); 45 os.precision(precEp);
46 os << " P= " << SixV[3] << " " << SixV[4] 46 os << " P= " << SixV[3] << " " << SixV[4] << " "
47 << SixV[5] << " "; // Momen 47 << SixV[5] << " "; // Momentum
48 os << " Pmag= " 48 os << " Pmag= "
49 << G4ThreeVector(SixV[3], SixV[4], Six 49 << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude
50 os << " Ekin= " << SixVec.fKineticEnergy 50 os << " Ekin= " << SixVec.fKineticEnergy ;
51 os.precision(precLen); 51 os.precision(precLen);
52 os << " l= " << SixVec.GetCurveLength(); 52 os << " l= " << SixVec.GetCurveLength();
53 os.precision(6); 53 os.precision(6);
54 os << " m0= " << SixVec.fRestMass_c2; 54 os << " m0= " << SixVec.fRestMass_c2;
55 os << " (Pdir-1)= " << SixVec.fMomentumD 55 os << " (Pdir-1)= " << SixVec.fMomentumDir.mag()-1.0;
56 if( SixVec.fLabTimeOfFlight > 0.0 ) 56 if( SixVec.fLabTimeOfFlight > 0.0 )
57 { 57 {
58 os.precision(precTime); 58 os.precision(precTime);
59 } 59 }
60 else 60 else
61 { 61 {
62 os.precision(3); 62 os.precision(3);
63 } 63 }
64 os << " t_lab= " << SixVec.fLabTimeOfF 64 os << " t_lab= " << SixVec.fLabTimeOfFlight;
65 os << " t_proper= " << SixVec.fProperTime 65 os << " t_proper= " << SixVec.fProperTimeOfFlight ;
66 G4ThreeVector pol= SixVec.GetPolarization 66 G4ThreeVector pol= SixVec.GetPolarization();
67 if( pol.mag2() > 0.0 ) 67 if( pol.mag2() > 0.0 )
68 { 68 {
69 os.precision(precSpin); 69 os.precision(precSpin);
70 os << " PolV= " << pol; // SixVec.GetP 70 os << " PolV= " << pol; // SixVec.GetPolarization();
71 } 71 }
72 else 72 else
73 { 73 {
74 os << " PolV= (0,0,0) "; 74 os << " PolV= (0,0,0) ";
75 } 75 }
76 os << " ) "; 76 os << " ) ";
77 os.precision(oldpr); 77 os.precision(oldpr);
78 return os; 78 return os;
79 } 79 }
80 80
81 G4FieldTrack::G4FieldTrack( const G4ThreeVecto 81 G4FieldTrack::G4FieldTrack( const G4ThreeVector& pPosition,
82 G4double LaboratoryTimeO 82 G4double LaboratoryTimeOfFlight,
83 const G4ThreeVector& pMomentumDirect 83 const G4ThreeVector& pMomentumDirection,
84 G4double kineticEnergy, 84 G4double kineticEnergy,
85 G4double restMass_c2, 85 G4double restMass_c2,
86 G4double charge, 86 G4double charge,
87 const G4ThreeVector& vecPolarization 87 const G4ThreeVector& vecPolarization,
88 G4double magnetic_dipole 88 G4double magnetic_dipole_moment,
89 G4double 89 G4double curve_length,
90 G4double 90 G4double pdgSpin )
91 : fDistanceAlongCurve(curve_length), 91 : fDistanceAlongCurve(curve_length),
92 fKineticEnergy(kineticEnergy), 92 fKineticEnergy(kineticEnergy),
93 fRestMass_c2(restMass_c2), 93 fRestMass_c2(restMass_c2),
94 fLabTimeOfFlight(LaboratoryTimeOfFlight), 94 fLabTimeOfFlight(LaboratoryTimeOfFlight),
95 fProperTimeOfFlight(0.), 95 fProperTimeOfFlight(0.),
96 // fMomentumDir(pMomentumDirection), 96 // fMomentumDir(pMomentumDirection),
97 fChargeState( charge, magnetic_dipole_mome 97 fChargeState( charge, magnetic_dipole_moment, pdgSpin )
98 // fChargeState( charge, magnetic_dipole_m 98 // fChargeState( charge, magnetic_dipole_moment ) ,
99 // fPDGSpin( pdgSpin ) 99 // fPDGSpin( pdgSpin )
100 { 100 {
101 UpdateFourMomentum( kineticEnergy, pMomentum 101 UpdateFourMomentum( kineticEnergy, pMomentumDirection );
102 // Sets momentum direction as well. 102 // Sets momentum direction as well.
103 103
104 SetPosition( pPosition ); 104 SetPosition( pPosition );
105 SetPolarization( vecPolarization ); 105 SetPolarization( vecPolarization );
106 } 106 }
107 107
108 G4FieldTrack::G4FieldTrack( const G4ThreeVecto 108 G4FieldTrack::G4FieldTrack( const G4ThreeVector& pPosition,
109 const G4ThreeVecto 109 const G4ThreeVector& pMomentumDirection,
110 G4double 110 G4double curve_length,
111 G4double 111 G4double kineticEnergy,
112 const G4double 112 const G4double restMass_c2,
113 G4double, 113 G4double, // velocity
114 G4double 114 G4double pLaboratoryTimeOfFlight,
115 G4double 115 G4double pProperTimeOfFlight,
116 const G4ThreeVecto 116 const G4ThreeVector* pPolarization,
117 G4double 117 G4double pdgSpin )
118 : fDistanceAlongCurve(curve_length), 118 : fDistanceAlongCurve(curve_length),
119 fKineticEnergy(kineticEnergy), 119 fKineticEnergy(kineticEnergy),
120 fRestMass_c2(restMass_c2), 120 fRestMass_c2(restMass_c2),
121 fLabTimeOfFlight(pLaboratoryTimeOfFlight), 121 fLabTimeOfFlight(pLaboratoryTimeOfFlight),
122 fProperTimeOfFlight(pProperTimeOfFlight), 122 fProperTimeOfFlight(pProperTimeOfFlight),
123 fChargeState( DBL_MAX, DBL_MAX, -1.0 ) // 123 fChargeState( DBL_MAX, DBL_MAX, -1.0 ) // charge not set
124 { 124 {
125 UpdateFourMomentum( kineticEnergy, pMomentum 125 UpdateFourMomentum( kineticEnergy, pMomentumDirection );
126 // Sets momentum direction as well. 126 // Sets momentum direction as well.
127 127
128 SetPosition( pPosition ); 128 SetPosition( pPosition );
129 fChargeState.SetPDGSpin( pdgSpin ); 129 fChargeState.SetPDGSpin( pdgSpin );
130 130
131 G4ThreeVector PolarVec(0.0, 0.0, 0.0); 131 G4ThreeVector PolarVec(0.0, 0.0, 0.0);
132 if( pPolarization != nullptr ) { PolarVec= << 132 if( pPolarization ) { PolarVec= *pPolarization; }
133 SetPolarization( PolarVec ); 133 SetPolarization( PolarVec );
134 } 134 }
135 135
136 G4FieldTrack::G4FieldTrack( char ) 136 G4FieldTrack::G4FieldTrack( char ) // Nothing is set !!
137 : fKineticEnergy(0.), fRestMass_c2(0.), fLab 137 : fKineticEnergy(0.), fRestMass_c2(0.), fLabTimeOfFlight(0.),
138 fProperTimeOfFlight(0.), fChargeState( DBL 138 fProperTimeOfFlight(0.), fChargeState( DBL_MAX , DBL_MAX, -1 )
139 { 139 {
140 G4ThreeVector Zero(0.0, 0.0, 0.0); 140 G4ThreeVector Zero(0.0, 0.0, 0.0);
141 SetCurvePnt( Zero, Zero, 0.0 ); 141 SetCurvePnt( Zero, Zero, 0.0 );
142 SetPolarization( Zero ); 142 SetPolarization( Zero );
143 // fInitialMomentumMag = 0.00; // Invalid 143 // fInitialMomentumMag = 0.00; // Invalid
144 // fLastMomentumMag = 0.0; 144 // fLastMomentumMag = 0.0;
145 } 145 }
146 146
147 void G4FieldTrack:: 147 void G4FieldTrack::
148 SetChargeAndMoments(G4double charge, 148 SetChargeAndMoments(G4double charge,
149 G4double magnetic_dipole_moment, // def 149 G4double magnetic_dipole_moment, // default = DBL_MAX
150 G4double electric_dipole_moment, // dit 150 G4double electric_dipole_moment, // ditto
151 G4double magnetic_charge ) // dit 151 G4double magnetic_charge ) // ditto
152 { 152 {
153 fChargeState.SetChargesAndMoments( charge, 153 fChargeState.SetChargesAndMoments( charge,
154 magnetic_ 154 magnetic_dipole_moment,
155 electric_ 155 electric_dipole_moment,
156 magnetic_ 156 magnetic_charge );
157 157
158 // NOTE: Leaves Spin unchanged ! 158 // NOTE: Leaves Spin unchanged !
159 // 159 //
160 // G4double pdgSpin= fChargeState.GetSpin(); 160 // G4double pdgSpin= fChargeState.GetSpin();
161 // New Property of ChargeState (not well doc 161 // New Property of ChargeState (not well documented! )
162 162
163 // IDEA: Improve the implementation using ha 163 // IDEA: Improve the implementation using handles
164 // -- and handle to the old one (which can 164 // -- and handle to the old one (which can be shared by other copies) and
165 // must not be left to hang loose 165 // must not be left to hang loose
166 // 166 //
167 // fpChargeState= new G4ChargeState( charge 167 // fpChargeState= new G4ChargeState( charge, magnetic_dipole_moment,
168 // electric_dipole_moment, magneti 168 // electric_dipole_moment, magnetic_charge );
169 } 169 }
170 170
171 // Load values from array 171 // Load values from array
172 // 172 //
173 // Note that momentum direction must-be/is nor 173 // Note that momentum direction must-be/is normalised
174 // 174 //
175 void G4FieldTrack::LoadFromArray(const G4doubl 175 void G4FieldTrack::LoadFromArray(const G4double valArrIn[ncompSVEC],
176 G4int n 176 G4int noVarsIntegrated)
177 { 177 {
178 // Fill the variables not integrated with ze 178 // Fill the variables not integrated with zero -- so it's clear !!
179 // 179 //
180 G4double valArr[ncompSVEC]; 180 G4double valArr[ncompSVEC];
181 for(G4int i=0; i<noVarsIntegrated; ++i) 181 for(G4int i=0; i<noVarsIntegrated; ++i)
182 { 182 {
183 valArr[i] = valArrIn[i]; 183 valArr[i] = valArrIn[i];
184 } 184 }
185 for(G4int i=noVarsIntegrated; i<ncompSVEC; + 185 for(G4int i=noVarsIntegrated; i<ncompSVEC; ++i)
186 { 186 {
187 valArr[i] = 0.0; 187 valArr[i] = 0.0;
188 } 188 }
189 189
190 SixVector[0] = valArr[0]; 190 SixVector[0] = valArr[0];
191 SixVector[1] = valArr[1]; 191 SixVector[1] = valArr[1];
192 SixVector[2] = valArr[2]; 192 SixVector[2] = valArr[2];
193 SixVector[3] = valArr[3]; 193 SixVector[3] = valArr[3];
194 SixVector[4] = valArr[4]; 194 SixVector[4] = valArr[4];
195 SixVector[5] = valArr[5]; 195 SixVector[5] = valArr[5];
196 196
197 G4ThreeVector Momentum(valArr[3],valArr[4],v 197 G4ThreeVector Momentum(valArr[3],valArr[4],valArr[5]);
198 198
199 G4double momentum_square= Momentum.mag2(); 199 G4double momentum_square= Momentum.mag2();
200 fMomentumDir= Momentum.unit(); 200 fMomentumDir= Momentum.unit();
201 201
202 fKineticEnergy = momentum_square 202 fKineticEnergy = momentum_square
203 / (std::sqrt(momentum_square+ 203 / (std::sqrt(momentum_square+fRestMass_c2*fRestMass_c2)
204 + fRestMass_c2 ); 204 + fRestMass_c2 );
205 // The above equation is stable for small 205 // The above equation is stable for small and large momenta
206 206
207 // The following components may or may not b 207 // The following components may or may not be
208 // integrated over -- integration is optiona 208 // integrated over -- integration is optional
209 // fKineticEnergy = valArr[6]; 209 // fKineticEnergy = valArr[6];
210 210
211 fLabTimeOfFlight = valArr[7]; 211 fLabTimeOfFlight = valArr[7];
212 fProperTimeOfFlight = valArr[8]; 212 fProperTimeOfFlight = valArr[8];
213 G4ThreeVector vecPolarization= G4ThreeVector 213 G4ThreeVector vecPolarization= G4ThreeVector(valArr[9],valArr[10],valArr[11]);
214 SetPolarization( vecPolarization ); 214 SetPolarization( vecPolarization );
215 215
216 // fMomentumDir=G4ThreeVector(valArr[13],val 216 // fMomentumDir=G4ThreeVector(valArr[13],valArr[14],valArr[15]);
217 // fDistanceAlongCurve= valArr[]; 217 // fDistanceAlongCurve= valArr[];
218 } 218 }
219 219