diff --git a/README.md b/README.md index b9dbda1..e9be977 100644 --- a/README.md +++ b/README.md @@ -59,10 +59,13 @@ Run with custom config file * secondsipmenabled - Enable the second SiPM # Changelog +<<<<<<< HEAD ## 2020-04-08 * SiPM position fix * Changed scintillator mother volume +======= +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 ## 2020-02-19 * Scintillator subtracted from it's coating * Scintillation process fix 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/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/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 a1e4759..4a33d88 100644 --- a/src/SiPMDetectorConstruction.cc +++ b/src/SiPMDetectorConstruction.cc @@ -165,7 +165,11 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() { z_pos_helper = coatingThickness + (scintillator_size.getZ() / 2); z_pos_helper = z_pos_helper - (container_sizeZ / 2); +<<<<<<< HEAD //posScint = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y +======= + posScint = G4ThreeVector(0, 0, z_pos_helper); //center on X and Y +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 } //------------------------------------------------------------------------------------------------------------------ @@ -199,10 +203,16 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() //---World definitions---------------------------------------------------------------------------------------------- G4Box* solidWorld = new G4Box( "World", //its name +<<<<<<< HEAD world_sizeX, world_sizeY, world_sizeZ); //its size //Prevent overlapping so the world size is doubled +======= + 0.5 * world_sizeX, + 0.5 * world_sizeY, + 0.5 * world_sizeZ); //its size +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 G4LogicalVolume* logicWorld = new G4LogicalVolume( solidWorld, //its solid world_mat, //its material @@ -250,14 +260,24 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() 0.5 * scintillator_size.getZ()); //its size G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid +<<<<<<< HEAD scint_mat, //its material +======= + scint_mat, //its material +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 "Scintillator"); //its name G4VSolid* solidCoating = new G4SubtractionSolid("ScintillatorCoating", solidScintCoating, +<<<<<<< HEAD solidScint/*, 0, G4ThreeVector(0, 0, 0)*/); +======= + solidScint, + 0, + G4ThreeVector(0, 0, 0)); +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 G4LogicalVolume* logicScintCoating = new G4LogicalVolume( solidCoating, //its solid scint_coating, //its material @@ -267,7 +287,11 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() posScint, //at position logicScint, //its logical volume "Scintillator", //its name +<<<<<<< HEAD logicContainer, //its mother volume +======= + logicContainer, //its mother volume +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 false, //no boolean operation 0, //copy lxenumber checkOverlaps); //overlaps checking @@ -366,7 +390,11 @@ 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); +<<<<<<< HEAD scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 1000./MeV); //50 volt +======= + scint_material_mpt -> AddConstProperty("SCINTILLATIONYIELD", 50000./MeV); //50 volt +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 scint_material_mpt -> AddConstProperty("RESOLUTIONSCALE", 1.0); scint_material_mpt -> AddConstProperty("FASTTIMECONSTANT", 0.01*ns); scint_material_mpt -> AddConstProperty("SLOWTIMECONSTANT", 1.*ns); @@ -379,7 +407,11 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() G4LogicalBorderSurface *ScintillatorSurface = new G4LogicalBorderSurface("Scintillator Surface", physScint, physScintCoating, OpScintillatorSurface); +<<<<<<< HEAD G4double reflectivityCoating[n] = {0.99, 0.99}; +======= + G4double reflectivityCoating[n] = {0.9, 0.9}; +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 G4double efficiencyCoating[n] = {0, 0}; G4MaterialPropertiesTable *ScintillatorToWolframMaterialPropertyTable = new G4MaterialPropertiesTable(); @@ -429,6 +461,7 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() { snprintf(s1, 30, "Container_x%d_y%d", i, j); logicContainer -> SetName(s1); +<<<<<<< HEAD physContainer[i][j] = new G4PVPlacement( 0, G4ThreeVector(i * container_sizeX + (container_sizeX / 2) - (world_sizeX / 2), j * container_sizeY + (container_sizeY / 2) - (world_sizeY / 2), container_sizeZ / 2), logicContainer, @@ -438,6 +471,17 @@ G4VPhysicalVolume* SiPMDetectorConstruction::Construct() copyNumber, //copy number checkOverlaps); copyNumber++; +======= + physContainer[x][y] = new G4PVPlacement(0, + G4ThreeVector(i * (container_sizeX / 2), j * (container_sizeY / 2), container_sizeZ / 2), + logicContainer, + s1, //its name + logicWorld, + false, + helper, //copy number + checkOverlaps); + helper++; +>>>>>>> 247f4aad085b756a9fdbb30ec68ccec1ad1b97b7 } } //------------------------------------------------------------------------------------------------------------------ diff --git a/src/SiPMParameters.cc b/src/SiPMParameters.cc index ad133a6..e908ac8 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 = " << 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 = " << 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 = " << scintillator_Dimension.getZ() << 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); }