SiPM position fix and changes on scintillator mother volume

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

config.conf

@@ -1,15 +1,22 @@
-pgpositionx=50
+pgpositionx=0
-pgpositiony=-5
+pgpositiony=-15
-pgpositionz=0
+pgpositionz=10
-pgmomentumx=0.5
+pgmomentumx=0
 pgmomentumy=1
 pgmomentumz=0
 sipmsizex=1
 sipmsizey=1
 sipmsizez=1
-scintillatorlength=40
+scintillatorsizex=1
-scintillatorradius=0.25
+scintillatorsizey=1
-xdivision=10
+scintillatorlength=20
-ydivision=10
+scintillatorisbox=1
-pgenergy=1
+coatingthickness=0.1
-numberofevents=10
+scintillatorradius=5
+xdivision=1
+ydivision=1
+pgenergy=100
+numberofevents=100
+lengthunit=cm
+firstsipmenabled=1
+secondsipmenabled=0

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.

include/SiPMParameters.hh

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

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
          }
      }
      //------------------------------------------------------------------------------------------------------------------

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)

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