From a0f17b90846cb45189ccb29b35151f75350c0152 Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Tue, 6 May 2014 11:11:12 +0100
Subject: [PATCH] * Adding missing SIliconScorer File * Fixing minor bug in
NPBeam and NPReaction
---
NPLib/DummyDetector/TDUMMYDetectorData.cxx | 21 +-
NPLib/Physics/NPBeam.cxx | 5 +-
NPLib/Physics/NPBeam.h | 14 +-
NPLib/Physics/NPReaction.cxx | 14 +-
NPLib/SSSD/TSSSDData.cxx | 29 +-
NPLib/VDetector/DetectorManager.cxx | 9 +-
.../include/EventGeneratorTwoBodyReaction.hh | 2 +-
NPSimulation/src/EventGeneratorBeam.cc | 3 +-
.../src/EventGeneratorTwoBodyReaction.cc | 6 +-
NPSimulation/src/ParticleStack.cc | 2 -
NPSimulation/src/SiliconScorers.cc | 244 ++++++++++++
NPSimulation/src/Target.cc | 347 +++++++++---------
12 files changed, 471 insertions(+), 225 deletions(-)
create mode 100644 NPSimulation/src/SiliconScorers.cc
diff --git a/NPLib/DummyDetector/TDUMMYDetectorData.cxx b/NPLib/DummyDetector/TDUMMYDetectorData.cxx
index 1d60e763a..5a5b47e67 100644
--- a/NPLib/DummyDetector/TDUMMYDetectorData.cxx
+++ b/NPLib/DummyDetector/TDUMMYDetectorData.cxx
@@ -24,30 +24,27 @@
#include "TDUMMYDetectorData.h"
ClassImp(TDUMMYDetectorData)
-TDUMMYDetectorData::TDUMMYDetectorData()
-{
+////////////////////////////////////////////////////////////////////////////////
+TDUMMYDetectorData::TDUMMYDetectorData(){
}
-
-TDUMMYDetectorData::~TDUMMYDetectorData()
-{
+////////////////////////////////////////////////////////////////////////////////
+TDUMMYDetectorData::~TDUMMYDetectorData(){
}
-
-void TDUMMYDetectorData::Clear()
-{
+////////////////////////////////////////////////////////////////////////////////
+void TDUMMYDetectorData::Clear(){
fDUMMYDetector_Energy.clear();
fDUMMYDetector_Number.clear();
fDUMMYDetector_Time.clear();
}
-
-void TDUMMYDetectorData::Dump() const
-{
- cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+////////////////////////////////////////////////////////////////////////////////
+void TDUMMYDetectorData::Dump() const{
+ cout << "XXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
for(unsigned short i = 0 ; i<fDUMMYDetector_Energy.size() ; i ++)
{
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index 5cc783910..c3bfebf35 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -309,8 +309,7 @@ void Beam::ReadConfigurationFile(string Path){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY )
-{
+void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY ){
X = Y = 1000000*cm;
if(fSigmaEnergy!=-1)
@@ -319,7 +318,6 @@ void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, doubl
E = fEnergyHist->GetRandom();
if(fSigmaX!=-1){
-
// Shoot within the target unless target size is null (no limit)
while(sqrt(X*X+Y*Y>fEffectiveTargetSize) || fEffectiveTargetSize == 0){
NPL::RandomGaussian2D(fMeanX, fMeanThetaX, fSigmaX, fSigmaThetaX, X, ThetaX);
@@ -333,7 +331,6 @@ void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, doubl
fYPhiYHist->GetRandom2(Y,PhiY);
}
}
-
Z = fTargetZ + fEffectiveTargetThickness*(gRandom->Uniform()-0.5);
}
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index 71bdb92c0..197e596d9 100755
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -115,14 +115,22 @@ namespace NPL{
double fEffectiveTargetThickness; // target thickness has seen by the beam
double fTargetAngle;
double fTargetZ;
-
- public: // Event Generation
- void GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY );
+
+ public:
void SetTargetSize(double TargetSize);
void SetTargetThickness(double TargetThickness);
void SetTargetAngle(double TargetAngle);
void SetTargetZ(double TargetZ) {fTargetZ = TargetZ;}
+ double GetTargetSize(){return fTargetSize;}
+ double GetTargetThickness(){return fTargetThickness;}
+ double GetTargetAngle(){return fTargetAngle;}
+ double GetTargetZ() {return fTargetZ;}
+ double GetTargetEffectiveThickness(){return fEffectiveTargetThickness;}
+ double GetTargetEffectiveTargetSize(){return fEffectiveTargetSize;}
+
+ public: // Event Generation
+ void GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY );
public: // Print private paremeter
void Print() const;
};
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index b67f2ccc6..8f8a6e303 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -214,12 +214,12 @@ double Reaction::ShootRandomThetaCM(){
// Those test are there for the tail event of the energy distribution
// In case the energy is outside the range of the 2D histo we take the
- // closest availabile CS
+ // closest available CS
if(Y->FindBin(fBeamEnergy) > Y->GetLast())
binY = Y->GetLast()-1;
else if(Y->FindBin(fBeamEnergy) < Y->GetFirst())
- binY = Y->GetFirst();
+ binY = Y->GetFirst()+1;
else
binY = Y->FindBin(fBeamEnergy);
@@ -248,11 +248,11 @@ void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
// case of inverse kinematics
double theta = fThetaCM;
if (m1 > m2) theta = M_PI - fThetaCM;
-
+
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3*sin(theta),0,pCM_3*cos(theta),ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
- fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
+ fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
fEnergyImpulsionLab_3.Boost(0,0,BetaCM);
fEnergyImpulsionLab_4 = fEnergyImpulsionCM_4;
fEnergyImpulsionLab_4.Boost(0,0,BetaCM);
@@ -517,13 +517,15 @@ void Reaction::ReadConfigurationFile(string Path){
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::initializePrecomputeVariable(){
+ if(fBeamEnergy < 0)
+ fBeamEnergy = 0 ;
+
m1 = fNuclei1->Mass();
m2 = fNuclei2->Mass();
m3 = fNuclei3->Mass() + fExcitation3;
m4 = fNuclei4->Mass() + fExcitation4;
s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1);
-
fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
@@ -537,7 +539,7 @@ void Reaction::initializePrecomputeVariable(){
pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
fImpulsionLab_1 = TVector3(0,0,sqrt(fBeamEnergy*fBeamEnergy + 2*fBeamEnergy*m1));
- fImpulsionLab_2 = TVector3(0,0,0);
+ fImpulsionLab_3 = TVector3(0,0,0);
fEnergyImpulsionLab_1 = TLorentzVector(fImpulsionLab_1,m1+fBeamEnergy);
fEnergyImpulsionLab_2 = TLorentzVector(fImpulsionLab_2,m2);
diff --git a/NPLib/SSSD/TSSSDData.cxx b/NPLib/SSSD/TSSSDData.cxx
index a89703fa5..3e881a81b 100644
--- a/NPLib/SSSD/TSSSDData.cxx
+++ b/NPLib/SSSD/TSSSDData.cxx
@@ -26,8 +26,8 @@ using namespace std;
ClassImp(TSSSDData)
-TSSSDData::TSSSDData()
-{
+////////////////////////////////////////////////////////////////////////////////
+TSSSDData::TSSSDData(){
// Default constructor
// SSSD
@@ -42,11 +42,12 @@ TSSSDData::TSSSDData()
}
-TSSSDData::~TSSSDData()
-{}
+////////////////////////////////////////////////////////////////////////////////
+TSSSDData::~TSSSDData(){
+}
-void TSSSDData::Clear()
-{
+////////////////////////////////////////////////////////////////////////////////
+void TSSSDData::Clear(){
// DSSD
// Energy
fSSSD_StripE_DetectorNbr.clear();
@@ -59,19 +60,23 @@ void TSSSDData::Clear()
}
-
-void TSSSDData::Dump() const
-{
- cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+////////////////////////////////////////////////////////////////////////////////
+void TSSSDData::Dump() const{
+ cout << "XXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
// SSSD
// Energy
cout << "SSSD_StripE_Mult = " << fSSSD_StripE_DetectorNbr.size() << endl;
for (UShort_t i = 0; i < fSSSD_StripE_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fSSSD_StripE_DetectorNbr[i] << " Strip: " << fSSSD_StripE_StripNbr[i] << " Energy: " << fSSSD_StripE_Energy[i] << endl;
+ cout << "DetNbr: " << fSSSD_StripE_DetectorNbr[i]
+ << " Strip: " << fSSSD_StripE_StripNbr[i]
+ << " Energy: " << fSSSD_StripE_Energy[i] << endl;
+
// Time
cout << "SSSD_StripXT_Mult = " << fSSSD_StripT_DetectorNbr.size() << endl;
for (UShort_t i = 0; i < fSSSD_StripT_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fSSSD_StripT_DetectorNbr[i] << " Strip: " << fSSSD_StripT_StripNbr[i] << " Time: " << fSSSD_StripT_Time[i] << endl;
+ cout << "DetNbr: " << fSSSD_StripT_DetectorNbr[i]
+ << " Strip: " << fSSSD_StripT_StripNbr[i]
+ << " Time: " << fSSSD_StripT_Time[i] << endl;
}
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index f39ca9215..b9b5ed31b 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -61,7 +61,7 @@
/////////////////////////////////////////////////////////////////////////////////////////////////
// Default Constructor
-DetectorManager::DetectorManager() {
+DetectorManager::DetectorManager(){
}
/////////////////////////////////////////////////////////////////////////////////////////////////
@@ -81,7 +81,7 @@ void DetectorManager::ReadConfigurationFile(string Path) {
string LineBuffer;
string DataBuffer;
- /////////Boolean////////////////////
+ //////// List of Detector ///////
bool AnnularS1 = false;
bool CATS = false;
bool Charissa = false;
@@ -89,7 +89,6 @@ void DetectorManager::ReadConfigurationFile(string Path) {
bool EXL = false;
bool Exogam = false;
bool GPDTracker = false;
- bool GeneralTarget = false;
bool HYD2Tracker = false;
bool IonisationChamber = false;
bool LaBr3 = false;
@@ -108,6 +107,10 @@ void DetectorManager::ReadConfigurationFile(string Path) {
bool TiaraHyball = false;
bool Trifoil = false;
bool W1 = false;
+ ////////////////////////////////
+ bool GeneralTarget = false;
+
+
//////////////////////////////////////////////////////////////////////////////////////////
ifstream ConfigFile;
ConfigFile.open(Path.c_str());
diff --git a/NPSimulation/include/EventGeneratorTwoBodyReaction.hh b/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
index a43006666..c5ad6db91 100644
--- a/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
+++ b/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
@@ -40,7 +40,7 @@
#include "NPReaction.h"
using namespace std;
-
using namespace CLHEP;
+using namespace CLHEP;
using namespace NPL ;
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index 0354bbeaf..02d1e4d5d 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -72,6 +72,7 @@ void EventGeneratorBeam::ReadConfiguration(string Path,int){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
+
//--------------write the DeDx Table -------------------
if( anEvent->GetEventID()==0){
// Define the particle to be shoot
@@ -98,7 +99,7 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
G4ThreeVector BeamPos(x0,y0,z0);
Beam_theta = BeamDir.theta() ;
Beam_phi = BeamDir.phi() ; Beam_phi *= 1;
- FinalBeamEnergy = m_Target->SlowDownBeam(m_particle, InitialBeamEnergy,z0,Beam_theta);
+ FinalBeamEnergy = m_Target->SlowDownBeam(m_particle, InitialBeamEnergy,z0-m_Beam->GetTargetZ(),Beam_theta);
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
diff --git a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
index 7cc5671b0..ed40af9e0 100644
--- a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
+++ b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
@@ -111,7 +111,6 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
}
}
-
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
@@ -135,9 +134,9 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
// Get the beam particle form the Particle Stack
Particle BeamParticle = m_ParticleStack->SearchAndRemoveParticle(m_BeamName);
m_Reaction->SetBeamEnergy(BeamParticle.GetParticleKineticEnergy());
+
G4double Beam_theta = BeamParticle.GetParticleMomentumDirection().theta();
G4double Beam_phi = BeamParticle.GetParticleMomentumDirection().phi();
-
//////////////////////////////////////////////////////////
///// Build rotation matrix to go from the incident //////
@@ -163,7 +162,7 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
// Shoot and Set a Random ThetaCM
G4double ThetaCM = m_Reaction->ShootRandomThetaCM();
G4double phi = RandFlat::shoot() * 2. * pi;
-
+
//////////////////////////////////////////////////
///// Momentum and angles from kinematics /////
///// Angles are in the beam frame /////
@@ -173,6 +172,7 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
// Compute Kinematic using previously defined ThetaCM
m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
+
// Momentum in beam frame for light particle
G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi),
sin(ThetaLight) * sin(phi),
diff --git a/NPSimulation/src/ParticleStack.cc b/NPSimulation/src/ParticleStack.cc
index 5f036fa36..39d60d06e 100644
--- a/NPSimulation/src/ParticleStack.cc
+++ b/NPSimulation/src/ParticleStack.cc
@@ -19,7 +19,6 @@
* Comment: *
* *
*****************************************************************************/
-
// NPS
#include"ParticleStack.hh"
@@ -236,7 +235,6 @@ void ParticleStack::ShootAllParticle(G4Event* anEvent){
m_particleGun->SetParticleMomentumDirection(m_ParticleStack[i].GetParticleMomentumDirection());
m_particleGun->SetParticlePosition(m_ParticleStack[i].GetParticlePosition());
m_particleGun->GeneratePrimaryVertex(anEvent);
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
m_InitialConditions-> SetParticleName ( m_ParticleStack[i].GetParticleDefinition()->GetParticleName()) ;
m_InitialConditions-> SetThetaCM ( m_ParticleStack[i].GetParticleThetaCM()/deg);
diff --git a/NPSimulation/src/SiliconScorers.cc b/NPSimulation/src/SiliconScorers.cc
new file mode 100644
index 000000000..a19a63bec
--- /dev/null
+++ b/NPSimulation/src/SiliconScorers.cc
@@ -0,0 +1,244 @@
+ /*****************************************************************************
+ * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : February 2013 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * File old the scorer specific to the Sharc Detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * This new type of scorer is aim to become the standard for DSSD,SSSD and *
+ * PAD detector (any Silicon Detector) *
+ *****************************************************************************/
+#include "SiliconScorers.hh"
+#include "G4UnitsTable.hh"
+using namespace SILICONSCORERS ;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Rectangle::PS_Silicon_Rectangle(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_StripPlaneLength = StripPlaneLength;
+ m_StripPlaneWidth = StripPlaneWidth;
+ m_NumberOfStripLength = NumberOfStripLength;
+ m_NumberOfStripWidth = NumberOfStripWidth;
+ m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength;
+ m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
+
+ m_Position = G4ThreeVector(-1000,-1000,-1000);
+ m_DetectorNumber = -1;
+ m_StripLengthNumber = -1;
+ m_StripWidthNumber = -1;
+ m_Index = -1 ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Rectangle::~PS_Silicon_Rectangle(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[10];
+ Infos[0] = aStep->GetTotalEnergyDeposit();
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+
+ m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position);
+
+ m_StripLengthNumber = (int)((m_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength ) + 1 ;
+ m_StripWidthNumber = (int)((m_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ;
+
+ m_StripLengthNumber = m_NumberOfStripLength - m_StripLengthNumber + 1 ;
+ m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ;
+
+ //Rare case where particle is close to edge of silicon plan
+ if (m_StripLengthNumber > m_NumberOfStripLength) m_StripLengthNumber = m_StripLengthNumber;
+ if (m_StripWidthNumber > m_NumberOfStripWidth) m_StripWidthNumber = m_StripWidthNumber;
+
+ Infos[7] = m_DetectorNumber;
+ Infos[8] = m_StripLengthNumber;
+ Infos[9] = m_StripWidthNumber;
+
+
+ m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripLengthNumber * 1e6 + m_StripWidthNumber * 1e9;
+
+ // Check if the particle has interact before, if yes, add up the energies.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Rectangle::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new G4THitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Rectangle::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Rectangle::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Rectangle::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Rectangle::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Annular::PS_Silicon_Annular(G4String name, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double PhiStart,G4double PhiStop, G4int NumberOfStripRing,G4int NumberOfStripSector,G4int NumberOfQuadrant,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+
+ m_StripPlaneInnerRadius = StripPlaneInnerRadius;
+ m_StripPlaneOuterRadius = StripPlaneOuterRadius;
+ m_PhiStart = PhiStart;
+ m_PhiStop = PhiStop;
+ m_NumberOfStripRing = NumberOfStripRing;
+ m_NumberOfStripSector = NumberOfStripSector;
+ m_NumberOfStripQuadrant = NumberOfQuadrant;
+
+ m_StripPitchRing = (m_StripPlaneOuterRadius-m_StripPlaneInnerRadius)/m_NumberOfStripRing;
+ m_StripPitchSector = (m_PhiStop-m_PhiStart)/m_NumberOfStripSector;
+ m_StripPitchQuadrant = (m_PhiStop-m_PhiStart)/m_NumberOfStripQuadrant;
+
+ m_Position = G4ThreeVector(-1000,-1000,-1000);
+ m_uz = G4ThreeVector(0,0,1);
+ m_DetectorNumber = -1;
+ m_StripRingNumber = -1;
+ m_StripSectorNumber = -1;
+ m_Index = -1 ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Annular::~PS_Silicon_Annular(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[11];
+ Infos[0] = aStep->GetTotalEnergyDeposit();
+
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+
+ m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position);
+ m_StripRingNumber = (int) ((m_Position.rho() - m_StripPlaneInnerRadius) / m_StripPitchRing ) + 1 ;
+
+ double phi = m_Position.phi()+M_PI;
+ m_StripSectorNumber = (int) ((phi - m_PhiStart) / m_StripPitchSector ) + 1 ;
+ m_StripQuadrantNumber = (int) ((phi - m_PhiStart) /m_StripPitchQuadrant) + 1 ;
+
+ //Rare case where particle is close to edge of silicon plan
+ if (m_StripRingNumber > m_NumberOfStripRing) m_StripRingNumber = m_NumberOfStripRing;
+ if (m_StripSectorNumber > m_NumberOfStripSector) m_StripSectorNumber = m_NumberOfStripSector;
+ if (m_StripQuadrantNumber > m_NumberOfStripQuadrant) m_StripQuadrantNumber = m_NumberOfStripQuadrant;
+
+ Infos[7] = m_DetectorNumber;
+ Infos[8] = m_StripRingNumber;
+ Infos[9] = m_StripSectorNumber;
+ Infos[10] = m_StripQuadrantNumber;
+
+
+ m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripRingNumber * 1e6 + m_NumberOfStripSector * 1e9;
+
+ // Check if the particle has interact before, if yes, add up the energies.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Annular::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new G4THitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Annular::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Annular::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Annular::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Annular::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index 8b80e8194..6883e3ba5 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -57,8 +57,7 @@ using namespace std;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Specific Method of this class
-Target::Target()
-{
+Target::Target(){
m_TargetType = true;
m_TargetThickness = 0 ;
m_TargetAngle = 0 ;
@@ -70,130 +69,126 @@ Target::Target()
m_TargetNbLayers = 50; // Number of steps by default
}
-G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Temperature, G4double Pressure)
-{
-
+G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Temperature, G4double Pressure){
G4Material* myMaterial;
-
+
if (MaterialName == "D2") {
G4double density = 0;
-
+
if(Pressure == 0 )
- if (Temperature == 0) {
+ if (Temperature == 0) {
density = 0.000083771* g / cm3;
G4cout << "CryoTarget temp set to 300K with P = 1bar" << G4endl;
}
-
- else if (Pressure == 1) {
- G4cout << "CryoTarget pressure set to 1 bar" << G4endl;
-
- if (Temperature == 24) {
- density = 0.0020182 * g / cm3;
- G4cout << "CryoTarget temp set to 24K" << G4endl;
- }
-
- else if (Temperature == 25) {
- density = 0.0019377 * g / cm3;
- G4cout << "CryoTarget temp set to 25K" << G4endl;
- }
-
- else if (Temperature == 26) {
- density = 0.001863 * g / cm3;
- G4cout << "CryoTarget temp set to 26K" << G4endl;
- }
-
- else if (Temperature == 30) {
- density = 0.00020475 * g / cm3;
- G4cout << "CryoTarget temp set to 30K" << G4endl;
- }
-
- else if (Temperature == 300) {
- density = 8.3771e-5* g / cm3;
- G4cout << "CryoTarget temp set to 30K" << G4endl;
- }
-
-
- else {
- G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
- }
- }
-
- else if (Pressure == 0.5) {
- G4cout << "CryoTarget pressure set to 0.5 bar" << G4endl;
-
- if (Temperature == 24) {
- density = 0.0010091 * g / cm3;
- G4cout << "CryoTarget temp set to 24K" << G4endl;
- }
-
- else if (Temperature == 25) {
- density = 0.00096875 * g / cm3;
- G4cout << "CryoTarget temp set to 25K" << G4endl;
- }
-
- else if (Temperature == 26) {
- density = 0.00093149 * g / cm3;
- G4cout << "CryoTarget temp set to 26K" << G4endl;
- }
-
-
- else if (Pressure == 0.7) {
- G4cout << "CryoTarget pressure set to 0.7 bar" << G4endl;
-
- if (Temperature == 26) {
- density = 0.0013125 * g / cm3;
+
+ else if (Pressure == 1) {
+ G4cout << "CryoTarget pressure set to 1 bar" << G4endl;
+
+ if (Temperature == 24) {
+ density = 0.0020182 * g / cm3;
+ G4cout << "CryoTarget temp set to 24K" << G4endl;
+ }
+
+ else if (Temperature == 25) {
+ density = 0.0019377 * g / cm3;
+ G4cout << "CryoTarget temp set to 25K" << G4endl;
+ }
+
+ else if (Temperature == 26) {
+ density = 0.001863 * g / cm3;
G4cout << "CryoTarget temp set to 26K" << G4endl;
}
+
+ else if (Temperature == 30) {
+ density = 0.00020475 * g / cm3;
+ G4cout << "CryoTarget temp set to 30K" << G4endl;
+ }
+
+ else if (Temperature == 300) {
+ density = 8.3771e-5* g / cm3;
+ G4cout << "CryoTarget temp set to 30K" << G4endl;
+ }
+
+
+ else {
+ G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
+ }
}
-
-
- else {
- G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
+
+ else if (Pressure == 0.5) {
+ G4cout << "CryoTarget pressure set to 0.5 bar" << G4endl;
+
+ if (Temperature == 24) {
+ density = 0.0010091 * g / cm3;
+ G4cout << "CryoTarget temp set to 24K" << G4endl;
+ }
+
+ else if (Temperature == 25) {
+ density = 0.00096875 * g / cm3;
+ G4cout << "CryoTarget temp set to 25K" << G4endl;
+ }
+
+ else if (Temperature == 26) {
+ density = 0.00093149 * g / cm3;
+ G4cout << "CryoTarget temp set to 26K" << G4endl;
+ }
+
+
+ else if (Pressure == 0.7) {
+ G4cout << "CryoTarget pressure set to 0.7 bar" << G4endl;
+
+ if (Temperature == 26) {
+ density = 0.0013125 * g / cm3;
+ G4cout << "CryoTarget temp set to 26K" << G4endl;
+ }
+ }
+
+
+ else {
+ G4cout << ">>> !!!!WARNING INVALID TEMP FOR CRYOGENIC TARGET!!!! <<<" << G4endl;
+ }
}
- }
-
-
-
+
G4Element* D = new G4Element("Deuteron" , "D" , 1., 2.0141*g / mole);
myMaterial = new G4Material("D2", density, 1, kStateGas, Temperature, Pressure);
myMaterial->AddElement(D, 2);
return(myMaterial);
}
-
+
else if (MaterialName == "D2_solid") {
G4Element* D = new G4Element("Deuteron", "D", 1., 2.0141*g / mole);
G4Material* myMaterial = new G4Material("D2_solid", 0.0715*g/cm3, 1);
myMaterial->AddElement(D, 2);
return(myMaterial);
}
-
+
else if (MaterialName == "H2_solid") {
G4Element* H = new G4Element("Hydrogen", "H", 1., 1.01*g / mole);
G4Material* myMaterial = new G4Material("H2_solid", 0.0715*g/cm3, 1);
myMaterial->AddElement(H, 2);
return(myMaterial);
}
-
+
else if (MaterialName == "Mylar") {
G4cout << "Mylar Material" << G4endl ;
G4Element* H = new G4Element("Hydrogen", "H" , 1. , 1.01 *g / mole);
G4Element* C = new G4Element("Carbon" , "C" , 6. , 12.011*g / mole);
G4Element* O = new G4Element("Oxygen" , "O" , 8. , 16.00 *g / mole);
-
+
G4Material* myMaterial = new G4Material("Mylar", 1.397*g / cm3, 3);
myMaterial->AddElement(C , 10);
myMaterial->AddElement(H , 8);
myMaterial->AddElement(O , 4);
return myMaterial;
}
-
+
else if (MaterialName == "Harvar") {
G4Element* Co = new G4Element("Cobalt" , "Co" , 27 , 58.933*g / mole);
G4Element* Cr = new G4Element("Cromium", "Cr" , 24 , 51.996*g / mole);
G4Element* Ni = new G4Element("Nickel" , "Ni" , 28 , 58.69*g / mole);
G4Element* Fe = new G4Element("Iron" , "Fe" , 26 , 55.847*g / mole);
G4Element* W = new G4Element("Teflon" , "W" , 74 , 183.5*g / mole);
-
+
G4Material* myMaterial = new G4Material("Havar", 8.3*g / cm3, 5);
myMaterial->AddElement(Co , 42);
myMaterial->AddElement(Cr , 20);
@@ -202,45 +197,45 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
myMaterial->AddElement(W , 1);
return myMaterial;
}
-
+
else if (MaterialName == "CD2") {
G4Element* C = new G4Element("Carbon" , "C" , 6. , 12.011*g / mole);
G4Element* D = new G4Element("Deuteron" , "D" , 1., 2.0141*g / mole);
-
+
G4Material* myMaterial = new G4Material("CD2", 1.06*g / cm3, 2);
myMaterial->AddElement(C , 1);
myMaterial->AddElement(D , 2);
return myMaterial;
}
-
+
else if (MaterialName == "CH2") {
G4Element* C = new G4Element("Carbon" , "C" , 6. , 12.011*g / mole);
G4Element* H = new G4Element("Hydrogen", "H" , 1. , 1.01 *g / mole);
-
+
G4Material* myMaterial = new G4Material("CH2", 0.93*g / cm3, 2);
myMaterial->AddElement(C , 1);
myMaterial->AddElement(H , 2);
return myMaterial;
}
-
+
else if (MaterialName == "208Pb") {
G4Element* Pb = new G4Element("Lead" , "Pb" , 82. , 207.2*g / mole);
-
+
G4Material* myMaterial = new G4Material("208Pb", 11.342*g / cm3, 1);
myMaterial->AddElement(Pb , 1);
return myMaterial;
}
-
+
else if (MaterialName == "Si") {
G4Material* myMaterial = new G4Material("Si", 14., 28.0855 * g / mole, 2.321 * g / cm3);
return myMaterial;
}
-
+
else if (MaterialName == "Al") {
G4Material* myMaterial = new G4Material("Aluminium", 13., 26.98 * g / mole, 2.702 * g / cm3);
return myMaterial;
}
-
+
else {
G4cout << "No Matching Material in the Target Library Default is Vacuum" << G4endl;
G4Element* N = new G4Element("Nitrogen", "N", 7., 14.01*g / mole);
@@ -260,30 +255,29 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void Target::ReadConfiguration(string Path)
-{
+void Target::ReadConfiguration(string Path){
ifstream ConfigFile;
ConfigFile.open(Path.c_str());
string LineBuffer;
string DataBuffer;
-
+
bool ReadingStatusTarget = false ;
bool ReadingStatusCryoTarget = false ;
-
+
bool check_Thickness = false ;
bool check_Radius = false ;
-// bool check_Angle = false ;
+ // bool check_Angle = false ;
bool check_Material = false ;
bool check_X = false ;
bool check_Y = false ;
bool check_Z = false ;
-// bool check_m_TargetNbLayers = false;
-
+ // bool check_m_TargetNbLayers = false;
+
bool check_Temperature = false ;
bool check_Pressure = false ;
bool check_WinThickness = false ;
bool check_WinMaterial = false ;
-
+
int VerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
while (!ConfigFile.eof()) {
@@ -298,251 +292,248 @@ void Target::ReadConfiguration(string Path)
m_TargetType = false ;
ReadingStatusCryoTarget = true ;
}
-
+
while (ReadingStatusTarget) {
ConfigFile >> DataBuffer;
-
+
//Search for comment Symbol %
if (DataBuffer.compare(0, 1, "%") == 0) { ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
-
+
else if (DataBuffer.compare(0, 10, "THICKNESS=") == 0) {
check_Thickness = true ;
ConfigFile >> DataBuffer;
m_TargetThickness = atof(DataBuffer.c_str()) * micrometer;
if(VerboseLevel==1) cout << "Target Thickness: " << m_TargetThickness / micrometer << " micrometer" << endl;
}
-
+
else if (DataBuffer.compare(0, 6, "ANGLE=") == 0) {
-// check_Angle = true ;
+ // check_Angle = true ;
ConfigFile >> DataBuffer;
m_TargetAngle = atof(DataBuffer.c_str()) * deg;
if(VerboseLevel==1) cout << "Target Angle: " << m_TargetAngle / deg << endl ;
}
-
+
else if (DataBuffer.compare(0, 7, "RADIUS=") == 0) {
check_Radius = true ;
ConfigFile >> DataBuffer;
m_TargetRadius = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << "Target Radius: " << m_TargetRadius / mm << " mm " << endl;
}
-
+
else if (DataBuffer.compare(0, 9, "MATERIAL=") == 0) {
check_Material = true ;
ConfigFile >> DataBuffer;
m_TargetMaterial = GetMaterialFromLibrary(DataBuffer);
if(VerboseLevel==1) cout << "Target Material: " << m_TargetMaterial << endl ;
}
-
+
else if (DataBuffer.compare(0, 2, "X=") == 0) {
check_X = true ;
ConfigFile >> DataBuffer;
m_TargetX = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << "Target coordinate (mm): ( " << m_TargetX / mm << " ; ";
}
-
+
else if (DataBuffer.compare(0, 2, "Y=") == 0) {
check_Y = true ;
ConfigFile >> DataBuffer;
m_TargetY = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << m_TargetY / mm << " ; ";
}
-
+
else if (DataBuffer.compare(0, 2, "Z=") == 0) {
check_Z = true ;
ConfigFile >> DataBuffer;
m_TargetZ = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << m_TargetZ / mm << " )" << endl ;
}
-
+
else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) {
-// check_m_TargetNbLayers = true ;
+ // check_m_TargetNbLayers = true ;
ConfigFile >> DataBuffer;
m_TargetNbLayers = atoi(DataBuffer.c_str());
if(VerboseLevel==1) cout << "Number of steps for slowing down the beam in target: " << m_TargetNbLayers << endl;
}
-
+
///////////////////////////////////////////////////
// If no Beam Token and no comment, toggle out
else{ReadingStatusTarget = false;
G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;
}
-
+
///////////////////////////////////////////////////
// If all Token found toggle out
if( check_Thickness && check_Radius && check_Material && check_X && check_Y && check_Z ){
ReadingStatusTarget = false ;
}
}
-
+
while(ReadingStatusCryoTarget){
-
+
ConfigFile >> DataBuffer;
-
+
//Search for comment Symbol %
if (DataBuffer.compare(0, 1, "%") == 0) {/*Do Nothing*/;}
-
+
else if (DataBuffer.compare(0, 10, "THICKNESS=") == 0) {
check_Thickness = true ;
ConfigFile >> DataBuffer;
m_TargetThickness = atof(DataBuffer.c_str()) * micrometer;
if(VerboseLevel==1) cout << "Target Thickness: " << m_TargetThickness / micrometer << "um" << endl ;
}
-
+
else if (DataBuffer.compare(0, 7, "RADIUS=") == 0) {
check_Radius = true ;
ConfigFile >> DataBuffer;
m_TargetRadius = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << "Target Radius: " << m_TargetRadius / mm << "mm" << endl ;
}
-
+
else if (DataBuffer.compare(0, 12, "TEMPERATURE=") == 0) {
check_Temperature = true ;
ConfigFile >> DataBuffer;
if(VerboseLevel==1) m_TargetTemperature = atof(DataBuffer.c_str());
}
-
+
else if (DataBuffer.compare(0, 9, "PRESSURE=") == 0) {
check_Pressure = true ;
ConfigFile >> DataBuffer;
m_TargetPressure = atof(DataBuffer.c_str());
}
-
+
else if (DataBuffer.compare(0, 9, "MATERIAL=") == 0) {
check_Material = true ;
ConfigFile >> DataBuffer;
m_TargetMaterial = GetMaterialFromLibrary(DataBuffer, m_TargetTemperature, m_TargetPressure);
if(VerboseLevel==1) cout << "Target Material: " << m_TargetMaterial << endl ;
}
-
+
else if (DataBuffer.compare(0, 17, "WINDOWSTHICKNESS=") == 0) {
check_WinThickness = true ;
ConfigFile >> DataBuffer;
m_WindowsThickness = atof(DataBuffer.c_str()) * micrometer;
if(VerboseLevel==1) cout << "Windows Thickness: " << m_WindowsThickness / micrometer << "um" << endl ;
}
-
+
else if (DataBuffer.compare(0, 16, "WINDOWSMATERIAL=") == 0) {
check_WinMaterial = true ;
ConfigFile >> DataBuffer;
m_WindowsMaterial = GetMaterialFromLibrary(DataBuffer);
if(VerboseLevel==1) cout << "Windows Material: " << m_WindowsMaterial << endl ;
}
-
+
else if (DataBuffer.compare(0, 2, "X=") == 0) {
check_X = true ;
ConfigFile >> DataBuffer;
m_TargetX = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << "Target coordinate (mm): ( " << m_TargetX / mm << " ; ";
}
-
+
else if (DataBuffer.compare(0, 2, "Y=") == 0) {
check_Y = true ;
ConfigFile >> DataBuffer;
m_TargetY = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << m_TargetY / mm << " ; ";
}
-
+
else if (DataBuffer.compare(0, 2, "Z=") == 0) {
check_Z = true ;
ConfigFile >> DataBuffer;
m_TargetZ = atof(DataBuffer.c_str()) * mm;
if(VerboseLevel==1) cout << m_TargetZ / mm << " )" << endl ;
}
-
+
else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) {
-// check_m_TargetNbLayers = true ;
+ // check_m_TargetNbLayers = true ;
ConfigFile >> DataBuffer;
m_TargetNbLayers = atoi(DataBuffer.c_str());
if(VerboseLevel==1) cout << "Number of steps for slowing down the beam in target: " << m_TargetNbLayers << endl;
}
-
+
///////////////////////////////////////////////////
// If no Beam Token and no comment, toggle out
else{
ReadingStatusCryoTarget = false;
G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;
}
-
+
///////////////////////////////////////////////////
// If all Token found toggle out
if( check_Thickness && check_Radius && check_Material && check_X && check_Y && check_Z && check_WinThickness && check_WinMaterial && check_Pressure && check_Temperature)
m_EffectiveThickness = m_TargetThickness / cos(m_TargetAngle);
ReadingStatusCryoTarget = false ;
-
+
}
-
+
}
// if the target as a null radius then no target exist
if(m_TargetRadius==0) {m_TargetThickness=0;m_TargetRadius=0.1*um;}
-
+
}
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
-void Target::ConstructDetector(G4LogicalVolume* world)
-{
-
-
+void Target::ConstructDetector(G4LogicalVolume* world){
if (m_TargetType) { // case of standard target
-
+
if (m_TargetThickness > 0) {
G4Tubs* solidTarget = new G4Tubs("solidTarget", 0, m_TargetRadius, 0.5*m_TargetThickness, 0*deg, 360*deg);
G4LogicalVolume* logicTarget = new G4LogicalVolume(solidTarget, m_TargetMaterial, "logicTarget");
-
+
// rotation of target
G4RotationMatrix *rotation = new G4RotationMatrix();
rotation->rotateY(m_TargetAngle);
-
+
new G4PVPlacement(rotation, G4ThreeVector(m_TargetX, m_TargetY, m_TargetZ), logicTarget, "Target", world, false, 0);
-
+
G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));//Blue
logicTarget->SetVisAttributes(TargetVisAtt);
}
}
-
+
else { // case of cryogenic target
-
+
if (m_TargetThickness > 0) {
G4Tubs* solidTarget = new G4Tubs("solidTarget", 0, m_TargetRadius, 0.5*m_TargetThickness, 0*deg, 360*deg);
G4LogicalVolume* logicTarget = new G4LogicalVolume(solidTarget, m_TargetMaterial, "logicTarget");
new G4PVPlacement(0, G4ThreeVector(m_TargetX, m_TargetY, m_TargetZ), logicTarget, "Target", world, false, 0);
-
+
G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));//Blue
logicTarget->SetVisAttributes(TargetVisAtt);
}
-
-
+
+
if (m_WindowsThickness > 0) {
G4ThreeVector TargetPos = G4ThreeVector(m_TargetX, m_TargetY, m_TargetZ);
-
+
G4Tubs* solidWindowsF =
- new G4Tubs("solidTargetWindowsF", 0, m_TargetRadius, 0.5*m_WindowsThickness, 0*deg, 360*deg);
+ new G4Tubs("solidTargetWindowsF", 0, m_TargetRadius, 0.5*m_WindowsThickness, 0*deg, 360*deg);
G4LogicalVolume* logicWindowsF = new G4LogicalVolume(solidWindowsF, m_WindowsMaterial, "logicTargetWindowsF");
-
+
G4Tubs* solidWindowsB =
- new G4Tubs("solidTargetWindowsB", 0, m_TargetRadius, 0.5*m_WindowsThickness, 0*deg, 360*deg);
+ new G4Tubs("solidTargetWindowsB", 0, m_TargetRadius, 0.5*m_WindowsThickness, 0*deg, 360*deg);
G4LogicalVolume* logicWindowsB = new G4LogicalVolume(solidWindowsB, m_WindowsMaterial, "logicTargetWindowsB");
-
+
new G4PVPlacement( 0 ,
- TargetPos + G4ThreeVector(0., 0., 0.5*(m_TargetThickness + m_WindowsThickness)) ,
- logicWindowsF ,
- "Target Window Front" ,
- world ,
- false, 0 );
-
+ TargetPos + G4ThreeVector(0., 0., 0.5*(m_TargetThickness + m_WindowsThickness)) ,
+ logicWindowsF ,
+ "Target Window Front" ,
+ world ,
+ false, 0 );
+
new G4PVPlacement( 0 ,
- TargetPos + G4ThreeVector(0., 0., -0.5*(m_TargetThickness + m_WindowsThickness)) ,
- logicWindowsB ,
- "Target Window Back" ,
- world ,
- false, 0 );
-
+ TargetPos + G4ThreeVector(0., 0., -0.5*(m_TargetThickness + m_WindowsThickness)) ,
+ logicWindowsB ,
+ "Target Window Back" ,
+ world ,
+ false, 0 );
+
G4VisAttributes* WindowsVisAtt = new G4VisAttributes(G4Colour(0.5, 1., 0.5));
logicWindowsF->SetVisAttributes(WindowsVisAtt);
logicWindowsB->SetVisAttributes(WindowsVisAtt);
}
}
-
+
}
// Add Detector branch to the EventTree.
@@ -560,7 +551,7 @@ void Target::ReadSensitive(const G4Event*)
G4double Target::SlowDownBeam(G4ParticleDefinition* Beam, G4double IncidentEnergy, G4double ZInteraction, G4double IncidentTheta){
G4double ThicknessBeforeInteraction = abs(ZInteraction - 0.5*m_EffectiveThickness) / cos(m_TargetAngle);
G4double dedx,de;
-
+
if(m_TargetType){
G4EmCalculator emCalculator;
if(m_TargetThickness!=0){
@@ -571,7 +562,7 @@ G4double Target::SlowDownBeam(G4ParticleDefinition* Beam, G4double IncidentEnerg
}
}
}
-
+
else{
G4EmCalculator emCalculator;
// Windows
@@ -581,7 +572,7 @@ G4double Target::SlowDownBeam(G4ParticleDefinition* Beam, G4double IncidentEnerg
de = dedx * m_TargetNbLayers * m_WindowsThickness / cos(IncidentTheta);
IncidentEnergy -= de;
}
-
+
// Target
if(m_TargetThickness!=0)
for (G4int i = 0; i < m_TargetNbLayers; i++){
@@ -598,11 +589,11 @@ void Target::RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double
if (SigmaX != 0) {
X = 2 * NumberOfSigma*SigmaX;
while (X > NumberOfSigma*SigmaX) X = RandGauss::shoot(MeanX, SigmaX);
-
+
double a = NumberOfSigma * SigmaX/2;
double b = NumberOfSigma * SigmaY/2;
double SigmaYPrim = b * sqrt(1 - X*X / (a*a));
-
+
SigmaYPrim = 2*SigmaYPrim / NumberOfSigma;
Y = RandGauss::shoot(MeanY, SigmaYPrim);
}
@@ -620,35 +611,35 @@ void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4doub
G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_TargetMaterial->GetName() + ".G4table";
ofstream File ;
File.open(Path) ;
-
+
if(!File) return ;
-
+
File << "Table from Geant4 generate using NPSimulation \t"
- << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_TargetMaterial->GetName() << endl ;
-
+ << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_TargetMaterial->GetName() << endl ;
+
G4EmCalculator emCalculator;
-
+
for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/1000.){
G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_TargetMaterial);
File << E/MeV << "\t" << dedx/(MeV/micrometer) << endl ;
}
-
+
File.close();
if(!m_TargetType){
G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + m_WindowsMaterial->GetName() + ".G4table";
File.open(Path) ;
if(!File) return ;
File << "Table from Geant4 generate using NPSimulation \t "
- << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_WindowsMaterial->GetName() << endl ;
-
+ << "Particle: " << Particle->GetParticleName() << "\tMaterial: " << m_WindowsMaterial->GetName() << endl ;
+
for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/10.){
// G4double dedx = emCalculator.ComputeTotalDEDX(E, Particle, m_WindowsMaterial);
G4double dedx = emCalculator.ComputeDEDX( E, Particle ,
- "ionIoni", m_WindowsMaterial);
+ "ionIoni", m_WindowsMaterial);
File << E/MeV << "\t" << dedx/(MeV/micrometer) << endl ;
}
}
File.close();
-
+
}
}
--
GitLab