diff --git a/include/ASCIIPrimaryGenerator.hh b/include/ASCIIPrimaryGenerator.hh deleted file mode 100644 index 4bb8a9b28bc8914dc5d439dc7be70bd60bc27d73..0000000000000000000000000000000000000000 --- a/include/ASCIIPrimaryGenerator.hh +++ /dev/null @@ -1,71 +0,0 @@ -// -// ******************************************************************** -// * 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 diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh index 72d3fb564a278c1cdac2934f06bc38f5911d0e7b..c06e061b99766dbd90a52709669842eda5077aa3 100644 --- a/include/DetectorConstruction.hh +++ b/include/DetectorConstruction.hh @@ -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*/ diff --git a/include/DetectorMessenger.hh b/include/DetectorMessenger.hh index 1bbc683fe0547e6e9a1b82263243df315f293304..da10fd7db21037f75d04d4a700b9aacf72d4eb33 100644 --- a/include/DetectorMessenger.hh +++ b/include/DetectorMessenger.hh @@ -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 diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh index 0cc68fba564630fcf17098a552d53e80949f0bd7..c7e8dd93fd30fed1fb03448bc8b1a14f09bdb7dd 100644 --- a/include/PrimaryGeneratorAction.hh +++ b/include/PrimaryGeneratorAction.hh @@ -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; diff --git a/src/ASCIIPrimaryGenerator.cc b/src/ASCIIPrimaryGenerator.cc deleted file mode 100644 index c1adf2a126f584b7c1d30ec776ba61b5be9bdc17..0000000000000000000000000000000000000000 --- a/src/ASCIIPrimaryGenerator.cc +++ /dev/null @@ -1,178 +0,0 @@ -// -// ******************************************************************** -// * 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() ); -} diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index 595d0e453800a70bc938995e0fb7ea0aba007992..1d205fa5f8057ef7b9066f2b89dca5a1a91d4041 100644 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -67,15 +67,11 @@ 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), @@ -88,15 +84,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 +113,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){ } if(!m_logicWorld) BuildWorld(); - if(!m_logicLightGuide) BuildTrapezoidLG(); PlaceGeometry(); m_ConstructionHasBeenDone = true; @@ -166,80 +152,109 @@ 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 + + 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_logicCore"); // 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_logicCore"); // 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 +277,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 +287,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_logicLightGuide, + "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_physHalfWorld, + // materials->AlSurface ) ); + // + m_Surface = new G4LogicalBorderSurface("AlSurface", + m_physHalfWorld, + m_physLightGuide, + materials->AlSurface ) ); + } /* @@ -361,22 +348,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 @@ -405,7 +376,7 @@ void DetectorConstruction::SetRotation(G4ThreeVector arg){ */ void DetectorConstruction::SetTranslation(G4ThreeVector arg){ if(m_ConstructionHasBeenDone){ - if(m_physLightGuide.size() == 1){ + if(m_physLG.size() == 1){ m_physLightGuide.back()->SetTranslation(arg); } else { m_LGpos = new G4ThreeVector(arg); @@ -453,40 +424,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_logicLightGuide; + delete m_physLightGuide; + 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('.')); @@ -540,44 +491,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 +546,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; -} diff --git a/src/DetectorMessenger.cc b/src/DetectorMessenger.cc index f80d887d2418cbcc3fa0c260020c1104a1082304..fb703b58f4d3936c473c6fb56bc295e233c24587 100644 --- a/src/DetectorMessenger.cc +++ b/src/DetectorMessenger.cc @@ -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); - } } diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index 749f98a2e445f9656a8cc4b06e6ce1995466b65b..a72754d7f6f801a4776ac1d231d994dfa43c07c7 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -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; } @@ -131,11 +129,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(); } }