Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-03-28 07:50:15

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file SAXSDetectorConstruction.cc
0027 /// \brief Implementation of the SAXSDetectorConstruction class
0028 
0029 #include "SAXSDetectorConstruction.hh"
0030 
0031 #include "SAXSSensitiveDetector.hh"
0032 
0033 #include "G4AssemblyVolume.hh"
0034 #include "G4Box.hh"
0035 #include "G4Colour.hh"
0036 #include "G4Cons.hh"
0037 #include "G4Element.hh"
0038 #include "G4ElementTable.hh"
0039 #include "G4ExtendedMaterial.hh"
0040 #include "G4GeometryManager.hh"
0041 #include "G4IntersectionSolid.hh"
0042 #include "G4LogicalVolume.hh"
0043 #include "G4LogicalVolumeStore.hh"
0044 #include "G4MIData.hh"
0045 #include "G4Material.hh"
0046 #include "G4MaterialTable.hh"
0047 #include "G4MultiFunctionalDetector.hh"
0048 #include "G4NistManager.hh"
0049 #include "G4PSDoseDeposit.hh"
0050 #include "G4PSEnergyDeposit.hh"
0051 #include "G4PVPlacement.hh"
0052 #include "G4PhysicalConstants.hh"
0053 #include "G4PhysicalVolumeStore.hh"
0054 #include "G4Polyhedra.hh"
0055 #include "G4RotationMatrix.hh"
0056 #include "G4RunManager.hh"
0057 #include "G4SDManager.hh"
0058 #include "G4SDParticleFilter.hh"
0059 #include "G4SolidStore.hh"
0060 #include "G4Sphere.hh"
0061 #include "G4SubtractionSolid.hh"
0062 #include "G4SystemOfUnits.hh"
0063 #include "G4ThreeVector.hh"
0064 #include "G4Trd.hh"
0065 #include "G4Tubs.hh"
0066 #include "G4UniformMagField.hh"
0067 #include "G4UnionSolid.hh"
0068 #include "G4VSensitiveDetector.hh"
0069 #include "G4VisAttributes.hh"
0070 #include "globals.hh"
0071 
0072 #include <cmath>
0073 #include <sstream>
0074 #include <vector>
0075 
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077 
0078 SAXSDetectorConstruction::SAXSDetectorConstruction() : G4VUserDetectorConstruction(), fWorldLogic(0)
0079 {
0080   G4cout << "### DetectorConstruction Instantiated ###" << G4endl;
0081 
0082   // instantiate the messenger (set methods will be called after construct)
0083   fMessenger = new SAXSDetectorConstructionMessenger(this);
0084 
0085   // set geometrical variables
0086   SetGeometricalVariables();
0087 
0088   // Initialization
0089   fPhantomMaterialIndex = 1;
0090 
0091   fComp0 = 0.0;  // components of the "Medical Material (MedMat)"
0092   fComp1 = 1.0;  // Through macro I can set one Medical Material only
0093   fComp2 = 0.0;
0094   fComp3 = 0.0;
0095 
0096   fCustomMatDensity = 1.00;  // g/mol
0097   fCustomMatHmassfract = 0.1119;
0098   fCustomMatCmassfract = 0.;
0099   fCustomMatNmassfract = 0.;
0100   fCustomMatOmassfract = 0.8881;
0101   fCustomMatNamassfract = 0.;
0102   fCustomMatPmassfract = 0.;
0103   fCustomMatSmassfract = 0.;
0104   fCustomMatClmassfract = 0.;
0105   fCustomMatKmassfract = 0.;
0106   fCustomMatCamassfract = 0.;
0107 
0108   fCustomMatFF = "";  // MIFF filename for a custom (extended) material
0109 
0110   fSensitiveVolume = 0;
0111 
0112   fIWantSlits = false;
0113 }
0114 
0115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0116 
0117 SAXSDetectorConstruction::~SAXSDetectorConstruction() {}
0118 
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120 
0121 void SAXSDetectorConstruction::DefineMaterials()
0122 {
0123   // Define the NIST manager
0124   G4NistManager* NistMan = G4NistManager::Instance();
0125 
0126   // Define the required elements for compounds
0127   G4Element* elH = NistMan->FindOrBuildElement("H");
0128   G4Element* elC = NistMan->FindOrBuildElement("C");
0129   G4Element* elN = NistMan->FindOrBuildElement("N");
0130   G4Element* elO = NistMan->FindOrBuildElement("O");
0131   G4Element* elNa = NistMan->FindOrBuildElement("Na");
0132   G4Element* elP = NistMan->FindOrBuildElement("P");
0133   G4Element* elS = NistMan->FindOrBuildElement("S");
0134   G4Element* elCl = NistMan->FindOrBuildElement("Cl");
0135   G4Element* elK = NistMan->FindOrBuildElement("K");
0136   G4Element* elCa = NistMan->FindOrBuildElement("Ca");
0137 
0138   // variable definition
0139   G4double d;  // density
0140   G4int nel;  // number of elements
0141   G4String matname;
0142 
0143   // Air
0144   d = 1.29 * mg / cm3;
0145   nel = 2;
0146   G4double tAir = 293.15 * CLHEP::kelvin;  // 20° Celsius
0147   G4double pAir = 1. * CLHEP::atmosphere;  // 1 atm
0148   fAir = new G4Material("Air", d, nel, kStateGas, tAir, pAir);
0149   fAir->AddElement(elN, 0.7);
0150   fAir->AddElement(elO, 0.3);
0151 
0152   // Fat (Tartari2002) (FF from Tartari2002)
0153   G4double d_Fat = 0.923 * g / cm3;
0154   nel = 3;
0155   matname = "Fat_MI";
0156   fFat = new G4Material(matname, d_Fat, nel);
0157   fFat->AddElement(elH, 0.119);
0158   fFat->AddElement(elC, 0.772);
0159   fFat->AddElement(elO, 0.109);
0160 
0161   // Water (FF from Tartari2002)
0162   G4double d_Water = 1. * g / cm3;
0163   nel = 2;
0164   matname = "Water_MI";
0165   fWater = new G4Material(matname, d_Water, nel);
0166   fWater->AddElement(elH, 2);
0167   fWater->AddElement(elO, 1);
0168 
0169   // BoneMatrix (Collagen) (FF from Tartari2002)
0170   G4double d_BoneMatrix = 1.263 * g / cm3;
0171   nel = 4;
0172   matname = "BoneMatrix_MI";
0173   fBoneMatrix = new G4Material(matname, d_BoneMatrix, nel);
0174   fBoneMatrix->AddElement(elH, 0.0344);
0175   fBoneMatrix->AddElement(elC, 0.7140);
0176   fBoneMatrix->AddElement(elN, 0.1827);
0177   fBoneMatrix->AddElement(elO, 0.0689);
0178 
0179   // Mineral (Hydroxyapatite) (Tartari2002) (FF from Tartari2002)
0180   G4double d_Mineral = 2.74 * g / cm3;
0181   nel = 4;
0182   matname = "Mineral_MI";
0183   fMineral = new G4Material(matname, d_Mineral, nel);
0184   fMineral->AddElement(elH, 0.002);
0185   fMineral->AddElement(elO, 0.414);
0186   fMineral->AddElement(elP, 0.185);
0187   fMineral->AddElement(elCa, 0.399);
0188 
0189   // Medical Material (compostion of Water, Fat, BoneMatrix and Mineral)
0190   G4double comp[] = {fComp0, fComp1, fComp2, fComp3};
0191   G4double d_MedMat =
0192     1 / (comp[0] / d_Fat + comp[1] / d_Water + comp[2] / d_BoneMatrix + comp[3] / d_Mineral);
0193   G4int n_MedMat = 0;
0194   for (size_t i = 0; i < 4; i++) {
0195     if (comp[i] > 0) n_MedMat++;
0196     if (comp[i] < 0 || comp[i] > 1) {
0197       G4String excep = "Error in Medical Material composition: comp[i]<0 or comp[i]>1";
0198       G4Exception("DetectorConstuction::DefineMaterials()", "dc0001", FatalException, excep);
0199       return;
0200     }
0201   }
0202   std::stringstream ss0, ss1, ss2, ss3;
0203   ss0 << comp[0];
0204   ss1 << comp[1];
0205   ss2 << comp[2];
0206   ss3 << comp[3];
0207   if (comp[0] == 0 || comp[0] == 1) ss0 << ".00";
0208   if (comp[1] == 0 || comp[1] == 1) ss1 << ".00";
0209   if (comp[2] == 0 || comp[2] == 1) ss2 << ".00";
0210   if (comp[3] == 0 || comp[3] == 1) ss3 << ".00";
0211   if (ss0.str().size() < 4) ss0 << "0";
0212   if (ss1.str().size() < 4) ss1 << "0";
0213   if (ss2.str().size() < 4) ss2 << "0";
0214   if (ss3.str().size() < 4) ss3 << "0";
0215   if (ss0.str().size() != 4 || ss1.str().size() != 4 || ss2.str().size() != 4
0216       || ss3.str().size() != 4)
0217   {
0218     G4String excep = "Error in MedMaterial composition: check the digits of the elements of comp";
0219     G4Exception("DetectorConstuction::DefineMaterials()", "dc0002", FatalException, excep);
0220     return;
0221   }
0222   matname = "MedMat_" + ss0.str() + "_" + ss1.str() + "_" + ss2.str() + "_" + ss3.str();
0223   fMedMat = new G4Material(matname, d_MedMat, n_MedMat);
0224   if (comp[0]) fMedMat->AddMaterial(fFat, comp[0]);
0225   if (comp[1]) fMedMat->AddMaterial(fWater, comp[1]);
0226   if (comp[2]) fMedMat->AddMaterial(fBoneMatrix, comp[2]);
0227   if (comp[3]) fMedMat->AddMaterial(fMineral, comp[3]);
0228   if (comp[0] + comp[1] + comp[2] + comp[3] != 1) {
0229     G4String excep = "Error in Medical Material composition: sum(comp) != 1";
0230     G4Exception("DetectorConstuction::DefineMaterials()", "dc0003", FatalException, excep);
0231     return;
0232   }
0233   // If the user wants to use more than one MedMat, he has to create the mix
0234   // and label the material properly, such as "MedMat_0.55_0.25_0.05_0.15".
0235   // Such a name enables the automatic form factor calculation.
0236 
0237   // PMMA (FF from Tartari2002)
0238   d = 1.18 * g / cm3;
0239   nel = 3;
0240   matname = "PMMA_MI";
0241   fPMMA = new G4Material(matname, d, nel);
0242   fPMMA->AddElement(elH, 8);
0243   fPMMA->AddElement(elC, 5);
0244   fPMMA->AddElement(elO, 2);
0245 
0246   // Adipose (Poletti2002) (FF from Poletti2002)
0247   d = 0.92 * g / cm3;
0248   nel = 4;
0249   matname = "adipose_MI";
0250   fAdipose = new G4Material(matname, d, nel);
0251   fAdipose->AddElement(elH, 0.124);
0252   fAdipose->AddElement(elC, 0.765);
0253   fAdipose->AddElement(elN, 0.004);
0254   fAdipose->AddElement(elO, 0.107);
0255 
0256   // Glandular (Poletti2002) (FF from Poletti2002)
0257   d = 1.04 * g / cm3;
0258   nel = 4;
0259   matname = "glandular_MI";
0260   fGlandular = new G4Material(matname, d, nel);
0261   fGlandular->AddElement(elH, 0.093);
0262   fGlandular->AddElement(elC, 0.184);
0263   fGlandular->AddElement(elN, 0.044);
0264   fGlandular->AddElement(elO, 0.679);
0265 
0266   // human breast 50/50 (ICRU44) (FF from Peplow1998)
0267   d = 0.96 * g / cm3;
0268   nel = 3;
0269   matname = "breast5050_MI";
0270   fBreast5050 = new G4Material(matname, d, nel);
0271   fBreast5050->AddElement(elH, 0.115);
0272   fBreast5050->AddElement(elC, 0.387);
0273   fBreast5050->AddElement(elO, 0.498);
0274 
0275   // Liver (ICRU46) (FF_pork_liver_Peplow1998)
0276   d = 1.06 * g / cm3;
0277   nel = 9;
0278   matname = "liver_MI";
0279   fliver = new G4Material(matname, d, nel);
0280   fliver->AddElement(elH, 0.102);
0281   fliver->AddElement(elC, 0.139);
0282   fliver->AddElement(elN, 0.030);
0283   fliver->AddElement(elO, 0.716);
0284   fliver->AddElement(elNa, 0.002);
0285   fliver->AddElement(elP, 0.003);
0286   fliver->AddElement(elS, 0.003);
0287   fliver->AddElement(elCl, 0.002);
0288   fliver->AddElement(elK, 0.003);
0289 
0290   // Kidney (ICRU46) (FF_pork_kidney_Peplow1998)
0291   d = 1.05 * g / cm3;
0292   nel = 10;
0293   matname = "kidney_MI";
0294   fkidney = new G4Material(matname, d, nel);
0295   fkidney->AddElement(elH, 0.103);
0296   fkidney->AddElement(elC, 0.132);
0297   fkidney->AddElement(elN, 0.030);
0298   fkidney->AddElement(elO, 0.724);
0299   fkidney->AddElement(elNa, 0.002);
0300   fkidney->AddElement(elP, 0.002);
0301   fkidney->AddElement(elS, 0.002);
0302   fkidney->AddElement(elCl, 0.002);
0303   fkidney->AddElement(elK, 0.002);
0304   fkidney->AddElement(elCa, 0.001);
0305 
0306   // Lexan (Polycarbonate) (FF from Peplow1998)
0307   d = 1.221 * g / cm3;
0308   nel = 3;
0309   matname = "Lexan_MI";
0310   fLexan = new G4Material(matname, d, nel);
0311   fLexan->AddElement(elH, 14);
0312   fLexan->AddElement(elC, 16);
0313   fLexan->AddElement(elO, 3);
0314 
0315   // Kapton (FF from Peplow1998)
0316   d = 1.42 * g / cm3;
0317   nel = 4;
0318   matname = "Kapton_MI";
0319   fKapton = new G4Material(matname, d, nel);
0320   fKapton->AddElement(elH, 28);
0321   fKapton->AddElement(elC, 35);
0322   fKapton->AddElement(elN, 2);
0323   fKapton->AddElement(elO, 7);
0324 
0325   // Carcinoma (muscle ICRU44) (FF from Kidane1999)
0326   d = 1.05 * g / cm3;  // check the density
0327   nel = 9;
0328   matname = "carcinoma_MI";
0329   fcarcinoma = new G4Material(matname, d, nel);
0330   fcarcinoma->AddElement(elH, 0.102);
0331   fcarcinoma->AddElement(elC, 0.143);
0332   fcarcinoma->AddElement(elN, 0.034);
0333   fcarcinoma->AddElement(elO, 0.710);
0334   fcarcinoma->AddElement(elNa, 0.001);
0335   fcarcinoma->AddElement(elP, 0.002);
0336   fcarcinoma->AddElement(elS, 0.003);
0337   fcarcinoma->AddElement(elCl, 0.001);
0338   fcarcinoma->AddElement(elK, 0.004);
0339 
0340   // Nylon (FF from Kosanetzky1987)
0341   d = 1.15 * g / cm3;
0342   nel = 4;
0343   matname = "Nylon_MI";
0344   fNylon = new G4Material(matname, d, nel);
0345   fNylon->AddElement(elH, 11);
0346   fNylon->AddElement(elC, 6);
0347   fNylon->AddElement(elN, 1);
0348   fNylon->AddElement(elO, 1);
0349 
0350   // Polyethylene (FF from Kosanetzky1987)
0351   d = 0.94 * g / cm3;  // MDPE => 0.92*g/cm3, HDPE => 0.94*g/cm3
0352   nel = 2;
0353   matname = "Polyethylene_MI";
0354   fPolyethylene = new G4Material(matname, d, nel);
0355   fPolyethylene->AddElement(elH, 4);
0356   fPolyethylene->AddElement(elC, 2);
0357 
0358   // Polystyrene (FF from Kosanetzky1987)
0359   d = 1.05 * g / cm3;
0360   nel = 2;
0361   matname = "Polystyrene_MI";
0362   fPolystyrene = new G4Material(matname, d, nel);
0363   fPolystyrene->AddElement(elH, 8);
0364   fPolystyrene->AddElement(elC, 8);
0365 
0366   // GrayMatter (DeFelici2008) (FF from DeFelici2008)
0367   d = 0.991 * g / cm3;
0368   nel = 3;
0369   matname = "grayMatter_MI";
0370   fGrayMatter = new G4Material(matname, d, nel);
0371   fGrayMatter->AddElement(elH, 0.1127);
0372   fGrayMatter->AddElement(elC, 0.0849);
0373   fGrayMatter->AddElement(elO, 0.8024);
0374 
0375   // WhiteMatter (DeFelici2008) (FF from DeFelici2008)
0376   d = 0.983 * g / cm3;
0377   nel = 3;
0378   matname = "whiteMatter_MI";
0379   fWhiteMatter = new G4Material(matname, d, nel);
0380   fWhiteMatter->AddElement(elH, 0.1134);
0381   fWhiteMatter->AddElement(elC, 0.1621);
0382   fWhiteMatter->AddElement(elO, 0.7245);
0383 
0384   // Blood (beef) (ICRU46) (FF from Peplow1998)
0385   d = 1.06 * g / cm3;
0386   nel = 9;
0387   matname = "blood_MI";
0388   fbeefBlood = new G4Material(matname, d, nel);
0389   fbeefBlood->AddElement(elH, 0.102);
0390   fbeefBlood->AddElement(elC, 0.11);
0391   fbeefBlood->AddElement(elN, 0.033);
0392   fbeefBlood->AddElement(elO, 0.746);
0393   fbeefBlood->AddElement(elNa, 0.001);
0394   fbeefBlood->AddElement(elP, 0.001);
0395   fbeefBlood->AddElement(elS, 0.002);
0396   fbeefBlood->AddElement(elCl, 0.003);
0397   fbeefBlood->AddElement(elK, 0.002);
0398 
0399   // Formaline (FF from Peplow1998)
0400   d = 1.083 * g / cm3;
0401   nel = 3;
0402   matname = "Formaline_MI";
0403   fFormaline = new G4Material(matname, d, nel);
0404   fFormaline->AddElement(elH, 2);
0405   fFormaline->AddElement(elC, 1);
0406   fFormaline->AddElement(elO, 1);
0407 
0408   // Acetone (FF from Cozzini2010)
0409   d = 0.7845 * g / cm3;
0410   nel = 3;
0411   matname = "Acetone_MI";
0412   fAcetone = new G4Material(matname, d, nel);
0413   fAcetone->AddElement(elH, 6);
0414   fAcetone->AddElement(elC, 3);
0415   fAcetone->AddElement(elO, 1);
0416 
0417   // Hperoxide (FF from Cozzini2010)
0418   d = 1.11 * g / cm3;
0419   nel = 2;
0420   matname = "Hperoxide_MI";
0421   fHperoxide = new G4Material(matname, d, nel);
0422   fHperoxide->AddElement(elH, 2);
0423   fHperoxide->AddElement(elO, 2);
0424 
0425   // CIRS30-70 (Poletti2002) (FF from Poletti2002)
0426   d = 0.97 * g / cm3;
0427   nel = 5;
0428   matname = "CIRS30-70_MI";
0429   fCIRS3070 = new G4Material(matname, d, nel);
0430   fCIRS3070->AddElement(elH, 0.1178);
0431   fCIRS3070->AddElement(elC, 0.7512);
0432   fCIRS3070->AddElement(elN, 0.0066);
0433   fCIRS3070->AddElement(elO, 0.1214);
0434   fCIRS3070->AddElement(elCa, 0.0030);
0435 
0436   // CIRS50-50 (Poletti2002) (FF from Poletti2002)
0437   d = 0.98 * g / cm3;
0438   nel = 5;
0439   matname = "CIRS50-50_MI";
0440   fCIRS5050 = new G4Material(matname, d, nel);
0441   fCIRS5050->AddElement(elH, 0.1110);
0442   fCIRS5050->AddElement(elC, 0.7274);
0443   fCIRS5050->AddElement(elN, 0.0104);
0444   fCIRS5050->AddElement(elO, 0.1482);
0445   fCIRS5050->AddElement(elCa, 0.0030);
0446 
0447   // CIRS70-30 (Poletti2002) (FF from Poletti2002)
0448   d = 1.01 * g / cm3;
0449   nel = 5;
0450   matname = "CIRS70-30_MI";
0451   fCIRS7030 = new G4Material(matname, d, nel);
0452   fCIRS7030->AddElement(elH, 0.1172);
0453   fCIRS7030->AddElement(elC, 0.7378);
0454   fCIRS7030->AddElement(elN, 0.0130);
0455   fCIRS7030->AddElement(elO, 0.1244);
0456   fCIRS7030->AddElement(elCa, 0.0076);
0457 
0458   // RMI454 (Poletti2002) (FF from Poletti2002)
0459   d = 0.98 * g / cm3;
0460   nel = 4;
0461   matname = "RMI454_MI";
0462   fRMI454 = new G4Material(matname, d, nel);
0463   fRMI454->AddElement(elH, 0.0924);
0464   fRMI454->AddElement(elC, 0.6935);
0465   fRMI454->AddElement(elN, 0.0198);
0466   fRMI454->AddElement(elO, 0.1943);
0467 
0468   // Bone (King2011 decomposition) (FF from King2011)
0469   d = 1.344 * g / cm3;
0470   nel = 6;
0471   matname = "bone_MI";
0472   fBone = new G4Material(matname, d, nel);
0473   fBone->AddElement(elH, 0.0582);
0474   fBone->AddElement(elC, 0.3055);
0475   fBone->AddElement(elN, 0.0347);
0476   fBone->AddElement(elO, 0.3856);
0477   fBone->AddElement(elP, 0.0684);
0478   fBone->AddElement(elCa, 0.1476);
0479 
0480   // FatLowX (Tartari2002) (FF_fat_Tartari2002_joint_lowXdata_ESRF2003)
0481   nel = 3;
0482   matname = "FatLowX_MI";
0483   ffatLowX = new G4Material(matname, d_Fat, nel);
0484   ffatLowX->AddElement(elH, 0.119);
0485   ffatLowX->AddElement(elC, 0.772);
0486   ffatLowX->AddElement(elO, 0.109);
0487 
0488   // BonematrixLowX (Collagen)
0489   //(Tartari2002) (FF_bonematrix_Tartari2002_joint_lowXdata)
0490   nel = 4;
0491   matname = "BoneMatrixLowX_MI";
0492   fbonematrixLowX = new G4Material(matname, d_BoneMatrix, nel);
0493   fbonematrixLowX->AddElement(elH, 0.0344);
0494   fbonematrixLowX->AddElement(elC, 0.7140);
0495   fbonematrixLowX->AddElement(elN, 0.1827);
0496   fbonematrixLowX->AddElement(elO, 0.0689);
0497 
0498   // dryBoneLowX
0499   //(Tartari2002) (FF_dryBone_Tartari2002_joint_lowXdata_ESRF2003)
0500   d = 2.06 * g / cm3;
0501   nel = 6;
0502   matname = "dryBoneLowX_MI";
0503   fdryBoneLowX = new G4Material(matname, d, nel);
0504   fdryBoneLowX->AddElement(elH, 0.0112);
0505   fdryBoneLowX->AddElement(elC, 0.2013);
0506   fdryBoneLowX->AddElement(elN, 0.0515);
0507   fdryBoneLowX->AddElement(elO, 0.3148);
0508   fdryBoneLowX->AddElement(elP, 0.1327);
0509   fdryBoneLowX->AddElement(elCa, 0.2885);
0510 
0511   // CustomMat (FF read from file)
0512   nel = 0;
0513   G4double sumMF = 0.;
0514   G4cout << "CustomMat composition: " << G4endl;
0515   if (fCustomMatHmassfract) {
0516     G4cout << "CustomMatHmassfract: " << fCustomMatHmassfract << G4endl;
0517     nel++;
0518     sumMF += fCustomMatHmassfract;
0519   }
0520   if (fCustomMatCmassfract) {
0521     G4cout << "CustomMatCmassfract: " << fCustomMatCmassfract << G4endl;
0522     nel++;
0523     sumMF += fCustomMatCmassfract;
0524   }
0525   if (fCustomMatNmassfract) {
0526     G4cout << "CustomMatNmassfract: " << fCustomMatNmassfract << G4endl;
0527     nel++;
0528     sumMF += fCustomMatNmassfract;
0529   }
0530   if (fCustomMatOmassfract) {
0531     G4cout << "CustomMatOmassfract: " << fCustomMatOmassfract << G4endl;
0532     nel++;
0533     sumMF += fCustomMatOmassfract;
0534   }
0535   if (fCustomMatNamassfract) {
0536     G4cout << "CustomMatNamassfract: " << fCustomMatNamassfract << G4endl;
0537     nel++;
0538     sumMF += fCustomMatNamassfract;
0539   }
0540   if (fCustomMatPmassfract) {
0541     G4cout << "CustomMatPmassfract: " << fCustomMatPmassfract << G4endl;
0542     nel++;
0543     sumMF += fCustomMatPmassfract;
0544   }
0545   if (fCustomMatSmassfract) {
0546     G4cout << "CustomMatSmassfract: " << fCustomMatSmassfract << G4endl;
0547     nel++;
0548     sumMF += fCustomMatSmassfract;
0549   }
0550   if (fCustomMatClmassfract) {
0551     G4cout << "CustomMatClmassfract: " << fCustomMatClmassfract << G4endl;
0552     nel++;
0553     sumMF += fCustomMatClmassfract;
0554   }
0555   if (fCustomMatKmassfract) {
0556     G4cout << "CustomMatKmassfract: " << fCustomMatKmassfract << G4endl;
0557     nel++;
0558     sumMF += fCustomMatKmassfract;
0559   }
0560   if (fCustomMatCamassfract) {
0561     G4cout << "CustomMatCamassfract: " << fCustomMatCamassfract << G4endl;
0562     nel++;
0563     sumMF += fCustomMatCamassfract;
0564   }
0565   if (sumMF == 0.) {
0566     // set a default material (water),
0567     // otherwiswe an error appears in the interactive mode
0568     fCustomMatDensity = 1.00;  // g/cm3
0569     fCustomMatHmassfract = 0.1119;
0570     G4cout << "CustomMat set, but not used!" << G4endl;
0571     G4cout << "CustomMatHmassfract: " << fCustomMatHmassfract << G4endl;
0572     fCustomMatOmassfract = 0.8881;
0573     G4cout << "CustomMatOmassfract: " << fCustomMatOmassfract << G4endl;
0574     nel = 2;
0575     sumMF = 1;
0576   }
0577   if (sumMF != 1.) {
0578     G4String excep = "Error in Custom Material composition: check elemental mass fractions";
0579     G4Exception("DetectorConstuction::DefineMaterials()", "dc0004", FatalException, excep);
0580     return;
0581   }
0582 
0583   d = fCustomMatDensity * g / cm3;
0584   matname = "CustomMat";
0585   // Notice: this is an extended material
0586   fCustomMat = new G4ExtendedMaterial(matname, d, nel);
0587   if (fCustomMatHmassfract) fCustomMat->AddElement(elH, fCustomMatHmassfract);
0588   if (fCustomMatCmassfract) fCustomMat->AddElement(elC, fCustomMatCmassfract);
0589   if (fCustomMatNmassfract) fCustomMat->AddElement(elN, fCustomMatNmassfract);
0590   if (fCustomMatOmassfract) fCustomMat->AddElement(elO, fCustomMatOmassfract);
0591   if (fCustomMatNamassfract) fCustomMat->AddElement(elNa, fCustomMatNamassfract);
0592   if (fCustomMatPmassfract) fCustomMat->AddElement(elP, fCustomMatPmassfract);
0593   if (fCustomMatSmassfract) fCustomMat->AddElement(elS, fCustomMatSmassfract);
0594   if (fCustomMatClmassfract) fCustomMat->AddElement(elCl, fCustomMatClmassfract);
0595   if (fCustomMatKmassfract) fCustomMat->AddElement(elK, fCustomMatKmassfract);
0596   if (fCustomMatCamassfract) fCustomMat->AddElement(elCa, fCustomMatCamassfract);
0597   // Register MI extension
0598   fCustomMat->RegisterExtension(std::unique_ptr<G4MIData>(new G4MIData("MI")));
0599   G4MIData* dataMICustomMat = (G4MIData*)fCustomMat->RetrieveExtension("MI");
0600   dataMICustomMat->SetFilenameFF(fCustomMatFF);
0601 
0602   // Nist Materials
0603   fLead = NistMan->FindOrBuildMaterial("G4_Pb");
0604   fTungsten = NistMan->FindOrBuildMaterial("G4_W");
0605   fGe = NistMan->FindOrBuildMaterial("G4_Ge");
0606 }
0607 
0608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0609 
0610 void SAXSDetectorConstruction::SetGeometricalVariables()
0611 {
0612   // World
0613   fWorldSize = 10000. * mm;
0614 
0615   // Phantom
0616   fPhantomDiameter = 50. * mm;
0617   fPhantomHeight = 15. * mm;
0618   fPhantomZ = 0. * mm;
0619 
0620   // setup angle (rad)
0621   fthetaSetup = 0.;
0622 
0623   // Slits
0624   fSlitSize = 50. * mm;
0625   fSlit1Thickness = 5. * mm;
0626   fSlit2Thickness = 5. * mm;
0627   fSlit3Thickness = 5. * mm;
0628   fSlit4Thickness = 5. * mm;
0629   fSlit1SampleDistance = 350. * mm;
0630   fSlit2SampleDistance = 50. * mm;
0631   fSlit3SampleDistance = 50. * mm;
0632   fSlit4SampleDistance = 450. * mm;
0633   fSlit1xAperture = 4. * mm;
0634   fSlit2xAperture = 4. * mm;
0635   fSlit3xAperture = 4. * mm;
0636   fSlit4xAperture = 4. * mm;
0637   fSlit1yAperture = 4. * mm;
0638   fSlit2yAperture = 4. * mm;
0639   fSlit3yAperture = 4. * mm;
0640   fSlit4yAperture = 4. * mm;
0641 
0642   // Detector
0643   fDetectorSize = 20. * mm;
0644   fDetectorThickness = 20. * mm;
0645   fDetectorSampleDistance = 500. * mm;
0646 
0647   // Shielding
0648   fShieldingThickness = 4. * mm;
0649 }
0650 
0651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0652 
0653 G4VPhysicalVolume* SAXSDetectorConstruction::Construct()
0654 {
0655   // define custom materials
0656   DefineMaterials();
0657 
0658   // World
0659   fWorldSolid = new G4Box("WorldSolid", fWorldSize * 0.5, fWorldSize * 0.5, fWorldSize * 0.5);
0660 
0661   fWorldLogic = new G4LogicalVolume(fWorldSolid, fAir, "WorldLogic");
0662 
0663   fWorldLogic->SetVisAttributes(G4VisAttributes::GetInvisible());
0664 
0665   fWorldPhysical =
0666     new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), fWorldLogic, "WorldPhysical", 0, false, 0);
0667 
0668   // choose the phantom material
0669   switch (fPhantomMaterialIndex) {
0670     case (1):
0671       fPhantomMaterial = fWater;
0672       break;
0673     case (2):
0674       fPhantomMaterial = fMedMat;
0675       break;
0676     case (3):
0677       fPhantomMaterial = fPMMA;
0678       break;
0679     case (4):
0680       fPhantomMaterial = fAdipose;
0681       break;
0682     case (5):
0683       fPhantomMaterial = fGlandular;
0684       break;
0685     case (6):
0686       fPhantomMaterial = fBreast5050;
0687       break;
0688     case (7):
0689       fPhantomMaterial = fcarcinoma;
0690       break;
0691     case (8):
0692       fPhantomMaterial = fkidney;
0693       break;
0694     case (9):
0695       fPhantomMaterial = fliver;
0696       break;
0697     case (10):
0698       fPhantomMaterial = fFat;
0699       break;
0700     case (11):
0701       fPhantomMaterial = fBoneMatrix;
0702       break;
0703     case (12):
0704       fPhantomMaterial = fMineral;
0705       break;
0706     case (13):
0707       fPhantomMaterial = fBone;
0708       break;
0709     case (14):
0710       fPhantomMaterial = ffatLowX;
0711       break;
0712     case (15):
0713       fPhantomMaterial = fbonematrixLowX;
0714       break;
0715     case (16):
0716       fPhantomMaterial = fdryBoneLowX;
0717       break;
0718     case (17):
0719       fPhantomMaterial = fLexan;
0720       break;
0721     case (18):
0722       fPhantomMaterial = fKapton;
0723       break;
0724     case (19):
0725       fPhantomMaterial = fNylon;
0726       break;
0727     case (20):
0728       fPhantomMaterial = fPolyethylene;
0729       break;
0730     case (21):
0731       fPhantomMaterial = fPolystyrene;
0732       break;
0733     case (22):
0734       fPhantomMaterial = fFormaline;
0735       break;
0736     case (23):
0737       fPhantomMaterial = fAcetone;
0738       break;
0739     case (24):
0740       fPhantomMaterial = fHperoxide;
0741       break;
0742     case (25):
0743       fPhantomMaterial = fCIRS3070;
0744       break;
0745     case (26):
0746       fPhantomMaterial = fCIRS5050;
0747       break;
0748     case (27):
0749       fPhantomMaterial = fCIRS7030;
0750       break;
0751     case (28):
0752       fPhantomMaterial = fRMI454;
0753       break;
0754     case (29):
0755       fPhantomMaterial = fAir;
0756       break;
0757     case (30):
0758       fPhantomMaterial = fCustomMat;
0759       break;
0760   }
0761 
0762   // Phantom (cylinder with axis orthogonal to the X-ray beam axis)
0763   G4Tubs* PhantomSolid = new G4Tubs("PhantomSolid", 0., fPhantomDiameter * 0.5,
0764                                     fPhantomHeight * 0.5, 0. * deg, 360. * deg);
0765 
0766   fPhantomLogic = new G4LogicalVolume(PhantomSolid, fPhantomMaterial, "PhantomLogic");
0767 
0768   G4VisAttributes* PhantomVisAttribute = new G4VisAttributes(G4Colour(1., 1., 1.));
0769   PhantomVisAttribute->SetForceSolid(true);
0770   fPhantomLogic->SetVisAttributes(PhantomVisAttribute);
0771 
0772   G4cout << "Phantom material: " << fPhantomMaterial->GetName() << G4endl;
0773   G4cout << "Phantom density: " << fPhantomMaterial->GetDensity() / (g / cm3) << " g/cm3" << G4endl;
0774   G4cout << "Phantom mass: " << fPhantomLogic->GetMass() / g << " g" << G4endl;
0775 
0776   G4double rotAngle = 90. * CLHEP::deg;
0777   G4RotationMatrix* PhantomRotationMatrix = new G4RotationMatrix(0., 0., 0.);
0778   PhantomRotationMatrix->rotateX(rotAngle);
0779 
0780   fPhantomPhysical = new G4PVPlacement(PhantomRotationMatrix, G4ThreeVector(0., 0., fPhantomZ),
0781                                        fPhantomLogic, "PhantomPhysical", fWorldLogic, false, 0);
0782 
0783   // setup rotation matrix (downstream of the phantom/sample)
0784   G4RotationMatrix* SetupRotationMatrix = new G4RotationMatrix();
0785   SetupRotationMatrix->rotateY(-fthetaSetup);
0786 
0787   // Slits
0788   G4Box* Slit1OutSolid =
0789     new G4Box("Slit1OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit1Thickness * 0.5);
0790 
0791   G4Box* Slit2OutSolid =
0792     new G4Box("Slit2OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit2Thickness * 0.5);
0793 
0794   G4Box* Slit3OutSolid =
0795     new G4Box("Slit3OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit3Thickness * 0.5);
0796 
0797   G4Box* Slit4OutSolid =
0798     new G4Box("Slit4OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit4Thickness * 0.5);
0799 
0800   G4Box* Hole1Solid =
0801     new G4Box("Hole1Solid", fSlit1xAperture * 0.5, fSlit1yAperture * 0.5, fSlit1Thickness * 0.51);
0802 
0803   G4Box* Hole2Solid =
0804     new G4Box("Hole2Solid", fSlit2xAperture * 0.5, fSlit2yAperture * 0.5, fSlit2Thickness * 0.51);
0805 
0806   G4Box* Hole3Solid =
0807     new G4Box("Hole3Solid", fSlit3xAperture * 0.5, fSlit3yAperture * 0.5, fSlit3Thickness * 0.51);
0808 
0809   G4Box* Hole4Solid =
0810     new G4Box("Hole4Solid", fSlit4xAperture * 0.5, fSlit4yAperture * 0.5, fSlit4Thickness * 0.51);
0811 
0812   G4SubtractionSolid* Slit1Solid = new G4SubtractionSolid("Slit1Solid", Slit1OutSolid, Hole1Solid);
0813 
0814   G4SubtractionSolid* Slit2Solid = new G4SubtractionSolid("Slit1Solid", Slit2OutSolid, Hole2Solid);
0815 
0816   G4SubtractionSolid* Slit3Solid = new G4SubtractionSolid("Slit3Solid", Slit3OutSolid, Hole3Solid);
0817 
0818   G4SubtractionSolid* Slit4Solid = new G4SubtractionSolid("Slit4Solid", Slit4OutSolid, Hole4Solid);
0819 
0820   fSlit1Logic = new G4LogicalVolume(Slit1Solid, fTungsten, "Slit1Logic");
0821   fSlit2Logic = new G4LogicalVolume(Slit2Solid, fTungsten, "Slit2Logic");
0822   fSlit3Logic = new G4LogicalVolume(Slit3Solid, fTungsten, "Slit3Logic");
0823   fSlit4Logic = new G4LogicalVolume(Slit4Solid, fTungsten, "Slit4Logic");
0824 
0825   if (fIWantSlits) {
0826     G4cout << "Slit material: Tungsten" << G4endl;
0827     G4cout << "Slit1 thickness: " << fSlit1Thickness / mm << " mm" << G4endl;
0828     G4cout << "Slit2 thickness: " << fSlit2Thickness / mm << " mm" << G4endl;
0829     G4cout << "Slit3 thickness: " << fSlit3Thickness / mm << " mm" << G4endl;
0830     G4cout << "Slit4 thickness: " << fSlit4Thickness / mm << " mm" << G4endl;
0831     G4cout << "Slit1 aperture: " << fSlit1xAperture / mm << " x " << fSlit1yAperture / mm << " mm2"
0832            << G4endl;
0833     G4cout << "Slit2 aperture: " << fSlit2xAperture / mm << " x " << fSlit2yAperture / mm << " mm2"
0834            << G4endl;
0835     G4cout << "Slit3 aperture: " << fSlit3xAperture / mm << " x " << fSlit3yAperture / mm << " mm2"
0836            << G4endl;
0837     G4cout << "Slit4 aperture: " << fSlit4xAperture / mm << " x " << fSlit4yAperture / mm << " mm2"
0838            << G4endl;
0839   }
0840 
0841   G4VisAttributes* SlitlVisAttribute = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));
0842   SlitlVisAttribute->SetForceSolid(true);
0843   fSlit1Logic->SetVisAttributes(SlitlVisAttribute);
0844   fSlit2Logic->SetVisAttributes(SlitlVisAttribute);
0845   fSlit3Logic->SetVisAttributes(SlitlVisAttribute);
0846   fSlit4Logic->SetVisAttributes(SlitlVisAttribute);
0847 
0848   G4double Slit1z = fPhantomZ - fSlit1SampleDistance;
0849   G4ThreeVector Slit1PositionVector = G4ThreeVector(0., 0., Slit1z);
0850 
0851   G4double Slit2z = fPhantomZ - fSlit2SampleDistance;
0852   G4ThreeVector Slit2PositionVector = G4ThreeVector(0., 0., Slit2z);
0853 
0854   G4double Slit3x = fSlit3SampleDistance * std::sin(fthetaSetup);
0855   G4double Slit3z = fPhantomZ + fSlit3SampleDistance * std::cos(fthetaSetup);
0856   G4ThreeVector Slit3PositionVector = G4ThreeVector(Slit3x, 0., Slit3z);
0857 
0858   G4double Slit4x = fSlit4SampleDistance * std::sin(fthetaSetup);
0859   G4double Slit4z = fPhantomZ + fSlit4SampleDistance * std::cos(fthetaSetup);
0860   G4ThreeVector Slit4PositionVector = G4ThreeVector(Slit4x, 0., Slit4z);
0861 
0862   if (fIWantSlits) {
0863     fSlit1Physical = new G4PVPlacement(0, Slit1PositionVector, fSlit1Logic, "Slit1Physical",
0864                                        fWorldLogic, false, 0);
0865 
0866     fSlit2Physical = new G4PVPlacement(0, Slit2PositionVector, fSlit2Logic, "Slit2Physical",
0867                                        fWorldLogic, false, 0);
0868 
0869     fSlit3Physical = new G4PVPlacement(SetupRotationMatrix, Slit3PositionVector, fSlit3Logic,
0870                                        "Slit3Physical", fWorldLogic, false, 0);
0871 
0872     fSlit4Physical = new G4PVPlacement(SetupRotationMatrix, Slit4PositionVector, fSlit4Logic,
0873                                        "Slit4Physical", fWorldLogic, false, 0);
0874   }
0875 
0876   // Detector (with shielding)
0877   G4Tubs* DetectorSolid = new G4Tubs("DetectorSolid", 0., fDetectorSize * 0.5,
0878                                      fDetectorThickness * 0.5, 0. * deg, 360. * deg);
0879 
0880   fDetectorLogic = new G4LogicalVolume(DetectorSolid, fGe, "DetectorLogic");
0881 
0882   G4VisAttributes* DetectorVisAttribute = new G4VisAttributes(G4Colour(0., 0.5, 0.));
0883   DetectorVisAttribute->SetForceSolid(true);
0884   fDetectorLogic->SetVisAttributes(DetectorVisAttribute);
0885 
0886   G4double Detx = fDetectorSampleDistance * std::sin(fthetaSetup);
0887   G4double Detz = fPhantomZ + fDetectorSampleDistance * std::cos(fthetaSetup);
0888   G4ThreeVector DetectorPositionVector = G4ThreeVector(Detx, 0., Detz);
0889 
0890   fDetectorPhysical = new G4PVPlacement(SetupRotationMatrix, DetectorPositionVector, fDetectorLogic,
0891                                         "DetectorPhysical", fWorldLogic, false, 0);
0892 
0893   // Shielding
0894   G4double ShieldingSize = fDetectorSize + 2 * fShieldingThickness;
0895 
0896   G4double margin = 2.;
0897   G4double ShieldingLength = fDetectorThickness + fShieldingThickness * margin;
0898 
0899   G4double ShieldingSampleDistance = fDetectorSampleDistance + fDetectorThickness * 0.5
0900                                      + fShieldingThickness - ShieldingLength * 0.5;
0901 
0902   G4Tubs* ShieldingSolid = new G4Tubs("ShieldingSolid", fDetectorSize * 0.5, ShieldingSize * 0.5,
0903                                       ShieldingLength * 0.5, 0. * deg, 360. * deg);
0904 
0905   fShieldingLogic = new G4LogicalVolume(ShieldingSolid, fLead, "ShieldingLogic");
0906 
0907   G4cout << "Shielding material: Lead" << G4endl;
0908   G4cout << "Shielding thickness: " << fShieldingThickness / mm << " mm" << G4endl;
0909 
0910   G4VisAttributes* ShieldingVisAttribute = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
0911   ShieldingVisAttribute->SetForceSolid(true);
0912   fShieldingLogic->SetVisAttributes(ShieldingVisAttribute);
0913 
0914   G4double Shieldx = ShieldingSampleDistance * std::sin(fthetaSetup);
0915   G4double Shieldz = fPhantomZ + ShieldingSampleDistance * std::cos(fthetaSetup);
0916   G4ThreeVector ShieldingPositionVector = G4ThreeVector(Shieldx, 0., Shieldz);
0917 
0918   G4double ShieldingBackSampleDistance =
0919     fDetectorSampleDistance + fDetectorThickness * 0.5 + fShieldingThickness * 0.5;
0920 
0921   G4Tubs* ShieldingBackSolid = new G4Tubs("ShieldingBackSolid", 0., fDetectorSize * 0.5,
0922                                           fShieldingThickness * 0.5, 0. * deg, 360. * deg);
0923 
0924   fShieldingBackLogic = new G4LogicalVolume(ShieldingBackSolid, fLead, "ShieldingBackLogic");
0925 
0926   fShieldingBackLogic->SetVisAttributes(ShieldingVisAttribute);
0927 
0928   G4double ShieldBackx = ShieldingBackSampleDistance * std::sin(fthetaSetup);
0929   G4double ShieldBackz = fPhantomZ + ShieldingBackSampleDistance * std::cos(fthetaSetup);
0930   G4ThreeVector ShieldingBackPositionVector = G4ThreeVector(ShieldBackx, 0., ShieldBackz);
0931 
0932   fShieldingPhysical =
0933     new G4PVPlacement(SetupRotationMatrix, ShieldingPositionVector, fShieldingLogic,
0934                       "ShieldingPhysical", fWorldLogic, false, 0);
0935 
0936   fShieldingBackPhysical =
0937     new G4PVPlacement(SetupRotationMatrix, ShieldingBackPositionVector, fShieldingBackLogic,
0938                       "ShieldingBackPhysical", fWorldLogic, false, 0);
0939 
0940   return fWorldPhysical;
0941 }
0942 
0943 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0944 
0945 void SAXSDetectorConstruction::ConstructSDandField()
0946 {
0947   // Sensitive Volume
0948   G4VSensitiveDetector* vDetector = new SAXSSensitiveDetector("det");
0949   G4SDManager::GetSDMpointer()->AddNewDetector(vDetector);
0950   fSensitiveVolume = fDetectorLogic;
0951   fSensitiveVolume->SetSensitiveDetector(vDetector);
0952 }
0953 
0954 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......