Skip to content
Snippets Groups Projects
Commit af122f3e authored by lheitz's avatar lheitz
Browse files

Add 3D possibility for Isotropic, modified DSSDScorer

parent 3c4da222
No related branches found
No related tags found
2 merge requests!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...,!23Np tool.2.dev
Hide whitespace changes
Inline Side-by-side
NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+ 37
− 13
View file @ af122f3e
......@@ -65,7 +65,9 @@ EventGeneratorIsotropic::SourceParameters::SourceParameters(){
m_z0 = 0 ;
m_SigmaX = 0 ;
m_SigmaY = 0 ;
m_SigmaZ = 0 ;
m_particle = NULL;
m_direction = 'z' ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -73,7 +75,7 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Isotropic");
m_Parameters.reserve(blocks.size());
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "\033[1;35m//// Isotropic reaction found " << endl;
cout << endl << "\033[1;35m//// Isotropic reaction found " << endl;
vector<string> token = {"EnergyLow","EnergyHigh","HalfOpenAngleMin","HalfOpenAngleMax","x0","y0","z0","Particle"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
......@@ -89,6 +91,8 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
vector<string> file = blocks[i]->GetVectorString("EnergyDistributionHist");
m_EnergyDistributionHist = NPL::Read1DProfile(file[0],file[1]);
}
if(blocks[i]->HasToken("Direction"))
it->m_direction=blocks[i]->GetString("Direction");
it->m_HalfOpenAngleMin =blocks[i]->GetDouble("HalfOpenAngleMin","deg");
it->m_HalfOpenAngleMax =blocks[i]->GetDouble("HalfOpenAngleMax","deg");
it->m_x0 =blocks[i]->GetDouble("x0","mm");
......@@ -96,7 +100,7 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
it->m_z0 =blocks[i]->GetDouble("z0","mm");
vector<string> particleName =blocks[i]->GetVectorString("Particle");
for(unsigned int j = 0 ; j < particleName.size() ; j++){
if(particleName[j]=="proton"){ it->m_particleName.push_back("1H") ;}
if(particleName[j]=="proton"){ it->m_particleName.push_back("1H") ;}
else if(particleName[j]=="deuton"){ it->m_particleName.push_back("2H") ; }
else if(particleName[j]=="triton"){ it->m_particleName.push_back("3H") ; }
else if(particleName[j]=="3He" || particleName[j]=="He3") { it->m_particleName.push_back("3He") ; }
......@@ -116,11 +120,13 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
it->m_SigmaX=blocks[i]->GetDouble("SigmaX","mm");
if(blocks[i]->HasToken("SigmaY"))
it->m_SigmaY=blocks[i]->GetDouble("SigmaY","mm");
if(blocks[i]->HasToken("SigmaZ"))
it->m_SigmaZ=blocks[i]->GetDouble("SigmaZ","mm");
if(blocks[i]->HasToken("Multiplicity"))
it->m_Multiplicty=blocks[i]->GetVectorInt("Multiplicity");
}
else{
cout << "ERROR: check your input file formatting \033[0m" << endl;
cout << "ERROR: check your input file formatting \033[0m" << endl;
exit(1);
}
}
......@@ -138,7 +144,7 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
fEnergyDist = new TF1("fDist", par.m_EnergyDistribution, par.m_EnergyLow, par.m_EnergyHigh);
}
}
}
}
}
......@@ -181,19 +187,37 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
particle_energy = fEnergyDist->GetRandom();
}
G4double x0, y0, z0;
G4double momentum_x, momentum_y, momentum_z;
x0 = RandGauss::shoot(par.m_x0,par.m_SigmaX);
y0 = RandGauss::shoot(par.m_y0,par.m_SigmaY);
z0 = RandGauss::shoot(par.m_z0,par.m_SigmaZ);
if(par.m_direction == 'z')
{
// Direction of particle, energy and laboratory angle
momentum_x = sin(theta) * cos(phi) ;
momentum_y = sin(theta) * sin(phi) ;
momentum_z = cos(theta) ;
// Direction of particle, energy and laboratory angle
G4double momentum_x = sin(theta) * cos(phi) ;
G4double momentum_y = sin(theta) * sin(phi) ;
G4double momentum_z = cos(theta) ;
G4double x0 = RandGauss::shoot(par.m_x0,par.m_SigmaX);
G4double y0 = RandGauss::shoot(par.m_y0,par.m_SigmaY);
Particle particle(par.m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, par.m_z0));
}
else if(par.m_direction == 'y')
{
// Direction of particle, energy and laboratory angle
momentum_z = sin(theta) * cos(phi) ;
momentum_x = sin(theta) * sin(phi) ;
momentum_y = cos(theta) ;
}
else // = 'x'
{
// Direction of particle, energy and laboratory angle
momentum_y = sin(theta) * cos(phi) ;
momentum_z = sin(theta) * sin(phi) ;
momentum_x = cos(theta) ;
}
Particle particle(par.m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, z0));
m_ParticleStack->AddParticleToStack(particle);
}
}
......
NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
+ 5
− 3
View file @ af122f3e
......@@ -45,19 +45,19 @@ class EventGeneratorIsotropic : public NPS::VEventGenerator{
public: // Constructor and destructor
EventGeneratorIsotropic() ;
virtual ~EventGeneratorIsotropic();
public: // Inherit from VEventGenerator Class
void ReadConfiguration(NPL::InputParser) ;
void GenerateEvent(G4Event*) ;
void InitializeRootOutput() ;
private: // Source parameter from input file
G4int event_ID;
struct SourceParameters {
SourceParameters() ;
G4double m_EnergyLow ; // Lower limit of energy range
G4double m_EnergyHigh ; // Upper limit of energy range
TString m_EnergyDistribution;
TString m_EnergyDistribution;
G4double m_HalfOpenAngleMin ; // Min Half open angle of the source
G4double m_HalfOpenAngleMax ; // Max Half open angle of the source
G4double m_x0 ; // Vertex Position X
......@@ -65,6 +65,8 @@ private: // Source parameter from input file
G4double m_z0 ; // Vertex Position Z
G4double m_SigmaX ;
G4double m_SigmaY ;
G4double m_SigmaZ ;
TString m_direction ;
G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
vector<G4double> m_ExcitationEnergy ; // Excitation energy of the emitted particle
vector<string> m_particleName ;
......
NPSimulation/Scorers/DSSDScorers.cc
+ 114
− 26
View file @ af122f3e
......@@ -22,6 +22,21 @@
#include "DSSDScorers.hh"
#include "G4UnitsTable.hh"
using namespace DSSDScorers;
/*
vector<DSSDData>::iterator DSSDDataVector::find(const unsigned int& index, const double& time , const double TimeThreshold = 0) {
for (vector<DSSDData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) {
G4bool checkIndex =(*it).GetIndex() == index;
G4bool checkTime = 1;
if (TimeThreshold>0)
{
checkTime = std::abs((*it).GetTime() - time) < TimeThreshold;
};
if (checkIndex & checkTime)
return it; //Autre possibilité : rajouter Temps en argument, + check *it.GetTime() == Temps)
}
return m_Data.end();
}
*/
vector<DSSDData>::iterator DSSDDataVector::find(const unsigned int& index) {
for (vector<DSSDData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) {
if ((*it).GetIndex() == index)
......@@ -29,10 +44,31 @@ vector<DSSDData>::iterator DSSDDataVector::find(const unsigned int& index) {
}
return m_Data.end();
}
/*
vector<DSSDData>::iterator DSSDDataVector::find(const unsigned int& index) {
for (vector<DSSDData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) {
G4bool checkIndex =(*it).GetIndex() == index;
if (checkIndex)
return it; //Autre possibilité : rajouter Temps en argument, + check *it.GetTime() == Temps)
}
return m_Data.end();
}
*/
/*
vector<DSSDData>::iterator DSSDDataVector::findTime(const unsigned int& index, const double& time, const double thresh = 1000* ns) {
for (vector<DSSDData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) {
G4bool checkIndex =(*it).GetIndex() == index;
G4bool checkTime = std::abs((*it).GetTime() - time) < thresh;
if ( checkIndex && checkTime)
return it; //Autre possibilité : rajouter Temps en argument, + check *it.GetTime() == Temps)
}
return m_Data.end();
}
*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PS_Images::PS_Images(G4String name, string imageFront, string imageBack, double scalingFront, double scalingBack,
double centerOffsetX, double centerOffsetY, unsigned int ignoreValue, G4int depth)
double centerOffsetX, double centerOffsetY, unsigned int ignoreValue, G4int depth,bool saveAll)
: G4VPrimitiveScorer(name, depth) {
m_ImageFront = new NPL::Image(imageFront, scalingFront, scalingFront);
m_ImageBack = new NPL::Image(imageBack, scalingBack, scalingBack);
......@@ -42,6 +78,7 @@ PS_Images::PS_Images(G4String name, string imageFront, string imageBack, double
m_CenterOffsetY = centerOffsetY;
m_IgnoreValue = ignoreValue;
m_Level = depth;
m_SaveAll = saveAll;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -49,44 +86,68 @@ G4bool PS_Images::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
// contain Energy Time, DetNbr, PixelFront and PixelBack
t_Energy = aStep->GetTotalEnergyDeposit();
G4Track* track = aStep->GetTrack();
// Get the track ID
t_trackID = track->GetTrackID();
if(t_Energy == 0)
{
return FALSE;
}
//std::cout << "In scrorer, energy = " << t_Energy << std::endl;
t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
t_DetectorNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
t_Position = aStep->GetPreStepPoint()->GetPosition();
// transforming the position to the local volume
t_Position =
aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position);
t_PixelFront = m_ImageFront->GetPixelAtCoordinate(t_Position.x(), t_Position.y());
t_PixelBack = m_ImageBack->GetPixelAtCoordinate(t_Position.x(), t_Position.y());
t_PixelFront = m_ImageFront->GetPixelAtCoordinate(t_Position.x()-m_CenterOffsetX, t_Position.y()-m_CenterOffsetY);
t_PixelBack = m_ImageBack->GetPixelAtCoordinate(t_Position.x()-m_CenterOffsetX, t_Position.y()-m_CenterOffsetY);
// If front and back are in inactive part of the wafer,
// nothing is added to the unordered_map
if (t_PixelFront == m_IgnoreValue && t_PixelBack == m_IgnoreValue)
return FALSE;
// Check if the particle has interact before, if yes, add up the energies.
vector<DSSDData>::iterator it;
it = m_HitFront.find(DSSDData::CalculateIndex(t_PixelFront, t_DetectorNbr));
if (it != m_HitFront.end()) {
it->Add(t_Energy);
if(m_SaveAll) // In that case, save every energy deposit
{
m_HitFront.Set(t_Energy, t_Time, t_PixelFront, t_DetectorNbr);
m_TrackId.push_back(t_trackID);
}
else {
m_HitFront.Set(t_Energy, t_Time, t_PixelFront, t_DetectorNbr);
else //Add up energies if necessary
{
if (it != m_HitFront.end()) {// Check if the particle has interact before, if yes, add up the energies.
it->Add(t_Energy);
}
else {
m_HitFront.Set(t_Energy, t_Time, t_PixelFront, t_DetectorNbr);
}
}
// Check if the particle has interact before, if yes, add up the energies.
it = m_HitBack.find(DSSDData::CalculateIndex(t_PixelBack, t_DetectorNbr));
if (it != m_HitBack.end()) {
it->Add(t_Energy);
if(m_SaveAll) // Save every step
{
m_HitBack.Set(t_Energy, t_Time, t_PixelBack, t_DetectorNbr);
}
else {
m_HitBack.Set(t_Energy, t_Time, t_PixelBack, t_DetectorNbr);
else // Add up energies if conditions are met
{
if (it != m_HitBack.end()) {// Check if the particle has interact before, if yes, add up the energies.
it->Add(t_Energy);
}
else {
m_HitBack.Set(t_Energy, t_Time, t_PixelBack, t_DetectorNbr);
}
}
return TRUE;
}
......@@ -100,6 +161,7 @@ void PS_Images::EndOfEvent(G4HCofThisEvent*) {}
void PS_Images::clear() {
m_HitFront.clear();
m_HitBack.clear();
m_TrackId.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_Images::GetARGBFront(unsigned int& i, unsigned int& a, unsigned int& r, unsigned int& g, unsigned int& b) {
......@@ -120,15 +182,20 @@ void PS_Images::GetARGBBack(unsigned int& i, unsigned int& a, unsigned int& r, u
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PS_Rectangle::PS_Rectangle(G4String name, G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth,
G4int NumberOfStripLength, G4int NumberOfStripWidth, G4int depth, G4String axis)
G4int NumberOfStripLength, G4int NumberOfStripWidth, G4int depth, G4String axis,
G4double TimeThreshold, G4double InterStripLength, G4double InterStripWidth)
: G4VPrimitiveScorer(name, depth) {
m_StripPlaneLength = StripPlaneLength;
m_StripPlaneWidth = StripPlaneWidth;
m_NumberOfStripLength = NumberOfStripLength;
m_NumberOfStripWidth = NumberOfStripWidth;
m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength;
m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength ;
m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
m_Level = Level;
m_TimeThreshold = TimeThreshold;
m_InterStripLength = InterStripLength;
m_InterStripWidth = InterStripWidth;
if (axis == "xy")
m_Axis = ps_xy;
else if (axis == "yz")
......@@ -142,7 +209,6 @@ PS_Rectangle::~PS_Rectangle() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool PS_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
// contain Energy Time, DetNbr, StripFront and StripBack
t_Energy = aStep->GetTotalEnergyDeposit();
t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
......@@ -152,11 +218,21 @@ G4bool PS_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
t_Position =
aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position);
G4bool checkInterStrip_x = 1;
G4bool checkInterStrip_y = 1;
if (m_Axis == ps_xy) {
t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
t_StripWidthNumber = (int)((t_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
}
G4double x0 =t_Position.x() + m_StripPlaneLength / 2.;
G4double middle_x = (2*t_StripLengthNumber-1) * m_StripPitchLength/2;
checkInterStrip_x =std::abs(x0-middle_x)<(m_StripPitchLength - m_InterStripLength)/2 ;
G4double y0 = t_Position.y() + m_StripPlaneWidth/2.;
G4double middle_y = (2*t_StripWidthNumber-1) * m_StripPitchWidth/2;
checkInterStrip_y =std::abs(y0-middle_y)<(m_StripPitchWidth - m_InterStripWidth)/2 ;
}
else if (m_Axis == ps_yz) {
t_StripLengthNumber = (int)((t_Position.y() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
......@@ -165,30 +241,42 @@ G4bool PS_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
}
// Rare case where particle is close to edge of silicon plan
if (t_StripLengthNumber > m_NumberOfStripLength)
t_StripLengthNumber = m_NumberOfStripLength;
if (t_StripWidthNumber > m_NumberOfStripWidth)
t_StripWidthNumber = m_NumberOfStripWidth;
// Check if the particle has interact before, if yes, add up the energies.
vector<DSSDData>::iterator it;
// Length
//it = m_HitLength.find(DSSDData::CalculateIndex(t_StripLengthNumber, t_DetectorNumber),t_Time,m_TimeThreshold);
it = m_HitLength.find(DSSDData::CalculateIndex(t_StripLengthNumber, t_DetectorNumber));
//if(checkInterStrip_x!=0)
//{
/*
if (it != m_HitLength.end()) {
it->Add(t_Energy);
}
else
*/
m_HitLength.Set(t_Energy, t_Time, t_StripLengthNumber, t_DetectorNumber);
// Width
//it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber),t_Time,m_TimeThreshold);
it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber));
//if(checkInterStrip_y!=0)
//{
// Width
//it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber));
/*
if (it != m_HitWidth.end()) {
it->Add(t_Energy);
}
else
*/
m_HitWidth.Set(t_Energy, t_Time, t_StripWidthNumber, t_DetectorNumber);
//}
return TRUE;
}
......@@ -356,7 +444,7 @@ G4bool PS_Resistive::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
else if (m_dir == "z")
P = (t_Position.z()) / (0.5 * m_StripPlaneLength);
else {
std::cout << "Error : Resistive strip DSSD scorer direction incorrect, should be x,y, or z " << std::endl;
G4cout << "Error : Resistive strip DSSD scorer direction incorrect, should be x,y, or z " << G4endl;
exit(1);
}
......
NPSimulation/Scorers/DSSDScorers.hh
+ 13
− 2
View file @ af122f3e
......@@ -87,7 +87,10 @@ namespace DSSDScorers {
vector<DSSDData> m_Data;
public:
//vector<DSSDData>::iterator find(const unsigned int& index, const double& time, const double TimeThreshold);
vector<DSSDData>::iterator find(const unsigned int& index);
//vector<DSSDData>::iterator findTime(const unsigned int& index, const double& time, const double thresh) ;
//const unsigned int& index, const double& time, const double TimeThreshold
inline void clear() { m_Data.clear(); };
inline vector<DSSDData>::iterator end() { return m_Data.end(); };
inline vector<DSSDData>::iterator begin() { return m_Data.begin(); };
......@@ -104,7 +107,7 @@ namespace DSSDScorers {
public: // with description
PS_Images(G4String name, string imageFront, string imageBack, double scalingFront, double scalingBack,
double centerOffsetX, double centerOffsetY, unsigned int ignoreValue, G4int depth = 0);
double centerOffsetX, double centerOffsetY, unsigned int ignoreValue, G4int depth = 0, bool saveAll = false);
~PS_Images(){};
protected: // with description
......@@ -125,6 +128,7 @@ namespace DSSDScorers {
double m_CenterOffsetX;
double m_CenterOffsetY;
unsigned int m_IgnoreValue;
bool m_SaveAll;
// Level at which to find the copy number linked to the detector number
G4int m_Level;
......@@ -132,11 +136,13 @@ namespace DSSDScorers {
private: // inherited from G4VPrimitiveScorer
DSSDDataVector m_HitFront;
DSSDDataVector m_HitBack;
std::vector<int> m_TrackId;
private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
G4ThreeVector t_Position;
double t_Energy;
double t_Time;
int t_trackID;
unsigned int t_DetectorNbr;
unsigned int t_PixelFront;
unsigned int t_PixelBack;
......@@ -152,6 +158,7 @@ namespace DSSDScorers {
inline unsigned int GetDetectorBack(const unsigned int& i) { return m_HitBack[i]->GetDetector(); };
inline double GetEnergyBack(const unsigned int& i) { return m_HitBack[i]->GetEnergy(); };
inline double GetTimeBack(const unsigned int& i) { return m_HitBack[i]->GetTime(); };
inline int GetTrackId(const unsigned int& i){return m_TrackId[i]; };
void GetARGBFront(unsigned int& i, unsigned int& a, unsigned int& r, unsigned int& g, unsigned int& b);
void GetARGBBack(unsigned int& i, unsigned int& a, unsigned int& r, unsigned int& g, unsigned int& b);
......@@ -162,7 +169,8 @@ namespace DSSDScorers {
public: // with description
PS_Rectangle(G4String name, G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth,
G4int NumberOfStripLength, G4int NumberOfStripWidth, G4int depth = 0, G4String axis = "xy");
G4int NumberOfStripLength, G4int NumberOfStripWidth, G4int depth = 0, G4String axis = "xy",
G4double TimeThreshold = 0, G4double InterStripLength = 0,G4double InterStripWidth = 0);
~PS_Rectangle();
private:
......@@ -186,6 +194,9 @@ namespace DSSDScorers {
unsigned int m_NumberOfStripWidth;
double m_StripPitchLength;
double m_StripPitchWidth;
double m_TimeThreshold;
double m_InterStripLength;
double m_InterStripWidth;
// Level at which to find the copy number linked to the detector number
int m_Level;
......
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