00001
00002 #ifndef BDSContinuousSR_h
00003 #define BDSContinuousSR_h 1
00004
00005 #include "BDSGlobalConstants.hh"
00006 #include "BDSMaterials.hh"
00007
00008 #include "G4ios.hh"
00009 #include "globals.hh"
00010 #include "G4VDiscreteProcess.hh"
00011 #include "G4Track.hh"
00012 #include "G4Step.hh"
00013 #include "G4Electron.hh"
00014 #include "G4Positron.hh"
00015 #include "G4Field.hh"
00016 #include "G4FieldManager.hh"
00017 #include "G4AffineTransform.hh"
00018 #include "G4NavigationHistory.hh"
00019
00020 class BDSContinuousSR : public G4VDiscreteProcess
00021 {
00022 public:
00023
00024 BDSContinuousSR(const G4String& processName = "contSR");
00025
00026 ~BDSContinuousSR();
00027
00028 G4bool IsApplicable(const G4ParticleDefinition&);
00029
00030 G4double GetMeanFreePath(const G4Track& track,
00031 G4double previousStepSize,
00032 G4ForceCondition* condition );
00033
00034 G4VParticleChange *PostStepDoIt(const G4Track& track,
00035 const G4Step& step);
00036
00037 G4double SynGenC(G4double xmin);
00038 G4double SynRadC(G4double x);
00039
00040 protected:
00041
00042 private:
00043
00044
00045 BDSContinuousSR & operator=(const BDSContinuousSR &right);
00046 BDSContinuousSR(const BDSContinuousSR&);
00047
00048 G4double nExpConst;
00049 G4double CritEngFac;
00050
00051 private:
00052 };
00053
00054 inline G4bool
00055 BDSContinuousSR::IsApplicable(const G4ParticleDefinition& particle)
00056 {
00057 return( (&particle == G4Electron::Electron())
00058 ||(&particle == G4Positron::Positron()) );
00059 }
00060
00061 inline G4double
00062 BDSContinuousSR::GetMeanFreePath(const G4Track& track,
00063 G4double,
00064 G4ForceCondition* ForceCondition)
00065 {
00066 *ForceCondition = NotForced ;
00067
00068 G4double rfact = 1.;
00069
00070
00071
00072
00073 G4double MeanFreePath;
00074 G4FieldManager* TheFieldManager=
00075 track.GetVolume()->GetLogicalVolume()->GetFieldManager();
00076
00077 if(track.GetTotalEnergy()<BDSGlobalConstants::Instance()->GetThresholdCutCharged())
00078 return DBL_MAX;
00079
00080
00081
00082
00083
00084 if(TheFieldManager)
00085 {
00086 const G4Field* pField = TheFieldManager->GetDetectorField() ;
00087 G4ThreeVector globPosition = track.GetPosition() ;
00088 G4double globPosVec[3], FieldValueVec[3]={0.} ;
00089 globPosVec[0] = globPosition.x() ;
00090 globPosVec[1] = globPosition.y() ;
00091 globPosVec[2] = globPosition.z() ;
00092
00093 if(pField)
00094 pField->GetFieldValue( globPosVec, FieldValueVec ) ;
00095
00096 G4double Blocal= FieldValueVec[1];
00097 if ( FieldValueVec[0]!=0.)
00098 Blocal=sqrt(Blocal*Blocal+FieldValueVec[0]*FieldValueVec[0]);
00099
00100
00101
00102 if(track.GetMaterial()==BDSMaterials::Instance()->GetMaterial(BDSGlobalConstants::Instance()->GetVacuumMaterial())
00103 && Blocal !=0 )
00104 {
00105 G4ThreeVector InitMag=track.GetMomentum();
00106
00107
00108
00109 G4AffineTransform tf(track.GetTouchable()->GetHistory()->GetTopTransform().Inverse());
00110 const G4RotationMatrix Rot=tf.NetRotation();
00111 const G4ThreeVector Trans=-tf.NetTranslation();
00112 InitMag=Rot*InitMag;
00113
00114
00115 G4double Rlocal=(InitMag.z()/CLHEP::GeV)/(0.299792458 * Blocal/CLHEP::tesla) *CLHEP::m;
00116
00117 MeanFreePath=
00118 fabs(Rlocal)/(track.GetTotalEnergy()*nExpConst);
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 return rfact * MeanFreePath;
00129 }
00130 else
00131 return DBL_MAX;
00132 }
00133 else
00134 return DBL_MAX;
00135
00136 }
00137
00138
00139
00140 #endif