Commit e96d7708 authored by Warren's avatar Warren
Browse files

works again for particle source

parent 1ec0a9d2
Pipeline #105951 passed with stages
in 13 minutes and 34 seconds
......@@ -72,7 +72,7 @@ namespace TACTIC_NS{
//const double ResoEnergy = 1.0*MeV ;
// const double cathode_radius = 12.*mm;
const double drift_radius = 50.*mm;
//const double anode_radius = 51.*mm;
const double anode_radius = 51.*mm;
const double active_length = 251.9*mm;
const double window_pos = 104.*mm; //from centre of TACTIC from https://elog.triumf.ca/Tactic/Documentation/18
const double window_radius = 12.*mm; //guess
......@@ -127,7 +127,7 @@ void TACTIC::AddDetector(double R, double Theta, double Phi, string Shape){
G4LogicalVolume* TACTIC::BuildCylindricalDetector(){
if(!m_CylindricalDetector){
//G4Material* Cu = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
G4Material* Cu = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
G4Material* Mylar = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("G4_Galactic");
......@@ -168,6 +168,7 @@ G4LogicalVolume* TACTIC::BuildCylindricalDetector(){
G4Tubs* world_volume = new G4Tubs("world_volume",0,TACTIC_NS::drift_radius+1*mm,TACTIC_NS::active_length*0.5+1*mm,0,360*deg);
G4Tubs* window = new G4Tubs("window",0,TACTIC_NS::window_radius,TACTIC_NS::window_width*0.5,0,360*deg);
//G4Tubs* anode = new G4Tubs("anode",TACTIC_NS::drift_radius,TACTIC_NS::anode_radius,TACTIC_NS::active_length*0.5,0,360*deg);
G4Tubs* gas_cathode = new G4Tubs("gas_cathode",0,TACTIC_NS::window_radius,TACTIC_NS::window_pos-TACTIC_NS::window_width*0.5,0,360*deg); //window pos doesn't need halving
G4Tubs* gas_drift = new G4Tubs("gas_drift",TACTIC_NS::window_radius,TACTIC_NS::drift_radius,TACTIC_NS::active_length*0.5,0,360*deg);
G4UnionSolid* gas_volume = new G4UnionSolid("gas_volume",gas_cathode,gas_drift);
......@@ -175,16 +176,19 @@ G4LogicalVolume* TACTIC::BuildCylindricalDetector(){
//G4Tubs* gas_volume = new G4Tubs("gas_volume",0,TACTIC_NS::drift_radius,TACTIC_NS::active_length*0.5,0,360*deg);
m_CylindricalDetector = new G4LogicalVolume(world_volume, Vacuum, "m_CylindricalDetector_log",0,0,0);
//anode_log = new G4LogicalVolume(anode, Cu, "anode_log", 0,0,0);
gas_volume_log = new G4LogicalVolume(gas_volume, TACTIC_gas, "gas_volume_log",0,0,0);
window_log = new G4LogicalVolume(window, Mylar, "window_log",0,0,0);
new G4PVPlacement(0,G4ThreeVector(0,0,0),gas_volume_log,"gas_volume_phys",m_CylindricalDetector,false,0);
//new G4PVPlacement(0,G4ThreeVector(0,0,0),anode_log,"anode_phys",m_CylindricalDetector,false,0,true);
new G4PVPlacement(0,G4ThreeVector(0,0,TACTIC_NS::window_pos),window_log,"window_phys1",m_CylindricalDetector,false,0,true);
new G4PVPlacement(0,G4ThreeVector(0,0,-TACTIC_NS::window_pos),window_log,"window_phys2",m_CylindricalDetector,false,0,true);
gas_volume_log->SetVisAttributes(m_VisGas);
m_CylindricalDetector->SetVisAttributes(m_VisGas);
//anode_log->SetVisAttributes(m_VisGas);
G4UserLimits *gas_volume_step = new G4UserLimits();
G4double maxStep = 0.1*mm;
gas_volume_step->SetMaxAllowedStep(maxStep);
......@@ -305,11 +309,11 @@ void TACTIC::ReadSensitive(const G4Event* event ){
m_Event->Clear();
ofstream file;
/*
file.open("signal.dat", std::ios::app);
file << "Event" << endl;
file.close();
*/
file.open("out.dat",std::ios::app);
G4THitsMap<G4double*>* LightHitMap;
......@@ -319,11 +323,11 @@ void TACTIC::ReadSensitive(const G4Event* event ){
for (Light_itr = LightHitMap->GetMap()->begin(); Light_itr != LightHitMap->GetMap()->end(); Light_itr++) {
G4double* Info = *(Light_itr->second);
file << event->GetEventID() << "\t";
//file << event->GetEventID() << "\t";
for(int s = 0; s<11; s++) {
file << Info[s] << "\t";
}
file << Info[11] << endl;
file << Info[12] << endl;
}
LightHitMap->clear();
......@@ -335,11 +339,11 @@ void TACTIC::ReadSensitive(const G4Event* event ){
for (Heavy_itr = HeavyHitMap->GetMap()->begin(); Heavy_itr != HeavyHitMap->GetMap()->end(); Heavy_itr++) {
G4double* Info = *(Heavy_itr->second);
file << event->GetEventID() << "\t";
//file << event->GetEventID() << "\t";
for(int s = 0; s<11; s++) {
file << Info[s] << "\t";
}
file << Info[11] << endl;
file << Info[12] << endl;
}
HeavyHitMap->clear();
......@@ -351,11 +355,11 @@ void TACTIC::ReadSensitive(const G4Event* event ){
for (Beam_itr = BeamHitMap->GetMap()->begin(); Beam_itr != BeamHitMap->GetMap()->end(); Beam_itr++) {
G4double* Info = *(Beam_itr->second);
file << event->GetEventID() << "\t";
//file << event->GetEventID() << "\t";
for(int s = 0; s<11; s++) {
file << Info[s] << "\t";
}
file << Info[11] << endl;
file << Info[12] << endl;
}
BeamHitMap->clear();
......@@ -374,34 +378,58 @@ void TACTIC::InitializeScorers() {
return ;
// Otherwise the scorer is initialised
// NEED SEPERATE SCORERS FOR EACH PARTICLE! OTHERWISE SCORER JUST RETURNS OUTPUT FOR ONE PARTICLE WHEN THERE IS OVERLAP
G4VPrimitiveScorer* LightScorer = new TACTICScorer::Gas_Scorer("LightScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* HeavyScorer = new TACTICScorer::Gas_Scorer("HeavyScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* BeamScorer = new TACTICScorer::Gas_Scorer("BeamScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4SDParticleFilter* LightFilter = new G4SDParticleFilter("LightFilter");
G4SDParticleFilter* HeavyFilter = new G4SDParticleFilter("HeavyFilter"); //For studying alpha source data
G4SDParticleFilter* BeamFilter = new G4SDParticleFilter("BeamFilter");
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
NPL::Reaction m_Reaction;
m_Reaction.ReadConfigurationFile(input);
if(input.GetAllBlocksWithToken("TwoBodyReaction").size()>0) {
// NEED SEPERATE SCORERS FOR EACH PARTICLE! OTHERWISE SCORER JUST RETURNS OUTPUT FOR ONE PARTICLE WHEN THERE IS OVERLAP
G4VPrimitiveScorer* LightScorer = new TACTICScorer::Gas_Scorer("LightScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* HeavyScorer = new TACTICScorer::Gas_Scorer("HeavyScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* BeamScorer = new TACTICScorer::Gas_Scorer("BeamScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4SDParticleFilter* LightFilter = new G4SDParticleFilter("LightFilter");
G4SDParticleFilter* HeavyFilter = new G4SDParticleFilter("HeavyFilter");
G4SDParticleFilter* BeamFilter = new G4SDParticleFilter("BeamFilter");
/*
cout << m_Reaction.GetParticle1()->GetName() << "\t" << m_Reaction.GetParticle1()->GetA() << "\t" << m_Reaction.GetParticle1()->GetZ() << endl;
cout << m_Reaction.GetParticle2()->GetName() << "\t" << m_Reaction.GetParticle2()->GetA() << "\t" << m_Reaction.GetParticle2()->GetZ() << endl;
cout << m_Reaction.GetParticle3()->GetName() << "\t" << m_Reaction.GetParticle3()->GetA() << "\t" << m_Reaction.GetParticle3()->GetZ() << endl;
cout << m_Reaction.GetParticle4()->GetName() << "\t" << m_Reaction.GetParticle4()->GetA() << "\t" << m_Reaction.GetParticle4()->GetZ() << endl;
*/
LightFilter->addIon(m_Reaction.GetParticle3()->GetZ(),m_Reaction.GetParticle3()->GetA());
HeavyFilter->addIon(m_Reaction.GetParticle4()->GetZ(),m_Reaction.GetParticle4()->GetA());
BeamFilter->addIon(m_Reaction.GetParticle1()->GetZ(),m_Reaction.GetParticle1()->GetA());
LightScorer->SetFilter(LightFilter);
HeavyScorer->SetFilter(HeavyFilter);
BeamScorer->SetFilter(BeamFilter);
LightFilter->addIon(m_Reaction.GetParticle3()->GetZ(),m_Reaction.GetParticle3()->GetA());
HeavyFilter->addIon(m_Reaction.GetParticle4()->GetZ(),m_Reaction.GetParticle4()->GetA());
BeamFilter->addIon(m_Reaction.GetParticle1()->GetZ(),m_Reaction.GetParticle1()->GetA());
LightScorer->SetFilter(LightFilter);
HeavyScorer->SetFilter(HeavyFilter);
BeamScorer->SetFilter(BeamFilter);
m_Scorer->RegisterPrimitive(LightScorer);
m_Scorer->RegisterPrimitive(HeavyScorer);
m_Scorer->RegisterPrimitive(BeamScorer);
m_Scorer->RegisterPrimitive(LightScorer);
m_Scorer->RegisterPrimitive(HeavyScorer);
m_Scorer->RegisterPrimitive(BeamScorer);
}
if(input.GetAllBlocksWithToken("Isotropic").size()>0) { //For alpha, proton or neutron source (obviously neutron wont do anything so this is spare).
G4VPrimitiveScorer* LightScorer = new TACTICScorer::Gas_Scorer("LightScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* HeavyScorer = new TACTICScorer::Gas_Scorer("HeavyScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4VPrimitiveScorer* BeamScorer = new TACTICScorer::Gas_Scorer("BeamScorer",1,TACTIC_NS::active_length,(int)TACTIC_NS::NumberOfStrips);
G4SDParticleFilter* LightFilter = new G4SDParticleFilter("LightFilter","alpha");
G4SDParticleFilter* HeavyFilter = new G4SDParticleFilter("HeavyFilter","proton");
G4SDParticleFilter* BeamFilter = new G4SDParticleFilter("BeamFilter","neutron");
LightScorer->SetFilter(LightFilter);
HeavyScorer->SetFilter(HeavyFilter);
BeamScorer->SetFilter(BeamFilter);
m_Scorer->RegisterPrimitive(LightScorer);
m_Scorer->RegisterPrimitive(HeavyScorer);
m_Scorer->RegisterPrimitive(BeamScorer);
}
G4SDManager::GetSDMpointer()->AddNewDetector(m_Scorer);
}
......
......@@ -5,9 +5,9 @@ TACTIC
% GasMaterial_1 = P10_gas
GasMaterial_1 = He_gas
GasMaterial_2 = G4_CARBON_DIOXIDE
GasFraction_1 = 100
GasFraction_2 = 0
GasFraction_1 = 90
GasFraction_2 = 10
Temperature = 293.15 kelvin
Pressure = 0.5 bar
Pressure = 0.2 bar
% Pressure = 1.0 bar
Active = gas
\ No newline at end of file
......@@ -2,6 +2,7 @@
#include "TACTICScorer.hh"
#include "G4UnitsTable.hh"
#include "G4RunManager.hh"
#ifdef USE_Garfield
#include "GARFDRIFT.h"
......@@ -24,31 +25,33 @@ Gas_Scorer::~Gas_Scorer(){}
G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
G4double* Infos = new G4double[12];
G4double* Infos = new G4double[13];
m_Position = aStep->GetPreStepPoint()->GetPosition();
Infos[0] = aStep->GetTrack()->GetTrackID();
Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
Infos[1] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();;
Infos[1] = aStep->GetTrack()->GetTrackID();
Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();;
Infos[2] = aStep->GetPreStepPoint()->GetGlobalTime();
Infos[3] = aStep->GetPreStepPoint()->GetKineticEnergy();
Infos[4] = aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit();
Infos[3] = aStep->GetPreStepPoint()->GetGlobalTime();
Infos[4] = aStep->GetPreStepPoint()->GetKineticEnergy();
Infos[5] = aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit();
m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength ) + 1; //Pad number
Infos[5] = m_SegmentNumber;
Infos[6] = m_Position.z();
Infos[7] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R
Infos[8] = aStep->GetTrack()->GetVertexPosition()[2];
Infos[9] = aStep->GetTrack()->GetVertexKineticEnergy();
Infos[6] = m_SegmentNumber;
Infos[7] = m_Position.z();
Infos[8] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R
Infos[9] = aStep->GetTrack()->GetVertexPosition()[2];
Infos[10] = aStep->GetTrack()->GetVertexKineticEnergy();
G4ThreeVector p_vec = aStep->GetTrack()->GetVertexMomentumDirection();
Infos[10] = acos(p_vec[2]/pow(pow(p_vec[0],2)+pow(p_vec[1],2)+pow(p_vec[2],2),0.5))/deg; //angle relative to z axis (theta);
Infos[11] = aStep->GetTrack()->GetTrackLength();
Infos[11] = acos(p_vec[2]/pow(pow(p_vec[0],2)+pow(p_vec[1],2)+pow(p_vec[2],2),0.5))/deg; //angle relative to z axis (theta);
Infos[12] = aStep->GetTrack()->GetTrackLength();
m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6;
if(isnan(Infos[9])) {
if(isnan(Infos[10])) {
aStep->GetTrack()->SetTrackStatus(fStopAndKill);
return 0;
}
......@@ -56,7 +59,7 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
#ifdef USE_Garfield
G4ThreeVector delta_Position = aStep->GetDeltaPosition();
GARFDRIFT(Infos[4]/eV, Infos[2], m_Position/cm, delta_Position/cm, Infos[7]/cm, Infos[5], Infos[1], m_ScorerLength/cm, m_SegmentLength/cm);
GARFDRIFT(Infos[5]/eV, Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0], Infos[2], Infos[11]);
#endif
map<G4int, G4double**>::iterator it;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment