In This Package:

DsG4OpRayleigh Class Reference

A slightly modified version of G4OpRayleigh. More...

#include <DsG4OpRayleigh.h>

List of all members.

Public Member Functions
	DsG4OpRayleigh (const G4String &processName="OpRayleigh", G4ProcessType type=fOptical)
	~DsG4OpRayleigh ()
G4bool	IsApplicable (const G4ParticleDefinition &aParticleType)
G4double	GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
G4PhysicsTable *	GetPhysicsTable () const
void	DumpPhysicsTable () const
Protected Attributes
G4PhysicsTable *	thePhysicsTable
Private Member Functions
void	BuildThePhysicsTable ()
G4PhysicsOrderedFreeVector *	RayleighAttenuationLengthGenerator (G4MaterialPropertiesTable *aMPT)
Private Attributes
G4bool	DefaultWater

---

Detailed Description

A slightly modified version of G4OpRayleigh.

It is modified to make the Rayleigh Scattering happen with different waters defined in /dd/Material/

This was taken from G4.9.1p1

zhangh@bnl.gov on 8th, July

Definition at line 88 of file DsG4OpRayleigh.h.

---

Constructor & Destructor Documentation

DsG4OpRayleigh::DsG4OpRayleigh	(	const G4String &	processName = "OpRayleigh",
		G4ProcessType	type = fOptical
	)

Definition at line 91 of file DsG4OpRayleigh.cc.

           : G4VDiscreteProcess(processName, type)
{
        SetProcessSubType(fOpRayleigh);

        thePhysicsTable = 0;

        DefaultWater = false;

        if (verboseLevel>0) {
           G4cout << GetProcessName() << " is created " << G4endl;
        }

        BuildThePhysicsTable();
}

DsG4OpRayleigh::~DsG4OpRayleigh

(

)

Definition at line 115 of file DsG4OpRayleigh.cc.

{
        if (thePhysicsTable!= 0) {
           thePhysicsTable->clearAndDestroy();
           delete thePhysicsTable;
        }
}

---

Member Function Documentation

G4bool DsG4OpRayleigh::IsApplicable

(

const G4ParticleDefinition &

aParticleType

)

[inline]

Definition at line 170 of file DsG4OpRayleigh.h.

{
  return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
}

G4double DsG4OpRayleigh::GetMeanFreePath	(	const G4Track &	aTrack,
		G4double	,
		G4ForceCondition *
	)

Definition at line 261 of file DsG4OpRayleigh.cc.

{
        const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
        const G4Material* aMaterial = aTrack.GetMaterial();

        G4double thePhotonEnergy = aParticle->GetTotalEnergy();

        G4double AttenuationLength = DBL_MAX;

        if ((strcmp(aMaterial->GetName(), "/dd/Materials/Water") == 0 ||
             strcmp(aMaterial->GetName(), "/dd/Materials/OwsWater") == 0 ||
             strcmp(aMaterial->GetName(), "/dd/Materials/IwsWater") == 0 )
            && DefaultWater){

           G4bool isOutRange;

           AttenuationLength =
                (*thePhysicsTable)(aMaterial->GetIndex())->
                           GetValue(thePhotonEnergy, isOutRange);
        }
        else {

           G4MaterialPropertiesTable* aMaterialPropertyTable =
                           aMaterial->GetMaterialPropertiesTable();

           if(aMaterialPropertyTable){
             G4MaterialPropertyVector* AttenuationLengthVector =
                   aMaterialPropertyTable->GetProperty("RAYLEIGH");
             if(AttenuationLengthVector){
               AttenuationLength = AttenuationLengthVector ->
                                    GetProperty(thePhotonEnergy);
             }
             else{
//               G4cout << "No Rayleigh scattering length specified" << G4endl;
             }
           }
           else{
//             G4cout << "No Rayleigh scattering length specified" << G4endl; 
           }
        }

        return AttenuationLength;
}

G4VParticleChange * DsG4OpRayleigh::PostStepDoIt	(	const G4Track &	aTrack,
		const G4Step &	aStep
	)

Definition at line 131 of file DsG4OpRayleigh.cc.

{
        aParticleChange.Initialize(aTrack);

        const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();

        if (verboseLevel>0) {
                G4cout << "Scattering Photon!" << G4endl;
                G4cout << "Old Momentum Direction: "
                     << aParticle->GetMomentumDirection() << G4endl;
                G4cout << "Old Polarization: "
                     << aParticle->GetPolarization() << G4endl;
        }

        // find polar angle w.r.t. old polarization vector

        G4double rand = G4UniformRand();

        G4double CosTheta = std::pow(rand, 1./3.);
        G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);

        if(G4UniformRand() < 0.5)CosTheta = -CosTheta;

        // find azimuthal angle w.r.t old polarization vector 

        rand = G4UniformRand();

        G4double Phi = twopi*rand;
        G4double SinPhi = std::sin(Phi); 
        G4double CosPhi = std::cos(Phi); 
        
        G4double unit_x = SinTheta * CosPhi; 
        G4double unit_y = SinTheta * SinPhi;  
        G4double unit_z = CosTheta; 
        
        G4ThreeVector NewPolarization (unit_x,unit_y,unit_z);

        // Rotate new polarization direction into global reference system 

        G4ThreeVector OldPolarization = aParticle->GetPolarization();
        OldPolarization = OldPolarization.unit();

        NewPolarization.rotateUz(OldPolarization);
        NewPolarization = NewPolarization.unit();
        
        // -- new momentum direction is normal to the new
        // polarization vector and in the same plane as the
        // old and new polarization vectors --

        G4ThreeVector NewMomentumDirection = 
                              OldPolarization - NewPolarization * CosTheta;

        if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection;
        NewMomentumDirection = NewMomentumDirection.unit();

        aParticleChange.ProposePolarization(NewPolarization);

        aParticleChange.ProposeMomentumDirection(NewMomentumDirection);

        if (verboseLevel>0) {
                G4cout << "New Polarization: " 
                     << NewPolarization << G4endl;
                G4cout << "Polarization Change: "
                     << *(aParticleChange.GetPolarization()) << G4endl;  
                G4cout << "New Momentum Direction: " 
                     << NewMomentumDirection << G4endl;
                G4cout << "Momentum Change: "
                     << *(aParticleChange.GetMomentumDirection()) << G4endl; 
        }

        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}

