401 lines
15 KiB
C++
401 lines
15 KiB
C++
|
|
#include "CloverDetectorConstruction.hh"
|
|
#include "CloverCrystalSD.hh"
|
|
#include "Detector_Parameters.hh"
|
|
|
|
#include "G4Material.hh"
|
|
#include "G4NistManager.hh"
|
|
|
|
#include "G4Box.hh"
|
|
#include "G4Tubs.hh"
|
|
#include "G4Sphere.hh"
|
|
#include "G4Trd.hh"
|
|
#include "G4UnionSolid.hh"
|
|
#include "G4IntersectionSolid.hh"
|
|
#include "G4SubtractionSolid.hh"
|
|
|
|
#include "G4LogicalVolume.hh"
|
|
#include "G4PVPlacement.hh"
|
|
#include "G4PVReplica.hh"
|
|
#include "G4GlobalMagFieldMessenger.hh"
|
|
#include "G4AutoDelete.hh"
|
|
|
|
#include "G4RotationMatrix.hh"
|
|
|
|
#include "G4SDManager.hh"
|
|
|
|
#include "G4VisAttributes.hh"
|
|
#include "G4Colour.hh"
|
|
|
|
#include "G4PhysicalConstants.hh"
|
|
#include "G4SystemOfUnits.hh"
|
|
|
|
#include <vector>
|
|
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4ThreadLocal
|
|
G4GlobalMagFieldMessenger* CloverDetectorConstruction::fMagFieldMessenger = 0;
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
CloverDetectorConstruction::CloverDetectorConstruction()
|
|
: G4VUserDetectorConstruction(),
|
|
fNumOfCrystal(0),
|
|
fLogicCrystal(NULL), fCrystalMaterial(NULL),
|
|
fCheckOverlaps(true)
|
|
{
|
|
fNumOfCrystal = Crystal_Num; //also need to change the fNDet in CloverEventAction.cc
|
|
fLogicCrystal = new G4LogicalVolume* [fNumOfCrystal];
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
CloverDetectorConstruction::~CloverDetectorConstruction()
|
|
{
|
|
delete [] fLogicCrystal;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4VPhysicalVolume* CloverDetectorConstruction::Construct()
|
|
{
|
|
// Define materials
|
|
DefineMaterials();
|
|
|
|
// Define volumes
|
|
return DefineVolumes();
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void CloverDetectorConstruction::DefineMaterials()
|
|
{
|
|
// Lead material defined using NIST Manager
|
|
auto nistManager = G4NistManager::Instance();
|
|
|
|
// Vacuum
|
|
new G4Material("Galactic", 1., 1.01*g/mole, universe_mean_density, kStateGas, 2.73*kelvin, 3.e-18*pascal);
|
|
|
|
// Air
|
|
nistManager->FindOrBuildMaterial("G4_AIR");
|
|
|
|
// Iron
|
|
nistManager->FindOrBuildMaterial("G4_Fe");
|
|
|
|
// Aluminium
|
|
nistManager->FindOrBuildMaterial("G4_Al");
|
|
|
|
// Ge
|
|
fCrystalMaterial = nistManager->FindOrBuildMaterial("G4_Ge");
|
|
|
|
// Print materials
|
|
G4cout << *(G4Material::GetMaterialTable()) << G4endl;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
|
|
{
|
|
// Geometry parameters
|
|
G4double crystalLength = 8.*cm;
|
|
G4double crystalRadius = 25.*mm;
|
|
G4double crystalZPos = Radius*mm + crystalLength/2.; //distance to target + distance to center of crystals.;
|
|
G4double cutXY = 46.0 * mm;
|
|
|
|
G4double worldSizeZ = 2*(crystalLength + crystalZPos);
|
|
G4double worldSizeXY = worldSizeZ;
|
|
/*
|
|
//pipe
|
|
G4double pipeOuterRadius = 100/2.* mm;
|
|
G4double pipeWallThickness = 5 * mm;
|
|
G4double pipeLength = 70* mm;
|
|
G4double pipeZPos = 0.*cm;
|
|
*/
|
|
|
|
//spherical chamber
|
|
G4double chamberOuterRadius = 150/2.* mm;
|
|
G4double chamberWallThickness = 5 * mm;
|
|
G4double chamberInnerRadius = chamberOuterRadius - chamberWallThickness;
|
|
G4double chamberZPos = 0.*cm;
|
|
|
|
|
|
//Clover casing
|
|
G4double caseXYInner = 2*(cutXY + 1.0) * mm;
|
|
G4double caseZInner = crystalLength/2 + crystalLength + 10*mm;
|
|
G4double caseXYWallThickness = 3.0 * mm;
|
|
G4double caseZWallThickness = 3.0 * mm;
|
|
|
|
|
|
//BGO Detector
|
|
|
|
G4double bgoXYInner = 10*mm + 2*(cutXY + 1.0) * mm;
|
|
G4double bgoZInner = crystalLength/2 + 4 * mm + crystalLength + 15*mm;
|
|
G4double bgoXYWallThickness = 103.5 * mm; //from schematic
|
|
G4double bgoanglecut0 = 53 * mm ;
|
|
G4double zbgo0 = 178 * mm;
|
|
|
|
G4double bgoanglecut1 = 24.935 * mm;
|
|
G4double zbgo1 = 80*mm;
|
|
|
|
G4double bgoanglecut2 = 10 *mm ;
|
|
G4double zbgo2 = 20*mm;
|
|
|
|
// Get materials
|
|
auto backgroundMaterial = G4Material::GetMaterial("G4_AIR");
|
|
//auto backgroundMaterial = G4Material::GetMaterial("Galactic");
|
|
|
|
// World
|
|
auto worldS = new G4Box("World", worldSizeXY/2, worldSizeXY/2, worldSizeZ);
|
|
|
|
auto worldLV = new G4LogicalVolume( worldS, backgroundMaterial, "World");
|
|
|
|
auto worldPV = new G4PVPlacement( 0, // no rotation
|
|
G4ThreeVector(), // at (0,0,0)
|
|
worldLV, // its logical volume
|
|
"World", // its name
|
|
0, // its mother volume
|
|
false, // no boolean operation
|
|
0, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
worldLV->SetVisAttributes (G4VisAttributes::GetInvisible());
|
|
|
|
// vaccum Pipe
|
|
|
|
// G4Tubs * pipe = new G4Tubs("pipe", pipeOuterRadius - pipeWallThickness , pipeOuterRadius, pipeLength, 0, 360*degree);
|
|
// G4Tubs * pipe0 = new G4Tubs("pipe0", 0 , pipeOuterRadius, pipeLength, 0, 360*degree);
|
|
|
|
G4Sphere *chamber = new G4Sphere("chamber",chamberInnerRadius,chamberOuterRadius,0,2*pi,0,pi);
|
|
|
|
G4RotationMatrix * rotX = new G4RotationMatrix(); rotX->rotateX(90*degree);
|
|
G4RotationMatrix * rotY = new G4RotationMatrix(); rotY->rotateY(90*degree);
|
|
|
|
G4LogicalVolume * chamberLV = new G4LogicalVolume(chamber, G4Material::GetMaterial("G4_Fe"),"targetchamber");
|
|
|
|
new G4PVPlacement( 0, // no rotation
|
|
G4ThreeVector(0, 0, chamberZPos), // at (0,0,0)
|
|
chamberLV, // its logical volume
|
|
"targetchamber", // its name
|
|
worldLV, // its mother volume
|
|
false, // no boolean operation
|
|
0, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
chamberLV->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,1.0)));
|
|
|
|
|
|
// Crystals
|
|
G4VisAttributes * crystalVisAtt= new G4VisAttributes(G4Colour(0.5,0.5,1.0));
|
|
crystalVisAtt->SetVisibility(true);
|
|
|
|
//define crystal shapes
|
|
|
|
G4Tubs * base = new G4Tubs("base", 0, crystalRadius, crystalLength/2., 0, 360*degree);
|
|
G4Box * cut = new G4Box("cut", cutXY/2. , cutXY/2. , crystalLength/2. * 1.2);
|
|
G4IntersectionSolid * crystalS = new G4IntersectionSolid("HPGe", base, cut);
|
|
|
|
G4Box* case1 = new G4Box("case1", caseXYInner/2., caseXYInner/2., caseZInner );
|
|
G4Box* case2 = new G4Box("case2", (caseXYInner + caseXYWallThickness)/2. , (caseXYInner + caseXYWallThickness)/2., caseZInner + caseZWallThickness*2.);
|
|
G4SubtractionSolid * casing = new G4SubtractionSolid("casing", case2, case1);
|
|
|
|
G4Box* bgo1 = new G4Box("bgo1", bgoXYInner/2., bgoXYInner/2., bgoZInner +zbgo1);
|
|
|
|
G4Trd* bgo2 = new G4Trd("bgo2", (bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut1 ,(bgoXYInner+ bgoXYWallThickness)/2.
|
|
, (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut1, (bgoXYInner + bgoXYWallThickness)/2., zbgo0/2);
|
|
|
|
G4Trd* bgo3 = new G4Trd("bgo3", (bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 + bgoanglecut2,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut1
|
|
, (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0 + bgoanglecut2, (bgoXYInner + bgoXYWallThickness)/2.- bgoanglecut1, zbgo1/2);
|
|
G4ThreeVector bgo3_v = G4ThreeVector(0.,0.,-zbgo0 + 50*mm);
|
|
|
|
G4RotationMatrix * bgo3_m = new G4RotationMatrix();
|
|
|
|
|
|
G4VSolid* bgo_u = new G4UnionSolid("bgo2 + bgo3",bgo2,bgo3,bgo3_m,bgo3_v);
|
|
//G4Box* bgo2 = new G4Box("bgo2", (bgoXYInner+ bgoXYWallThickness)/2., (bgoXYInner + bgoXYWallThickness)/2., bgoZInner + bgoZWallThickness*2.);
|
|
|
|
|
|
G4Trd* bgowindow = new G4Trd("bgowindow", (bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 ,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 +bgoanglecut2
|
|
, (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0, (bgoXYInner + bgoXYWallThickness)/2. -bgoanglecut0 +bgoanglecut2, zbgo2/2);
|
|
|
|
|
|
|
|
G4SubtractionSolid * bgo = new G4SubtractionSolid("bgo", bgo_u, bgo1);
|
|
|
|
// creating the subtraction solid that will form the actual bgo sheild.
|
|
|
|
|
|
|
|
double dTheta,dPsi,zmod;
|
|
|
|
for( G4int i = 0 ; i < fNumOfCrystal ; i++){
|
|
G4String name = "HPGe"+ std::to_string(i);
|
|
G4cout << " crystal name : " << name << G4endl;
|
|
G4double phi = 360./4. * degree;
|
|
|
|
G4double rho = 0;
|
|
rho = (cutXY/2 + 0.3* mm) / sin(phi/2.);
|
|
|
|
int iClover = i/4; //How to get each group of crystals.
|
|
int jCrystal = i%4; //How to renumber each group as 0,1,2,3.
|
|
|
|
fLogicCrystal[i] = new G4LogicalVolume(crystalS, fCrystalMaterial, "CrystalLV");
|
|
|
|
|
|
auto caseLV = new G4LogicalVolume * [fNumOfCrystal/4];
|
|
auto bgoLV = new G4LogicalVolume * [fNumOfCrystal/4];
|
|
auto bgowindowLV = new G4LogicalVolume * [fNumOfCrystal/4];
|
|
|
|
fLogicCrystal[i]->SetVisAttributes(crystalVisAtt);
|
|
|
|
|
|
if(int j = iClover){
|
|
|
|
dTheta = Clover_Location[j][1];
|
|
|
|
dPsi = Clover_Location[j][2];
|
|
|
|
zmod = Clover_Location[j][3];
|
|
|
|
}
|
|
|
|
|
|
G4ThreeVector pos = G4ThreeVector(rho * cos(jCrystal* phi + phi/2.), rho * sin(jCrystal * phi + phi/2.), crystalZPos+zmod);
|
|
|
|
// if(iClover != DetNum){return 0;}
|
|
|
|
|
|
G4ThreeVector case_pos = G4ThreeVector(0,0,88*mm+crystalZPos); //The crystalZpos is the radius of the detector,
|
|
//bgo uses an identical rotation matrix as the case,
|
|
//you can also add in zmod if you want to shift the shielding.
|
|
G4ThreeVector bgo_pos = G4ThreeVector(0,0,115*mm+crystalZPos);
|
|
G4ThreeVector bgowindow_pos = G4ThreeVector(0,0,-63*mm + crystalZPos);
|
|
|
|
//Set rotation matrix settings for clovers.
|
|
|
|
G4RotationMatrix * bgowindowmat = new G4RotationMatrix();
|
|
G4RotationMatrix * bgomat = new G4RotationMatrix();
|
|
G4RotationMatrix * rotmat = new G4RotationMatrix();
|
|
G4RotationMatrix * crystalmat= new G4RotationMatrix();
|
|
|
|
//x = sin(theta)cos(phi) || y = sin(theta)sin(phi) || z = cos(theta)
|
|
|
|
|
|
for(int j = 0; j <= iClover;j++){
|
|
if(j == iClover){
|
|
|
|
// rotation matrix must undo the operations you perform with the position vector.
|
|
bgowindowmat->rotateZ(-dPsi*degree);
|
|
bgowindowmat->rotateX(-dTheta*degree);
|
|
|
|
bgomat->rotateZ(-dPsi*degree);
|
|
bgomat->rotateX(-dTheta*degree);
|
|
|
|
rotmat->rotateZ(-dPsi*degree);
|
|
rotmat->rotateX(-dTheta*degree);
|
|
|
|
crystalmat->rotateZ(-dPsi*degree);
|
|
crystalmat->rotateX(-dTheta*degree);
|
|
|
|
//rotate position vector based on iClover.
|
|
|
|
pos.rotateX(dTheta*degree);
|
|
pos.rotateZ(dPsi *degree);
|
|
|
|
case_pos.rotateX(dTheta *degree);
|
|
case_pos.rotateZ(dPsi*degree);
|
|
|
|
bgo_pos.rotateX(dTheta *degree);
|
|
bgo_pos.rotateZ(dPsi*degree);
|
|
|
|
bgowindow_pos.rotateX(dTheta *degree);
|
|
bgowindow_pos.rotateZ(dPsi*degree);
|
|
}
|
|
}
|
|
new G4PVPlacement( crystalmat, // orientation matrix
|
|
pos, // its position vector
|
|
fLogicCrystal[i], // its logical volume
|
|
name, // its name
|
|
worldLV, // its mother volume
|
|
false, // no boolean operation
|
|
i, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
|
|
if( jCrystal == 0 ) {
|
|
caseLV[iClover] = new G4LogicalVolume(casing, G4Material::GetMaterial("G4_Al"), "Case");
|
|
bgoLV[iClover] = new G4LogicalVolume(bgo , G4Material::GetMaterial("G4_Al"), "BGO ");
|
|
bgowindowLV[iClover] = new G4LogicalVolume(bgowindow, G4Material::GetMaterial("G4_Al"), "BGO WINDOW ");
|
|
|
|
new G4PVPlacement( rotmat , // case orientation
|
|
case_pos, // case location
|
|
caseLV[i/4], // its logical volume
|
|
"Casing", // its name
|
|
worldLV, // its mother volume
|
|
false, // no boolean operation
|
|
iClover, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
new G4PVPlacement( bgomat , // bgo orientation
|
|
bgo_pos, // bgo location
|
|
bgoLV[i/4], // its logical volume
|
|
"Bgo", // its name
|
|
worldLV, // its mother volume
|
|
false, // no boolean operation
|
|
iClover, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
new G4PVPlacement( bgowindowmat , // bgo orientation
|
|
bgowindow_pos, // bgo location
|
|
bgowindowLV[i/4], // its logical volume
|
|
"Bgowindow", // its name
|
|
worldLV, // its mother volume
|
|
false, // no boolean operation
|
|
iClover, // copy number
|
|
fCheckOverlaps); // checking overlaps
|
|
|
|
caseLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,0.0)));
|
|
bgoLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.,1.0,0.0)));
|
|
bgowindowLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.0,1.0,1.0)));
|
|
}
|
|
//Set a i+1%4 for each
|
|
//
|
|
|
|
}
|
|
|
|
// Al casing
|
|
|
|
|
|
|
|
// Always return the physical World
|
|
return worldPV;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void CloverDetectorConstruction::ConstructSDandField()
|
|
{
|
|
|
|
G4cout << "********************* CloverDetectorConstruction::ConstructSDandField() " << G4endl;
|
|
G4SDManager::GetSDMpointer()->SetVerboseLevel(1);
|
|
|
|
// Sensitive detectors
|
|
CloverCrystalSD * crystalSD = new CloverCrystalSD("Crystal", "CrystalHitsCollection", fNumOfCrystal);
|
|
G4SDManager::GetSDMpointer()->AddNewDetector(crystalSD);
|
|
|
|
//Set crystalSD to all logical volumn under the same name of "CrystalLV"
|
|
SetSensitiveDetector("CrystalLV",crystalSD, true);
|
|
|
|
// Create global magnetic field messenger.
|
|
// Uniform magnetic field is then created automatically if
|
|
// the field value is not zero.
|
|
G4ThreeVector fieldValue;
|
|
fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
|
|
fMagFieldMessenger->SetVerboseLevel(1);
|
|
|
|
// Register the field messenger for deleting
|
|
G4AutoDelete::Register(fMagFieldMessenger);
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|