diff --git a/u4/CMakeLists.txt b/u4/CMakeLists.txt index 8290c5942a..4990168eec 100644 --- a/u4/CMakeLists.txt +++ b/u4/CMakeLists.txt @@ -24,8 +24,6 @@ set(SOURCES U4Cerenkov_Debug.cc U4Hit_Debug.cc - ShimG4OpAbsorption.cc - ShimG4OpRayleigh.cc Local_G4Cerenkov_modified.cc Local_DsG4Scintillation.cc @@ -88,9 +86,6 @@ set(HEADERS U4Cerenkov_Debug.hh U4Hit_Debug.hh - ShimG4OpAbsorption.hh - ShimG4OpRayleigh.hh - Local_G4Cerenkov_modified.hh Local_DsG4Scintillation.hh U4Physics.hh diff --git a/u4/Local_DsG4Scintillation.cc b/u4/Local_DsG4Scintillation.cc index d17cd9e03f..db4b3c1c51 100644 --- a/u4/Local_DsG4Scintillation.cc +++ b/u4/Local_DsG4Scintillation.cc @@ -84,22 +84,6 @@ #include "G4Electron.hh" #include "globals.hh" -#ifdef WITH_G4OPTICKS -#include "G4Opticks.hh" -#include "CGenstep.hh" -#include "CTrack.hh" -#include "CPhotonInfo.hh" -#include "SLOG.hh" -#endif - - - -#ifdef WITH_G4OPTICKS -//const plog::Severity Local_DsG4Scintillation::LEVEL = SLOG::EnvLevel("Local_DsG4Scintillation", "DEBUG") ; -const plog::Severity Local_DsG4Scintillation::LEVEL = error ; -#endif - - #include "U4Stack.h" #include "SEvt.hh" @@ -582,28 +566,6 @@ Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep } -//-------------------------------------------------// - -#ifdef WITH_G4OPTICKS - - CPho ancestor = CPhotonInfo::Get(&aTrack); - int ancestor_id = ancestor.get_id() ; - - - /** - ancestor_id:-1 - normal case, meaning that aTrack was not a photon - so the generation loop will yield "primary" photons - with id : gs.offset + i - - ancestor_id>-1 - aTrack is a photon that may be about to reemit (Num=0 or 1) - ancestor_id is the absolute id of the primary parent photon, - this id is retained thru any subsequent remission secondary generations - **/ - -#endif - #ifdef STANDALONE U4::GenPhotonAncestor(&aTrack); #endif @@ -639,18 +601,6 @@ Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep G4int NumPhoton = NumVec[scnt] ; - -#ifdef WITH_G4OPTICKS - if(flagReemission) assert( NumPhoton == 0 || NumPhoton == 1); // expecting only 0 or 1 remission photons - bool is_opticks_genstep = NumPhoton > 0 && !flagReemission ; - - CGenstep gs ; - if(is_opticks_genstep && (m_opticksMode & 1)) - { - gs = G4Opticks::Get()->collectGenstep_Local_DsG4Scintillation_r4695( &aTrack, &aStep, NumPhoton, scnt, ScintillationTime); - } -#endif - #ifdef STANDALONE if(flagReemission) assert( NumPhoton == 0 || NumPhoton == 1); // expecting only 0 or 1 remission photons bool is_opticks_genstep = NumPhoton > 0 && !flagReemission ; @@ -665,9 +615,6 @@ Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep { for(G4int i = 0 ; i < NumPhoton ; i++) { -#ifdef WITH_G4OPTICKS - G4Opticks::Get()->setAlignIndex( ancestor_id > -1 ? ancestor_id : gs.offset + i ); // aka photon_id -#endif #ifdef STANDALONE U4::GenPhotonBegin(i); #endif @@ -845,11 +792,6 @@ Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep if ( verboseLevel > 0 ) { G4cout << " aSecondaryTrack->SetWeight( " << weight<< " ) ; aSecondaryTrack->GetWeight() = " << aSecondaryTrack->GetWeight() << G4endl;} -#ifdef WITH_G4OPTICKS - aSecondaryTrack->SetUserInformation(CPhotonInfo::MakeScintillation(gs, i, ancestor )); - G4Opticks::Get()->setAlignIndex(-1); -#endif - #ifdef STANDALONE U4::GenPhotonEnd(i, aSecondaryTrack); #endif @@ -875,56 +817,6 @@ Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep return change ; } - - -#ifdef WITH_G4OPTICKS -G4MaterialPropertyVector* Local_DsG4Scintillation::getMaterialProperty(const char* name, G4int materialIndex) -{ - const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); - G4int numOfMaterials = G4Material::GetNumberOfMaterials(); - assert( materialIndex < numOfMaterials ); - - G4Material* aMaterial = (*theMaterialTable)[materialIndex]; - G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable(); - G4MaterialPropertyVector* prop = aMaterialPropertiesTable ? aMaterialPropertiesTable->GetProperty(name) : nullptr ; - return prop ; -} - -G4PhysicsOrderedFreeVector* Local_DsG4Scintillation::getScintillationIntegral(G4int scnt, G4int materialIndex) const -{ - G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL; - - if ( scnt == 0 ){ - ScintillationIntegral = - (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex)); - } - else{ - ScintillationIntegral = - (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex)); - } - - return ScintillationIntegral ; -} - - -G4double Local_DsG4Scintillation::getSampledEnergy(G4int scnt, G4int materialIndex) const -{ - G4PhysicsOrderedFreeVector* ScintillationIntegral = getScintillationIntegral(scnt, materialIndex); - G4double CIIvalue = G4UniformRand()*ScintillationIntegral->GetMaxValue(); - G4double sampledEnergy = ScintillationIntegral->GetEnergy(CIIvalue); - return sampledEnergy ; -} - -G4double Local_DsG4Scintillation::getSampledWavelength(G4int scnt, G4int materialIndex) const -{ - G4double sampledEnergy = getSampledEnergy(scnt, materialIndex ); - G4double wavelength = h_Planck*c_light/sampledEnergy ; - G4double wavelength_nm = wavelength/nm ; - return wavelength_nm ; -} -#endif - - // BuildThePhysicsTable for the scintillation process // -------------------------------------------------- // diff --git a/u4/Local_DsG4Scintillation.hh b/u4/Local_DsG4Scintillation.hh index e64ac93575..7500b34519 100644 --- a/u4/Local_DsG4Scintillation.hh +++ b/u4/Local_DsG4Scintillation.hh @@ -70,10 +70,6 @@ #define Local_DsG4Scintillation_h 1 -#ifdef WITH_G4OPTICKS -#include "plog/Severity.h" -#endif - #include "globals.hh" #include "templates.hh" #include "Randomize.hh" @@ -111,11 +107,6 @@ class G4UIdirectory; class Local_DsG4Scintillation : public G4VRestDiscreteProcess, public G4UImessenger { //too lazy to create another UImessenger class -private: -#ifdef WITH_G4OPTICKS - static const plog::Severity LEVEL ; -#endif - public: static const bool FLOAT ; static const int PIDX ; @@ -194,13 +185,6 @@ public: // With description void SetUseFastMu300nsTrick(const G4bool fastMu300nsTrick); G4bool GetUseFastMu300nsTrick() const; -#ifdef WITH_G4OPTICKS - G4MaterialPropertyVector* getMaterialProperty(const char* name, G4int materialIndex) ; - G4PhysicsOrderedFreeVector* getScintillationIntegral(G4int scnt, G4int materialIndex) const; - G4double getSampledEnergy(G4int scnt, G4int materialIndex) const ; - G4double getSampledWavelength(G4int scnt, G4int materialIndex) const ; -#endif - void SetScintillationExcitationRatio(const G4double excitationratio); // Called to set the scintillation exciation ratio, needed when // the scintillation level excitation is different for different diff --git a/u4/Local_G4Cerenkov_modified.cc b/u4/Local_G4Cerenkov_modified.cc index 8224872287..9bb4aeb8f5 100644 --- a/u4/Local_G4Cerenkov_modified.cc +++ b/u4/Local_G4Cerenkov_modified.cc @@ -75,11 +75,6 @@ #include "Local_G4Cerenkov_modified.hh" -#ifdef INSTRUMENTED -#include "OpticksDebug.hh" -#include "OpticksRandom.hh" -#endif - #ifdef STANDALONE #include "SLOG.hh" #include "U4.hh" @@ -112,9 +107,6 @@ Local_G4Cerenkov_modified::Local_G4Cerenkov_modified(const G4String& processName fMaxBetaChange(0.0), fMaxPhotons(0), fStackingFlag(true), -#ifdef INSTRUMENTED - override_fNumPhotons(0), -#endif fNumPhotons(0) { SetProcessSubType(fCerenkov); @@ -243,13 +235,6 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack fNumPhotons = (G4int) G4Poisson(MeanNumberOfPhotons); -#ifdef INSTRUMENTED - if( override_fNumPhotons > 0 ) - { - fNumPhotons = override_fNumPhotons ; - } -#endif - // calculate the fNumPhotons1 and fNumPhotons2 { @@ -317,19 +302,6 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack fNumPhotons2 = MeanNumberOfPhotons2; -#ifdef INSTRUMENTED - par->append( BetaInverse, "BetaInverse" ); - par->append( beta , "beta" ); - par->append( Pmin , "Pmin" ); - par->append( Pmax , "Pmax" ); - - par->append( nMax , "nMax" ); - par->append( maxCos , "maxCos" ); - par->append( maxSin2 , "maxSin2" ); - par->append( fNumPhotons, "fNumPhotons" ); -#endif - - if (MeanNumberOfPhotons1 <= 0.0 or MeanNumberOfPhotons2<=0.0) { @@ -411,36 +383,9 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack G4double sin2Theta(0.); #endif -#ifdef INSTRUMENTED - unsigned head_count = 0 ; - unsigned tail_count = 0 ; - unsigned continue_count = 0 ; - unsigned condition_count = 0 ; - int seqidx = -1 ; - if(rnd) - { - rnd->setSequenceIndex(i); - seqidx = rnd->getSequenceIndex(); - - if(i < 10) std::cout - << " i " << std::setw(6) << i - << " seqidx " << std::setw(7) << seqidx - << " Pmin/eV " << std::fixed << std::setw(10) << std::setprecision(5) << Pmin/eV - << " Pmax/eV " << std::fixed << std::setw(10) << std::setprecision(5) << Pmax/eV - << " dp/eV " << std::fixed << std::setw(10) << std::setprecision(5) << dp/eV - << " maxSin2 " << std::fixed << std::setw(10) << std::setprecision(5) << maxSin2 - << std::endl - ; - - } -#endif - // sample an energy do { -#ifdef INSTRUMENTED - head_count += 1 ; -#endif rand0 = G4UniformRand(); sampledEnergy = Pmin + rand0 * dp; sampledRI = Rindex->Value(sampledEnergy); @@ -450,17 +395,10 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack #else // check if n(E) > 1/beta if (sampledRI < BetaInverse) { -#ifdef INSTRUMENTED - continue_count += 1 ; -#endif continue; } #endif - -#ifdef INSTRUMENTED - tail_count += 1 ; -#endif cosTheta = BetaInverse / sampledRI; @@ -483,46 +421,7 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack #endif // Loop checking, 07-Aug-2015, Vladimir Ivanchenko -#ifdef INSTRUMENTED - - if( i < 10 ) std::cout - << " tc " << std::setw(6) << tail_count - << " u0 " << std::fixed << std::setw(10) << std::setprecision(5) << rand0 - << " eV " << std::fixed << std::setw(10) << std::setprecision(5) << sampledEnergy/eV - << " ri " << std::fixed << std::setw(10) << std::setprecision(5) << sampledRI - << " ct " << std::fixed << std::setw(10) << std::setprecision(5) << cosTheta - << " s2 " << std::fixed << std::setw(10) << std::setprecision(5) << sin2Theta - << " rand1*maxSin2 " << std::fixed << std::setw(10) << std::setprecision(5) << rand1*maxSin2 - << " rand1*maxSin2 - sin2Theta " << std::fixed << std::setw(10) << std::setprecision(5) << rand1*maxSin2 - sin2Theta - << " loop " << ( rand1*maxSin2 > sin2Theta ? "Y" : "N" ) - << std::endl - ; - - - } while ( looping_condition(condition_count) && rand1*maxSin2 > sin2Theta ); -#else } while (rand1*maxSin2 > sin2Theta); -#endif - -#ifdef INSTRUMENTED - G4double sampledEnergy_eV = sampledEnergy/eV ; - G4double sampledWavelength_nm = h_Planck*c_light/sampledEnergy/nm ; - - gen->append( sampledEnergy_eV , "sampledEnergy" ); - gen->append( sampledWavelength_nm , "sampledWavelength" ); - gen->append( sampledRI , "sampledRI" ); - gen->append( cosTheta , "cosTheta" ); - - gen->append( sin2Theta , "sin2Theta" ); - gen->append( head_count , tail_count, "head_tail" ); - gen->append( continue_count , condition_count, "continue_condition" ); - gen->append( BetaInverse , "BetaInverse" ); - - if(rnd) - { - rnd->setSequenceIndex(-1); - } -#endif @@ -630,17 +529,6 @@ G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack } -#ifdef INSTRUMENTED -bool Local_G4Cerenkov_modified::looping_condition(unsigned& count) -{ - count += 1 ; - return true ; -} - - -#endif - - // BuildThePhysicsTable for the Cerenkov process // --------------------------------------------- // @@ -872,17 +760,6 @@ G4double const G4double Rfact = 369.81/(eV * cm); -#ifdef X_INSTRUMENTED - std::cout - << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons" - << " Rfact " << std::fixed << std::setw(10) << std::setprecision(5) << Rfact - << " eV " << std::fixed << std::setw(10) << std::setprecision(7) << eV - << " cm " << std::fixed << std::setw(10) << std::setprecision(5) << cm - << " charge " << std::fixed << std::setw(10) << std::setprecision(5) << charge - << std::endl - ; -#endif - if(beta <= 0.0)return 0.0; G4double BetaInverse = 1./beta; @@ -1133,40 +1010,6 @@ G4double G4double NumPhotons = Rfact * charge/eplus * charge/eplus * (dp1 - ge1 * BetaInverse*BetaInverse); - -#ifdef X_INSTRUMENTED - std::cout - << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons" - << " BetaInverse " << std::fixed << std::setw(10) << std::setprecision(5) << BetaInverse - << " maxRI " << std::fixed << std::setw(10) << std::setprecision(5) << maxRI - << " minRI " << std::fixed << std::setw(10) << std::setprecision(5) << minRI - << " cross_num " << cross_num - << " dp1 " << std::fixed << std::setw(10) << std::setprecision(5) << dp1 - << " dp1/eV " << std::fixed << std::setw(10) << std::setprecision(5) << dp1/eV - << " ge1 " << std::fixed << std::setw(10) << std::setprecision(5) << ge1 - << " NumPhotons " << std::fixed << std::setw(10) << std::setprecision(5) << NumPhotons - << std::endl - ; - - for(int i=0 ; i < cross_num ; i++) - { - - G4bool isOutRange; - G4double cai0 = CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i+0], isOutRange); - G4double cai1 = CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i+1], isOutRange); - - std::cout - << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons" - << " the_energies_threshold[2*i+0]/eV " << std::fixed << std::setw(10) << std::setprecision(5) << the_energies_threshold[2*i+0]/eV - << " the_energies_threshold[2*i+1]/eV " << std::fixed << std::setw(10) << std::setprecision(5) << the_energies_threshold[2*i+1]/eV - << " cai0 " << std::fixed << std::setw(20) << std::setprecision(10) << cai0 - << " cai1 " << std::fixed << std::setw(20) << std::setprecision(10) << cai1 - << std::endl - ; - } -#endif - - return NumPhotons; } @@ -1282,19 +1125,6 @@ G4double Local_G4Cerenkov_modified::GetAverageNumberOfPhotons_s2(const G4double const G4double Rfact = 369.81/(eV * cm); G4double NumPhotons = Rfact * charge/eplus * charge/eplus * s2integral ; -#ifdef X_INSTRUMENTED - std::cout - << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons_s2" - << " Rfact " << std::fixed << std::setw(10) << std::setprecision(5) << Rfact - << " eV " << std::fixed << std::setw(10) << std::setprecision(7) << eV - << " cm " << std::fixed << std::setw(10) << std::setprecision(5) << cm - << " mm " << std::fixed << std::setw(10) << std::setprecision(5) << mm - << " charge " << std::fixed << std::setw(10) << std::setprecision(5) << charge - << " s2integral " << std::fixed << std::setw(10) << std::setprecision(5) << s2integral - << " NumPhotons " << std::fixed << std::setw(10) << std::setprecision(5) << NumPhotons - << std::endl - ; -#endif return NumPhotons ; } @@ -1310,4 +1140,3 @@ void Local_G4Cerenkov_modified::DumpPhysicsTable() const } } - diff --git a/u4/Local_G4Cerenkov_modified.hh b/u4/Local_G4Cerenkov_modified.hh index 387565c67e..6d7b002744 100644 --- a/u4/Local_G4Cerenkov_modified.hh +++ b/u4/Local_G4Cerenkov_modified.hh @@ -78,25 +78,10 @@ // Class Definition ///////////////////// -#ifdef INSTRUMENTED -template struct OpticksDebug ; -struct OpticksRandom ; -#endif - - class Local_G4Cerenkov_modified : public G4VProcess { public: -#ifdef INSTRUMENTED - template friend struct Local_G4Cerenkov_modifiedTest ; - bool looping_condition(unsigned& count); - OpticksDebug* gen ; - OpticksDebug* par ; - OpticksRandom* rnd ; -#endif - - //////////////////////////////// // Constructors and Destructor //////////////////////////////// @@ -252,12 +237,6 @@ private: G4bool fStackingFlag; -#ifdef INSTRUMENTED -public: - G4int override_fNumPhotons ; -private: -#endif - G4int fNumPhotons; G4int fNumPhotons1; // mean G4int fNumPhotons2; // mean @@ -325,5 +304,3 @@ G4PhysicsTable* Local_G4Cerenkov_modified::GetPhysicsTable() const } #endif /* Local_G4Cerenkov_modified_h */ - - diff --git a/u4/ShimG4OpAbsorption.cc b/u4/ShimG4OpAbsorption.cc deleted file mode 100644 index c608a0b156..0000000000 --- a/u4/ShimG4OpAbsorption.cc +++ /dev/null @@ -1,147 +0,0 @@ -#include - -#include "G4SystemOfUnits.hh" - -#include "ShimG4OpAbsorption.hh" -#include "SEvt.hh" -#include "SLOG.hh" -#include "ssys.h" -#include "U4UniformRand.h" -#include "U4Stack.h" - - -const plog::Severity ShimG4OpAbsorption::LEVEL = SLOG::EnvLevel("ShimG4OpAbsorption", "DEBUG") ; - - -ShimG4OpAbsorption::ShimG4OpAbsorption(const G4String& processName, G4ProcessType type ) - : - G4OpAbsorption(processName, type) -{ -} - -ShimG4OpAbsorption::~ShimG4OpAbsorption(){} - - - - -//#ifdef DEBUG_TAG -/** -ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft --------------------------------------------------------- - -Shim makes process classname appear in SBacktrace.h enabling U4Random::flat/U4Stack::Classify - -**/ - -const bool ShimG4OpAbsorption::FLOAT = ssys::getenvbool("ShimG4OpAbsorption_FLOAT") ; -const int ShimG4OpAbsorption::PIDX = ssys::getenvint("PIDX",-1); -const bool ShimG4OpAbsorption::PIDX_ENABLED = ssys::getenvbool("ShimG4OpAbsorption__PIDX_ENABLED") ; - - -// void ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft(){ G4VProcess::ResetNumberOfInteractionLengthLeft(); } - void ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft() -{ - G4double u = G4UniformRand() ; - - LOG(LEVEL) - << U4UniformRand::Desc(u, SEvt::UU ) - ; - -#ifndef PRODUCTION - SEvt::AddTag( 1, U4Stack_AbsorptionDiscreteReset, u ); -#endif - - if(FLOAT) - { - float f = -1.f*std::log( float(u) ) ; - theNumberOfInteractionLengthLeft = f ; - } - else - { - theNumberOfInteractionLengthLeft = -1.*G4Log(u) ; - } - theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; - -} - - G4double ShimG4OpAbsorption::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition* ) -{ - G4double len = G4OpAbsorption::GetMeanFreePath( aTrack, 0, 0 ); - //std::cout << "ShimG4OpAbsorption::GetMeanFreePath " << len << std::endl ; - - //std::raise(SIGINT); - - return len ; -} - - - G4double ShimG4OpAbsorption::PostStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4ForceCondition* condition) -{ - if ( (previousStepSize < 0.0) || (theNumberOfInteractionLengthLeft<=0.0)) { - // beggining of tracking (or just after DoIt of this process) - ResetNumberOfInteractionLengthLeft(); - } else if ( previousStepSize > 0.0) { - // subtract NumberOfInteractionLengthLeft - SubtractNumberOfInteractionLengthLeft(previousStepSize); - } else { - // zero step - // DO NOTHING - } - - // condition is set to "Not Forced" - *condition = NotForced; - - // get mean free path - currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition); - - G4double value; - if (currentInteractionLength 1){ - G4cout << "G4VDiscreteProcess::PostStepGetPhysicalInteractionLength "; - G4cout << "[ " << GetProcessName() << "]" <DumpInfo(); - G4cout << " in Material " << track.GetMaterial()->GetName() < 0.0) { - // subtract NumberOfInteractionLengthLeft - SubtractNumberOfInteractionLengthLeft(previousStepSize); - } else { - // zero step - // DO NOTHING - } - - // condition is set to "Not Forced" - *condition = NotForced; - - // get mean free path - currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition); - - G4double value; - if (currentInteractionLength 1){ - G4cout << "G4VDiscreteProcess::PostStepGetPhysicalInteractionLength "; - G4cout << "[ " << GetProcessName() << "]" <DumpInfo(); - G4cout << " in Material " << track.GetMaterial()->GetName() <0) { - G4cout << "Scattering Photon!" << G4endl; - G4cout << "Old Momentum Direction: " - << aParticle->GetMomentumDirection() << G4endl; - G4cout << "Old Polarization: " - << aParticle->GetPolarization() << G4endl; - } - - G4double cosTheta; - G4ThreeVector OldMomentumDirection, NewMomentumDirection; - G4ThreeVector OldPolarization, NewPolarization; - - G4double rand, constant; - G4double CosTheta, SinTheta, SinPhi, CosPhi, unit_x, unit_y, unit_z; - - G4double u ; - G4double u_loopexit ; - - - do { - // Try to simulate the scattered photon momentum direction - // w.r.t. the initial photon momentum direction - - CosTheta = G4UniformRand(); - -#ifndef PRODUCTION -#ifdef DEBUG_TAG - SEvt::AddTag( 1, U4Stack_RayleighScatter, CosTheta ); // 0 -#endif -#endif - - SinTheta = std::sqrt(1.-CosTheta*CosTheta); - // consider for the angle 90-180 degrees - - u = G4UniformRand() ; - -#ifndef PRODUCTION -#ifdef DEBUG_TAG - SEvt::AddTag( 1, U4Stack_RayleighScatter, u ); // 1 -#endif -#endif - - if (u < 0.5) CosTheta = -CosTheta; - - // simulate the phi angle - - u = G4UniformRand() ; - -#ifndef PRODUCTION -#ifdef DEBUG_TAG - SEvt::AddTag( 1, U4Stack_RayleighScatter, u ); // 2 -#endif -#endif - - rand = twopi*u; - SinPhi = std::sin(rand); - CosPhi = std::cos(rand); - - // start constructing the new momentum direction - unit_x = SinTheta * CosPhi; - unit_y = SinTheta * SinPhi; - unit_z = CosTheta; - NewMomentumDirection.set (unit_x,unit_y,unit_z); - - // Rotate the new momentum direction into global reference system - OldMomentumDirection = aParticle->GetMomentumDirection(); - OldMomentumDirection = OldMomentumDirection.unit(); - NewMomentumDirection.rotateUz(OldMomentumDirection); - NewMomentumDirection = NewMomentumDirection.unit(); - - // calculate the new polarization direction - // The new polarization needs to be in the same plane as the new - // momentum direction and the old polarization direction - OldPolarization = aParticle->GetPolarization(); - constant = -NewMomentumDirection.dot(OldPolarization); - - NewPolarization = OldPolarization + constant*NewMomentumDirection; - NewPolarization = NewPolarization.unit(); - - // There is a corner case, where the Newmomentum direction - // is the same as oldpolariztion direction: - // random generate the azimuthal angle w.r.t. Newmomentum direction - - u = G4UniformRand() ; - -#ifndef PRODUCTION -#ifdef DEBUG_TAG - SEvt::AddTag( 1, U4Stack_RayleighScatter, u ); // 3 -#endif -#endif - - if (NewPolarization.mag() == 0.) { - rand = u*twopi; - NewPolarization.set(std::cos(rand),std::sin(rand),0.); - NewPolarization.rotateUz(NewMomentumDirection); - } else { - // There are two directions which are perpendicular - // to the new momentum direction - if (u < 0.5) NewPolarization = -NewPolarization; - } - - // simulate according to the distribution cos^2(theta) - cosTheta = NewPolarization.dot(OldPolarization); - - u_loopexit = G4UniformRand() ; - -#ifndef PRODUCTION -#ifdef DEBUG_TAG - SEvt::AddTag( 1, U4Stack_RayleighScatter, u_loopexit ); // 4 -#endif -#endif - - // Loop checking, 13-Aug-2015, Peter Gumplinger - } while (std::pow(cosTheta,2) < u_loopexit ); - - 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); -} - - - - diff --git a/u4/ShimG4OpRayleigh.hh b/u4/ShimG4OpRayleigh.hh deleted file mode 100644 index d3114c183a..0000000000 --- a/u4/ShimG4OpRayleigh.hh +++ /dev/null @@ -1,42 +0,0 @@ -#pragma once -/** -ShimG4OpRayleigh -=================== - -class G4OpRayleigh : public G4VDiscreteProcess - -**/ - -#include "plog/Severity.h" -#include "G4OpRayleigh.hh" -#include "U4_API_EXPORT.hh" - -class U4_API ShimG4OpRayleigh : public G4OpRayleigh -{ - static const plog::Severity LEVEL ; - public: - ShimG4OpRayleigh(); - virtual ~ShimG4OpRayleigh(); - - static const bool FLOAT ; - static const int PIDX ; - static const bool PIDX_ENABLED ; - - void ResetNumberOfInteractionLengthLeft(); - G4double GetMeanFreePath(const G4Track& aTrack, - G4double , - G4ForceCondition* ); - - G4double PostStepGetPhysicalInteractionLength( - const G4Track& track, - G4double previousStepSize, - G4ForceCondition* condition - ); - - - G4VParticleChange* PostStepDoIt(const G4Track& aTrack, - const G4Step& aStep); - -}; - - diff --git a/u4/U4Physics.cc b/u4/U4Physics.cc index 9d0caf334d..734d0e239b 100644 --- a/u4/U4Physics.cc +++ b/u4/U4Physics.cc @@ -20,9 +20,9 @@ Boundary class changes need to match in all the below:: #include "G4ProcessManager.hh" #include "G4FastSimulationManagerProcess.hh" - +#include "G4OpAbsorption.hh" #include "G4OpBoundaryProcess.hh" - +#include "G4OpRayleigh.hh" #include "SLOG.hh" const plog::Severity U4Physics::LEVEL = SLOG::EnvLevel("U4Physics", "DEBUG") ; @@ -153,9 +153,6 @@ void U4Physics::ConstructEM() #include "Local_G4Cerenkov_modified.hh" #include "Local_DsG4Scintillation.hh" -#include "ShimG4OpAbsorption.hh" -#include "ShimG4OpRayleigh.hh" - std::string U4Physics::desc() const { @@ -177,14 +174,7 @@ std::string U4Physics::desc() const std::string U4Physics::Desc() // static { - std::stringstream ss ; -#ifdef DEBUG_TAG - ss << ( ShimG4OpAbsorption::FLOAT ? "ShimG4OpAbsorption_FLOAT" : "ShimG4OpAbsorption_ORIGINAL" ) ; - ss << "_" ; - ss << ( ShimG4OpRayleigh::FLOAT ? "ShimG4OpRayleigh_FLOAT" : "ShimG4OpRayleigh_ORIGINAL" ) ; -#endif - std::string str = ss.str(); - return str ; + return {}; } @@ -260,20 +250,12 @@ void U4Physics::ConstructOp() if(OpAbsorption_DISABLE == 0) { -#ifdef DEBUG_TAG - fAbsorption = new ShimG4OpAbsorption(); -#else fAbsorption = new G4OpAbsorption(); -#endif } if(OpRayleigh_DISABLE == 0) { -#ifdef DEBUG_TAG - fRayleigh = new ShimG4OpRayleigh(); -#else fRayleigh = new G4OpRayleigh(); -#endif } if(OpBoundaryProcess_DISABLE == 0) diff --git a/u4/U4Physics.hh b/u4/U4Physics.hh index 0f5754c4cc..19c198e3d4 100644 --- a/u4/U4Physics.hh +++ b/u4/U4Physics.hh @@ -17,13 +17,8 @@ This is intended solely for use from U4AppTest class Local_G4Cerenkov_modified ; class Local_DsG4Scintillation ; -#ifdef DEBUG_TAG -class ShimG4OpAbsorption ; -class ShimG4OpRayleigh ; -#else class G4OpAbsorption ; class G4OpRayleigh ; -#endif class G4VProcess ; class G4ProcessManager ; @@ -40,13 +35,8 @@ struct U4_API U4Physics : public G4VUserPhysicsList Local_G4Cerenkov_modified* fCerenkov ; Local_DsG4Scintillation* fScintillation ; -#ifdef DEBUG_TAG - ShimG4OpAbsorption* fAbsorption ; - ShimG4OpRayleigh* fRayleigh ; -#else G4OpAbsorption* fAbsorption ; G4OpRayleigh* fRayleigh ; -#endif G4VProcess* fBoundary ; G4FastSimulationManagerProcess* fFastSim ; @@ -80,5 +70,3 @@ struct U4_API U4Physics : public G4VUserPhysicsList int OpBoundaryProcess_LASTPOST = 0 ; int FastSim_ENABLE = 0 ; }; - - diff --git a/u4/U4Stack.h b/u4/U4Stack.h index 64dbef6e13..a929237e31 100644 --- a/u4/U4Stack.h +++ b/u4/U4Stack.h @@ -47,8 +47,8 @@ struct U4Stack static unsigned TagToStack(unsigned tag); static constexpr const char* Unclassified_ = "Unclassified" ; // 0 - static constexpr const char* RestDiscreteReset_ = "RestDiscreteReset" ; // 1: not used as use Shims to idenify processes in SBacktrace - static constexpr const char* DiscreteReset_ = "DiscreteReset" ; // 2 : G4OpAbsoption G4OpRayleigh but now use Shims to identify + static constexpr const char* RestDiscreteReset_ = "RestDiscreteReset"; // 1 + static constexpr const char* DiscreteReset_ = "DiscreteReset"; // 2: G4OpAbsorption/G4OpRayleigh share this reset path static constexpr const char* ScintDiscreteReset_ = "ScintDiscreteReset" ; // 3 static constexpr const char* BoundaryDiscreteReset_ = "BoundaryDiscreteReset" ; // 4 static constexpr const char* RayleighDiscreteReset_ = "RayleighDiscreteReset" ; // 5 @@ -171,4 +171,3 @@ inline unsigned U4Stack::TagToStack(unsigned tag) } return stack ; } - diff --git a/u4/U4StackAuto.h b/u4/U4StackAuto.h index 1abcfeb07c..1c478cc0e9 100644 --- a/u4/U4StackAuto.h +++ b/u4/U4StackAuto.h @@ -59,68 +59,6 @@ G4SteppingManager::DefinePhysicalStepLength G4SteppingManager::Stepping )" ; - - static constexpr const char* RayleighDiscreteReset = R"( -U4Random::flat -G4VProcess::ResetNumberOfInteractionLengthLeft -ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft -G4VDiscreteProcess::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)" ; - - static constexpr const char* ShimRayleighDiscreteReset_ = "ShimRayleighDiscreteReset" ; // 5 - static constexpr const char* ShimRayleighDiscreteReset = R"( -U4Random::flat -ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft -G4VDiscreteProcess::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)" ; - - static constexpr const char* Shim2RayleighDiscreteReset_ = "Shim2RayleighDiscreteReset" ; // 5 - static constexpr const char* Shim2RayleighDiscreteReset = R"( -U4Random::flat -ShimG4OpRayleigh::ResetNumberOfInteractionLengthLeft -ShimG4OpRayleigh::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)" ; - - - - static constexpr const char* AbsorptionDiscreteReset = R"( -U4Random::flat -G4VProcess::ResetNumberOfInteractionLengthLeft -ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft -G4VDiscreteProcess::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)" ; - static constexpr const char* ShimAbsorptionDiscreteReset_ = "ShimAbsorptionDiscreteReset" ; // 6 - static constexpr const char* ShimAbsorptionDiscreteReset = R"( -U4Random::flat -ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft -G4VDiscreteProcess::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)" ; - - static constexpr const char* Shim2AbsorptionDiscreteReset_ = "Shim2AbsorptionDiscreteReset" ; // 6 - static constexpr const char* Shim2AbsorptionDiscreteReset = R"( -U4Random::flat -ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft -ShimG4OpAbsorption::PostStepGetPhysicalInteractionLength -G4VProcess::PostStepGPIL -G4SteppingManager::DefinePhysicalStepLength -G4SteppingManager::Stepping -)"; - static constexpr const char *BoundaryBurn_SurfaceReflectTransmitAbsorb = R"( U4Random::flat G4OpBoundaryProcess::PostStepDoIt @@ -173,14 +111,6 @@ inline unsigned U4StackAuto::Classify(const char* summary) if(strstr(summary, BoundaryDiscreteReset)) stack = U4Stack_BoundaryDiscreteReset ; if(strstr(summary, BoundaryDiscreteReset2)) stack = U4Stack_BoundaryDiscreteReset ; - if(strstr(summary, RayleighDiscreteReset)) stack = U4Stack_RayleighDiscreteReset ; - if(strstr(summary, ShimRayleighDiscreteReset)) stack = U4Stack_RayleighDiscreteReset ; - if(strstr(summary, Shim2RayleighDiscreteReset)) stack = U4Stack_RayleighDiscreteReset ; - - if(strstr(summary, AbsorptionDiscreteReset)) stack = U4Stack_AbsorptionDiscreteReset ; - if(strstr(summary, ShimAbsorptionDiscreteReset)) stack = U4Stack_AbsorptionDiscreteReset ; - if(strstr(summary, Shim2AbsorptionDiscreteReset)) stack = U4Stack_AbsorptionDiscreteReset ; - if(strstr(summary, BoundaryBurn_SurfaceReflectTransmitAbsorb)) stack = U4Stack_BoundaryBurn_SurfaceReflectTransmitAbsorb ; if(strstr(summary, BoundaryDiDiTransCoeff)) stack = U4Stack_BoundaryDiDiTransCoeff ; if(strstr(summary, AbsorptionEffDetect)) stack = U4Stack_AbsorptionEffDetect ; @@ -192,4 +122,3 @@ inline bool U4StackAuto::IsClassified(unsigned stack) { return stack != U4Stack_Unclassified ; } - diff --git a/u4/tests/U4SimulateTest.sh b/u4/tests/U4SimulateTest.sh index cb7bca0f02..5f9dd78e6f 100755 --- a/u4/tests/U4SimulateTest.sh +++ b/u4/tests/U4SimulateTest.sh @@ -144,8 +144,6 @@ loglevel(){ export junoPMTOpticalModelSimple=INFO #export SEvt=INFO ## thus is exceedingly verbose export SEventConfig=INFO - export ShimG4OpAbsorption=INFO - export ShimG4OpRayleigh=INFO } diff --git a/u4/u4s.sh b/u4/u4s.sh index 800b7c57b5..f790bcb213 100755 --- a/u4/u4s.sh +++ b/u4/u4s.sh @@ -97,16 +97,7 @@ fi export A_FOLD # A_FOLD is needed for loading "$A_FOLD/sframe.npy" export A_CFBASE # A_CFBASE needed for loading "$A_CFBASE/CSGFoundry/SSim" -export ShimG4OpAbsorption_FLOAT=1 -export ShimG4OpRayleigh_FLOAT=1 - -# cf U4Physics::Desc -physdesc="" -[ -n "$ShimG4OpAbsorption_FLOAT" ] && physdesc="${physdesc}ShimG4OpAbsorption_FLOAT" -[ -z "$ShimG4OpAbsorption_FLOAT" ] && physdesc="${physdesc}ShimG4OpAbsorption_ORIGINAL" -physdesc="${physdesc}_" -[ -n "$ShimG4OpRayleigh_FLOAT" ] && physdesc="${physdesc}ShimG4OpRayleigh_FLOAT" -[ -z "$ShimG4OpRayleigh_FLOAT" ] && physdesc="${physdesc}ShimG4OpRayleigh_ORIGINAL" +physdesc="G4OpAbsorption_G4OpRayleigh" #sel=PIDX_0_ sel=ALL