divaldo
/
SiPM
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

SPACE

This commit is contained in:
David Baranyai 2020-02-18 15:14:15 +01:00
parent 9b2d6fcf08
commit 11540e0bf4
6 changed files with 553 additions and 244 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
README.md
View File

@ -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

27
config.conf
View File

@ -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
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

50
include/SiPMParameters.hh
View File

@ -17,6 +17,7 @@
#include <string>
#include <sstream>
#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 */

511
src/SiPMDetectorConstruction.cc
View File

@ -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();
// ------------- Materials -------------
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
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);
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);
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
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);
/*
* Wolfram material
*/
G4Material *wolfram = nist -> FindOrBuildMaterial("G4_W");
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
//
G4ThreeVector sipm_size = parameters.GetSiPMSize();
//---World definitions----------------------------------------------------------------------------------------------
G4Box* solidWorld = new G4Box( "World", //its name
world_sizeX,
world_sizeY,
world_sizeZ); //its size
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");
G4LogicalVolume* logicWorld = new G4LogicalVolume( solidWorld, //its solid
world_mat, //its material
"World"); //its name
G4Box* solidWorld =
new G4Box("World", //its name
/*0.5**/world_sizeX, /*0.5**/world_sizeY, /*0.5**/world_sizeZ); //its size
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
//------------------------------------------------------------------------------------------------------------------
G4LogicalVolume* logicWorld =
new G4LogicalVolume(solidWorld, //its solid
world_mat, //its material
"World"); //its name
//---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
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);
logicContainer -> SetVisAttributes(containerVisAtt);
logicContainer -> SetVisAttributes(G4VisAttributes::GetInvisible());
//------------------------------------------------------------------------------------------------------------------
G4double sizeX = sipm_size.getX()*cm;
G4double sizeY = sipm_size.getY()*cm;
G4double sipm_width = sipm_size.getZ()*cm;
//---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
//
// Sipm0
//
G4Material* sipm0_mat = air; //nist->FindOrBuildMaterial("G4_Si");
G4LogicalVolume* logicScintCoating = new G4LogicalVolume( solidScintCoating, //its solid
scint_coating, //its material
"ScintillatorCoating"); //its name
// 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);
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
//
// 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
//---Scintillator definitions---------------------------------------------------------------------------------------
/*
* Scintillator
*/
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
G4double scint_radius = parameters.GetScintillatorRadius()*cm;
G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid
scint_mat, //its material
"Scintillator"); //its name
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);
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
*/
G4LogicalVolume *logicScint =
new G4LogicalVolume(solidScint, //its solid
scint_mat, //its material
"Scintillator"); //its name
G4Box* solidScint = new G4Box( "Scintillator", //its name
0.5 * scintillator_size.getX(),
0.5 * scintillator_size.getY(),
0.5 * scintillator_size.getZ()); //its size
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
G4LogicalVolume* logicScint = new G4LogicalVolume( solidScint, //its solid
scint_mat, //its material
"Scintillator"); //its name
G4Colour scintColour( 1.0, 0, 0.0);
G4VisAttributes* scintVisAtt = new G4VisAttributes( scintColour );
logicScint -> SetVisAttributes(scintVisAtt);
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
//
// Sipm1
//
G4Colour scintColour(1.0, 0, 0.0);
G4VisAttributes* scintVisAtt = new G4VisAttributes(scintColour);
logicScint->SetVisAttributes(scintVisAtt);
//------------------------------------------------------------------------------------------------------------------
G4double sipm2_pos_z = (0.5*sipm_width) + scint_sizeZ;
//---SiPM0 definitions----------------------------------------------------------------------------------------------
G4Colour sipmColour(0.0, 1.0, 0.0);
G4Box* solidSipm0;
G4LogicalVolume* logicSipm0;
G4VPhysicalVolume* physSipm0;
G4VisAttributes* sipmVisAtt;
G4Material* sipm1_mat = air; //nist->FindOrBuildMaterial("G4_Si");
G4ThreeVector pos_sipm1 = G4ThreeVector(0, 0*cm, sipm_width+sipm2_pos_z-container_sizeZ*0.5);
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;
G4LogicalVolume* logicSipm1 =
new G4LogicalVolume(solidSipm1, //its solid
sipm1_mat, //its material
"Sipm1"); //its name
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
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 = new G4LogicalVolume( solidSipm1, //its solid
sipm1_mat, //its material
"Sipm1"); //its name
logicSipm1->SetVisAttributes(sipmVisAtt);
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
//----------------------------------General values--------------------------------------------------------
logicSipm1->SetVisAttributes(sipmVisAtt);
}
//------------------------------------------------------------------------------------------------------------------
//---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;
G4LogicalBorderSurface *SurfacefromScintillatorToSipm1_logical2 = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm1", physScint, physSipm1, SurfacefromScintillatorToSipm);
if (parameters.FirstSiPMEnabled())
{
SurfacefromScintillatorToSipm0_logical = new G4LogicalBorderSurface("SurfacefromScintillatorToSipm0", physScint, physSipm0, SurfacefromScintillatorToSipm);
}
if (parameters.SecondSiPMEnabled())
{
SurfacefromScintillatorToSipm1_logical = 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;
}

86
src/SiPMParameters.cc
View File

@ -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)
@ -167,6 +169,90 @@ void SiPMParameters::StoreConfigValues(std::string key1, std::string 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)
{

9
src/SiPMPrimaryGeneratorAction.cc
View File

@ -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 &parameters = 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);
}