Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • develop
  • master
  • revert-39cd0316
3 results

Target

Select target project
  • lphe-lab-tools/SciFiMatG4
  • bleverin/SciFiMatG4_v2
2 results
Select Git revision
  • develop
  • master
  • revert-39cd0316
3 results
Show changes
Commits on Source (16)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1073 additions and 2131 deletions
README.md
View file @ 07fea035
......@@ -6,35 +6,41 @@ This is based on code from Mirco Deckenhoff and modified by Blake Leverington, P
_Installation of GEANT4 help is at the bottom of this readme. Check the cmake options for compatibility._
> $ git clone https://$USER@gitlab.cern.ch/lhcb-scifi/SciFiMatG4_v2
> $ cd SciFiMatG4_v2
> $ git clone https://gitlab.cern.ch/lphe-lab-tools/SciFiMatG4.git
> $ cd SciFiMatG4
Edit the environment script setup.sh to point to your Geant and ROOT installation and `source` it to add it to your environment:
> $ source setup.sh //(required before running the software too if you open a new terminal.)
Now to cun cmake in the build folder:
Now to run cmake in the build folder:
> mkdir -p SciFiSim-build //it is safe to delete this folder if you have problems during compilation
> $ mkdir -p SciFiSim-build //it is safe to delete this folder if you have problems during compilation
> cd SciFiSim-build
> $ cd SciFiSim-build
> cmake ../SciFiSim
> cmake --build . // or make clean; make
> $ cmake ../SciFiSim
> $ cmake --build . // or make clean; make // or just make
To run the Geant4 simulation in batch mode:
> cd ../SimulationData
> $ cd SimulationData/test/ //(or any other subdirectory of /SimulationData that contains a .mac file)
> ../SciFiSim-build/scifiSim muongun.mac
> $ ../../SciFiSim-build/scifiSim test.mac
To run the simulation from within the GUI (for visualization only) very slow if optical photon physics are enabled:
> $ ../SciFiSim-build/scifiSim //opens the gui
> /control/execute/ vis.mac
> $ cd SimulationData/visualisation/
> $ ../../SciFiSim-build/scifiSim //opens the gui
Then in the Session box, type:
> $ /control/execute/ vis.mac
......@@ -64,7 +70,7 @@ There are 5 different TTrees created.
- InitialParticle(Not sure what the difference is)
- EnergyTrack (track length in fibres and energy deposited)
## SciFiSimG4_v2 requires installation of Geant4 first
## SciFiSimG4 requires installation of Geant4 first
> $ wget http://geant4.web.cern.ch/geant4/support/source/geant4.10.07.p01.tar.gz
......
SciFiSim/CMakeLists.txt
View file @ 07fea035
......@@ -2,7 +2,7 @@
#----------------------------------------------------------------------------
# Setup the project
cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
cmake_minimum_required(VERSION 3.7 FATAL_ERROR)
project(scifiSim)
#----------------------------------------------------------------------------
......@@ -24,12 +24,20 @@ endif()
include(${Geant4_USE_FILE})
include_directories(${PROJECT_SOURCE_DIR}/include)
#include("FindROOT.cmake")
set(CMAKE_MODULE_PATH$ENV{ROOTSYS}/etc/cmake)
find_package(ROOT REQUIRED)
list(APPEND libs $ENV{ROOTSYS}/lib)
include_directories($ENV{ROOTSYS}/include)
message(status $ENV{ROOTSYS}/include)
#----------------------------------------------------------------------------
# ROOT setup
#----------------------------------------------------------------------------
# If ROOT is installed in a non-standard prefix, set CMAKE_PREFIX_PATH when calling cmake:
# cmake .. -DCMAKE_PREFIX_PATH="/opt/homebrew"
find_package(ROOT REQUIRED COMPONENTS Core RIO Tree)
message(STATUS "ROOT include dirs: ${ROOT_INCLUDE_DIRS}")
message(STATUS "ROOT libraries: ${ROOT_LIBRARIES}")
find_package(nlohmann_json REQUIRED)
#----------------------------------------------------------------------------
......@@ -43,7 +51,8 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
# Add the executable, and link it to the Geant4 libraries
#
add_executable(scifiSim scifiSim.cc ${sources} ${headers})
target_link_libraries(scifiSim ${Geant4_LIBRARIES} ${ROOT_LIBRARIES})
target_include_directories(scifiSim PRIVATE ${ROOT_INCLUDE_DIRS})
target_link_libraries(scifiSim ${Geant4_LIBRARIES} ${ROOT_LIBRARIES} nlohmann_json::nlohmann_json)
#----------------------------------------------------------------------------
# Copy all scripts to the build directory, i.e. the directory in which we
......
SciFiSim/include/Analysis.hh
View file @ 07fea035
......@@ -35,10 +35,12 @@ public:
void EndOfRun();
void Close();
void LoadPDETable();
double GetPDE(double wavelength);
// Adds a new row to the tree : detectedPhotons
void FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb, Double_t xPixel,
Double_t yPixel, Double_t energy,
void FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb,
Double_t xPixel, Double_t yPixel, Double_t energy,
Double_t time, Double_t length, Double_t absTime,
Double_t x, Double_t y, Double_t z,
Double_t px, Double_t py, Double_t pz,
......@@ -53,7 +55,8 @@ public:
// Adds a new to to the tree : initialParticle
void FillInitialParticle(Double_t runID, Double_t eventID,
Double_t energy, Double_t xMom, Double_t yMom, Double_t zMom);
Double_t energy, Double_t xPos, Double_t yPos, Double_t zPos, // Modified to output Pos information
Double_t xMom, Double_t yMom, Double_t zMom);
// Adds a new to to the tree : primaryParticleTrack
void FillPrimaryParticleTrack(Double_t runID, Double_t eventID, Double_t xPos, Double_t yPos, Double_t zPos);
......@@ -99,7 +102,7 @@ private:
TFile* dataFile;
char fileName[40];
char fileName[256];
TH1F * h_energy;
TH1F * h_trackL;
......@@ -110,10 +113,13 @@ private:
TTree* initialParticle;
TTree* energyTrack;
std::vector<std::pair<double, double>> pdeTable;
bool pdeLoaded = false;
Float_t runIDBuffer;
Float_t eventIDBuffer;
Float_t detNumbBuffer;
Float_t xPixelBuffer;
Float_t yPixelBuffer;
Float_t xPosBuffer;
......@@ -121,6 +127,8 @@ private:
Float_t zPosBuffer;
Float_t energyBuffer;
Float_t wavelengthBuffer;
Float_t detProbBuffer;
Bool_t ifDetBuffer;
Float_t timeBuffer;
Float_t lengthBuffer;
Float_t xMomBuffer;
......
SciFiSim/include/DetectorConstruction.hh
View file @ 07fea035
......@@ -33,7 +33,7 @@ private:
/* ++ methods ++ */
void DefineMaterials();
void DefineMaterialProperties();
void ConstructFiber(); // Constructs and places the fibres
void ConstructFibre(G4double, G4LogicalVolume*); // Constructs and places the fibres
void ConstructFiberSheet(); // Constructs and places the fibres
G4ThreeVector objectPos(G4int, G4int); // for positioning of multiple fibres
G4ThreeVector objectPos(G4int, G4int, G4double);// for positioning of detector
......@@ -59,6 +59,7 @@ private:
G4Material *PMMA;
G4Material *PMMA2;
G4Material *Epoxy;
G4Material *PolyimideBlack;
G4Material *Glue;
G4Material *TiO2;
G4Material *Abs_plastic;
......@@ -81,8 +82,8 @@ private:
// epoxybox
G4LogicalVolume* epoxyLog;
G4VPhysicalVolume* epoxyPhy;
// G4LogicalVolume* fibreContainerLog;
// G4VPhysicalVolume* epoxyPhy;
// abs plastic
G4LogicalVolume* absLog;
......
SciFiSim/include/Parameters.hh
View file @ 07fea035
......@@ -36,6 +36,16 @@ public:
G4double RandomSeed();
G4double RandomNumber();
G4double FibreLength();
G4double AirGapZ1();
G4double AirGapL1();
G4double AirGapZ2();
G4double AirGapL2();
G4int NLayers();
G4bool IfMIEHG();
G4double MIEHG_FORWARD_RATIO();
G4double MIEHG_FORWARD();
G4double MIEHG_BACKWARD();
std::string SipmDetectionEfficiencyFileName();
G4double SemiAxisZ();
G4double SemiAxisY();
G4double ProbabilityOfPhotonLossAtSurface();
......@@ -56,6 +66,14 @@ public:
G4double* Intensity;
G4int NumberOfInterpolatedPoints();
char* MieScatteringFileName();
G4int NumberOfMieEnergies();
G4double* MieEnergy;
G4double* MieMFP;
G4double* MieForward;
G4double* MieBackward;
G4double* MieG;
char* WlsAbsSpectrumFileName();
G4int NumberOfWlsAbsEnergies();
G4double* WlsAbsEnergy;
......@@ -110,6 +128,16 @@ private:
G4double randomSeed;
G4double randomNumber;
G4double fibreLength;
G4double airGapZ1;
G4double airGapL1;
G4double airGapZ2;
G4double airGapL2;
G4int nLayers;
G4bool ifMIEHG;
G4double mIEHG_FORWARD_RATIO;
G4double mIEHG_FORWARD;
G4double mIEHG_BACKWARD;
std::string sipmDetectionEfficiencyFileName;
G4double semiAxisZ;
G4double semiAxisY;
G4double probabilityOfPhotonLossAtSurface;
......@@ -130,6 +158,9 @@ private:
// G4double[] intensity;
G4int numberOfInterpolatedPoints;
char mieScatteringFileName[256];
G4int numberOfMieEnergies;
char wlsAbsSpectrumFileName[256];
G4int numberOfWlsAbsEnergies;
......
SciFiSim/python/CalculateTiO2EpoxyMIE.py 0 → 100644
View file @ 07fea035
import numpy as np
import miepython
import pandas as pd
results = []
df = pd.read_csv("../../SimulationData/parameterFiles/TPBD_emission.csv", sep=r"\s+", header=None)
df.columns = ["Energy_eV", "OtherValue"]
print(df["Energy_eV"])
df["Wavelength_nm"] = 1240.0 / df["Energy_eV"]
# radius_nm = 150 # TiO₂ particle radius
# n_particle = 2.7 # TiO₂ index (rutile)
# n_host = 1.5 # epoxy index
# Constants for geometry
radius_nm = 150 # TiO2 particle radius: 150 nm
volume_particle_cm3 = (4/3) * np.pi * (radius_nm*1e-7)**3
density_tio2 = 4.23 # g/cm³, from internet
particle_mass = density_tio2 * volume_particle_cm3 # g, mass of individual particles (each particle can contain large number of molecules)
density_epoxy = 1.347 # g/cm³
mass_fraction_tio2 = 0.2 # fraction of TiO₂ in epoxy by mass
mass_density_tio2 = mass_fraction_tio2 * density_epoxy
number_density_tio2_cm3 = mass_density_tio2 / particle_mass
number_density_tio2_m3 = number_density_tio2_cm3 * 1e6
# print (number_density*1e-12)
ref_index_tio2 = 2.7
ref_index_glue = 1.59
# Output storage
mfp_m = [] # mean free path in m, corresponding to MIEHG value in Geant4
forward_ratio = []
backward_ratio = []
g_values = []
for wl_nm in df["Wavelength_nm"]:
d = 2 * radius_nm
m = ref_index_tio2 / ref_index_glue # relative index: TiO₂ / glue (https://miepython.readthedocs.io/en/latest/01_basics.html, equation 3)
# Since material thin enough, we can ignore attenuation (light will not be lost inside the epoxy). Therefore ignore the imaginary part.
qext, qsca, qback, g = miepython.efficiencies(m, d, wl_nm) # extinction efficiency, scattering efficiency, backscatter efficiency, scattering anisotropy
# not using qext
forward_ratio.append(1 - qback/qsca)
backward_ratio.append(qback/qsca)
mfp_value_m = 1 / (number_density_tio2_m3 * np.pi * (radius_nm*1e-9)**2 * qsca)
mfp_m.append(mfp_value_m)
g_values.append(g)
df_out = pd.DataFrame({
"Energy_eV": df["Energy_eV"],
"mfp_m": mfp_m,
"forward_ratio": forward_ratio,
"backward_ratio": backward_ratio,
"g_values": g_values,
})
df_out.to_csv("mie_results.csv", index=False, sep=" ")
print(df_out)
SciFiSim/python/ConvertToJson.py 0 → 100644
View file @ 07fea035
import json
import argparse
def dat_to_json(dat_path, json_path):
with open(dat_path, "r") as f:
lines = [line.strip() for line in f if line.strip()]
param_dict = {}
i = 0
try:
param_dict["randomSeed"] = float(lines[i]); i += 1
param_dict["randomNumber"] = float(lines[i]); i += 1
param_dict["fibreLength"] = float(lines[i]); i += 1
param_dict["semiAxisZ"] = float(lines[i]); i += 1
param_dict["semiAxisY"] = float(lines[i]); i += 1
param_dict["triggerX"] = float(lines[i]); i += 1
param_dict["triggerY"] = float(lines[i]); i += 1
param_dict["triggerZ"] = float(lines[i]); i += 1
param_dict["triggerXPos"] = float(lines[i]); i += 1
param_dict["triggerZPos"] = float(lines[i]); i += 1
param_dict["probabilityOfPhotonLossAtSurface"] = float(lines[i]); i += 1
param_dict["placeMirror"] = bool(int(lines[i])); i += 1
param_dict["mirrorReflectivity"] = float(lines[i]); i += 1
param_dict["detectorMaterial"] = bool(int(lines[i])); i += 1
param_dict["emissionSpectrumFileName"] = lines[i]; i += 1
param_dict["numberOfInterpolatedPoints"] = int(lines[i]); i += 1
param_dict["wlsAbsSpectrumFileName"] = lines[i]; i += 1
param_dict["wlsEmissionSpectrumFileName"] = lines[i]; i += 1
param_dict["scintillationYield"] = float(lines[i]); i += 1
param_dict["resolutionScale"] = float(lines[i]); i += 1
param_dict["decayTimeFast"] = float(lines[i]); i += 1
param_dict["decayTimeSlow"] = float(lines[i]); i += 1
param_dict["yieldRatio"] = float(lines[i]); i += 1
param_dict["birksConstant"] = float(lines[i]); i += 1
param_dict["wlsDecayTime"] = float(lines[i]); i += 1
param_dict["refractiveIndexVacuum"] = lines[i]; i += 1
param_dict["refractiveIndexCore"] = lines[i]; i += 1
param_dict["refractiveIndexClad1"] = lines[i]; i += 1
param_dict["refractiveIndexClad2"] = lines[i]; i += 1
param_dict["absorptionCore"] = lines[i]; i += 1
param_dict["absorptionClad1"] = lines[i]; i += 1
param_dict["absorptionClad2"] = lines[i]; i += 1
param_dict["absorptionFromIrradiationCore"] = lines[i]; i += 1
param_dict["absorptionFromIrradiationClad1"] = lines[i]; i += 1
param_dict["absorptionFromIrradiationClad2"] = lines[i]; i += 1
param_dict["sectionsFileName"] = lines[i]; i += 1
param_dict["rayleighCore"] = lines[i]; i += 1
param_dict["rayleighClad1"] = lines[i]; i += 1
param_dict["rayleighClad2"] = lines[i]; i += 1
except IndexError:
print("Error: .dat file appears to be incomplete or malformed.")
return
with open(json_path, "w") as f:
json.dump(param_dict, f, indent=4)
print(f"Converted to JSON and saved to: {json_path}")
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Convert .dat parameter file to JSON.")
parser.add_argument("input", help="Path to the input .dat file")
parser.add_argument("output", help="Path to the output .json file")
args = parser.parse_args()
dat_to_json(args.input, args.output)
SciFiSim/python/mie_results.csv 0 → 100644
View file @ 07fea035
Energy_eV mfp_m forward_ratio backward_ratio g_values
1.9 2.2460430870153308e-06 0.8933484240653992 0.10665157593460084 0.530794021100411
1.922 2.1514546541869238e-06 0.9110701230038827 0.08892987699611726 0.5413462052073424
2.046 1.6821091070319383e-06 0.9655116988039657 0.034488301196034306 0.582352897241809
2.17 1.347506954097682e-06 0.9469358966904547 0.05306410330954524 0.5884489355956701
2.294 1.1608421298397205e-06 0.8913354482450572 0.10866455175494284 0.5717611512067021
2.418 1.0853427044436351e-06 0.8376634341476806 0.16233656585231934 0.5526453749501974
2.542 1.0636579754251157e-06 0.8018418094844373 0.19815819051556272 0.544293511540108
2.591 1.0589042946323193e-06 0.7927559846738559 0.20724401532614406 0.5450622275747696
2.641 1.0526255155471e-06 0.7861961169674857 0.21380388303251427 0.548274150527059
2.678 1.04593533850466e-06 0.7830954345739758 0.21690456542602427 0.5521789908036471
2.715 1.0367955992384133e-06 0.7815655378970024 0.21843446210299758 0.5573217169679914
2.765 1.0196614733983262e-06 0.7822630780041915 0.2177369219958085 0.5661089375890631
2.814 9.967629020827126e-07 0.786484187348927 0.21351581265107297 0.5765678273186746
2.914 9.304195721185947e-07 0.8077541470332248 0.19224585296677518 0.6022554604621921
3.038 8.235897878384737e-07 0.8536671892978943 0.14633281070210571 0.6364630684535498
3.1 7.725069295226393e-07 0.876151509226725 0.12384849077327499 0.651096427459337
3.162 7.318834931130283e-07 0.8906479507209527 0.1093520492790473 0.66224223863101
3.224 7.056431800145182e-07 0.8930801736482404 0.10691982635175955 0.6692979933572701
3.286 6.93944772684968e-07 0.8826661102539248 0.11733388974607523 0.6723548254590311
3.4 7.010389549125751e-07 0.8374085910543141 0.16259140894568586 0.6690912442098049
3.5 7.236809589045797e-07 0.7824829682671686 0.2175170317328314 0.6587267936629412
3.6 7.495085113781051e-07 0.7229059784941927 0.2770940215058073 0.6435259460902165
3.7 7.687717102409724e-07 0.6610266026165954 0.33897339738340454 0.6263889084099602
3.8 7.712025325626651e-07 0.5930689353880929 0.40693106461190703 0.6113082727942503
3.9 7.472942659917688e-07 0.5175831549180008 0.48241684508199917 0.6035925457665996
4.0 6.989474938482231e-07 0.4578126759243978 0.5421873240756022 0.6089759629233108
4.1 6.539925121086237e-07 0.4630000164839789 0.5369999835160211 0.6271817505581787
4.2 6.438481635059883e-07 0.5357346132327996 0.4642653867672004 0.6473111509136431
4.3 6.665966601529686e-07 0.6128280966395372 0.38717190336046275 0.6596121912170165
4.4275 7.174726042432396e-07 0.6450591605400392 0.3549408394599608 0.6614169772701893
SciFiSim/python/setup.sh 0 → 100644
View file @ 07fea035
source /Users/isiral/IsmetPython/bin/activate
SciFiSim/src/Analysis.cc
View file @ 07fea035
// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler
// modified by Bastian Rössler, Lai Gui
#include "Analysis.hh"
#include "Parameters.hh"
#include "G4Threading.hh"
#include "Randomize.hh"
#include "G4GeneralParticleSource.hh"
#include <ctime>
#include <fstream>
#include <filesystem>
#include <string>
#include <sstream>
#include <iostream>
#include <cmath>
#include <regex>
#include <cstdlib>
Analysis* Analysis::singleton = 0;
Analysis::Analysis()
{
// Set file name for simulation output
G4int threadNumber;
std::ifstream inFile;
inFile.open("threadNumber.txt");
inFile >> threadNumber;
inFile.close();
// Create output filename
const char* jobIdEnv = std::getenv("SLURM_JOB_ID");
std::string jobId = jobIdEnv ? jobIdEnv : "unknown";
sprintf(fileName, "outFile_%s.root", jobId.c_str());
std::cout << "Output will be written to: " << fileName << std::endl;
sprintf(fileName,"outFile_%i.root", threadNumber);
////////////////////////////////////////////////////////////////////////////////////////
detectedPhotons = new TTree("DetectedPhotons","Photons detected at the detector strip.");
detectedPhotons->Branch("runID", &runIDBuffer, "runID/F"); //modifed by me
detectedPhotons->Branch("eventID", &eventIDBuffer, "eventID/F"); //modifed by me
detectedPhotons->Branch("detNumb", &detNumbBuffer, "detNumb/F"); //modifed by me
detectedPhotons->Branch("runID", &runIDBuffer, "runID/F");
detectedPhotons->Branch("eventID", &eventIDBuffer, "eventID/F");
detectedPhotons->Branch("detNumb", &detNumbBuffer, "detNumb/F");
detectedPhotons->Branch("xPixel", &xPixelBuffer, "xPixel/F"); //modifed by me
detectedPhotons->Branch("yPixel", &yPixelBuffer, "yPixel/F"); //modifed by me
detectedPhotons->Branch("xPixel", &xPixelBuffer, "xPixel/F");
detectedPhotons->Branch("yPixel", &yPixelBuffer, "yPixel/F");
detectedPhotons->Branch("energy",&energyBuffer,"energy/F");
detectedPhotons->Branch("wavelength",&wavelengthBuffer,"wavelength/F");
detectedPhotons->Branch("detProb", &detProbBuffer, "detProb/F");
detectedPhotons->Branch("ifDet", &ifDetBuffer, "ifDet/O");
detectedPhotons->Branch("localtime",&timeBuffer,"localtime/F");
detectedPhotons->Branch("abstime",&absTimeBuffer,"abstime/F");
detectedPhotons->Branch("length",&lengthBuffer,"length/F");
......@@ -102,11 +115,14 @@ Analysis::Analysis()
initialParticle->Branch("runID", &runIDBuffer, "runID/F");
initialParticle->Branch("eventID", &eventIDBuffer, "eventID/F");
initialParticle->Branch("energy", &energyBuffer, "energy/F");
initialParticle->Branch("xPos", &xPosBuffer, "xPos/F");
initialParticle->Branch("yPos", &yPosBuffer, "yPos/F");
initialParticle->Branch("zPos", &zPosBuffer, "zPos/F");
initialParticle->Branch("xMom", &xMomBuffer, "xMom/F");
initialParticle->Branch("yMom", &yMomBuffer, "yMom/F");
initialParticle->Branch("zMom", &zMomBuffer, "zMom/F");
}
......@@ -168,13 +184,18 @@ void Analysis::PrepareNewEvent(const G4Event* anEvent)
lengthInOuterCladding.resize(0);
}
void Analysis::FillInitialParticle(Double_t runID, Double_t eventID, Double_t energy, Double_t xMom, Double_t yMom, Double_t zMom)
void Analysis::FillInitialParticle(Double_t runID, Double_t eventID, Double_t energy, Double_t xPos, Double_t yPos, Double_t zPos, Double_t xMom, Double_t yMom, Double_t zMom)
{
if(initialParticle != NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
energyBuffer = energy;
xPosBuffer = xPos; // Modified to output Pos information
yPosBuffer = yPos;
zPosBuffer = zPos;
xMomBuffer = xMom;
yMomBuffer = yMom;
zMomBuffer = zMom;
......@@ -239,8 +260,53 @@ void Analysis::FillPrimaryParticleTrack(Double_t runID, Double_t eventID,
}
}
void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb, Double_t xPixel,
Double_t yPixel, Double_t energy,
void Analysis::LoadPDETable() {
if (pdeLoaded) return;
std::ifstream infile(Parameters::GetInstance()->SipmDetectionEfficiencyFileName()); // or your actual filename
if (!infile.is_open()) {
std::cerr << "Failed to open PDE table!" << std::endl;
return;
}
std::string line;
std::getline(infile, line); // skip header line
double lambda, pde;
while (std::getline(infile, line)) {
std::istringstream iss(line);
if (iss >> lambda >> pde) {
pdeTable.emplace_back(lambda, pde);
}
}
pdeLoaded = true;
}
double Analysis::GetPDE(double wavelength_nm) {
LoadPDETable();
if (pdeTable.empty()) return 0.0;
double min_diff = std::numeric_limits<double>::max();
double closest_pde = 0.0;
for (const auto& entry : pdeTable) {
double diff = std::abs(entry.first - wavelength_nm);
if (diff < min_diff) {
min_diff = diff;
closest_pde = entry.second;
}
}
return closest_pde;
}
void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb,
Double_t xPixel, Double_t yPixel, Double_t energy,
Double_t time, Double_t length, Double_t absTime,
Double_t x, Double_t y, Double_t z,
Double_t px, Double_t py, Double_t pz,
......@@ -254,6 +320,7 @@ void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t de
runIDBuffer = runID;
eventIDBuffer = eventID;
detNumbBuffer = detNumb;
xPixelBuffer = xPixel;
yPixelBuffer = yPixel;
......@@ -262,6 +329,8 @@ void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t de
lengthBuffer = length;
wavelengthBuffer = 1239.842/energy;
detProbBuffer = GetPDE(wavelengthBuffer) / 100.0;
ifDetBuffer = (G4UniformRand() < detProbBuffer);
absTimeBuffer = absTime;
xBuffer = x;
......
SciFiSim/src/Analysis.cc~ deleted 100644 → 0
View file @ 6f8f2267
// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler
#include "Analysis.hh"
#include "Parameters.hh"
#include <ctime>
#include <fstream>
Analysis* Analysis::singleton = 0;
Analysis::Analysis()
{
// Set file name for simulation output
G4int threadNumber;
std::ifstream inFile;
inFile.open("threadNumber.txt");
inFile >> threadNumber;
inFile.close();
sprintf(fileName,"outFile_%i.root", threadNumber);
detectedPhotons = new TTree("DetectedPhotons","Photons detected at the detector strip.");
detectedPhotons->Branch("runID", &runIDBuffer, "runID/F"); //modifed by me
detectedPhotons->Branch("eventID", &eventIDBuffer, "eventID/F"); //modifed by me
detectedPhotons->Branch("detNumb", &detNumbBuffer, "detNumb/F"); //modifed by me
detectedPhotons->Branch("xPixel", &xPixelBuffer, "xPixel/F"); //modifed by me
detectedPhotons->Branch("yPixel", &yPixelBuffer, "yPixel/F"); //modifed by me
detectedPhotons->Branch("energy",&energyBuffer,"energy/F");
detectedPhotons->Branch("wavelength",&wavelengthBuffer,"wavelength/F");
detectedPhotons->Branch("localtime",&timeBuffer,"localtime/F");
detectedPhotons->Branch("abstime",&absTimeBuffer,"abstime/F");
detectedPhotons->Branch("length",&lengthBuffer,"length/F");
detectedPhotons->Branch("x",&xBuffer,"x/F");
detectedPhotons->Branch("y",&yBuffer,"y/F");
detectedPhotons->Branch("z",&zBuffer,"z/F");
detectedPhotons->Branch("px",&pxBuffer,"px/F");
detectedPhotons->Branch("py",&pyBuffer,"py/F");
detectedPhotons->Branch("pz",&pzBuffer,"pz/F");
detectedPhotons->Branch("vertexX",&vertexXBuffer,"vertexX/F");
detectedPhotons->Branch("vertexY",&vertexYBuffer,"vertexY/F");
detectedPhotons->Branch("vertexZ",&vertexZBuffer,"vertexZ/F");
detectedPhotons->Branch("vertexPx",&vertexPxBuffer,"vertexPx/F");
detectedPhotons->Branch("vertexPy",&vertexPyBuffer,"vertexPy/F");
detectedPhotons->Branch("vertexPz",&vertexPzBuffer,"vertexPz/F");
detectedPhotons->Branch("gpsPosX",&gpsPositionX,"gpsPosX/F");
detectedPhotons->Branch("gpsPosY",&gpsPositionY,"gpsPosY/F");
detectedPhotons->Branch("gpsPosZ",&gpsPositionZ,"gpsPosZ/F");
detectedPhotons->Branch("gpsDirX",&gpsDirectionX,"gpsPosX/F");
detectedPhotons->Branch("gpsDirY",&gpsDirectionY,"gpsPosY/F");
detectedPhotons->Branch("gpsDirZ",&gpsDirectionZ,"gpsPosZ/F");
detectedPhotons->Branch("runId",&runIdBuffer,"runId/I");
detectedPhotons->Branch("eventId",&eventIdBuffer,"eventId/I");
detectedPhotons->Branch("trackId",&trackIdBuffer,"trackId/I");
detectedPhotons->Branch("creatorProcess",&creatorProcessBuffer,"creatorProcess/I");
detectedPhotons->Branch("parentId",&parentIdBuffer,"parentId/I");
detectedPhotons->Branch("reflMirr",&reflMirrBuffer,"reflMirr/I");
detectedPhotons->Branch("reflSurf",&reflSurfBuffer,"reflSurf/I");
detectedPhotons->Branch("reflTotalCladClad",&reflTotalCladCladBuffer,"reflTotalCladClad/I");
detectedPhotons->Branch("reflTotalCoreClad",&reflTotalCoreCladBuffer,"reflTotalCoreClad/I");
detectedPhotons->Branch("reflFresnelCladClad",&reflFresnelCladCladBuffer,"reflFresnelCladClad/I");
detectedPhotons->Branch("reflFresnelCoreClad",&reflFresnelCoreCladBuffer,"reflFresnelCoreClad/I");
detectedPhotons->Branch("refracCladClad",&refracCladCladBuffer,"refracCladClad/I");
detectedPhotons->Branch("refracCoreClad",&refracCoreCladBuffer,"refracCoreClad/I");
detectedPhotons->Branch("rayleighScatterings",&rayleighScatteringsBuffer,"rayleighScatterings/I");
detectedPhotons->Branch("lengthInCore",&lengthInCoreBuffer,"lengthInCore/F");
detectedPhotons->Branch("lengthInInnerCladding",&lengthInInnerCladdingBuffer,"lengthInInnerCladding/F");
detectedPhotons->Branch("lengthInOuterCladding",&lengthInOuterCladdingBuffer,"lengthInOuterCladding/F");
trigger = new TTree("Trigger","Hits in the trigger.");
trigger->Branch("runID", &runIDBuffer, "runID/F");
trigger->Branch("eventID", &eventIDBuffer, "eventID/F");
trigger->Branch("edep", &edepBuffer, "edep/F");
trigger->Branch("xPos", &xPosBuffer, "xPos/F");
trigger->Branch("yPos", &yPosBuffer, "yPos/F");
trigger->Branch("zPos", &zPosBuffer, "zPos/F");
energyTrack = new TTree("EnergyTrack", "Deposited energy and tracklength in core of fibre.");
energyTrack->Branch("runID", &runIDBuffer, "runID/F");
energyTrack->Branch("eventID", &eventIDBuffer, "eventID/F");
energyTrack->Branch("edep", &energyBuffer, "edep/F");
energyTrack->Branch("trackL", &lengthBuffer, "trackL/F");
h_energy = new TH1F("h_energy","h_energy",200,0.0,1.0);
h_trackL = new TH1F("h_trackL","h_trackL",400,0.0,2.0);
primaryParticleTrack = new TTree("PrimaryParticleTrack", "Track of the primary particles");
primaryParticleTrack->Branch("runID", &runIDBuffer, "runID/F");
primaryParticleTrack->Branch("eventID", &eventIDBuffer, "eventID/F");
primaryParticleTrack->Branch("xPos", &xPosBuffer, "xPos/F");
primaryParticleTrack->Branch("yPos", &yPosBuffer, "yPos/F");
primaryParticleTrack->Branch("zPos", &zPosBuffer, "zPos/F");
initialParticle = new TTree("InitialParticle", "Initial Particle");
initialParticle->Branch("runID", &runIDBuffer, "runID/F");
initialParticle->Branch("eventID", &eventIDBuffer, "eventID/F");
initialParticle->Branch("energy", &energyBuffer, "energy/F");
initialParticle->Branch("xMom", &xMomBuffer, "xMom/F");
initialParticle->Branch("yMom", &yMomBuffer, "yMom/F");
initialParticle->Branch("zMom", &zMomBuffer, "zMom/F");
}
void Analysis::PrepareNewRun(const G4Run* aRun)
{
runIdBuffer = aRun->GetRunID();
}
void Analysis::EndOfRun()
{
dataFile = new TFile(fileName,"update");
//detectedPhotons->Write("",TObject::kOverwrite);
//trigger->Write("", TObject::kOverwrite);
// primaryParticleTrack->Write("", TObject::kOverwrite);
//initialParticle->Write("", TObject::kOverwrite);
// energyTrack->Write("", TObject::kOverwrite);
h_energy->Write("", TObject::kOverwrite);
h_trackL->Write("", TObject::kOverwrite);
G4cout << "Data written to file " << fileName << G4endl;
dataFile->Close();
G4cout << "Data file closed: " << fileName << G4endl;
delete dataFile;
}
void Analysis::Close()
{
delete detectedPhotons;
delete trigger;
delete initialParticle;
delete energyTrack;
}
void Analysis::PrepareNewEvent(const G4Event* anEvent)
{
eventIdBuffer = anEvent->GetEventID();
reflectionsAtMirror.resize(0);
reflectionsAtFibreSurface.resize(0);
reflectionsTotalAtCladCladInterface.resize(0);
reflectionsTotalAtCoreCladInterface.resize(0);
reflectionsFresnelAtCladCladInterface.resize(0);
reflectionsFresnelAtCoreCladInterface.resize(0);
refractionsAtCladCladInterface.resize(0);
refractionsAtCoreCladInterface.resize(0);
rayleighScatterings.resize(0);
lengthInCore.resize(0);
lengthInInnerCladding.resize(0);
lengthInOuterCladding.resize(0);
}
void Analysis::FillInitialParticle(Double_t runID, Double_t eventID, Double_t energy, Double_t xMom, Double_t yMom, Double_t zMom)
{
if(initialParticle != NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
energyBuffer = energy;
xMomBuffer = xMom;
yMomBuffer = yMom;
zMomBuffer = zMom;
// initialParticle->Fill();
}
else
G4cout << "Nullpointer to initialParticle tree." << G4endl;
}
void Analysis::FillTrigger(Double_t runID, Double_t eventID, Double_t edep, Double_t xPos, Double_t yPos, Double_t zPos)
{
if(trigger != NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
edepBuffer = edep;
xPosBuffer = xPos;
yPosBuffer = yPos;
zPosBuffer = zPos;
// trigger->Fill();
}
else
{
G4cout << "Nullpointer to trigger tree." << G4endl;
}
}
void Analysis::FillEnergyTrack(Double_t runID, Double_t eventID, Double_t energy, Double_t trackL)
{
if(energyTrack != NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
energyBuffer = energy;
lengthBuffer = trackL;
// energyTrack->Fill();
h_energy->Fill(energy);
h_trackL->Fill(trackL);
}
else
{
G4cout << "Nullpointer to energyTrack tree" << G4endl;
}
}
void Analysis::FillPrimaryParticleTrack(Double_t runID, Double_t eventID,
Double_t xPos, Double_t yPos, Double_t zPos)
{
if(primaryParticleTrack != NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
xPosBuffer = xPos;
yPosBuffer = yPos;
zPosBuffer = zPos;
// primaryParticleTrack->Fill();
}
else
{
G4cout << "Nullpointer to primaryParticleTrack tree" << G4endl;
}
}
void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb, Double_t xPixel,
Double_t yPixel, Double_t energy,
Double_t time, Double_t length, Double_t absTime,
Double_t x, Double_t y, Double_t z,
Double_t px, Double_t py, Double_t pz,
Double_t vertexX, Double_t vertexY, Double_t vertexZ,
Double_t vertexPx, Double_t vertexPy, Double_t vertexPz,
Int_t trackId,
Int_t creatorProcess, Int_t parentId)
{
if(detectedPhotons!=NULL)
{
runIDBuffer = runID;
eventIDBuffer = eventID;
detNumbBuffer = detNumb;
xPixelBuffer = xPixel;
yPixelBuffer = yPixel;
energyBuffer = energy;
timeBuffer = time;
lengthBuffer = length;
wavelengthBuffer = 1239.842/energy;
absTimeBuffer = absTime;
xBuffer = x;
yBuffer = y;
zBuffer = z;
pxBuffer = px;
pyBuffer = py;
pzBuffer = pz;
vertexXBuffer = vertexX;
vertexYBuffer = vertexY;
vertexZBuffer = vertexZ;
vertexPxBuffer = vertexPx;
vertexPyBuffer = vertexPy;
vertexPzBuffer = vertexPz;
trackIdBuffer = trackId;
creatorProcessBuffer = creatorProcess;
parentIdBuffer = parentId;
reflMirrBuffer = reflectionsAtMirror[trackId];
reflSurfBuffer = reflectionsAtFibreSurface[trackId];
reflTotalCladCladBuffer = reflectionsTotalAtCladCladInterface[trackId];
reflTotalCoreCladBuffer = reflectionsTotalAtCoreCladInterface[trackId];
reflFresnelCladCladBuffer = reflectionsFresnelAtCladCladInterface[trackId];
reflFresnelCoreCladBuffer = reflectionsFresnelAtCoreCladInterface[trackId];
refracCladCladBuffer = refractionsAtCladCladInterface[trackId];
refracCoreCladBuffer = refractionsAtCoreCladInterface[trackId];
rayleighScatteringsBuffer = rayleighScatterings[trackId];
lengthInCoreBuffer = lengthInCore[trackId];
lengthInInnerCladdingBuffer = lengthInInnerCladding[trackId];
lengthInOuterCladdingBuffer = lengthInOuterCladding[trackId];
// detectedPhotons->Fill();
}
else
{
G4cout << G4endl << "Pointer to detectedPhotons tree is NULL!" << G4endl;
}
}
char* Analysis::FileName()
{
return fileName;
}
void Analysis::SetGpsPosition(G4ThreeVector position)
{
gpsPositionX = position[0];
gpsPositionY = position[1];
gpsPositionZ = position[2];
}
void Analysis::IncreaseReflectionsAtMirror(Int_t trackId)
{
reflectionsAtMirror[trackId]++;
}
void Analysis::SetGpsDirection(G4ThreeVector direction)
{
gpsDirectionX = direction[0];
gpsDirectionY = direction[1];
gpsDirectionZ = direction[2];
}
void Analysis::IncreaseReflectionsAtFibreSurface(Int_t trackId)
{
reflectionsAtFibreSurface[trackId]++;
}
void Analysis::IncreaseTotalReflectionsAtCladCladInterface(Int_t trackId)
{
reflectionsTotalAtCladCladInterface[trackId]++;
}
void Analysis::IncreaseTotalReflectionsAtCoreCladInterface(Int_t trackId)
{
reflectionsTotalAtCoreCladInterface[trackId]++;
}
void Analysis::IncreaseFresnelReflectionsAtCladCladInterface(Int_t trackId)
{
reflectionsFresnelAtCladCladInterface[trackId]++;
}
void Analysis::IncreaseFresnelReflectionsAtCoreCladInterface(Int_t trackId)
{
reflectionsFresnelAtCoreCladInterface[trackId]++;
}
void Analysis::IncreaseRefractionsAtCladCladInterface(Int_t trackId)
{
refractionsAtCladCladInterface[trackId]++;
}
void Analysis::IncreaseRefractionsAtCoreCladInterface(Int_t trackId)
{
refractionsAtCoreCladInterface[trackId]++;
}
void Analysis::IncreaseRayleighScatterings(Int_t trackId)
{
rayleighScatterings[trackId]++;
}
void Analysis::IncreaseReflectionRefractionAndScatteringVectors(Int_t trackId)
{
for(int i=reflectionsAtMirror.size(); i<trackId+1; i++)
reflectionsAtMirror.push_back(0);
for(int i=reflectionsAtFibreSurface.size(); i<trackId+1; i++)
reflectionsAtFibreSurface.push_back(0);
for(int i=reflectionsTotalAtCladCladInterface.size(); i<trackId+1; i++)
reflectionsTotalAtCladCladInterface.push_back(0);
for(int i=reflectionsTotalAtCoreCladInterface.size(); i<trackId+1; i++)
reflectionsTotalAtCoreCladInterface.push_back(0);
for(int i=reflectionsFresnelAtCladCladInterface.size(); i<trackId+1; i++)
reflectionsFresnelAtCladCladInterface.push_back(0);
for(int i=reflectionsFresnelAtCoreCladInterface.size(); i<trackId+1; i++)
reflectionsFresnelAtCoreCladInterface.push_back(0);
for(int i=refractionsAtCladCladInterface.size(); i<trackId+1; i++)
refractionsAtCladCladInterface.push_back(0);
for(int i=refractionsAtCoreCladInterface.size(); i<trackId+1; i++)
refractionsAtCoreCladInterface.push_back(0);
for(int i=rayleighScatterings.size(); i<trackId+1; i++)
rayleighScatterings.push_back(0);
}
void Analysis::IncreaseLengthInCore(Int_t trackId, Float_t lengthValue)
{
lengthInCore[trackId] += lengthValue;
}
void Analysis::IncreaseLengthInInnerCladding(Int_t trackId, Float_t lengthValue)
{
lengthInInnerCladding[trackId] += lengthValue;
}
void Analysis::IncreaseLengthInOuterCladding(Int_t trackId, Float_t lengthValue)
{
lengthInOuterCladding[trackId] += lengthValue;
}
void Analysis::IncreaseLengthVectors(Int_t trackId)
{
for(int i=lengthInCore.size(); i<trackId+1; i++)
lengthInCore.push_back(0);
for(int i=lengthInInnerCladding.size(); i<trackId+1; i++)
lengthInInnerCladding.push_back(0);
for(int i=lengthInOuterCladding.size(); i<trackId+1; i++)
lengthInOuterCladding.push_back(0);
}
SciFiSim/src/DetectorConstruction.cc
View file @ 07fea035
This diff is collapsed.
SciFiSim/src/DetectorConstruction.cc~ deleted 100755 → 0
View file @ 6f8f2267
This diff is collapsed.
SciFiSim/src/Parameters.cc
View file @ 07fea035
This diff is collapsed.
SciFiSim/src/PhysicsList.cc
View file @ 07fea035
......@@ -11,6 +11,7 @@
#include "G4EmStandardPhysics.hh"
#include "G4DecayPhysics.hh"
#include "G4OpticalPhysics.hh"
#include "G4OpticalParameters.hh"
#include "G4ProcessManager.hh"
#include "G4ParticleTypes.hh"
#include "G4hIonisation.hh"
......@@ -27,21 +28,37 @@ PhysicsList::PhysicsList() : G4VModularPhysicsList()
// * EM Physics
RegisterPhysics( new G4EmStandardPhysics() );
// * Optical Physics
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
RegisterPhysics( opticalPhysics );
// adjust some parameters for the optical physics
opticalPhysics->SetWLSTimeProfile("exponential");
auto* optParams = G4OpticalParameters::Instance();
// optParams->SetWLSProcess(true);
// optParams->SetRayleigh(true);
optParams->SetWLSTimeProfile("exponential");
// NOTE: SetScintillationYieldFactor() and SetScintillationExcitationRatio()
// no longer exist in v11. Instead you must define material properties:
// in your DetectorConstruction (or wherever you build the MPT):
//
// mpt->AddConstProperty("SCINTILLATIONYIELD", 1.0);
// mpt->AddConstProperty("SCINTILLATIONYIELDRATIO", Parameters::GetInstance()->YieldRatio());
//
// For particle-dependent yields use the ALPHASCINTILLATIONYIELD, PROTONSCINTILLATIONYIELD, etc. keys .
// If you need particle-type–specific yields:
optParams->SetScintByParticleType(false);
opticalPhysics->SetScintillationYieldFactor(1.0);
opticalPhysics->SetScintillationExcitationRatio(Parameters::GetInstance()->YieldRatio());
// Cerenkov parameters: the methods were renamed in v11
optParams->SetCerenkovMaxPhotonsPerStep(100); // was SetMaxNumPhotonsPerStep
optParams->SetCerenkovMaxBetaChange(10.0); // was SetMaxBetaChangePerStep
opticalPhysics->SetMaxNumPhotonsPerStep(100);
opticalPhysics->SetMaxBetaChangePerStep(10.0);
// Optional: track optical secondaries first
// optParams->SetScintTrackSecondariesFirst(true);
// opticalPhysics->SetTrackSecondariesFirst(true);
}
......
SciFiSim/src/PhysicsList.cc~ deleted 100755 → 0
View file @ 6f8f2267
// Written by Mirco DECKENHOFF based on
//
// $Id: PhysicsList.cc,v 1.17 2009/11/10 05:16:23 gum Exp $
//
//
// modified by Peter Stromberger
// modified by Bastian Rössler
#include "PhysicsList.hh"
#include "G4EmStandardPhysics.hh"
#include "G4DecayPhysics.hh"
#include "G4OpticalPhysics.hh"
#include "G4ProcessManager.hh"
#include "G4ParticleTypes.hh"
#include "G4hIonisation.hh"
#include "Parameters.hh"
#include "G4SystemOfUnits.hh"
PhysicsList::PhysicsList() : G4VModularPhysicsList()
{
// default cut value (1.0mm)
defaultCutValue = 1.0*mm;
// * EM Physics
RegisterPhysics( new G4EmStandardPhysics() );
/*
// * Optical Physics
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
RegisterPhysics( opticalPhysics );
// adjust some parameters for the optical physics
opticalPhysics->SetWLSTimeProfile("exponential");
opticalPhysics->SetScintillationYieldFactor(1.0);
opticalPhysics->SetScintillationExcitationRatio(Parameters::GetInstance()->YieldRatio());
opticalPhysics->SetMaxNumPhotonsPerStep(100);
opticalPhysics->SetMaxBetaChangePerStep(10.0);
*/
// opticalPhysics->SetTrackSecondariesFirst(true);
}
void PhysicsList::ConstructParticle()
{
// Constructs all paricles
G4VModularPhysicsList::ConstructParticle();
}
void PhysicsList::SetCuts()
{
// " G4VUserPhysicsList::SetCutsWithDefault" method sets
// the default cut value for all particle types
//
SetCutsWithDefault();
if (verboseLevel>0)
DumpCutValuesTable();
}
SciFiSim/src/PrimaryGeneratorAction.cc
View file @ 07fea035
// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler
// modified by Bastian Rössler, Lai Gui
#include "PrimaryGeneratorAction.hh"
#include "G4Event.hh"
......@@ -26,11 +26,17 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
G4double runID = (G4double) G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
G4double eventID = (G4double) anEvent->GetEventID();
G4double energy = (G4double) gun->GetParticleEnergy();
G4ThreeVector position = gun->GetParticlePosition();
G4double xPos = position.x();
G4double yPos = position.y();
G4double zPos = position.z();
G4double xMom = (G4double) gun->GetParticleMomentumDirection()[0];
G4double yMom = (G4double) gun->GetParticleMomentumDirection()[1];
G4double zMom = (G4double) gun->GetParticleMomentumDirection()[2];
Analysis::GetInstance()->FillInitialParticle(runID, eventID, energy, xMom, yMom, zMom);
Analysis::GetInstance()->FillInitialParticle(runID, eventID, energy, xPos, yPos, zPos, xMom, yMom, zMom);
// Necessery?
Analysis::GetInstance()->SetGpsPosition(gun->GetParticlePosition());
......
SciFiSim/src/SensitiveDetector.cc
View file @ 07fea035
......@@ -55,31 +55,31 @@ G4bool SensitiveDetector::ProcessHits(G4Step * step, G4TouchableHistory *)
G4double yPixel = G4UIcommand::ConvertToDouble((G4String)
(step->GetPreStepPoint()->GetPhysicalVolume()->GetName()).substr(21,2));
Analysis::GetInstance()->FillDetectedPhotons(myRunID,
myEventID,
myDetectorID,
xPixel,
yPixel,
step->GetTrack()->GetTotalEnergy()*1e6,
step->GetTrack()->GetLocalTime()-step->GetDeltaTime(),
step->GetTrack()->GetTrackLength()-step->GetStepLength(),
step->GetTrack()->GetGlobalTime()-step->GetDeltaTime(),
step->GetTrack()->GetPosition()[0]-step->GetDeltaPosition()[0],
step->GetTrack()->GetTotalEnergy()*1e6, // energy
step->GetTrack()->GetLocalTime()-step->GetDeltaTime(), // local time
step->GetTrack()->GetTrackLength()-step->GetStepLength(), // length
step->GetTrack()->GetGlobalTime()-step->GetDeltaTime(), // absolute time
step->GetTrack()->GetPosition()[0]-step->GetDeltaPosition()[0], // x y z
step->GetTrack()->GetPosition()[1]-step->GetDeltaPosition()[1],
step->GetTrack()->GetPosition()[2]-step->GetDeltaPosition()[2],
step->GetDeltaPosition()[0]/step->GetStepLength(),
step->GetDeltaPosition()[1]/step->GetStepLength(),
step->GetDeltaPosition()[2]/step->GetStepLength(),
step->GetTrack()->GetVertexPosition()[0],
step->GetTrack()->GetVertexPosition()[0], // vertex x y z
step->GetTrack()->GetVertexPosition()[1],
step->GetTrack()->GetVertexPosition()[2],
step->GetTrack()->GetVertexMomentumDirection()[0],
step->GetTrack()->GetVertexMomentumDirection()[0], // vertex px py pz
step->GetTrack()->GetVertexMomentumDirection()[1],
step->GetTrack()->GetVertexMomentumDirection()[2],
step->GetTrack()->GetTrackID(),
step->GetTrack()->GetTrackID(), // trackId
creatorProcessId,
step->GetTrack()->GetParentID());
step->GetTrack()->GetParentID()); // parentId
}
// Also kill reflected photons to avoid many reflections without detection
......
SciFiSim/src/SteppingAction.cc
View file @ 07fea035
......@@ -57,25 +57,29 @@ void SteppingAction::UserSteppingAction( const G4Step * step )
{
if(postPhysicalVolumeName == "EpoxyBox" && prePhysicalVolumeName.substr(0,9) == "Cladding2")
{
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() == 1)
G4double probLoss = Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface();
if (probLoss == 1.0)
{
// Always lose the photon
step->GetTrack()->SetTrackStatus(fStopAndKill);
}
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() < 1
&& Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() > 0)
else
{
G4double randomNumber = CLHEP::RandFlat::shoot();
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() < randomNumber)
if (randomNumber < probLoss)
{
// Kill photon
step->GetTrack()->SetTrackStatus(fStopAndKill);
}
else
{
// Survives, count reflection
if (step->GetStepLength() > 0)
Analysis::GetInstance()->IncreaseReflectionsAtFibreSurface(trackId);
}
else
step->GetTrack()->SetTrackStatus(fStopAndKill);
}
if (step->GetStepLength() > 0)
Analysis::GetInstance()->IncreaseReflectionsAtFibreSurface(trackId);
}
if(postPhysicalVolumeName == "Mirror")
Analysis::GetInstance()->IncreaseReflectionsAtMirror(trackId);
......
SciFiSim/src/SteppingAction.cc~ deleted 100755 → 0
View file @ 6f8f2267
// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler
#include "SteppingAction.hh"
#include "G4Step.hh"
#include "G4VTouchable.hh"
#include "G4SteppingManager.hh"
#include "G4UnitsTable.hh"
#include "G4ParticleTypes.hh"
#include "G4ProcessManager.hh"
#include "G4VProcess.hh"
#include "G4ProcessVector.hh"
#include "EventAction.hh"
#include <string>
#include "Parameters.hh"
#include "Analysis.hh"
#include "G4RunManager.hh"
#include "Randomize.hh"
#include "G4SystemOfUnits.hh"
SteppingAction::SteppingAction(EventAction* eventAction, PrimaryGeneratorAction* p)
{
creatorProcess = -1;
fEventAction = eventAction;
pGA = p;
}
SteppingAction::~SteppingAction(){}
void SteppingAction::UserSteppingAction( const G4Step * step )
{
primDefinition = (pGA->GetGun())->GetParticleDefinition();
primParticle = (primDefinition->GetParticleName()).c_str();
/* ++ SteppingAction for optical photons ++ */
if(step->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
{
G4int trackId = step->GetTrack()->GetTrackID();
G4double stepLength = step->GetStepLength();
G4String prePhysicalVolumeName = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
G4String postPhysicalVolumeName = "";
// PostStep does only exist if the particle is not exiting the EpoxyBox volume!
if(step->GetTrack()->GetTrackStatus() != fStopAndKill)
postPhysicalVolumeName = step->GetPostStepPoint()->GetPhysicalVolume()->GetName();
// Kill photons at fibre surface
// Count the reflections and refractions
if(step->GetPostStepPoint()->GetStepStatus()==fGeomBoundary)
{
if(postPhysicalVolumeName == "EpoxyBox" && prePhysicalVolumeName.substr(0,9) == "Cladding2")
{
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() == 1)
step->GetTrack()->SetTrackStatus(fStopAndKill);
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() < 1
&& Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() > 0)
{
G4double randomNumber = CLHEP::RandFlat::shoot();
if(Parameters::GetInstance()->ProbabilityOfPhotonLossAtSurface() < randomNumber)
{
if (step->GetStepLength() > 0)
Analysis::GetInstance()->IncreaseReflectionsAtFibreSurface(trackId);
}
else
step->GetTrack()->SetTrackStatus(fStopAndKill);
}
if (step->GetStepLength() > 0)
Analysis::GetInstance()->IncreaseReflectionsAtFibreSurface(trackId);
}
if(postPhysicalVolumeName == "Mirror")
Analysis::GetInstance()->IncreaseReflectionsAtMirror(trackId);
// Retrieve status of the photon boundary process
G4OpBoundaryProcessStatus theBoundaryStatus = Undefined;
G4ProcessManager* OpManager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
if(OpManager)
{
G4int MAXofPostStepLoops = OpManager->GetPostStepProcessVector()->entries();
G4ProcessVector* fPostStepDoItVector = OpManager->GetPostStepProcessVector(typeDoIt);
for(G4int i = 0; i<MAXofPostStepLoops; i++)
{
G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i];
opBoundaryProcess = dynamic_cast<G4OpBoundaryProcess*>(fCurrentProcess);
if(opBoundaryProcess)
{
theBoundaryStatus = opBoundaryProcess->GetStatus();
break;
}
}
}
if(theBoundaryStatus != Undefined)
{
if((prePhysicalVolumeName.substr(0,9) == "Cladding1" && postPhysicalVolumeName.substr(0,9) == "Cladding2")
|| (prePhysicalVolumeName.substr(0,9) == "Cladding2" && postPhysicalVolumeName.substr(0,9) == "Cladding1"))
{
if(theBoundaryStatus == TotalInternalReflection)
Analysis::GetInstance()->IncreaseTotalReflectionsAtCladCladInterface(trackId);
if(theBoundaryStatus == FresnelReflection)
Analysis::GetInstance()->IncreaseFresnelReflectionsAtCladCladInterface(trackId);
if(theBoundaryStatus == FresnelRefraction)
Analysis::GetInstance()->IncreaseRefractionsAtCladCladInterface(trackId);
}
if ((prePhysicalVolumeName.substr(0,4) == "Core" && postPhysicalVolumeName.substr(0,9) == "Cladding1")
|| (prePhysicalVolumeName.substr(0,9) == "Cladding1" && postPhysicalVolumeName.substr(0,9) == "Core"))
{
if(theBoundaryStatus == TotalInternalReflection)
Analysis::GetInstance()->IncreaseTotalReflectionsAtCoreCladInterface(trackId);
if(theBoundaryStatus == FresnelReflection)
Analysis::GetInstance()->IncreaseFresnelReflectionsAtCoreCladInterface(trackId);
if(theBoundaryStatus == FresnelRefraction)
Analysis::GetInstance()->IncreaseRefractionsAtCoreCladInterface(trackId);
}
}
}
else // Check for Rayleigh scattering
if (step->GetDeltaMomentum()[0] != 0 || step->GetDeltaMomentum()[1] != 0 || step->GetDeltaMomentum()[2] != 0)
Analysis::GetInstance()->IncreaseRayleighScatterings(trackId);
// Store length per volume
if (prePhysicalVolumeName.substr(0,4) == "Core")
Analysis::GetInstance()->IncreaseLengthInCore(trackId,stepLength);
if (prePhysicalVolumeName.substr(0,9) == "Cladding1")
Analysis::GetInstance()->IncreaseLengthInInnerCladding(trackId,stepLength);
if (prePhysicalVolumeName.substr(0,9) == "Cladding2")
Analysis::GetInstance()->IncreaseLengthInOuterCladding(trackId,stepLength);
// Store photons killed by OpWLS process, if not primary particle
if(step->GetTrack()->GetParentID() > 0 && step->GetTrack()->GetTrackStatus() == fStopAndKill && step->GetSecondary()->size() == 1)
{
G4int creatorProcessId = 0;
if(step->GetTrack()->GetCreatorProcess()->GetProcessName() == "Scintillation")
creatorProcessId = 1;
if(step->GetTrack()->GetCreatorProcess()->GetProcessName() == "Cerenkov")
creatorProcessId = 2;
if(step->GetTrack()->GetCreatorProcess()->GetProcessName() == "OpWLS")
creatorProcessId = 3;
}
// Radiation damage
dose = h_dose->GetBinContent(1,1 );
dose = 0;
std::cout << step->GetTrack()->GetPosition().x() << " " << step->GetTrack()->GetPosition().y() << " " << step->GetTrack()->GetPosition().z() << " " << dose << "\n ";
G4float attlength = 100.*1./(dose * 1.5*exp(9.412 -0.0256721 * 1240./(step->GetTrack()->GetTotalEnergy()/eV)) );
double xf = G4UniformRand();
//std::cout << xf << " " << attlength<< " " << 1240./(step->GetTrack()->GetTotalEnergy()/eV) << "\n";
if(xf>(exp(-stepLength*cm/attlength)) ){ //randomly kill photon based on attenuation length
step->GetTrack()->SetTrackStatus(fStopAndKill);}
}
/* ++ End of SteppingAction of optical photons. ++ */
/* ++ SteppingAction for electrons and muon ++ */
if((step->GetTrack()->GetDefinition() == primDefinition))
{
G4double runID = (G4double) G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
G4double eventID = (G4double) G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
Analysis::GetInstance()->FillPrimaryParticleTrack(runID,
eventID,
step->GetTrack()->GetPosition()[0],
step->GetTrack()->GetPosition()[1],
step->GetTrack()->GetPosition()[2]);
G4String prePhysicalVolumeName = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
if (prePhysicalVolumeName.substr(0,4) == "Core")
{
G4double edep = step->GetTotalEnergyDeposit();
G4double stepLength = step->GetStepLength();
fEventAction->AddCore(edep, stepLength);
}
if (prePhysicalVolumeName == "Trigger")
{
G4double edep = step->GetTotalEnergyDeposit();
Analysis::GetInstance()->FillTrigger(runID,
eventID,
edep,
step->GetTrack()->GetPosition()[0],
step->GetTrack()->GetPosition()[1],
step->GetTrack()->GetPosition()[2]);
}
}
/* ++ End of SteppingAction for electrons and muon ++ */
}