Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/error_propagation/src/G4ErrorSurfaceTrajState.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // ------------------------------------------------------------
 28 //      GEANT 4 class implementation file
 29 // ------------------------------------------------------------
 30 //
 31 
 32 #include "G4ErrorSurfaceTrajState.hh"
 33 #include "G4ErrorPropagatorData.hh"
 34 
 35 #include "G4PhysicalConstants.hh"
 36 #include "G4SystemOfUnits.hh"
 37 #include "G4Field.hh"
 38 #include "G4FieldManager.hh"
 39 #include "G4TransportationManager.hh"
 40 
 41 #include "G4ErrorMatrix.hh"
 42 
 43 #include <iomanip>
 44 
 45 //------------------------------------------------------------------------
 46 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajState(
 47   const G4String& partType, const G4Point3D& pos, const G4Vector3D& mom,
 48   const G4Vector3D& vecU, const G4Vector3D& vecV, const G4ErrorTrajErr& errmat)
 49   : G4ErrorTrajState(partType, pos, mom, errmat)
 50 {
 51   Init();
 52   fTrajParam = G4ErrorSurfaceTrajParam(pos, mom, vecU, vecV);
 53 }
 54 
 55 //------------------------------------------------------------------------
 56 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajState(const G4String& partType,
 57                                                  const G4Point3D& pos,
 58                                                  const G4Vector3D& mom,
 59                                                  const G4Plane3D& plane,
 60                                                  const G4ErrorTrajErr& errmat)
 61   : G4ErrorTrajState(partType, pos, mom, errmat)
 62 {
 63   Init();
 64   fTrajParam = G4ErrorSurfaceTrajParam(pos, mom, plane);
 65 }
 66 
 67 //------------------------------------------------------------------------
 68 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajState(G4ErrorFreeTrajState& tpSC,
 69                                                  const G4Plane3D& plane)
 70   : G4ErrorTrajState(tpSC.GetParticleType(), tpSC.GetPosition(),
 71                      tpSC.GetMomentum())
 72 {
 73   //  fParticleType = tpSC.GetParticleType();
 74   //  fPosition = tpSC.GetPosition();
 75   //  fMomentum = tpSC.GetMomentum();
 76   fTrajParam = G4ErrorSurfaceTrajParam(fPosition, fMomentum, plane);
 77   Init();
 78 
 79   //----- Get the error matrix in SC coordinates
 80   BuildErrorMatrix(tpSC, GetVectorV(), GetVectorW());
 81 }
 82 
 83 //------------------------------------------------------------------------
 84 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajState(G4ErrorFreeTrajState& tpSC,
 85                                                  const G4Vector3D& vecU,
 86                                                  const G4Vector3D& vecV,
 87                                                  G4ErrorMatrix& transfM)
 88   : G4ErrorTrajState(tpSC.GetParticleType(), tpSC.GetPosition(),
 89                      tpSC.GetMomentum())
 90 {
 91   Init();  // needed to define charge sign
 92   fTrajParam = G4ErrorSurfaceTrajParam(fPosition, fMomentum, vecU, vecV);
 93   //----- Get the error matrix in SC coordinates
 94   transfM = BuildErrorMatrix(tpSC, vecU, vecV);
 95 }
 96 
 97 //------------------------------------------------------------------------
 98 G4ErrorMatrix G4ErrorSurfaceTrajState::BuildErrorMatrix(
 99   G4ErrorFreeTrajState& tpSC, const G4Vector3D&, const G4Vector3D&)
