Skip to content

GitLab

Commit 224cb773 authored by shengy3's avatar shengy3
Browse files

merge from master

parents 58b32833 698a7799
Hide whitespace changes
Inline Side-by-side
MonteCarlo/.vscode/settings.json 0 → 100644
View file @ 224cb773
{
"files.associations": {
"analysis.C": "cpp"
}
}
\ No newline at end of file
MonteCarlo/include/ReactorTestFiber.hh 0 → 100644
View file @ 224cb773
//
// ********************************************************************
// * 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. *
// ********************************************************************
//
// Author: Sheng Yang
#ifndef ReactorTestFiber_h
#define ReactorTestFiber_h 1
#include "globals.hh"
#include "G4PVPlacement.hh"
#include "G4SubtractionSolid.hh"
#include "Materials.hh"
#include <vector>
class G4VPhysicalVolume;
class G4LogicalVolume;
class G4VSolid;
class G4Material;
class SharedData;
/// Detector construction class to define materials and geometry.
class ReactorTestFiber
{
public:
ReactorTestFiber(const int,
const G4ThreeVector&,
G4LogicalVolume*, SharedData*);
ReactorTestFiber();
~ReactorTestFiber();
virtual void Construct();
virtual void DefineMaterials();
// virtual void DefineBorderProperties();
virtual void ConstructDetector();
protected:
const int m_modNum;
const G4ThreeVector m_pos;
G4LogicalVolume* m_logicMother;
SharedData* m_sd;
Materials* materials;
protected:
G4Material* m_matQuartz;
G4Material* m_Al;
G4Material* m_Poly;
G4Material* m_Air;
G4Material* m_PMMA;
G4Material* m_Grease;
// G4LogicalBorderSurface* quartz_grease_LBS[64];
// G4LogicalBorderSurface* quartz_grease_LBS[64];
// G4LogicalBorderSurface* grease_clad_LBS[64];
// G4LogicalBorderSurface* clad_core_LBS[64];
// G4LogicalSkinSurface* greaseLSS[4];
// G4LogicalSkinSurface* pmmaLSS[4];
// G4LogicalSkinSurface* polyLSS[4];
// G4LogicalSkinSurface* alum1LSS;
// G4LogicalSkinSurface* alum2LSS;
// G4LogicalSkinSurface* quartzLSS;
G4OpticalSurface* m_photonDetSurface;
bool m_simCherenkov;
G4SubtractionSolid* m_tile;
G4VSolid* m_tile_no_fiber_hole;
G4VSolid* m_fiber_subtract;
G4LogicalVolume* m_tileLogical;
G4VPhysicalVolume* m_tilePhysical[4][4];
G4VSolid* m_foilV[4];
G4LogicalVolume* m_foilVLogical[4];
G4VPhysicalVolume* m_foilVPhysical[4][4];
G4VSolid* m_foilVfront;
G4LogicalVolume* m_foilVfrontLogical;
G4VPhysicalVolume* m_foilVfrontPhysical[4][4];
G4VSolid* m_foilH;
G4LogicalVolume* m_foilHLogical;
G4VPhysicalVolume* m_foilHPhysical[4][4];
G4SubtractionSolid* m_foilHtop_hole;
G4VSolid* m_foilHtop;
G4LogicalVolume* m_foilHtopLogical;
G4VPhysicalVolume* m_foilHtopPhysical[4];
G4VSolid* m_AlcaseV;
G4LogicalVolume* m_AlcaseVLogical;
G4VPhysicalVolume* m_AlcaseVPhysical[5];
G4VSolid* m_Alcase;
G4LogicalVolume* m_AlcaseLogical;
G4VPhysicalVolume* m_AlcasePhysical[2];
G4VSolid* m_fiber[4];
G4LogicalVolume* m_fiberLogical[64];//4
G4VPhysicalVolume* m_fiberPhysical[64];
G4VSolid* m_fiberclad[4];
G4LogicalVolume* m_fibercladLogical[4];
G4VPhysicalVolume* m_fibercladPhysical[64];
G4VSolid* m_fibergrease[4];
G4LogicalVolume* m_fibergreaseLogical[4];
G4VPhysicalVolume* m_fibergreasePhysical[64];
G4VSolid* m_air_detect;
G4LogicalVolume* m_air_detect_Logical[64];
G4VPhysicalVolume* m_air_detectPhysical[64];
//test_setup
G4VSolid* m_test_tile;
G4LogicalVolume* m_test_tileLogical;
G4VPhysicalVolume* m_test_tilePhysical;
G4VSolid* m_test_alum;
G4LogicalVolume* m_test_alumLogical;
G4VPhysicalVolume* m_test_alumPhysical;
G4VSolid* m_test_wls;
G4LogicalVolume* m_test_wlsLogical;
G4VPhysicalVolume* m_test_wlsPhysical;
G4VSolid* m_test_PD;
G4LogicalVolume* m_test_PDLogical;
G4VPhysicalVolume* m_test_PDPhysical;
G4VSolid* m_test_clad;
G4LogicalVolume* m_test_cladLogical;
G4VPhysicalVolume* m_test_cladPhysical[2];
G4VSolid* m_test_grease;
G4LogicalVolume* m_test_greaseLogical;
G4VPhysicalVolume* m_test_greasePhysical[2];
G4VSolid* m_test_block;
G4LogicalVolume* m_test_blockLogical;
G4VPhysicalVolume* m_test_blockPhysical[2];
//pan flute rpd
G4VSolid* m_PFrpd[4];
G4LogicalVolume* m_PFrpdLogical[512];
G4VPhysicalVolume* m_PFrpdPhysical[512];
G4int m_PFrpd_cnt;
G4VSolid* m_PFrpd_foil[4];
G4LogicalVolume* m_PFrpd_foilLogical[512];
G4VPhysicalVolume* m_PFrpd_foilPhysical[512];
G4VSolid* m_PFdetec;
G4LogicalVolume* m_PFdetecLogical[512];
G4VPhysicalVolume* m_PFdetecPhysical[512];
//rpd booleans
bool rpd_comp[8];
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
MonteCarlo/src/FiberSD.cc
View file @ 224cb773
......@@ -143,6 +143,8 @@ if (config->GetValue("OPTICAL_ON",false) == 1){
newHit->setNCherenkovs (capturedPhotons);
// if (capturedPhotons != 0 ) std::cout << " capturedPhotons " << capturedPhotons << std::endl;
if( particle->GetPDGEncoding() == 0 ) fiberCollection->insert (newHit ); // only want to record photons
return true;
......
MonteCarlo/src/ModTypeRPD.cc
View file @ 224cb773
......@@ -140,6 +140,27 @@ void ModTypeRPD::DefineMaterials()
void ModTypeRPD::ConstructDetector()
{
bool pan_flute = true;
bool test_tile = false;
if(pan_flute || test_tile){
for(int i =0;i<8;i++){
rpd_comp[i]=false;
}
}
else for(int i =0;i<8;i++){
rpd_comp[i]=true;
}
// RPD components
bool tile_flag = rpd_comp[0]; // quartz tiles
bool fiber_flag = rpd_comp[1]; // readout fibers
bool cladding_flag = rpd_comp[2]; // fiber cladding
bool grease_flag = rpd_comp[3]; // optical grease
bool foil_flag = rpd_comp[4]; // aluminum foil
bool fr_bck_flag = rpd_comp[5]; // front and back aluminum casing
bool vertical_flag = rpd_comp[6]; // veritcal aluminum dividers
bool air_detec_flag = rpd_comp[7]; // photo detectors
bool make_Mike_happy = false;
// Get Config
......
MonteCarlo/src/ReactorTestFiber.cc 0 → 100644
View file @ 224cb773
//
// ********************************************************************
// * 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. *
// ********************************************************************
//
// Sheng Yang
// For an explanation of the hierarchy scheme see: https://twiki.cern.ch/twiki/bin/view/Atlas/ZdcSimulation#Geometry_Implementation_Develope
#include "ReactorTestFiber.hh"
#include "QuartzSD.hh"
#include "RpdSD.hh"
#include "FiberSD.hh"
#include "CherenkovSD.hh"
#include "SharedData.hh"
#include "G4GeometryManager.hh"
#include "G4SolidStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpticalSurface.hh"
#include "G4MaterialTable.hh"
#include "G4UserLimits.hh"
#include "G4RunManager.hh"
#include "G4SDManager.hh"
#include "G4NistManager.hh"
#include "G4CSGSolid.hh"
#include "G4Box.hh"
#include "G4Para.hh"
#include "G4Tubs.hh"
#include "G4Cons.hh"
#include "G4Orb.hh"
#include "G4Sphere.hh"
#include "G4Trd.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include <iostream>
#include <stdio.h>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ReactorTestFiber::ReactorTestFiber(const int cn,const G4ThreeVector& pos,
G4LogicalVolume* mother, SharedData* sd)
: m_modNum( cn ), m_pos( pos ), m_logicMother( mother ),
m_sd( sd ),
m_matQuartz(0),
m_PFrpd_cnt(0)
{
materials = Materials::getInstance();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ReactorTestFiber::ReactorTestFiber()
: m_modNum( 0 ), m_pos(G4ThreeVector()), m_logicMother(NULL),
m_sd(NULL),
m_matQuartz(0),
m_PFrpd_cnt(0)
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ReactorTestFiber::~ReactorTestFiber()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ReactorTestFiber::Construct(){
DefineMaterials();
ConstructDetector();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ReactorTestFiber::DefineMaterials()
{
materials->UseOpticalMaterials(true);
materials->DefineOpticalProperties();
//----------------------------------------------
// Define Materials
//----------------------------------------------
//Quartz
m_matQuartz = materials->pQuartz;
// m_quartzOS = materials->quartzOS;
// m_opSurface = materials->opSurface;
//Aluminum
m_Al = materials->Al;
//Air
m_Air = materials->Air;
//Polyethylene/Clad
m_Poly = materials->Polyethylene;
//Optical Grease
m_Grease = materials->Grease;
//Wavelength Shifter/Fiber
m_PMMA = materials->PMMA;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ReactorTestFiber::ConstructDetector()
{
// Get Config
TEnv* config = m_sd->GetConfig();
//retrieve RPD parameters
float tileX = config->GetValue("tileXsize",20);
float tileY = config->GetValue("tileYsize",20);
float fiber_diam = config->GetValue("mod6CoreDiameter",1);
float min_wall_thickness = config->GetValue("min_wall_thickness",1.0);
bool checkOverlaps = config->GetValue("checkOverlaps",false);
char name[256];
float RPD_centerX = m_pos.getX();
float RPD_centerY = m_pos.getY();
float RPD_centerZ = m_pos.getZ();
int n_rows = 1;
int n_columns = 1;
int n_cycles_per_tile = (tileX/fiber_diam)/n_rows; //Divide to round down to a whole number
int n_fibers_per_tile = n_cycles_per_tile*2*n_rows; //The pattern will be repeated twice per tile
//If you asked for a
if(.707*min_wall_thickness < fiber_diam){
min_wall_thickness = 1.414*fiber_diam;
std::cout << "Minimum wall thickness entered is too small!!!" << std::endl;
std::cout << "Using " << min_wall_thickness << "mm instead" << std::endl;
}
// If the remaining space in x leaves less than 1mm for each wall, remove a cycle and calculate
// the new wall thickness
calculate_wall_thickness:
float wall_thickness = ((tileX - n_fibers_per_tile*fiber_diam)/n_fibers_per_tile)*mm;
if( wall_thickness < min_wall_thickness*mm ){
--n_cycles_per_tile;
n_fibers_per_tile = n_cycles_per_tile*n_rows;
goto calculate_wall_thickness;
}
// Distance on center of fibers
float pitch = wall_thickness + fiber_diam;
// Distance in X and Z from one fiber to another diagonally from it
float offset = pitch/2;
// Distance between fiber and aluminum sheath
float gap = 0.01*mm;
// Will be filled with fiber height for each row
float fiber_height;
// Distance from the top of the RPD area of interest to the readout
float distance_to_readout = 5.0*mm;
// Thickness of the readout puck
float readout_thickness = 0.5*mm;
// Positions of the current fiber for pattern1 and pattern2
float posx1, posz1, posy, posx2, posz2;
// Count the number of fibers as we go along
m_PFrpd_cnt = 0;
//create some rotation matrices
G4RotationMatrix* stripRotation = new G4RotationMatrix();
stripRotation->rotateX(90.*deg);
G4RotationMatrix* nullRotation = new G4RotationMatrix();
// The SD caps that will go on the fibers. Each fiber will get a separate SD,
// so we wait to create the logical volume in the loop
m_PFdetec = new G4Tubs( "m_PFdetec",
0.0*mm,
(fiber_diam/2.0)*mm,
(readout_thickness/2.0)*mm ,
0.0*deg,
360.0*deg);
for(int row = 0; row < n_rows; row++){
fiber_height = tileY*(row+1) + distance_to_readout;
// The fiber solid. One length for each row
sprintf(name,"m_PFrpd_%d", row);
m_PFrpd[row] = new G4Tubs(name,
0.0*mm,
(fiber_diam/2.0)*mm,
fiber_height*mm/2.0 ,
0.0*deg,
360.0*deg);
//Cylindrical foil dividers
sprintf(name,"m_PFfoil_%d", row);
m_PFrpd_foil[row] = new G4Tubs(name,
(fiber_diam/2.0 + gap/2.0)*mm,
(fiber_diam/2.0 + gap/2.0 + wall_thickness/4)*mm,
fiber_height*mm/2.0 ,
0.0*deg,
360.0*deg);
sprintf(name,"m_PFfoil_log_%d",row);
m_PFrpd_foilLogical[row] = new G4LogicalVolume(m_PFrpd_foil[row]
,m_Al,
name);
m_PFrpd_foilLogical[row]->SetVisAttributes( G4Colour(1.0,0.0,0.0,0.3) );// G4Colour(1.0,0.0,0.0,0.3) or G4VisAttributes::Invisible
for(int col = 0; col < n_columns; col++){
//Now we're in the realm of working on a single tile
//we have to cycle through the two patterns until the tile
//is filled in X
for(int cycle = 0; cycle < n_cycles_per_tile; cycle++){
for(int fiber = 0; fiber < n_columns; fiber++){
// !!!!!!!Position calculations assume an even number of rows and columns!!!!!!! //////
//Start at RPD center + tile width* number of tiles + cycle number * cycle width + stack number in cycle * pitch
posx1 = RPD_centerX + tileX*((n_columns/2) - col ) - (cycle*n_columns*pitch) - fiber*pitch - pitch/4; //ARIC ADDED - pitch/4
posx2 = posx1 - offset;
// if(row==0){
// std::cout << "Rod# = " << m_PFrpd_cnt << ", posx1 = " << posx1 << std::endl;
// std::cout << "Rod# = " << m_PFrpd_cnt + 1 << ", posx2 = " << posx2 << std::endl;
// }
//Start at Z center - half the stack depth + moving front to back and looping back to front
posz1 = RPD_centerZ - pitch*n_columns/2 + pitch*((row + fiber )%4);
posz2 = RPD_centerZ - pitch*n_columns/2 + pitch*((row + fiber + 2)%4) + offset; //Pattern is offset by 2
//Start at RPDY center + distance to bottom of top tile + half the fiber height
posy = RPD_centerY + ((n_rows/2 - 1) - row)*tileY + fiber_height/2;
//----------------------- Place the rods -----------------------//
sprintf(name,"m_PFrpd_log_%d_%d_%d_%d_0",row,col,cycle,fiber);
m_PFrpdLogical[m_PFrpd_cnt] =
new G4LogicalVolume(m_PFrpd[row],
m_matQuartz,
name);
m_PFrpdLogical[m_PFrpd_cnt]->SetUserLimits(new G4UserLimits(DBL_MAX,DBL_MAX,10*ms));
sprintf(name,"m_PFrpd_pyhs_%d_%d_%d_%d_0",row,col,cycle,fiber);
m_PFrpdPhysical[m_PFrpd_cnt] =
new G4PVPlacement(stripRotation,
G4ThreeVector(posx1*mm, posy*mm, posz1*mm),
m_PFrpdLogical[m_PFrpd_cnt],
name,
m_logicMother,
false,
m_PFrpd_cnt,
checkOverlaps);
sprintf(name,"m_PFrpd_log_%d_%d_%d_%d_1",row,col,cycle,fiber);
m_PFrpdLogical[m_PFrpd_cnt+1] =
new G4LogicalVolume(m_PFrpd[row],
m_matQuartz,
name);
m_PFrpdLogical[m_PFrpd_cnt+1]->SetUserLimits(new G4UserLimits(DBL_MAX,DBL_MAX,10*ms));
sprintf(name,"m_PFrpd_pyhs_%d_%d_%d_%d_1",row,col,cycle,fiber);
m_PFrpdPhysical[m_PFrpd_cnt+1] =
new G4PVPlacement(stripRotation,
G4ThreeVector(posx2*mm, posy*mm, posz2*mm),
m_PFrpdLogical[m_PFrpd_cnt+1],
name,
m_logicMother,
false,
m_PFrpd_cnt+1,
checkOverlaps);
if(row==0){
m_PFrpdLogical[m_PFrpd_cnt]->SetVisAttributes( G4Colour::Cyan() );
m_PFrpdLogical[m_PFrpd_cnt+1]->SetVisAttributes( G4Colour::Cyan() );}
if(row==1){
m_PFrpdLogical[m_PFrpd_cnt]->SetVisAttributes( G4Colour::Red() );
m_PFrpdLogical[m_PFrpd_cnt+1]->SetVisAttributes( G4Colour::Red() );}
if(row==2){
m_PFrpdLogical[m_PFrpd_cnt]->SetVisAttributes( G4Colour::Green() );
m_PFrpdLogical[m_PFrpd_cnt+1]->SetVisAttributes( G4Colour::Green() );}
if(row==3){
m_PFrpdLogical[m_PFrpd_cnt]->SetVisAttributes( G4Colour::Magenta() );
m_PFrpdLogical[m_PFrpd_cnt+1]->SetVisAttributes( G4Colour::Magenta() );}
//----------------------- Place the walls/gap filler -----------------------//
sprintf(name,"m_PFfoil_pyhs_%d_%d_%d_%d_0",row,col,cycle,fiber);
m_PFrpd_foilPhysical[m_PFrpd_cnt] =
new G4PVPlacement(stripRotation,
G4ThreeVector(posx1*mm, posy*mm, posz1*mm),
m_PFrpd_foilLogical[row],
name,
m_logicMother,
false,
m_PFrpd_cnt,
checkOverlaps);
sprintf(name,"m_PFfoil_pyhs_%d_%d_%d_%d_1",row,col,cycle,fiber);
m_PFrpd_foilPhysical[m_PFrpd_cnt+1] =
new G4PVPlacement(stripRotation,
G4ThreeVector(posx2*mm, posy*mm, posz2*mm),
m_PFrpd_foilLogical[row],
name,
m_logicMother,
false,
m_PFrpd_cnt+1,
checkOverlaps);
//----------------------- Make and place the photodector volumes -----------------------//
sprintf(name,"m_PFdetec_log_%d_%d_%d_%d_0",row,col,cycle,fiber);
m_PFdetecLogical[m_PFrpd_cnt] =
new G4LogicalVolume(m_PFdetec,
m_matQuartz,
name);
m_PFdetecLogical[m_PFrpd_cnt]->SetVisAttributes(G4Colour(0.0,6.0,4.0,0.3));//G4Colour(1.0,0.0,0.0,0.3)
sprintf(name,"photo_det_phys_%d_%d_%d_%d_0",row,col,cycle,fiber);
m_PFdetecPhysical[m_PFrpd_cnt] =
new G4PVPlacement(nullRotation,
G4ThreeVector(0*mm, 0*mm, ((fiber_height/2-readout_thickness/2))*mm*mm),
m_PFdetecLogical[m_PFrpd_cnt],
name,
m_PFrpdLogical[m_PFrpd_cnt],
false,
m_PFrpd_cnt,
checkOverlaps);
sprintf(name,"m_PFdetec_log_%d_%d_%d_%d_1",row,col,cycle,fiber);
m_PFdetecLogical[m_PFrpd_cnt+1] =
new G4LogicalVolume(m_PFdetec,
m_matQuartz,