eic-opticks/u4/ShimG4OpRayleigh.cc

0001 #include "G4SystemOfUnits.hh"
0002 #include "G4PhysicalConstants.hh"
0003 #include "ShimG4OpRayleigh.hh"
0004 #include "SLOG.hh"
0005 #include "ssys.h"
0006
0007 const plog::Severity ShimG4OpRayleigh::LEVEL = SLOG::EnvLevel("ShimG4OpRayleigh", "DEBUG");
0008
0009
0010 ShimG4OpRayleigh::ShimG4OpRayleigh()
0011     :
0012     G4OpRayleigh("OpRayleigh",fOptical)
0013 {
0014 }
0015
0016 ShimG4OpRayleigh::~ShimG4OpRayleigh()
0017 {
0018 }
0019
0020
0021
0022 #include "U4Stack.h"
0023 #include "U4UniformRand.h"
0024 #include "SEvt.hh"
0025
0026
0027 const bool ShimG4OpRayleigh::FLOAT = ssys::getenvbool("ShimG4OpRayleigh__FLOAT") ;
0028 const int  ShimG4OpRayleigh::PIDX  = ssys::getenvint("PIDX",-1);
0029 const bool ShimG4OpRayleigh::PIDX_ENABLED = ssys::getenvbool("ShimG4OpRayleigh__PIDX_ENABLED") ;
0030
0031 /**
0032 ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft
0033 -------------------------------------------------------
0034
0035 Shim makes process classname appear in SBacktrace.h enabling U4Random::flat/U4Stack::Classify
0036
0037
0038 **/
0039
0040 // void ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft(){ G4VProcess::ResetNumberOfInteractionLengthLeft(); }
0041  void ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft()
0042 {
0043     G4double u = G4UniformRand() ;
0044
0045     LOG(LEVEL)
0046         << U4UniformRand::Desc(u, SEvt::UU )
0047         ;
0048
0049 #ifndef PRODUCTION
0050 #ifdef DEBUG_TAG
0051     SEvt::AddTag( 1, U4Stack_RayleighDiscreteReset, u );
0052 #endif
0053 #endif
0054
0055
0056     if(FLOAT)
0057     {
0058         float f =  -1.f*std::log( float(u) ) ;
0059         theNumberOfInteractionLengthLeft = f ;
0060     }
0061     else
0062     {
0063         theNumberOfInteractionLengthLeft = -1.*G4Log(u) ;
0064     }
0065     theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
0066 }
0067
0068
0069  G4double ShimG4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,  G4double, G4ForceCondition* )
0070 {
0071      G4double len = G4OpRayleigh::GetMeanFreePath( aTrack, 0, 0 );
0072      //std::cout << "ShimG4OpRayleigh::GetMeanFreePath " << len << std::endl ;
0073      return len ;
0074 }
0075
0076  G4double ShimG4OpRayleigh::PostStepGetPhysicalInteractionLength(const G4Track& track, G4double   previousStepSize, G4ForceCondition* condition)
0077 {
0078   if ( (previousStepSize < 0.0) || (theNumberOfInteractionLengthLeft<=0.0)) {
0079     // beggining of tracking (or just after DoIt of this process)
0080     ResetNumberOfInteractionLengthLeft();
0081   } else if ( previousStepSize > 0.0) {
0082     // subtract NumberOfInteractionLengthLeft
0083     SubtractNumberOfInteractionLengthLeft(previousStepSize);
0084   } else {
0085     // zero step
0086     //  DO NOTHING
0087   }
0088
0089   // condition is set to "Not Forced"
0090   *condition = NotForced;
0091
0092   // get mean free path
0093   currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition);
0094
0095   G4double value;
0096   if (currentInteractionLength <DBL_MAX) {
0097
0098      if( FLOAT )
0099      {
0100           float fvalue = float(theNumberOfInteractionLengthLeft) * float(currentInteractionLength) ;
0101           value = fvalue ;
0102      }
0103      else
0104      {
0105           value = theNumberOfInteractionLengthLeft * currentInteractionLength ;
0106      }
0107
0108
0109       if(track.GetTrackID() - 1  == PIDX && PIDX_ENABLED)
0110       {
0111            std::cout
0112                << "ShimG4OpRayleigh::PostStepGetPhysicalInteractionLength"
0113                << " PIDX " << PIDX
0114                << " currentInteractionLength " << std::setw(10) << std::fixed << std::setprecision(7) << currentInteractionLength
0115                << " theNumberOfInteractionLengthLeft " << std::setw(10) << std::fixed << std::setprecision(7) << theNumberOfInteractionLengthLeft
0116                << " value " << std::setw(10) << std::fixed << std::setprecision(7) << value
0117                << std::endl
0118                ;
0119       }
0120
0121
0122   } else {
0123     value = DBL_MAX;
0124   }
0125 #ifdef G4VERBOSE
0126   if (verboseLevel>1){
0127     G4cout << "G4VDiscreteProcess::PostStepGetPhysicalInteractionLength ";
0128     G4cout << "[ " << GetProcessName() << "]" <<G4endl;
0129     track.GetDynamicParticle()->DumpInfo();
0130     G4cout << " in Material  " <<  track.GetMaterial()->GetName() <<G4endl;
0131     G4cout << "InteractionLength= " << value/cm <<"[cm] " <<G4endl;
0132   }
0133 #endif
0134   return value;
0135
0136 }
0137
0138
0139
0140
0141 G4VParticleChange* ShimG4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
0142 {
0143         aParticleChange.Initialize(aTrack);
0144
0145         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0146
0147         if (verboseLevel>0) {
0148                 G4cout << "Scattering Photon!" << G4endl;
0149                 G4cout << "Old Momentum Direction: "
0150                        << aParticle->GetMomentumDirection() << G4endl;
0151                 G4cout << "Old Polarization: "
0152                        << aParticle->GetPolarization() << G4endl;
0153         }
0154
0155         G4double cosTheta;
0156         G4ThreeVector OldMomentumDirection, NewMomentumDirection;
0157         G4ThreeVector OldPolarization, NewPolarization;
0158
0159         G4double rand, constant;
0160         G4double CosTheta, SinTheta, SinPhi, CosPhi, unit_x, unit_y, unit_z;
0161
0162         G4double u ;
0163         G4double u_loopexit ;
0164
0165
0166         do {
0167            // Try to simulate the scattered photon momentum direction
0168            // w.r.t. the initial photon momentum direction
0169
0170            CosTheta = G4UniformRand();
0171
0172 #ifndef PRODUCTION
0173 #ifdef DEBUG_TAG
0174            SEvt::AddTag( 1, U4Stack_RayleighScatter, CosTheta );   // 0
0175 #endif
0176 #endif
0177
0178            SinTheta = std::sqrt(1.-CosTheta*CosTheta);
0179            // consider for the angle 90-180 degrees
0180
0181            u = G4UniformRand() ;
0182
0183 #ifndef PRODUCTION
0184 #ifdef DEBUG_TAG
0185            SEvt::AddTag( 1, U4Stack_RayleighScatter, u );      // 1
0186 #endif
0187 #endif
0188
0189            if (u < 0.5) CosTheta = -CosTheta;
0190
0191            // simulate the phi angle
0192
0193            u = G4UniformRand() ;
0194
0195 #ifndef PRODUCTION
0196 #ifdef DEBUG_TAG
0197            SEvt::AddTag( 1, U4Stack_RayleighScatter, u );    // 2
0198 #endif
0199 #endif
0200
0201            rand = twopi*u;
0202            SinPhi = std::sin(rand);
0203            CosPhi = std::cos(rand);
0204
0205            // start constructing the new momentum direction
0206        unit_x = SinTheta * CosPhi;
0207        unit_y = SinTheta * SinPhi;
0208        unit_z = CosTheta;
0209        NewMomentumDirection.set (unit_x,unit_y,unit_z);
0210
0211            // Rotate the new momentum direction into global reference system
0212            OldMomentumDirection = aParticle->GetMomentumDirection();
0213            OldMomentumDirection = OldMomentumDirection.unit();
0214            NewMomentumDirection.rotateUz(OldMomentumDirection);
0215            NewMomentumDirection = NewMomentumDirection.unit();
0216
0217            // calculate the new polarization direction
0218            // The new polarization needs to be in the same plane as the new
0219            // momentum direction and the old polarization direction
0220            OldPolarization = aParticle->GetPolarization();
0221            constant = -NewMomentumDirection.dot(OldPolarization);
0222
0223            NewPolarization = OldPolarization + constant*NewMomentumDirection;
0224            NewPolarization = NewPolarization.unit();
0225
0226            // There is a corner case, where the Newmomentum direction
0227            // is the same as oldpolariztion direction:
0228            // random generate the azimuthal angle w.r.t. Newmomentum direction
0229
0230            u = G4UniformRand() ;
0231
0232 #ifndef PRODUCTION
0233 #ifdef DEBUG_TAG
0234            SEvt::AddTag( 1, U4Stack_RayleighScatter, u );   // 3
0235 #endif
0236 #endif
0237
0238            if (NewPolarization.mag() == 0.) {
0239               rand = u*twopi;
0240               NewPolarization.set(std::cos(rand),std::sin(rand),0.);
0241               NewPolarization.rotateUz(NewMomentumDirection);
0242            } else {
0243               // There are two directions which are perpendicular
0244               // to the new momentum direction
0245               if (u < 0.5) NewPolarization = -NewPolarization;
0246            }
0247
0248        // simulate according to the distribution cos^2(theta)
0249            cosTheta = NewPolarization.dot(OldPolarization);
0250
0251            u_loopexit = G4UniformRand() ;
0252
0253 #ifndef PRODUCTION
0254 #ifdef DEBUG_TAG
0255            SEvt::AddTag( 1, U4Stack_RayleighScatter, u_loopexit );   // 4
0256 #endif
0257 #endif
0258
0259           // Loop checking, 13-Aug-2015, Peter Gumplinger
0260         } while (std::pow(cosTheta,2) < u_loopexit );
0261
0262         aParticleChange.ProposePolarization(NewPolarization);
0263         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
0264
0265         if (verboseLevel>0) {
0266                 G4cout << "New Polarization: "
0267                      << NewPolarization << G4endl;
0268                 G4cout << "Polarization Change: "
0269                      << *(aParticleChange.GetPolarization()) << G4endl;
0270                 G4cout << "New Momentum Direction: "
0271                      << NewMomentumDirection << G4endl;
0272                 G4cout << "Momentum Change: "
0273                      << *(aParticleChange.GetMomentumDirection()) << G4endl;
0274         }
0275
0276         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0277 }
0278
0279
0280
0281