Skip to content
Snippets Groups Projects
Commit e9e9ce6a authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* importing change in nebula simulation from Louis Lemair

parent e14a5361
No related branches found
No related tags found
No related merge requests found
Pipeline #232422 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPSimulation/Detectors/Nebula/Nebula.cc
+ 140
− 46
View file @ e9e9ce6a
......@@ -33,6 +33,7 @@
//G4 various object
#include "G4Material.hh"
#include "G4MaterialPropertiesTable.hh"
#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
#include "G4VisAttributes.hh"
......@@ -42,6 +43,7 @@
#include "Nebula.hh"
#include "PlasticBar.hh"
#include "InteractionScorers.hh"
#include "ProcessScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
#include "NPSDetectorFactory.hh"
......@@ -56,19 +58,24 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Nebula_NS{
// Energy and time Resolution
const double EnergyThreshold = 0.1*MeV;
const double ResoTime = 4.5*ns ;
const double ResoEnergy = 1.0*MeV ;
const double ModuleWidth = 120*mm ;
const double ModuleLength = 120*mm ;
const double ModuleHeight = 1800*mm ;
const double InterModule = 1*mm ;
const double VetoWidth = 320*mm ;
const double VetoLength = 10*mm ;
const double VetoHeight = 1900*mm ;
const double InterVeto = 1*mm ;
const int VetoPerWall = 12;
const double WallToVeto = 10*cm;
const double LightThreshold = 0.1*MeV;
const double ResoTime = 0.75/2.355*ns; //0.75
const double ResoEnergy = 0.1/2.355*MeV;
const double ResoLight = 0.1/2.355*MeV;
const double ResoPosition = 1.0*um; //1.0
const double ModuleWidth = 120*mm ;
const double ModuleLength = 120*mm ;
const double ModuleHeight = 1800*mm ;
const double InterModule = 1*mm ;
const double VetoWidth = 320*mm ;
const double VetoLength = 10*mm ;
const double VetoHeight = 1900*mm ;
const double InterVeto = 1*mm ;
const int VetoPerWall = 12;
const int VetoPerExpand = 6;
const double WallToVeto = 10*cm;
const double MaterialIndex = 1.58;
const double Attenuation = 6680*mm;
const string Material = "BC400";
}
......@@ -85,7 +92,8 @@ Nebula::Nebula(){
// RGB Color + Transparency
m_VisModule = new G4VisAttributes(G4Colour(0.3, 0.3, 1, 1));
m_VisModule = new G4VisAttributes(G4Colour(0.263, 0.682, 0.639, 1));
//m_VisModule = new G4VisAttributes(G4Colour(0.145, 0.384, 0.596, 1));
m_VisVeto = new G4VisAttributes(G4Colour(0.4, 0.4, 0.4, 0.2));
m_VisPMT = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1, 1));
m_VisFrame = new G4VisAttributes(G4Colour(0, 0.3, 1, 0.5));
......@@ -131,7 +139,6 @@ G4LogicalVolume* Nebula::BuildVeto(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -162,6 +169,9 @@ void Nebula::ReadConfiguration(NPL::InputParser parser){
exit(1);
}
}
std::for_each(m_NbrModule.begin(), m_NbrModule.end(), [&] (int n) {
m_TotalModule += n;
});
}
......@@ -180,30 +190,41 @@ void Nebula::ConstructDetector(G4LogicalVolume* world){
new G4PVPlacement(G4Transform3D(*Rot,m_Pos[i]+Offset),
BuildModule(),
"NebulaModule",world,false,nbrM++);
}
}
if(m_HasVeto[i]){
for (unsigned short m = 0 ; m < Nebula_NS::VetoPerWall ; m++) {
G4RotationMatrix* Rot = new G4RotationMatrix();
double offset = (Nebula_NS::VetoWidth+Nebula_NS::InterVeto)*(-Nebula_NS::VetoPerWall*0.5+m)+Nebula_NS::VetoWidth*0.5;
G4ThreeVector Offset(offset,0,-Nebula_NS::WallToVeto);
new G4PVPlacement(G4Transform3D(*Rot,m_Pos[i]+Offset),
BuildVeto(),
"NebulaVeto",world,false,nbrV++);
}
if(m_NbrModule[i] > 15){
for (unsigned short m = 0 ; m < Nebula_NS::VetoPerWall ; m++) {
G4RotationMatrix* Rot = new G4RotationMatrix();
double offset = (Nebula_NS::VetoWidth+Nebula_NS::InterVeto)*(-Nebula_NS::VetoPerWall*0.5+m)+Nebula_NS::VetoWidth*0.5;
G4ThreeVector Offset(offset,0,-Nebula_NS::WallToVeto);
new G4PVPlacement(G4Transform3D(*Rot,m_Pos[i]+Offset),
BuildVeto(),
"NebulaVeto",world,false,nbrV++);
}
}
}
else{
for (unsigned short m = 0 ; m < Nebula_NS::VetoPerExpand ; m++) {
G4RotationMatrix* Rot = new G4RotationMatrix();
double offset = (Nebula_NS::VetoWidth+Nebula_NS::InterVeto)*(-Nebula_NS::VetoPerExpand*0.5+m)+Nebula_NS::VetoWidth*0.5;
G4ThreeVector Offset(offset,0,-Nebula_NS::WallToVeto);
new G4PVPlacement(G4Transform3D(*Rot,m_Pos[i]+Offset),
BuildVeto(),
"NebulaVeto",world,false,nbrV++);
}
}
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
// Called After DetectorConstruction::AddDetector Method
void Nebula::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
if(!pTree->FindBranch("Nebula")){
pTree->Branch("Nebula", "TNebulaData", &m_Event) ;
pTree->Branch("Nebula", "TNebulaData",&m_Event) ;
}
pTree->SetBranchAddress("Nebula", &m_Event) ;
}
......@@ -215,22 +236,91 @@ void Nebula::ReadSensitive(const G4Event* ){
m_Event->Clear();
///////////
// Module scorer
PlasticBar::PS_PlasticBar* Scorer= (PlasticBar::PS_PlasticBar*) m_ModuleScorer->GetPrimitive(0);
unsigned int size = Scorer->GetMult();
for(unsigned int i = 0 ; i < size ; i++){
vector<unsigned int> level = Scorer->GetLevel(i);
double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Nebula_NS::ResoEnergy);
if(Energy>Nebula_NS::EnergyThreshold){
double Time = RandGauss::shoot(Scorer->GetTime(i),Nebula_NS::ResoTime);
// PlasticBar scorer
PlasticBar::PS_PlasticBar* PlasticScorer_Module = (PlasticBar::PS_PlasticBar*) m_ModuleScorer->GetPrimitive(0);
PlasticBar::PS_PlasticBar* PlasticScorer_Veto = (PlasticBar::PS_PlasticBar*) m_VetoScorer->GetPrimitive(0);
// Should we put a ProcessScorer here to get the info if the particle is first neutron and give it to NebulaData ?
double Time_up, Time_down;
double Energy_tmp, Light_tmp;
//////////// TRIAL TO GET THE OPTICAL INDEX FROM MATERIAL PROPERTIES /////////////
//Trying to get Optical Index from Material directly
//const G4Material* aMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(Nebula_NS::Material);
//G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
//if(!aMaterialPropertiesTable->PropertyExists("RINDEX")){
// MaterialIndex = !aMaterialPropertiesTable->GetConstProperty("RINDEX");
//}
//else{
// MaterialIndex = 0;
//}
//cout << MaterialManager::getInstance()->GetMaterialFromLibrary(Nebula_NS::Material)->GetMaterialPropertiesTable()->GetMaterialPropertyNames()[0] << endl;
//////////////////////////////////////////////////////////////////////////////////
///////////////////////////////// MODULE SCORER //////////////////////////////////
unsigned int ModuleHits_size = PlasticScorer_Module->GetMult();
for(unsigned int i = 0 ; i < ModuleHits_size ; i++){
vector<unsigned int> level = PlasticScorer_Module->GetLevel(i);
Energy_tmp = PlasticScorer_Module->GetEnergy(i);
Light_tmp = PlasticScorer_Module->GetLight(i);
Energy = RandGauss::shoot(Energy_tmp, Energy_tmp*Nebula_NS::ResoEnergy);
Light = RandGauss::shoot(Light_tmp, Light_tmp*Nebula_NS::ResoLight);
if(Light>Nebula_NS::LightThreshold){
double Time = RandGauss::shoot(PlasticScorer_Module->GetTime(i),Nebula_NS::ResoTime);
//cout << "Time is " << Time << endl;
double Position = RandGauss::shoot(PlasticScorer_Module->GetPosition(i),Nebula_NS::ResoPosition);
//cout << "Position is " << Position << endl;
int DetectorNbr = level[0];
m_Event->SetChargeUp(DetectorNbr,Energy/2);
m_Event->SetChargeDown(DetectorNbr,Energy/2);
m_Event->SetTimeUp(DetectorNbr,Time/2);
m_Event->SetTimeDown(DetectorNbr,Time/2);
//cout << "Detector ID: " << DetectorNbr << endl;
m_Event->SetChargeUp(DetectorNbr,Light*exp(-(Nebula_NS::ModuleHeight/2-Position)/Nebula_NS::Attenuation));
m_Event->SetChargeDown(DetectorNbr,Light*exp(-(Nebula_NS::ModuleHeight/2+Position)/Nebula_NS::Attenuation));
Time_up = (Nebula_NS::ModuleHeight/2-Position)/(c_light/Nebula_NS::MaterialIndex) + Time;
//cout << "Time_up is " << Time_up << endl;
m_Event->SetTimeUp(DetectorNbr,Time_up);
Time_down = (Nebula_NS::ModuleHeight/2+Position)/(c_light/Nebula_NS::MaterialIndex) + Time;
//cout << "Time_down is " << Time_down << endl;
m_Event->SetTimeDown(DetectorNbr,Time_down);
}
}
//cout << endl;
///////////////////////////////// VETO SCORER //////////////////////////////////
unsigned int VetoHits_size = PlasticScorer_Veto->GetMult();
for(unsigned int i = 0 ; i < VetoHits_size ; i++){
vector<unsigned int> level = PlasticScorer_Veto->GetLevel(i);
Energy_tmp = PlasticScorer_Veto->GetEnergy(i);
Light_tmp = PlasticScorer_Veto->GetLight(i);
Energy = RandGauss::shoot(Energy_tmp, Energy_tmp*Nebula_NS::ResoEnergy);
Light = RandGauss::shoot(Light_tmp, Light_tmp*Nebula_NS::ResoLight);
if(Light>Nebula_NS::LightThreshold){
double Time = RandGauss::shoot(PlasticScorer_Veto->GetTime(i),Nebula_NS::ResoTime);
//cout << "Time is " << Time << endl;
double Position = RandGauss::shoot(PlasticScorer_Veto->GetPosition(i),Nebula_NS::ResoPosition);
//cout << "Position is " << Position << endl;
int DetectorNbr = level[0] + m_TotalModule;
//cout << "Veto ID: " << DetectorNbr << endl;
m_Event->SetChargeUp(DetectorNbr,Light*exp(-(Nebula_NS::VetoHeight/2-Position)/Nebula_NS::Attenuation));
m_Event->SetChargeDown(DetectorNbr,Light*exp(-(Nebula_NS::VetoHeight/2+Position)/Nebula_NS::Attenuation));
Time_up = (Nebula_NS::VetoHeight/2-Position)/(c_light/Nebula_NS::MaterialIndex) + Time;
//cout << "Time_up is " << Time_up << endl;
m_Event->SetTimeUp(DetectorNbr,Time_up);
Time_down = (Nebula_NS::VetoHeight/2+Position)/(c_light/Nebula_NS::MaterialIndex) + Time;
//cout << "Time_down is " << Time_down << endl;
m_Event->SetTimeDown(DetectorNbr,Time_down);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -247,19 +337,23 @@ void Nebula::InitializeScorers() {
// Otherwise the scorer is initialise
// Module
vector<int> level; level.push_back(0);
G4VPrimitiveScorer* ModulePlasticBar= new PlasticBar::PS_PlasticBar("ModulePlasticBar",level, 0) ;
G4VPrimitiveScorer* ModuleInteraction= new InteractionScorers::PS_Interactions("ModuleInteraction",ms_InterCoord, 0) ;
G4VPrimitiveScorer* ModulePlasticBar= new PlasticBar::PS_PlasticBar("ModulePlasticBar",level, 0);
G4VPrimitiveScorer* ModuleInteraction= new InteractionScorers::PS_Interactions("ModuleInteraction",ms_InterCoord, 0);
G4VPrimitiveScorer* ModuleProcess= new ProcessScorers::PS_Process("ModuleProcess", 0);
//and register it to the multifunctionnal detector
m_ModuleScorer->RegisterPrimitive(ModulePlasticBar);
m_ModuleScorer->RegisterPrimitive(ModuleInteraction);
m_ModuleScorer->RegisterPrimitive(ModuleProcess);
G4SDManager::GetSDMpointer()->AddNewDetector(m_ModuleScorer) ;
// Veto
G4VPrimitiveScorer* VetoPlasticBar= new PlasticBar::PS_PlasticBar("VetoPlasticBar",level, 0) ;
G4VPrimitiveScorer* VetoInteraction= new InteractionScorers::PS_Interactions("VetoInteraction",ms_InterCoord, 0) ;
G4VPrimitiveScorer* VetoPlasticBar= new PlasticBar::PS_PlasticBar("VetoPlasticBar",level, 0);
G4VPrimitiveScorer* VetoInteraction= new InteractionScorers::PS_Interactions("VetoInteraction",ms_InterCoord, 0);
G4VPrimitiveScorer* VetoProcess= new ProcessScorers::PS_Process("ModuleProcess", 0);
//and register it to the multifunctionnal detector
m_VetoScorer->RegisterPrimitive(VetoPlasticBar);
m_VetoScorer->RegisterPrimitive(VetoInteraction);
m_VetoScorer->RegisterPrimitive(VetoProcess);
G4SDManager::GetSDMpointer()->AddNewDetector(m_VetoScorer) ;
......
NPSimulation/Detectors/Nebula/Nebula.hh
+ 3
− 0
View file @ e9e9ce6a
......@@ -58,6 +58,8 @@ class Nebula : public NPS::VDetector{
private:
G4LogicalVolume* m_Module;
G4LogicalVolume* m_Veto;
double Energy;
double Light;
////////////////////////////////////////////////////
////// Inherite from NPS::VDetector class /////////
......@@ -99,6 +101,7 @@ class Nebula : public NPS::VDetector{
// Detector Coordinate
vector<G4ThreeVector> m_Pos;
vector<int> m_NbrModule;
int m_TotalModule = 0;
vector<bool> m_HasVeto;
vector<bool> m_HasFrame;
......
NPSimulation/Scorers/PlasticBar.cc
+ 114
− 27
View file @ e9e9ce6a
......@@ -20,8 +20,11 @@
* PAD detector (any Silicon Detector) *
*****************************************************************************/
#include "PlasticBar.hh"
#include "G4SteppingManager.hh"
#include "G4UnitsTable.hh"
#include "RootOutput.h"
using namespace PlasticBar;
using namespace std;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
unsigned int
......@@ -37,6 +40,53 @@ PlasticBarData::CalculateIndex(const vector<unsigned int>& level) {
return result;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double
PS_PlasticBar::EnergyToLight(double Energy, int Z){
//Get Energy to Light conversion factor depending on the particle doing the ionisation
double a0 = t_Convertor[Z][0];
double a1 = t_Convertor[Z][1];
double a2 = t_Convertor[Z][2];
double a3 = t_Convertor[Z][3];
double Light = a0*Energy - a1*( 1. - std::exp(-a2*std::pow(Energy,a3)) );
return( std::max(Light ,0.) );
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void
PS_PlasticBar::Compute_Light(){
std::vector<PlasticBarData>::iterator it = m_Data.begin();
std::vector<std::array<double,7>>::iterator it_Energy = t_Energy_by_ChargeNumber.begin();
std::vector<std::array<double,7>>::iterator it_Light = t_Light_by_ChargeNumber.begin();
for(; it_Energy < t_Energy_by_ChargeNumber.end() ; it_Energy++, it_Light++, it++){
//cout << "Barre ID : " << (*it).GetIndex() << endl;
double Light = 0, Energy = 0;
double ConvertedLight = 0;
//cout << " Energy Light Z " << endl;
for(int i=0 ; i<7 ; i++){
Energy += (*it_Energy)[i];
if(t_DoConversion){
ConvertedLight = EnergyToLight( (*it_Energy)[i], i);
(*it_Light)[i] = ConvertedLight;
Light += ConvertedLight;
t_TotalLight_by_Z[i] += ConvertedLight;
}
//cout << " " << Energy;
//cout << " " << Light;
//cout << " " << i << endl;
}
//cout << endl;
(*it).SetEnergy(Energy);
if(t_DoConversion)
(*it).SetLight(Light);
else
(*it).SetLight(Energy);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
vector<PlasticBarData>::iterator
PlasticBarDataVector::find(const unsigned int& index) {
......@@ -49,10 +99,14 @@ PlasticBarDataVector::find(const unsigned int& index) {
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PS_PlasticBar::PS_PlasticBar(G4String name, vector<G4int> NestingLevel,
G4int depth)
PS_PlasticBar::PS_PlasticBar(G4String name, vector<G4int> NestingLevel, G4int depth)
: G4VPrimitiveScorer(name, depth) {
m_NestingLevel = NestingLevel;
auto tree = RootOutput::getInstance()->GetTree();
tree->Branch("PlasticBar_Time", &t_Time);
tree->Branch("PlasticBar_Position", &t_Position);
tree->Branch("PlasticBar_TotalEnergy_by_Z", &t_TotalEnergy_by_Z);
tree->Branch("PlasticBar_TotalLight_by_Z", &t_TotalLight_by_Z);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -60,40 +114,73 @@ PS_PlasticBar::~PS_PlasticBar() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool PS_PlasticBar::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
G4String particlename
= aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
cout << particlename << " has left ";
// Contain Energy, Time + as many copy number as nested volume
unsigned int mysize = m_NestingLevel.size();
t_Energy = aStep->GetTotalEnergyDeposit();
cout << t_Energy << " MeV";
t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
cout << " after " << t_Time << " ms" << endl;
t_Level.clear();
for (unsigned int i = 0; i < mysize; i++) {
t_Level.push_back(
aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(
m_NestingLevel[i]));
}
// Check if the particle has interact before, if yes, add up the energies.
vector<PlasticBarData>::iterator it;
it = m_Data.find(PlasticBarData::CalculateIndex(t_Level));
if (it != m_Data.end()) {
it->Add(t_Energy);
} else {
m_Data.Set(t_Energy, t_Time, t_Level);
}
return TRUE;
// Contain Energy, Time + as many copy number as nested volume
unsigned int mysize = m_NestingLevel.size();
// first lets consider a step going inside the bar, the interaction point can be associated to the position of the end of the step (which helps avoiding having all the steps at the entrance of the bar)
G4String processname = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
if(processname!="Transportation"){
G4String particlename = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
int Z = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();
//under the energy cut, neutrons can deposit energy via hadIon. But the conversion cant be done since we dont know the recoil particle. We chose to neglect their effect. (validated because deposited energy is less than 1MeV, which is the PM threshold)
if(particlename!="neutron" && Z<7){
t_Energy = aStep->GetTotalEnergyDeposit();
t_Time = aStep->GetPostStepPoint()->GetGlobalTime();
t_Position = aStep->GetPostStepPoint()->GetPosition()[1];
t_Level.clear();
for (unsigned int i = 0; i < mysize; i++) {
t_Level.push_back(
aStep->GetPostStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i])
);
}
if(t_Time < 150){
// Check if the bar has been hit before, if not, set the first interaction point.
std::vector<PlasticBarData>::iterator it;
it = m_Data.find(PlasticBarData::CalculateIndex(t_Level));
if (it == m_Data.end()) {
m_Data.Set(0.0, 0.0, t_Position, t_Time, t_Level);
AddEntry(); //Adds null entry to t_Energy_by_ChargeNumber
auto lastEnergyEntry = t_Energy_by_ChargeNumber.end() - 1; //iterator
(*lastEnergyEntry)[Z] += t_Energy;
t_TotalEnergy_by_Z[Z] += t_Energy;
}
else {
t_Energy_by_ChargeNumber[it - m_Data.begin()][Z] += t_Energy;
t_TotalEnergy_by_Z[Z] += t_Energy;
}
}
}
//if(it->GetLight()>160){
// cout << "{!!!} Deposited Light : " << it->GetLight() << endl;
//}
//it->Add_Light(t_Light);
//it->Set_Light(EnergyToLight(it->GetEnergy(),Z));
//cout << particlename << " of Z = " << Z << " has left :" << endl;
//cout << t_Energy << " MeV" << endl;
//cout << EnergyToLight(it->GetEnergy(),Z)*MeV << " MeVEE" << endl;
//now, if the next interaction is a Transportation, and there is an energy deposition, it means that we are observing the last step. To avoid having all the positions on the exit of the bar, lets consider the PreStepPosition instead.
}
return TRUE;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_PlasticBar::Initialize(G4HCofThisEvent*) { clear(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_PlasticBar::EndOfEvent(G4HCofThisEvent*) {}
void PS_PlasticBar::EndOfEvent(G4HCofThisEvent*) {
Compute_Light();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_PlasticBar::clear() {
t_Energy_by_ChargeNumber.clear();
t_Light_by_ChargeNumber.clear();
m_Data.clear();
t_Level.clear();
}
......
NPSimulation/Scorers/PlasticBar.hh
+ 93
− 41
View file @ e9e9ce6a
......@@ -29,40 +29,55 @@
//#include "NPSecondaries.hh"
#include <map>
using namespace std;
using namespace CLHEP;
namespace PlasticBar {
// Hold One hit info
class PlasticBarData{
public:
PlasticBarData(const double& Energy,const double& Time,const vector<unsigned int>& Nesting){
PlasticBarData(const double& Energy, const double& Light, const double& Position, const double& Time,const std::vector<unsigned int>& Nesting){
m_Index=CalculateIndex(Nesting);
m_Level=Nesting;
m_Energy=Energy;
m_Light=Light;
m_Position=Position;
m_Time=Time;
};
};
~PlasticBarData(){};
private:
unsigned int m_Index;
vector<unsigned int> m_Level;
std::vector<unsigned int> m_Level;
double m_Energy;
double m_Light;
double m_Position;
double m_Time;
public:
static unsigned int CalculateIndex(const vector<unsigned int>& Nesting);
static unsigned int CalculateIndex(const std::vector<unsigned int>& Nesting);
public:
inline unsigned int GetIndex() const {return m_Index;}
inline vector<unsigned int> GetLevel() const {return m_Level;};
inline unsigned int GetIndex() const {return m_Index;};
inline std::vector<unsigned int> GetLevel() const {return m_Level;};
inline double GetEnergy() const {return m_Energy;};
inline double GetLight() const {return m_Light;};
inline double GetPosition() const {return m_Position;};
inline double GetTime() const {return m_Time;};
public:
void Add(const double& Energy) {m_Energy+=Energy;};
inline void SetEnergy(const double& Energy) {m_Energy = Energy;};
inline void SetLight(const double& Light) {m_Light = Light;};
};
// Manage a vector of DSSD hit
class PlasticBarDataVector{
public:
......@@ -70,56 +85,93 @@ namespace PlasticBar {
~PlasticBarDataVector(){};
private:
vector<PlasticBarData> m_Data;
std::vector<PlasticBarData> m_Data;
public:
vector<PlasticBarData>::iterator find(const unsigned int& index) ;
std::vector<PlasticBarData>::iterator find(const unsigned int& index);
inline void clear(){m_Data.clear();} ;
inline vector<PlasticBarData>::iterator end() {return m_Data.end();};
inline vector<PlasticBarData>::iterator begin() {return m_Data.begin();};
inline std::vector<PlasticBarData>::iterator end() {return m_Data.end();};
inline std::vector<PlasticBarData>::iterator begin() {return m_Data.begin();};
inline unsigned int size() {return m_Data.size();};
inline void Add(const unsigned int& index,const double& Energy) {find(index)->Add(Energy);};
inline void Set(const double& Energy, const double& Time, const vector<unsigned int>& Nesting) {m_Data.push_back(PlasticBarData(Energy,Time,Nesting));};
PlasticBarData* operator[](const unsigned int& i){return &m_Data[i];};
inline void Set(const double& Energy, const double& Light, const double& Position, const double& Time, const std::vector<unsigned int>& Nesting) {
m_Data.push_back(PlasticBarData(Energy, Light, Position, Time, Nesting));
};
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class PS_PlasticBar : public G4VPrimitiveScorer{
public: // with description
PS_PlasticBar(G4String name, vector<G4int> NestingLevel,G4int depth=0);
~PS_PlasticBar();
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class PS_PlasticBar : public G4VPrimitiveScorer{
public: // with description
PS_PlasticBar(G4String name, std::vector<G4int>NestingLevel, G4int depth=0);
~PS_PlasticBar();
protected: // with description
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
double EnergyToLight(double Energy, int Z);
void Compute_Light();
public:
void Initialize(G4HCofThisEvent*);
void EndOfEvent(G4HCofThisEvent*);
void clear();
void DrawAll();
void PrintAll();
private: // How much level of volume nesting should be considered
// Give the list of the nesting level at which the copy number should be return.
// 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
vector<G4int> m_NestingLevel;
void Initialize(G4HCofThisEvent*);
void EndOfEvent(G4HCofThisEvent*);
void clear();
void DrawAll();
void PrintAll();
private:
// How much level of volume nesting should be considered
// Give the list of the nesting level at which the copy number should be return.
// 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
std::vector<G4int> m_NestingLevel;
private:
PlasticBarDataVector m_Data;
double t_Energy;
double t_Time;
vector<unsigned int> t_Level;
double t_ParticleMultiplicity;
PlasticBarDataVector m_Data;
private:
std::vector< std::array<double,7> > t_Energy_by_ChargeNumber;
std::vector< std::array<double,7> > t_Light_by_ChargeNumber;
inline void AddEntry() {
t_Energy_by_ChargeNumber.push_back({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
t_Light_by_ChargeNumber.push_back({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
};
std::array<double,7> t_TotalEnergy_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
std::array<double,7> t_TotalLight_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
private:
double t_Energy;
double t_Light;
double t_Position;
double t_Time;
std::vector<unsigned int> t_Level;
private:
bool t_DoConversion = 1;
const double t_Convertor[7][4] = {
{1.0 , 0.0 , 0.0 , 0.0 }, //electron
{0.90 , 7.55 , 0.099 , 0.74}, //proton
{0.78 , 39.3 , 0.022 , 0.91}, //alpha
{0.61 , 57.1 , 0.012 , 0.95}, //7Li
{0.44 , 45.5 , 0.010 , 0.92}, //9Be
{0.35 , 34.5 , 0.011 , 0.91}, //11B
{0.30 , 25.6 , 0.013 , 0.91}, //12C
};
public:
inline unsigned int GetMult() {return m_Data.size();};
inline double GetEnergy(const unsigned int& i) {return m_Data[i]->GetEnergy();};
inline double GetLight(const unsigned int& i) {return m_Data[i]->GetLight();};
inline double GetPosition(const unsigned int& i) {return m_Data[i]->GetPosition();};
inline double GetTime(const unsigned int& i) {return m_Data[i]->GetTime();};
inline vector<unsigned int> GetLevel(const unsigned int& i) {return m_Data[i]->GetLevel();};
};
inline std::vector<unsigned int> GetLevel(const unsigned int& i) {
return m_Data[i]->GetLevel();
};
};
}
#endif
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