100 {
101   G4double sclambda = tpSC.GetParameters().GetLambda();
102   G4double scphi    = tpSC.GetParameters().GetPhi();
103   if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() ==
104      G4ErrorMode_PropBackwards)
105   {
106     sclambda *= -1;
107     scphi += CLHEP::pi;
108   }
109   G4double cosLambda = std::cos(sclambda);
110   G4double sinLambda = std::sin(sclambda);
111   G4double sinPhi    = std::sin(scphi);
112   G4double cosPhi    = std::cos(scphi);
113 
114 #ifdef G4EVERBOSE
115   if(iverbose >= 4)
116     G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi "
117            << scphi << G4endl;
118 #endif
119 
120   G4ThreeVector vTN(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda);
121   G4ThreeVector vUN(-sinPhi, cosPhi, 0.);
122   G4ThreeVector vVN(-vTN.z() * vUN.y(), vTN.z() * vUN.x(), cosLambda);
123 
124 #ifdef G4EVERBOSE
125   if(iverbose >= 4)
126     G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN
127            << G4endl;
128 #endif
129   G4double UJ = vUN * GetVectorV();
130   G4double UK = vUN * GetVectorW();
131   G4double VJ = vVN * GetVectorV();
132   G4double VK = vVN * GetVectorW();
133 
134   //--- Get transformation first
135   G4ErrorMatrix transfM(5, 5, 0);
136   //--- Get magnetic field
137   const G4Field* field = G4TransportationManager::GetTransportationManager()
138                            ->GetFieldManager()
139                            ->GetDetectorField();
140 
141   G4Vector3D vectorU = GetVectorV().cross(GetVectorW());
142   G4double T1R       = 1. / (vTN * vectorU);
143 
144 #ifdef G4EVERBOSE
145   if(iverbose >= 4)
146     G4cout << "surf vectors:U,V,W " << vectorU << " " << GetVectorV() << " "
147            << GetVectorW() << "  T1R:" << T1R << G4endl;
148 #endif
149 
150   if(fCharge != 0 && field)
151   {
152     G4double pos[3];
153     pos[0] = fPosition.x() * cm;
154     pos[1] = fPosition.y() * cm;
155     pos[2] = fPosition.z() * cm;
156     G4double Hd[3];
157     field->GetFieldValue(pos, Hd);
158     G4ThreeVector H =
159       G4ThreeVector(Hd[0], Hd[1], Hd[2]) / tesla * 10.;  // in kilogauus
160     G4double magH  = H.mag();
161     G4double invP  = 1. / (fMomentum.mag() / GeV);
162     G4double magHM = magH * invP;
163     if(magH != 0.)
164     {
165       G4double magHM2 = fCharge / magH;
166       G4double Q      = -magHM * c_light / (km / ns);
167 #ifdef G4EVERBOSE
168       if(iverbose >= 4)
169         G4cout << GeV << " Q " << Q << " magHM " << magHM << " c_light/(km/ns) "
170                << c_light / (km / ns) << G4endl;
171 #endif
172 
173       G4double sinz = -H * vUN * magHM2;
174       G4double cosz = H * vVN * magHM2;
175       G4double T3R  = Q * std::pow(T1R, 3);
176       G4double UI   = vUN * vectorU;
177       G4double VI   = vVN * vectorU;
178 #ifdef G4EVERBOSE
179       if(iverbose >= 4)
180       {
181         G4cout << " T1R " << T1R << " T3R " << T3R << G4endl;
182         G4cout << " UI " << UI << " VI " << VI << " vectorU " << vectorU
183                << G4endl;
184         G4cout << " UJ " << UJ << " VJ " << VJ << G4endl;
185         G4cout << " UK " << UK << " VK " << VK << G4endl;
186       }
187 #endif
188 
189       transfM[1][3] = -UI * (VK * cosz - UK * sinz) * T3R;
190       transfM[1][4] = -VI * (VK * cosz - UK * sinz) * T3R;
191       transfM[2][3] = UI * (VJ * cosz - UJ * sinz) * T3R;
192       transfM[2][4] = VI * (VJ * cosz - UJ * sinz) * T3R;
193     }
194   }
195 
196   G4double T2R  = T1R * T1R;
197   transfM[0][0] = 1.;
198   transfM[1][1] = -UK * T2R;
199   transfM[1][2] = VK * cosLambda * T2R;
200   transfM[2][1] = UJ * T2R;
201   transfM[2][2] = -VJ * cosLambda * T2R;
202   transfM[3][3] = VK * T1R;
203   transfM[3][4] = -UK * T1R;
204   transfM[4][3] = -VJ * T1R;
205   transfM[4][4] = UJ * T1R;
206 
207 #ifdef G4EVERBOSE
208   if(iverbose >= 4)
209     G4cout << " SC2SD transf matrix A= " << transfM << G4endl;
210 #endif
211   fError = G4ErrorTrajErr(tpSC.GetError().similarity(transfM));
212 
213 #ifdef G4EVERBOSE
214   if(iverbose >= 1)
215     G4cout << "G4EP: error matrix SC2SD " << fError << G4endl;
216   if(iverbose >= 4)
217     G4cout << "G4ErrorSurfaceTrajState from SC " << *this << G4endl;
218 #endif
219 
220   return transfM;  // want to use trasnfM for the reverse transform
221 }
222 
223 //------------------------------------------------------------------------
224 void G4ErrorSurfaceTrajState::Init()
225 {
226   theTSType = G4eTS_OS;
227   BuildCharge();
228 }
229 
230 //------------------------------------------------------------------------
231 void G4ErrorSurfaceTrajState::Dump(std::ostream& out) const { out << *this; }
232 
233 //------------------------------------------------------------------------
234 std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts)
235 {
236   std::ios::fmtflags oldFlags = out.flags();
237   out.setf(std::ios::fixed, std::ios::floatfield);
238 
239   ts.DumpPosMomError(out);
240 
241   out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl;
242   out.flags(oldFlags);
243   return out;
244 }
245