git clone https://github.com/saraheno/AlanHorstGEANTTile.git
cd AlanHostGEANTTile
source g4env.csh
mkdir build
cd build
cmake -DWITH_GEANT4_UIVIS=ON -DGeant4_DIR=$G4LIB ../UMDSRDGEStudy
make
cp -r ../UMDSRDGEStudy-utils/* .
./LYSim at the prompt. do the commands in photontest.mac
look at Analysis.root
code originally clonedfrom Alan Horst
We study the light output and collection efficiency of scintillators being considered for an update at the hadronic calorimeter of the CMS detector. Experimental data is useful for determining properties of different tiles, however, simulations using Geant4 could provide another more efficient way of testing different geometries and tiles. Here, I test a commercial scintillator by Eljen Technology the EJ-200 tile. It has a polyvinyltoluene (PVT) base, an emission peak of 425nm, a σ shaped groove with a wavelength shifting fiber inside connected to a Hamamatsu R6091 PMT and is wrapped in Tyvek.
-Fiber Geometry Update Our tile simulation was in need of a modified geometry for the EJ-200 material. This fiber was unlike the fiber used originally in the simulation. The main section I altered was the output corner, where the fiber leaves the tile, and where the opposite end of the fiber ends inside the tile. Originally, the output of the fiber was connected to the end of the fiber in a T like style on the output corner. This is not true of the EJ-200 tile, in which the end and output are not connected at this corner, nor is the end of the fiber straight. The EJ-200 fiber looks closer to a σ, where the output is not connected to the end of the fiber.
-Unirradiated & Irradiated Absorption Lengths The experimental researchers provided absorption lengths from tests with irradiated and unirradiated tiles. Using these experimental results, I was able to test uniformity of photon absorption in the tile based on position. This involved creating a function to record the photon source position to gather position based data on the amount of photons being detected by the unirradiated and irradiated tiles. Although partially used to compare the unirradiated tiles to the irradiated tiles, this process was also used to test the accuracy of the simulation itself.
-Uniformity Testing To test uniformity, I needed a method of creating a changing particle source to test specific points on the tile. Using the position commands and a basic python loop, I could run multiple jobs at different points. I also added code to the Analysis and PrimaryGeneratorAction source files and headers to record the source position for later use in data visualization and analysis.
-Light Yield Comparison using Cosmic Ray Muons To test the energy deposit of muons in the EJ-200 tile, I created a planar source of muons shot into the tile, and created histograms that displayed the deposited energy and the initial total energy of each muon from the source. I used two different source energy distributions: one from a sea level muon experimental measurement, and one using the Geant4 CDG (cosmic diffuse gamma ray) distribution.
The primary goal of this project is to modify/extend and verify out existing GEANT4 setup, including updating the tile geometry of the scintillator tile to the commercial EJ-200 geometry, by altering the tile geometry and the fiber groove. This change was made so I could more realistically describe real measurements, such as the cosmic ray, alpha source, and test beam measurements, which could then be used to explain behaviors in the real tiles. After validating these results, others can further develop this simulation based on my work to morph/extrapolate into more geometries for other tiles. By comparing/contrasting different types/geometries of tiles, one could easily find tiles best suited for use in the real detector based on radiation damage, longevity, and accuracy. These tiles are crucial to the measurements made at CMS, and better tiles will mean more accurate measurement with less noise. Using GEANT4 will be a straight forward method of finding tiles that best suit the detectors.
Much of this project also involved python, for submitting jobs and automatically replacing strings (coordinates, random seeds, root names, etc.) in macros used for giving Geant4 commands for each job. I also used python and root extensively for analysis and graphing of results.
Contributions to the main source code include: - A new fiber core geometry, with curved readout corner, and ends that do not connect at the readout corner. - A function for outputting the particle source position at the beginning of the run (useful for uniformity testing). - A function for outputting initial particle energy. (Useful for checking emission spectra) - New histograms (Source initial energy, Step energy)
- Each section is dictated by an uppercase title surrounded by ~~~~~ as seen below for the
analysis.cc
- **This is a quick summary of what was added**
- What follows is the section(s) of code including what was added. It does not include the
entire version of said .cc or header file.
ANALYSIS.CC
Added position variable G4ThreeVector (pos) to record particle source position for each run.
void Analysis::EndOfRun(const G4Run*) { outputfile.open(fOutputFileName.c_str(), ofstream::out | ofstream::app); G4ThreeVector pos = generatorAction->GetSourcePosition(); G4double detEff = (PhotonCount > 0 ? (G4double)HitCount/(G4double)PhotonCount: 0.0); G4cout << "Efficiency in this run is " << detEff << G4endl; if (outputfile.is_open()) { //outputfile << "#Mu_tile [cm^-1]\tMu_fiber [cm^-1]\tEfficiency" << G4endl; outputfile << inducedMuTile << "\t" << inducedMuFiber << "\t" << detEff << "\t" << pos.x() << "\t" << pos.y() << G4endl; } else { G4cout << "Output file not open" << G4endl; //G4cout << "#Mu_tile [cm^-1]\tMu_fiber [cm^-1]\tEfficiency" << G4endl; G4cout << inducedMuTile << "\t" << inducedMuFiber << "\t" << detEff << "\t" << pos.x() << "\t" << pos.y() << G4endl; } outputfile.close();
// Save histograms
G4AnalysisManager* man = G4AnalysisManager::Instance();
man->Write();
man->CloseFile();
}
ANALYSIS.HH
Added pointers to LYSimPrimaryGeneratorAction.cc
class Analysis { public:
};
//Set Pointer to Generator Action
void SetGeneratorAction(LYSimPrimaryGeneratorAction* genaction){
generatorAction = genaction;
};
private: //Pointer to GeneratorAction class for access to source properties LYSimPrimaryGeneratorAction* generatorAction;
LYSIMRUNACTION.CC
**Added Total/Step Muon energy histograms
void LYSimRunAction::BeginOfRunAction(const G4Run* aRun) { Analysis::GetInstance()->PrepareNewRun(aRun); G4AnalysisManager* man = G4AnalysisManager::Instance(); outFileName = Analysis::GetInstance()->GetROOTFileName(); G4cout << "Output filename: " << outFileName << G4endl; man->OpenFile(outFileName.c_str()); man->SetFirstHistoId(1);
// Create histogram(s) (avoid non-integer bins)
if (pDetectorConstruction->GetDetectorType()==1) {
man->CreateH1("h1","Optical photons energy [eV]", //histoID,histo name
100,0.,10.); //bins' number, xmin, xmax
man->CreateH1("h2","Number of detected photons per event",
100,0.,100.); //bins' number, xmin, xmax
man->CreateH1("h3","Total optical photons energy deposited per event [eV]",
500,0.,500.); //bins' number, xmin, xmax
man->CreateH1("h4","Muon energy deposited per step [keV]",
50,0.,5.); //bins' number, xmin, xmax
man->CreateH1("h5","Total Muon total energy[MeV]",
100,0.,100.); //bins' number, xmin, xmax
} else {
man->CreateH1("h1","Optical photons energy [eV]", //histoID,histo name
100,0.,5.); //bins' number, xmin, xmax
man->CreateH1("h2","Number of detected photons per event",
250,0.,250.); //bins' number, xmin, xmax
man->CreateH1("h3","Total optical photons energy deposited per event [eV]",
100,0.,5.); //bins' number, xmin, xmax
man->CreateH1("h4","Muon energy deposited per step [keV]",
50,0.,5.); //bins' number, xmin, xmax
man->CreateH1("h5","Total Muon total energy[MeV]",
100,0.,100.); //bins' number, xmin, xmax
}
}
LYSIMEVENTACTION.CC
Added initial kenetic energy function, filled h5 histogram
void LYSimEventAction::BeginOfEventAction(const G4Event* anEvent ) { G4AnalysisManager* man = G4AnalysisManager::Instance(); G4PrimaryParticle* primary = anEvent->GetPrimaryVertex(0)->GetPrimary(0); //G4cout << G4endl // << ">>> Event " << anEvent->GetEventID() << " >>> Simulation truth : " // << primary->GetG4code()->GetParticleName() // << " " << primary->GetMomentum() << G4endl; G4PrimaryVertex* primaryVertex = anEvent->GetPrimaryVertex(); G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary(); G4double ke = primaryParticle->GetKineticEnergy(); man->FillH1(5,ke); //total energy of source if ( anEvent->GetEventID() % 100 == 0 ) { G4cout<<"Starting Event: "<GetEventID()<<G4endl; } Analysis::GetInstance()->PrepareNewEvent(anEvent); //Retrieve the ID for the hit collection //if ( hitsCollID == -1 ) // { // G4SDManager * SDman = G4SDManager::GetSDMpointer(); // hitsCollID = SDman->GetCollectionID(hitsCollName); // } }
LYSIMSCINTILLATION.CC
Added h4 histogram fill with step energy function
LYSimScintillation::LYSimScintillation(const G4String &processName, G4ProcessType type) : G4Scintillation(processName, type) { G4AnalysisManager::Instance(); }
G4VParticleChange* LYSimScintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) { G4VParticleChange result = G4Scintillation::PostStepDoIt(aTrack, aStep); G4double depenergy = aStep.GetTotalEnergyDeposit(); G4ParticleDefinition particleType = aTrack.GetDefinition(); G4String particleName = particleType->GetParticleName(); G4AnalysisManager* man = G4AnalysisManager::Instance(); if (depenergy > 0.0) { G4ThreeVector pos = aStep.GetPreStepPoint()->GetPosition(); G4cout << "[LYSim] Energy Deposited keV: " << depenergy << G4endl; man->FillH1(4,depenergy/keV); //energy deposited per step
// if (outFile.is_open()) { // outFile << "# scintillating: " << depenergy/keV << " keV of " // << std::setprecision(4) // << std::setw(10) << depenergy/keV << " keV of " // << std::setw(9) << aStep.GetPreStepPoint()->GetKineticEnergy()/keV << " keV [kine] " // << std::setw(10) << aStep.GetPreStepPoint()->GetTotalEnergy()/keV << " keV [total] " // << "deposited at (" // << std::fixed // << std::setw(7) << pos.x()/mm << " mm," // << std::setw(7) << pos.y()/mm << " mm," // << std::setw(7) << pos.z()/mm << " mm) " // << "by parent ID " << std::setw(5) << aTrack.GetTrackID() << " : " // << std::setw(5) << particleName << " " // << "producing " << result->GetNumberOfSecondaries() << " optical photons" // << std::endl // << std::resetiosflags(std::ios::fixed); // } else { G4cout << "scintillating: " << std::setprecision(4) << std::setw(10) << depenergy/keV << " keV of " << std::setw(9) << aStep.GetPreStepPoint()->GetKineticEnergy()/keV << " keV [kine] " << std::setw(10) << aStep.GetPreStepPoint()->GetTotalEnergy()/keV << " keV [total] " << "deposited at (" << std::fixed << std::setw(7) << pos.x()/mm << " mm," << std::setw(7) << pos.y()/mm << " mm," << std::setw(7) << pos.z()/mm << " mm) " << "by parent ID " << std::setw(5) << aTrack.GetTrackID() << " : " << std::setw(5) << particleName << " " << "producing " << result->GetNumberOfSecondaries() << " optical photons" << G4endl << std::resetiosflags(std::ios::fixed); //} } return result; }
LYSIMPRIMARYGENERATORACTION.CC
Added particleSource function to provide pos G4ThreeVector for Analysis use.
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
LYSimPrimaryGeneratorAction::LYSimPrimaryGeneratorAction(LYSimDetectorConstruction* det) : PhotonEnergy(2.95eV), GammaEnergy(660keV), BetaEnergy(511*keV) { fDetector = det; particleSource = new G4GeneralParticleSource();
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = G4OpticalPhoton::OpticalPhotonDefinition();
particleSource->SetParticleDefinition(particle);
particleSource->SetParticleTime(0.0*ns);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
LYSimPrimaryGeneratorAction::~LYSimPrimaryGeneratorAction() { delete particleSource; }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void LYSimPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { if (particleSource->GetParticleDefinition()->GetParticleName() == "opticalphoton") { SetOptPhotonPolar(); } G4ThreeVector pos = particleSource->GetParticlePosition();
particleSource->GeneratePrimaryVertex(anEvent);
//Analysis
//Analysis::GetInstance()->AddPhotonCount(1);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
const G4ThreeVector LYSimPrimaryGeneratorAction::GetSourcePosition() { G4ThreeVector pos = particleSource->GetParticlePosition(); return pos; }
//....oooOO-1OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
LYSIMPRIMARYGENERATORACTION.HH
** Added functions/constants for use in Analysis.cc for particle source position**
class LYSimPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction { public: LYSimPrimaryGeneratorAction(LYSimDetectorConstruction*); ~LYSimPrimaryGeneratorAction();
public: void GeneratePrimaries(G4Event*);
void SetOptPhotonPolar();
void SetOptPhotonPolar(G4double);
const G4ThreeVector GetSourcePosition();
private:
G4GeneralParticleSource* particleSource;
LYSimDetectorConstruction* fDetector;
G4double GammaEnergy;
G4double BetaEnergy;
G4double PhotonEnergy;
};
LYSIMDETECTORCONSTRUCTION.CC
Added new geometry based on previously defined lengths for simplicity
//Curved Fiber -- Alan Horst
G4VSolid* LYSimDetectorConstruction::ConstructFiberSolid(const G4String& name, G4double radiusI, G4double radiusO, G4double bendRadius, G4double distance, G4double angle1, G4double angle2, G4int readoutCorner) { G4double l1 = (bendRadius + distance); G4double l2 = (bendRadius + distance);
/*
(6) ____________________________ (7) / /| (2)/_2(3)/ | | | | | | | | | | | | | 3 | 1 | | | | | | | | | | |(4) | |(5) (0)|/0(1)|/
Basic corner labels for tile shown in (). Corners (0,1,2,3) are consistent for centers, curved fiber components.
Straight component labels are labeled along the sides of the tile Ex: _0_
Previously Defined Lengths:
- radiusI => Inner torus radius
- radiusO => Outer torus radius
- bendRadius => Torus swept radius
- distance => Distance from tile edge to fiber center
Readout Corner can be at corner (0) or corner (1)
Currently, only square tile shape is supported, the other fiber geometry has trapezoidal options.
*/
G4ThreeVector center[4]; //Center of curvature for tori at each corner
center[0] = 0.5*(corners[0] + corners[4]) + G4ThreeVector((l1), bendRadius + distance, 0);
center[1] = 0.5*(corners[1] + corners[5]) + G4ThreeVector(-1*(l2), bendRadius + distance, 0);
center[2] = 0.5*(corners[2] + corners[6]) + G4ThreeVector((l1), -bendRadius - distance, 0);
center[3] = 0.5*(corners[3] + corners[7]) + G4ThreeVector(-1*(l2), -bendRadius - distance, 0);
G4ThreeVector endpointsA[4];
G4ThreeVector endpointsB[4];
G4ThreeVector centersStraight[4];
G4double lengthsStraight[4];
endpointsA[0] = center[0] + G4ThreeVector(0, -bendRadius, 0);
endpointsA[1] = center[1] + G4ThreeVector(0, -bendRadius, 0);
endpointsA[2] = center[2] + G4ThreeVector(0, +bendRadius, 0);
endpointsA[3] = center[3] + G4ThreeVector(0, +bendRadius, 0);
endpointsB[0] = center[0] + G4ThreeVector(-bendRadius, 0, 0);
endpointsB[1] = center[1] + G4ThreeVector(+bendRadius, 0, 0);
endpointsB[2] = center[2] + G4ThreeVector(-bendRadius, 0, 0);
endpointsB[3] = center[3] + G4ThreeVector(+bendRadius, 0, 0);
centersStraight[0] = (endpointsA[0] + endpointsA[1]) / 2;
centersStraight[1] = (endpointsB[1] + endpointsB[3]) / 2;
centersStraight[2] = (endpointsA[3] + endpointsA[2]) / 2;
centersStraight[3] = (endpointsB[2] + endpointsB[0]) / 2;
lengthsStraight[0] = (endpointsA[0] - endpointsA[1]).mag()-0.0*mm;
lengthsStraight[1] = (endpointsB[1] - endpointsB[3]).mag()-0.0*mm;
lengthsStraight[2] = (endpointsA[3] - endpointsA[2]).mag()-0.0*mm;
lengthsStraight[3] = (endpointsB[2] - endpointsB[0]).mag()-0.0*mm;
G4RotationMatrix* rotStraight0 = new G4RotationMatrix;
G4RotationMatrix* rotStraight1 = new G4RotationMatrix;
G4RotationMatrix* rotStraight2 = new G4RotationMatrix;
G4RotationMatrix* rotStraight3 = new G4RotationMatrix;
rotStraight0->rotateY(pi/2*rad); //z to x
rotStraight0->invert();
rotStraight1->rotateY(pi/2*rad); //z to x
rotStraight1->rotateZ(pi/2*rad);
rotStraight1->invert();
rotStraight2->rotateY(pi/2*rad); //z to x
rotStraight2->invert();
rotStraight3->rotateY(pi/2*rad); //z to x
rotStraight3->rotateZ(pi/2*rad);
rotStraight3->invert();
G4double lengthToEdge0, lengthToEdge1;
G4double lengthReadoutSection0, lengthReadoutSection1;
G4ThreeVector centerReadoutSection0, centerReadoutSection1;
lengthToEdge0 = (endpointsB[0].y() - corners[0].y());
lengthToEdge1 = (endpointsB[1].y() - corners[1].y());
lengthReadoutSection0 = 2*lengthToEdge0;
lengthReadoutSection1 = 2*lengthToEdge1;
readout0 = endpointsB[0] + lengthReadoutSection0 * G4ThreeVector(0, -1, 0);
readout1 = endpointsB[1] + lengthReadoutSection1 * G4ThreeVector(0, -1, 0);
centerReadoutSection0 = 0.5*(endpointsB[0] + readout0);
centerReadoutSection1 = 0.5*(endpointsB[1] + readout1);
G4VSolid* solidFiberCurved0 =
new G4Torus(name+"CurvedSection0",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
1.0*pi, //G4double pSPhi,
0.5*pi); //G4double pDPhi)
G4VSolid* solidFiberCurved1 =
new G4Torus(name+"CurvedSection1",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
1.5*pi + 0.0000, //G4double pSPhi,
0.5*pi); //G4double pDPhi)
G4VSolid* solidFiberCurved2 =
new G4Torus(name+"CurvedSection2",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
0.5*pi + 0.0000, //G4double pSPhi,
0.5*pi); //G4double pDPhi)
G4VSolid* solidFiberCurved3 =
new G4Torus(name+"CurvedSection3",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
0.0*pi, //G4double pSPhi,
0.5*pi); //G4double pDPhi)
G4VSolid* solidFiberStraight0 =
new G4Tubs(name+"StraightSection0",
radiusI,
radiusO,
0.5*lengthsStraight[0],
0,
2.*pi);
G4VSolid* solidFiberStraight1 =
new G4Tubs(name+"StraightSection1",
radiusI,
radiusO,
0.5*lengthsStraight[1],
0,
2.*pi);
G4VSolid* solidFiberStraight2 =
new G4Tubs(name+"StraightSection2",
radiusI,
radiusO,
0.5*lengthsStraight[2],
0,
2.*pi);
G4VSolid* solidFiberStraight3 =
new G4Tubs(name+"StraightSection3",
radiusI,
radiusO,
0.5*lengthsStraight[3],
0,
2.*pi);
G4VSolid* solidReadoutSection0 =
new G4Tubs(name+"ReadoutSection0",
radiusI,
radiusO,
0.5*lengthReadoutSection0,
0,
2.*pi);
G4VSolid* solidReadoutSection1 =
new G4Tubs(name+"ReadoutSection1",
radiusI,
radiusO,
0.5*lengthReadoutSection1,
0,
2.*pi);
G4Box* solidFiberBase =
new G4Box(name+"FiberBase",
0.5*world_sizeXY, 0.5*world_sizeXY, 0.5*world_sizeZ);
G4VSolid* solidPrev = solidFiberBase;
G4RotationMatrix* rotIdentity = new G4RotationMatrix;
if(readoutCorner != 0)
{
G4VSolid* solidFiberComponent0 =
new G4IntersectionSolid(name+"Component0",
solidPrev,
solidFiberCurved0,
rotIdentity,
center[0]);
solidPrev = solidFiberComponent0;
}
//Straight components
G4VSolid* solidFiberComponent4 =
new G4UnionSolid(name+"Component4",
solidPrev,
solidFiberStraight0,
rotStraight0,
centersStraight[0]);
solidPrev = solidFiberComponent4;
if(readoutCorner != 1)
{
G4VSolid* solidFiberComponent1 =
new G4UnionSolid(name+"Component1",
solidPrev,
solidFiberCurved1,
rotIdentity,
center[1]);
solidPrev = solidFiberComponent1;
}
G4VSolid* solidFiberComponent5 =
new G4UnionSolid(name+"Component5",
solidPrev,
solidFiberStraight1,
rotStraight1,
centersStraight[1]);
solidPrev = solidFiberComponent5;
G4VSolid* solidFiberComponent2 =
new G4UnionSolid(name+"Component2",
solidPrev,
solidFiberCurved2,
rotIdentity,
center[2]);
solidPrev = solidFiberComponent2;
G4VSolid* solidFiberComponent6 =
new G4UnionSolid(name+"Component6",
solidPrev,
solidFiberStraight2,
rotStraight2,
centersStraight[2]);
solidPrev = solidFiberComponent6;
G4VSolid* solidFiberComponent3 =
new G4UnionSolid(name+"Component3",
solidPrev,
solidFiberCurved3,
rotIdentity,
center[3]);
solidPrev = solidFiberComponent3;
G4VSolid* solidFiberComponent7 =
new G4UnionSolid(name+"Component7",
solidPrev,
solidFiberStraight3,
rotStraight3,
centersStraight[3]);
solidPrev = solidFiberComponent7;
if(readoutCorner == 0)
{
G4VSolid* solidFiberComponent8 =
new G4UnionSolid(name+"Component8",
solidPrev,
solidReadoutSection0,
rotStraight3,
centerReadoutSection0);
solidPrev = solidFiberComponent8;
//Curved Readout
G4VSolid* solidFiberCurved0R =
new G4Torus(name+"CurvedSection0R",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
1.15*pi, //G4double pSPhi,
0.35*pi); //G4double pDPhi)
G4VSolid* solidFiberComponent9 =
new G4UnionSolid(name+"Component9",
solidPrev,
solidFiberCurved0R,
rotIdentity,
center[0]);
solidPrev = solidFiberComponent9;
}
else if(readoutCorner == 1)
{
G4VSolid* solidFiberComponent8 =
new G4UnionSolid(name+"Component8",
solidPrev,
solidReadoutSection1,
rotStraight1,
centerReadoutSection1);
solidPrev = solidFiberComponent8;
//Curved Readout
G4VSolid* solidFiberCurved1R =
new G4Torus(name+"CurvedSection1R",
radiusI, //G4double pRmin,
radiusO, //G4double pRmax,
bendRadius, //G4double pRtor,
1.5*pi, //G4double pSPhi,
0.3*pi); //G4double pDPhi)
G4VSolid* solidFiberComponent9 =
new G4UnionSolid(name+"Component9",
solidPrev,
solidFiberCurved1R,
rotIdentity,
center[1]);
solidPrev = solidFiberComponent9;
}
return solidPrev;
}