From 11540e0bf486f53e991a93dbc616882dcf257442 Mon Sep 17 00:00:00 2001 From: David Baranyai Date: Tue, 18 Feb 2020 15:14:15 +0100 Subject: [PATCH 1/2] SPACE --- README.md | 56 ++- config.conf | 27 +- include/SiPMParameters.hh | 52 ++- src/SiPMDetectorConstruction.cc | 567 ++++++++++++++++++------------ src/SiPMParameters.cc | 86 +++++ src/SiPMPrimaryGeneratorAction.cc | 9 +- 6 files changed, 553 insertions(+), 244 deletions(-) diff --git a/README.md b/README.md index 60f7acb..61ba08e 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ The output file is a ROOT file. Can be checked by TBrowser. * CERN Geant4 * CERN Root (tested on 6.19/01) * Linux or MacOS (should work on Windows, not tested) -* Also tested on WLS Ubuntu (X11 needed) +* Also tested on WLS Ubuntu (X11 (for e.g. VcXsrv) and OpenGL needed) ## Building Before building, be sure that you changed the macroPath to the right directory in sipm.cc. @@ -20,7 +20,53 @@ cmake ../SiPM make -jN (where N is the number of jobs to run simultaneously) ``` -## Updates -* 2020/02/06 - SiPMAnalisys' and SiPMParameters' GetInstance functions are returning a static reference instead of pointer -* 2020/02/06 - G4Mutex replaced with std::mutex -* 2020/02/06 - Input config file and output data file name can be changed with the GetInstance functions on the first call +## Running +Run with default parameters (not always works) +``` +./sipm +``` +Run with default config file (config.conf) +``` +./sipm -df +``` +Run with custom config file +``` +./sipm -f config_file_name.conf +``` + +## Config file parameters +* pgpositionx - Particle Gun X position +* pgpositiony - Particle Gun Y position +* pgpositionz - Particle Gun Z position +* pgmomentumx - Particle Gun X momentum +* pgmomentumy - Particle Gun Y momentum +* pgmomentumz - Particle Gun Z momentum +* sipmsizex - SiPM X Size +* sipmsizey - SiPM Y Size +* sipmsizez - SiPM Z Size +* scintillatorsizex - Scintillator X Size +* scintillatorsizey - Scintillator Y Size +* scintillatorlength - Scintillator Z Size +* scintillatorisbox - Scintillator can be a tube or a box +* coatingthickness - Scintillator coating thickness +* scintillatorradius - Radius of the scintillator. Only used when it is a tube. +* xdivision - How many detectors should be placed along the X axis +* ydivision - How many detectors should be placed along the Y axis +* pgenergy - Particle Gun energy +* numberofevents - Number of events +* lengthunit - Size unit in detector construction (mm | cm | m) +* firstsipmenabled - Enable the first SiPM +* secondsipmenabled - Enable the second SiPM + +# Changelog +## 2020-02-18 +* Both SiPMs can be enabled or disabled +* Added an option to change between a box and tube scintillator (only box tested yet) +* Added new parameters +* Cleaned up the detector construction +* Small changes according to the new parameters + +## 2020-02-06 +* SiPMAnalisys' and SiPMParameters' GetInstance functions are returning a static reference instead of pointer +* G4Mutex replaced with std::mutex +* Input config file and output data file name can be changed with the GetInstance functions on the first call diff --git a/config.conf b/config.conf index 38dd80c..4831330 100644 --- a/config.conf +++ b/config.conf @@ -1,15 +1,22 @@ -pgpositionx=50 -pgpositiony=-5 -pgpositionz=0 -pgmomentumx=0.5 +pgpositionx=0 +pgpositiony=-15 +pgpositionz=10 +pgmomentumx=0 pgmomentumy=1 pgmomentumz=0 sipmsizex=1 sipmsizey=1 sipmsizez=1 -scintillatorlength=40 -scintillatorradius=0.25 -xdivision=10 -ydivision=10 -pgenergy=1 -numberofevents=10 \ No newline at end of file +scintillatorsizex=1 +scintillatorsizey=1 +scintillatorlength=20 +scintillatorisbox=1 +coatingthickness=0.1 +scintillatorradius=5 +xdivision=1 +ydivision=1 +pgenergy=100 +numberofevents=100 +lengthunit=cm +firstsipmenabled=1 +secondsipmenabled=0 \ No newline at end of file diff --git a/include/SiPMParameters.hh b/include/SiPMParameters.hh index 71d981e..ee1985e 100644 --- a/include/SiPMParameters.hh +++ b/include/SiPMParameters.hh @@ -17,6 +17,7 @@ #include #include #include "G4ThreeVector.hh" +#include "G4SystemOfUnits.hh" class SiPMParameters { @@ -24,6 +25,30 @@ public: static SiPMParameters& GetInstance(const std::string& config_file_name = "config.conf"); ~SiPMParameters(); + //---Units--------------------------------------------------------------------------------------------- + enum lengthUnit { mm, cm, m }; + + void SetGlobalLengthUnit(lengthUnit l) { globalLengthUnit = l; } + double GetLengthMultiplier() //returns the correct length multiplier for the global length unit (use it instead of *cm/*mm/*m) + { + switch (globalLengthUnit) + { + case mm: + return millimeter; + break; + case cm: + return centimeter; + break; + case m: + return meter; + break; + default: + return centimeter; + break; + } + } + + //---Config file--------------------------------------------------------------------------------------- void ParseConfigFile(); //Read the default config file void ParseConfigFile(std::string config_file1); //Read another config file @@ -44,10 +69,18 @@ public: void SetSiPMSize(G4ThreeVector sipm_s) { sipm_Dimension = sipm_s; } void SetSiPMSize(G4double x, G4double y, G4double z) { sipm_Dimension.set(x, y, z); } G4ThreeVector GetSiPMSize() { return sipm_Dimension; } + void FirstSiPMEnabled(bool enabled) { sipm1Enabled = enabled; } + bool FirstSiPMEnabled() { return sipm1Enabled; } + void SecondSiPMEnabled(bool enabled) { sipm2Enabled = enabled; } + bool SecondSiPMEnabled() { return sipm2Enabled; } //---Scintillator parameters--------------------------------------------------------------------------- - void SetScintillatorLength(G4double sc_l) { scintillator_length = sc_l; } + void SetScintillatorLength(G4double sc_l) { scintillator_length = sc_l; scintillator_Dimension.setZ(sc_l); } G4double GetScintillatorLength() { return scintillator_length; } + void SetScintillatorSize(G4double x, G4double y, G4double z) { scintillator_Dimension.set(x, y, z); } + G4ThreeVector GetScintillatorSize() { return scintillator_Dimension; } + double CoatingThickness() { return cThickness; } + void CoatingThickness(double t) { cThickness = t; } //---Division parameters------------------------------------------------------------------------------- void SetDivision(G4int x, G4int y) { x_division = x; y_division = y; }; @@ -58,7 +91,11 @@ public: //---Radius parameters--------------------------------------------------------------------------------- void SetScintillatorRadius(G4double sc_r) { scint_radius = sc_r; } - G4double GetScintillatorRadius() { return scint_radius; } + G4double GetScintillatorRadius() { return scintIsBox ? 0.0 : scint_radius; } //If it is a box, return 0 + + //---Scintillator Tube vs box parameters--------------------------------------------------------------- + void ScintillatorIsBox(bool isBox) { scintIsBox = isBox; } + bool ScintillatorIsBox() { return scintIsBox; } //---Number of Events---------------------------------------------------------------------------------- void SetNumberOfEvents(G4int noe1) { numberofevents = noe1; } @@ -85,16 +122,21 @@ private: G4ThreeVector sipm_Dimension; //applies to both (in cm) G4double scintillator_length; //the size same as sipm + G4double scint_radius; + G4ThreeVector scintillator_Dimension; G4int x_division; G4int y_division; - G4double scint_radius; - G4int numberofevents; + bool scintIsBox = true; + bool sipm1Enabled = true; + bool sipm2Enabled = false; + + lengthUnit globalLengthUnit = cm; //default length is cm - + double cThickness = 1; }; #endif /* Parameters_hh */ diff --git a/src/SiPMDetectorConstruction.cc b/src/SiPMDetectorConstruction.cc index 563185a..846fd54 100644 --- a/src/SiPMDetectorConstruction.cc +++ b/src/SiPMDetectorConstruction.cc @@ -4,19 +4,22 @@ // // Created by Baranyai David on 2018. 08. 22.. // +// 2020.02.07 - SiPM size fixed +// - Scintillator is now a box #include "SiPMDetectorConstruction.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -SiPMDetectorConstruction::SiPMDetectorConstruction() -: G4VUserDetectorConstruction() -{ } +SiPMDetectorConstruction::SiPMDetectorConstruction() : G4VUserDetectorConstruction() +{ +} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... SiPMDetectorConstruction::~SiPMDetectorConstruction() -{ } +{ +} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -28,229 +31,341 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //Get the parameters instance SiPMParameters& parameters = SiPMParameters::GetInstance(); - // Option to switch on/off checking of volumes overlaps G4bool checkOverlaps = true; + + //---Object size calculations----------------------------------------------------------------------------------- + //Calculate everything with real lengths + double lengthMultiplier = parameters.GetLengthMultiplier(); + + G4ThreeVector sipm_size = parameters.GetSiPMSize(); + sipm_size.setX(sipm_size.getX() * lengthMultiplier); + sipm_size.setY(sipm_size.getY() * lengthMultiplier); + sipm_size.setZ(sipm_size.getZ() * lengthMultiplier); + + G4ThreeVector scintillator_size = parameters.GetScintillatorSize(); + scintillator_size.setX(scintillator_size.getX() * lengthMultiplier); + scintillator_size.setY(scintillator_size.getY() * lengthMultiplier); + scintillator_size.setZ(scintillator_size.getZ() * lengthMultiplier); + + //Scintillator radius | only used if the shape of the scintillator is set to tube + G4double scint_radius = parameters.GetScintillatorRadius() * lengthMultiplier; + + G4double coatingThickness = parameters.CoatingThickness() * lengthMultiplier; + + //Scintillator size with wolfram + G4ThreeVector coated_scintillator_size = scintillator_size; //scintillator size already calculated with the length mulitiplier + coated_scintillator_size.setX(coated_scintillator_size.getX() + coatingThickness); + coated_scintillator_size.setY(coated_scintillator_size.getY() + coatingThickness); + coated_scintillator_size.setZ(coated_scintillator_size.getZ() + (coatingThickness * 2)); //both sides + + //World size + int sipm_size_multiplier = 0; + if(parameters.FirstSiPMEnabled() || parameters.SecondSiPMEnabled()) sipm_size_multiplier = 1; //if at least one sipm is enabled, then increase the world size by sipm size * 1 + else if (parameters.FirstSiPMEnabled() && parameters.SecondSiPMEnabled()) sipm_size_multiplier = 2; //if both sipm is enabled, then increase the world size by sipm size * 2 - - // ------------- Materials ------------- - + G4double world_sizeZ = 0; + if (coatingThickness > sipm_size.getZ()) //if wolfram thicker than sipm, then increase the worldZ by the wolfram thickness + { + world_sizeZ = coated_scintillator_size.getZ(); + } + else + { + world_sizeZ = sipm_size_multiplier * sipm_size.getZ() + scintillator_size.getZ(); + } + + G4double world_sizeX = parameters.GetXDivison() * coated_scintillator_size.getX(); + G4double world_sizeY = parameters.GetYDivison() * coated_scintillator_size.getY(); + + //Container sizes + G4double container_sizeX = coated_scintillator_size.getX(); + G4double container_sizeY = coated_scintillator_size.getY(); + G4double container_sizeZ = world_sizeZ; + + //------------------------------------------------------------------------------------------------------------------ + + //---Object position calculations----------------------------------------------------------------------------------- + G4ThreeVector pos_sipm0 = G4ThreeVector(); + G4ThreeVector pos_sipm1 = G4ThreeVector(); + G4ThreeVector posScint = G4ThreeVector(); + G4ThreeVector posScintCoating = G4ThreeVector(); + double z_pos_helper = 0; + + //World position + G4ThreeVector posWorld = G4ThreeVector(); //at (0, 0, 0) + + //SiPM0 position + if (parameters.FirstSiPMEnabled()) + { + if (coatingThickness > sipm_size.getZ()) + { + z_pos_helper = coatingThickness - (sipm_size.getZ() / 2); + } + else + { + z_pos_helper = sipm_size.getZ() / 2; + } + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + pos_sipm0 = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + + //SiPM1 position + if (parameters.SecondSiPMEnabled()) + { + if (coatingThickness > sipm_size.getZ()) + { + z_pos_helper = coatingThickness + scintillator_size.getZ() + (sipm_size.getZ() / 2); + } + else + { + z_pos_helper = sipm_size.getZ() + scintillator_size.getZ() + sipm_size.getZ() / 2; + } + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + pos_sipm1 = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + + //Coating position + if (parameters.FirstSiPMEnabled()) + { + if (coatingThickness > sipm_size.getZ()) + { + z_pos_helper = coated_scintillator_size.getZ() / 2; + } + else + { + z_pos_helper = sipm_size.getZ() - coatingThickness + (coated_scintillator_size.getZ() / 2); + } + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + posScintCoating = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + else + { + z_pos_helper = coated_scintillator_size.getZ() / 2; + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + posScintCoating = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + + //Scintillator position + if (parameters.FirstSiPMEnabled()) + { + if (coatingThickness > sipm_size.getZ()) + { + z_pos_helper = coatingThickness + (scintillator_size.getZ() / 2); + } + else + { + z_pos_helper = sipm_size.getZ() + (scintillator_size.getZ() / 2); + } + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + posScint = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + else + { + z_pos_helper = coatingThickness + (scintillator_size.getZ() / 2); + z_pos_helper = z_pos_helper - (container_sizeZ / 2); + posScint = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y + } + + //------------------------------------------------------------------------------------------------------------------ + + //---Material definitions------------------------------------------------------------------------------------------- G4double a, z, density; G4int nelements; - + // Air - // - G4Element* N = new G4Element("Nitrogen", "N", z=7 , a=14.01*g/mole); - G4Element* O = new G4Element("Oxygen" , "O", z=8 , a=16.00*g/mole); - - G4Material* air = new G4Material("Air", density=1.29*mg/cm3, nelements=2); - air->AddElement(N, 70.*perCent); - air->AddElement(O, 30.*perCent); - - // Water - // - G4Element* H = new G4Element("Hydrogen", "H", z=1 , a=1.01*g/mole); - - G4Material* water = new G4Material("Water", density= 1.0*g/cm3, nelements=2); + G4Element* N = new G4Element("Nitrogen", "N", z = 7, a = 14.01 * g / mole); + G4Element* O = new G4Element("Oxygen", "O", z = 8, a = 16.00 * g / mole); + G4Material* air = new G4Material("Air", density = 1.29 * mg / cm3, nelements = 2); + air->AddElement(N, 70. * perCent); + air->AddElement(O, 30. * perCent); + + //Water + G4Element* H = new G4Element("Hydrogen", "H", z = 1, a = 1.01 * g / mole); + G4Material* water = new G4Material("Water", density = 1.0 * g / cm3, nelements = 2); water->AddElement(H, 2); water->AddElement(O, 1); + + G4Material* wolfram = nist->FindOrBuildMaterial("G4_W"); + + //Object materials + G4Material* world_mat = air; //nist->FindOrBuildMaterial("G4_AIR"); + G4Material* scint_mat = water; //nist->FindOrBuildMaterial("G4_Si"); + G4Material* scint_coating = wolfram; //nist->FindOrBuildMaterial("G4_W"); + G4Material* sipm0_mat = air; //nist->FindOrBuildMaterial("G4_Si"); + G4Material* sipm1_mat = air; //nist->FindOrBuildMaterial("G4_Si"); + //------------------------------------------------------------------------------------------------------------------ + + //---World definitions---------------------------------------------------------------------------------------------- + G4Box* solidWorld = new G4Box( "World", //its name + world_sizeX, + world_sizeY, + world_sizeZ); //its size - /* - * Wolfram material - */ - G4Material *wolfram = nist -> FindOrBuildMaterial("G4_W"); + G4LogicalVolume* logicWorld = new G4LogicalVolume( solidWorld, //its solid + world_mat, //its material + "World"); //its name + G4VPhysicalVolume* physWorld = new G4PVPlacement( 0, //no rotation + posWorld, //at (0,0,0) + logicWorld, //its logical volume + "World", //its name + 0, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + //------------------------------------------------------------------------------------------------------------------ + + //---Container definitions------------------------------------------------------------------------------------------ + //Place a container which contains everything (sipms, scintillator, coating) for G4PVPlacement + G4Box *solidContainer = new G4Box( "Container", //name + container_sizeX * 0.5, + container_sizeY * 0.5, + container_sizeZ * 0.5); //size - // - // World - // - G4ThreeVector sipm_size = parameters.GetSiPMSize(); - - G4double world_sizeX = parameters.GetXDivison() * sipm_size.getX() * cm; //2*m; - G4double world_sizeY = parameters.GetYDivison() * sipm_size.getY() * cm; //2*m; - G4double world_sizeZ = 2*sipm_size.getZ() + parameters.GetScintillatorLength(); - G4Material* world_mat = air; //nist->FindOrBuildMaterial("G4_AIR"); - - G4Box* solidWorld = - new G4Box("World", //its name - /*0.5**/world_sizeX, /*0.5**/world_sizeY, /*0.5**/world_sizeZ); //its size - - G4LogicalVolume* logicWorld = - new G4LogicalVolume(solidWorld, //its solid - world_mat, //its material - "World"); //its name - - G4VPhysicalVolume* physWorld = - new G4PVPlacement(0, //no rotation - G4ThreeVector(), //at (0,0,0) - logicWorld, //its logical volume - "World", //its name - 0, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - - //Place a container which contains everything for G4Replica - G4double container_sizeX = sipm_size.getX()*cm; - G4double container_sizeY = sipm_size.getY()*cm; - G4double container_sizeZ = (sipm_size.getZ()*2 + parameters.GetScintillatorLength())*cm; - - G4Box *solidContainer = - new G4Box("Container", container_sizeX*0.5, container_sizeY*0.5, container_sizeZ*0.5); - - G4LogicalVolume *logicContainer = - new G4LogicalVolume(solidContainer, world_mat, "Container"); + G4LogicalVolume *logicContainer = new G4LogicalVolume( solidContainer, //its solid + world_mat, //its material + "Container"); //its name G4Colour containerColour( 1.0, 1.0, 0.0); G4VisAttributes* containerVisAtt = new G4VisAttributes( containerColour ); - //logicContainer -> SetVisAttributes(containerVisAtt); - - G4double sizeX = sipm_size.getX()*cm; - G4double sizeY = sipm_size.getY()*cm; - G4double sipm_width = sipm_size.getZ()*cm; - - // - // Sipm0 - // - G4Material* sipm0_mat = air; //nist->FindOrBuildMaterial("G4_Si"); - - // sipm0 shape - G4double sipm0_sizeX = sizeX; - G4double sipm0_sizeY = sizeY; - G4double sipm0_sizeZ = sipm_width; - - G4ThreeVector pos_sipm0 = G4ThreeVector(0, 0*cm, (sipm_width/2)-container_sizeZ*0.5); - - G4Box* solidSipm0 = - new G4Box("Sipm0", //its name - 0.5*sipm0_sizeX, 0.5*sipm0_sizeY, - 0.5*sipm0_sizeZ); //its size - - G4LogicalVolume* logicSipm0 = - new G4LogicalVolume(solidSipm0, //its solid - sipm0_mat, //its material - "Sipm0"); //its name - - G4VPhysicalVolume *physSipm0 = - new G4PVPlacement(0, //no rotation - pos_sipm0, //at position - logicSipm0, //its logical volume - "Sipm0", //its name - logicContainer, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - - G4Colour sipmColour( 0.0, 1.0, 0.0); - G4VisAttributes* sipmVisAtt = new G4VisAttributes( sipmColour ); - logicSipm0->SetVisAttributes(sipmVisAtt); - - - // - // Box - // Changed the scintillator construction so now the box is Wolfram whic contains the scintillator - // - // - G4Material* scint_mat = water; //nist->FindOrBuildMaterial("G4_Si"); - - // box shape - G4double scint_sizeX = sizeX; - G4double scint_sizeY = sizeY; - G4double scint_sizeZ = parameters.GetScintillatorLength() * cm; - - G4double z_pos = sipm_width + (scint_sizeZ*0.5); - - G4ThreeVector posScint = G4ThreeVector(0, 0*cm, z_pos-container_sizeZ*0.5); - - G4Box* solidScint_W = - new G4Box("Scintillator_W", //its name - 0.5*scint_sizeX, 0.5*scint_sizeY, - 0.5*scint_sizeZ); //its size - - G4LogicalVolume* logicScint_W = - new G4LogicalVolume(solidScint_W, //its solid - wolfram, //its material - "Scintillator_W"); //its name - - G4VPhysicalVolume *physScint_W = - new G4PVPlacement(0, //no rotation - posScint, //at position - logicScint_W, //its logical volume - "Scintillator_W", //its name - logicContainer, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - + logicContainer -> SetVisAttributes(containerVisAtt); + logicContainer -> SetVisAttributes(G4VisAttributes::GetInvisible()); + //------------------------------------------------------------------------------------------------------------------ + + //---Scintillator coating definitions------------------------------------------------------------------------------- + G4Box* solidScintCoating = new G4Box( "ScintillatorCoating", //its name + 0.5 * coated_scintillator_size.getX(), + 0.5 * coated_scintillator_size.getY(), + 0.5 * coated_scintillator_size.getZ()); //its size + + G4LogicalVolume* logicScintCoating = new G4LogicalVolume( solidScintCoating, //its solid + scint_coating, //its material + "ScintillatorCoating"); //its name + + G4VPhysicalVolume* physScintCoating = new G4PVPlacement(0, //no rotation + posScintCoating, //at position + logicScintCoating, //its logical volume + "ScintillatorCoating", //its name + logicContainer, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + + + //---Scintillator definitions--------------------------------------------------------------------------------------- + /* - * Scintillator - */ - - G4double scint_radius = parameters.GetScintillatorRadius()*cm; - - G4Tubs * solidScint = new G4Tubs("tube", 0, scint_radius, 0.5*(scint_sizeZ+(0.5*mm)), 0, 2*CLHEP::pi); //name, inner R, outter R, Half length in Z, starting angle, angle of the segment in rad - new G4Box("Scintillator", - 0.5*scint_sizeX, - 0.5*scint_sizeY, - 0.5*scint_sizeZ); - - G4LogicalVolume *logicScint = - new G4LogicalVolume(solidScint, //its solid - scint_mat, //its material - "Scintillator"); //its name - - G4VPhysicalVolume *physScint = - new G4PVPlacement(0, //no rotation - G4ThreeVector(0,0,0), //at position (as the mother volume is not the world, maybe this will be the right place) - logicScint, //its logical volume - "Scintillator", //its name - logicScint_W, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - - G4Colour scintColour( 1.0, 0, 0.0); - G4VisAttributes* scintVisAtt = new G4VisAttributes( scintColour ); - logicScint -> SetVisAttributes(scintVisAtt); - - // - // Sipm1 - // - - G4double sipm2_pos_z = (0.5*sipm_width) + scint_sizeZ; - - G4Material* sipm1_mat = air; //nist->FindOrBuildMaterial("G4_Si"); - G4ThreeVector pos_sipm1 = G4ThreeVector(0, 0*cm, sipm_width+sipm2_pos_z-container_sizeZ*0.5); - + G4Tubs* solidScint = new G4Tubs("tube", //name + 0, //inner radius + scint_radius, //outter radius + scintillator_size.getZ(), //half length in Z + 0, //starting angle + 2 * CLHEP::pi); //angle of the segment in rad + + G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid + scint_mat, //its material + "Scintillator"); //its name + + G4VPhysicalVolume* physScint = new G4PVPlacement( 0, //no rotation + posScint, //at position + logicScint, //its logical volume + "Scintillator", //its name + logicScintCoating, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + */ + + G4Box* solidScint = new G4Box( "Scintillator", //its name + 0.5 * scintillator_size.getX(), + 0.5 * scintillator_size.getY(), + 0.5 * scintillator_size.getZ()); //its size + + G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid + scint_mat, //its material + "Scintillator"); //its name + + G4VPhysicalVolume* physScint = new G4PVPlacement( 0, //no rotation + posScint, //at position + logicScint, //its logical volume + "Scintillator", //its name + logicContainer, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + + G4Colour scintColour(1.0, 0, 0.0); + G4VisAttributes* scintVisAtt = new G4VisAttributes(scintColour); + logicScint->SetVisAttributes(scintVisAtt); + //------------------------------------------------------------------------------------------------------------------ + + //---SiPM0 definitions---------------------------------------------------------------------------------------------- + G4Colour sipmColour(0.0, 1.0, 0.0); + G4Box* solidSipm0; + G4LogicalVolume* logicSipm0; + G4VPhysicalVolume* physSipm0; + G4VisAttributes* sipmVisAtt; + + if (parameters.FirstSiPMEnabled()) + { + solidSipm0 = new G4Box( "Sipm0", //its name + 0.5 * sipm_size.getX(), + 0.5 * sipm_size.getY(), + 0.5 * sipm_size.getZ()); //its size + + logicSipm0 = new G4LogicalVolume( solidSipm0, //its solid + sipm0_mat, //its material + "Sipm0"); //its name + + physSipm0 = new G4PVPlacement( 0, //no rotation + pos_sipm0, //at position + logicSipm0, //its logical volume + "Sipm0", //its name + logicContainer, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + sipmVisAtt = new G4VisAttributes(sipmColour); + logicSipm0->SetVisAttributes(sipmVisAtt); + } + //------------------------------------------------------------------------------------------------------------------ + + //---SiPM1 definitions---------------------------------------------------------------------------------------------- // sipm1 shape -> same as sipm0 - G4Box* solidSipm1 = - new G4Box("Sipm1", //its name - 0.5*sipm0_sizeX, 0.5*sipm0_sizeY, - 0.5*sipm0_sizeZ); //its size + G4Box* solidSipm1; + G4LogicalVolume* logicSipm1; + G4VPhysicalVolume* physSipm1; + + if (parameters.SecondSiPMEnabled()) + { + solidSipm1 = new G4Box( "Sipm1", //its name + 0.5 * sipm_size.getX(), + 0.5 * sipm_size.getY(), + 0.5 * sipm_size.getZ()); //its size + + logicSipm1 = new G4LogicalVolume( solidSipm1, //its solid + sipm1_mat, //its material + "Sipm1"); //its name + + physSipm1 = new G4PVPlacement( 0, //no rotation + pos_sipm1, //at position + logicSipm1, //its logical volume + "Sipm1", //its name + logicContainer, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + + logicSipm1->SetVisAttributes(sipmVisAtt); + } + //------------------------------------------------------------------------------------------------------------------ - G4LogicalVolume* logicSipm1 = - new G4LogicalVolume(solidSipm1, //its solid - sipm1_mat, //its material - "Sipm1"); //its name - - G4VPhysicalVolume *physSipm1 = - new G4PVPlacement(0, //no rotation - pos_sipm1, //at position - logicSipm1, //its logical volume - "Sipm1", //its name - logicContainer, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - - logicSipm1->SetVisAttributes(sipmVisAtt); - - //----------------------------------General values-------------------------------------------------------- + //---General values------------------------------------------------------------------------------------------------- const G4int n = 2; - G4double pp[n] = {2.0*eV, 3.0*eV}; //distribution of optical photons produced in eV - - //---------------------------------General Material Settings---------------------------------------------- + G4double pp[n] = {2.0 * eV, 3.0 * eV}; //distribution of optical photons produced in eV + //------------------------------------------------------------------------------------------------------------------ + //---General Material Settings-------------------------------------------------------------------------------------- //material of scintillator G4double rind_scintillator[n] = {1.59, 1.57}; //refraction index G4double absl[n] = {35.*m, 35.*m}; //absorption length @@ -263,7 +378,7 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() scint_material_mpt -> AddProperty("ABSLENGTH", pp, absl, n); scint_material_mpt -> AddProperty("SLOWCOMPONENT", pp, slow, n); scint_material_mpt -> AddProperty("FASTCOMPONENT", pp, fast, n); - scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 50./MeV); //50 volt + scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 500./MeV); //50 volt scint_material_mpt -> AddConstProperty("RESOLUTIONSCALE", 1.0); scint_material_mpt -> AddConstProperty("FASTTIMECONSTANT", 0.01*ns); scint_material_mpt -> AddConstProperty("SLOWTIMECONSTANT", 1.*ns); @@ -274,14 +389,14 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //Surface of scintillator to wolfram G4OpticalSurface *OpScintillatorSurface = new G4OpticalSurface("Scintillator Surface to Wolfram",glisur, polished, dielectric_metal); - G4LogicalBorderSurface *ScintillatorSurface = new G4LogicalBorderSurface("Scintillator Surface", physScint, physScint_W, OpScintillatorSurface); + G4LogicalBorderSurface *ScintillatorSurface = new G4LogicalBorderSurface("Scintillator Surface", physScint, physScintCoating, OpScintillatorSurface); - G4double reflectivity_W[n] = {0.9, 0.9}; - G4double efficiency_W[n] = {0, 0}; + G4double reflectivityCoating[n] = {0.9, 0.9}; + G4double efficiencyCoating[n] = {0, 0}; G4MaterialPropertiesTable *ScintillatorToWolframMaterialPropertyTable = new G4MaterialPropertiesTable(); - ScintillatorToWolframMaterialPropertyTable -> AddProperty("REFLECTIVITY", pp, reflectivity_W, n); - ScintillatorToWolframMaterialPropertyTable -> AddProperty("EFFICIENCY", pp, efficiency_W, n); + ScintillatorToWolframMaterialPropertyTable -> AddProperty("REFLECTIVITY", pp, reflectivityCoating, n); + ScintillatorToWolframMaterialPropertyTable -> AddProperty("EFFICIENCY", pp, efficiencyCoating, n); OpScintillatorSurface -> SetMaterialPropertiesTable(ScintillatorToWolframMaterialPropertyTable); @@ -289,12 +404,19 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //Surface from scintillator to sipm0 and sipm1 G4OpticalSurface *SurfacefromScintillatorToSipm = new G4OpticalSurface("SurfacefromScintillatorToSipm", glisur, polished, dielectric_metal); - G4LogicalBorderSurface *SurfacefromScintillatorToSipm0_logical = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm0", physScint, physSipm0, SurfacefromScintillatorToSipm); + G4LogicalBorderSurface* SurfacefromScintillatorToSipm0_logical; + G4LogicalBorderSurface* SurfacefromScintillatorToSipm1_logical; + + if (parameters.FirstSiPMEnabled()) + { + SurfacefromScintillatorToSipm0_logical = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm0", physScint, physSipm0, SurfacefromScintillatorToSipm); + } + + if (parameters.SecondSiPMEnabled()) + { + SurfacefromScintillatorToSipm1_logical = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm1", physScint, physSipm1, SurfacefromScintillatorToSipm); + } - G4LogicalBorderSurface *SurfacefromScintillatorToSipm1_logical2 = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm1", physScint, physSipm1, SurfacefromScintillatorToSipm); - - - //---------------------------------- G4double reflectivity[n] = {1, 1}; //ha nem 1,1 akkor nem éri el a sipm-et az aki eléri a scint végét. G4MaterialPropertiesTable *ScintillatorMaterialPropertyTable = new G4MaterialPropertiesTable(); @@ -302,6 +424,9 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() SurfacefromScintillatorToSipm -> SetMaterialPropertiesTable(ScintillatorMaterialPropertyTable); + //------------------------------------------------------------------------------------------------------------------ + + //---Container placement-------------------------------------------------------------------------------------------- //Using G4PVPlacement instead of replica or others int x = parameters.GetXDivison(); @@ -317,7 +442,7 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() snprintf(s1, 30, "Container_x%d_y%d", i, j); logicContainer -> SetName(s1); physContainer[x][y] = new G4PVPlacement(0, - G4ThreeVector(i*container_sizeX, j*container_sizeY, 0), + G4ThreeVector(i * (container_sizeX / 2), j * (container_sizeY / 2), container_sizeZ / 2), logicContainer, s1, //its name logicWorld, @@ -327,6 +452,6 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() helper++; } } - + //------------------------------------------------------------------------------------------------------------------ return physWorld; } diff --git a/src/SiPMParameters.cc b/src/SiPMParameters.cc index ad133a6..4d4119a 100644 --- a/src/SiPMParameters.cc +++ b/src/SiPMParameters.cc @@ -70,6 +70,7 @@ void SiPMParameters::ParseConfigFile() } } CheckValues(); + conf_loaded = true; } catch (char param) { @@ -147,6 +148,7 @@ void SiPMParameters::StoreConfigValues(std::string key1, std::string value1) else if(key1.compare("scintillatorlength") == 0) { scintillator_length = std::stod(value1); + scintillator_Dimension.setZ(std::stod(value1)); std::cout << "Scintillator length parsed from config file! Value = " << scintillator_length << std::endl; } else if(key1.compare("scintillatorradius") == 0) @@ -166,6 +168,90 @@ void SiPMParameters::StoreConfigValues(std::string key1, std::string value1) y_division = std::stod(value1); std::cout << "Y division parsed from config file! Value = " << y_division << std::endl; } + + + else if (key1.compare("scintillatorsizex") == 0) + { + scintillator_Dimension.setX(std::stod(value1)); + std::cout << "Scintillator X size set from config file! Value = " << y_division << std::endl; + } + else if (key1.compare("scintillatorsizey") == 0) + { + scintillator_Dimension.setY(std::stod(value1)); + std::cout << "Scintillator Y size set from config file! Value = " << y_division << std::endl; + } + else if (key1.compare("scintillatorsizez") == 0) + { + scintillator_Dimension.setZ(std::stod(value1)); + std::cout << "Scintillator Z size set from config file! Value = " << y_division << std::endl; + } + + else if (key1.compare("scintillatorisbox") == 0) + { + if (std::stod(value1) != 0) + { + scintIsBox = true; + std::cout << "Using box scintillator." << std::endl; + } + else + { + scintIsBox = false; + std::cout << "Using tube scintillator." << std::endl; + } + } + + else if (key1.compare("firstsipmenabled") == 0) + { + if (std::stod(value1) != 0) + { + sipm1Enabled = true; + std::cout << "First SiPM enabled." << std::endl; + } + else + { + sipm1Enabled = false; + std::cout << "First SiPM disabled." << std::endl; + } + } + + else if (key1.compare("secondsipmenabled") == 0) + { + if (std::stod(value1) != 0) + { + sipm2Enabled = true; + std::cout << "Second SiPM enabled." << std::endl; + } + else + { + sipm2Enabled = false; + std::cout << "Second SiPM disabled." << std::endl; + } + } + + else if (key1.compare("coatingthickness") == 0) + { + cThickness = std::stod(value1); + std::cout << "Coating thickness set to " << cThickness << std::endl; + } + + else if (key1.compare("lengthunit") == 0) + { + if (value1.compare("mm") == 0) + { + std::cout << "Length unit is mm." << std::endl; + globalLengthUnit = mm; + } + else if (value1.compare("cm") == 0) + { + std::cout << "Length unit is cm." << std::endl; + globalLengthUnit = cm; + } + else if (value1.compare("m") == 0) + { + std::cout << "Length unit is m." << std::endl; + globalLengthUnit = m; + } + } //---Number of events--------------------------------------------------------------------------------------------------------- else if(key1.compare("numberofevents") == 0) diff --git a/src/SiPMPrimaryGeneratorAction.cc b/src/SiPMPrimaryGeneratorAction.cc index 695a7e6..552af37 100644 --- a/src/SiPMPrimaryGeneratorAction.cc +++ b/src/SiPMPrimaryGeneratorAction.cc @@ -16,10 +16,10 @@ SiPMPrimaryGeneratorAction::SiPMPrimaryGeneratorAction() : G4VUserPrimaryGenerat // default particle kinematic G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); G4String particleName; - G4ParticleDefinition* particle = particleTable->FindParticle(particleName="mu+"); + G4ParticleDefinition* particle = particleTable->FindParticle(particleName="gamma"); fParticleGun->SetParticleDefinition(particle); fParticleGun->SetParticleMomentumDirection(parameters.GetParticleGunMomentumDirection()); - fParticleGun->SetParticleEnergy(parameters.GetParticleGunEnergy()*GeV); //1GeV + fParticleGun->SetParticleEnergy(parameters.GetParticleGunEnergy()*keV); //1GeV } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -34,6 +34,9 @@ SiPMPrimaryGeneratorAction::~SiPMPrimaryGeneratorAction() void SiPMPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { SiPMParameters ¶meters = SiPMParameters::GetInstance(); - fParticleGun->SetParticlePosition(parameters.GetParticleGunPosition()); + G4ThreeVector pgPos = parameters.GetParticleGunPosition(); + double lengthMultiplier = parameters.GetLengthMultiplier(); + pgPos.set(pgPos.getX() * lengthMultiplier, pgPos.getY() * lengthMultiplier, pgPos.getZ() * lengthMultiplier); + fParticleGun->SetParticlePosition(pgPos); fParticleGun->GeneratePrimaryVertex(anEvent); } From 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 Mon Sep 17 00:00:00 2001 From: David Baranyai Date: Wed, 19 Feb 2020 11:51:42 +0100 Subject: [PATCH 2/2] Fixes and minor changes --- README.md | 6 +++ include/SiPMDetectorConstruction.hh | 1 + src/SiPMDetectorConstruction.cc | 71 ++++++++++++----------------- src/SiPMParameters.cc | 6 +-- 4 files changed, 39 insertions(+), 45 deletions(-) diff --git a/README.md b/README.md index 61ba08e..58626ad 100644 --- a/README.md +++ b/README.md @@ -59,6 +59,12 @@ Run with custom config file * secondsipmenabled - Enable the second SiPM # Changelog +## 2020-02-19 +* Scintillator subtracted from it's coating +* Scintillation process fix +* Size and position fixes +* Print fixes in Parameters + ## 2020-02-18 * Both SiPMs can be enabled or disabled * Added an option to change between a box and tube scintillator (only box tested yet) diff --git a/include/SiPMDetectorConstruction.hh b/include/SiPMDetectorConstruction.hh index 481a1c1..03ca0ea 100644 --- a/include/SiPMDetectorConstruction.hh +++ b/include/SiPMDetectorConstruction.hh @@ -31,6 +31,7 @@ #include "G4PVParameterised.hh" #include "SiPMParameterisation.hh" #include "SiPMParameters.hh" +#include "G4SubtractionSolid.hh" /// Detector construction class to define materials and geometry. diff --git a/src/SiPMDetectorConstruction.cc b/src/SiPMDetectorConstruction.cc index 846fd54..a7c4bfe 100644 --- a/src/SiPMDetectorConstruction.cc +++ b/src/SiPMDetectorConstruction.cc @@ -55,8 +55,8 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //Scintillator size with wolfram G4ThreeVector coated_scintillator_size = scintillator_size; //scintillator size already calculated with the length mulitiplier - coated_scintillator_size.setX(coated_scintillator_size.getX() + coatingThickness); - coated_scintillator_size.setY(coated_scintillator_size.getY() + coatingThickness); + coated_scintillator_size.setX(coated_scintillator_size.getX() + (coatingThickness * 2)); + coated_scintillator_size.setY(coated_scintillator_size.getY() + (coatingThickness * 2)); coated_scintillator_size.setZ(coated_scintillator_size.getZ() + (coatingThickness * 2)); //both sides //World size @@ -82,6 +82,8 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() G4double container_sizeY = coated_scintillator_size.getY(); G4double container_sizeZ = world_sizeZ; + + //------------------------------------------------------------------------------------------------------------------ //---Object position calculations----------------------------------------------------------------------------------- @@ -196,10 +198,10 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //------------------------------------------------------------------------------------------------------------------ //---World definitions---------------------------------------------------------------------------------------------- - G4Box* solidWorld = new G4Box( "World", //its name - world_sizeX, - world_sizeY, - world_sizeZ); //its size + G4Box* solidWorld = new G4Box( "World", //its name + 0.5 * world_sizeX, + 0.5 * world_sizeY, + 0.5 * world_sizeZ); //its size G4LogicalVolume* logicWorld = new G4LogicalVolume( solidWorld, //its solid world_mat, //its material @@ -238,44 +240,9 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() 0.5 * coated_scintillator_size.getY(), 0.5 * coated_scintillator_size.getZ()); //its size - G4LogicalVolume* logicScintCoating = new G4LogicalVolume( solidScintCoating, //its solid - scint_coating, //its material - "ScintillatorCoating"); //its name - - G4VPhysicalVolume* physScintCoating = new G4PVPlacement(0, //no rotation - posScintCoating, //at position - logicScintCoating, //its logical volume - "ScintillatorCoating", //its name - logicContainer, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - //---Scintillator definitions--------------------------------------------------------------------------------------- - /* - G4Tubs* solidScint = new G4Tubs("tube", //name - 0, //inner radius - scint_radius, //outter radius - scintillator_size.getZ(), //half length in Z - 0, //starting angle - 2 * CLHEP::pi); //angle of the segment in rad - - G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid - scint_mat, //its material - "Scintillator"); //its name - - G4VPhysicalVolume* physScint = new G4PVPlacement( 0, //no rotation - posScint, //at position - logicScint, //its logical volume - "Scintillator", //its name - logicScintCoating, //its mother volume - false, //no boolean operation - 0, //copy lxenumber - checkOverlaps); //overlaps checking - */ - G4Box* solidScint = new G4Box( "Scintillator", //its name 0.5 * scintillator_size.getX(), 0.5 * scintillator_size.getY(), @@ -285,6 +252,16 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() scint_mat, //its material "Scintillator"); //its name + G4VSolid* solidCoating = new G4SubtractionSolid("ScintillatorCoating", + solidScintCoating, + solidScint, + 0, + G4ThreeVector(0, 0, 0)); + + G4LogicalVolume* logicScintCoating = new G4LogicalVolume( solidCoating, //its solid + scint_coating, //its material + "ScintillatorCoating"); //its name + G4VPhysicalVolume* physScint = new G4PVPlacement( 0, //no rotation posScint, //at position logicScint, //its logical volume @@ -294,6 +271,16 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() 0, //copy lxenumber checkOverlaps); //overlaps checking + G4VPhysicalVolume* physScintCoating = new G4PVPlacement(0, //no rotation + posScintCoating, //at position + logicScintCoating, //its logical volume + "ScintillatorCoating", //its name + logicContainer, //its mother volume + false, //no boolean operation + 0, //copy lxenumber + checkOverlaps); //overlaps checking + + G4Colour scintColour(1.0, 0, 0.0); G4VisAttributes* scintVisAtt = new G4VisAttributes(scintColour); logicScint->SetVisAttributes(scintVisAtt); @@ -378,7 +365,7 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() scint_material_mpt -> AddProperty("ABSLENGTH", pp, absl, n); scint_material_mpt -> AddProperty("SLOWCOMPONENT", pp, slow, n); scint_material_mpt -> AddProperty("FASTCOMPONENT", pp, fast, n); - scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 500./MeV); //50 volt + scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 50000./MeV); //50 volt scint_material_mpt -> AddConstProperty("RESOLUTIONSCALE", 1.0); scint_material_mpt -> AddConstProperty("FASTTIMECONSTANT", 0.01*ns); scint_material_mpt -> AddConstProperty("SLOWTIMECONSTANT", 1.*ns); diff --git a/src/SiPMParameters.cc b/src/SiPMParameters.cc index 4d4119a..e908ac8 100644 --- a/src/SiPMParameters.cc +++ b/src/SiPMParameters.cc @@ -173,17 +173,17 @@ void SiPMParameters::StoreConfigValues(std::string key1, std::string value1) else if (key1.compare("scintillatorsizex") == 0) { scintillator_Dimension.setX(std::stod(value1)); - std::cout << "Scintillator X size set from config file! Value = " << y_division << std::endl; + std::cout << "Scintillator X size set from config file! Value = " << scintillator_Dimension.getX() << std::endl; } else if (key1.compare("scintillatorsizey") == 0) { scintillator_Dimension.setY(std::stod(value1)); - std::cout << "Scintillator Y size set from config file! Value = " << y_division << std::endl; + std::cout << "Scintillator Y size set from config file! Value = " << scintillator_Dimension.getY() << std::endl; } else if (key1.compare("scintillatorsizez") == 0) { scintillator_Dimension.setZ(std::stod(value1)); - std::cout << "Scintillator Z size set from config file! Value = " << y_division << std::endl; + std::cout << "Scintillator Z size set from config file! Value = " << scintillator_Dimension.getZ() << std::endl; } else if (key1.compare("scintillatorisbox") == 0)