#include "RE02DetectorConstruction.hh" #include "G4PSEnergyDeposit3D.hh" #include "G4PSNofStep3D.hh" #include "G4PSCellFlux3D.hh" #include "G4PSPassageCellFlux3D.hh" #include "G4PSFlatSurfaceFlux3D.hh" #include "G4PSFlatSurfaceCurrent3D.hh" #include "G4SDParticleWithEnergyFilter.hh" #include "G4SDParticleFilter.hh" #include "G4SDChargedFilter.hh" #include "G4NistManager.hh" #include "G4Material.hh" #include "G4Box.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4SDManager.hh" #include "G4PVParameterised.hh" #include "G4PhantomParameterisation.hh" #include "G4VisAttributes.hh" #include "G4Colour.hh" #include "G4SystemOfUnits.hh" #include "G4ios.hh" RE02DetectorConstruction::RE02DetectorConstruction() : G4VUserDetectorConstruction() { // Default size of water phantom,and segmentation. fPhantomSize.setX(200.*mm); fPhantomSize.setY(200.*mm); fPhantomSize.setZ(400.*mm); fNx = fNy = fNz = 100; fInsertLead = TRUE; } RE02DetectorConstruction::~RE02DetectorConstruction() {;} G4VPhysicalVolume* RE02DetectorConstruction::Construct() { //===================== // Material Definitions //===================== G4NistManager* NISTman = G4NistManager::Instance(); G4Material* air = NISTman->FindOrBuildMaterial("G4_AIR"); G4Material* water = NISTman->FindOrBuildMaterial("G4_WATER"); G4Material* lead = NISTman->FindOrBuildMaterial("G4_AIR"); // Print all the materials defined. G4cout << G4endl << "The materials defined are : " << G4endl << G4endl; G4cout << *(G4Material::GetMaterialTable()) << G4endl; std::vector phantomMat = { air , water }; //============================================================================ // Definitions of Solids, Logical Volumes, Physical Volumes //============================================================================ //------------------------------------------------------------------ // World Volume //------------------------------------------------------------------ G4ThreeVector worldSize = G4ThreeVector(200*cm, 200*cm, 200*cm); G4Box * solidWorld = new G4Box("world", worldSize.x()/2., worldSize.y()/2., worldSize.z()/2.); G4LogicalVolume * logicWorld = new G4LogicalVolume(solidWorld, air, "World", 0, 0, 0); // Must place the World Physical volume unrotated at (0,0,0). G4VPhysicalVolume * physiWorld = new G4PVPlacement(0, // no rotation G4ThreeVector(), // at (0,0,0) logicWorld, // its logical volume "World", // its name 0, // its mother volume false, // no boolean operations 0); // copy number //------------------------------------------------------------------ // Constructing the container for the phantom //------------------------------------------------------------------ G4ThreeVector phantomSize = fPhantomSize; G4Box * solidPhantom = new G4Box("phantom", phantomSize.x()/2., phantomSize.y()/2., phantomSize.z()/2.); G4LogicalVolume * logicPhantom = new G4LogicalVolume(solidPhantom, water, "Phantom", 0, 0, 0); G4RotationMatrix* rot = new G4RotationMatrix(); //rot->rotateY(30.*deg); G4ThreeVector positionPhantom; G4VPhysicalVolume * physiPhantom = new G4PVPlacement(rot, // no rotation positionPhantom, // at (x,y,z) logicPhantom, // its logical volume "Phantom", // its name logicWorld, // its mother volume false, // no boolean operations 0); // copy number // - Default number of segmentation is defined at constructor. G4int nxCells = fNx; G4int nyCells = fNy; G4int nzCells = fNz; // This is setting the pixel size G4ThreeVector sensSize; sensSize.setX(phantomSize.x()/(G4double)nxCells); sensSize.setY(phantomSize.y()/(G4double)nyCells); sensSize.setZ(phantomSize.z()/(G4double)nzCells); //------------------------------------------------------------------ // Constructing the Voxel and parameterisation //------------------------------------------------------------------ G4String zVoxName("phantomSens"); G4VSolid* solVoxel = new G4Box(zVoxName,sensSize.x()/2.,sensSize.y()/2.,sensSize.z()/2.); fLVPhantomSens = new G4LogicalVolume(solVoxel,water,zVoxName); // Parameterisation for transformation of voxels. G4PhantomParameterisation* params = new G4PhantomParameterisation(); params->SetVoxelDimensions(sensSize.x() / 2., sensSize.y() / 2., sensSize.z() / 2.); params->SetNoVoxels(fNx, fNy, fNz); params->SetMaterials(phantomMat); // Get fills container params->CheckVoxelsFillContainer(solidPhantom->GetXHalfLength(), solidPhantom->GetYHalfLength(), solidPhantom->GetZHalfLength()); // set the material IDs as just a patterened list in x and y size_t* materialIDs = new size_t[fNx * fNy * fNz]; G4int bx, by, bz, bw; for (G4int index = 0; index < fNx * fNy * fNz; index++) { // Get indexes bz = index / (fNy * fNx); // x in N(A) bw = index % (fNy * fNx); // w in N(BC) by = bw / (fNx); // y in N(B) bx = bw % (fNx); // z in N(C) // If conditions are true then change material to one, else set to zero if ((bx % 100 == 0) || (by % 50 == 0)) { materialIDs[index] = 1; } else { materialIDs[index] = 0; } } params->SetMaterialIndices(materialIDs); params->SetSkipEqualMaterials(1); params->BuildContainerSolid(physiPhantom); // Create parameterisation G4PVParameterised* physiPhantomSens = new G4PVParameterised("PhantomSens", // their name fLVPhantomSens, // their logical volume physiPhantom, // Mother logical volume kUndefined, // Are placed along this axis fNx * fNy * fNz, // Number of cells params); // Parameterisation. // Optimization flag is avaiable for, // kUndefined, kXAxis, kYAxis, kZAxis. // // physiPhantomSens-> physiPhantomSens->SetRegularStructureId(1); return physiWorld; } void RE02DetectorConstruction::ConstructSDandField() { //================================================ // Sensitive detectors : MultiFunctionalDetector //================================================ // Sensitive Detector Manager. G4SDManager* pSDman = G4SDManager::GetSDMpointer(); // Sensitive Detector Name G4String phantomSDname = "PhantomSD"; //------------------------ // MultiFunctionalDetector //------------------------ // Define MultiFunctionalDetector with name. G4MultiFunctionalDetector* mFDet = new G4MultiFunctionalDetector(phantomSDname); pSDman->AddNewDetector( mFDet ); // Register SD to SDManager. fLVPhantomSens->SetSensitiveDetector(mFDet); // Assign SD to the logical volume. //--------------------------------------- // SDFilter : Sensitive Detector Filters //--------------------------------------- // Particle Filter for Primitive Scorer with filter name(fltName) // and particle name(particleName), // or particle names are given by add("particle name"); method. G4String fltName,particleName; //-- proton filter G4SDParticleFilter* protonFilter = new G4SDParticleFilter(fltName="protonFilter", particleName="proton"); //-- electron filter G4SDParticleFilter* electronFilter = new G4SDParticleFilter(fltName="electronFilter"); electronFilter->add(particleName="e+"); // accept electrons. electronFilter->add(particleName="e-"); // accept positorons. //-- charged particle filter G4SDChargedFilter* chargedFilter = new G4SDChargedFilter(fltName="chargedFilter"); //------------------------ // PS : Primitive Scorers //------------------------ // Primitive Scorers are used with SDFilters according to your purpose. //-- Primitive Scorer for Energy Deposit. // Total, by protons, by electrons. G4String psName; G4PSEnergyDeposit3D * scorer0 = new G4PSEnergyDeposit3D(psName="totalEDep", fNx,fNy,fNz); G4PSEnergyDeposit3D * scorer1 = new G4PSEnergyDeposit3D(psName="protonEDep",fNx,fNy,fNz); scorer1->SetFilter(protonFilter); //-- Number of Steps for protons G4PSNofStep3D * scorer2 = new G4PSNofStep3D(psName="protonNStep",fNx,fNy,fNz); scorer2->SetFilter(protonFilter); //-- CellFlux for charged particles G4PSPassageCellFlux3D * scorer3 = new G4PSPassageCellFlux3D(psName="chargedPassCellFlux", fNx,fNy,fNz); G4PSCellFlux3D * scorer4 = new G4PSCellFlux3D(psName="chargedCellFlux", fNx,fNy,fNz); G4PSFlatSurfaceFlux3D * scorer5 = new G4PSFlatSurfaceFlux3D(psName="chargedSurfFlux", fFlux_InOut,fNx,fNy,fNz); scorer3->SetFilter(chargedFilter); scorer4->SetFilter(chargedFilter); scorer5->SetFilter(chargedFilter); //------------------------------------------------------------ // Register primitive scorers to MultiFunctionalDetector //------------------------------------------------------------ mFDet->RegisterPrimitive(scorer0); mFDet->RegisterPrimitive(scorer1); mFDet->RegisterPrimitive(scorer2); mFDet->RegisterPrimitive(scorer3); mFDet->RegisterPrimitive(scorer4); mFDet->RegisterPrimitive(scorer5); //======================== // More additional Primitive Scoreres //======================== //--- Surface Current for gamma with energy bin. // This example creates four primitive scorers. // 4 bins with energy --- Primitive Scorer Name // 1. to 10 KeV, gammaSurfCurr000 // 10 keV to 100 KeV, gammaSurfCurr001 // 100 keV to 1 MeV, gammaSurfCurr002 // 1 MeV to 10 MeV. gammaSurfCurr003 // for ( G4int i = 0; i < 4; i++){ std::ostringstream name; name << "gammaSurfCurr" << std::setfill('0') << std::setw(3) << i; G4String psgName = name.str(); G4double kmin = std::pow(10.,(G4double)i)*keV; G4double kmax = std::pow(10.,(G4double)(i+1))*keV; //-- Particle with kinetic energy filter. G4SDParticleWithEnergyFilter* pkinEFilter = new G4SDParticleWithEnergyFilter(fltName="gammaE filter",kmin,kmax); pkinEFilter->add("gamma"); // Accept only gamma. pkinEFilter->show(); // Show accepting condition to stdout. //-- Surface Current Scorer which scores number of tracks in unit area. G4PSFlatSurfaceCurrent3D * scorer = new G4PSFlatSurfaceCurrent3D(psgName,fCurrent_InOut,fNx,fNy,fNz); scorer->SetFilter(pkinEFilter); // Assign filter. mFDet->RegisterPrimitive(scorer); // Register it to MultiFunctionalDetector. } }