From 616679a0c069cc3367d0eefc15ad7c4eb1aa4297 Mon Sep 17 00:00:00 2001 From: Mirald Tuzi Date: Tue, 8 Feb 2022 15:54:44 +0100 Subject: [PATCH 1/2] update README with instructions to set up configurations for dichroic filter --- README.md | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/README.md b/README.md index 27711a3..9b1c895 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,43 @@ # G4OpSim Geant4-based optical simulation. + +# Configuring a dichroic filter in GEANT4 + +Within the detector construction section of the code, one defines the contstruction of the dichroic filter. + +## General remarks + +* For the dichroic filter, defining a G4OpticalSurface object is necessary: + * `G4OpticalSurface* filtersurf_opsurf = new G4OpticalSurface(filtersurf_name, dichroic, polished, dielectric_dichroic, 1.);` + * `filtersurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); //material may vary` + * `new G4LogicalSkinSurface("FILTER_SURFACE", filter_logic_vol, filtersurf_opsurf);` +* To characterise the transmittance of the dichroic filter, it is necessary to define a configuration file that defines transmittance $T(\lambda, \theta)$ as a function of the wavelength $\lambda$ of the photon and the incident angle $\theta$ (_**defined with respect to to the normal vector of the plane**_). A template will be shown in the next section. After the file is set up, its path needs to be called as an environment variable in the detector class method where one defines the dichroic filter + * setenv("G4DICHROICDATA", "data/dichroic_data", 1); + * The environment variable name has to be as written above + +## Configuration of Transmittance + +The following variables are used for simplification: + +* `G4VectorType`: a value telling GEANT4 what type of vector to pick from an enum. For the dichroic filter, we set this to 4 +* `m `: Number of wavelength samples +* `n `: Number of angle samples +* `wl_i`: $i$-th wavelength value ($1 + ... + ... + ... +. . . +. . . +. . . + ... +``` + +Remarks on restrictions: +* $m,n\geq2$, otherwise GEANT4 does not accept the configuration file +* There can be formatting issues; the file has been tested that it works if the encoding is UTF-8 and there is always exactly one white space between values in a row and no extra whitespaces otherwise and no extra empty lines at the end of the file From 212dcd8beb26229116347d2cbac80146a65b5b39 Mon Sep 17 00:00:00 2001 From: mhw0i Date: Tue, 8 Feb 2022 16:01:05 +0100 Subject: [PATCH 2/2] Configuration of Dichroic filter according to ASAHI model data and rudimentary data generation for post analysis --- src/Analysis.h | 6 + src/CMakeLists.txt | 18 +- src/DetectorConstruction.cpp | 1136 +++++++++++++++-------------- src/DetectorConstruction.h | 81 +- src/EventAction.cpp | 58 +- src/EventAction.h | 59 +- src/GenericPhotosensor.cpp | 238 +++--- src/GenericPhotosensor.h | 76 +- src/Materials.cpp | 188 ++--- src/Materials.h | 54 +- src/OpticalHit.cpp | 132 ++-- src/OpticalHit.h | 134 ++-- src/OpticalMaterialProperties.cpp | 726 +++++++++--------- src/OpticalMaterialProperties.h | 82 +-- src/OpticalSD.cpp | 80 +- src/OpticalSD.h | 60 +- src/PhysicsList.cpp | 68 +- src/PhysicsList.h | 42 +- src/PrimaryGeneration.cpp | 183 +++-- src/PrimaryGeneration.h | 54 +- src/RunAction.cpp | 126 +++- src/RunAction.h | 54 +- src/SteppingAction.cpp | 324 +++++--- src/SteppingAction.h | 61 +- src/TrackingAction.cpp | 30 +- src/TrackingAction.h | 54 +- src/params.h | 3 + 27 files changed, 2247 insertions(+), 1880 deletions(-) create mode 100644 src/Analysis.h create mode 100644 src/params.h diff --git a/src/Analysis.h b/src/Analysis.h new file mode 100644 index 0000000..cccc062 --- /dev/null +++ b/src/Analysis.h @@ -0,0 +1,6 @@ +#include "g4root.hh" +//#include "g4csv.hh" +//#include "g4xml.hh" + +//#include "G4GenericAnalysisManager.hh" +//using G4AnalysisManager = G4GenericAnalysisManager; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b59de97..fe2ce75 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,9 +1,9 @@ -## ----------------------------------------------------------------------------- -## G4OpSim | src/CMakeLists.txt -## -## * Author: -## * Creation date: 10 February 2020 -## ----------------------------------------------------------------------------- - -file(GLOB SRCS ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) -add_library(${CMAKE_PROJECT_NAME}_SRC OBJECT ${SRCS}) +## ----------------------------------------------------------------------------- +## G4OpSim | src/CMakeLists.txt +## +## * Author: +## * Creation date: 10 February 2020 +## ----------------------------------------------------------------------------- + +file(GLOB SRCS ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) +add_library(${CMAKE_PROJECT_NAME}_SRC OBJECT ${SRCS}) diff --git a/src/DetectorConstruction.cpp b/src/DetectorConstruction.cpp index 709a395..d42c29c 100644 --- a/src/DetectorConstruction.cpp +++ b/src/DetectorConstruction.cpp @@ -1,533 +1,603 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | DetectorConstruction.cpp -// -// * Author: -// * Creation date: 10 February 2020 -// ----------------------------------------------------------------------------- - -#include "DetectorConstruction.h" - -#include "GenericPhotosensor.h" -#include "Materials.h" -#include "OpticalMaterialProperties.h" -#include "OpticalSD.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -DetectorConstruction::DetectorConstruction(): - G4VUserDetectorConstruction(), - world_size_(10.*m), - plate_width_ (112.0*mm), // X - plate_thickn_( 4.0*mm), // Y - plate_length_(491.5*mm), // Z - foil_thickn_(0.165*mm), - num_phsensors(48) -{ -} - - -DetectorConstruction::~DetectorConstruction() -{ -} - - -G4VPhysicalVolume* DetectorConstruction::Construct() -{ - // WORLD /////////////////////////////////////////////////////////// - // Sphere of liquid argon that contains all other volumes. - - const G4String world_name = "WORLD"; - - G4Material* LAr = G4NistManager::Instance()->FindOrBuildMaterial("G4_lAr"); - LAr->SetMaterialPropertiesTable(OpticalMaterialProperties::LAr()); - - G4Sphere* world_solid_vol = - new G4Sphere(world_name, 0., world_size_/2., 0., 360.*deg, 0., 180.*deg); - - G4LogicalVolume* world_logic_vol = - new G4LogicalVolume(world_solid_vol, LAr, world_name); - world_logic_vol->SetVisAttributes(G4VisAttributes::Invisible); - - G4VPhysicalVolume* world_phys_vol = - new G4PVPlacement(nullptr, G4ThreeVector(), - world_logic_vol, world_name, nullptr, false, 0, true); - - //////////////////////////////////////////////////////////////////// - - ConstructWLSPlate(world_phys_vol); - ConstructPhotosensors(world_phys_vol); - ConstructReflectiveFoils(world_phys_vol); - - return world_phys_vol; -} - - -void DetectorConstruction::ConstructWLSPlate(G4VPhysicalVolume* world_phys_vol) const -{ - // WLS PLATE /////////////////////////////////////////////////////// - - Assert(world_phys_vol, "DetectorConstruction::ConstructWLSPlate()"); - - const G4String plate_name = "WLS_PLATE"; - - G4Box* plate_solid_vol = - new G4Box(plate_name, plate_width_/2., plate_thickn_/2., plate_length_/2.); - - G4Material* pvt = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); - pvt->SetMaterialPropertiesTable(OpticalMaterialProperties::PVT()); - - G4LogicalVolume* plate_logic_vol = - new G4LogicalVolume(plate_solid_vol, pvt, plate_name); - - new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,0.), - plate_logic_vol, plate_name, world_phys_vol->GetLogicalVolume(), - false, 0, true); -} - - -void DetectorConstruction::ConstructPhotosensors(G4VPhysicalVolume* world_phys_vol) const -{ - Assert(world_phys_vol, "DetectorConstruction::ConstructPhotosensors()"); - - GenericPhotosensor photosensor_geom; - photosensor_geom.Construct(); - G4LogicalVolume* photosensor_logic_vol = photosensor_geom.GetLogicalVolume(); - - if (!photosensor_logic_vol) { - G4Exception("DetectorConstruction::ConstructPhotosensors()", "nullptr", - FatalException, " Description : Null pointer to logical volume."); - } - - G4int phsensor_id = 0; - - G4RotationMatrix* rot = new G4RotationMatrix(); - rot->rotateY(-90*deg); - - for (G4int i=0; iGetName(), - world_phys_vol->GetLogicalVolume(), - false, phsensor_id, true); - } - - G4RotationMatrix* rot2 = new G4RotationMatrix(); - rot2->rotateY(90*deg); - - for (G4int i=0; iGetName(), - world_phys_vol->GetLogicalVolume(), - false, phsensor_id, true); - } - - // ////////////////////////////////////////////////////////// - // // describing photosensors position - // const G4int nsens = 48; //n SiPMs per arapuca cell - // G4int nsides = 2; // n sides of the sell covered with SiPMs - // const G4double sensSpacing = plate_length/(nsens/2)*mm; - // G4double firstPosition = sensSpacing/2; - // G4ThreeVector SiPMpos; - // G4double Zpos; - // G4double Xpos; - // G4int nsensStart; - // G4String SiPM_name = "SiPM"; - // G4RotationMatrix* rotationMatrix = new G4RotationMatrix(); - // - // G4VPhysicalVolume* photosens_phys_vol[nsens]; - // - // for (G4int iside = 0; iside < nsides; iside ++) - // { - // if(iside == 0) - // { - // Xpos = plate_width/2 + 6*mm; - // rotationMatrix->rotateY(90.*deg); //they have to be facing de plate - // nsensStart = 0; - // //rotationMatrix->print(std::cout); - // - // } - // else if(iside == 1) - // { - // Xpos = -(plate_width/2+6*mm); - // rotationMatrix->rotateY(180.*deg); //rotations apply sequently, so this rotation is 90º+180º - // nsensStart = nsens/2; - // //rotationMatrix->print(std::cout); - // } - // for (G4int isens = 0; isens < nsens/2; isens++) - // { - // Zpos = -plate_length/2 + firstPosition + isens * sensSpacing; - // // Not sure if the sensors shuld be plased in the world logic word - // // When putting them in plate_logic_vol it produces overlap errors - // photosens_phys_vol[nsensStart+isens]= - // new G4PVPlacement(rotationMatrix,G4ThreeVector(Xpos,0.,Zpos), - // this->GenericPhotosensor(),SiPM_name, world_logic_vol, - // false, isens, true); - // } - // } - // - // return world_phys_vol; - -} - -void DetectorConstruction::ConstructReflectiveFoils(G4VPhysicalVolume* world_phys_vol) const -{ - // REFLECTIVE FOILS /////////////////////////////////////////////////////// - - Assert(world_phys_vol, "DetectorConstruction::ConstructWLSPlate()"); - - //there are three reflective foils, one in the bottom, two in the sides - //first create the logical volume and then create the surface with LAr - const G4String bottom_foil_name = "REF_BOTTOM_FOIL"; - const G4String side_foil_name = "REF_SIDE_FOIL"; - - //for the moment, same area as WLS plate. Should be checked. - //thickness from manufacturer, should be checked. - //https://www.isoltronic.ch/assets/of-m-vikuiti-esr-app-guide.pdf - - G4Box* bottom_foil_solid_vol = - new G4Box(bottom_foil_name, plate_width_/2. ,foil_thickn_/2., plate_length_/2.); - - G4Box* side_foil_solid_vol = - new G4Box(side_foil_name, plate_width_/2., plate_thickn_/2., foil_thickn_/2.); - - G4Material* plexiglass = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLEXIGLASS"); - - G4LogicalVolume* bottom_foil_logic_vol = - new G4LogicalVolume(bottom_foil_solid_vol, plexiglass, bottom_foil_name); - - G4LogicalVolume* side_foil_logic_vol = - new G4LogicalVolume(side_foil_solid_vol, plexiglass, side_foil_name); - - //the position is still unknown. For the moment put 1mm, should be checked - G4ThreeVector bottom_foil_pos(0, -plate_thickn_/2 - foil_thickn_/2 - 1*mm, 0); - new G4PVPlacement(nullptr, bottom_foil_pos, - bottom_foil_logic_vol, bottom_foil_name, - world_phys_vol->GetLogicalVolume(), - false, 0, true); - - G4ThreeVector side_foil_pos(0, 0, plate_length_/2 + foil_thickn_/2 + 1*mm); - new G4PVPlacement(nullptr, side_foil_pos, - side_foil_logic_vol, side_foil_name, - world_phys_vol->GetLogicalVolume(), - false, 0, true); - - new G4PVPlacement(nullptr, -side_foil_pos, - side_foil_logic_vol, side_foil_name, - world_phys_vol->GetLogicalVolume(), - false, 0, true); - - //now create the surface - const G4String refsurf_name = "REF_SURFACE"; - G4OpticalSurface* refsurf_opsurf = - new G4OpticalSurface(refsurf_name, unified, polishedfrontpainted, dielectric_dielectric, 1); - - refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::VIKUITI()); - new G4LogicalSkinSurface("REF_FOIL_SURFACE",bottom_foil_logic_vol,refsurf_opsurf); - new G4LogicalSkinSurface("REF_FOIL_SURFACE",side_foil_logic_vol,refsurf_opsurf); -} - -void DetectorConstruction::Assert(G4VPhysicalVolume* ptr, const G4String& origin) const -{ - if (!ptr) { - G4ExceptionDescription ed; - ed << " Description : World physical volume is not defined."; - G4Exception(origin, "nullptr", FatalException, ed); - } -} - - - - - // // REFLECTIVE FOIL ///////////////////////////////////////////// - // - // //there are three reflective foils, one in the bottom, two in the sides - // //first create the logical volume and then create the surface with LAr - // const G4String bottom_foil_name = "REF_BOTTOM_FOIL"; - // const G4String side_foil_name = "REF_SIDE_FOIL"; - // - // //for the moment, same area as WLS plate. Should be checked. - // //thickness from manufacturer, should be checked. - // //https://www.isoltronic.ch/assets/of-m-vikuiti-esr-app-guide.pdf - // const G4double bottom_foil_width = 112.0*mm; // X - // const G4double bottom_foil_thickn = 0.165*mm; // Y - // const G4double bottom_foil_length = 491.5*mm; // Z - // - // const G4double side_foil_width = 112.0*mm; // X - // const G4double side_foil_thickn = 4.0*mm; - // const G4double side_foil_length = 0.165*mm; // Y - // - // G4Box* bottom_foil_solid_vol = - // new G4Box(bottom_foil_name, bottom_foil_width/2., - // bottom_foil_thickn/2., bottom_foil_length/2.); - // - // G4Box* side_foil_solid_vol = - // new G4Box(side_foil_name, side_foil_width/2., - // side_foil_thickn/2., side_foil_length/2.); - // - // G4Material* plexiglass = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLEXIGLASS"); - // - // G4LogicalVolume* bottom_foil_logic_vol = new - // G4LogicalVolume(bottom_foil_solid_vol, plexiglass, bottom_foil_name); - // - // G4LogicalVolume* side_foil_logic_vol = new - // G4LogicalVolume(side_foil_solid_vol, plexiglass, side_foil_name); - // - // //the position is still unknown. For the moment put 3mm, should be checked - // const G4double bottom_foil_posy = -bottom_foil_thickn/2. - 3*mm; - // const G4double side_foil_posz = bottom_foil_length/2 + 3*mm; - // - // G4VPhysicalVolume* bottom_foil_phys_vol = - // new G4PVPlacement(nullptr, G4ThreeVector(0.,bottom_foil_posy,0.), - // bottom_foil_logic_vol, bottom_foil_name, world_logic_vol, - // false, 0, true); - // - // G4VPhysicalVolume* side1_foil_phys_vol = - // new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,side_foil_posz), - // side_foil_logic_vol, side_foil_name, world_logic_vol, - // false, 0, true); - // - // G4VPhysicalVolume* side2_foil_phys_vol = - // new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,-side_foil_posz), - // side_foil_logic_vol, side_foil_name, world_logic_vol, - // false, 0, true); - // - // //now create the surface - // const G4String refsurf_name = "REF_SURFACE"; - // G4OpticalSurface* refsurf_opsurf = new - // G4OpticalSurface(refsurf_name, unified, polishedfrontpainted, dielectric_dielectric, 1); - // - // refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::VIKUITI()); - // new G4LogicalSkinSurface("REF_FOIL_SURFACE",bottom_foil_logic_vol,refsurf_opsurf); - // new G4LogicalSkinSurface("REF_FOIL_SURFACE",side_foil_logic_vol,refsurf_opsurf); - // - // - // // DICHROIC FILTER /////////////////////////////////////// - // - // const G4String filter_name = "DICHROIC_FILTER"; - // - // //dimensions should be checked - // const G4double filter_width = 112.0*mm; // X - // const G4double filter_thickn = 1.0*mm; // Y - // const G4double filter_length = 491.5*mm; // Z - // - // G4Box* filter_solid_vol = - // new G4Box(filter_name, filter_width/2., filter_thickn/2., filter_length/2.); - // - // G4Material* filter_material = G4NistManager::Instance()->FindOrBuildMaterial("G4_SILICON_DIOXIDE"); - // filter_material->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); - // - // G4LogicalVolume* filter_logic_vol = - // new G4LogicalVolume(filter_solid_vol,filter_material,filter_name); - // - // G4double filter_wls_gap = 1.5*mm; - // G4double filter_ypos = plate_thickn/2. + filter_wls_gap + filter_thickn/2.; - // - // G4VPhysicalVolume* filter_phys_vol = - // new G4PVPlacement(nullptr, G4ThreeVector(0.,filter_ypos,0.), - // filter_logic_vol, filter_name, world_logic_vol, - // false, 0, true); - // - // //now create the surface. It is a dichroic filter, so we need more things - // setenv("G4DICHROICDATA","data/dichroic_data",1); - // const G4String filtersurf_name = "FILTER_SURFACE"; - // G4OpticalSurface* filtersurf_opsurf = new - // G4OpticalSurface(filtersurf_name, dichroic, polished, dielectric_dichroic, 1); - // - // refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); - // new G4LogicalSkinSurface("FILTER_SURFACE",filter_logic_vol,filtersurf_opsurf); - // - // - // // PTP LAYER /////////////////////////////////////// - // - // const G4String ptp_name = "PTP_LAYER"; - // - // //real dimensions not known. To be reviewed - // const G4double ptp_width = 112.0*mm; // X - // const G4double ptp_thickn = 0.002*mm; // Y from https://arxiv.org/pdf/1912.09191.pdf - // const G4double ptp_length = 491.5*mm; // Z - // - // G4Box* ptp_solid_vol = - // new G4Box(ptp_name, ptp_width/2., ptp_thickn/2., ptp_length/2.); - // - // G4Material* ptp = G4NistManager::Instance()->FindOrBuildMaterial("G4_TERPHENYL"); - // ptp->SetMaterialPropertiesTable(OpticalMaterialProperties::PTP()); - // - // G4LogicalVolume* ptp_logic_vol = new G4LogicalVolume(ptp_solid_vol, ptp, ptp_name); - // - // //position still not known, to be reviewed - // const G4double ptp_posy = plate_thickn/2 + filter_wls_gap + filter_thickn + ptp_thickn/2; - // //1.5*mm from https://indico.fnal.gov/event/45283/contributions/195721/attachments/133823/165234/X-ARAPUCA_Cuts_Study.pdf - // //and 1*mm is the assumed filter width - // G4VPhysicalVolume* ptp_phys_vol = - // new G4PVPlacement(nullptr, G4ThreeVector(0.,ptp_posy,0.), - // ptp_logic_vol, ptp_name, world_logic_vol, - // false, 0, true); - // - // - // ////////////////////////////////////////////////////////// - // // describing photosensors position - // const G4int nsens = 48; //n SiPMs per arapuca cell - // G4int nsides = 2; // n sides of the sell covered with SiPMs - // const G4double sensSpacing = plate_length/(nsens/2)*mm; - // G4double firstPosition = sensSpacing/2; - // G4ThreeVector SiPMpos; - // G4double Zpos; - // G4double Xpos; - // G4int nsensStart; - // G4String SiPM_name = "SiPM"; - // G4RotationMatrix* rotationMatrix = new G4RotationMatrix(); - // - // G4VPhysicalVolume* photosens_phys_vol[nsens]; - // - // for (G4int iside = 0; iside < nsides; iside ++) - // { - // if(iside == 0) - // { - // Xpos = plate_width/2 + 6*mm; - // rotationMatrix->rotateY(90.*deg); //they have to be facing de plate - // nsensStart = 0; - // //rotationMatrix->print(std::cout); - // - // } - // else if(iside == 1) - // { - // Xpos = -(plate_width/2+6*mm); - // rotationMatrix->rotateY(180.*deg); //rotations apply sequently, so this rotation is 90º+180º - // nsensStart = nsens/2; - // //rotationMatrix->print(std::cout); - // } - // for (G4int isens = 0; isens < nsens/2; isens++) - // { - // Zpos = -plate_length/2 + firstPosition + isens * sensSpacing; - // // Not sure if the sensors shuld be plased in the world logic word - // // When putting them in plate_logic_vol it produces overlap errors - // photosens_phys_vol[nsensStart+isens]= - // new G4PVPlacement(rotationMatrix,G4ThreeVector(Xpos,0.,Zpos), - // this->GenericPhotosensor(),SiPM_name, world_logic_vol, - // false, isens, true); - // } - // } - // - // return world_phys_vol; -// } - - - - - - - - -// G4LogicalVolume* DetectorConstruction::GenericPhotosensor() -// { -// // PHOTOSENSOR ENCASING ////////////////////////////////// -// -// const G4double width_ = 6.*mm; -// const G4double height_ = 6.*mm; -// const G4double thickness_ = 2.*mm; -// -// G4String name = "PHOTOSENSOR"; -// -// G4Box* encasing_solid_vol = -// new G4Box(name, width_/2., height_/2., thickness_/2.); -// -// G4LogicalVolume* encasing_logic_vol = -// new G4LogicalVolume(encasing_solid_vol, Materials::FR4(), name); -// -// // OPTICAL WINDOW //////////////////////////////////////// -// -// name = "PHOTOSENSOR_WINDOW"; -// -// G4double window_thickness = thickness_/4.; -// -// G4Box* window_solid_vol = -// new G4Box(name, width_/2., height_/2., window_thickness/2.); -// -// -// G4LogicalVolume* window_logic_vol = -// new G4LogicalVolume(window_solid_vol, -// Materials::OpticalSilicone(), -// name); -// -// //THE SIPM IS CREATED FACING +Z DIRECTION -// G4double zpos = thickness_/2. - window_thickness/2.; -// -// new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), -// window_logic_vol, name, encasing_logic_vol, -// false, 0, false); -// -// // PHOTOSENSITIVE AREA ///////////////////////////////////////////// -// -// name = "PHOTOSENSOR_SENSAREA"; -// -// G4double sensarea_thickness = 0.1*mm; -// -// G4Box* sensarea_solid_vol = -// new G4Box(name, width_/2., height_/2., sensarea_thickness/2.); -// -// G4LogicalVolume* sensarea_logic_vol = -// new G4LogicalVolume(sensarea_solid_vol, -// G4NistManager::Instance()->FindOrBuildMaterial("G4_Si"), -// name); -// -// zpos = thickness_/2. - window_thickness - sensarea_thickness/2.; -// -// new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), -// sensarea_logic_vol, name, encasing_logic_vol, -// false, 0, false); -// -// // OPTICAL PROPERTIES ////////////////////////////////////////////// -// -// name = "PHOTOSENSOR_OPSURF"; -// -// G4double energy[] = {0.2*eV, 11.5*eV}; -// G4double reflectivity[] = {0.0 , 0.0}; -// G4double efficiency[] = {1.0 , 1.0}; -// -// G4MaterialPropertiesTable* photosensor_mpt = new G4MaterialPropertiesTable(); -// photosensor_mpt->AddProperty("REFLECTIVITY", energy, reflectivity, 2); -// photosensor_mpt->AddProperty("EFFICIENCY", energy, efficiency, 2); -// -// G4OpticalSurface* photosensor_opsurf = -// new G4OpticalSurface(name, unified, polished, dielectric_metal); -// photosensor_opsurf->SetMaterialPropertiesTable(photosensor_mpt); -// new G4LogicalSkinSurface(name, sensarea_logic_vol, photosensor_opsurf); -// -// // SENSITIVE DETECTOR ////////////////////////////////////////////// -// -// OpticalSD* sensdet = new OpticalSD("/GENERIC_PHOTOSENSOR/SiPM"); -// G4SDManager::GetSDMpointer()->AddNewDetector(sensdet); -// sensarea_logic_vol->SetSensitiveDetector(sensdet); -// -// //////////////////////////////////////////////////////////////////// -// -// return encasing_logic_vol; -// } +// ----------------------------------------------------------------------------- +// G4OpSim | DetectorConstruction.cpp +// +// * Author: +// * Creation date: 10 February 2020 +// ----------------------------------------------------------------------------- + +#include "DetectorConstruction.h" + +#include "GenericPhotosensor.h" +#include "Materials.h" +#include "OpticalMaterialProperties.h" +#include "OpticalSD.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +DetectorConstruction::DetectorConstruction(): + G4VUserDetectorConstruction(), + //world_size_(10.*m), + //plate_width_ (112.0*mm), // X + //plate_thickn_( 4.0*mm), // Y + //plate_length_(491.5*mm), // Z + //foil_thickn_(0.165*mm), + //num_phsensors(48) + world_size_(10. * m), + plate_width_(1000.0 * mm), // X 112.0 + plate_thickn_(4.0 * mm), // Y + plate_length_(1000.0 * mm), // Z 491.5 + foil_thickn_(0.165 * mm), // Z&X + foil_gap_(1. * mm), + filter_thickn_(1. * mm), // Y + filter_gap_(1.5 * mm), + ptp_thickn_(0.002 * mm), // Y + num_phsensors(48) +{ +} + + +DetectorConstruction::~DetectorConstruction() +{ +} + + +G4VPhysicalVolume* DetectorConstruction::Construct() +{ + // WORLD /////////////////////////////////////////////////////////// + // Sphere of liquid argon that contains all other volumes. + + const G4String world_name = "WORLD"; + + G4Material* LAr = G4NistManager::Instance()->FindOrBuildMaterial("G4_lAr"); + LAr->SetMaterialPropertiesTable(OpticalMaterialProperties::LAr()); + + G4Sphere* world_solid_vol = + new G4Sphere(world_name, 0., world_size_/2., 0., 360.*deg, 0., 180.*deg); + + G4LogicalVolume* world_logic_vol = + new G4LogicalVolume(world_solid_vol, LAr, world_name); + world_logic_vol->SetVisAttributes(G4VisAttributes::Invisible); + + G4VPhysicalVolume* world_phys_vol = + new G4PVPlacement(nullptr, G4ThreeVector(), + world_logic_vol, world_name, nullptr, false, 0, true); + + //////////////////////////////////////////////////////////////////// + + //ConstructWLSPlate(world_phys_vol); + //ConstructPhotosensors(world_phys_vol); + //ConstructReflectiveFoils(world_phys_vol); + ConstructDichroicFilter(world_phys_vol); + + return world_phys_vol; +} + + +void DetectorConstruction::ConstructWLSPlate(G4VPhysicalVolume* world_phys_vol) const +{ + // WLS PLATE /////////////////////////////////////////////////////// + + Assert(world_phys_vol, "DetectorConstruction::ConstructWLSPlate()"); + + const G4String plate_name = "WLS_PLATE"; + + G4Box* plate_solid_vol = + new G4Box(plate_name, plate_width_/2., plate_thickn_/2., plate_length_/2.); + + G4Material* pvt = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); + pvt->SetMaterialPropertiesTable(OpticalMaterialProperties::PVT()); + + G4LogicalVolume* plate_logic_vol = + new G4LogicalVolume(plate_solid_vol, pvt, plate_name); + + new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,0.), + plate_logic_vol, plate_name, world_phys_vol->GetLogicalVolume(), + false, 0, true); +} + + +void DetectorConstruction::ConstructPhotosensors(G4VPhysicalVolume* world_phys_vol) const +{ + Assert(world_phys_vol, "DetectorConstruction::ConstructPhotosensors()"); + + GenericPhotosensor photosensor_geom; + photosensor_geom.Construct(); + G4LogicalVolume* photosensor_logic_vol = photosensor_geom.GetLogicalVolume(); + + if (!photosensor_logic_vol) { + G4Exception("DetectorConstruction::ConstructPhotosensors()", "nullptr", + FatalException, " Description : Null pointer to logical volume."); + } + + G4int phsensor_id = 0; + + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateY(-90*deg); + + for (G4int i=0; iGetName(), + world_phys_vol->GetLogicalVolume(), + false, phsensor_id, true); + } + + G4RotationMatrix* rot2 = new G4RotationMatrix(); + rot2->rotateY(90*deg); + + for (G4int i=0; iGetName(), + world_phys_vol->GetLogicalVolume(), + false, phsensor_id, true); + } + + // ////////////////////////////////////////////////////////// + // // describing photosensors position + // const G4int nsens = 48; //n SiPMs per arapuca cell + // G4int nsides = 2; // n sides of the sell covered with SiPMs + // const G4double sensSpacing = plate_length/(nsens/2)*mm; + // G4double firstPosition = sensSpacing/2; + // G4ThreeVector SiPMpos; + // G4double Zpos; + // G4double Xpos; + // G4int nsensStart; + // G4String SiPM_name = "SiPM"; + // G4RotationMatrix* rotationMatrix = new G4RotationMatrix(); + // + // G4VPhysicalVolume* photosens_phys_vol[nsens]; + // + // for (G4int iside = 0; iside < nsides; iside ++) + // { + // if(iside == 0) + // { + // Xpos = plate_width/2 + 6*mm; + // rotationMatrix->rotateY(90.*deg); //they have to be facing de plate + // nsensStart = 0; + // //rotationMatrix->print(std::cout); + // + // } + // else if(iside == 1) + // { + // Xpos = -(plate_width/2+6*mm); + // rotationMatrix->rotateY(180.*deg); //rotations apply sequently, so this rotation is 90º+180º + // nsensStart = nsens/2; + // //rotationMatrix->print(std::cout); + // } + // for (G4int isens = 0; isens < nsens/2; isens++) + // { + // Zpos = -plate_length/2 + firstPosition + isens * sensSpacing; + // // Not sure if the sensors shuld be plased in the world logic word + // // When putting them in plate_logic_vol it produces overlap errors + // photosens_phys_vol[nsensStart+isens]= + // new G4PVPlacement(rotationMatrix,G4ThreeVector(Xpos,0.,Zpos), + // this->GenericPhotosensor(),SiPM_name, world_logic_vol, + // false, isens, true); + // } + // } + // + // return world_phys_vol; + +} + +void DetectorConstruction::ConstructReflectiveFoils(G4VPhysicalVolume* world_phys_vol) const +{ + // REFLECTIVE FOILS /////////////////////////////////////////////////////// + + Assert(world_phys_vol, "DetectorConstruction::ConstructWLSPlate()"); + + //there are three reflective foils, one in the bottom, two in the sides + //first create the logical volume and then create the surface with LAr + const G4String bottom_foil_name = "REF_BOTTOM_FOIL"; + const G4String side_foil_name = "REF_SIDE_FOIL"; + + //for the moment, same area as WLS plate. Should be checked. + //thickness from manufacturer, should be checked. + //https://www.isoltronic.ch/assets/of-m-vikuiti-esr-app-guide.pdf + + G4Box* bottom_foil_solid_vol = + new G4Box(bottom_foil_name, plate_width_/2. ,foil_thickn_/2., plate_length_/2.); + + G4Box* side_foil_solid_vol = + new G4Box(side_foil_name, plate_width_/2., plate_thickn_/2., foil_thickn_/2.); + + G4Material* plexiglass = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLEXIGLASS"); + + G4LogicalVolume* bottom_foil_logic_vol = + new G4LogicalVolume(bottom_foil_solid_vol, plexiglass, bottom_foil_name); + + G4LogicalVolume* side_foil_logic_vol = + new G4LogicalVolume(side_foil_solid_vol, plexiglass, side_foil_name); + + //the position is still unknown. For the moment put 1mm, should be checked + G4ThreeVector bottom_foil_pos(0, -plate_thickn_ / 2 - foil_thickn_ / 2 - foil_gap_, 0); + new G4PVPlacement(nullptr, bottom_foil_pos, + bottom_foil_logic_vol, bottom_foil_name, + world_phys_vol->GetLogicalVolume(), + false, 0, true); + + G4ThreeVector side_foil_pos(0, 0, plate_length_ / 2 + foil_thickn_ / 2 + foil_gap_); + new G4PVPlacement(nullptr, side_foil_pos, + side_foil_logic_vol, side_foil_name, + world_phys_vol->GetLogicalVolume(), + false, 0, true); + + new G4PVPlacement(nullptr, -side_foil_pos, + side_foil_logic_vol, side_foil_name, + world_phys_vol->GetLogicalVolume(), + false, 0, true); + + //now create the surface + const G4String refsurf_name = "REF_SURFACE"; + G4OpticalSurface* refsurf_opsurf = + new G4OpticalSurface(refsurf_name, unified, polishedfrontpainted, dielectric_dielectric, 1); + + refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::VIKUITI()); + new G4LogicalSkinSurface("REF_FOIL_SURFACE",bottom_foil_logic_vol,refsurf_opsurf); + new G4LogicalSkinSurface("REF_FOIL_SURFACE",side_foil_logic_vol,refsurf_opsurf); +} + +void DetectorConstruction::ConstructDichroicFilter(G4VPhysicalVolume* world_phys_vol) const //vol2 +{ + // DICHROIC FILTER /////////////////////////////////////// + const G4String filter_name = "DICHROIC_FILTER"; + + //dimensions should be checked + //1.5*mm from https://indico.fnal.gov/event/45283/contributions/195721/attachments/133823/165234/X-ARAPUCA_Cuts_Study.pdf + //and 1*mm is the assumed filter width + G4Box* filter_solid_vol = + new G4Box(filter_name, plate_width_ / 2., filter_thickn_ / 2., plate_length_ / 2.); //width: 112 mm, thickness: 1mm, length: 491.5 mm, + + //G4Material* filter_material = G4NistManager::Instance()->FindOrBuildMaterial("G4_SILICON_DIOXIDE"); + G4Material* filter_material = G4NistManager::Instance()->FindOrBuildMaterial("G4_lAr"); + filter_material->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); + + G4LogicalVolume* filter_logic_vol = + new G4LogicalVolume(filter_solid_vol, filter_material, filter_name); + + G4ThreeVector filter_pos(0, plate_thickn_ / 2 + filter_gap_ + filter_thickn_ / 2, 0); + G4VPhysicalVolume* filter_phys_vol = new G4PVPlacement(nullptr, filter_pos, // vol1 + filter_logic_vol, filter_name, + world_phys_vol->GetLogicalVolume(), + false, 0, true); + + //now create the surface. It is a dichroic filter, so we need more things... DONE + setenv("G4DICHROICDATA", "data/dichroic_data", 1); + const G4String filtersurf_name = "FILTER_SURFACE"; + G4OpticalSurface* filtersurf_opsurf = + new G4OpticalSurface(filtersurf_name, dichroic, polished, dielectric_dichroic, 1.); + + //filtersurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); + filtersurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::LAr()); + //new G4LogicalSkinSurface("FILTER_SURFACE", filter_logic_vol, filtersurf_opsurf); // replace with border, takes two physical, not logical volumes + new G4LogicalBorderSurface("FILTER_SURFACE", world_phys_vol, filter_phys_vol, filtersurf_opsurf); //physical volume of media photon passes through, in order of how photon passes through them +} + +void DetectorConstruction::ConstructPTPLayer(G4VPhysicalVolume* world_phys_vol) const +{ + // PTP LAYER /////////////////////////////////////// + const G4String ptp_name = "PTP_LAYER"; + + //real dimensions not known. To be reviewed + G4Box* ptp_solid_vol = + new G4Box(ptp_name, plate_width_ / 2., ptp_thickn_ / 2., plate_length_ / 2.); + + G4Material* ptp = G4NistManager::Instance()->FindOrBuildMaterial("G4_TERPHENYL"); + ptp->SetMaterialPropertiesTable(OpticalMaterialProperties::PTP()); + + G4LogicalVolume* ptp_logic_vol = new G4LogicalVolume(ptp_solid_vol, ptp, ptp_name); + + //position still not known, to be reviewed + G4ThreeVector ptp_pos(0, plate_thickn_ / 2 + filter_gap_ + filter_thickn_ / 2 + ptp_thickn_ / 2, 0); + new G4PVPlacement(nullptr, ptp_pos, + ptp_logic_vol, ptp_name, + world_phys_vol->GetLogicalVolume(), + false, 0, true); + +} + +void DetectorConstruction::Assert(G4VPhysicalVolume* ptr, const G4String& origin) const +{ + if (!ptr) { + G4ExceptionDescription ed; + ed << " Description : World physical volume is not defined."; + G4Exception(origin, "nullptr", FatalException, ed); + } +} + + + + + // // REFLECTIVE FOIL ///////////////////////////////////////////// + // + // //there are three reflective foils, one in the bottom, two in the sides + // //first create the logical volume and then create the surface with LAr + // const G4String bottom_foil_name = "REF_BOTTOM_FOIL"; + // const G4String side_foil_name = "REF_SIDE_FOIL"; + // + // //for the moment, same area as WLS plate. Should be checked. + // //thickness from manufacturer, should be checked. + // //https://www.isoltronic.ch/assets/of-m-vikuiti-esr-app-guide.pdf + // const G4double bottom_foil_width = 112.0*mm; // X + // const G4double bottom_foil_thickn = 0.165*mm; // Y + // const G4double bottom_foil_length = 491.5*mm; // Z + // + // const G4double side_foil_width = 112.0*mm; // X + // const G4double side_foil_thickn = 4.0*mm; + // const G4double side_foil_length = 0.165*mm; // Y + // + // G4Box* bottom_foil_solid_vol = + // new G4Box(bottom_foil_name, bottom_foil_width/2., + // bottom_foil_thickn/2., bottom_foil_length/2.); + // + // G4Box* side_foil_solid_vol = + // new G4Box(side_foil_name, side_foil_width/2., + // side_foil_thickn/2., side_foil_length/2.); + // + // G4Material* plexiglass = G4NistManager::Instance()->FindOrBuildMaterial("G4_PLEXIGLASS"); + // + // G4LogicalVolume* bottom_foil_logic_vol = new + // G4LogicalVolume(bottom_foil_solid_vol, plexiglass, bottom_foil_name); + // + // G4LogicalVolume* side_foil_logic_vol = new + // G4LogicalVolume(side_foil_solid_vol, plexiglass, side_foil_name); + // + // //the position is still unknown. For the moment put 3mm, should be checked + // const G4double bottom_foil_posy = -bottom_foil_thickn/2. - 3*mm; + // const G4double side_foil_posz = bottom_foil_length/2 + 3*mm; + // + // G4VPhysicalVolume* bottom_foil_phys_vol = + // new G4PVPlacement(nullptr, G4ThreeVector(0.,bottom_foil_posy,0.), + // bottom_foil_logic_vol, bottom_foil_name, world_logic_vol, + // false, 0, true); + // + // G4VPhysicalVolume* side1_foil_phys_vol = + // new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,side_foil_posz), + // side_foil_logic_vol, side_foil_name, world_logic_vol, + // false, 0, true); + // + // G4VPhysicalVolume* side2_foil_phys_vol = + // new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,-side_foil_posz), + // side_foil_logic_vol, side_foil_name, world_logic_vol, + // false, 0, true); + // + // //now create the surface + // const G4String refsurf_name = "REF_SURFACE"; + // G4OpticalSurface* refsurf_opsurf = new + // G4OpticalSurface(refsurf_name, unified, polishedfrontpainted, dielectric_dielectric, 1); + // + // refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::VIKUITI()); + // new G4LogicalSkinSurface("REF_FOIL_SURFACE",bottom_foil_logic_vol,refsurf_opsurf); + // new G4LogicalSkinSurface("REF_FOIL_SURFACE",side_foil_logic_vol,refsurf_opsurf); + // + // + // // DICHROIC FILTER /////////////////////////////////////// + // + // const G4String filter_name = "DICHROIC_FILTER"; + // + // //dimensions should be checked + // const G4double filter_width = 112.0*mm; // X + // const G4double filter_thickn = 1.0*mm; // Y + // const G4double filter_length = 491.5*mm; // Z + // + // G4Box* filter_solid_vol = + // new G4Box(filter_name, filter_width/2., filter_thickn/2., filter_length/2.); + // + // G4Material* filter_material = G4NistManager::Instance()->FindOrBuildMaterial("G4_SILICON_DIOXIDE"); + // filter_material->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); + // + // G4LogicalVolume* filter_logic_vol = + // new G4LogicalVolume(filter_solid_vol,filter_material,filter_name); + // + // G4double filter_wls_gap = 1.5*mm; + // G4double filter_ypos = plate_thickn/2. + filter_wls_gap + filter_thickn/2.; + // + // G4VPhysicalVolume* filter_phys_vol = + // new G4PVPlacement(nullptr, G4ThreeVector(0.,filter_ypos,0.), + // filter_logic_vol, filter_name, world_logic_vol, + // false, 0, true); + // + // //now create the surface. It is a dichroic filter, so we need more things + // setenv("G4DICHROICDATA","data/dichroic_data",1); + // const G4String filtersurf_name = "FILTER_SURFACE"; + // G4OpticalSurface* filtersurf_opsurf = new + // G4OpticalSurface(filtersurf_name, dichroic, polished, dielectric_dichroic, 1); + // + // refsurf_opsurf->SetMaterialPropertiesTable(OpticalMaterialProperties::FusedSilica()); + // new G4LogicalSkinSurface("FILTER_SURFACE",filter_logic_vol,filtersurf_opsurf); + // + // + // // PTP LAYER /////////////////////////////////////// + // + // const G4String ptp_name = "PTP_LAYER"; + // + // //real dimensions not known. To be reviewed + // const G4double ptp_width = 112.0*mm; // X + // const G4double ptp_thickn = 0.002*mm; // Y from https://arxiv.org/pdf/1912.09191.pdf + // const G4double ptp_length = 491.5*mm; // Z + // + // G4Box* ptp_solid_vol = + // new G4Box(ptp_name, ptp_width/2., ptp_thickn/2., ptp_length/2.); + // + // G4Material* ptp = G4NistManager::Instance()->FindOrBuildMaterial("G4_TERPHENYL"); + // ptp->SetMaterialPropertiesTable(OpticalMaterialProperties::PTP()); + // + // G4LogicalVolume* ptp_logic_vol = new G4LogicalVolume(ptp_solid_vol, ptp, ptp_name); + // + // //position still not known, to be reviewed + // const G4double ptp_posy = plate_thickn/2 + filter_wls_gap + filter_thickn + ptp_thickn/2; + // //1.5*mm from https://indico.fnal.gov/event/45283/contributions/195721/attachments/133823/165234/X-ARAPUCA_Cuts_Study.pdf + // //and 1*mm is the assumed filter width + // G4VPhysicalVolume* ptp_phys_vol = + // new G4PVPlacement(nullptr, G4ThreeVector(0.,ptp_posy,0.), + // ptp_logic_vol, ptp_name, world_logic_vol, + // false, 0, true); + // + // + // ////////////////////////////////////////////////////////// + // // describing photosensors position + // const G4int nsens = 48; //n SiPMs per arapuca cell + // G4int nsides = 2; // n sides of the sell covered with SiPMs + // const G4double sensSpacing = plate_length/(nsens/2)*mm; + // G4double firstPosition = sensSpacing/2; + // G4ThreeVector SiPMpos; + // G4double Zpos; + // G4double Xpos; + // G4int nsensStart; + // G4String SiPM_name = "SiPM"; + // G4RotationMatrix* rotationMatrix = new G4RotationMatrix(); + // + // G4VPhysicalVolume* photosens_phys_vol[nsens]; + // + // for (G4int iside = 0; iside < nsides; iside ++) + // { + // if(iside == 0) + // { + // Xpos = plate_width/2 + 6*mm; + // rotationMatrix->rotateY(90.*deg); //they have to be facing de plate + // nsensStart = 0; + // //rotationMatrix->print(std::cout); + // + // } + // else if(iside == 1) + // { + // Xpos = -(plate_width/2+6*mm); + // rotationMatrix->rotateY(180.*deg); //rotations apply sequently, so this rotation is 90º+180º + // nsensStart = nsens/2; + // //rotationMatrix->print(std::cout); + // } + // for (G4int isens = 0; isens < nsens/2; isens++) + // { + // Zpos = -plate_length/2 + firstPosition + isens * sensSpacing; + // // Not sure if the sensors shuld be plased in the world logic word + // // When putting them in plate_logic_vol it produces overlap errors + // photosens_phys_vol[nsensStart+isens]= + // new G4PVPlacement(rotationMatrix,G4ThreeVector(Xpos,0.,Zpos), + // this->GenericPhotosensor(),SiPM_name, world_logic_vol, + // false, isens, true); + // } + // } + // + // return world_phys_vol; +// } + + + + + + + + +// G4LogicalVolume* DetectorConstruction::GenericPhotosensor() +// { +// // PHOTOSENSOR ENCASING ////////////////////////////////// +// +// const G4double width_ = 6.*mm; +// const G4double height_ = 6.*mm; +// const G4double thickness_ = 2.*mm; +// +// G4String name = "PHOTOSENSOR"; +// +// G4Box* encasing_solid_vol = +// new G4Box(name, width_/2., height_/2., thickness_/2.); +// +// G4LogicalVolume* encasing_logic_vol = +// new G4LogicalVolume(encasing_solid_vol, Materials::FR4(), name); +// +// // OPTICAL WINDOW //////////////////////////////////////// +// +// name = "PHOTOSENSOR_WINDOW"; +// +// G4double window_thickness = thickness_/4.; +// +// G4Box* window_solid_vol = +// new G4Box(name, width_/2., height_/2., window_thickness/2.); +// +// +// G4LogicalVolume* window_logic_vol = +// new G4LogicalVolume(window_solid_vol, +// Materials::OpticalSilicone(), +// name); +// +// //THE SIPM IS CREATED FACING +Z DIRECTION +// G4double zpos = thickness_/2. - window_thickness/2.; +// +// new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), +// window_logic_vol, name, encasing_logic_vol, +// false, 0, false); +// +// // PHOTOSENSITIVE AREA ///////////////////////////////////////////// +// +// name = "PHOTOSENSOR_SENSAREA"; +// +// G4double sensarea_thickness = 0.1*mm; +// +// G4Box* sensarea_solid_vol = +// new G4Box(name, width_/2., height_/2., sensarea_thickness/2.); +// +// G4LogicalVolume* sensarea_logic_vol = +// new G4LogicalVolume(sensarea_solid_vol, +// G4NistManager::Instance()->FindOrBuildMaterial("G4_Si"), +// name); +// +// zpos = thickness_/2. - window_thickness - sensarea_thickness/2.; +// +// new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), +// sensarea_logic_vol, name, encasing_logic_vol, +// false, 0, false); +// +// // OPTICAL PROPERTIES ////////////////////////////////////////////// +// +// name = "PHOTOSENSOR_OPSURF"; +// +// G4double energy[] = {0.2*eV, 11.5*eV}; +// G4double reflectivity[] = {0.0 , 0.0}; +// G4double efficiency[] = {1.0 , 1.0}; +// +// G4MaterialPropertiesTable* photosensor_mpt = new G4MaterialPropertiesTable(); +// photosensor_mpt->AddProperty("REFLECTIVITY", energy, reflectivity, 2); +// photosensor_mpt->AddProperty("EFFICIENCY", energy, efficiency, 2); +// +// G4OpticalSurface* photosensor_opsurf = +// new G4OpticalSurface(name, unified, polished, dielectric_metal); +// photosensor_opsurf->SetMaterialPropertiesTable(photosensor_mpt); +// new G4LogicalSkinSurface(name, sensarea_logic_vol, photosensor_opsurf); +// +// // SENSITIVE DETECTOR ////////////////////////////////////////////// +// +// OpticalSD* sensdet = new OpticalSD("/GENERIC_PHOTOSENSOR/SiPM"); +// G4SDManager::GetSDMpointer()->AddNewDetector(sensdet); +// sensarea_logic_vol->SetSensitiveDetector(sensdet); +// +// //////////////////////////////////////////////////////////////////// +// +// return encasing_logic_vol; +// } diff --git a/src/DetectorConstruction.h b/src/DetectorConstruction.h index 003fb19..0f6599c 100644 --- a/src/DetectorConstruction.h +++ b/src/DetectorConstruction.h @@ -1,38 +1,43 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | DetectorConstruction.h -// -// * Author: -// * Creation date: 10 February 2020 -// ----------------------------------------------------------------------------- - -#ifndef DETECTOR_CONSTRUCTION_H -#define DETECTOR_CONSTRUCTION_H - -#include - -class G4Material; -class G4LogicalVolume; - - -class DetectorConstruction: public G4VUserDetectorConstruction -{ -public: - DetectorConstruction(); - ~DetectorConstruction(); - G4VPhysicalVolume* Construct() override; -private: - void ConstructWorld(G4VPhysicalVolume&); - void ConstructWLSPlate(G4VPhysicalVolume*) const; - void ConstructPhotosensors(G4VPhysicalVolume*) const; - void ConstructReflectiveFoils(G4VPhysicalVolume*) const; - - void Assert(G4VPhysicalVolume*, const G4String&) const; - -private: - const G4double world_size_; - const G4double plate_width_, plate_thickn_, plate_length_; - const G4double foil_thickn_; - const G4int num_phsensors; -}; - -#endif +// ----------------------------------------------------------------------------- +// G4OpSim | DetectorConstruction.h +// +// * Author: +// * Creation date: 10 February 2020 +// ----------------------------------------------------------------------------- + +#ifndef DETECTOR_CONSTRUCTION_H +#define DETECTOR_CONSTRUCTION_H + +#include + +class G4Material; +class G4LogicalVolume; + + +class DetectorConstruction: public G4VUserDetectorConstruction +{ +public: + DetectorConstruction(); + ~DetectorConstruction(); + G4VPhysicalVolume* Construct() override; +private: + void ConstructWorld(G4VPhysicalVolume&); + void ConstructWLSPlate(G4VPhysicalVolume*) const; + void ConstructPhotosensors(G4VPhysicalVolume*) const; + void ConstructReflectiveFoils(G4VPhysicalVolume*) const; + void ConstructDichroicFilter(G4VPhysicalVolume*) const; + void ConstructPTPLayer(G4VPhysicalVolume*) const; + + void Assert(G4VPhysicalVolume*, const G4String&) const; + +private: + const G4double world_size_; + const G4double plate_width_, plate_thickn_, plate_length_; + //const G4double foil_thickn_; + const G4double foil_thickn_, foil_gap_; + const G4double filter_thickn_, filter_gap_; + const G4double ptp_thickn_; + const G4int num_phsensors; +}; + +#endif diff --git a/src/EventAction.cpp b/src/EventAction.cpp index 6f24841..cddd93d 100644 --- a/src/EventAction.cpp +++ b/src/EventAction.cpp @@ -1,18 +1,40 @@ -// ----------------------------------------------------------------------------- -// G4Basic | EventAction.cpp -// -// User run action class. -// ----------------------------------------------------------------------------- - -#include "EventAction.h" - -#include - - -void EventAction::BeginOfEventAction(const G4Event*) -{ -} - -void EventAction::EndOfEventAction(const G4Event*) -{ -} +// ----------------------------------------------------------------------------- +// G4Basic | EventAction.cpp +// +// User run action class. +// ----------------------------------------------------------------------------- + +//#include "params.h" +#include "EventAction.h" +#include "Analysis.h" + +#include +#include +#include "G4Event.hh" +#include "G4UnitsTable.hh" +#include +#include + +using CLHEP::h_Planck; +using CLHEP::c_light; + + +EventAction::EventAction() + : G4UserEventAction(), + fInitEnerg(6 * eV) + +{} + +EventAction::~EventAction() +{} + +void EventAction::BeginOfEventAction(const G4Event*) +{ +} + +void EventAction::EndOfEventAction(const G4Event* anEvent) +{ + //G4PrimaryVertex* primaryVertex = anEvent->GetPrimaryVertex(); + //G4double initEnerg = primaryVertex->GetPrimary()->GetTotalEnergy(); + //G4cout << G4BestUnit((h_Planck * c_light) / (initEnerg), "Length") << G4endl; +} diff --git a/src/EventAction.h b/src/EventAction.h index 0e5ac5c..49afcad 100644 --- a/src/EventAction.h +++ b/src/EventAction.h @@ -1,27 +1,32 @@ -// ----------------------------------------------------------------------------- -// G4Basic | EventAction.h -// -// User run action class. -// ----------------------------------------------------------------------------- - -#ifndef EVENT_ACTION_H -#define EVENT_ACTION_H - -#include - -class G4Event; - - -class EventAction: public G4UserEventAction -{ -public: - EventAction(); - virtual ~EventAction(); - virtual void BeginOfEventAction(const G4Event*); - virtual void EndOfEventAction(const G4Event*); -}; - -inline EventAction::EventAction() {} -inline EventAction::~EventAction() {} - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | EventAction.h +// +// User run action class. +// ----------------------------------------------------------------------------- + +#ifndef EVENT_ACTION_H +#define EVENT_ACTION_H + +#include +#include + +class G4Event; +class G4Step; + + +class EventAction: public G4UserEventAction +{ +public: + EventAction(); + virtual ~EventAction(); + virtual void BeginOfEventAction(const G4Event*); + virtual void EndOfEventAction(const G4Event*); + +private: + G4double fInitEnerg; +}; + +//inline EventAction::EventAction() {} +//inline EventAction::~EventAction() {} + +#endif diff --git a/src/GenericPhotosensor.cpp b/src/GenericPhotosensor.cpp index a13d186..0f433c6 100644 --- a/src/GenericPhotosensor.cpp +++ b/src/GenericPhotosensor.cpp @@ -1,119 +1,119 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | GenericPhotosensor.cpp -// -// * Author: -// * Creation date: 10 February 2020 -// ----------------------------------------------------------------------------- - -#include "GenericPhotosensor.h" -#include "OpticalMaterialProperties.h" - -#include "Materials.h" -#include "OpticalSD.h" - -#include -#include -#include -#include -#include -#include -#include -#include - - -GenericPhotosensor::GenericPhotosensor(): - width_ (6.0*mm), - height_ (6.0*mm), - thickness_(2.0*mm), - logvol_(nullptr) -{ -} - - -GenericPhotosensor::~GenericPhotosensor() -{ -} - - -void GenericPhotosensor::Construct() -{ - // PHOTOSENSOR ENCASING ////////////////////////////////// - - G4String name = "PHOTOSENSOR"; - - G4Box* encasing_solid_vol = - new G4Box(name, width_/2., height_/2., thickness_/2.); - - G4LogicalVolume* encasing_logic_vol = - new G4LogicalVolume(encasing_solid_vol, Materials::FR4(), name); - - logvol_ = encasing_logic_vol; - - // OPTICAL WINDOW //////////////////////////////////////// - - name = "PHOTOSENSOR_WINDOW"; - - G4double window_thickness = thickness_/4.; - - G4Material* op_sil = Materials::OpticalSilicone(); - op_sil->SetMaterialPropertiesTable(OpticalMaterialProperties::GlassEpoxy()); - - G4Box* window_solid_vol = - new G4Box(name, width_/2., height_/2., window_thickness/2.); - - G4LogicalVolume* window_logic_vol = - new G4LogicalVolume(window_solid_vol, - op_sil, - name); - - G4double zpos = thickness_/2. - window_thickness/2.; - - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), - window_logic_vol, name, encasing_logic_vol, - false, 0, true); - - // PHOTOSENSITIVE AREA ///////////////////////////////////////////// - - name = "PHOTOSENSOR_SENSAREA"; - - G4double sensarea_thickness = 0.1*mm; - - G4Box* sensarea_solid_vol = - new G4Box(name, width_/2., height_/2., sensarea_thickness/2.); - - G4LogicalVolume* sensarea_logic_vol = - new G4LogicalVolume(sensarea_solid_vol, - G4NistManager::Instance()->FindOrBuildMaterial("G4_Si"), - name); - - zpos = thickness_/2. - window_thickness - sensarea_thickness/2.; - - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), - sensarea_logic_vol, name, encasing_logic_vol, - false, 0, false); - - // OPTICAL PROPERTIES ////////////////////////////////////////////// - - name = "PHOTOSENSOR_OPSURF"; - - G4double energy[] = {0.2*eV, 11.5*eV}; - G4double reflectivity[] = {0.0 , 0.0}; - G4double efficiency[] = {1.0 , 1.0}; - - G4MaterialPropertiesTable* photosensor_mpt = new G4MaterialPropertiesTable(); - photosensor_mpt->AddProperty("REFLECTIVITY", energy, reflectivity, 2); - photosensor_mpt->AddProperty("EFFICIENCY", energy, efficiency, 2); - - G4OpticalSurface* photosensor_opsurf = - new G4OpticalSurface(name, unified, polished, dielectric_metal); - photosensor_opsurf->SetMaterialPropertiesTable(photosensor_mpt); - new G4LogicalSkinSurface(name, sensarea_logic_vol, photosensor_opsurf); - - // SENSITIVE DETECTOR ////////////////////////////////////////////// - - OpticalSD* sensdet = new OpticalSD("/GENERIC_PHOTOSENSOR/SiPM"); - G4SDManager::GetSDMpointer()->AddNewDetector(sensdet); - sensarea_logic_vol->SetSensitiveDetector(sensdet); - - //////////////////////////////////////////////////////////////////// -} +// ----------------------------------------------------------------------------- +// G4OpSim | GenericPhotosensor.cpp +// +// * Author: +// * Creation date: 10 February 2020 +// ----------------------------------------------------------------------------- + +#include "GenericPhotosensor.h" +#include "OpticalMaterialProperties.h" + +#include "Materials.h" +#include "OpticalSD.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + +GenericPhotosensor::GenericPhotosensor(): + width_ (6.0*mm), + height_ (6.0*mm), + thickness_(2.0*mm), + logvol_(nullptr) +{ +} + + +GenericPhotosensor::~GenericPhotosensor() +{ +} + + +void GenericPhotosensor::Construct() +{ + // PHOTOSENSOR ENCASING ////////////////////////////////// + + G4String name = "PHOTOSENSOR"; + + G4Box* encasing_solid_vol = + new G4Box(name, width_/2., height_/2., thickness_/2.); + + G4LogicalVolume* encasing_logic_vol = + new G4LogicalVolume(encasing_solid_vol, Materials::FR4(), name); + + logvol_ = encasing_logic_vol; + + // OPTICAL WINDOW //////////////////////////////////////// + + name = "PHOTOSENSOR_WINDOW"; + + G4double window_thickness = thickness_/4.; + + G4Material* op_sil = Materials::OpticalSilicone(); + op_sil->SetMaterialPropertiesTable(OpticalMaterialProperties::GlassEpoxy()); + + G4Box* window_solid_vol = + new G4Box(name, width_/2., height_/2., window_thickness/2.); + + G4LogicalVolume* window_logic_vol = + new G4LogicalVolume(window_solid_vol, + op_sil, + name); + + G4double zpos = thickness_/2. - window_thickness/2.; + + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), + window_logic_vol, name, encasing_logic_vol, + false, 0, true); + + // PHOTOSENSITIVE AREA ///////////////////////////////////////////// + + name = "PHOTOSENSOR_SENSAREA"; + + G4double sensarea_thickness = 0.1*mm; + + G4Box* sensarea_solid_vol = + new G4Box(name, width_/2., height_/2., sensarea_thickness/2.); + + G4LogicalVolume* sensarea_logic_vol = + new G4LogicalVolume(sensarea_solid_vol, + G4NistManager::Instance()->FindOrBuildMaterial("G4_Si"), + name); + + zpos = thickness_/2. - window_thickness - sensarea_thickness/2.; + + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., zpos), + sensarea_logic_vol, name, encasing_logic_vol, + false, 0, false); + + // OPTICAL PROPERTIES ////////////////////////////////////////////// + + name = "PHOTOSENSOR_OPSURF"; + + G4double energy[] = {0.2*eV, 11.5*eV}; + G4double reflectivity[] = {0.0 , 0.0}; + G4double efficiency[] = {1.0 , 1.0}; + + G4MaterialPropertiesTable* photosensor_mpt = new G4MaterialPropertiesTable(); + photosensor_mpt->AddProperty("REFLECTIVITY", energy, reflectivity, 2); + photosensor_mpt->AddProperty("EFFICIENCY", energy, efficiency, 2); + + G4OpticalSurface* photosensor_opsurf = + new G4OpticalSurface(name, unified, polished, dielectric_metal); + photosensor_opsurf->SetMaterialPropertiesTable(photosensor_mpt); + new G4LogicalSkinSurface(name, sensarea_logic_vol, photosensor_opsurf); + + // SENSITIVE DETECTOR ////////////////////////////////////////////// + + OpticalSD* sensdet = new OpticalSD("/GENERIC_PHOTOSENSOR/SiPM"); + G4SDManager::GetSDMpointer()->AddNewDetector(sensdet); + sensarea_logic_vol->SetSensitiveDetector(sensdet); + + //////////////////////////////////////////////////////////////////// +} diff --git a/src/GenericPhotosensor.h b/src/GenericPhotosensor.h index 4424e27..cdb0482 100644 --- a/src/GenericPhotosensor.h +++ b/src/GenericPhotosensor.h @@ -1,38 +1,38 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | GenericPhotosensor.h -// -// * Author: -// * Creation date: 10 February 2020 -// ----------------------------------------------------------------------------- - -#ifndef GENERIC_PHOTOSENSOR_H -#define GENERIC_PHOTOSENSOR_H - -#include - -class G4LogicalVolume; - - -class GenericPhotosensor -{ -public: - GenericPhotosensor(); - ~GenericPhotosensor(); - void Construct(); - G4LogicalVolume* GetLogicalVolume(); - - G4double GetWidth() const; - G4double GetHeight() const; - G4double GetThickness() const; - -private: - G4double width_, height_, thickness_; - G4LogicalVolume* logvol_; -}; - -inline G4LogicalVolume* GenericPhotosensor::GetLogicalVolume() { return logvol_; } -inline G4double GenericPhotosensor::GetWidth() const { return width_; } -inline G4double GenericPhotosensor::GetHeight() const { return height_; } -inline G4double GenericPhotosensor::GetThickness() const { return thickness_; } - -#endif +// ----------------------------------------------------------------------------- +// G4OpSim | GenericPhotosensor.h +// +// * Author: +// * Creation date: 10 February 2020 +// ----------------------------------------------------------------------------- + +#ifndef GENERIC_PHOTOSENSOR_H +#define GENERIC_PHOTOSENSOR_H + +#include + +class G4LogicalVolume; + + +class GenericPhotosensor +{ +public: + GenericPhotosensor(); + ~GenericPhotosensor(); + void Construct(); + G4LogicalVolume* GetLogicalVolume(); + + G4double GetWidth() const; + G4double GetHeight() const; + G4double GetThickness() const; + +private: + G4double width_, height_, thickness_; + G4LogicalVolume* logvol_; +}; + +inline G4LogicalVolume* GenericPhotosensor::GetLogicalVolume() { return logvol_; } +inline G4double GenericPhotosensor::GetWidth() const { return width_; } +inline G4double GenericPhotosensor::GetHeight() const { return height_; } +inline G4double GenericPhotosensor::GetThickness() const { return thickness_; } + +#endif diff --git a/src/Materials.cpp b/src/Materials.cpp index e273945..65becaf 100644 --- a/src/Materials.cpp +++ b/src/Materials.cpp @@ -1,94 +1,94 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | Materials.cpp -// -// Definitions of materials used in the geometry. -// * Author: -// * Creation date: 26 March 2020 -// ----------------------------------------------------------------------------- - -#include "Materials.h" - -#include -#include -#include - - -G4Material* Materials::FR4() -{ - // FR-4 is a composite material widely used for printed circuits boards. - // It consists of woven fiberglass cloth with an epoxy resin binder that is - // flame resistant. Typical proportions are 60% fused silica and 40% epoxy. - - const G4String name = "FR4"; - G4Material* mat = G4Material::GetMaterial(name, false); - - if (!mat) { - mat = new G4Material(name, 1.850*g/cm3, 2, kStateSolid); - mat->AddMaterial(Materials::FusedSilica(), 0.60); - mat->AddMaterial(Materials::Epoxy(), 0.40); - } - - return mat; -} - - -G4Material* Materials::OpticalSilicone() -{ - // Silicone resin with a methyl group - // (https://en.wikipedia.org/wiki/Silicone_resin) - - const G4String name = "OPTICAL_SILICONE"; - - G4Material* mat = G4Material::GetMaterial(name, false); - - if (!mat) { - G4NistManager* nist = G4NistManager::Instance(); - - G4Element* H = nist->FindOrBuildElement("H"); - G4Element* C = nist->FindOrBuildElement("C"); - G4Element* O = nist->FindOrBuildElement("O"); - G4Element* Si = nist->FindOrBuildElement("Si"); - - mat = new G4Material(name, 1.05*g/cm3, 4, kStateSolid); - mat->AddElement(H, 3); - mat->AddElement(C, 1); - mat->AddElement(Si, 1); - mat->AddElement(O, 1); - } - - return mat; -} - - -G4Material* Materials::Epoxy() -{ - // Definition taken from the Geant4 advanced example "Composite Calorimeter" - // (Geant4/examples/advanced/composite_calorimeter/dataglobal/material.cms). - - const G4String name = "EPOXY"; - - G4Material* mat = G4Material::GetMaterial(name, false); - - if (!mat) { - G4NistManager* nist = G4NistManager::Instance(); - - G4Element* H = nist->FindOrBuildElement("H"); - G4Element* C = nist->FindOrBuildElement("C"); - G4Element* O = nist->FindOrBuildElement("O"); - - mat = new G4Material(name, 1.3*g/cm3, 3); - mat->AddElement(H, 44); - mat->AddElement(C, 15); - mat->AddElement(O, 7); - } - - return mat; -} - - -G4Material* Materials::FusedSilica() -{ - G4Material* mat = - G4NistManager::Instance()->FindOrBuildMaterial("G4_SILICON_DIOXIDE"); - return mat; -} +// ----------------------------------------------------------------------------- +// G4OpSim | Materials.cpp +// +// Definitions of materials used in the geometry. +// * Author: +// * Creation date: 26 March 2020 +// ----------------------------------------------------------------------------- + +#include "Materials.h" + +#include +#include +#include + + +G4Material* Materials::FR4() +{ + // FR-4 is a composite material widely used for printed circuits boards. + // It consists of woven fiberglass cloth with an epoxy resin binder that is + // flame resistant. Typical proportions are 60% fused silica and 40% epoxy. + + const G4String name = "FR4"; + G4Material* mat = G4Material::GetMaterial(name, false); + + if (!mat) { + mat = new G4Material(name, 1.850*g/cm3, 2, kStateSolid); + mat->AddMaterial(Materials::FusedSilica(), 0.60); + mat->AddMaterial(Materials::Epoxy(), 0.40); + } + + return mat; +} + + +G4Material* Materials::OpticalSilicone() +{ + // Silicone resin with a methyl group + // (https://en.wikipedia.org/wiki/Silicone_resin) + + const G4String name = "OPTICAL_SILICONE"; + + G4Material* mat = G4Material::GetMaterial(name, false); + + if (!mat) { + G4NistManager* nist = G4NistManager::Instance(); + + G4Element* H = nist->FindOrBuildElement("H"); + G4Element* C = nist->FindOrBuildElement("C"); + G4Element* O = nist->FindOrBuildElement("O"); + G4Element* Si = nist->FindOrBuildElement("Si"); + + mat = new G4Material(name, 1.05*g/cm3, 4, kStateSolid); + mat->AddElement(H, 3); + mat->AddElement(C, 1); + mat->AddElement(Si, 1); + mat->AddElement(O, 1); + } + + return mat; +} + + +G4Material* Materials::Epoxy() +{ + // Definition taken from the Geant4 advanced example "Composite Calorimeter" + // (Geant4/examples/advanced/composite_calorimeter/dataglobal/material.cms). + + const G4String name = "EPOXY"; + + G4Material* mat = G4Material::GetMaterial(name, false); + + if (!mat) { + G4NistManager* nist = G4NistManager::Instance(); + + G4Element* H = nist->FindOrBuildElement("H"); + G4Element* C = nist->FindOrBuildElement("C"); + G4Element* O = nist->FindOrBuildElement("O"); + + mat = new G4Material(name, 1.3*g/cm3, 3); + mat->AddElement(H, 44); + mat->AddElement(C, 15); + mat->AddElement(O, 7); + } + + return mat; +} + + +G4Material* Materials::FusedSilica() +{ + G4Material* mat = + G4NistManager::Instance()->FindOrBuildMaterial("G4_SILICON_DIOXIDE"); + return mat; +} diff --git a/src/Materials.h b/src/Materials.h index 4a08d2e..5ec4650 100644 --- a/src/Materials.h +++ b/src/Materials.h @@ -1,27 +1,27 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | Materials.h -// -// Definitions of materials used in the geometry. -// * Author: -// * Creation date: 26 March 2020 -// ----------------------------------------------------------------------------- - -#ifndef MATERIALS_H_ -#define MATERIALS_H_ - -class G4Material; - - -namespace Materials { - - G4Material* FR4(); - - G4Material* OpticalSilicone(); - - G4Material* Epoxy(); - - G4Material* FusedSilica(); - -} // end namespace Materials - -#endif +// ----------------------------------------------------------------------------- +// G4OpSim | Materials.h +// +// Definitions of materials used in the geometry. +// * Author: +// * Creation date: 26 March 2020 +// ----------------------------------------------------------------------------- + +#ifndef MATERIALS_H_ +#define MATERIALS_H_ + +class G4Material; + + +namespace Materials { + + G4Material* FR4(); + + G4Material* OpticalSilicone(); + + G4Material* Epoxy(); + + G4Material* FusedSilica(); + +} // end namespace Materials + +#endif diff --git a/src/OpticalHit.cpp b/src/OpticalHit.cpp index 2ff8bb1..148b9ff 100644 --- a/src/OpticalHit.cpp +++ b/src/OpticalHit.cpp @@ -1,66 +1,66 @@ -// ----------------------------------------------------------------------------- -// G4Basic | OpticalHit.h -// -// TODO: Class description -// * Author: -// * Creation date: 5 Nov 2020 -// ----------------------------------------------------------------------------- - -#include "OpticalHit.h" - - -G4Allocator OpticalHitAllocator; - - -OpticalHit::OpticalHit(): - G4VHit(), - sensor_id_(-1), - time_bin_width_(0.) -{ -} - - -OpticalHit::~OpticalHit() -{ -} - - -OpticalHit::OpticalHit(const OpticalHit& other): G4VHit() -{ - *this = other; -} - - -const OpticalHit& OpticalHit::operator=(const OpticalHit& other) -{ - sensor_id_ = other.sensor_id_; - time_bin_width_ = other.time_bin_width_; - wvf_ = other.wvf_; - - return *this; -} - - -G4int OpticalHit::operator==(const OpticalHit& other) const -{ - return (this == &other) ? 1 : 0; -} - - -void OpticalHit::SetTimeBinWidth(G4double bin_size) -{ - if (wvf_.size() == 0) { - time_bin_width_ = bin_size; - } - else { - G4String msg = "A OpticalHit cannot be rebinned once it has been filled."; - G4Exception("[OpticalHit]", "SetTimeBinWidth()", JustWarning, msg); - } -} - - -void OpticalHit::Fill(G4double time, G4int counts) -{ - G4double time_bin = floor(time/time_bin_width_) * time_bin_width_; - wvf_[time_bin] += counts; -} +// ----------------------------------------------------------------------------- +// G4Basic | OpticalHit.h +// +// TODO: Class description +// * Author: +// * Creation date: 5 Nov 2020 +// ----------------------------------------------------------------------------- + +#include "OpticalHit.h" + + +G4Allocator OpticalHitAllocator; + + +OpticalHit::OpticalHit(): + G4VHit(), + sensor_id_(-1), + time_bin_width_(0.) +{ +} + + +OpticalHit::~OpticalHit() +{ +} + + +OpticalHit::OpticalHit(const OpticalHit& other): G4VHit() +{ + *this = other; +} + + +const OpticalHit& OpticalHit::operator=(const OpticalHit& other) +{ + sensor_id_ = other.sensor_id_; + time_bin_width_ = other.time_bin_width_; + wvf_ = other.wvf_; + + return *this; +} + + +G4int OpticalHit::operator==(const OpticalHit& other) const +{ + return (this == &other) ? 1 : 0; +} + + +void OpticalHit::SetTimeBinWidth(G4double bin_size) +{ + if (wvf_.size() == 0) { + time_bin_width_ = bin_size; + } + else { + G4String msg = "A OpticalHit cannot be rebinned once it has been filled."; + G4Exception("[OpticalHit]", "SetTimeBinWidth()", JustWarning, msg); + } +} + + +void OpticalHit::Fill(G4double time, G4int counts) +{ + G4double time_bin = floor(time/time_bin_width_) * time_bin_width_; + wvf_[time_bin] += counts; +} diff --git a/src/OpticalHit.h b/src/OpticalHit.h index c71e512..f0dc3d7 100644 --- a/src/OpticalHit.h +++ b/src/OpticalHit.h @@ -1,67 +1,67 @@ -// ----------------------------------------------------------------------------- -// G4Basic | OpticalHit.h -// -// TODO: Class description -// * Author: -// * Creation date: 5 Nov 2020 -// ----------------------------------------------------------------------------- - -#ifndef OPTICAL_HIT_H -#define OPTICAL_HIT_H - -#include -#include -#include -#include - - - -class OpticalHit: public G4VHit -{ -public: - OpticalHit(); - OpticalHit(const OpticalHit&); - ~OpticalHit(); - - const OpticalHit& operator=(const OpticalHit&); - G4int operator==(const OpticalHit&) const; - - void* operator new(size_t); - void operator delete(void*); - - G4int GetSensorID() const; - void SetSensorID(G4int); - - G4double GetTimeBinWidth() const; - void SetTimeBinWidth(G4double); - - void Fill(G4double time, G4int counts=1); - - const std::map& GetWaveform() const; - -private: - G4int sensor_id_; - G4double time_bin_width_; - std::map wvf_; -}; - -////////////////////////////////////////////////////////////////////// - -typedef G4THitsCollection OpticalHitCollection; -extern G4Allocator OpticalHitAllocator; - -inline void* OpticalHit::operator new(size_t) -{ return ((void*) OpticalHitAllocator.MallocSingle()); } - -inline void OpticalHit::operator delete(void* hit) -{ OpticalHitAllocator.FreeSingle((OpticalHit*) hit); } - -inline G4int OpticalHit::GetSensorID() const { return sensor_id_; } -inline void OpticalHit::SetSensorID(G4int id) { sensor_id_ = id; } - -inline G4double OpticalHit::GetTimeBinWidth() const { return time_bin_width_; } - -inline const std::map& OpticalHit::GetWaveform() const -{ return wvf_; } - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | OpticalHit.h +// +// TODO: Class description +// * Author: +// * Creation date: 5 Nov 2020 +// ----------------------------------------------------------------------------- + +#ifndef OPTICAL_HIT_H +#define OPTICAL_HIT_H + +#include +#include +#include +#include + + + +class OpticalHit: public G4VHit +{ +public: + OpticalHit(); + OpticalHit(const OpticalHit&); + ~OpticalHit(); + + const OpticalHit& operator=(const OpticalHit&); + G4int operator==(const OpticalHit&) const; + + void* operator new(size_t); + void operator delete(void*); + + G4int GetSensorID() const; + void SetSensorID(G4int); + + G4double GetTimeBinWidth() const; + void SetTimeBinWidth(G4double); + + void Fill(G4double time, G4int counts=1); + + const std::map& GetWaveform() const; + +private: + G4int sensor_id_; + G4double time_bin_width_; + std::map wvf_; +}; + +////////////////////////////////////////////////////////////////////// + +typedef G4THitsCollection OpticalHitCollection; +extern G4Allocator OpticalHitAllocator; + +inline void* OpticalHit::operator new(size_t) +{ return ((void*) OpticalHitAllocator.MallocSingle()); } + +inline void OpticalHit::operator delete(void* hit) +{ OpticalHitAllocator.FreeSingle((OpticalHit*) hit); } + +inline G4int OpticalHit::GetSensorID() const { return sensor_id_; } +inline void OpticalHit::SetSensorID(G4int id) { sensor_id_ = id; } + +inline G4double OpticalHit::GetTimeBinWidth() const { return time_bin_width_; } + +inline const std::map& OpticalHit::GetWaveform() const +{ return wvf_; } + +#endif diff --git a/src/OpticalMaterialProperties.cpp b/src/OpticalMaterialProperties.cpp index 436743d..607357f 100644 --- a/src/OpticalMaterialProperties.cpp +++ b/src/OpticalMaterialProperties.cpp @@ -1,361 +1,365 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | OpticalMaterialProperties.cpp -// -// * Author: -// * Creation date: 26 March 2020 -// ----------------------------------------------------------------------------- - -#include "OpticalMaterialProperties.h" - -#include - -#include - -#include - -using CLHEP::h_Planck; -using CLHEP::c_light; - - -G4MaterialPropertiesTable* OpticalMaterialProperties::Vacuum() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - - G4double energies[] = {OpticalMaterialProperties::energy_min, - OpticalMaterialProperties::energy_max}; - - // Refractive index (RINDEX) - G4double rindex[] = {1., 1.}; - mpt->AddProperty("RINDEX", energies, rindex, 2); - - // Absorption length (ABSLENGTH) - G4double abslength[] = {OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max}; - mpt->AddProperty("ABSLENGTH", energies, abslength, 2); - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::LAr() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - - G4double energy_step = 0.05 * eV; - G4double energy_curr = OpticalMaterialProperties::energy_min; - std::vector energies; - - do { - energies.push_back(energy_curr); - energy_curr += energy_step; - } while(energy_curr < OpticalMaterialProperties::energy_max); - - // Refractive index (RINDEX) - // https://arxiv.org/pdf/2002.09346.pdf - G4double a0 = 0.335; - G4double aUV = 0.099; - G4double aIR = 0.008; - G4double lambdaUV = 106.6 * nm; - G4double lambdaIR = 908.3 * nm; - std::vector rindex; - - for (auto energy: energies) { - - G4double wavelength = h_Planck * c_light / energy; - G4double x = a0 + aUV * pow (wavelength, 2) / (pow (wavelength, 2) - pow (lambdaUV, 2)) + aIR * pow (wavelength, 2) / (pow (wavelength, 2) - pow (lambdaIR, 2)); - G4double rindex_element = sqrt (1 + 3 * x / (3 - x)); - rindex.push_back(rindex_element); - } - - mpt->AddProperty("RINDEX", energies.data(), rindex.data(), energies.size()); - - // Absorption length (ABSLENGTH) - G4double energies_lim[] = {OpticalMaterialProperties::energy_min, - OpticalMaterialProperties::energy_max}; - G4double abslength[] = {OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max}; - mpt->AddProperty("ABSLENGTH", energies_lim, abslength, 2); - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::PVT() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - - G4double energy_step = 0.05 * eV; - G4double energy_curr = OpticalMaterialProperties::energy_min; - std::vector energies; - - do { - energies.push_back(energy_curr); - energy_curr += energy_step; - } while(energy_curr < 7.5 * eV); // Put 7.5 eV instead of energy_max = 11.3 eV in order to avoid divergence. - - // Refractive index (RINDEX) - // The refractive index is fitted to a Sellmeier function: https://en.wikipedia.org/wiki/Sellmeier_equation - G4double A = 1.421; - G4double B1 = 0.9944; - G4double C1 = 26250 * pow (nm, 2); - std::vector rindex; - - for (auto energy: energies) { - - G4double wavelength = h_Planck * c_light / energy; - G4double rindex_element = sqrt (A + B1 * pow (wavelength, 2) / (pow (wavelength, 2) - C1));; - rindex.push_back(rindex_element); - } - - mpt->AddProperty("RINDEX", energies.data(), rindex.data(), energies.size()); - - //Absorption length (ABSLENGTH) - //assuming just WLS can absorve - G4double energies_lim[] = {OpticalMaterialProperties::energy_min, - OpticalMaterialProperties::energy_max}; - G4double abslength[] = {OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max}; - mpt->AddProperty("ABSLENGTH", energies_lim, abslength, 2); - - //WLS absorption WLSABSLENGTH - //assuming maximum absorption length for extremal value - G4double WLS_abslength[] = {OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max, - 879.673*mm ,406.487*mm ,296.786*mm,165.772*mm, - 118.88*mm ,96.5343*mm ,70.4532*mm,38.5933*mm, - 26.714 *mm ,16.9549*mm ,11.2558*mm,9.49174*mm, - 8.78259*mm ,7.31249*mm ,6.59717*mm,6.05747*mm, - 5.24668*mm ,4.38833*mm ,3.73978*mm,2.71100*mm, - 2.1048*mm ,1.97013*mm ,1.20762*mm,0.71176*mm, - 0.803417*mm,0.803417*mm,0.97114*mm,1.63289*mm}; - G4double WLS_abs_energy[] = {OpticalMaterialProperties::energy_min, - 2.90388, - 2.9723 *eV,2.98595*eV,2.9987 *eV,3.00867*eV, - 3.01317*eV,3.01997*eV,3.02443*eV,3.03555*eV, - 3.04441*eV,3.0508*eV ,3.05709*eV,3.06822*eV, - 3.07265*eV,3.07916*eV,3.08353*eV,3.08565*eV, - 3.09449*eV,3.10326*eV,3.11208*eV,3.12862*eV, - 3.14538*eV,3.15101*eV,3.20314*eV,3.25895*eV, - 3.27264*eV,3.88859*eV,3.91598*eV,4.01202*eV}; - G4int WLS_abs_entries = 30; - mpt->AddProperty("WLSABSLENGTH", WLS_abs_energy, WLS_abslength, WLS_abs_entries); - - // Emision spectrum (WLSCOMPONENT) - // from https://arxiv.org/pdf/1912.09191.pdf EJ286 - G4double WLS_emi_energy[] = {2.32899*eV,2.38369*eV,2.46565*eV,2.57017*eV,2.61099*eV, - 2.64714*eV,2.69664*eV,2.72479*eV,2.77871*eV,2.79583*eV, - 2.84050*eV,2.85033*eV,2.86809*eV,2.89408*eV,2.91221*eV, - 2.94121*eV,2.97986*eV,2.98148*eV,3.00011*eV,3.01669*eV, - 3.01794*eV,3.03617*eV,3.04968*eV,3.04979*eV,3.06040*eV, - 3.07350*eV,3.07475*eV,3.08767*eV,3.11679*eV,3.12031*eV}; - G4double WLS_emi_Spectrum[] = {0.00293557,0.0122252,0.0439192,0.121378,0.166870, - 0.23135100,0.3592200,0.4271160,0.539329,0.579575, - 0.75539500,0.7947390,0.8865420,0.973837,0.995832, - 0.97616000,0.8592200,0.8414600,0.771788,0.658673, - 0.63818100,0.5706950,0.4619520,0.442280,0.378345, - 0.27670600,0.2562140,0.1865420,0.074247,0.056214}; - G4int WLS_emi_entries = 30; - mpt->AddProperty("WLSCOMPONENT", WLS_emi_energy, WLS_emi_Spectrum, WLS_emi_entries); - - //time that the WLS takes to emmit the absorved photon - mpt->AddConstProperty("WLSTIMECONSTANT", 1. * ns); - //mpt->AddConstProperty("WLSMEANNUMBERPHOTONS", 1); - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::PTP() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - - // Literature is not very extensive about this material - // Refractive index (RINDEX) - // https://en.wikipedia.org/wiki/Terphenyl it is not the best reference ever, i know. - G4double rindex_energy[] = {OpticalMaterialProperties::energy_min, - OpticalMaterialProperties::energy_max}; - //TODO: look for RINDEX dependence - G4double rindex[] = {1.65,1.65}; - int rindex_entries = 2; - - mpt->AddProperty("RINDEX", rindex_energy, rindex, rindex_entries); - - // Absorption length (WLSABSLENGTH) - // for the moment it is assumed that all LAr photons are absorved and then all emited photons pass through - //TODO: look for real spectrum - G4double abslength_energy[] = {OpticalMaterialProperties::energy_min, 5.49*eV, - 5.5*eV, OpticalMaterialProperties::energy_max}; - G4double abslength[] = {OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_min, - OpticalMaterialProperties::abslength_min}; - int abslength_entries = 4; - mpt->AddProperty("WLSABSLENGTH", abslength_energy, abslength, abslength_entries); - - //PTP emision spsectra "WLSCOMPONENT" - //https://deepblue.lib.umich.edu/bitstream/handle/2027.42/30880/0000545.pdf;jsessionid=574463A6129BB3C95B63594EBC262D68?sequence=1 - G4double emi_energy[] = {3.06894*eV,3.15885*eV,3.21398*eV,3.30108*eV,3.32320*eV, - 3.36117*eV,3.41334*eV,3.43137*eV,3.45221*eV,3.46007*eV, - 3.47011*eV,3.47800*eV,3.50550*eV,3.51258*eV,3.52198*eV, - 3.52792*eV,3.53843*eV,3.55350*eV,3.56298*eV,3.57339*eV, - 3.59301*eV,3.60393*eV,3.62111*eV,3.63872*eV,3.64298*eV, - 3.65721*eV,3.66756*eV,3.67709*eV,3.68821*eV,3.72512*eV}; - G4double emi_Spectrum[] = {0.0226077,0.101296,0.199985,0.378345,0.506214, - 0.6566240,0.712772,0.752936,0.864903,0.931214, - 0.9730180,0.989821,0.966870,0.922608,0.872608, - 0.8275260,0.782444,0.738181,0.742280,0.785996, - 0.8144110,0.763045,0.681078,0.517143,0.432444, - 0.3308040,0.245012,0.170695,0.082170,0.006214}; - G4int emi_entries = 30; - - mpt->AddProperty("WLSCOMPONENT", emi_energy, emi_Spectrum, emi_entries); - - mpt->AddConstProperty("WLSTIMECONSTANT", 1. * ns); - //mpt->AddConstProperty("WLSMEANNUMBERPHOTONS", 1); - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::VIKUITI() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - //info from David Warner 11/2019 talk - //it is said that the reflective foils are coated with TPB, - //shall we care about this? - - //reflection (REFLECTIVITY) - //from this talk https://indico.fnal.gov/event/24273/contributions/188657/attachments/130083/158244/DUNE_60Review1.pdf - //we can see depedences in angle for large wavelengths > 800 nm. - //For shorter wavelength almost no angle dependence, so we are using 60º curves. - //We are only considering two options: the photon is either reflected or bulk-absorved. - //If we want to consider transmission to the outer world we need to know refractive index. - const G4int entries = 50; - G4double energy[entries] = {h_Planck * c_light / (351.408*nm),h_Planck * c_light / (360.963*nm), - h_Planck * c_light / (368.310*nm),h_Planck * c_light / (372.535*nm), - h_Planck * c_light / (375.704*nm),h_Planck * c_light / (376.761*nm), - h_Planck * c_light / (378.873*nm),h_Planck * c_light / (377.817*nm), - h_Planck * c_light / (379.930*nm),h_Planck * c_light / (382.042*nm), - h_Planck * c_light / (383.099*nm),h_Planck * c_light / (385.211*nm), - h_Planck * c_light / (387.621*nm),h_Planck * c_light / (393.662*nm), - h_Planck * c_light / (401.987*nm),h_Planck * c_light / (424.296*nm), - h_Planck * c_light / (451.312*nm),h_Planck * c_light / (476.056*nm), - h_Planck * c_light / (496.127*nm),h_Planck * c_light / (528.873*nm), - h_Planck * c_light / (555.282*nm),h_Planck * c_light / (591.197*nm), - h_Planck * c_light / (632.394*nm),h_Planck * c_light / (665.141*nm), - h_Planck * c_light / (691.231*nm),h_Planck * c_light / (720.826*nm), - h_Planck * c_light / (726.189*nm),h_Planck * c_light / (759.155*nm), - h_Planck * c_light / (796.127*nm),h_Planck * c_light / (837.324*nm), - h_Planck * c_light / (883.803*nm),h_Planck * c_light / (901.460*nm), - h_Planck * c_light / (911.037*nm),h_Planck * c_light / (923.967*nm), - h_Planck * c_light / (925.883*nm),h_Planck * c_light / (933.545*nm), - h_Planck * c_light / (934.981*nm),h_Planck * c_light / (940.845*nm), - h_Planck * c_light / (946.953*nm),h_Planck * c_light / (948.869*nm), - h_Planck * c_light / (957.010*nm),h_Planck * c_light / (958.925*nm), - h_Planck * c_light / (968.503*nm),h_Planck * c_light / (1008.45*nm), - h_Planck * c_light / (1038.03*nm),h_Planck * c_light / (1060.21*nm), - h_Planck * c_light / (1048.59*nm),h_Planck * c_light / (1103.52*nm), - h_Planck * c_light / (1154.23*nm),h_Planck * c_light / (1186.97*nm)}; - - G4double reflectivity[entries] = {0.120930,0.132038,0.125581,0.153488,0.206977, - 0.272093,0.479070,0.376744,0.548837,0.623256, - 0.702326,0.781395,0.848983,0.932558,0.970015, - 0.979070,0.971596,0.983721,0.983721,0.979070, - 0.988372,0.986047,0.981395,0.988372,0.970542, - 0.969909,0.973705,0.990698,0.990698,0.993023, - 0.995349,0.961053,0.912205,0.788151,0.764430, - 0.658897,0.638444,0.553488,0.449727,0.420735, - 0.310914,0.287368,0.216732,0.186047,0.169767, - 0.172093,0.162791,0.167442,0.158140,0.174419}; - - mpt->AddProperty("REFLECTIVITY" , energy, reflectivity , entries); - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::FusedSilica() -{ - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - //refractive index RINDEX - //data from https://www.filmetrics.com/refractive-index-database/SiO2 - ///Fused-Silica-Silica-Silicon-Dioxide-Thermal-Oxide-ThermalOxide - //it can be seen from https://nvlpubs.nist.gov/nistpubs/jres/75A/jresv75An4p279_A1b.pdf - //that there is almost no dependence with temperature. - //since PTP emits at 3.45 eV (~360 nm) and WLS emits at 2.9 eV (~430 nm) - //no need to go to lower wavelengths - G4double energy[] = {h_Planck * c_light / (210*nm) ,h_Planck * c_light / (220*nm), - h_Planck * c_light / (240*nm) ,h_Planck * c_light / (260*nm), - h_Planck * c_light / (280*nm) ,h_Planck * c_light / (300*nm), - h_Planck * c_light / (320*nm) ,h_Planck * c_light / (340*nm), - h_Planck * c_light / (360*nm) ,h_Planck * c_light / (380*nm), - h_Planck * c_light / (400*nm) ,h_Planck * c_light / (420*nm), - h_Planck * c_light / (440*nm) ,h_Planck * c_light / (460*nm), - h_Planck * c_light / (480*nm) ,h_Planck * c_light / (500*nm), - h_Planck * c_light / (520*nm) ,h_Planck * c_light / (540*nm), - h_Planck * c_light / (560*nm) ,h_Planck * c_light / (580*nm), - h_Planck * c_light / (600*nm) ,h_Planck * c_light / (650*nm), - h_Planck * c_light / (700*nm) ,h_Planck * c_light / (750*nm), - h_Planck * c_light / (800*nm) ,h_Planck * c_light / (850*nm), - h_Planck * c_light / (900*nm) ,h_Planck * c_light / (1000*nm), - h_Planck * c_light / (1100*nm),h_Planck * c_light / (1200*nm)}; - G4double rindex[] = {1.5384,1.5285,1.5133,1.5024,1.4942, - 1.4878,1.4827,1.4787,1.4753,1.4725, - 1.4701,1.4681,1.4663,1.4648,1.4635, - 1.4623,1.4613,1.4603,1.4595,1.4587, - 1.4580,1.4565,1.4553,1.4542,1.4533, - 1.4525,1.4518,1.4504,1.4492,1.4481}; - G4int entries = 30; - - mpt->AddProperty("RINDEX", energy, rindex, entries); - - //data for the dichroic filter is stored in /data/dichroic_data - - - return mpt; -} - -G4MaterialPropertiesTable* OpticalMaterialProperties::GlassEpoxy() -{ - // Optical properties of Optorez 1330 glass epoxy. - // Obtained from http://refractiveindex.info and - // https://www.zeonex.com/Optics.aspx.html#glass-like - G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); - // REFRACTIVE INDEX - // The range is chosen to be up to ~10.7 eV because Sellmeier's equation - // for fused silica is valid only in that range - const G4int ri_entries = 200; - G4double eWidth = - (OpticalMaterialProperties::energy_max - OpticalMaterialProperties::energy_min) / ri_entries; - G4double ri_energy[ri_entries]; - for (int i = 0; i < ri_entries; i++){ - ri_energy[i] = OpticalMaterialProperties::energy_min + i * eWidth; - } - G4double rIndex[ri_entries]; - for (int i = 0; i < ri_entries; i++) { - G4double lambda = h_Planck*c_light/ri_energy[i]*1000; // in micron - G4double n2 = 2.291142 - 3.311944E-2*pow(lambda,2) - 1.630099E-2*pow(lambda,-2) + - 7.265983E-3*pow(lambda,-4) - 6.806145E-4*pow(lambda,-6) + - 1.960732E-5*pow(lambda,-8); - rIndex[i] = sqrt(n2); - // G4cout << "* GlassEpoxy rIndex: " << std::setw(5) - // << ri_energy[i]/eV << " eV -> " << rIndex[i] << G4endl; - } - assert(sizeof(rIndex) == sizeof(ri_energy)); - mpt->AddProperty("RINDEX", ri_energy, rIndex, ri_entries); - // ABSORPTION LENGTH - G4double abs_energy[] = { - OpticalMaterialProperties::energy_max, - 2.001 * eV, 2.132 * eV, 2.735 * eV, 2.908 * eV, - 3.119 * eV, 3.320 * eV, 3.476 * eV, 3.588 * eV, - 3.749 * eV, 3.869 * eV, 3.973 * eV, 4.120 * eV, - 4.224 * eV, 4.320 * eV, 4.420 * eV, 5.018 * eV - }; - const G4int abs_entries = sizeof(abs_energy) / sizeof(G4double); - G4double absLength[] = { - OpticalMaterialProperties::abslength_max, - OpticalMaterialProperties::abslength_max, - 326.00 * mm, 117.68 * mm, 85.89 * mm, - 50.93 * mm, 31.25 * mm , 17.19 * mm, 10.46 * mm, - 5.26 * mm , 3.77 * mm , 2.69 * mm, 1.94 * mm, - 1.33 * mm , 0.73 * mm , 0.32 * mm, 0.10 * mm - }; - assert(sizeof(absLength) == sizeof(abs_energy)); - mpt->AddProperty("ABSLENGTH", abs_energy, absLength, abs_entries); - return mpt; -} +// ----------------------------------------------------------------------------- +// G4OpSim | OpticalMaterialProperties.cpp +// +// * Author: +// * Creation date: 26 March 2020 +// ----------------------------------------------------------------------------- + +#include "OpticalMaterialProperties.h" + +#include + +#include + +#include + +using CLHEP::h_Planck; +using CLHEP::c_light; + + +G4MaterialPropertiesTable* OpticalMaterialProperties::Vacuum() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + + G4double energies[] = {OpticalMaterialProperties::energy_min, + OpticalMaterialProperties::energy_max}; + + // Refractive index (RINDEX) + G4double rindex[] = {1., 1.}; + mpt->AddProperty("RINDEX", energies, rindex, 2); + + // Absorption length (ABSLENGTH) + G4double abslength[] = {OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max}; + mpt->AddProperty("ABSLENGTH", energies, abslength, 2); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::LAr() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + + G4double energy_step = 0.05 * eV; + G4double energy_curr = OpticalMaterialProperties::energy_min; + std::vector energies; + + do { + energies.push_back(energy_curr); + energy_curr += energy_step; + } while(energy_curr < OpticalMaterialProperties::energy_max); + + // Refractive index (RINDEX) + // https://arxiv.org/pdf/2002.09346.pdf + G4double a0 = 0.335; + G4double aUV = 0.099; + G4double aIR = 0.008; + G4double lambdaUV = 106.6 * nm; + G4double lambdaIR = 908.3 * nm; + std::vector rindex; + + for (auto energy: energies) { + + G4double wavelength = h_Planck * c_light / energy; + G4double x = a0 + aUV * pow (wavelength, 2) / (pow (wavelength, 2) - pow (lambdaUV, 2)) + aIR * pow (wavelength, 2) / (pow (wavelength, 2) - pow (lambdaIR, 2)); + G4double rindex_element = sqrt (1 + 3 * x / (3 - x)); + rindex.push_back(rindex_element); + } + + mpt->AddProperty("RINDEX", energies.data(), rindex.data(), energies.size()); + + // Absorption length (ABSLENGTH) + G4double energies_lim[] = {OpticalMaterialProperties::energy_min, + OpticalMaterialProperties::energy_max}; + G4double abslength[] = {OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max}; + mpt->AddProperty("ABSLENGTH", energies_lim, abslength, 2); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::PVT() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + + G4double energy_step = 0.05 * eV; + G4double energy_curr = OpticalMaterialProperties::energy_min; + std::vector energies; + + do { + energies.push_back(energy_curr); + energy_curr += energy_step; + } while(energy_curr < 7.5 * eV); // Put 7.5 eV instead of energy_max = 11.3 eV in order to avoid divergence. + + // Refractive index (RINDEX) + // The refractive index is fitted to a Sellmeier function: https://en.wikipedia.org/wiki/Sellmeier_equation + G4double A = 1.421; + G4double B1 = 0.9944; + G4double C1 = 26250 * pow (nm, 2); + std::vector rindex; + + for (auto energy: energies) { + + G4double wavelength = h_Planck * c_light / energy; + G4double rindex_element = sqrt (A + B1 * pow (wavelength, 2) / (pow (wavelength, 2) - C1));; + rindex.push_back(rindex_element); + } + + mpt->AddProperty("RINDEX", energies.data(), rindex.data(), energies.size()); + + //Absorption length (ABSLENGTH) + //assuming just WLS can absorve + G4double energies_lim[] = {OpticalMaterialProperties::energy_min, + OpticalMaterialProperties::energy_max}; + G4double abslength[] = {OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max}; + mpt->AddProperty("ABSLENGTH", energies_lim, abslength, 2); + + //WLS absorption WLSABSLENGTH + //assuming maximum absorption length for extremal value + G4double WLS_abslength[] = {OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max, + 879.673*mm ,406.487*mm ,296.786*mm,165.772*mm, + 118.88*mm ,96.5343*mm ,70.4532*mm,38.5933*mm, + 26.714 *mm ,16.9549*mm ,11.2558*mm,9.49174*mm, + 8.78259*mm ,7.31249*mm ,6.59717*mm,6.05747*mm, + 5.24668*mm ,4.38833*mm ,3.73978*mm,2.71100*mm, + 2.1048*mm ,1.97013*mm ,1.20762*mm,0.71176*mm, + 0.803417*mm,0.803417*mm,0.97114*mm,1.63289*mm}; + G4double WLS_abs_energy[] = {OpticalMaterialProperties::energy_min, + 2.90388, + 2.9723 *eV,2.98595*eV,2.9987 *eV,3.00867*eV, + 3.01317*eV,3.01997*eV,3.02443*eV,3.03555*eV, + 3.04441*eV,3.0508*eV ,3.05709*eV,3.06822*eV, + 3.07265*eV,3.07916*eV,3.08353*eV,3.08565*eV, + 3.09449*eV,3.10326*eV,3.11208*eV,3.12862*eV, + 3.14538*eV,3.15101*eV,3.20314*eV,3.25895*eV, + 3.27264*eV,3.88859*eV,3.91598*eV,4.01202*eV}; + G4int WLS_abs_entries = 30; + mpt->AddProperty("WLSABSLENGTH", WLS_abs_energy, WLS_abslength, WLS_abs_entries); + + // Emision spectrum (WLSCOMPONENT) + // from https://arxiv.org/pdf/1912.09191.pdf EJ286 + G4double WLS_emi_energy[] = {2.32899*eV,2.38369*eV,2.46565*eV,2.57017*eV,2.61099*eV, + 2.64714*eV,2.69664*eV,2.72479*eV,2.77871*eV,2.79583*eV, + 2.84050*eV,2.85033*eV,2.86809*eV,2.89408*eV,2.91221*eV, + 2.94121*eV,2.97986*eV,2.98148*eV,3.00011*eV,3.01669*eV, + 3.01794*eV,3.03617*eV,3.04968*eV,3.04979*eV,3.06040*eV, + 3.07350*eV,3.07475*eV,3.08767*eV,3.11679*eV,3.12031*eV}; + G4double WLS_emi_Spectrum[] = {0.00293557,0.0122252,0.0439192,0.121378,0.166870, + 0.23135100,0.3592200,0.4271160,0.539329,0.579575, + 0.75539500,0.7947390,0.8865420,0.973837,0.995832, + 0.97616000,0.8592200,0.8414600,0.771788,0.658673, + 0.63818100,0.5706950,0.4619520,0.442280,0.378345, + 0.27670600,0.2562140,0.1865420,0.074247,0.056214}; + G4int WLS_emi_entries = 30; + mpt->AddProperty("WLSCOMPONENT", WLS_emi_energy, WLS_emi_Spectrum, WLS_emi_entries); + + //time that the WLS takes to emmit the absorved photon + mpt->AddConstProperty("WLSTIMECONSTANT", 1. * ns); + //mpt->AddConstProperty("WLSMEANNUMBERPHOTONS", 1); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::PTP() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + + // Literature is not very extensive about this material + // Refractive index (RINDEX) + // https://en.wikipedia.org/wiki/Terphenyl it is not the best reference ever, i know. + G4double rindex_energy[] = {OpticalMaterialProperties::energy_min, + OpticalMaterialProperties::energy_max}; + //TODO: look for RINDEX dependence + G4double rindex[] = {1.65,1.65}; + int rindex_entries = 2; + + mpt->AddProperty("RINDEX", rindex_energy, rindex, rindex_entries); + + // Absorption length (WLSABSLENGTH) + // for the moment it is assumed that all LAr photons are absorved and then all emited photons pass through + //TODO: look for real spectrum + G4double abslength_energy[] = {OpticalMaterialProperties::energy_min, 5.49*eV, + 5.5*eV, OpticalMaterialProperties::energy_max}; + G4double abslength[] = {OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_min, + OpticalMaterialProperties::abslength_min}; + int abslength_entries = 4; + mpt->AddProperty("WLSABSLENGTH", abslength_energy, abslength, abslength_entries); + + //PTP emision spsectra "WLSCOMPONENT" + //https://deepblue.lib.umich.edu/bitstream/handle/2027.42/30880/0000545.pdf;jsessionid=574463A6129BB3C95B63594EBC262D68?sequence=1 + G4double emi_energy[] = {3.06894*eV,3.15885*eV,3.21398*eV,3.30108*eV,3.32320*eV, + 3.36117*eV,3.41334*eV,3.43137*eV,3.45221*eV,3.46007*eV, + 3.47011*eV,3.47800*eV,3.50550*eV,3.51258*eV,3.52198*eV, + 3.52792*eV,3.53843*eV,3.55350*eV,3.56298*eV,3.57339*eV, + 3.59301*eV,3.60393*eV,3.62111*eV,3.63872*eV,3.64298*eV, + 3.65721*eV,3.66756*eV,3.67709*eV,3.68821*eV,3.72512*eV}; + G4double emi_Spectrum[] = {0.0226077,0.101296,0.199985,0.378345,0.506214, + 0.6566240,0.712772,0.752936,0.864903,0.931214, + 0.9730180,0.989821,0.966870,0.922608,0.872608, + 0.8275260,0.782444,0.738181,0.742280,0.785996, + 0.8144110,0.763045,0.681078,0.517143,0.432444, + 0.3308040,0.245012,0.170695,0.082170,0.006214}; + G4int emi_entries = 30; + + mpt->AddProperty("WLSCOMPONENT", emi_energy, emi_Spectrum, emi_entries); + + mpt->AddConstProperty("WLSTIMECONSTANT", 1. * ns); + //mpt->AddConstProperty("WLSMEANNUMBERPHOTONS", 1); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::VIKUITI() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + //info from David Warner 11/2019 talk + //it is said that the reflective foils are coated with TPB, + //shall we care about this? + + //reflection (REFLECTIVITY) + //from this talk https://indico.fnal.gov/event/24273/contributions/188657/attachments/130083/158244/DUNE_60Review1.pdf + //we can see depedences in angle for large wavelengths > 800 nm. + //For shorter wavelength almost no angle dependence, so we are using 60º curves. + //We are only considering two options: the photon is either reflected or bulk-absorved. + //If we want to consider transmission to the outer world we need to know refractive index. + const G4int entries = 50; + G4double energy[entries] = {h_Planck * c_light / (351.408*nm),h_Planck * c_light / (360.963*nm), + h_Planck * c_light / (368.310*nm),h_Planck * c_light / (372.535*nm), + h_Planck * c_light / (375.704*nm),h_Planck * c_light / (376.761*nm), + h_Planck * c_light / (378.873*nm),h_Planck * c_light / (377.817*nm), + h_Planck * c_light / (379.930*nm),h_Planck * c_light / (382.042*nm), + h_Planck * c_light / (383.099*nm),h_Planck * c_light / (385.211*nm), + h_Planck * c_light / (387.621*nm),h_Planck * c_light / (393.662*nm), + h_Planck * c_light / (401.987*nm),h_Planck * c_light / (424.296*nm), + h_Planck * c_light / (451.312*nm),h_Planck * c_light / (476.056*nm), + h_Planck * c_light / (496.127*nm),h_Planck * c_light / (528.873*nm), + h_Planck * c_light / (555.282*nm),h_Planck * c_light / (591.197*nm), + h_Planck * c_light / (632.394*nm),h_Planck * c_light / (665.141*nm), + h_Planck * c_light / (691.231*nm),h_Planck * c_light / (720.826*nm), + h_Planck * c_light / (726.189*nm),h_Planck * c_light / (759.155*nm), + h_Planck * c_light / (796.127*nm),h_Planck * c_light / (837.324*nm), + h_Planck * c_light / (883.803*nm),h_Planck * c_light / (901.460*nm), + h_Planck * c_light / (911.037*nm),h_Planck * c_light / (923.967*nm), + h_Planck * c_light / (925.883*nm),h_Planck * c_light / (933.545*nm), + h_Planck * c_light / (934.981*nm),h_Planck * c_light / (940.845*nm), + h_Planck * c_light / (946.953*nm),h_Planck * c_light / (948.869*nm), + h_Planck * c_light / (957.010*nm),h_Planck * c_light / (958.925*nm), + h_Planck * c_light / (968.503*nm),h_Planck * c_light / (1008.45*nm), + h_Planck * c_light / (1038.03*nm),h_Planck * c_light / (1060.21*nm), + h_Planck * c_light / (1048.59*nm),h_Planck * c_light / (1103.52*nm), + h_Planck * c_light / (1154.23*nm),h_Planck * c_light / (1186.97*nm)}; + + G4double reflectivity[entries] = {0.120930,0.132038,0.125581,0.153488,0.206977, + 0.272093,0.479070,0.376744,0.548837,0.623256, + 0.702326,0.781395,0.848983,0.932558,0.970015, + 0.979070,0.971596,0.983721,0.983721,0.979070, + 0.988372,0.986047,0.981395,0.988372,0.970542, + 0.969909,0.973705,0.990698,0.990698,0.993023, + 0.995349,0.961053,0.912205,0.788151,0.764430, + 0.658897,0.638444,0.553488,0.449727,0.420735, + 0.310914,0.287368,0.216732,0.186047,0.169767, + 0.172093,0.162791,0.167442,0.158140,0.174419}; + + mpt->AddProperty("REFLECTIVITY" , energy, reflectivity , entries); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::FusedSilica() +{ + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + //refractive index RINDEX + //data from https://www.filmetrics.com/refractive-index-database/SiO2 + ///Fused-Silica-Silica-Silicon-Dioxide-Thermal-Oxide-ThermalOxide + //it can be seen from https://nvlpubs.nist.gov/nistpubs/jres/75A/jresv75An4p279_A1b.pdf + //that there is almost no dependence with temperature. + //since PTP emits at 3.45 eV (~360 nm) and WLS emits at 2.9 eV (~430 nm) + //no need to go to lower wavelengths + G4double energy[] = {h_Planck * c_light / (210*nm) ,h_Planck * c_light / (220*nm), + h_Planck * c_light / (240*nm) ,h_Planck * c_light / (260*nm), + h_Planck * c_light / (280*nm) ,h_Planck * c_light / (300*nm), + h_Planck * c_light / (320*nm) ,h_Planck * c_light / (340*nm), + h_Planck * c_light / (360*nm) ,h_Planck * c_light / (380*nm), + h_Planck * c_light / (400*nm) ,h_Planck * c_light / (420*nm), + h_Planck * c_light / (440*nm) ,h_Planck * c_light / (460*nm), + h_Planck * c_light / (480*nm) ,h_Planck * c_light / (500*nm), + h_Planck * c_light / (520*nm) ,h_Planck * c_light / (540*nm), + h_Planck * c_light / (560*nm) ,h_Planck * c_light / (580*nm), + h_Planck * c_light / (600*nm) ,h_Planck * c_light / (650*nm), + h_Planck * c_light / (700*nm) ,h_Planck * c_light / (750*nm), + h_Planck * c_light / (800*nm) ,h_Planck * c_light / (850*nm), + h_Planck * c_light / (900*nm) ,h_Planck * c_light / (1000*nm), + h_Planck * c_light / (1100*nm),h_Planck * c_light / (1200*nm)}; + /*G4double rindex[] = {1.5384,1.5285,1.5133,1.5024,1.4942, + 1.4878,1.4827,1.4787,1.4753,1.4725, + 1.4701,1.4681,1.4663,1.4648,1.4635, + 1.4623,1.4613,1.4603,1.4595,1.4587, + 1.4580,1.4565,1.4553,1.4542,1.4533, + 1.4525,1.4518,1.4504,1.4492,1.4481};*/ + G4double rindex[] = {1.,1.,1.,1.,1., + 1.,1.,1.,1.,1., + 1.,1.,1.,1.,1., + 1.,1.,1.,1.,1., + 1.,1.,1.,1.,1., + 1.,1.,1.,1.,1.}; + + G4int entries = 30; + + mpt->AddProperty("RINDEX", energy, rindex, entries); + + return mpt; +} + +G4MaterialPropertiesTable* OpticalMaterialProperties::GlassEpoxy() +{ + // Optical properties of Optorez 1330 glass epoxy. + // Obtained from http://refractiveindex.info and + // https://www.zeonex.com/Optics.aspx.html#glass-like + G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); + // REFRACTIVE INDEX + // The range is chosen to be up to ~10.7 eV because Sellmeier's equation + // for fused silica is valid only in that range + const G4int ri_entries = 200; + G4double eWidth = + (OpticalMaterialProperties::energy_max - OpticalMaterialProperties::energy_min) / ri_entries; + G4double ri_energy[ri_entries]; + for (int i = 0; i < ri_entries; i++){ + ri_energy[i] = OpticalMaterialProperties::energy_min + i * eWidth; + } + G4double rIndex[ri_entries]; + for (int i = 0; i < ri_entries; i++) { + G4double lambda = h_Planck*c_light/ri_energy[i]*1000; // in micron + G4double n2 = 2.291142 - 3.311944E-2*pow(lambda,2) - 1.630099E-2*pow(lambda,-2) + + 7.265983E-3*pow(lambda,-4) - 6.806145E-4*pow(lambda,-6) + + 1.960732E-5*pow(lambda,-8); + rIndex[i] = sqrt(n2); + // G4cout << "* GlassEpoxy rIndex: " << std::setw(5) + // << ri_energy[i]/eV << " eV -> " << rIndex[i] << G4endl; + } + assert(sizeof(rIndex) == sizeof(ri_energy)); + mpt->AddProperty("RINDEX", ri_energy, rIndex, ri_entries); + // ABSORPTION LENGTH + G4double abs_energy[] = { + OpticalMaterialProperties::energy_max, + 2.001 * eV, 2.132 * eV, 2.735 * eV, 2.908 * eV, + 3.119 * eV, 3.320 * eV, 3.476 * eV, 3.588 * eV, + 3.749 * eV, 3.869 * eV, 3.973 * eV, 4.120 * eV, + 4.224 * eV, 4.320 * eV, 4.420 * eV, 5.018 * eV + }; + const G4int abs_entries = sizeof(abs_energy) / sizeof(G4double); + G4double absLength[] = { + OpticalMaterialProperties::abslength_max, + OpticalMaterialProperties::abslength_max, + 326.00 * mm, 117.68 * mm, 85.89 * mm, + 50.93 * mm, 31.25 * mm , 17.19 * mm, 10.46 * mm, + 5.26 * mm , 3.77 * mm , 2.69 * mm, 1.94 * mm, + 1.33 * mm , 0.73 * mm , 0.32 * mm, 0.10 * mm + }; + assert(sizeof(absLength) == sizeof(abs_energy)); + mpt->AddProperty("ABSLENGTH", abs_energy, absLength, abs_entries); + return mpt; +} diff --git a/src/OpticalMaterialProperties.h b/src/OpticalMaterialProperties.h index 4a2a8c0..8f60ef8 100644 --- a/src/OpticalMaterialProperties.h +++ b/src/OpticalMaterialProperties.h @@ -1,41 +1,41 @@ -// ----------------------------------------------------------------------------- -// G4OpSim | OpticalMaterialProperties.h -// -// * Author: -// * Creation date: 26 March 2020 -// ----------------------------------------------------------------------------- - -#ifndef OPTICAL_MATERIAL_PROPERTIES_H -#define OPTICAL_MATERIAL_PROPERTIES_H - -#include -#include - -class G4MaterialPropertiesTable; - - -namespace OpticalMaterialProperties { - - inline const G4double energy_min = 2.0 * eV; // 620 nm - inline const G4double energy_max = 11.3 * eV; // 110 nm - - inline const G4double abslength_min = 1. * nm; - inline const G4double abslength_max = 1.E4 * m; - - G4MaterialPropertiesTable* Vacuum(); - - G4MaterialPropertiesTable* LAr(); - - G4MaterialPropertiesTable* PVT(); - - G4MaterialPropertiesTable* PTP(); - - G4MaterialPropertiesTable* VIKUITI(); - - G4MaterialPropertiesTable* FusedSilica(); - - G4MaterialPropertiesTable* GlassEpoxy(); - -} // end namespace - -#endif +// ----------------------------------------------------------------------------- +// G4OpSim | OpticalMaterialProperties.h +// +// * Author: +// * Creation date: 26 March 2020 +// ----------------------------------------------------------------------------- + +#ifndef OPTICAL_MATERIAL_PROPERTIES_H +#define OPTICAL_MATERIAL_PROPERTIES_H + +#include +#include + +class G4MaterialPropertiesTable; + + +namespace OpticalMaterialProperties { + + inline const G4double energy_min = 2.0 * eV; // 620 nm + inline const G4double energy_max = 11.3 * eV; // 110 nm + + inline const G4double abslength_min = 1. * nm; + inline const G4double abslength_max = 1.E4 * m; + + G4MaterialPropertiesTable* Vacuum(); + + G4MaterialPropertiesTable* LAr(); + + G4MaterialPropertiesTable* PVT(); + + G4MaterialPropertiesTable* PTP(); + + G4MaterialPropertiesTable* VIKUITI(); + + G4MaterialPropertiesTable* FusedSilica(); + + G4MaterialPropertiesTable* GlassEpoxy(); + +} // end namespace + +#endif diff --git a/src/OpticalSD.cpp b/src/OpticalSD.cpp index 2e5054f..064be80 100644 --- a/src/OpticalSD.cpp +++ b/src/OpticalSD.cpp @@ -1,40 +1,40 @@ -// ----------------------------------------------------------------------------- -// G4Basic | OpticalSD.cpp -// -// Sensitive detector for optical sensors (e.g. SiPMs). -// * Author: -// * Creation date: 5 Nov 2020 -// ----------------------------------------------------------------------------- - -#include "OpticalSD.h" -#include "OpticalHit.h" - -#include - - - -OpticalSD::OpticalSD(const G4String& sdname): - G4VSensitiveDetector(sdname) -{ - collectionName.insert("Optical"); -} - - -OpticalSD::~OpticalSD() -{ -} - - -void OpticalSD::Initialize(G4HCofThisEvent* hce) -{ -} - - -G4bool OpticalSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) -{ -} - - -void OpticalSD::EndOfEvent(G4HCofThisEvent*) -{ -} +// ----------------------------------------------------------------------------- +// G4Basic | OpticalSD.cpp +// +// Sensitive detector for optical sensors (e.g. SiPMs). +// * Author: +// * Creation date: 5 Nov 2020 +// ----------------------------------------------------------------------------- + +#include "OpticalSD.h" +#include "OpticalHit.h" + +#include + + + +OpticalSD::OpticalSD(const G4String& sdname): + G4VSensitiveDetector(sdname) +{ + collectionName.insert("Optical"); +} + + +OpticalSD::~OpticalSD() +{ +} + + +void OpticalSD::Initialize(G4HCofThisEvent* hce) +{ +} + + +G4bool OpticalSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ +} + + +void OpticalSD::EndOfEvent(G4HCofThisEvent*) +{ +} diff --git a/src/OpticalSD.h b/src/OpticalSD.h index 3804ff4..804fc0d 100644 --- a/src/OpticalSD.h +++ b/src/OpticalSD.h @@ -1,30 +1,30 @@ -// ----------------------------------------------------------------------------- -// G4Basic | OpticalSD.h -// -// Sensitive detector for optical sensors (e.g. SiPMs). -// * Author: -// * Creation date: 5 Nov 2020 -// ----------------------------------------------------------------------------- - -#ifndef OPTICAL_SD_H -#define OPTICAL_SD_H - -#include -#include "OpticalHit.h" - - -class OpticalSD: public G4VSensitiveDetector -{ -public: - OpticalSD(const G4String&); - ~OpticalSD(); - - void Initialize(G4HCofThisEvent*) override; - G4bool ProcessHits(G4Step*, G4TouchableHistory*) override; - void EndOfEvent(G4HCofThisEvent*) override; - -private: - OpticalHitCollection* hc_; -}; - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | OpticalSD.h +// +// Sensitive detector for optical sensors (e.g. SiPMs). +// * Author: +// * Creation date: 5 Nov 2020 +// ----------------------------------------------------------------------------- + +#ifndef OPTICAL_SD_H +#define OPTICAL_SD_H + +#include +#include "OpticalHit.h" + + +class OpticalSD: public G4VSensitiveDetector +{ +public: + OpticalSD(const G4String&); + ~OpticalSD(); + + void Initialize(G4HCofThisEvent*) override; + G4bool ProcessHits(G4Step*, G4TouchableHistory*) override; + void EndOfEvent(G4HCofThisEvent*) override; + +private: + OpticalHitCollection* hc_; +}; + +#endif diff --git a/src/PhysicsList.cpp b/src/PhysicsList.cpp index 5b021a4..b77dc2d 100644 --- a/src/PhysicsList.cpp +++ b/src/PhysicsList.cpp @@ -1,34 +1,34 @@ -// ----------------------------------------------------------------------------- -// G4Basic | G4Basic.cpp -// -// -// ----------------------------------------------------------------------------- - -#include "PhysicsList.h" - -#include -#include -#include -#include -#include - - -PhysicsList::PhysicsList(): G4VModularPhysicsList() -{ - RegisterPhysics(new G4EmStandardPhysics_option4()); - RegisterPhysics(new G4DecayPhysics()); - RegisterPhysics(new G4RadioactiveDecayPhysics()); - RegisterPhysics(new G4StepLimiterPhysics()); - RegisterPhysics(new G4OpticalPhysics()); -} - - -PhysicsList::~PhysicsList() -{ -} - - -void PhysicsList::SetCuts() -{ - G4VUserPhysicsList::SetCuts(); -} +// ----------------------------------------------------------------------------- +// G4Basic | G4Basic.cpp +// +// +// ----------------------------------------------------------------------------- + +#include "PhysicsList.h" + +#include +#include +#include +#include +#include + + +PhysicsList::PhysicsList(): G4VModularPhysicsList() +{ + RegisterPhysics(new G4EmStandardPhysics_option4()); + RegisterPhysics(new G4DecayPhysics()); + RegisterPhysics(new G4RadioactiveDecayPhysics()); + RegisterPhysics(new G4StepLimiterPhysics()); + RegisterPhysics(new G4OpticalPhysics()); +} + + +PhysicsList::~PhysicsList() +{ +} + + +void PhysicsList::SetCuts() +{ + G4VUserPhysicsList::SetCuts(); +} diff --git a/src/PhysicsList.h b/src/PhysicsList.h index 54d6625..02c1e88 100644 --- a/src/PhysicsList.h +++ b/src/PhysicsList.h @@ -1,21 +1,21 @@ -// ----------------------------------------------------------------------------- -// G4Basic | PhysicsList.h -// -// -// ----------------------------------------------------------------------------- - -#ifndef PHYSICS_LIST_H -#define PHYSICS_LIST_H - -#include - - -class PhysicsList: public G4VModularPhysicsList -{ -public: - PhysicsList(); - virtual ~PhysicsList(); - virtual void SetCuts(); -}; - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | PhysicsList.h +// +// +// ----------------------------------------------------------------------------- + +#ifndef PHYSICS_LIST_H +#define PHYSICS_LIST_H + +#include + + +class PhysicsList: public G4VModularPhysicsList +{ +public: + PhysicsList(); + virtual ~PhysicsList(); + virtual void SetCuts(); +}; + +#endif diff --git a/src/PrimaryGeneration.cpp b/src/PrimaryGeneration.cpp index df28698..f8339fb 100644 --- a/src/PrimaryGeneration.cpp +++ b/src/PrimaryGeneration.cpp @@ -1,79 +1,104 @@ -// ----------------------------------------------------------------------------- -// G4Basic | PrimaryGeneration.cpp -// -// -// ----------------------------------------------------------------------------- - -#include "PrimaryGeneration.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -PrimaryGeneration::PrimaryGeneration(): - G4VUserPrimaryGeneratorAction(), - kinetic_energy_(6*eV) -{ -} - - -PrimaryGeneration::~PrimaryGeneration() -{ -} - - -void PrimaryGeneration::GeneratePrimaries(G4Event* event) -{ - // Generate a random momentum direction for the photon - - G4double cost = 1. - 2.*G4UniformRand(); - G4double sint = std::sqrt((1.-cost)*(1.+cost)); - - G4double phi = twopi*G4UniformRand(); - G4double sinp = std::sin(phi); - G4double cosp = std::cos(phi); - - G4double px = sint*cosp; - G4double py = sint*sinp; - G4double pz = cost; - - //G4ThreeVector momentum(px,py,pz); - G4ThreeVector momentum(0.,-1,0.); - //G4ThreeVector momentum(0.,-0.99,sqrt(1-0.99*0.99)); - - // Determine the polarization of the photon - - G4double sx = cost*cosp; - G4double sy = cost*sinp; - G4double sz = -sint; - - G4ThreeVector polarization(sx, sy, sz); - - G4ThreeVector perp = momentum.cross(polarization); - - phi = twopi*G4UniformRand(); - sinp = std::sin(phi); - cosp = std::cos(phi); - - polarization = cosp * polarization + sinp * perp; - polarization = polarization.unit(); - - // Create a new photon - - G4PrimaryParticle* particle = new G4PrimaryParticle(G4OpticalPhoton::Definition()); - particle->SetMomentumDirection(momentum); - particle->SetPolarization(polarization); - particle->SetKineticEnergy(kinetic_energy_); - - G4PrimaryVertex* vertex = new G4PrimaryVertex(G4ThreeVector(0.,10.*cm,0.), 0.); - vertex->SetPrimary(particle); - - event->AddPrimaryVertex(vertex); -} +// ----------------------------------------------------------------------------- +// G4Basic | PrimaryGeneration.cpp +// +// +// ----------------------------------------------------------------------------- + +#include +#include +#include + +#include "params.h" +#include "PrimaryGeneration.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using CLHEP::h_Planck; +using CLHEP::c_light; + +PrimaryGeneration::PrimaryGeneration(): + G4VUserPrimaryGeneratorAction(), + kinetic_energy_(6*eV) +{ +} + + +PrimaryGeneration::~PrimaryGeneration() +{ +} + + +void PrimaryGeneration::GeneratePrimaries(G4Event* event) +{ + // Generate a random momentum direction for the photon + + //G4double cost = 1. - 2.*G4UniformRand(); + //G4double sint = std::sqrt((1.-cost)*(1.+cost)); + //G4double theta = twopi / 4. + twopi / 360 * 30 + twopi / 360 * 15 * G4UniformRand(); //isotropic distribution, careful with angles for which the photon misses the ARAPUCA + G4double theta; // = twopi / 4. + twopi / 360. * 45.; + G4double randpicker = G4UniformRand(); + if (randpicker < 0.5) theta = twopi / 4. + twopi / 360 * 30; + else theta = twopi / 4. + twopi / 360 * 45; + G4double cost = std::cos(theta); + G4double sint = std::sin(theta); + + G4double phi = twopi / 4.; //* G4UniformRand(); + G4double sinp = std::sin(phi); //twopi/2. + twopi/360*30+ twopi/360*15*G4UniformRand(); + G4double cosp = std::cos(phi); + + G4double wavelengths[29] = { 250, 295, 300, 305, 310, 315, 320, 325, 330, 335, 340, 345, 350, 355, 360, 365, 370, 375, 380, 385, 390, 395, 400, 405, 410, 415, 420, 425, 500 }; + //srand(time(NULL)); + int randPicker = int((G4UniformRand()*1000000.)) % sizeof(wavelengths); + //G4double wavelengths_sample[3] = {300, 350, 400}; // simulation tests with smaller sample size + //int randPicker_sample = int((G4UniformRand()*100000.))%3; + + G4double px = sint*cosp; + G4double py = sint*sinp; + G4double pz = cost; + + //G4double Edelta = h_Planck * c_light * (1 / (200 * nm) - 1 / (500 * nm) ); // used in combination with a minimal energy in order to simulate a better approximation of a continuous spectrum + + init_Ekin = h_Planck * c_light / (wavelengths[randPicker] * nm); + + G4ThreeVector momentum(px,pz,py); + //G4ThreeVector momentum(0.,-1.,0.); //vertical incidence onto the ARAPUCA + + //momentum.unit() for normalisation + + // Determine the polarization of the photon + + G4double sx = cost*cosp; + G4double sy = cost*sinp; + G4double sz = -sint; + + G4ThreeVector polarization(sx, sy, sz); + + G4ThreeVector perp = momentum.cross(polarization); + + phi = twopi*G4UniformRand(); + sinp = std::sin(phi); + cosp = std::cos(phi); + + polarization = cosp * polarization + sinp * perp; + polarization = polarization.unit(); + + // Create a new photon + + G4PrimaryParticle* particle = new G4PrimaryParticle(G4OpticalPhoton::Definition()); + particle->SetMomentumDirection(momentum); + particle->SetPolarization(polarization); + particle->SetKineticEnergy(init_Ekin); + + G4PrimaryVertex* vertex = new G4PrimaryVertex(G4ThreeVector(0.,10.*cm,0.), 0.); + vertex->SetPrimary(particle); + + event->AddPrimaryVertex(vertex); +} diff --git a/src/PrimaryGeneration.h b/src/PrimaryGeneration.h index 084f03b..9fba3b0 100644 --- a/src/PrimaryGeneration.h +++ b/src/PrimaryGeneration.h @@ -1,27 +1,27 @@ -// ----------------------------------------------------------------------------- -// G4Basic | PrimaryGeneration.cpp -// -// -// ----------------------------------------------------------------------------- - -#ifndef PRIMARY_GENERATION_H -#define PRIMARY_GENERATION_H - -#include -#include - -class G4ParticleDefinition; - - -class PrimaryGeneration: public G4VUserPrimaryGeneratorAction -{ -public: - PrimaryGeneration(); - virtual ~PrimaryGeneration(); - virtual void GeneratePrimaries(G4Event*); - -private: - G4double kinetic_energy_; -}; - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | PrimaryGeneration.cpp +// +// +// ----------------------------------------------------------------------------- + +#ifndef PRIMARY_GENERATION_H +#define PRIMARY_GENERATION_H + +#include +#include + +class G4ParticleDefinition; + + +class PrimaryGeneration: public G4VUserPrimaryGeneratorAction +{ +public: + PrimaryGeneration(); + virtual ~PrimaryGeneration(); + virtual void GeneratePrimaries(G4Event*); + +private: + G4double kinetic_energy_; +}; + +#endif diff --git a/src/RunAction.cpp b/src/RunAction.cpp index 019b47e..ff60773 100644 --- a/src/RunAction.cpp +++ b/src/RunAction.cpp @@ -1,23 +1,103 @@ -// ----------------------------------------------------------------------------- -// G4Basic | RunAction.cpp -// -// User run action class. -// ----------------------------------------------------------------------------- - -#include "RunAction.h" - -#include - - -void RunAction::BeginOfRunAction(const G4Run* run) -{ - G4cout << "------------------------------------------------------------\n" - << "Run ID " << run->GetRunID() << G4endl; -} - -void RunAction::EndOfRunAction(const G4Run*) -{ - G4cout << "End of run." - << "------------------------------------------------------------" - << G4endl; -} +// ----------------------------------------------------------------------------- +// G4Basic | RunAction.cpp +// +// User run action class. +// ----------------------------------------------------------------------------- + +#include "RunAction.h" +#include "Analysis.h" + +#include +#include +#include +#include +#include + + +using CLHEP::h_Planck; +using CLHEP::c_light; +using CLHEP::twopi; + + +RunAction::RunAction() + : G4UserRunAction() +{ + // set printing event number per each event + G4RunManager::GetRunManager()->SetPrintProgress(1); + + // Create analysis manager + // The choice of analysis technology is done via selectin of a namespace + // in B4Analysis.hh + auto analysisManager = G4AnalysisManager::Instance(); + G4cout << "Using " << analysisManager->GetType() << G4endl; + + // Create directories + //analysisManager->SetHistoDirectoryName("histograms"); + //analysisManager->SetNtupleDirectoryName("ntuple"); + analysisManager->SetVerboseLevel(1); + + // Book histograms, ntuple + // + // Creating histograms + analysisManager->CreateH1("trm30", "Transmittance at 30 deg", 100, 150. * nm, 550. * nm); + analysisManager->CreateH1("trm45", "Transmittance at 45 deg", 100, 150. * nm, 550. * nm); + analysisManager->CreateH1("ypos", "y-pos of final particle pos", 100, -1000 * cm, 1000. * cm); + analysisManager->CreateH2("angvswl", "angle vs wavelength", 200, 150. * nm, 550. * nm, 200, 25, 50); + + // Creating ntuple + // + analysisManager->CreateNtuple("trm", "dichroic simulation"); + analysisManager->CreateNtupleDColumn("trm30"); + analysisManager->CreateNtupleDColumn("trm45"); // this and trm30 is for quick check of distributions, below is full information that can be used for post simulation analysis + analysisManager->CreateNtupleDColumn("ypos_end"); + analysisManager->CreateNtupleDColumn("angle"); + analysisManager->CreateNtupleDColumn("wavelength"); + analysisManager->FinishNtuple(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +RunAction::~RunAction() +{ + delete G4AnalysisManager::Instance(); +} + +void RunAction::BeginOfRunAction(const G4Run* run) +{ + G4cout << "------------------------------------------------------------\n" + << "Run ID " << run->GetRunID() << G4endl; + + //inform the runManager to save random number seed + //G4RunManager::GetRunManager()->SetRandomNumberStore(true); + + //TODO: initialise counters here so the numbers are reset for each run + + // Get analysis manager + auto analysisManager = G4AnalysisManager::Instance(); + + // Open an output file + // + G4String fileName = "dichroicTransmittance"; + analysisManager->OpenFile(fileName); + +} + +void RunAction::EndOfRunAction(const G4Run*) +{ + + // print histogram statistics + // + auto analysisManager = G4AnalysisManager::Instance(); + + G4cout << " Trm at 30 deg : mean = " + << analysisManager->GetH1(0)->mean() << G4endl; + + G4cout << "End of run." + << "------------------------------------------------------------" + << G4endl; + + // save histograms & ntuple + // + analysisManager->Write(); + analysisManager->CloseFile(); +} diff --git a/src/RunAction.h b/src/RunAction.h index 2f76116..f4994c6 100644 --- a/src/RunAction.h +++ b/src/RunAction.h @@ -1,27 +1,27 @@ -// ----------------------------------------------------------------------------- -// G4Basic | RunAction.h -// -// User run action class. -// ----------------------------------------------------------------------------- - -#ifndef RUN_ACTION_H -#define RUN_ACTION_H - -#include - -class G4Run; - - -class RunAction: public G4UserRunAction -{ -public: - RunAction(); - virtual ~RunAction(); - virtual void BeginOfRunAction(const G4Run*); - virtual void EndOfRunAction(const G4Run*); -}; - -inline RunAction::RunAction() {} -inline RunAction::~RunAction() {} - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | RunAction.h +// +// User run action class. +// ----------------------------------------------------------------------------- + +#ifndef RUN_ACTION_H +#define RUN_ACTION_H + +#include + +class G4Run; + + +class RunAction: public G4UserRunAction +{ +public: + RunAction(); + virtual ~RunAction(); + virtual void BeginOfRunAction(const G4Run*); + virtual void EndOfRunAction(const G4Run*); +}; + +//inline RunAction::RunAction() {} +//inline RunAction::~RunAction() {} + +#endif diff --git a/src/SteppingAction.cpp b/src/SteppingAction.cpp index b9101b3..31a2844 100644 --- a/src/SteppingAction.cpp +++ b/src/SteppingAction.cpp @@ -1,94 +1,230 @@ -// ----------------------------------------------------------------------------- -// G4Basic | SteppingAction.cpp -// -// User stepping action class. -// ----------------------------------------------------------------------------- - -#include "SteppingAction.h" - -#include -#include -#include -#include -#include -#include - - -SteppingAction::SteppingAction(): - G4UserSteppingAction(), counter(0) -{ - G4cout << "SteppingAction::SteppingAction" << G4endl; -} - - -SteppingAction::~SteppingAction() -{ - G4cout << "SteppingAction::~SteppingAction" << G4endl; - G4cout << "counter: " << counter << G4endl; -} - - -void SteppingAction::UserSteppingAction(const G4Step* step) -{ - if (step->GetTrack()->GetParentID() == 0) return; - - G4ParticleDefinition* pdef = step->GetTrack()->GetDefinition(); - - //Check whether the track is an optical photon - if (pdef != G4OpticalPhoton::Definition()) return; - - auto step_number = step->GetTrack()->GetCurrentStepNumber(); - auto volume_name = step->GetTrack()->GetVolume()->GetName(); - - G4cout << "step_number: " << step_number << ", volume_name: " << volume_name << G4endl; - - if (step_number==3 && volume_name=="WLS_PLATE") ++counter; - - /* - // example of information one can access about optical photons - - G4Track* track = step->GetTrack(); - G4int pid = track->GetParentID(); - G4int tid = track->GetTrackID(); - G4StepPoint* point1 = step->GetPreStepPoint(); - G4StepPoint* point2 = step->GetPostStepPoint(); - G4TouchableHandle touch1 = point1->GetTouchableHandle(); - G4TouchableHandle touch2 = point2->GetTouchableHandle(); - G4String vol1name = touch1->GetVolume()->GetName(); - G4String vol2name = touch2->GetVolume()->GetName(); - - G4String proc_name = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); - G4int copy_no = step->GetPostStepPoint()->GetTouchable()->GetReplicaNumber(1); - */ - - /* - // Retrieve the pointer to the optical boundary process. - // We do this only once per run defining our local pointer as static. - static G4OpBoundaryProcess* boundary = 0; - - if (!boundary) { // the pointer is not defined yet - // Get the list of processes defined for the optical photon - // and loop through it to find the optical boundary process. - G4ProcessVector* pv = pdef->GetProcessManager()->GetProcessList(); - for (G4int i=0; isize(); i++) { - if ((*pv)[i]->GetProcessName() == "OpBoundary") { - boundary = (G4OpBoundaryProcess*) (*pv)[i]; - break; - } - } - } - - if (step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) { - if (boundary->GetStatus() == Detection ){ - G4String detector_name = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName(); - G4cout << "##### Sensitive Volume: " << detector_name << G4endl; - - // detectorCounts::iterator it = my_counts_.find(detector_name); - // if (it != my_counts_.end()) my_counts_[it->first] += 1; - // else my_counts_[detector_name] = 1; - } - } - */ - return; - -} +// ----------------------------------------------------------------------------- +// G4Basic | SteppingAction.cpp +// +// User stepping action class. +// ----------------------------------------------------------------------------- + +#include + +#include "SteppingAction.h" +#include "EventAction.h" +#include "Analysis.h" + +#include +#include +//#include +#include +#include +#include +#include +#include +#include "G4UnitsTable.hh" +#include + + +using CLHEP::h_Planck; +using CLHEP::c_light; +using CLHEP::twopi; + +SteppingAction::SteppingAction( + EventAction* eventAction): + G4UserSteppingAction(), + fEventAction(eventAction), + counter(0), + ypos(-99999. * cm) +{ + G4cout << "SteppingAction::SteppingAction" << G4endl; + wl30_counters = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + wl30_counters_passed = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + wl45_counters = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + wl45_counters_passed = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; +} + + +SteppingAction::~SteppingAction() +{ + //G4double valuetoprint; + G4cout << "SteppingAction::~SteppingAction" << G4endl; + G4cout << "counter: " << counter << G4endl; + G4cout << "final y position: " << ypos << G4endl; + G4cout << "Final Transmittance fraction 30 deg: ["; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + //valuetoprint = 0.; + if (i==28) { + G4cout << ((wl30_counters[i]!=0) ? (double)wl30_counters_passed[i] / (double)wl30_counters[i] * 100. : 0.) << "]"; //wl30_counters[i] << "]"; + continue; + } + G4cout << ((wl30_counters[i]!=0) ? (double)wl30_counters_passed[i] / (double)wl30_counters[i] * 100. : 0.) << ", "; //wl30_counters[i] << ", "; + } + G4cout << G4endl; + + G4cout << "Final Transmittance fraction 45 deg: ["; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + if (i==28) { + G4cout << ((wl30_counters[i] != 0) ? (double)wl45_counters_passed[i] / (double)wl45_counters[i] * 100. : 0.) << "]"; //wl45_counters[i] << "] "; + continue; + } + G4cout << ((wl30_counters[i] != 0) ? (double)wl45_counters_passed[i] / (double)wl45_counters[i] * 100. : 0.) << ", "; //wl45_counters[i] << ", "; + } + G4cout << G4endl; + +} + + +void SteppingAction::UserSteppingAction(const G4Step* step) +{ + //if (step->GetTrack()->GetParentID() == 0) return; + + auto analysisManager = G4AnalysisManager::Instance(); + G4ParticleDefinition* pdef = step->GetTrack()->GetDefinition(); + + //Check whether the track is an optical photon + if (pdef != G4OpticalPhoton::Definition()) return; + + auto step_number = step->GetTrack()->GetCurrentStepNumber(); + auto volume_name = step->GetTrack()->GetVolume()->GetName(); + + std::vector theWavelengths = { 250, 295, 300, 305, 310, 315, 320, 325, 330, 335, 340, 345, 350, 355, 360, 365, 370, 375, 380, 385, 390, 395, 400, 405, 410, 415, 420, 425, 500 }; + G4double ekin_init = step->GetTrack()->GetVertexKineticEnergy(); + G4double wl_init = (h_Planck * c_light) / (ekin_init); + G4double momx = (step->GetTrack()->GetVertexMomentumDirection())[0]; + G4double momy = (step->GetTrack()->GetVertexMomentumDirection())[1]; + G4double momz = (step->GetTrack()->GetVertexMomentumDirection())[2]; + + G4StepPoint* postStep = step->GetPostStepPoint(); + + ypos = (postStep->GetPosition())[1]; + + + G4double theAngle = std::abs(std::atan(momy / (std::sqrt(momx * momx + momz * momz)))); + theAngle = theAngle / twopi * 360; + G4cout << "step_number: " << step_number << ", volume_name: " << volume_name << ", incident angle: " << (int)theAngle << ", current count for 250 nm: " << wl30_counters[0] << G4endl; + + G4cout << "Initial photon energy: " << G4BestUnit(wl_init, "Length") << G4endl; + G4cout << "ypos: " << ypos << G4endl; + + + G4double min_ang = 30.5; + G4double max_ang = 44.5; + G4cout << "incident wavelength without units: " << (int)((h_Planck * c_light) / (ekin_init) / nm) << G4endl; + if (step_number == 3) { + for (int i = 0; i < (int)wl30_counters.size(); i++) { + if ((int)((h_Planck * c_light) / (ekin_init) / nm) == theWavelengths[i]) { + //G4cout << "wave length entry: " << theWavelengths[i] << G4endl; + if (theAngle < min_ang) wl30_counters[i]++; + if (theAngle > max_ang) wl45_counters[i]++; + //G4cout << "current wl_counter for wavelength " << theWavelengths[i] << " nm: " << wl30_counters[i] << G4endl; + if (ypos < 0.) { + if (theAngle < min_ang) wl30_counters_passed[i]++; + if (theAngle > max_ang) wl45_counters_passed[i]++; + } + } + } + G4cout << "Temp counts 30 deg: "; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << wl30_counters[i] << " "; + } + G4cout << G4endl; + + G4cout << "Temp Transmittance counts 30 deg: "; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << wl30_counters_passed[i] << " "; + } + G4cout << G4endl; + + G4cout << "Temp Transmittance fraction 30 deg: ["; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << (double)wl30_counters_passed[i]/(double)wl30_counters[i]*100. << ", "; + } + G4cout << "]" << G4endl; + + G4cout << "Temp counts 45 deg: "; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << wl45_counters[i] << " "; + } + G4cout << G4endl; + + G4cout << "Temp Transmittance counts 45 deg: "; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << wl45_counters_passed[i] << " "; + } + G4cout << G4endl; + + G4cout << "Temp Transmittance fraction 45 deg: ["; + for (int i = 0; i < (int)wl30_counters.size(); i++) { + G4cout << (double)wl45_counters_passed[i] / (double)wl45_counters[i] * 100. << ", "; + } + G4cout << "]" << G4endl; + } + + + if (step_number==3 && volume_name=="WLS_PLATE") ++counter; + + if (step_number == 3) { + if (ypos < 0.) { + if (theAngle < min_ang) { + analysisManager->FillH1(0, wl_init); // transmittance + analysisManager->FillNtupleDColumn(0, wl_init); + } + if (theAngle > max_ang) { + analysisManager->FillH1(1, wl_init); // transmittance + analysisManager->FillNtupleDColumn(1, wl_init); + } + analysisManager->FillH1(2, ypos); // position distribution + analysisManager->FillH2(0, wl_init, theAngle); + } + analysisManager->FillNtupleDColumn(2, ypos); + analysisManager->FillNtupleDColumn(3, theAngle); + analysisManager->FillNtupleDColumn(4, wl_init); + analysisManager->AddNtupleRow(); + } + + + // example of information one can access about optical photons + /* + G4Track* track = step->GetTrack(); + G4int pid = track->GetParentID(); + G4int tid = track->GetTrackID(); + G4StepPoint* point1 = step->GetPreStepPoint(); + G4StepPoint* point2 = step->GetPostStepPoint(); + G4TouchableHandle touch1 = point1->GetTouchableHandle(); + G4TouchableHandle touch2 = point2->GetTouchableHandle(); + G4String vol1name = touch1->GetVolume()->GetName(); + G4String vol2name = touch2->GetVolume()->GetName(); + + G4String proc_name = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); + G4int copy_no = step->GetPostStepPoint()->GetTouchable()->GetReplicaNumber(1);*/ + + + /* + // Retrieve the pointer to the optical boundary process. + // We do this only once per run defining our local pointer as static. + static G4OpBoundaryProcess* boundary = 0; + + if (!boundary) { // the pointer is not defined yet + // Get the list of processes defined for the optical photon + // and loop through it to find the optical boundary process. + G4ProcessVector* pv = pdef->GetProcessManager()->GetProcessList(); + for (G4int i=0; isize(); i++) { + if ((*pv)[i]->GetProcessName() == "OpBoundary") { + boundary = (G4OpBoundaryProcess*) (*pv)[i]; + break; + } + } + } + + if (step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) { + if (boundary->GetStatus() == Detection ){ + G4String detector_name = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName(); + G4cout << "##### Sensitive Volume: " << detector_name << G4endl; + + // detectorCounts::iterator it = my_counts_.find(detector_name); + // if (it != my_counts_.end()) my_counts_[it->first] += 1; + // else my_counts_[detector_name] = 1; + } + } + */ + + return; + +} diff --git a/src/SteppingAction.h b/src/SteppingAction.h index 388a221..9a63703 100644 --- a/src/SteppingAction.h +++ b/src/SteppingAction.h @@ -1,25 +1,36 @@ -// ----------------------------------------------------------------------------- -// G4Basic | SteppingAction.h -// -// User stepping action class. -// ----------------------------------------------------------------------------- - -#ifndef STEPPING_ACTION_H -#define STEPPING_ACTION_H - -#include - -class G4Step; - - -class SteppingAction: public G4UserSteppingAction -{ -public: - SteppingAction(); - virtual ~SteppingAction(); - virtual void UserSteppingAction(const G4Step*); -private: - int counter; -}; - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | SteppingAction.h +// +// User stepping action class. +// ----------------------------------------------------------------------------- + +#ifndef STEPPING_ACTION_H +#define STEPPING_ACTION_H + +#include +#include + +#include + +class G4Step; +class EventAction; + + +class SteppingAction: public G4UserSteppingAction +{ +public: + SteppingAction(EventAction* eventAction); + virtual ~SteppingAction(); + virtual void UserSteppingAction(const G4Step*); +private: + EventAction* fEventAction; + int counter; + std::vector wl30_counters; + std::vector wl30_counters_passed; + std::vector wl45_counters; + std::vector wl45_counters_passed; + G4double ypos; + +}; + +#endif diff --git a/src/TrackingAction.cpp b/src/TrackingAction.cpp index 2873248..69d3130 100644 --- a/src/TrackingAction.cpp +++ b/src/TrackingAction.cpp @@ -1,15 +1,15 @@ -// ----------------------------------------------------------------------------- -// G4Basic | TrackingAction.cpp -// -// User run action class. -// ----------------------------------------------------------------------------- - -#include "TrackingAction.h" - - -void TrackingAction::PreUserTrackingAction(const G4Track*) -{} - - -void TrackingAction::PostUserTrackingAction(const G4Track*) -{} +// ----------------------------------------------------------------------------- +// G4Basic | TrackingAction.cpp +// +// User run action class. +// ----------------------------------------------------------------------------- + +#include "TrackingAction.h" + + +void TrackingAction::PreUserTrackingAction(const G4Track*) +{} + + +void TrackingAction::PostUserTrackingAction(const G4Track*) +{} diff --git a/src/TrackingAction.h b/src/TrackingAction.h index 58d0803..b5f68c4 100644 --- a/src/TrackingAction.h +++ b/src/TrackingAction.h @@ -1,27 +1,27 @@ -// ----------------------------------------------------------------------------- -// G4Basic | TrackingAction.h -// -// User run action class. -// ----------------------------------------------------------------------------- - -#ifndef TRACKING_ACTION_H -#define TRACKING_ACTION_H - -#include - -class G4Track; - - -class TrackingAction: public G4UserTrackingAction -{ -public: - TrackingAction(); - ~TrackingAction(); - void PreUserTrackingAction(const G4Track*); - void PostUserTrackingAction(const G4Track*); -}; - -inline TrackingAction::TrackingAction() {} -inline TrackingAction::~TrackingAction() {} - -#endif +// ----------------------------------------------------------------------------- +// G4Basic | TrackingAction.h +// +// User run action class. +// ----------------------------------------------------------------------------- + +#ifndef TRACKING_ACTION_H +#define TRACKING_ACTION_H + +#include + +class G4Track; + + +class TrackingAction: public G4UserTrackingAction +{ +public: + TrackingAction(); + ~TrackingAction(); + void PreUserTrackingAction(const G4Track*); + void PostUserTrackingAction(const G4Track*); +}; + +inline TrackingAction::TrackingAction() {} +inline TrackingAction::~TrackingAction() {} + +#endif diff --git a/src/params.h b/src/params.h new file mode 100644 index 0000000..aec196b --- /dev/null +++ b/src/params.h @@ -0,0 +1,3 @@ +#include + +G4double init_Ekin; \ No newline at end of file