From 302d0717b5ec1aba8cbabe387561644fd4fd31f0 Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Wed, 29 May 2013 22:42:08 +0000
Subject: [PATCH] * Fix some issue in the scoring of Sharc * - In the case of
multiinteraction in the same volume, the scorer was keeping only the last
energy deposit, instead of summing them - The threshold has been move in the
Read sensitive and is no longer in the scorer, as the energy deposit can be
down in multiple step each of them below threshold and be ignored if the
threshold is applyed individually.
---
.../DetectorConfiguration/AnnularS1.detector | 4 +-
Inputs/DetectorConfiguration/Sharc.detector | 56 +++---
Inputs/EventGenerator/10He.reaction | 30 ++--
.../pMakePhysicalTree/DataProcessor.cxx | 12 +-
NPAnalysis/pMakePhysicalTree/RunToTreat.txt | 2 +-
.../InitialConditions/TInitialConditions.cxx | 10 --
NPLib/InitialConditions/TInitialConditions.h | 10 +-
NPLib/Physics/NPBeam.cxx | 32 ++++
NPLib/Physics/NPBeam.h | 1 +
NPLib/Physics/NPNucleus.cxx | 6 +
NPLib/Physics/NPReaction.cxx | 141 ++++++++++++----
NPLib/Physics/NPReaction.h | 14 +-
NPLib/Sharc/TSharcPhysics.cxx | 11 +-
NPSimulation/include/Sharc.hh | 7 +-
NPSimulation/include/SharcScorers.hh | 19 ++-
NPSimulation/src/DetectorConstruction.cc | 27 ++-
NPSimulation/src/Sharc.cc | 159 +++++++++---------
NPSimulation/src/SharcScorers.cc | 58 ++++---
NPSimulation/src/Target.cc | 2 +-
NPSimulation/vis.mac | 4 +-
20 files changed, 381 insertions(+), 224 deletions(-)
diff --git a/Inputs/DetectorConfiguration/AnnularS1.detector b/Inputs/DetectorConfiguration/AnnularS1.detector
index 145d0667d..1a910320f 100644
--- a/Inputs/DetectorConfiguration/AnnularS1.detector
+++ b/Inputs/DetectorConfiguration/AnnularS1.detector
@@ -11,8 +11,8 @@ Target
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AnnularS1
Z= 100
- VIS= all
+ VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AnnularS1
Z= -100
- VIS= all
+ VIS= all
diff --git a/Inputs/DetectorConfiguration/Sharc.detector b/Inputs/DetectorConfiguration/Sharc.detector
index 2ede44d27..3884cf6fa 100644
--- a/Inputs/DetectorConfiguration/Sharc.detector
+++ b/Inputs/DetectorConfiguration/Sharc.detector
@@ -1,8 +1,10 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GeneralTarget
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%0.5mg/cm2
Target
- THICKNESS= 0.001
+ THICKNESS= 4.8
RADIUS= 5
MATERIAL= CD2
ANGLE= 0
@@ -18,30 +20,30 @@ Sharc
Z= -66
R= 0
Phi= 0
- ThicknessDector= 100
+ ThicknessDector= 500
SharcQQQ
Z= -66
R= 0
Phi= 90
- ThicknessDector= 100
+ ThicknessDector= 500
SharcQQQ
Z= -66
R= 0
Phi= 180
- ThicknessDector= 100
+ ThicknessDector= 500
SharcQQQ
Z= -66
R= 0
Phi= 270
- ThicknessDector= 100
+ ThicknessDector= 500
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream Box
SharcBOX
Z= -30.775
- ThicknessDector1= 100
- ThicknessDector2= 100
- ThicknessDector3= 100
- ThicknessDector4= 100
+ ThicknessDector1= 500
+ ThicknessDector2= 500
+ ThicknessDector3= 500
+ ThicknessDector4= 500
ThicknessPAD1= 0
ThicknessPAD2= 0
ThicknessPAD3= 0
@@ -50,34 +52,34 @@ Sharc
%Down Stream Box
SharcBOX
Z= 34.3
- ThicknessDector1= 100
- ThicknessDector2= 100
- ThicknessDector3= 100
- ThicknessDector4= 100
- ThicknessPAD1= 1000
- ThicknessPAD2= 1000
- ThicknessPAD3= 1000
- ThicknessPAD4= 1000
+ ThicknessDector1= 1000
+ ThicknessDector2= 1000
+ ThicknessDector3= 1000
+ ThicknessDector4= 1000
+ ThicknessPAD1= 0
+ ThicknessPAD2= 0
+ ThicknessPAD3= 0
+ ThicknessPAD4= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %Downstream CD
+%Downstream CD
%SharcQQQ
- Z= 66
+ Z= 70
R= 0
Phi= 0
- ThicknessDector= 100
+ ThicknessDector= 500
%SharcQQQ
- Z= 66
+ Z= 70
R= 0
Phi= 90
- ThicknessDector= 100
- % SharcQQQ
- Z= 66
+ ThicknessDector= 500
+ %SharcQQQ
+ Z= 70
R= 0
Phi= 180
- ThicknessDector= 100
+ ThicknessDector= 500
%SharcQQQ
- Z= 66
+ Z= 70
R= 0
Phi= 270
- ThicknessDector= 100
+ ThicknessDector= 500
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index d453d5fa3..a4fedfb40 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -13,13 +13,13 @@
ExcitationEnergy= 1.4
Beam
- Particle= 7Li
- Energy= 20
- SigmaEnergy= 20
- SigmaThetaX= 25
- SigmaPhiY= 25
- SigmaX= 5
- SigmaY= 5
+ Particle= 11Li
+ Energy= 550
+ SigmaEnergy= 0
+ SigmaThetaX= 1
+ SigmaPhiY= 1
+ SigmaX= 1
+ SigmaY= 1
MeanThetaX= 0
MeanPhiY= 0
MeanX= 0
@@ -31,19 +31,19 @@ Beam
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TwoBodyReaction
- Beam= 9Li
+ Beam= 11Li
Target= 2H
Light= 3He
- Heavy= 8He
+ Heavy= 10He
ExcitationEnergyLight= 0.0
- ExcitationEnergyHeavy= 0.0
+ ExcitationEnergyHeavy= 3.0
CrossSectionPath= flat.txt CS10He
- ShootLight= 0
+ ShootLight= 1
ShootHeavy= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%ParticleDecay 10He
+ParticleDecay 10He
Daughter= 9He n
ExcitationEnergy= 0 0
DifferentialCrossSection= 11Li(d,3He)10He.txt particle9He
@@ -51,9 +51,9 @@ TwoBodyReaction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GammaDecay 10He
+GammaDecay 10He
Cascade
- BranchingRatio= 100
+ BranchingRatio= 30
Energies= 0.1
DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma9He
Cascade
@@ -62,7 +62,7 @@ TwoBodyReaction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%ParticleDecay 9He
+ParticleDecay 9He
Daughter= 8He n
DifferentialCrossSection= flat.txt
shoot= 1 1
diff --git a/NPAnalysis/pMakePhysicalTree/DataProcessor.cxx b/NPAnalysis/pMakePhysicalTree/DataProcessor.cxx
index dd2769ae6..e005c9956 100644
--- a/NPAnalysis/pMakePhysicalTree/DataProcessor.cxx
+++ b/NPAnalysis/pMakePhysicalTree/DataProcessor.cxx
@@ -6,7 +6,7 @@ ClassImp(DataProcessor);
DataProcessor::DataProcessor(){
m_InputChain = 0 ;
m_OutputTree = 0 ;
- m_ProofFile = 0 ;
+ m_ProofFile = 0 ;
m_OutputFile = 0 ;
}
@@ -17,8 +17,6 @@ DataProcessor::~DataProcessor(){
//_____________________________________________________________________________
void DataProcessor::SlaveBegin(TTree*){
-
-
TString NPLPath = gSystem->Getenv("NPTOOL");
gROOT->ProcessLine(Form(".x %s/NPLib/scripts/NPToolLogon.C+", NPLPath.Data()));
@@ -47,9 +45,9 @@ void DataProcessor::SlaveBegin(TTree*){
m_OutputTree = RootOutput::getInstance()->GetTree();
//Merging via file:
- UInt_t opt = TProofOutputFile::kRegister | TProofOutputFile::kOverwrite | TProofOutputFile::kVerify;
- m_ProofFile = new TProofOutputFile("Local.root",TProofOutputFile::kDataset, opt );
- m_OutputFile = m_ProofFile->OpenFile("RECREATE");
+ UInt_t opt = TProofOutputFile::kRegister | TProofOutputFile::kOverwrite | TProofOutputFile::kVerify;
+ m_ProofFile = new TProofOutputFile("Local.root",TProofOutputFile::kDataset, opt );
+ m_OutputFile = m_ProofFile->OpenFile("RECREATE");
m_OutputTree->SetDirectory(m_OutputFile);
m_OutputTree->AutoSave();
@@ -98,7 +96,7 @@ void DataProcessor::SlaveTerminate(){
m_OutputFile->cd();
m_OutputTree->Write();
m_OutputFile->Flush();
- GetOutputList()->Add(m_ProofFile);
+ GetOutputList()->Add(m_ProofFile);
m_OutputFile->Close();
RootOutput::getInstance()->Destroy();
RootInput::getInstance()->Destroy();
diff --git a/NPAnalysis/pMakePhysicalTree/RunToTreat.txt b/NPAnalysis/pMakePhysicalTree/RunToTreat.txt
index 4d110fe07..c77677c0e 100644
--- a/NPAnalysis/pMakePhysicalTree/RunToTreat.txt
+++ b/NPAnalysis/pMakePhysicalTree/RunToTreat.txt
@@ -1,4 +1,4 @@
TTreeName
S1107Data
RootFileName
- /Volumes/S1107/RootData/data2253*.root
+ /Volumes/S1107/RootData/data2*.root
diff --git a/NPLib/InitialConditions/TInitialConditions.cxx b/NPLib/InitialConditions/TInitialConditions.cxx
index 8baf85729..774247afb 100644
--- a/NPLib/InitialConditions/TInitialConditions.cxx
+++ b/NPLib/InitialConditions/TInitialConditions.cxx
@@ -113,16 +113,6 @@ TVector3 TInitialConditions::GetParticleDirection (const int &i) const {
}
-double TInitialConditions::GetThetaLab_WorldFrame (const int &i) const {
- return (GetParticleDirection(i).Angle(TVector3(0,0,1)))/deg;
-
-}
-
-double TInitialConditions::GetThetaLab_IncidentFrame (const int &i) const{
- return (GetParticleDirection(i).Angle(GetBeamDirection()))/deg;
-}
-
-
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 2e4eb9cd0..5d059424b 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -124,8 +124,14 @@ public:
TVector3 GetBeamDirection () const ;
TVector3 GetParticleDirection (const int &i) const ;
- double GetThetaLab_WorldFrame (const int &i) const ;
- double GetThetaLab_IncidentFrame (const int &i) const ;
+
+ double GetThetaLab_WorldFrame (const int &i) const {
+ return (GetParticleDirection(i).Angle(TVector3(0,0,1)))/deg;
+ }
+
+ double GetThetaLab_IncidentFrame (const int &i) const{
+ return (GetParticleDirection(i).Angle(GetBeamDirection()))/deg;
+ }
unsigned int GetEmittedMult() const {return fIC_Particle_Name.size();}
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index 1269606f7..5cc783910 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -78,6 +78,38 @@ Beam::Beam(){
fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+Beam::Beam(string isotope){
+ SetUp(isotope);
+ fEnergy = 0;
+ fSigmaEnergy = -1 ;
+ fMeanX = 0 ;
+ fMeanY = 0 ;
+ fSigmaX = -1;
+ fSigmaY = -0;
+ fMeanThetaX = 0 ;
+ fMeanPhiY = 0 ;
+ fSigmaThetaX = -1 ;
+ fSigmaPhiY = -1 ;
+ fTargetSize = 0 ;
+ fEffectiveTargetSize = 0 ;
+ fTargetThickness = 0 ;
+ fEffectiveTargetThickness = 0 ;
+ fTargetAngle = 0 ;
+ fTargetZ = 0 ;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+
+ // case of user given distribution
+ // do that to avoid warning from multiple Hist with same name...
+ int offset = 0;
+ while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
+ ++offset;
+
+ fEnergyHist = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
+ fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
+ fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Beam::~Beam(){
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index 8c1659eb1..1c712f01f 100755
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -45,6 +45,7 @@ namespace NPL{
public: // Constructors and Destructors
Beam();
+ Beam(string);
~Beam();
public: // Various Method
diff --git a/NPLib/Physics/NPNucleus.cxx b/NPLib/Physics/NPNucleus.cxx
index e16ebd175..89c5a9031 100644
--- a/NPLib/Physics/NPNucleus.cxx
+++ b/NPLib/Physics/NPNucleus.cxx
@@ -69,6 +69,12 @@ void Nucleus::SetUp(string isotope){
//----------- Constructor Using nubtab03.asc by name----------
// open the file to read and check if it is open
+ // Replace the p,d,t,a by there standard name:
+ if(isotope=="p") isotope="1H";
+ else if(isotope=="d") isotope="2H";
+ else if(isotope=="t") isotope="3H";
+ else if(isotope=="a") isotope="4He";
+
const char* Isotope = isotope.c_str();
ifstream inFile;
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index 3a153d3bc..f1086e9bc 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -58,6 +58,11 @@ ClassImp(Reaction)
Reaction::Reaction(){
//------------- Default Constructor -------------
+ // Need to be done before initializePrecomputeVariable
+ fKineLine3 = 0 ;
+ fKineLine4 = 0 ;
+
+ //
fNuclei1 = new Beam();
fNuclei2 = new Nucleus();
fNuclei3 = new Nucleus();
@@ -78,8 +83,71 @@ Reaction::Reaction(){
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
fshoot3=true;
fshoot4=true;
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// This constructor aim to provide a fast way to instantiate a reaction without input file
+// The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
+Reaction::Reaction(string reaction){
+ // Instantiate the parameter to default
+ // Analyse the reaction and extract revelant information
+ string A,b,c,D,E;
+ unsigned int i=0;
+ for(; i < reaction.length() ; i++){
+ if(reaction.compare(i,1,"(")!=0) A.push_back(reaction[i]);
+ else break;
+ }
+
+ i++;
+ for(; i < reaction.length() ; i++){
+ if(reaction.compare(i,1,",")!=0) b.push_back(reaction[i]);
+ else break;
+ }
+
+ i++;
+ for(; i < reaction.length() ; i++){
+ if(reaction.compare(i,1,")")!=0) c.push_back(reaction[i]);
+ else break;
+ }
+
+ i++;
+ for(; i < reaction.length() ; i++){
+ if(reaction.compare(i,1,"@")!=0) D.push_back(reaction[i]);
+ else break;
+ }
+
+ i++;
+ for(; i < reaction.length() ; i++){
+ E.push_back(reaction[i]);
+ }
+
+ fKineLine3 = 0 ;
+ fKineLine4 = 0 ;
+ fNuclei1 = new Beam(A);
+ fNuclei2 = new Nucleus(b);
+ fNuclei3 = new Nucleus(c);
+ fNuclei4 = new Nucleus(D);
+ fBeamEnergy = atof(E.c_str());
+ fThetaCM = 0;
+ fExcitation3 = 0;
+ fExcitation4 = 0;
+ fQValue = 0;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ initializePrecomputeVariable();
+
+ // do that to avoid warning from multiple Hist with same name... int offset = 0;
+ int offset = 0;
+ while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
+ ++offset;
+
+ fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
+ fshoot3=true;
+ fshoot4=true;
+
+
+
+ initializePrecomputeVariable();
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Reaction::~Reaction(){
//------------- Default Destructor ------------
@@ -363,6 +431,10 @@ void Reaction::ReadConfigurationFile(string Path){
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::initializePrecomputeVariable(){
+ // delete the previously calculated kinematical line:
+ if(fKineLine3!=0) {delete fKineLine3 ; fKineLine3 = 0;}
+ if(fKineLine4!=0) {delete fKineLine4 ; fKineLine4 = 0;}
+
m1 = fNuclei1->Mass();
m2 = fNuclei2->Mass();
m3 = fNuclei3->Mass() + fExcitation3;
@@ -416,46 +488,51 @@ void Reaction::SetNuclei3(double EnergyLab, double ThetaLab){
////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetKinematicLine3(){
- int size = 360;
- double x[size];
- double y[size];
- double theta3,E3,theta4,E4;
-
- for (int i = 0; i < size; ++i){
- SetThetaCM(((double)i)/2*deg);
- KineRelativistic(theta3, E3, theta4, E4);
- fNuclei3->SetKineticEnergy(E3);
-
- x[i] = theta3/deg;
- y[i] = E3;
- }
- TGraph* KineLine3 = new TGraph(size,x,y);
- KineLine3->SetTitle("Kinematic Line of particule 3");
- return(KineLine3);
+ if(fKineLine3==0){
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4;
+
+ for (int i = 0; i < size; ++i){
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei3->SetKineticEnergy(E3);
+
+ x[i] = theta3/deg;
+ y[i] = E3;
+ }
+
+ fKineLine3 = new TGraph(size,x,y);
+ // fKineLine3->SetTitle("Kinematic Line of particule 3");
+ }
+ return(fKineLine3);
}
////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetKinematicLine4(){
- int size = 360;
- double x[size];
- double y[size];
- double theta3,E3,theta4,E4;
-
- for (int i = 0; i < size; ++i)
+ if(fKineLine4==0){
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4;
+
+ for (int i = 0; i < size; ++i)
{
- SetThetaCM(((double)i)/2*deg);
- KineRelativistic(theta3, E3, theta4, E4);
- fNuclei4->SetKineticEnergy(E4);
-
- x[i] = theta4/deg;
- y[i] = E4;
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei4->SetKineticEnergy(E4);
+
+ x[i] = theta4/deg;
+ y[i] = E4;
}
- TGraph* KineLine4= new TGraph(size,x,y);
- KineLine4->SetTitle("Kinematic Line of particule 4");
- return(KineLine4);
+ fKineLine4= new TGraph(size,x,y);
+ // fKineLine4->SetTitle("Kinematic Line of particule 4");
+ }
+ return(fKineLine4);
}
////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 969e6ba3e..f419e9595 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -58,6 +58,9 @@ namespace NPL{
public: // Constructors and Destructors
Reaction();
+ // This constructor aim to provide a fast way to instantiate a reaction without input file
+ // The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
+ Reaction(string);
virtual ~Reaction();
public: // Various Method
@@ -70,11 +73,14 @@ namespace NPL{
bool fshoot3;
bool fshoot4;
+ private: // use to display the kinematical line
+ TGraph* fKineLine3 ;
+ TGraph* fKineLine4 ;
private:
- Beam *fNuclei1; // Beam
- Nucleus *fNuclei2; // Target
- Nucleus *fNuclei3; // Light ejectile
- Nucleus *fNuclei4; // Heavy ejectile
+ Beam* fNuclei1; // Beam
+ Nucleus* fNuclei2; // Target
+ Nucleus* fNuclei3; // Light ejectile
+ Nucleus* fNuclei4; // Heavy ejectile
double fQValue; // Q-value in MeV
double fBeamEnergy; // Beam energy in MeV
double fThetaCM; // Center-of-mass angle in radian
diff --git a/NPLib/Sharc/TSharcPhysics.cxx b/NPLib/Sharc/TSharcPhysics.cxx
index d2d76932d..040b268ca 100644
--- a/NPLib/Sharc/TSharcPhysics.cxx
+++ b/NPLib/Sharc/TSharcPhysics.cxx
@@ -695,14 +695,14 @@ void TSharcPhysics::AddBoxDetector(double Z)
m_NumberOfDetector++;
if(Z<0){// Up Stream
- if(i==0) {U=TVector3(1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-36.,42.5,-56.) ;}
+ if(i==0) {U=TVector3(1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-36.,42.5,-56.) ;}
else if(i==1) {U=TVector3(0,1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-42.5,-36.,-56.) ;}
else if(i==2) {U=TVector3(-1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(36.,-42.5,-56.) ;}
else if(i==3) {U=TVector3(0,-1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(42.5,36.,-56.) ;}
}
if(Z>0){//Down Stream
- if(i==0) {U=TVector3(-1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(36.,40.5,60.) ;}
+ if(i==0) {U=TVector3(-1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(36.,40.5,60.) ;}
else if(i==1) {U=TVector3(0,-1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-40.5,36.,60.) ;}
else if(i==2) {U=TVector3(1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-36.,-40.5,60.) ;}
else if(i==3) {U=TVector3(0,1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(40.5,-36.,60.) ;}
@@ -742,13 +742,12 @@ void TSharcPhysics::AddBoxDetector(double Z)
void TSharcPhysics::AddQQQDetector( double R,double Phi,double Z){
- double QQQ_R_Min = 9.+R-R;
- double QQQ_R_Max = 41.0;
+ double QQQ_R_Min = 9.+R;
+ double QQQ_R_Max = 41.0+R;
double QQQ_Phi_Min = 2.0*M_PI/180. ;
double QQQ_Phi_Max = 83.6*M_PI/180. ;
Phi= Phi*M_PI/180.;
- Z= -63.5;
int QQQ_Radial_NumberOfStrip = 16 ;
int QQQ_Sector_NumberOfStrip = 24 ;
@@ -776,7 +775,7 @@ void TSharcPhysics::AddQQQDetector( double R,double Phi,double Z){
for(int b = 0 ; b < QQQ_Sector_NumberOfStrip ; b++){
StripCenter = Strip_1_1;
- StripCenter.SetY(QQQ_R_Max-StripPitchRadial*f);
+ StripCenter.SetY(QQQ_R_Max-f*StripPitchRadial);
StripCenter.SetZ(Z);
StripCenter.RotateZ(Phi+QQQ_Phi_Min+b*StripPitchSector);
lineX.push_back( StripCenter.X() );
diff --git a/NPSimulation/include/Sharc.hh b/NPSimulation/include/Sharc.hh
index 8a8fafe3b..c2f14653e 100644
--- a/NPSimulation/include/Sharc.hh
+++ b/NPSimulation/include/Sharc.hh
@@ -46,7 +46,10 @@ namespace SHARC
// Energy and time Resolution
const G4double ResoTime = 0 ;
const G4double ResoEnergy = 0.035*MeV ;// = zzkeV of Resolution // Unit is MeV/2.35
- const G4double EnergyThreshold = 0.4*MeV;
+ //const G4double EnergyThreshold = 0.4*MeV;
+
+ // Change for TRex simulations
+ const G4double EnergyThreshold = 0.5*MeV;
// Geometry
// BOX //
@@ -124,7 +127,7 @@ namespace SHARC
// QQQ Wafer
const G4double QQQ_Wafer_Outer_Radius = 42.63*mm;
- const G4double QQQ_Wafer_Inner_Radius = 7.8*mm;
+ const G4double QQQ_Wafer_Inner_Radius = 9.0*mm;
const G4double QQQ_Wafer_Starting_Phi = 8*deg;
const G4double QQQ_Wafer_Stopping_Phi = 162*deg;
const G4int QQQ_Wafer_NumberOf_RadialStrip = 16 ;
diff --git a/NPSimulation/include/SharcScorers.hh b/NPSimulation/include/SharcScorers.hh
index 6bae9a21d..08c59c9e9 100644
--- a/NPSimulation/include/SharcScorers.hh
+++ b/NPSimulation/include/SharcScorers.hh
@@ -26,15 +26,17 @@
*****************************************************************************/
#include "G4VPrimitiveScorer.hh"
#include "G4THitsMap.hh"
+
+#include <map>
+using namespace std;
+
namespace SHARC {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- // This Threshold is used in the above scorer. Any energy deposit under this threshold will not create an entry.
- const double TriggerThreshold = 0.1*keV ;
class PS_Silicon_Rectangle : public G4VPrimitiveScorer{
public: // with description
- PS_Silicon_Rectangle(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4double TriggerThreshold,G4int depth=0);
+ PS_Silicon_Rectangle(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4int depth=0);
~PS_Silicon_Rectangle();
protected: // with description
@@ -74,7 +76,7 @@ namespace SHARC {
class PS_Silicon_Annular : public G4VPrimitiveScorer{
public: // with description
- PS_Silicon_Annular(G4String name, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double DeltaTheta, G4int NumberOfStripRadial,G4int NumberOfStripTheta,G4double TriggerThreshold,G4int depth=0);
+ PS_Silicon_Annular(G4String name, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double PhiStart,G4double PhiStop, G4int NumberOfStripRadial,G4int NumberOfStripTheta,G4int depth=0);
~PS_Silicon_Annular();
protected: // with description
@@ -93,7 +95,8 @@ namespace SHARC {
private: // Geometry of the detector
G4double m_StripPlaneInnerRadius;
G4double m_StripPlaneOuterRadius;
- G4double m_DeltaTheta;
+ G4double m_PhiStart;
+ G4double m_PhiStop;
G4int m_NumberOfStripRadial;
G4int m_NumberOfStripTheta;
G4double m_StripPitchRadial;
@@ -102,13 +105,13 @@ namespace SHARC {
private: // inherited from G4VPrimitiveScorer
G4int HCID;
G4THitsMap<G4double*>* EvtMap;
-
+
private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
G4ThreeVector m_Position ;
- G4ThreeVector m_uz;
+ G4ThreeVector m_uz ;
G4int m_DetectorNumber ;
G4int m_StripRadialNumber ;
- G4int m_StripThetaNumber ;
+ G4int m_StripThetaNumber ;
G4int m_Index ;
};
diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc
index 2e4531c59..42d212417 100644
--- a/NPSimulation/src/DetectorConstruction.cc
+++ b/NPSimulation/src/DetectorConstruction.cc
@@ -63,6 +63,7 @@
#include "ThinSi.hh"
#include "Sharc.hh"
#include "Shield.hh"
+#include "Tigress.hh"
#include "W1.hh"
// STL
@@ -187,7 +188,8 @@ void DetectorConstruction::ReadConfigurationFile(string Path){
bool cParis = false;
bool cS1 = false;
bool cSharc = false;
- bool cShield = false;
+ bool cShield = false;
+ bool cTigress = false;
bool cW1 = false;
bool cHelios = false;
@@ -463,9 +465,30 @@ void DetectorConstruction::ReadConfigurationFile(string Path){
AddDetector(myDetector);
#endif
}
+
+ ////////////////////////////////////////////
+ ////////// Search for Tigress ///////////
+ ////////////////////////////////////////////
+ else if (LineBuffer.compare(0,7, "Tigress") == 0 && cTigress == false) {
+#ifdef INC_TIGRESS
+ cTigress = true ;
+ if(VerboseLevel==1) cout << endl << "//////// Tigress ////////" << endl << endl ;
+
+ // Instantiate the new array as a VDetector Object
+ VDetector* myDetector = new Tigress();
+
+ // Read Position of detector
+ ConfigFile.close();
+ myDetector->ReadConfiguration(Path);
+ ConfigFile.open(Path.c_str());
+
+ // Add array to the VDetector Vector
+ AddDetector(myDetector);
+#endif
+ }
////////////////////////////////////////////
- ////////// Search for Plastic ///////////
+ ////////// Search for Plastic ///////////
////////////////////////////////////////////
else if (LineBuffer.compare(0, 19, "ScintillatorPlastic") == 0 && cPlastic == false) {
#ifdef INC_PLASTIC
diff --git a/NPSimulation/src/Sharc.cc b/NPSimulation/src/Sharc.cc
index c1d34f433..09e0ec913 100644
--- a/NPSimulation/src/Sharc.cc
+++ b/NPSimulation/src/Sharc.cc
@@ -200,7 +200,7 @@ void Sharc::ReadConfiguration(string Path){
check_Thickness = true;
ConfigFile >> DataBuffer ;
Thickness= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessDetector= " << Thickness/um << "mm" << endl;
+ if(VerboseLevel==1) cout << " ThicknessDetector= " << Thickness/um << "um" << endl;
}
///////////////////////////////////////////////////
@@ -245,56 +245,56 @@ void Sharc::ReadConfiguration(string Path){
check_Thickness1 = true;
ConfigFile >> DataBuffer ;
Thickness1= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessDetector1= " << Thickness1/um << "mm" << endl;
+ if(VerboseLevel==1) cout << " ThicknessDetector1= " << Thickness1/um << "um" << endl;
}
else if (DataBuffer == "ThicknessDector2=") {
check_Thickness2 = true;
ConfigFile >> DataBuffer ;
Thickness2= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessDetector2= " << Thickness2/um << "mm" << endl;
+ if(VerboseLevel==1) cout << " ThicknessDetector2= " << Thickness2/um << "um" << endl;
}
else if (DataBuffer == "ThicknessDector3=") {
check_Thickness3 = true;
ConfigFile >> DataBuffer ;
Thickness3= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessDetector3= " << Thickness3/um << "mm" << endl;
+ if(VerboseLevel==1) cout << " ThicknessDetector3= " << Thickness3/um << "um" << endl;
}
else if (DataBuffer == "ThicknessDector4=") {
check_Thickness4 = true;
ConfigFile >> DataBuffer ;
Thickness4= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessDetector4= " << Thickness4/um << "mm" << endl;
+ if(VerboseLevel==1) cout << " ThicknessDetector4= " << Thickness4/um << "um" << endl;
}
else if (DataBuffer == "ThicknessPAD1=") {
check_PAD1 = true;
ConfigFile >> DataBuffer ;
ThicknessPAD1= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessPAD1= " << ThicknessPAD1 << endl;
+ if(VerboseLevel==1) cout << " ThicknessPAD1= " << ThicknessPAD1<< "um" << endl;
}
else if (DataBuffer == "ThicknessPAD2=") {
check_PAD2 = true;
ConfigFile >> DataBuffer ;
ThicknessPAD2= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessPAD2= " << ThicknessPAD2 << endl;
+ if(VerboseLevel==1) cout << " ThicknessPAD2= " << ThicknessPAD2<< "um" << endl;
}
else if (DataBuffer == "ThicknessPAD3=") {
check_PAD3 = true;
ConfigFile >> DataBuffer ;
ThicknessPAD3= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessPAD3= " << ThicknessPAD3 << endl;
+ if(VerboseLevel==1) cout << " ThicknessPAD3= " << ThicknessPAD3<< "um" << endl;
}
else if (DataBuffer == "ThicknessPAD4=") {
check_PAD4 = true;
ConfigFile >> DataBuffer ;
ThicknessPAD4= atof(DataBuffer.c_str())*um;
- if(VerboseLevel==1) cout << " ThicknessPAD4= " << ThicknessPAD4 << endl;
+ if(VerboseLevel==1) cout << " ThicknessPAD4= " << ThicknessPAD4<< "um" << endl;
}
///////////////////////////////////////////////////
@@ -633,8 +633,8 @@ void Sharc::ReadSensitive(const G4Event* event){
///////////
// BOX
- G4THitsMap<G4double*>* BOXHitMap;
- std::map<G4int, G4double**>::iterator BOX_itr;
+ G4THitsMap<G4double*>* BOXHitMap;
+ std::map<G4int, G4double**>::iterator BOX_itr;
G4int BOXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Sharc_BOXScorer/SharcBOX");
BOXHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(BOXCollectionID));
@@ -644,34 +644,35 @@ void Sharc::ReadSensitive(const G4Event* event){
G4double* Info = *(BOX_itr->second);
- double Energy = Info[0];
- double Time = Info[1];
- int DetNbr = (int) Info[2];
- int StripFront = (int) Info[3];
- int StripBack = (int) Info[4];
-
- m_Event->SetFront_DetectorNbr(DetNbr);
- m_Event->SetFront_StripNbr(StripFront);
- m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
- m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
- m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
-
- m_Event->SetBack_DetectorNbr(DetNbr);
- m_Event->SetBack_StripNbr(StripBack);
- m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
- m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
- m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
-
-
- // Interraction Coordinates
- ms_InterCoord->SetDetectedPositionX(Info[5]) ;
- ms_InterCoord->SetDetectedPositionY(Info[6]) ;
- ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
- ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
- ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
+ double Energy = Info[0];
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+ int StripFront = (int) Info[3];
+ int StripBack = (int) Info[4];
+
+ m_Event->SetFront_DetectorNbr(DetNbr);
+ m_Event->SetFront_StripNbr(StripFront);
+ m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ m_Event->SetBack_DetectorNbr(DetNbr);
+ m_Event->SetBack_StripNbr(StripBack);
+ m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ // Interraction Coordinates
+ ms_InterCoord->SetDetectedPositionX(Info[5]) ;
+ ms_InterCoord->SetDetectedPositionY(Info[6]) ;
+ ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
+ ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
+ ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
+
+ }
}
-
// clear map for next event
BOXHitMap->clear();
@@ -689,14 +690,15 @@ void Sharc::ReadSensitive(const G4Event* event){
G4double* Info = *(PAD_itr->second);
double Energy = Info[0];
- double Time = Info[1];
- int DetNbr = (int) Info[2];
-
- m_Event->SetPAD_DetectorNbr(DetNbr);
- m_Event->SetPAD_Energy(RandGauss::shoot(Energy, ResoEnergy));
- m_Event->SetPAD_TimeCFD(RandGauss::shoot(Time, ResoTime));
- m_Event->SetPAD_TimeLED(RandGauss::shoot(Time, ResoTime));
-
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+
+ m_Event->SetPAD_DetectorNbr(DetNbr);
+ m_Event->SetPAD_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetPAD_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetPAD_TimeLED(RandGauss::shoot(Time, ResoTime));
+ }
}
// clear map for next event
@@ -704,8 +706,8 @@ void Sharc::ReadSensitive(const G4Event* event){
///////////
// QQQ
- G4THitsMap<G4double*>* QQQHitMap;
- std::map<G4int, G4double**>::iterator QQQ_itr;
+ G4THitsMap<G4double*>* QQQHitMap;
+ std::map<G4int, G4double**>::iterator QQQ_itr;
G4int QQQCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Sharc_QQQScorer/SharcQQQ");
QQQHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(QQQCollectionID));
@@ -716,30 +718,31 @@ void Sharc::ReadSensitive(const G4Event* event){
G4double* Info = *(QQQ_itr->second);
double Energy = Info[0];
- double Time = Info[1];
- int DetNbr = (int) Info[2];
- int StripFront = (int) Info[3];
- int StripBack = (int) Info[4];
-
- m_Event->SetFront_DetectorNbr(DetNbr);
- m_Event->SetFront_StripNbr(StripFront);
- m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
- m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
- m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
-
- m_Event->SetBack_DetectorNbr(DetNbr);
- m_Event->SetBack_StripNbr(StripBack);
- m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
- m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
- m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
-
- // Interraction Coordinates
- ms_InterCoord->SetDetectedPositionX(Info[5]) ;
- ms_InterCoord->SetDetectedPositionY(Info[6]) ;
- ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
- ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
- ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
-
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+ int StripFront = (int) Info[3];
+ int StripBack = (int) Info[4];
+
+ m_Event->SetFront_DetectorNbr(DetNbr);
+ m_Event->SetFront_StripNbr(StripFront);
+ m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ m_Event->SetBack_DetectorNbr(DetNbr);
+ m_Event->SetBack_StripNbr(StripBack);
+ m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ // Interraction Coordinates
+ ms_InterCoord->SetDetectedPositionX(Info[5]) ;
+ ms_InterCoord->SetDetectedPositionY(Info[6]) ;
+ ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
+ ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
+ ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
+ }
}
// clear map for next event
@@ -760,27 +763,23 @@ void Sharc::InitializeScorers(){
BOX_Wafer_Length,
BOX_Wafer_Width,
BOX_Wafer_Back_NumberOfStrip ,
- BOX_Wafer_Front_NumberOfStrip,
- EnergyThreshold);
+ BOX_Wafer_Front_NumberOfStrip);
G4VPrimitiveScorer* PADScorer =
new SHARC::PS_Silicon_Rectangle("SharcPAD",
PAD_Wafer_Length,
PAD_Wafer_Width,
1 ,
- 1,
- EnergyThreshold);
+ 1);
G4VPrimitiveScorer* QQQScorer =
new SHARC::PS_Silicon_Annular("SharcQQQ",
QQQ_Wafer_Inner_Radius,
QQQ_Wafer_Outer_Radius,
- QQQ_Wafer_Stopping_Phi-QQQ_Wafer_Starting_Phi,
+ QQQ_Wafer_Starting_Phi,
+ QQQ_Wafer_Stopping_Phi,
QQQ_Wafer_NumberOf_RadialStrip,
- QQQ_Wafer_NumberOf_AnnularStrip,
- EnergyThreshold);
-
-
+ QQQ_Wafer_NumberOf_AnnularStrip);
//and register it to the multifunctionnal detector
m_BOXScorer->RegisterPrimitive(BOXScorer);
diff --git a/NPSimulation/src/SharcScorers.cc b/NPSimulation/src/SharcScorers.cc
index 8cda2113c..62612c224 100644
--- a/NPSimulation/src/SharcScorers.cc
+++ b/NPSimulation/src/SharcScorers.cc
@@ -24,8 +24,8 @@
using namespace SHARC ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_Silicon_Rectangle::PS_Silicon_Rectangle(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4double TriggerThreshold,G4int depth)
-:G4VPrimitiveScorer(name, depth), m_TriggerThreshold(TriggerThreshold),HCID(-1){
+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;
@@ -51,11 +51,7 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
G4double* EnergyAndTime = new G4double[10];
EnergyAndTime[0] = aStep->GetTotalEnergyDeposit();
- // If the energy deposit is below the threshold, the deposit is ignored
- if (EnergyAndTime[0] < m_TriggerThreshold){
- delete EnergyAndTime;
- return FALSE;
- }
+
EnergyAndTime[1] = aStep->GetPreStepPoint()->GetGlobalTime();
@@ -74,6 +70,10 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
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 ;
+
+
EnergyAndTime[2] = m_DetectorNumber;
EnergyAndTime[3] = m_StripLengthNumber;
EnergyAndTime[4] = m_StripWidthNumber;
@@ -84,6 +84,15 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
if (m_StripWidthNumber == m_NumberOfStripWidth+1) m_StripWidthNumber = 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);
+ EnergyAndTime[0]+=dummy[0];
+ }
+
EvtMap->set(m_Index, EnergyAndTime);
return TRUE;
}
@@ -126,15 +135,16 @@ void PS_Silicon_Rectangle::PrintAll(){
//....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 DeltaTheta, G4int NumberOfStripRadial,G4int NumberOfStripTheta,G4double TriggerThreshold,G4int depth)
-:G4VPrimitiveScorer(name, depth), m_TriggerThreshold(TriggerThreshold),HCID(-1){
+PS_Silicon_Annular::PS_Silicon_Annular(G4String name, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double PhiStart,G4double PhiStop, G4int NumberOfStripRadial,G4int NumberOfStripTheta,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
m_StripPlaneInnerRadius = StripPlaneInnerRadius;
m_StripPlaneOuterRadius = StripPlaneOuterRadius;
- m_DeltaTheta = DeltaTheta;
+ m_PhiStart = PhiStart;
+ m_PhiStop = PhiStop;
m_NumberOfStripRadial = NumberOfStripRadial;
m_NumberOfStripTheta = NumberOfStripTheta;
- m_StripPitchRadial = DeltaTheta/m_NumberOfStripTheta;
+ m_StripPitchRadial = (m_PhiStop-m_PhiStart)/m_NumberOfStripTheta;
m_StripPitchTheta = (m_StripPlaneOuterRadius-m_StripPlaneInnerRadius)/m_NumberOfStripRadial;
m_Position = G4ThreeVector(-1000,-1000,-1000);
@@ -151,33 +161,26 @@ 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* EnergyAndTime = new G4double[10];
EnergyAndTime[0] = aStep->GetTotalEnergyDeposit();
- // If the energy deposit is below the threshold, the deposit is ignored
- if (EnergyAndTime[0] < m_TriggerThreshold){
- delete EnergyAndTime;
- return FALSE;
- }
-
EnergyAndTime[1] = aStep->GetPreStepPoint()->GetGlobalTime();
m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0);
m_Position = aStep->GetPreStepPoint()->GetPosition();
-
+
// Interaction coordinates (used to fill the InteractionCoordinates branch)
EnergyAndTime[5] = m_Position.x();
EnergyAndTime[6] = m_Position.y();
EnergyAndTime[7] = m_Position.z();
EnergyAndTime[8] = m_Position.theta();
EnergyAndTime[9] = m_Position.phi();
-
+
m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position);
-
- m_StripRadialNumber = (int)((m_Position.angle(m_uz) - m_DeltaTheta*0.5) / m_StripPitchRadial ) + 1 ;
- m_StripThetaNumber = (int)((m_Position.rho()-m_StripPlaneInnerRadius) / m_StripPitchTheta ) + 1 ;
+ m_StripRadialNumber = (int) ((m_StripPlaneOuterRadius - m_Position.rho()) / m_StripPitchTheta ) + 1 ;
+
+ m_StripThetaNumber = (int) ((m_Position.phi() - m_PhiStart*0.5) / m_StripPitchRadial ) + 1 ;
EnergyAndTime[2] = m_DetectorNumber;
EnergyAndTime[3] = m_StripRadialNumber;
@@ -188,6 +191,15 @@ G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){
if (m_StripThetaNumber == m_NumberOfStripTheta+1) m_StripThetaNumber = m_NumberOfStripTheta;
m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripRadialNumber * 1e6 + m_NumberOfStripTheta * 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);
+ EnergyAndTime[0]+=dummy[0];
+ }
+
EvtMap->set(m_Index, EnergyAndTime);
return TRUE;
}
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index b935a3742..1ccdaa85d 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -616,7 +616,7 @@ void Target::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4doub
G4EmCalculator emCalculator;
- for (G4double E=Emin; E < Emax; E+=(Emax-Emin)/10.){
+ 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 ;
}
diff --git a/NPSimulation/vis.mac b/NPSimulation/vis.mac
index 96d2fc4a6..9150d12e3 100644
--- a/NPSimulation/vis.mac
+++ b/NPSimulation/vis.mac
@@ -40,7 +40,7 @@
/vis/viewer/set/style surface
#
# Draw coordinate axes:
-/vis/scene/add/axes 0 0 0 20 cm
+#/vis/scene/add/axes 0 0 0 20 cm
#
# Draw smooth trajectories at end of event, showing trajectory points
# as markers 2 pixels wide:
@@ -67,7 +67,7 @@
#/vis/modeling/trajectories/drawByParticleID-0/set e- blue
#
# To superimpose all of the events from a given run:
-#/vis/scene/endOfEventAction accumulate
+/vis/scene/endOfEventAction accumulate
#
# Re-establish auto refreshing and verbosity:
/vis/viewer/set/autoRefresh true
--
GitLab