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-patch-03 $ >> 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