source: trunk/source/error_propagation/src/G4ErrorFreeTrajState.cc @ 1254

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

update geant4.9.3 tag
File size: 30.9 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// $Id: G4ErrorFreeTrajState.cc,v 1.8 2009/05/14 13:53:06 arce Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// ------------------------------------------------------------
30//      GEANT 4 class implementation file
31// ------------------------------------------------------------
32//
33#include "G4ErrorFreeTrajState.hh"
34#include "G4ErrorFreeTrajParam.hh"
35#include "G4ErrorSurfaceTrajState.hh"
36
37#include "G4ErrorMatrix.hh"
38#include <iomanip>
39
40#include "G4Field.hh"
41#include "G4FieldManager.hh"
42#include "G4TransportationManager.hh"
43#include "G4GeometryTolerance.hh"
44#include "G4Material.hh"
45#include "G4ErrorPropagatorData.hh"
46
47//------------------------------------------------------------------------
48G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partType, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partType, pos, mom, errmat )
49{
50  fTrajParam = G4ErrorFreeTrajParam( pos, mom );
51  Init();
52}
53
54
55//------------------------------------------------------------------------
56G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() )
57{
58  //  G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
59  // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane 
60  //  G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
61  //  G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
62
63  fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum );
64  Init();
65
66  //----- Get the error matrix in SC coordinates
67  G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
68  G4double mom = fMomentum.mag();
69  G4double mom2 = fMomentum.mag2();
70  G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPV()*tpSDparam.GetPV()) );
71  G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 );
72  G4Vector3D vectorU = tpSDparam.GetVectorV().cross(  tpSDparam.GetVectorW() );
73  G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW();
74
75#ifdef G4EVERBOSE
76   if( iverbose >= 5){
77     G4double pc2 = std::asin( vTN.z() );
78     G4double pc3 = std::atan (vTN.y()/vTN.x());
79 
80     G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() <<  " diff " << pc2-GetParameters().GetLambda() << G4endl;
81     G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() <<  " diff " << pc3-GetParameters().GetPhi() << G4endl;
82   }
83#endif
84
85  //--- Get the unit vectors perp to P
86  G4double cosl = std::cos( GetParameters().GetLambda() ); 
87  if (cosl < 1.E-30) cosl = 1.E-30;
88  G4double cosl1 = 1./cosl;
89  G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. );
90  G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl );
91
92  G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.);
93  G4Vector3D vVperp = vUperp.cross( fMomentum );
94  vUperp *= 1./vUperp.mag();
95  vVperp *= 1./vVperp.mag();
96
97#ifdef G4EVERBOSE
98   if( iverbose >= 5){
99     G4cout << " CHECK: vUN " << vUN << " = " << vUperp <<  " diff " << (vUN-vUperp).mag() << G4endl;
100     G4cout << " CHECK: vVN " << vVN << " = " << vVperp <<  " diff " << (vVN-vVperp).mag() << G4endl;
101   }
102#endif
103
104  //get the dot products of vectors perpendicular to direction and vector defining SD plane
105  G4double dUU = vUperp * tpSD.GetVectorV();
106  G4double dUV = vUperp * tpSD.GetVectorW();
107  G4double dVU = vVperp * tpSD.GetVectorV();
108  G4double dVV = vVperp * tpSD.GetVectorW();
109
110
111  //--- Get transformation first
112  G4ErrorMatrix transfM(5, 5, 1 );
113  //--- Get magnetic field
114  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
115  G4ThreeVector dir = fTrajParam.GetDirection();
116  G4double invCosTheta = 1./std::cos( dir.theta() );
117
118  if( fCharge != 0 
119&& field ) {
120    G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm;
121    G4double h1[3];
122    field->GetFieldValue( pos1, h1 );
123    G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.;
124    G4double magHPre = HPre.mag();
125    G4double invP = 1./fMomentum.mag();
126    G4double magHPreM = magHPre * invP;
127    if( magHPre != 0. ) {
128      G4double magHPreM2 = fCharge / magHPre;
129
130      G4double Q = -magHPreM * c_light;
131      G4double sinz = -HPre*vUperp * magHPreM2;
132      G4double cosz =  HPre*vVperp * magHPreM2;
133
134      transfM[1][3] = -Q*dir.y()*sinz;
135      transfM[1][4] = -Q*dir.z()*sinz;
136      transfM[2][3] = -Q*dir.y()*cosz*invCosTheta;
137      transfM[2][4] = -Q*dir.z()*cosz*invCosTheta;
138    }
139  }
140
141  transfM[0][0] = 1.;
142  transfM[1][1] = dir.x()*dVU;
143  transfM[1][2] = dir.x()*dVV;
144  transfM[2][1] = dir.x()*dUU*invCosTheta;
145  transfM[2][2] = dir.x()*dUV*invCosTheta;
146  transfM[3][3] = dUU;
147  transfM[3][4] = dUV;
148  transfM[4][3] = dVU;
149  transfM[4][4] = dVV;
150
151  fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) );
152
153#ifdef G4EVERBOSE
154  if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl;
155  if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
156#endif
157}
158
159
160//------------------------------------------------------------------------
161void G4ErrorFreeTrajState::Init()
162{
163  theTSType = G4eTS_FREE;
164  BuildCharge();
165  theTransfMat = G4ErrorMatrix(5,5,0);
166  //-  theFirstStep = true;
167}
168
169//------------------------------------------------------------------------
170void G4ErrorFreeTrajState::Dump( std::ostream& out ) const
171{
172  out << *this;
173}
174
175//------------------------------------------------------------------------
176G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack )
177{
178  G4int ierr = 0;
179  fTrajParam.Update( aTrack );
180  UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() );
181  return ierr;
182
183}
184
185
186//------------------------------------------------------------------------
187std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
188{
189  out.setf(std::ios::fixed,std::ios::floatfield);
190
191 
192  ts.DumpPosMomError( out );
193 
194  out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
195
196  return out;
197
198}
199
200
201//------------------------------------------------------------------------
202G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack )
203{
204  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
205  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.;
206
207  G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
208
209  if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0;
210 
211#ifdef G4EVERBOSE
212  if( iverbose >= 2 )G4cout << "  G4ErrorFreeTrajState::PropagateError " << G4endl;
213#endif
214
215  // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
216  G4Point3D vposPost = aTrack->GetPosition()/cm;
217  G4Vector3D vpPost = aTrack->GetMomentum()/GeV;
218  //  G4Point3D vposPre = fPosition/cm;
219  //  G4Vector3D vpPre = fMomentum/GeV;
220  G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm;
221  G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV;
222  //correct to avoid propagation along Z
223  if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV );
224  if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV );
225
226  G4double pPre = vpPre.mag();
227  G4double pPost = vpPost.mag();
228#ifdef G4EVERBOSE
229  if( iverbose >= 2 ) {
230    G4cout << "G4EP: vposPre " << vposPre << G4endl
231              << "G4EP: vposPost " << vposPost << G4endl;
232    G4cout << "G4EP: vpPre " << vpPre << G4endl
233              << "G4EP: vpPost " << vpPost << G4endl;
234    G4cout << " err start step " << fError << G4endl;
235    G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
236  }
237#endif
238
239  if( pPre == 0. || pPost == 0 ) return 2;
240  G4double pInvPre = 1./pPre;
241  G4double pInvPost = 1./pPost;
242  G4double deltaPInv = pInvPost - pInvPre;
243
244  G4Vector3D vpPreNorm = vpPre * pInvPre;
245  G4Vector3D vpPostNorm = vpPost * pInvPost;
246  //  if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl;
247  //return if propagation along Z?? 
248  if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4;
249  G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre
250  G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost
251  G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() );
252
253#ifdef G4EVERBOSE
254  if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
255#endif
256  //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
257  //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
258  G4ErrorMatrix transf(5, 5, 0 );
259
260  transf[3][2] = stepLengthCm * sinpPost;
261  transf[4][1] = stepLengthCm;
262  for( size_t ii=0;ii < 5; ii++ ){
263    transf[ii][ii] = 1.;
264  }
265#ifdef G4EVERBOSE
266  if( iverbose >= 2 ) {
267    G4cout << "G4EP: transf matrix neutral " << transf;
268  }
269#endif
270
271  //  charge X propagation direction
272  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
273  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) {
274    charge *= -1.; 
275  }
276  //  G4cout << " charge " << charge << G4endl;
277  //t check if particle has charge
278  //t  if( charge == 0 ) goto 45;
279  // check if the magnetic field is = 0.
280
281  //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed)
282  G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm;
283  G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm;
284  G4double h1[3], h2[3];
285
286  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
287  if( !field ) return 0; //goto 45
288 
289  // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
290  if( charge != 0. && field ) {
291
292    field->GetFieldValue( pos1, h1 );
293    field->GetFieldValue( pos2, h2 );
294    G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss)
295    G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.;
296    G4double magHPre = HPre.mag();
297    G4double magHPost = HPost.mag();
298#ifdef G4EVERBOSE
299    if( iverbose >= 2 ) G4cout << "G4EP: HPre " << HPre << G4endl
300                            << "G4EP: HPost " << HPost << G4endl;
301#endif
302   
303  if( magHPre + magHPost != 0. ) {
304     
305   //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
306    G4double gam;
307    if( magHPost != 0. ){ 
308      gam = HPost * vpPostNorm / magHPost;
309    }else {
310      gam = HPre * vpPreNorm / magHPre;
311    }
312   
313    // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory
314    G4double alphaSqr = 1. - gam * gam;
315    G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2();
316    G4double delhp6Sqr = 300.*300.; 
317#ifdef G4EVERBOSE
318    if( iverbose >= 2 ) G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr << " diffHSqr " << diffHSqr << G4endl;
319#endif
320    if( diffHSqr * alphaSqr > delhp6Sqr ) return 3;
321
322
323    //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
324    G4double pInvAver = 1./(pInvPre + pInvPost );
325    G4double CFACT8 = 2.997925E-4; 
326    //G4double HAver
327    G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 );
328    G4double HAver = vHAverNorm.mag();
329    G4double invHAver = 1./HAver;
330    vHAverNorm *= invHAver;
331#ifdef G4EVERBOSE
332    if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl;
333#endif
334
335    G4double pAver = (pPre+pPost)*0.5;
336    G4double QAver = -HAver/pAver;
337    G4double thetaAver = QAver * stepLengthCm;
338    G4double sinThetaAver = std::sin(thetaAver);
339    G4double cosThetaAver = std::cos(thetaAver);
340    G4double gamma = vHAverNorm * vpPostNorm;
341    G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm );
342   
343#ifdef G4EVERBOSE
344    if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << G4endl;
345#endif
346    G4double AU = 1./vpPreNorm.perp();
347    //t  G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
348    G4ThreeVector vUPre( -AU*vpPreNorm.y(), 
349                      AU*vpPreNorm.x(), 
350                      0. );
351    G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(), 
352                      vpPreNorm.z()*vUPre.x(), 
353                      vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() );
354   
355    //
356    AU = 1./vpPostNorm.perp();
357    //t  G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
358    G4ThreeVector vUPost( -AU*vpPostNorm.y(), 
359                       AU*vpPostNorm.x(), 
360                       0. );
361    G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(), 
362                       vpPostNorm.z()*vUPost.x(), 
363                       vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() );
364#ifdef G4EVERBOSE
365    //-    G4cout << " vpPostNorm " << vpPostNorm << G4endl;
366    if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
367#endif
368    G4Point3D deltaPos( vposPre - vposPost );
369
370    // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
371    // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
372    // * *** TAKEN INTO ACCOUNT
373   
374    G4double QP = QAver * pAver; // = -HAver
375#ifdef G4EVERBOSE
376    if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl;
377#endif
378    G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() );
379    G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() );
380    G4double OMcosThetaAver = 1. - cosThetaAver;
381#ifdef G4EVERBOSE
382    if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl;
383#endif
384    G4double TMSINT = thetaAver - sinThetaAver;
385#ifdef G4EVERBOSE
386    if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
387#endif
388   
389    G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(),
390                          vHAverNorm.z() * vUPre.x(),
391                          vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() );
392#ifdef G4EVERBOSE
393    //    if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl;
394#endif
395    G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
396                          vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
397                          vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() );
398#ifdef G4EVERBOSE
399    if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
400#endif
401   
402    //------------------- COMPUTE MATRIX
403    //---------- 1/P
404   
405    transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm)
406      +2.*deltaPInv*pAver;
407   
408    transf[0][1] =  -deltaPInv/thetaAver*
409      ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) +
410        sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
411        OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) );
412   
413    transf[0][2] =  -sinpPre*deltaPInv/thetaAver*
414      ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            ) +
415        sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            ) +
416        OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) );
417   
418    transf[0][3] =  -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            );
419   
420    transf[0][4] =  -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
421   
422    // ***   Lambda
423    transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())
424      *(1.+deltaPInv*pAver);
425#ifdef G4EVERBOSE
426     if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0]  << " " <<  -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z()
427      << " " << deltaPInv<< " " << pAver << G4endl;
428#endif
429   
430    transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
431      sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
432      OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
433      (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
434      ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
435            OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
436            TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
437   
438    transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y()            ) +
439      sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
440      OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            )*
441      (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
442      ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            ) +
443            OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y()             ) -
444            TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            ) );
445    transf[1][2] = sinpPre*transf[1][2];
446   
447    transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            );
448   
449    transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
450   
451    // ***   Phi
452   
453    transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv
454      *(1.+deltaPInv*pAver);
455#ifdef G4EVERBOSE
456   if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv
457         <<" "<<deltaPInv<<" "<<pAver<< G4endl;
458#endif
459    transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y()            ) +
460      sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y()             ) +
461      OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
462      (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y()            ) +
463      ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
464            OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
465            TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
466    transf[2][1] = sinpPostInv*transf[2][1];
467   
468    transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y()            ) +
469      sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y()             ) +
470      OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            )*
471      (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y()            ) +
472      ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            ) +
473            OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y()             ) -
474            TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            ) );
475    transf[2][2] = sinpPostInv*sinpPre*transf[2][2];
476   
477    transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y()            )*sinpPostInv;
478#ifdef G4EVERBOSE
479    if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl;
480#endif
481   
482    transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv;
483   
484    // ***   Yt
485   
486    transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() )
487      *(1.+deltaPInv*pAver);
488#ifdef G4EVERBOSE
489   if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y() 
490      <<" "<<deltaPInv<<" "<<pAver<<G4endl;
491#endif
492
493    transf[3][1] = (   sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y()            ) +
494                       OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y()             ) +
495                       TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y()            )*
496                       (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
497   
498    transf[3][2] = (   sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y()            ) +
499                       OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y()             ) +
500                       TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y()            )*
501                       (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            ) )*sinpPre/QAver;
502#ifdef G4EVERBOSE
503   if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<<
504                       OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<<
505                       TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<<
506      vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl;
507#endif
508   
509    transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y()            );
510   
511    transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y()            );
512
513    // ***   Zt
514    transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z())
515      *(1.+deltaPInv*pAver);
516   
517    transf[4][1] = (   sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
518                       OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
519                       TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
520                       (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
521#ifdef G4EVERBOSE
522    if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl;
523#endif
524   
525    transf[4][2] = (   sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y()            ) +
526                       OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
527                       TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
528                       (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y()            ) )*sinpPre/QAver;
529
530    transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y()  );
531
532    transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()); 
533    //   if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x()  <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
534
535 
536#ifdef G4EVERBOSE
537    if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl;
538#endif
539    /*    for( G4int ii=0;ii<5;ii++){
540      for( G4int jj=0;jj<5;jj++){
541        G4cout << transf[ii][jj] << " ";
542      }
543      G4cout << G4endl;
544      } */
545   }
546  }
547  // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
548  /*  if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl;
549  if( theFirstStep ) {
550    theTransfMat = transf;
551    theFirstStep = false;
552  }else{
553    theTransfMat = theTransfMat * transf;
554    if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl;
555  }
556  */
557    theTransfMat = transf;
558#ifdef G4EVERBOSE
559    if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
560    if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl
561                                  << " transf matrix " << theTransfMat.T() << G4endl;
562#endif
563   
564    fError = fError.similarity(theTransfMat).T();
565    //-    fError = transf * fError * transf.T();
566#ifdef G4EVERBOSE
567    if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl;
568#endif
569
570    //? S = B*S*BT S.similarity(B)
571    //? R = S
572    // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
573   
574    PropagateErrorMSC( aTrack );
575   
576    PropagateErrorIoni( aTrack );
577   
578    return 0;
579}
580
581
582//------------------------------------------------------------------------
583G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack )
584{ 
585  G4ThreeVector vpPre = aTrack->GetMomentum()/GeV;
586  G4double pPre = vpPre.mag();
587  G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV);
588  G4double  stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
589
590  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
591  G4double effZ, effA;
592  CalculateEffectiveZandA( mate, effZ, effA );
593
594#ifdef G4EVERBOSE
595  if( iverbose >= 4 ) G4cout << "material " << mate->GetName() 
596                     //<< " " << mate->GetZ() << " "  << mate->GetA()
597                        << " " << effZ << " " << effA
598                        << " "  << mate->GetDensity()/g*mole << " " << mate->GetRadlen()/cm << " " << mate->GetNuclearInterLength()/cm << G4endl;
599#endif
600
601  G4double RI = stepLengthCm / (mate->GetRadlen()/cm);
602#ifdef G4EVERBOSE
603  if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI " << RI << " stepLengthCm " << stepLengthCm << " radlen " << (mate->GetRadlen()/cm) << " " << RI*1.e10 << G4endl;
604#endif
605  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
606  G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta );
607#ifdef G4EVERBOSE
608  if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl;
609#endif
610  G4double S1 = DD*stepLengthCm*stepLengthCm/3.;
611  G4double S2 = DD;
612  G4double S3 = DD*stepLengthCm/2.;
613
614  G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre;
615#ifdef G4EVERBOSE
616  if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 "  << S2 << " S3 "  << S3 << " CLA " << CLA << G4endl;
617#endif
618  fError[1][1] += S2;
619  fError[1][4] -= S3;
620  fError[2][2] += S2/CLA/CLA;
621  fError[2][3] += S3/CLA;
622  fError[3][3] += S1;
623  fError[4][4] += S1;
624
625#ifdef G4EVERBOSE
626  if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
627#endif
628
629  return 0;
630}
631
632
633//------------------------------------------------------------------------
634void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA )
635{ 
636  effZ = 0.;
637  effA = 0.;
638  G4int ii, nelem = mate->GetNumberOfElements();
639  const G4double* fracVec = mate->GetFractionVector();
640  for(ii=0; ii < nelem; ii++ ) {
641    effZ += mate->GetElement( ii )->GetZ() * fracVec[ii];
642    effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole;
643  }
644
645}
646
647
648//------------------------------------------------------------------------
649G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack )
650{ 
651  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
652  G4double DEDX2;
653  if( stepLengthCm < 1.E-7 ) {
654    DEDX2=0.;
655  }
656  //  *     Calculate xi factor (KeV).
657  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
658  G4double effZ, effA;
659  CalculateEffectiveZandA( mate, effZ, effA );
660
661  G4double Etot = aTrack->GetTotalEnergy()/GeV;
662  G4double beta = aTrack->GetMomentum().mag()/GeV / Etot;
663  G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
664  G4double gamma = Etot / mass;
665 
666  // *     Calculate xi factor (KeV).
667  G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) / 
668    (effA*beta*beta);
669
670#ifdef G4EVERBOSE
671  if( iverbose >= 2 ){
672    G4cout << "G4EP:IONI: XI " << XI << " beta " << beta << " gamma " << gamma << G4endl;
673    G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl;
674  }
675#endif
676  // *     Maximum energy transfer to atomic electron (KeV).
677  G4double eta = beta*gamma;
678  G4double etasq = eta*eta;
679  G4double eMass = 0.51099906/GeV;
680  G4double massRatio = eMass / mass;
681  G4double F1 = 2*eMass*etasq;
682  G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
683  G4double Emax = 1.E+6*F1/F2;
684
685  //  * *** and now sigma**2  in GeV
686  G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12;
687#ifdef G4EVERBOSE
688  if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX2 " << dedxSq << " emass " << eMass << " Emax " << Emax << G4endl;
689#endif
690
691  //  if( iverbose >= 2 ) G4cout << "G4EP:IONI: Etot " << Etot << " DEDX2 " << dedxSq << " emass " << eMass << G4endl;
692 
693  G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag();
694  pPre6 = std::pow(pPre6, 6 );
695  //Apply it to error
696  fError[0][0] += Etot*Etot*dedxSq / pPre6;
697#ifdef G4EVERBOSE
698  if( iverbose >= 2 ) G4cout << "G4:IONI getot " << Etot << " dedx2 " << dedxSq << " p " << pPre6 << G4endl;
699  if( iverbose >= 2 ) G4cout << "G4EP:IONI: error_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl;
700#endif
701
702  return 0;
703}
704
Note: See TracBrowser for help on using the repository browser.