G4PhysicsTable * DsG4OpRayleigh::GetPhysicsTable

(

)

const [inline]

Definition at line 189 of file DsG4OpRayleigh.h.

{
  return thePhysicsTable;
}

void DsG4OpRayleigh::DumpPhysicsTable

(

)

const [inline]

Definition at line 176 of file DsG4OpRayleigh.h.

{
        G4int PhysicsTableSize = thePhysicsTable->entries();
        G4PhysicsOrderedFreeVector *v;

        for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
        {
                v = (G4PhysicsOrderedFreeVector*)(*thePhysicsTable)[i];
                v->DumpValues();
        }
}

void DsG4OpRayleigh::BuildThePhysicsTable

(

)

[private]

Definition at line 207 of file DsG4OpRayleigh.cc.

{
//      Builds a table of scattering lengths for each material

        if (thePhysicsTable) return;

        const G4MaterialTable* theMaterialTable=
                               G4Material::GetMaterialTable();
        G4int numOfMaterials = G4Material::GetNumberOfMaterials();

        // create a new physics table

        thePhysicsTable = new G4PhysicsTable(numOfMaterials);

        // loop for materials

        for (G4int i=0 ; i < numOfMaterials; i++)
        {
            G4PhysicsOrderedFreeVector* ScatteringLengths =
                                new G4PhysicsOrderedFreeVector();

            G4MaterialPropertiesTable *aMaterialPropertiesTable =
                         (*theMaterialTable)[i]->GetMaterialPropertiesTable();
                                                                                
            if(aMaterialPropertiesTable){

              G4MaterialPropertyVector* AttenuationLengthVector =
                            aMaterialPropertiesTable->GetProperty("RAYLEIGH");

              if(!AttenuationLengthVector){

                if ((*theMaterialTable)[i]->GetName() == "/dd/Materials/Water" || 
                    (*theMaterialTable)[i]->GetName() == "/dd/Materials/OwsWater" ||
                    (*theMaterialTable)[i]->GetName() == "/dd/Materials/IwsWater"
                    )
                {
                   // Call utility routine to Generate
                   // Rayleigh Scattering Lengths

                   DefaultWater = true;

                   ScatteringLengths =
                   RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
                }
              }
            }

            thePhysicsTable->insertAt(i,ScatteringLengths);
        } 
}

G4PhysicsOrderedFreeVector * DsG4OpRayleigh::RayleighAttenuationLengthGenerator

(

G4MaterialPropertiesTable *

aMPT

)

[private]

Definition at line 312 of file DsG4OpRayleigh.cc.

{
        // Physical Constants

        // isothermal compressibility of water
        G4double betat = 7.658e-23*m3/MeV;

        // K Boltzman
        G4double kboltz = 8.61739e-11*MeV/kelvin;

        // Temperature of water is 10 degrees celsius
        // conversion to kelvin:
        // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
        G4double temp = 283.15*kelvin;

        // Retrieve vectors for refraction index
        // and photon energy from the material properties table

        G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");

        G4double refsq;
        G4double e;
        G4double xlambda;
        G4double c1, c2, c3, c4;
        G4double Dist;
        G4double refraction_index;

        G4PhysicsOrderedFreeVector *RayleighScatteringLengths = 
                                new G4PhysicsOrderedFreeVector();

        if (Rindex ) {

           Rindex->ResetIterator();

           while (++(*Rindex)) {

                e = (Rindex->GetPhotonEnergy());

                refraction_index = Rindex->GetProperty();
                refsq = refraction_index*refraction_index;
                xlambda = h_Planck*c_light/e;

                if (verboseLevel>0) {
                        G4cout << Rindex->GetPhotonEnergy() << " MeV\t";
                        G4cout << xlambda << " mm\t";
                }

                c1 = 1 / (6.0 * pi);
                c2 = std::pow((2.0 * pi / xlambda), 4);
                c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
                c4 = betat * temp * kboltz;

                Dist = 1.0 / (c1*c2*c3*c4);

                if (verboseLevel>0) {
                        G4cout << Dist << " mm" << G4endl;
                }
                RayleighScatteringLengths->
                        InsertValues(Rindex->GetPhotonEnergy(), Dist);
           }

        }

        return RayleighScatteringLengths;
}

---

Member Data Documentation

G4PhysicsTable* DsG4OpRayleigh::thePhysicsTable [protected]

Definition at line 154 of file DsG4OpRayleigh.h.

G4bool DsG4OpRayleigh::DefaultWater [private]

Definition at line 161 of file DsG4OpRayleigh.h.

---

The documentation for this class was generated from the following files:

• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Simulation/DetSim/src/DsG4OpRayleigh.h
• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Simulation/DetSim/src/DsG4OpRayleigh.cc

---