Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:21:14

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