#include #include #define THICKNESS 0.001 // define the thickness of air and aluminum in each forward calorimeter #include "LGDetectorConstruction.hh" #include "G4Element.hh" #include "G4ElementTable.hh" #include "G4Material.hh" #include "G4MaterialTable.hh" #include "G4NistManager.hh" #include "G4Box.hh" #include "G4Tubs.hh" #include "G4Trd.hh" #include "G4Cons.hh" #include "G4Torus.hh" #include "G4LogicalVolume.hh" #include "G4ThreeVector.hh" #include "G4RotationMatrix.hh" #include "G4Transform3D.hh" #include "G4PVPlacement.hh" #include "G4LogicalBorderSurface.hh" #include "G4OpBoundaryProcess.hh" #include "G4SubtractionSolid.hh" #include "G4ios.hh" LGDetectorConstruction::LGDetectorConstruction(G4double zLG, G4double BendR, G4double BendA): zTapperedSection(zLG), TorusBendRadius(BendR), TorusBendAngle(BendA) { //--- 1. LightGuide: You must be carefule because all the length means half of the actual length. x_LGrec = 2.54*cm/2.; //thicknes light guide rectangular cross section (paddle thickness) y_LGrec = 6.0*cm/2.; //hight light guide rectangular cross section (paddle width) z_LGrec = 3.0*cm/2.; //length or rectangular part x_Scin = x_LGrec; // part of scintillator to be used to generate light distribution at light guide entrance y_Scin = y_LGrec; z_Scin = 30.*cm; r_LGcyl = 2.5*cm; //radius of cylindrical part of light guide z_LGcyl = 7.*cm/2.; //length of cylindrical part of light guide x1_LGtrap = x_LGrec; y1_LGtrap = y_LGrec; x2_LGtrap = r_LGcyl; y2_LGtrap = x2_LGtrap; z_LGtrap = zTapperedSection*cm/2.; r1m_LGcone = sqrt(x1_LGtrap*x1_LGtrap + y1_LGtrap*y1_LGtrap); r1M_LGcone = y_LGrec+2.*cm; r2m_LGcone = r_LGcyl; r2M_LGcone = r_LGcyl+2.*cm; z_LGcone = z_LGtrap; // -------- 4. World x_World = 200*cm; y_World = 200*cm; z_World = 400*cm; } LGDetectorConstruction::~LGDetectorConstruction() {;} G4VPhysicalVolume* LGDetectorConstruction::Construct() { // ---------------------------------------------------------------------- // -------- 1. Elements and Materials // ---------------------------------------------------------------------- G4double a, z, density; // a: mass of a mole, z: mean number of protons G4int n_Elements; // n_Elements: the number of kinds of the elements G4int ncomponents, natoms; G4String name, symbol; G4NistManager* nist_Manager = G4NistManager::Instance(); // -------- 1.2. Materials for various parts of the detector G4Material* air = nist_Manager->FindOrBuildMaterial("G4_AIR"); // Air G4Material* Al = nist_Manager->FindOrBuildMaterial("G4_Al"); // Aluminum G4Material* PlexiGlass = nist_Manager->FindOrBuildMaterial("G4_PLEXIGLASS"); // Plexiglass //G4Material* pmtGlass = nist_Manager->FindOrBuildMaterial("G4_SILICON_DIOXIDE");// PMTWindow // for pmtGlass see lower // define Elements G4Element *elH, *elC; a = 1.01*g/mole; elH = new G4Element(name="Hydrogen",symbol="H" , z= 1., a); a = 12.01*g/mole; elC = new G4Element(name="Carbon" ,symbol="C" , z= 6., a); //G4Material* Hydrogen = nist_Manager->FindOrBuildMaterial("G4_H"); //G4Material* Carbon = nist_Manager->FindOrBuildMaterial("G4_C"); G4Material* Oxigen = nist_Manager->FindOrBuildMaterial("G4_O"); G4Material* Silicon = nist_Manager->FindOrBuildMaterial("G4_Si"); G4Material* Potassium = nist_Manager->FindOrBuildMaterial("G4_K"); G4Material* Sodium = nist_Manager->FindOrBuildMaterial("G4_Na"); //G4Material* Lead = nist_Manager->FindOrBuildMaterial("G4_Pb"); G4Material* Boron = nist_Manager->FindOrBuildMaterial("G4_B"); G4Material* Aluminum = nist_Manager->FindOrBuildMaterial("G4_Al"); // data for ELJEN 200 scintillator density = 1.023*g/cm3; G4Material* Scintillator = new G4Material(name="Scintillator", density, ncomponents=2); Scintillator->AddElement(elC, natoms=469); Scintillator->AddElement(elH, natoms=517); //Fraction by weight for Borosicilcat PMT window glass //source http://en.wikipedia.org/wiki/Borosilicate_glass //Element Atomicnumber Fraction //---------------------------------------- //B 5 0.040064 //O 8 0.539562 //Na 11 0.028191 //Al 13 0.011644 //Si 14 0.377220 //K 19 0.003321 G4Material* pmtGlass = new G4Material("Borosicilate",density=2.23*g/cm3,n_Elements=6); pmtGlass->AddMaterial(Oxigen, 0.539562); pmtGlass->AddMaterial(Boron, 0.040064); pmtGlass->AddMaterial(Sodium, 0.028191); pmtGlass->AddMaterial(Aluminum, 0.011644); pmtGlass->AddMaterial(Silicon, 0.377220); pmtGlass->AddMaterial(Potassium, 0.003321); // ---------------------------------------------------------------------- // -------- 2. Generate and Add Material Property Tables // ---------------------------------------------------------------------- const G4int num = 4; G4double eRange_Photon[num] = {0.5*eV,1.24*eV, 2.27*eV, 3.76*eV}; G4double absLength_PlexiGlass[num] = {500.*cm,500.*cm, 500.*cm, 500.*cm}; G4double absLength_PMTGlass[num] = {200.*cm,200.*cm, 200.*cm, 200.*cm}; G4double index_Air[num] = {1.,1., 1., 1.}; // Index of refraction of air G4double index_Scintillator[num] = {1.58,1.58, 1.58, 1.58}; G4double absLength_Scintillator[num] = {400.*cm,400.*cm, 400.*cm, 400.*cm}; // source http://en.wikipedia.org/wiki/Borosilicate_glass //G4double index_PMTGlass[num] = {1.52,1.52, 1.52, 1.52}; G4double index_PMTGlass[num] = {1.48,1.48, 1.48, 1.48}; // in original code // Index of refraction of the PlexiGlass light guide G4double index_PlexiGlass[num] = {1.51,1.51, 1.51, 1.51}; G4double eff_Cathode[num] = {1.,1., 1., 1.}; // eff = efficiency G4double ref_Cathode[num] = {0.,0., 0., 0.}; // ref = reflectivity // -------- 2.1. Table - Air MPT= Material Property Table G4MaterialPropertiesTable* air_MPT = new G4MaterialPropertiesTable(); air_MPT->AddProperty("RINDEX", eRange_Photon, index_Air, num); air->SetMaterialPropertiesTable(air_MPT); G4MaterialPropertiesTable* Scintillator_MPT = new G4MaterialPropertiesTable(); Scintillator_MPT->AddProperty("RINDEX", eRange_Photon, index_Scintillator, num); Scintillator_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_Scintillator, num); Scintillator->SetMaterialPropertiesTable(Scintillator_MPT); // -------- 2.4. Table - plexiglass G4MaterialPropertiesTable* PlexiGlass_MPT = new G4MaterialPropertiesTable(); PlexiGlass_MPT->AddProperty("RINDEX", eRange_Photon, index_PlexiGlass, num); PlexiGlass_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_PlexiGlass, num); PlexiGlass->SetMaterialPropertiesTable(PlexiGlass_MPT); // -------- 2.5. Table - PMTglass G4MaterialPropertiesTable* pmtGlass_MPT = new G4MaterialPropertiesTable(); pmtGlass_MPT->AddProperty("RINDEX", eRange_Photon, index_PMTGlass, num); pmtGlass_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_PMTGlass, num); pmtGlass->SetMaterialPropertiesTable(pmtGlass_MPT); // -------- 2.6. Table - PMTcathode G4MaterialPropertiesTable* cathode_MPT = new G4MaterialPropertiesTable(); cathode_MPT->AddProperty("EFFICIENCY", eRange_Photon, eff_Cathode, num); cathode_MPT->AddProperty("REFLECTIVITY", eRange_Photon, ref_Cathode, num); Al->SetMaterialPropertiesTable(cathode_MPT); // ---------------------------------------------------------------------- // -------- 3. Volumes // ---------------------------------------------------------------------- // -------- 3.1. World Volume // ---------------------------------------------------------------------- // Definiton for the world volume G4Box* world_box = new G4Box("World", x_World, y_World, z_World); // Definiton of the logical volume of the world volume filled with the air G4LogicalVolume* world_log = new G4LogicalVolume(world_box, air, "World"); // Placement of the world volume G4VPhysicalVolume* world_phys = new G4PVPlacement(0, G4ThreeVector(), world_log, "World", 0, false, 0); // ---------------------------------------------------------------------- // -------- 3.2. Scintillator with light guide // ---------------------------------------------------------------------- G4Box* LGBox_box = new G4Box("LGBox_box", x_LGrec, y_LGrec, z_LGrec); G4LogicalVolume* LGBox_log = new G4LogicalVolume(LGBox_box, PlexiGlass, "LGBox_log"); G4VPhysicalVolume* LGBox_phys = NULL; G4double zpos = 0; LGBox_phys = new G4PVPlacement(0, G4ThreeVector(0., 0., zpos), LGBox_log, "LGBox", world_log, false, 0); G4Box* ScinBox_box = new G4Box("ScinBox_box", x_Scin, y_Scin, z_Scin); G4LogicalVolume* ScinBox_log = new G4LogicalVolume(ScinBox_box, PlexiGlass,"ScinBox_log"); G4VPhysicalVolume* ScinBox_phys = NULL; ScinBox_phys = new G4PVPlacement(0, G4ThreeVector(0., 0., -z_Scin-z_LGrec), ScinBox_log, "ScinBox", world_log, false, 0); G4Trd* LGTrd_trd = new G4Trd("LGTrd", x1_LGtrap, x2_LGtrap, y1_LGtrap, y2_LGtrap, z_LGtrap); G4Cons* LGCon_con = new G4Cons("LGCon", r1m_LGcone, r1M_LGcone, r2m_LGcone, r2M_LGcone, z_LGcone, 0.*deg, 360.*deg); G4VSolid* TapperedSection = new G4SubtractionSolid("TapperedSection", LGTrd_trd, LGCon_con); G4LogicalVolume* LGTappered_log = new G4LogicalVolume(TapperedSection, PlexiGlass, "LGTappered_log"); G4double TorusBend = TorusBendAngle*deg; G4double TorusRad = TorusBendRadius*cm; zpos += (z_LGrec + z_LGcone); G4VPhysicalVolume* LGTappered_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,zpos), LGTappered_log, "LGTappered", world_log, false, 0); G4Torus* TorusSection_tor = new G4Torus("TorusSection", 0., r_LGcyl, TorusRad, 0., TorusBend); G4LogicalVolume* TorusSection_log = new G4LogicalVolume(TorusSection_tor, PlexiGlass, "TorusSection_log"); G4RotationMatrix* RotX = new G4RotationMatrix(); RotX->rotateX(-90.*deg); zpos += z_LGtrap; G4VPhysicalVolume*TorusSection_phys= new G4PVPlacement(RotX, G4ThreeVector(-TorusRad,0.,zpos), TorusSection_log, "TorusSection", world_log, false, 0); G4Tubs* LGTub_tubs = new G4Tubs("LGTub_tubs", 0., r_LGcyl, z_LGcyl, 0.*deg, 360.*deg); G4LogicalVolume* LGTub_log = new G4LogicalVolume(LGTub_tubs, PlexiGlass, "LGTub_log"); //zpos += (z_LGtrap + z_LGcyl); G4RotationMatrix* RotY = new G4RotationMatrix(); RotY->rotateY(TorusBend); G4double transz1 = TorusRad*sin(TorusBend); G4double transz2 = z_LGcyl*sin(M_PI/2.-TorusBend); G4double transx1 = TorusRad - TorusRad*cos(TorusBend); G4double transx2 = z_LGcyl*cos(M_PI/2.-TorusBend); G4double dx = transx1 + transx2; G4double dz = transz1 + transz2; zpos += dz; G4VPhysicalVolume* LGTub_phys = new G4PVPlacement(RotY, G4ThreeVector(-dx, 0., zpos), LGTub_log, "LGTub", world_log, false, 0); // pmtglass G4double PMTGlassThickness = 1.0*mm; G4Tubs* pmtGlass_tubs = new G4Tubs("PMTGlass_tubs", 0., r_LGcyl, PMTGlassThickness, 0.*deg, 360.*deg); G4LogicalVolume* pmtGlass_log = new G4LogicalVolume(pmtGlass_tubs, pmtGlass, "PMTGlass_log"); dx += transx2; transx1 = PMTGlassThickness * cos(M_PI/2.-TorusBend); dx += transx1; zpos += transz2; transz1 = PMTGlassThickness * sin(M_PI/2.-TorusBend); zpos += transz1; G4VPhysicalVolume* pmtGlass_phys = new G4PVPlacement(RotY, G4ThreeVector(-dx, 0, zpos), pmtGlass_log, "PMTGlass", world_log, false, 0); // pmt cathode G4double PMTCathodeThickness = 0.001*mm; G4Tubs* pmtCathode_tubs = new G4Tubs("PMTCathode_tubs", 0., r_LGcyl*0.99, PMTCathodeThickness, 0.*deg, 360.*deg); G4LogicalVolume* pmtCathode_log = new G4LogicalVolume(pmtCathode_tubs, Aluminum, "PMTCathode_log"); dx += transx1; transx2 = PMTCathodeThickness * cos(M_PI/2.-TorusBend); dx += transx2; zpos += transz1; transz2 = PMTCathodeThickness * sin(M_PI/2.-TorusBend); zpos += transz2; G4VPhysicalVolume* pmtCathode_phys = NULL; pmtCathode_phys = new G4PVPlacement(RotY, G4ThreeVector(-dx, 0, zpos), pmtCathode_log, "PMTCathode", world_log, false, 0); // ---------------------------------------------------------------------- // -------- 4. Surfaces // ---------------------------------------------------------------------- // -------- 4.1. Generate and Add Surface Property Tables // ---------------------------------------------------------------------- G4double ref_PlexiGlass[num] = {0.95,0.95, 0.95, 0.95}; G4double eff_PlexiGlass[num] = {0.05,0.05, 0.05, 0.05}; G4double specular_Lobe[num] = {0.045,0.045, 0.045, 0.045}; G4double specular_Spike[num] = {0.95,.95, 0.95, 0.95}; G4double back_Scatter[num] = {0.005,0.005, 0.005, 0.005}; // surface properties from Chris Cude simulation //G4double specular_Lobe[num] = {0.18,0.18, 0.18, 0.18}; //G4double specular_Spike[num] = {0.7954,.7954, 0.7954, 0.7954}; //G4double back_Scatter[num] = {0.0246,0.0246, 0.0246, 0.0246}; // specular_lobe + specular_spike + back_scatter + defuse_lobe = 1 // ---------------------------------------------------------------------- // -------- 4.2. Surfaces // ---------------------------------------------------------------------- // -------- 4.2.1. Surface of the leadglass in contact with air G4OpticalSurface* PlexiGlassAir_Surf = new G4OpticalSurface("PlexiGlassAirSurface"); PlexiGlassAir_Surf->SetType(dielectric_dielectric); //PlexiGlassAir_Surf->SetFinish(polished); PlexiGlassAir_Surf->SetFinish(ground); PlexiGlassAir_Surf->SetModel(unified); // becaus of option "ground" use a sigma for the angle alpha needs // to be set. A sigma of 10 would indicate a very rough surface. PlexiGlassAir_Surf->SetSigmaAlpha(0.1); // ---------------------------------------------------------------------- // -------- 4.3. Property Tables // ---------------------------------------------------------------------- // plexiglass to air interface G4LogicalBorderSurface* plexiGlassAir_Surf_log1 = NULL; G4LogicalBorderSurface* plexiGlassAir_Surf_log2 = NULL; G4LogicalBorderSurface* plexiGlassAir_Surf_log3 = NULL; G4LogicalBorderSurface* plexiGlassAir_Surf_log4 = NULL; plexiGlassAir_Surf_log1 = new G4LogicalBorderSurface("PlexiGlassAirSurface1", LGBox_phys, world_phys, PlexiGlassAir_Surf); plexiGlassAir_Surf_log2 = new G4LogicalBorderSurface("PlexiGlassAirSurface2", LGTappered_phys, world_phys, PlexiGlassAir_Surf); plexiGlassAir_Surf_log3 = new G4LogicalBorderSurface("PlexiGlassAirSurface3", TorusSection_phys, world_phys, PlexiGlassAir_Surf); plexiGlassAir_Surf_log4 = new G4LogicalBorderSurface("PlexiGlassAirSurface4", LGTub_phys, world_phys, PlexiGlassAir_Surf); // -------- 4.1.4. LeadGlass(LeadGlass_Air) G4MaterialPropertiesTable* PlexiGlass_Air_MPT = new G4MaterialPropertiesTable(); PlexiGlass_Air_MPT->AddProperty("SPECULARLOBECONSTANT", eRange_Photon, specular_Lobe, num); PlexiGlass_Air_MPT->AddProperty("SPECULARSPIKECONSTANT", eRange_Photon, specular_Spike, num); PlexiGlass_Air_MPT->AddProperty("BACKSCATTERCONSTANT", eRange_Photon, back_Scatter, num); PlexiGlassAir_Surf->SetMaterialPropertiesTable(PlexiGlass_Air_MPT); return world_phys; // return the physical volume of the world volume }