diff options
Diffstat (limited to 'src/G4TankIceSD.cxx')
| -rw-r--r-- | src/G4TankIceSD.cxx | 166 |
1 files changed, 0 insertions, 166 deletions
diff --git a/src/G4TankIceSD.cxx b/src/G4TankIceSD.cxx deleted file mode 100644 index 2974837..0000000 --- a/src/G4TankIceSD.cxx +++ /dev/null @@ -1,166 +0,0 @@ -#include "G4TankIceSD.h" -#include "G4BeamTestTank.h" - -#include <G4Step.hh> -#include <G4HCofThisEvent.hh> -#include <G4TouchableHistory.hh> -#include <G4ios.hh> -#include <iostream> -#include <ios> -#include <fstream> - -#include <G4VProcess.hh> -#include <G4OpticalPhoton.hh> -#include "G4Poisson.hh" - -/* #include <icetray/I3Units.h> */ - -#include "G4Version.hh" - -#if G4VERSION_NUMBER >= 950 -// The G4MaterialPropertyVector is gone since 4.9.5. -// It has been typedef'd to a G4UnorderedPhysicsVector -// with a different interface. Try to support both -// versions with an ifdef. -#define MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR -#endif - - -G4TankIceSD::G4TankIceSD(G4String name/* , const std::map<OMKey, G4ThreeVector>& domPositions */) - : G4VSensitiveDetector(name)/* , domPositions_(domPositions) */ -{} - - -G4TankIceSD::~G4TankIceSD() {} - - -void G4TankIceSD::Initialize(G4HCofThisEvent* HCE) -{ - // std::map<OMKey, G4ThreeVector>::const_iterator om_iter; - // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) - // { - sumEdep_ = 0; - cogTime_ = 0; - cherenkovCounter_ = 0; - cherenkovCounterWeight_ = 0; - // } -} - - -G4bool G4TankIceSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) -{ - // Get energy deposit/time --> CHECK WHAT THESE QUANTITIES ACTUALLY MEAN !!!!!!!!!!!!!!!!!!!!!!!! - G4double edep = aStep->GetTotalEnergyDeposit(); - G4double time = aStep->GetPreStepPoint()->GetGlobalTime(); - G4double cherenkovNumber = GetCerenkovNumber(aStep); - - if(edep<=0 && cherenkovNumber<=0) return false; - - // Get Position relative to ice center - G4ThreeVector preStepPoint = aStep->GetPreStepPoint()->GetPosition(); - G4ThreeVector postStepPoint = aStep->GetPostStepPoint()->GetPosition(); - G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle(); - G4ThreeVector worldPosition = (preStepPoint + postStepPoint) / 2.0; - G4ThreeVector localPosition = theTouchable->GetHistory()-> - GetTopTransform().TransformPoint(worldPosition); - - // std::map<OMKey, G4ThreeVector>::const_iterator om_iter; - // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) - // { - G4double radius = sqrt(pow(domPositions_.x() - localPosition.x(), 2) + - pow(domPositions_.y() - localPosition.y(), 2)); - G4double height = (domPositions_.z() - localPosition.z()); - - sumEdep_ += edep; - cogTime_ += edep*time; - cherenkovCounterWeight_ += GetProbability(radius, height) * cherenkovNumber; - cherenkovCounter_ += cherenkovNumber; - // } - return true; -} - - -void G4TankIceSD::EndOfEvent(G4HCofThisEvent*) -{ - // std::map<OMKey, G4ThreeVector>::const_iterator om_iter; - // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) - // { - if(sumEdep_>0) - { - cogTime_ /= sumEdep_; - } - // } -} - - -G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep) -{ - // same function as in "source/processes/electromagnetic/xrays/src/G4Cerenkov.cc" but without - // cerenkov angle and only for an energy independent refraction index - G4Track* aTrack = aStep->GetTrack(); - const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle(); - const G4Material* aMaterial = aTrack->GetMaterial(); - - G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint(); - G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint(); - - G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable(); - if (!aMaterialPropertiesTable) return 0.0; - - G4MaterialPropertyVector* Rindex = aMaterialPropertiesTable->GetProperty("RINDEX"); - if (!Rindex) return 0.0; - - const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm); // TODO(shivesh): check this - const G4double charge = aParticle->GetDefinition()->GetPDGCharge(); - const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.; - if (beta <= 0.0) return 0.0; - G4double BetaInverse = 1. / beta; - - // Min and Max photon energies -#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR - G4double Pmin = Rindex->GetMinLowEdgeEnergy(); - G4double Pmax = Rindex->GetMaxLowEdgeEnergy(); -#else - G4double Pmin = Rindex->GetMinPhotonEnergy(); - G4double Pmax = Rindex->GetMaxPhotonEnergy(); -#endif - - // Min and Max Refraction Indices -#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR - G4double nMin = Rindex->GetMinValue(); - G4double nMax = Rindex->GetMaxValue(); -#else - G4double nMin = Rindex->GetMinProperty(); - G4double nMax = Rindex->GetMaxProperty(); -#endif - if (nMin!=nMax) - { - /* log_error("Support for energy dependent refraction index not yet implemented!"); */ - G4cout << "Support for energy dependent refraction index not yet implemented!"; - return 0.0; - } - - // If n(Pmax) < 1/Beta -- no photons generated - if (nMax < BetaInverse) return 0.0; - - G4double MeanNumberOfPhotons = Rfact * charge/CLHEP::eplus * charge/CLHEP::eplus * (Pmax - Pmin) * - (1. - BetaInverse * BetaInverse / nMin / nMin); - - G4double step_length = aStep->GetStepLength(); - - return MeanNumberOfPhotons * step_length; -} - - -// TODO(shivesh): what is this -G4double G4TankIceSD::GetProbability(G4double radius, G4double height) -{ - height = 0.90 - height / CLHEP::m; - radius = radius / CLHEP::m; - G4double p0 = 0.0340 + 0.0120 * height; - G4double p1 = 0.0400 + 0.000622 * exp(-1.13 + 8.03 * height); - G4double p2 = 0.917 + 2.34 * exp(-1.82 + 3.12 * height); - G4double p3 = -0.00325 - 0.00199 * height - 0.0121 * height*height; - - return (p0 + p3 * radius) + (p1 - p0) * exp(-p2*p2 * radius*radius); -} |