Skip to content

GitLab

Commit 6b557cd8 authored by Warren's avatar Warren
Browse files

updates in preperation for running on university system

parent 31d48849
Pipeline #123599 passed with stages
in 3 minutes and 35 seconds
Show whitespace changes
Inline Side-by-side
Showing with 102 additions and 45 deletions
NPSimulation/Detectors/TACTIC/CMakeLists.txt
View file @ 6b557cd8
...@@ -12,7 +12,7 @@ endif() ...@@ -12,7 +12,7 @@ endif()
add_library(NPSTACTIC SHARED TACTIC.cc TACTICScorer.cc) add_library(NPSTACTIC SHARED TACTIC.cc TACTICScorer.cc)
if(WITH_GARFIELD) if(WITH_GARFIELD)
target_link_libraries(NPSTACTIC NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPTACTIC /Users/warren/Garfield/install/lib/libGarfield.0.3.0.dylib /Users/warren/Garfield/install/lib/libGarfield.dylib) target_link_libraries(NPSTACTIC NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPTACTIC /Users/warren/garfieldpp/install/lib/libGarfield.0.3.0.dylib /Users/warren/garfieldpp/install/lib/libGarfield.dylib)
else() else()
target_link_libraries(NPSTACTIC NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPTACTIC) target_link_libraries(NPSTACTIC NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPTACTIC)
endif() endif()
NPSimulation/Detectors/TACTIC/GARFDRIFT.h
View file @ 6b557cd8
...@@ -9,9 +9,14 @@ ...@@ -9,9 +9,14 @@
#include "Garfield/ViewMedium.hh" #include "Garfield/ViewMedium.hh"
double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector delta_pos, double R, int Pad_start, double ID, double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector delta_pos, double R, int Pad_start, double ID,
double ScoLength, double SegLength, double event) { double ScoLength, double SegLength, double event, double raw_energy_check) {
ofstream file;
/*
file.open("garf_E_check.dat", std::ios::app);
file << raw_energy_check << "\t" << R << endl;
file.close();
*/
const double rWire = 1.2; const double rWire = 1.2;
const double rTube = 5.; const double rTube = 5.;
const double lTube = 25.19; const double lTube = 25.19;
...@@ -24,17 +29,22 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect ...@@ -24,17 +29,22 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect
//const double vWire = -681.0312; //For V_total = 1400V //const double vWire = -681.0312; //For V_total = 1400V
//const double vWire = -730.; //const double vWire = -730.;
const double vWire = -1000.; const double vWire = -1000.;
//const double vWire = -584.7139; // ND run 4230 (2008)
const double vTube = 0.; const double vTube = 0.;
double x_start, y_start, z_start; double x_start, y_start, z_start;
ofstream file; //double production_bias = 0.01;
double production_bias = 0.001; double production_bias = 0.001;
double t_end, x_end, y_end, z_end; double t_end, x_end, y_end, z_end;
int status;
double Pad_end; double Pad_end;
//double w_value = 35.; //for HeCO2 --guess //double w_value = 35.; //for HeCO2 --guess
//double w_value = 26.31; //for argon //double w_value = 26.31; //for argon
double w_value = 41.1;// for He double w_value = 41.1;// for He
int electrons = (int)(energy/w_value*production_bias);//*production_bias; //Need to carry excess energy over later. energy = energy*production_bias;
//int electrons = (int)(energy/w_value*production_bias);
int electrons = (int)(energy/w_value);
double energy_excess = 0.; double energy_excess = 0.;
int sector;
vector<double> pad_vec, time_vec; vector<double> pad_vec, time_vec;
Garfield::MediumMagboltz* gas = new Garfield::MediumMagboltz(); Garfield::MediumMagboltz* gas = new Garfield::MediumMagboltz();
...@@ -45,17 +55,25 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect ...@@ -45,17 +55,25 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect
Garfield::Sensor* sensor = new Garfield::Sensor(); Garfield::Sensor* sensor = new Garfield::Sensor();
Garfield::AvalancheMC* drift = new Garfield::AvalancheMC(); Garfield::AvalancheMC* drift = new Garfield::AvalancheMC();
cout << "\n" << electrons << "\n" << endl; //cout << "\n" << electrons << "\n" << endl;
if(R>rWire && electrons>0) { energy_excess = (energy - electrons*w_value)/production_bias;
// if(R>rWire) {
energy_excess = (energy/w_value*production_bias - electrons)*w_value/production_bias; // energy_excess = (energy/w_value*production_bias - electrons)*w_value/production_bias;
// energy_excess = (energy/w_value - electrons)*w_value;
if(R>rWire && electrons>0) {
//gas->LoadGasFile("he_90_co2_10_200mbar.gas"); //gas->LoadGasFile("he_90_co2_10_200mbar.gas");
//gas->LoadGasFile("ar_90_ch4_10.gas"); //atomic fraction in this case (P10). //gas->LoadGasFile("ar_90_ch4_10.gas"); //atomic fraction in this case (P10).
//gas->LoadGasFile("he_90_co2_10_500mbar.gas"); //mass fraction in this case //gas->LoadGasFile("he_90_co2_10_500mbar.gas"); //mass fraction in this case
//gas->LoadGasFile("he_90_co2_10_1bar.gas"); //gas->LoadGasFile("he_90_co2_10_1bar.gas");
gas->LoadGasFile("he_90_co2_10_100mbar.gas"); gas->LoadGasFile("he_90_co2_10_100mbar.gas");
//gas->LoadGasFile("he_90_co2_10_600mbar.gas");
gas->Initialise(true); gas->Initialise(true);
...@@ -67,6 +85,8 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect ...@@ -67,6 +85,8 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect
sensor->AddComponent(cmp); sensor->AddComponent(cmp);
drift->SetSensor(sensor); drift->SetSensor(sensor);
//drift->DisableAttachment();
for(int e=0;e<electrons;e++) { for(int e=0;e<electrons;e++) {
double randomize = (double)std::rand() / (double)RAND_MAX; double randomize = (double)std::rand() / (double)RAND_MAX;
...@@ -74,23 +94,51 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect ...@@ -74,23 +94,51 @@ double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVect
x_start = start.x() + delta_pos.x()*randomize; x_start = start.x() + delta_pos.x()*randomize;
y_start = start.y() + delta_pos.y()*randomize; y_start = start.y() + delta_pos.y()*randomize;
z_start = start.z() + delta_pos.z()*randomize; z_start = start.z() + delta_pos.z()*randomize;
/*
file.open("garf_electrons_start_RZ.dat", std::ios::app);
file << pow(pow(x_start,2)+pow(y_start,2),0.5) << "\t" << z_start << endl;
file.close();
*/
drift->DriftElectron(x_start, y_start, z_start, t_start); drift->DriftElectron(x_start, y_start, z_start, t_start);
int np = drift->GetNumberOfDriftLinePoints(); int np = drift->GetNumberOfDriftLinePoints();
//drift->GetEndPoint(x_end, y_end, z_end, t_end, status);
drift->GetDriftLinePoint(np-1, x_end, y_end, z_end, t_end); //np-1 is last DriftLineEndPoint drift->GetDriftLinePoint(np-1, x_end, y_end, z_end, t_end); //np-1 is last DriftLineEndPoint
//int np = drift->GetNuberOfElectronEndPoints();
//drift->GetElectronEndPoint(np-1, x_end, y_end, z_end, t_end);
Pad_end = (int)((z_end + ScoLength / 2.) / SegLength ) + 1; //new Pad number Pad_end = (int)((z_end + ScoLength / 2.) / SegLength ) + 1; //new Pad number
if(pow(pow(x_end,2)+pow(y_end,2),0.5)>4.9) { //reached the anode if(x_end > 0 and y_end > 0) {
if(abs(x_end) > abs(y_end)) sector = 0;
if(abs(x_end) < abs(y_end)) sector = 1;
}
if(x_end < 0 and y_end > 0) {
if(abs(x_end) < abs(y_end)) sector = 2;
if(abs(x_end) > abs(y_end)) sector = 3;
}
if(x_end < 0 and y_end < 0) {
if(abs(x_end) > abs(y_end)) sector = 4;
if(abs(x_end) < abs(y_end)) sector = 5;
}
if(x_end > 0 and y_end < 0) {
if(abs(x_end) < abs(y_end)) sector = 6;
if(abs(x_end) > abs(y_end)) sector = 7;
}
//if(pow(pow(x_end,2)+pow(y_end,2),0.5)>4.9) { //reached the anode
file.open("signal.dat", std::ios::app); file.open("signal.dat", std::ios::app);
file << event << "\t" << ID << "\t" << Pad_start << "\t" << Pad_end << "\t" << t_end << endl; file << event << "\t" << ID << "\t" << sector << "\t" << Pad_end << "\t" << t_end << "\t" << pow(pow(x_start,2)+pow(y_start,2),0.5) << "\t" << pow(pow(x_end,2)+pow(y_end,2),0.5) << endl;
//file << event << "\t" << ID << "\t" << sector << "\t" << Pad_end << "\t" << t_end
file.close(); file.close();
} //}
} }
} }
//}
delete gas; delete mediumView; delete geo; delete tube; delete cmp; delete sensor; delete drift; //otherwise these are held in memory and cause a killed: 9 crash. delete gas; delete mediumView; delete geo; delete tube; delete cmp; delete sensor; delete drift; //otherwise these are held in memory and cause a killed: 9 crash.
return energy_excess; return energy_excess;
......
NPSimulation/Detectors/TACTIC/TACTICScorer.cc
View file @ 6b557cd8
...@@ -3,11 +3,14 @@ ...@@ -3,11 +3,14 @@
#include "TACTICScorer.hh" #include "TACTICScorer.hh"
#include "G4UnitsTable.hh" #include "G4UnitsTable.hh"
#include "G4RunManager.hh" #include "G4RunManager.hh"
#include "TACTIC.hh"
#ifdef USE_Garfield #ifdef USE_Garfield
#include "GARFDRIFT.h" #include "GARFDRIFT.h"
#endif #endif
double excess;
using namespace TACTICScorer; using namespace TACTICScorer;
Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int NumberOfSegments, G4int depth, G4double p0, G4double p1, G4double p2, G4double p3) //what do level and depth do? Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int NumberOfSegments, G4int depth, G4double p0, G4double p1, G4double p2, G4double p3) //what do level and depth do?
...@@ -30,6 +33,7 @@ Gas_Scorer::~Gas_Scorer(){} ...@@ -30,6 +33,7 @@ Gas_Scorer::~Gas_Scorer(){}
G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
G4double* Infos = new G4double[15]; G4double* Infos = new G4double[15];
//bool first_step = true;
m_Position = aStep->GetPreStepPoint()->GetPosition(); m_Position = aStep->GetPreStepPoint()->GetPosition();
Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(); Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
...@@ -44,6 +48,7 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ ...@@ -44,6 +48,7 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength ) + 1; //Pad number m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength ) + 1; //Pad number
Infos[6] = m_SegmentNumber; Infos[6] = m_SegmentNumber;
//prepad = Infos[6];
Infos[7] = m_Position.z(); Infos[7] = m_Position.z();
Infos[8] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R Infos[8] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R
Infos[9] = aStep->GetTrack()->GetVertexPosition()[2]; Infos[9] = aStep->GetTrack()->GetVertexPosition()[2];
...@@ -53,6 +58,10 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ ...@@ -53,6 +58,10 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
Infos[12] = aStep->GetTrack()->GetTrackLength(); Infos[12] = aStep->GetTrack()->GetTrackLength();
Infos[13] = m_p0 + m_p1*Infos[8] + m_p2*Infos[8]*Infos[8] + m_p3*Infos[8]*Infos[8]*Infos[8]; Infos[13] = m_p0 + m_p1*Infos[8] + m_p2*Infos[8]*Infos[8] + m_p3*Infos[8]*Infos[8]*Infos[8];
//Infos[14] = excess;
G4ThreeVector delta_Position = aStep->GetDeltaPosition();
m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6; m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6;
...@@ -61,31 +70,27 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ ...@@ -61,31 +70,27 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
return 0; return 0;
} }
if(aStep->IsFirstStepInVolume() == true) excess = 0.;
map<G4int, G4double**>::iterator it; map<G4int, G4double**>::iterator it;
it= EvtMap->GetMap()->find(m_Index); it= EvtMap->GetMap()->find(m_Index);
if(it!=EvtMap->GetMap()->end()){ if(it!=EvtMap->GetMap()->end()){
G4double* dummy = *(it->second); G4double* dummy = *(it->second);
if(Infos[1]==dummy[1]) Infos[5]+=dummy[5]; //accumulate ionisation energy deposit to get total accross pad if(Infos[1]==dummy[1]) Infos[5]+=dummy[5]; //accumulate ionisation energy deposit to get total accross pad
delete dummy;
}
#ifdef USE_Garfield #ifdef USE_Garfield
G4ThreeVector delta_Position = aStep->GetDeltaPosition();
G4double excess;
if(aStep->IsFirstStepInVolume() == 1) excess = 0.;
if(aStep->IsFirstStepInVolume() == 0) excess = dummy[14]/eV;
if(excess > 1.e06) excess = 0.; //If dummy[12] is not defined returns random value > 1.e06 ? Infos[14] = GARFDRIFT(((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV+excess), Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0], (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV)*eV;
/*
file.open("excess_test.dat",std::ios::app);
file << Infos[6] << "\t" << "\t" << aStep->IsFirstStepInVolume() << "\t" << excess << "\t" << Infos[14]/eV << "\t" << (int)((((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV+excess) / 41.1)*0.01) << "\t" << Infos[8] << endl;
file.close();
*/
excess = Infos[14]/eV;
if(Infos[8] > 12.) Infos[14] = GARFDRIFT((Infos[5]/eV+excess), Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0])*eV;
#endif #endif
if(Infos[8] < 12.) Infos[14] = 0;
if(Infos[14]>0.) cout << "Infos[14] " << Infos[14]/eV << endl;
delete dummy;
}
EvtMap->set(m_Index, Infos); EvtMap->set(m_Index, Infos);
......
NPSimulation/Detectors/TACTIC/TACTICScorer.hh
View file @ 6b557cd8
...@@ -36,6 +36,7 @@ namespace TACTICScorer { ...@@ -36,6 +36,7 @@ namespace TACTICScorer {
void DrawAll(); void DrawAll();
void PrintAll(); void PrintAll();
private: // Geometry of the detector private: // Geometry of the detector
G4double m_ScorerLength; G4double m_ScorerLength;
G4int m_NumberOfSegments; G4int m_NumberOfSegments;
...@@ -46,6 +47,9 @@ namespace TACTICScorer { ...@@ -46,6 +47,9 @@ namespace TACTICScorer {
G4double m_p1; G4double m_p1;
G4double m_p2; G4double m_p2;
G4double m_p3; G4double m_p3;
ofstream file;
// G4double excess;
private: // inherited from G4VPrimitiveScorer private: // inherited from G4VPrimitiveScorer
G4int HCID; G4int HCID;
......
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