Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • clantz2/zdclg
  • shengy3/zdclg
  • dmaclea2/zdclg
3 results
Show changes
Commits on Source (5)
Hide whitespace changes
Inline Side-by-side
Showing
with 185 additions and 672 deletions
geometry.mac
View file @ 5d1a82b9
/lightGuide/worldVolume 0.25 0.25 0.25 m
/lightGuide/envelope 89.75 113. 165. mm
/lightGuide/model/PMTDiameter 53 mm
/lightGuide/model/nSegmentsX 1
/lightGuide/model/nSegmentsZ 1
#----- Set the light guide to be used and position it -----
#If no model is selected, the program will default to a
#simplified version of the 2018 testbeam light guide
/lightGuide/model/CADmodel ../models/run2.stl
/lightGuide/model/translate 200 -125 0 mm
/lightGuide/model/translatePMT 0 -74.5 0 mm
#/lightGuide/model/CADmodel ../models/run2.stl
#/lightGuide/model/translate 200 -125 0 mm
#/lightGuide/model/translatePMT 0 -74.5 0 mm
#/lightGuide/model/CADmodel ../models/2020WC.stl
#/lightGuide/model/rotate 0 0 90
......
include/ASCIIPrimaryGenerator.hh deleted 100644 → 0
View file @ 878512fb
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
/// \file ASCIIPrimaryGenerator.hh
/// \brief Definition of the ASCIIPrimaryGenerator class
//
//
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef ASCIIPrimaryGenerator_h
#define ASCIIPrimaryGenerator_h 1
#include "G4VPrimaryGenerator.hh"
#include "lgAnalysis.hh"
#include <fstream>
#include <vector>
class G4Event;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class ASCIIPrimaryGenerator : public G4VPrimaryGenerator
{
public:
ASCIIPrimaryGenerator();
~ASCIIPrimaryGenerator();
virtual void SetInputFile(G4String _name);
virtual void GetNextEvent( );
inline G4int GetnEvents(){return fnEvents;}
public:
virtual void GeneratePrimaryVertex(G4Event*);
private:
std::vector< G4ThreeVector >* fPositionVec;
std::vector< G4ThreeVector >* fMomentumVec;
std::vector< G4double >* fEnergyVec;
G4int fEventNo;
G4int fnEvents = 0;
std::ifstream fInputFile;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
include/DetectorConstruction.hh
View file @ 5d1a82b9
......@@ -55,39 +55,29 @@ class DetectorConstruction : public G4VUserDetectorConstruction
virtual G4VPhysicalVolume* Construct();
void BuildWorld ();
void BuildTrapezoidLG ();
void BuildPMT ();
void BuildDefaultLG ();
void PlaceGeometry ();
void SetWorldVolume (G4ThreeVector arg);
void SetEnvelope (G4ThreeVector arg);
void SetRotation (G4ThreeVector arg);
void SetTranslation (G4ThreeVector arg);
void SetPMTTranslation (G4ThreeVector arg);
void SetPMTDiameter (G4double arg);
void SetLGthickness (G4double arg);
void UseCADModel (G4String fileName);
void OutputToGDML (G4String name);
void SetNSegmentsX (G4int arg);
void SetNSegmentsZ (G4int arg);
void SetSurfaceModel (const G4OpticalSurfaceModel model);
void SetSurfaceFinish (const G4OpticalSurfaceFinish finish);
void SetSurfaceType (const G4SurfaceType type);
void SetSurfaceSigmaAlpha (G4double v);
void AddSurfaceMPV (const char* c, G4MaterialPropertyVector* mpv);
void AddGasMPV (const char* c, G4MaterialPropertyVector* mpv);
private:
G4ThreeVector* m_worldDim;
G4ThreeVector* m_LGenvelope;
G4ThreeVector* m_LGpos;
G4ThreeVector* m_pmtPos;
G4RotationMatrix* m_rotation;
G4double m_pmtDia;
G4double m_PMTthickness;
G4double m_thickness;
G4int m_nSegmentsX;
G4int m_nSegmentsZ;
G4bool m_ConstructionHasBeenDone;
G4bool m_UsingCADmodel;
......@@ -95,29 +85,32 @@ class DetectorConstruction : public G4VUserDetectorConstruction
G4LogicalVolume* m_logicWorld;
G4VPhysicalVolume* m_physWorld;
G4Box* m_solidHalfWorld;
G4LogicalVolume* m_logicHalfWorld;
G4VPhysicalVolume* m_physHalfWorld;
//The light guide will probably be built from boolean solids
G4LogicalVolume* m_logicLG;
G4VPhysicalVolume* m_physLG;
G4Trd* m_inner;
G4Trd* m_outter;
G4SubtractionSolid* m_LightGuide;
G4LogicalVolume* m_logicLightGuide;
std::vector< G4VPhysicalVolume* > m_physLightGuide;
G4Tubs* m_fiberCore;
G4LogicalVolume* m_logicCore;
std::vector< G4VPhysicalVolume* > m_physCore;
G4Tubs* m_solidPMT;
G4LogicalVolume* m_logicPMT;
std::vector< G4VPhysicalVolume* > m_physPMT;
G4Tubs* m_fiberCladding;
G4LogicalVolume* m_logicCladding;
std::vector< G4VPhysicalVolume* > m_physCladding;
std::vector< G4LogicalBorderSurface* > m_Surfvec;
G4Tubs* m_fiberBuffer;
G4LogicalVolume* m_logicBuffer;
std::vector< G4VPhysicalVolume* > m_physBuffer;
Materials* materials;
G4Material* m_filler;
G4MaterialPropertiesTable* m_GasMPT;
G4Tubs* m_solidPMT;
G4LogicalVolume* m_logicPMT;
G4VPhysicalVolume* m_physPMT;
G4GDMLParser m_Parser;
G4RunManager* m_runMan;
DetectorMessenger* m_DetectorMessenger;
G4LogicalBorderSurface* m_Surface;
Materials* materials;
G4RunManager* m_runMan;
DetectorMessenger* m_DetectorMessenger;
};
#endif /*DetectorConstruction_h*/
include/DetectorMessenger.hh
View file @ 5d1a82b9
......@@ -64,15 +64,10 @@ class DetectorMessenger: public G4UImessenger{
// the model
G4UIcmdWithAString* fModelCmd;
G4UIcmdWith3VectorAndUnit* fWorldVolumeCmd;
G4UIcmdWith3VectorAndUnit* fEnvelopeCmd;
G4UIcmdWith3VectorAndUnit* fModelRotationCmd;
G4UIcmdWith3VectorAndUnit* fModelTranslationCmd;
G4UIcmdWith3VectorAndUnit* fPMTTranslationCmd;
G4UIcmdWithADoubleAndUnit* fPMTDiameterCmd;
G4UIcmdWithADoubleAndUnit* fLGThicknessCmd;
G4UIcmdWithAString* fOutputModelCmd;
G4UIcmdWithAnInteger* fNsegmentsXCmd;
G4UIcmdWithAnInteger* fNsegmentsZCmd;
// the surface
G4UIcmdWithAString* fSurfaceModelCmd;
......@@ -81,8 +76,5 @@ class DetectorMessenger: public G4UImessenger{
G4UIcmdWithADouble* fSurfaceSigmaAlphaCmd;
G4UIcmdWithAString* fSurfaceMatPropVectorCmd;
// the gas
G4UIcmdWithAString* fGasPropVectorCmd;
};
#endif
include/PrimaryGeneratorAction.hh
View file @ 5d1a82b9
......@@ -33,7 +33,6 @@
#include "G4VUserPrimaryGeneratorAction.hh"
#include "PrimaryGeneratorMessenger.hh"
#include "ASCIIPrimaryGenerator.hh"
#include "G4GeneralParticleSource.hh"
......@@ -57,10 +56,8 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
private:
G4GeneralParticleSource* fParticleGun;
ASCIIPrimaryGenerator* fASCIIParticleGun;
PrimaryGeneratorMessenger* fMessenger;
G4int fnEvents;
G4bool fUseASCIIInput;
G4bool fUseRootInput;
TFile* fInputFile;
TTree* fInputTree;
......
src/ASCIIPrimaryGenerator.cc deleted 100644 → 0
View file @ 878512fb
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
/// \file ASCIIPrimaryGenerator.cc
/// \brief Implementation of the PrimaryGenerator1 class
//
//
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "ASCIIPrimaryGenerator.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4PrimaryParticle.hh"
#include "G4PrimaryVertex.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
/*
*
*/
ASCIIPrimaryGenerator::ASCIIPrimaryGenerator()
: G4VPrimaryGenerator(),fEventNo(0)
{
fPositionVec = new std::vector< G4ThreeVector >;
fMomentumVec = new std::vector< G4ThreeVector >;
fEnergyVec = new std::vector< G4double >;
}
/*
*
*/
ASCIIPrimaryGenerator::~ASCIIPrimaryGenerator()
{ }
/*
*
*/
void ASCIIPrimaryGenerator::GeneratePrimaryVertex(G4Event* event)
{
G4ParticleDefinition* particleDefinition=
G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
std::vector< G4PrimaryVertex* > vertexVec;
std::vector< G4PrimaryParticle* > particleVec;
fPositionVec->clear();
fMomentumVec->clear();
fEnergyVec->clear();
GetNextEvent();
// Pass the original event number through to the ouput
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->FillNtupleIColumn(7, fEventNo );
//For each photon read in, add a new particle and vertex to the event
int nPhotons = fPositionVec->size();
for(int i = 0; i < nPhotons; i++){
particleVec.push_back( new G4PrimaryParticle(particleDefinition) );
particleVec.back()->SetMomentumDirection( fMomentumVec->at(i) );
particleVec.back()->SetKineticEnergy( fEnergyVec->at(i)*eV );
vertexVec.push_back( new G4PrimaryVertex( fPositionVec->at(i) , 0*s ) );
vertexVec.back()->SetPrimary( particleVec.back() );
event->AddPrimaryVertex( vertexVec.back() );
}
}
/*
* Fill member vectors with position, momentum and energy of all
* photons in the next event stored in the input ascii file
*/
void ASCIIPrimaryGenerator::GetNextEvent(){
std::string str;
double x,y,z,px,py,pz,energy;
// The first line of each event should be "Event #"
// If it's not, there's an error in the file formatting
// or we've reached the end of the file
getline(fInputFile,str);
if(str[0] == 'E'){
fEventNo = std::atoi( str.substr( 6, str.size() ).c_str() );
}else{
G4cout << "Error reading input event. Searching for next event..." << G4endl;
while( !fInputFile.eof() ){
getline(fInputFile,str);
if(str[0] == 'E'){
fEventNo = std::atoi( str.substr( 6, str.size() ).c_str() );
break;
}
}
if( fInputFile.eof() ){
G4cout << "Reached end of input file" << G4endl;
fEventNo = -1;
return;
}
}
// Loop through the event and read all photons.
// Format is:
// (Vertex) V, x, y, z
// (Momentum) P, px, py, pz, energy
// Then push all of those to their respective member vectors
getline(fInputFile,str);
while( !fInputFile.eof() ){
//
//Get Vertex
if(str[0] != 'V'){
G4cout << "Error reading input vertex at position " << fInputFile.tellg() << G4endl;
return;
}
x = std::atof( str.substr(2,19).c_str() );
y = std::atof( str.substr(20,37).c_str() );
z = std::atof( str.substr(38,55).c_str() );
fPositionVec->push_back( G4ThreeVector(x*mm,y*mm,z*mm) );
//
//Get momentum and energy
getline(fInputFile,str);
if(str[0] != 'P'){
G4cout << "Error reading input momentum at position " << fInputFile.tellg() << G4endl;
return;
}
px = std::atof( str.substr(2,19).c_str() );
py = std::atof( str.substr(20,37).c_str() );
pz = std::atof( str.substr(38,55).c_str() );
fMomentumVec->push_back( G4ThreeVector(px,py,pz) );
energy = std::atof( str.substr(56,73).c_str() );
fEnergyVec->push_back( energy );
//
//The next line is either another Vertex, or the end of the event
getline(fInputFile,str);
if( str[0] == 'E' ) return;
}
}
/*
* Open the input file and get the total number of events from it
*/
void ASCIIPrimaryGenerator::SetInputFile(G4String _name){
std::string str;
fInputFile.open( _name.c_str() );
getline(fInputFile,str);
fnEvents = std::atoi( str.substr(15,str.size() ).c_str() );
}
src/DetectorConstruction.cc
View file @ 5d1a82b9
......@@ -67,18 +67,16 @@
DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),
m_worldDim( new G4ThreeVector(0.25*m,0.5*m,0.25*m) ),
m_LGenvelope( new G4ThreeVector(89.75*mm,113.*mm,165.*mm) ),
m_LGpos( new G4ThreeVector() ),
m_pmtPos( new G4ThreeVector() ),
m_rotation( new G4RotationMatrix() ),
m_pmtDia(65.*mm),
m_PMTthickness(0.25*mm),
m_thickness(1.0*mm),
m_nSegmentsX(1),
m_nSegmentsZ(1),
m_ConstructionHasBeenDone(false),
m_UsingCADmodel(false),
m_logicWorld(0),
m_physLG(0),
m_physPMT(0),
m_DetectorMessenger(nullptr)
{
materials = Materials::getInstance();
......@@ -88,15 +86,6 @@ DetectorConstruction::DetectorConstruction()
m_DetectorMessenger = new DetectorMessenger(this);
m_runMan = G4RunManager::GetRunManager();
//Set default values
G4double PhotonEnergy[2] = {0.5*eV,8.0*eV};
G4double RefractiveIndexAir[2] = {1.0,1.0};
m_filler = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
m_GasMPT = new G4MaterialPropertiesTable();
m_GasMPT->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexAir, 2);
m_filler->SetMaterialPropertiesTable(m_GasMPT);
}
......@@ -126,7 +115,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
}
if(!m_logicWorld) BuildWorld();
if(!m_logicLightGuide) BuildTrapezoidLG();
if(!m_physLG) BuildDefaultLG();
PlaceGeometry();
m_ConstructionHasBeenDone = true;
......@@ -166,80 +155,111 @@ void DetectorConstruction::BuildWorld(){
G4VisAttributes* boxVisAtt_world= new G4VisAttributes(G4VisAttributes::Invisible);
m_logicWorld ->SetVisAttributes(boxVisAtt_world);
//----------------- Define Half of the world -----------------//
// This volume will be filled with gas with a different refractive
// index than the rest of the world and will hold the light guide
// and PMT
m_solidHalfWorld =
new G4Box("HalfWorld", //name
m_worldDim->x(), //sizeX
m_worldDim->y()/2, //sizeY
m_worldDim->z()); //sizeZ
m_logicHalfWorld =
new G4LogicalVolume(m_solidHalfWorld, //solid
m_filler, //material
"HalfWorld"); //name
m_physHalfWorld =
new G4PVPlacement(0, //no rotation
G4ThreeVector(0,m_worldDim->y()/2,0), //at (0,0,0)
m_logicHalfWorld, //logical volume
"HalfWorld", //name
m_logicWorld, //mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking
G4VisAttributes* boxVisAtt_half_world = new G4VisAttributes(G4Colour(0.0,0.0,1.0,0.1)); // or G4Colour(0.0,0.0,1.0,0.1)
m_logicHalfWorld ->SetVisAttributes(boxVisAtt_half_world);
}
/*
* Make the trapezoidal light guide logical volume
* Make the default light guide
*/
void DetectorConstruction::BuildTrapezoidLG( ){
//Divide the light guide envelope (cumulative) by the number
//of light guides that will occupy it
G4double xSize = (m_LGenvelope->x()/2.)/m_nSegmentsX;
G4double zSize = (m_LGenvelope->z()/2.)/m_nSegmentsZ;
//Determine the size of the top of the trapezoid so the
//square opening will be inscribed in the window of the PMT
G4double PMTwindow = .707*m_pmtDia/2.;
m_rotation->rotateX(90*deg);
//Aluminum outter
m_outter =
new G4Trd("BasicLightGuide",
xSize,
PMTwindow + m_thickness,
zSize,
PMTwindow + m_thickness,
m_LGenvelope->y()/2.);
//Air inner
m_inner =
new G4Trd("AirVolume",
xSize - m_thickness,
PMTwindow,
zSize - m_thickness,
PMTwindow,
m_LGenvelope->y()/2. + 0.1*mm);
//Subtract material to hollow out the light guide
m_LightGuide =
new G4SubtractionSolid("BasicLightGuide",
m_outter,
m_inner);
m_logicLightGuide =
new G4LogicalVolume(m_LightGuide,
materials->Al,
"BasicLightGuide");
void DetectorConstruction::BuildDefaultLG( ){
// Build the light guide here up until the logical volume
// Maybe build the aluminum portion into a single CAD solid and add
// fibers and air here
UseCADModel("../models/run2.stl");
bool CHECK_OVERLAPS = false;
//Change these variables to what they really should be
double fiberLength = 200*mm;
double core_diam = 0.610*mm;
double cladding_diam = 0.650*mm;
double buffer_diam = 0.710*mm;
//Core
m_fiberCore = new G4Tubs("fiber_core", // name
0.0*mm, // inner radius
core_diam/2.0,// outter radius
fiberLength, // length
0.0*deg, // sweep starting angle
360.0*deg); // sweep stopping angle
m_logicCore = new G4LogicalVolume( m_fiberCore, // solid to be created from
materials->pQuartz, // material
"fiber_logicCore"); // name
//Cladding
m_fiberCladding = new G4Tubs("fiber_cladding", // name
core_diam/2.0, // inner radius
cladding_diam/2.0, // outter radius
fiberLength, // length
0.0*deg, // sweep starting angle
360.0*deg); // sweep stopping angle
m_logicCladding = new G4LogicalVolume( m_fiberCladding, // solid to be created from
materials->pQuartz, // material
"fiber_logicCladding"); // name
//Buffer
m_fiberBuffer = new G4Tubs("fiber_buffer", // name
cladding_diam/2.0,// inner radius
buffer_diam/2.0, // outter radius
fiberLength, // length
0.0*deg, // sweep starting angle
360.0*deg); // sweep stopping angle
m_logicBuffer = new G4LogicalVolume( m_fiberBuffer, // solid to be created from
materials->pQuartz, // material
"fiber_logicBuffer"); // name
// Place the fibers
G4RotationMatrix* fiberRotation = new G4RotationMatrix();
fiberRotation->rotateX(90.*deg);
char name[32];
for(int i = 0; i < 16; i++){
// Find a good way to calculate the positions. I'm assuming only x varies here
double x = -20*mm + 0.8*i*mm;
//Core
sprintf(name,"fiberCore_%d", i);
m_physCore.push_back(
new G4PVPlacement(fiberRotation, //Rotation
G4ThreeVector( x , 0.0 , 0.0 ), //Position
m_logicCore, //Logical volume this is built from
name, //Name
m_logicWorld, //Mother volume
false, //If there are boolean operations
i, //Copy number
CHECK_OVERLAPS) ); //Check for overlaps with other geometry
//Cladding
sprintf(name,"fiberCladding_%d", i);
m_physCladding.push_back(
new G4PVPlacement(fiberRotation, //Rotation
G4ThreeVector( x , 0.0 , 0.0 ), //Position
m_logicCladding, //Logical volume this is built from
name, //Name
m_logicWorld, //Mother volume
false, //If there are boolean operations
i, //Copy number
CHECK_OVERLAPS) ); //Check for overlaps with other geometry
//Core
sprintf(name,"fiberBuffer_%d", i);
m_physBuffer.push_back(
new G4PVPlacement(fiberRotation, //Rotation
G4ThreeVector( x , 0.0 , 0.0 ), //Position
m_logicBuffer, //Logical volume this is built from
name, //Name
m_logicWorld, //Mother volume
false, //If there are boolean operations
i, //Copy number
CHECK_OVERLAPS) ); //Check for overlaps with other geometry
}
}
......@@ -262,7 +282,7 @@ void DetectorConstruction::BuildPMT(){
m_logicPMT =
new G4LogicalVolume(m_solidPMT, //solid
m_filler, //material
materials->Air, //material
"PMT"); //name
G4VisAttributes* VisAtt_PMT = new G4VisAttributes(G4Colour(1.0,1.0,0.6,0.7));
......@@ -272,75 +292,47 @@ void DetectorConstruction::BuildPMT(){
}
/*
*
* Place the geometry. This is done in it's own step in case
* the user has specified a CAD model to be used
*/
void DetectorConstruction::PlaceGeometry(){
G4double xPos = 0.;
G4double yPos = m_LGenvelope->y()/2. + m_LGpos->y() - m_worldDim->y()/2.;
G4double zPos = 0.;
G4double PMTx = 0.;
G4double PMTy = m_LGenvelope->y() + m_pmtPos->y() + m_PMTthickness - m_worldDim->y()/2.;
G4double PMTz = 0.;
G4RotationMatrix * PMTrot = new G4RotationMatrix();
PMTrot->rotateX(90*deg);
char name[40];
for(G4int xIndex = 0; xIndex < m_nSegmentsX; xIndex++ ){
//Start from -x, add user specified offset, add 1/2 the width of one LG, then a full width for each
xPos = - m_LGenvelope->x()/2. + m_LGpos->x() + (0.5 + xIndex)*m_LGenvelope->x()/m_nSegmentsX;
//Centered on the light guide plus any user specified offset
PMTx = xPos + m_pmtPos->x();
for(G4int zIndex = 0; zIndex < m_nSegmentsZ; zIndex++ ){
//Start from -z, add user specified offset, add 1/2 the length of one LG, then a full length for each
zPos = - m_LGenvelope->z()/2. + m_LGpos->z() + (0.5 + zIndex)*m_LGenvelope->z()/m_nSegmentsZ;
//Centered on the light guide plus any user specified offset
PMTz = zPos + m_pmtPos->z();
//----------------- Place the LightGuide -----------------//
sprintf(name,"LightGuide%d_%d",xIndex,zIndex);
m_physLightGuide.push_back(
new G4PVPlacement(m_rotation,
G4ThreeVector(xPos, yPos, zPos),
m_logicLightGuide,
name,
m_logicHalfWorld,
false,
0) );
//----------------- Place the PMT -----------------//
sprintf(name,"PMT%d_%d",xIndex,zIndex);
m_physPMT.push_back(
new G4PVPlacement(PMTrot,
G4ThreeVector(PMTx, PMTy, PMTz),
m_logicPMT,
name,
m_logicHalfWorld,
false,
0) );
//----------------- Define Optical Borders -----------------//
sprintf(name,"AlSurface%d_%d",xIndex,zIndex);
m_Surfvec.push_back(
new G4LogicalBorderSurface(name,
m_physLightGuide.back(),
m_physHalfWorld,
materials->AlSurface ) );
m_Surfvec.push_back(
new G4LogicalBorderSurface(name,
m_physHalfWorld,
m_physLightGuide.back(),
materials->AlSurface ) );
}//end y loop
}//end x loop
//----------------- Place the LightGuide -----------------//
m_physLG = new G4PVPlacement(m_rotation,
*m_LGpos,
m_logicLG,
"LightGuide",
m_logicWorld,
false,
0);
//----------------- Place the PMT -----------------//
m_physPMT = new G4PVPlacement(new G4RotationMatrix(),
*m_pmtPos,
m_logicPMT,
"PMT",
m_logicWorld,
false,
0);
//----------------- Define Optical Borders -----------------//
// One of these is correct. The order of the logical volumes is different
// between the two. This is specific to the direction of light i.e. from
// volume a to volume b or vice versa
// m_Surface = new G4LogicalBorderSurface("AlSurface",
// m_physLightGuide,
// m_physWorld,
// materials->AlSurface ) );
//
m_Surface = new G4LogicalBorderSurface("AlSurface",
m_physWorld,
m_physLG,
materials->AlSurface );
}
/*
......@@ -361,22 +353,6 @@ void DetectorConstruction::SetWorldVolume(G4ThreeVector arg){
}
}
/*
* Sets the envelope the light guide is constructed in
* If the geometry has already been constructed, reconstruct
* it with the new envelope
*/
void DetectorConstruction::SetEnvelope(G4ThreeVector arg){
if(m_LGenvelope) delete m_LGenvelope;
m_LGenvelope = new G4ThreeVector(arg);
if(m_ConstructionHasBeenDone){
Construct();
m_runMan->GeometryHasBeenModified();
G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
}
}
/*
* If the world already exists, rotate the light guide
* otherwise set the rotation for the light guide that will be built
......@@ -390,9 +366,7 @@ void DetectorConstruction::SetRotation(G4ThreeVector arg){
m_rotation->rotateZ(arg.z());
if(m_ConstructionHasBeenDone){
for(uint i = 0; i < m_physLightGuide.size(); i++){
m_physLightGuide[i]->SetRotation(m_rotation);
}
m_physLG->SetRotation(m_rotation);
m_runMan->GeometryHasBeenModified();
G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
}
......@@ -405,12 +379,7 @@ void DetectorConstruction::SetRotation(G4ThreeVector arg){
*/
void DetectorConstruction::SetTranslation(G4ThreeVector arg){
if(m_ConstructionHasBeenDone){
if(m_physLightGuide.size() == 1){
m_physLightGuide.back()->SetTranslation(arg);
} else {
m_LGpos = new G4ThreeVector(arg);
Construct();
}
m_physLG->SetTranslation(arg);
m_runMan->GeometryHasBeenModified();
G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
} else {
......@@ -426,12 +395,7 @@ void DetectorConstruction::SetTranslation(G4ThreeVector arg){
*/
void DetectorConstruction::SetPMTTranslation(G4ThreeVector arg){
if(m_ConstructionHasBeenDone){
if(m_physPMT.size() == 1){
m_physPMT.back()->SetTranslation(arg);
} else {
m_LGpos = new G4ThreeVector(arg);
Construct();
}
m_physPMT->SetTranslation(arg);
m_runMan->GeometryHasBeenModified();
G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
} else {
......@@ -453,40 +417,20 @@ void DetectorConstruction::SetPMTDiameter(G4double arg){
}
}
/*
* Set the thickness of the trapezoidal light guide shell
* If the geometry has already been constructed, reconstruct
* it with the new thickness
*/
void DetectorConstruction::SetLGthickness(G4double arg){
m_thickness = arg;
if(m_ConstructionHasBeenDone){
Construct();
m_runMan->GeometryHasBeenModified();
G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
}
}
/*
* Builds a light guide logical volume from CAD model file
*/
void DetectorConstruction::UseCADModel(G4String fileName){
//Delete anything we are going to redefine
if(m_logicLightGuide) delete m_logicLightGuide;
float xSize=0.0, ySize=0.0, zSize=0.0;
for(uint i = 0; i < m_physLightGuide.size(); i++){
delete m_physLightGuide[i];
delete m_physPMT[i];
delete m_Surfvec[2*i];
delete m_Surfvec[2*i+1];
if(m_ConstructionHasBeenDone){
delete m_logicLG;
delete m_physLG;
delete m_physPMT;
delete m_Surface;
}
m_physLightGuide.clear();
m_physPMT.clear();
m_Surfvec.clear();
float xSize=0.0, ySize=0.0, zSize=0.0;
G4String fileType = fileName.substr( fileName.last('.') + 1, fileName.size() - fileName.last('.'));
......@@ -497,24 +441,20 @@ void DetectorConstruction::UseCADModel(G4String fileName){
mesh->SetOffset( G4ThreeVector(-20*cm, 0, 0) );
mesh->SetReverse(false);
m_logicLightGuide =
m_logicLG =
new G4LogicalVolume(mesh->GetSolid(), //solid
materials->Al, //material
"LightGuide"); //name
}else{
G4cout << "File type not supported" << G4endl;
}
if(fileType == "gdml"){
m_Parser.Read(fileName);
}
if(fileType == "step"){
m_logicLightGuide = m_Parser.ParseST(fileName,materials->Air,materials->Al);
}
if(m_logicLightGuide !=0 ){
if(m_logicLG !=0 ){
//Print the extent to the console to help the user position the light guide
char message[64];
G4VisExtent extent = m_logicLightGuide->GetSolid()->GetExtent();
G4VisExtent extent = m_logicLG->GetSolid()->GetExtent();
xSize = extent.GetXmax() - extent.GetXmin();
ySize = extent.GetYmax() - extent.GetYmin();
zSize = extent.GetZmax() - extent.GetZmin();
......@@ -529,7 +469,6 @@ void DetectorConstruction::UseCADModel(G4String fileName){
}
if(m_pmtPos == 0) m_pmtPos = new G4ThreeVector( );
SetEnvelope( G4ThreeVector( m_nSegmentsX*xSize, ySize, m_nSegmentsZ*zSize ) );
if(m_ConstructionHasBeenDone){
PlaceGeometry();
......@@ -540,44 +479,6 @@ void DetectorConstruction::UseCADModel(G4String fileName){
}
/*
* Output the model of the light guide to GDML
*/
void DetectorConstruction::OutputToGDML(G4String fileName){
if(m_physLightGuide[0] != 0){
//For some reason this only works as a pointer
G4GDMLParser* gdml = new G4GDMLParser();
gdml->Write(fileName.c_str(),m_physLightGuide.back());
delete gdml;
}else{
G4cout << "No physical light guide defined..." << G4endl;
G4cout << "Can't write model to GDML" << G4endl;
}
}
/*
* Set the segmentation of the light guide envelope in x
* Reconstruct the geometry if necessary
*/
void DetectorConstruction::SetNSegmentsX(G4int arg){
m_nSegmentsX = arg;
if(m_ConstructionHasBeenDone){
Construct();
}
}
/*
* Set the segmentation of the light guide envelope in Z
* Reconstruct the geometry if necessary
*/
void DetectorConstruction::SetNSegmentsZ(G4int arg){
m_nSegmentsZ = arg;
if(m_ConstructionHasBeenDone){
Construct();
}
}
/*
* Set the surface model for the light guide
*/
......@@ -633,15 +534,3 @@ void DetectorConstruction::AddSurfaceMPV(const char* c, G4MaterialPropertyVector
materials->GetMPTArray().at(1)->DumpTable();
G4cout << "............." << G4endl;
}
/*
* Set the material property table of the volume
* which contains the light guide
*/
void DetectorConstruction::AddGasMPV(const char* c, G4MaterialPropertyVector* mpv){
mpv->SetSpline(true);
m_GasMPT->AddProperty(c, mpv);
G4cout << "The MPT for the gas is now: " << G4endl;
m_GasMPT->DumpTable();
G4cout << "............." << G4endl;
}
src/DetectorMessenger.cc
View file @ 5d1a82b9
......@@ -69,15 +69,6 @@ DetectorMessenger::DetectorMessenger(DetectorConstruction * Det)
fWorldVolumeCmd->SetDefaultValue(G4ThreeVector(0.25,0.25,0.5));
fWorldVolumeCmd->SetDefaultUnit("m");
fEnvelopeCmd =
new G4UIcmdWith3VectorAndUnit("/lightGuide/envelope", this);
fEnvelopeCmd->SetGuidance("Set the size light guide envelope");
fEnvelopeCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fEnvelopeCmd->SetToBeBroadcasted(false);
fEnvelopeCmd->SetParameterName("X","Y","Z",true);
fEnvelopeCmd->SetDefaultValue(G4ThreeVector(89.75*CLHEP::mm,113.*CLHEP::mm,339.*CLHEP::mm));
fEnvelopeCmd->SetDefaultUnit("m");
//Model commands
fModelCmd = new G4UIcmdWithAString("/lightGuide/model/CADmodel", this);
fModelCmd->SetGuidance("CAD model to be used");
......@@ -121,36 +112,6 @@ DetectorMessenger::DetectorMessenger(DetectorConstruction * Det)
fPMTDiameterCmd->SetDefaultValue( 65.0*CLHEP::mm );
fPMTDiameterCmd->SetDefaultUnit("mm");
fLGThicknessCmd =
new G4UIcmdWithADoubleAndUnit("/lightGuide/model/SkinThickness",this);
fLGThicknessCmd->SetGuidance("Set thickness of the light guide skin");
fLGThicknessCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fLGThicknessCmd->SetToBeBroadcasted(false);
fLGThicknessCmd->SetParameterName("thickness",true);
fLGThicknessCmd->SetDefaultValue( 1.0*CLHEP::mm );
fLGThicknessCmd->SetDefaultUnit("mm");
fOutputModelCmd = new G4UIcmdWithAString("/lightGuide/model/OutputModel",this);
fOutputModelCmd->SetGuidance("Creates a .gdml file of the light guide with the given name");
fOutputModelCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fOutputModelCmd->SetToBeBroadcasted(false);
fNsegmentsXCmd =
new G4UIcmdWithAnInteger("/lightGuide/model/nSegmentsX",this);
fNsegmentsXCmd->SetGuidance("Set segmentation of the light guide envelope in x");
fNsegmentsXCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fNsegmentsXCmd->SetToBeBroadcasted(false);
fNsegmentsXCmd->SetParameterName("nSegmentsX",true);
fNsegmentsXCmd->SetDefaultValue( 1 );
fNsegmentsZCmd =
new G4UIcmdWithAnInteger("/lightGuide/model/nSegmentsZ",this);
fNsegmentsZCmd->SetGuidance("Set segmentation of the light guide envelope in z");
fNsegmentsZCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fNsegmentsZCmd->SetToBeBroadcasted(false);
fNsegmentsZCmd->SetParameterName("nSegmentsZ",true);
fNsegmentsZCmd->SetDefaultValue( 1 );
//Surface commands
fSurfaceTypeCmd = new G4UIcmdWithAString("/lightGuide/surface/Type", this);
fSurfaceTypeCmd->SetGuidance("Surface type.");
......@@ -182,13 +143,6 @@ DetectorMessenger::DetectorMessenger(DetectorConstruction * Det)
fSurfaceMatPropVectorCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fSurfaceMatPropVectorCmd->SetToBeBroadcasted(false);
//Gas commands
fGasPropVectorCmd =
new G4UIcmdWithAString("/lightGuide/gasProperty", this);
fGasPropVectorCmd->SetGuidance("Set fill gas property vector");
fGasPropVectorCmd->AvailableForStates(G4State_PreInit, G4State_Init, G4State_Idle);
fGasPropVectorCmd->SetToBeBroadcasted(false);
}
/*
......@@ -197,21 +151,15 @@ DetectorMessenger::DetectorMessenger(DetectorConstruction * Det)
DetectorMessenger::~DetectorMessenger(){
delete fModelCmd;
delete fWorldVolumeCmd;
delete fEnvelopeCmd;
delete fModelRotationCmd;
delete fModelTranslationCmd;
delete fPMTTranslationCmd;
delete fPMTDiameterCmd;
delete fLGThicknessCmd;
delete fOutputModelCmd;
delete fNsegmentsXCmd;
delete fNsegmentsZCmd;
delete fSurfaceModelCmd;
delete fSurfaceFinishCmd;
delete fSurfaceTypeCmd;
delete fSurfaceSigmaAlphaCmd;
delete fSurfaceMatPropVectorCmd;
delete fGasPropVectorCmd;
}
/*
......@@ -227,10 +175,6 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
else if(command == fWorldVolumeCmd){
fDetector->SetWorldVolume(fWorldVolumeCmd->GetNew3VectorValue(newValue));
}
// LIGHT GUIDE ENVELOPE
else if(command == fEnvelopeCmd){
fDetector->SetEnvelope(fEnvelopeCmd->GetNew3VectorValue(newValue));
}
// MODEL ROTATION
else if(command == fModelRotationCmd){
fDetector->SetRotation(fModelRotationCmd->GetNew3VectorValue(newValue));
......@@ -247,22 +191,6 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
else if(command == fPMTDiameterCmd){
fDetector->SetPMTDiameter(fPMTDiameterCmd->GetNewDoubleValue(newValue));
}
// PMT DIAMETER
else if(command == fLGThicknessCmd){
fDetector->SetLGthickness(fLGThicknessCmd->GetNewDoubleValue(newValue));
}
// OUTPUT GEOMETRY TO GDML
else if(command == fOutputModelCmd){
fDetector->OutputToGDML(newValue);
}
// X SEGMENTATION
else if(command == fNsegmentsXCmd){
fDetector->SetNSegmentsX(fNsegmentsXCmd->GetNewIntValue(newValue));
}
// Z SEGMENTATION
else if(command == fNsegmentsZCmd){
fDetector->SetNSegmentsZ(fNsegmentsZCmd->GetNewIntValue(newValue));
}
// FINISH
else if (command == fSurfaceFinishCmd) {
if (newValue == "polished") {
......@@ -461,27 +389,4 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
fDetector->SetSurfaceSigmaAlpha(
G4UIcmdWithADouble::GetNewDoubleValue(newValue));
}
// Gas property
else if (command == fGasPropVectorCmd) {
// Convert string to physics vector
// string format is property name, then pairs of energy, value
// space delimited
G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector();
G4cout << newValue << G4endl;
std::istringstream instring(newValue);
G4String prop;
instring >> prop;
while (instring) {
G4String tmp;
instring >> tmp;
if (tmp == "") { break; }
G4double en = G4UIcommand::ConvertToDouble(tmp);
instring >> tmp;
G4double val;
val = G4UIcommand::ConvertToDouble(tmp);
mpv->InsertValues(en, val);
}
const char* c = prop.c_str();
fDetector->AddGasMPV(c, mpv);
}
}
src/PrimaryGeneratorAction.cc
View file @ 5d1a82b9
......@@ -43,10 +43,9 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction()
{
fParticleGun = new G4GeneralParticleSource();
fASCIIParticleGun = new ASCIIPrimaryGenerator();
fMessenger = new PrimaryGeneratorMessenger( this );
fUseRootInput = fUseASCIIInput = false;
fUseRootInput = false;
x = z = px = py = pz = Energy = time = 0;
evNo = 0;
}
......@@ -57,7 +56,6 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
delete fParticleGun;
delete fASCIIParticleGun;
delete fMessenger;
}
......@@ -66,9 +64,7 @@ PrimaryGeneratorAction::~PrimaryGeneratorAction()
*/
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
if(fUseASCIIInput){
fASCIIParticleGun->GeneratePrimaryVertex(anEvent);
}else if(fUseRootInput){
if(fUseRootInput){
GeneratePrimariesFromRootFile(anEvent);
}else{
fParticleGun->GeneratePrimaryVertex(anEvent);
......@@ -101,7 +97,7 @@ void PrimaryGeneratorAction::GeneratePrimariesFromRootFile(G4Event* anEvent){
}
void PrimaryGeneratorAction::SetInputFile(G4String _name){
if(fUseRootInput || fUseASCIIInput){
if(fUseRootInput){
G4cerr << "WARNING: Input file defined twice. Using first definition." << G4endl;
return;
}
......@@ -131,11 +127,5 @@ void PrimaryGeneratorAction::SetInputFile(G4String _name){
fInputTree->SetBranchAddress("Energy",&Energy);
fInputTree->SetBranchAddress("Time",&time);
}else if(check.contains(".txt")){
fUseASCIIInput = true;
G4cout << "Using " << _name << " as source" << G4endl;
fASCIIParticleGun->SetInputFile(_name);
fnEvents = fASCIIParticleGun->GetnEvents();
}
}