diff --git a/Inputs/DetectorConfiguration/AnnularS1.detector b/Inputs/DetectorConfiguration/AnnularS1.detector
index 145d0667d91151ae2d72c5224ae92b75eb9cfccf..1a910320f40a5ebaf9471c6096553266ef0b5940 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 2ede44d27e097a2cb889741571e43ce67d79de7c..3884cf6fa9852ce4aa72d92b25cc4b42ed12ec9d 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 d453d5fa38b465a95b50d1d399cc950e0fc63787..a4fedfb4000e66f2d84374576e285bdf2c6d6ed0 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 dd2769ae600be11e9df019763b72afed5741d330..e005c995602d8984cdc8443da4f7cae81ef15047 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 4d110fe07df6045069832b0236653227109e352f..c77677c0e4a5e1295ae3f772bbe0ee372849e0e6 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 8baf85729801de3b6286872394b024687240fc00..774247afbb726e1e31b1f81ee3118ed3eef73f8f 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 2e4eb9cd0350837f5dd0bd87fab0ec7799f2f0f4..5d059424b7e7ef5eeb70c4cb744d18fa54e5cbb5 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 1269606f76b21b3478fce4a2ed7ecaf68c077489..5cc7839100563b2a8b3b762a01011b7c52456b4e 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 8c1659eb14fd4e440fb2a3981fb50d6825e5e11a..1c712f01f57d66cab29690f3f71d9437a9b891e7 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 e16ebd17540353bbd9889d8a6e9c17426378f9d0..89c5a90311ef9aa20b3dcbb1b12a1456d0e24d12 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 3a153d3bca608fc4f773eadd552d2bef69dceb1d..f1086e9bc86a6666e39230159c5e53a2301328a8 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 969e6ba3edfcb0ec13906c006b0a19e61b5b51a8..f419e9595739f53fa83bcc77aa2905b8b1ace38a 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 d2d76932d2441b73e5c3ac2fcde944ea829e3625..040b268ca76245b91cf69c91d0926f3baa3a386c 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 8a8fafe3bf25a7f985e4d27c928b2786bf0a6faf..c2f14653e963e63e086d58fa3a5b1245e01e4708 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 6bae9a21d304a677e3c6abefefc0adce5a2c25b7..08c59c9e9729e921364182517eb94e5f3508d53d 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 2e4531c593ff07f39bd8a97b18f04a2778f222bd..42d212417c1f709a2fa40d68ca345f0e48fa2263 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 c1d34f433d1f0c8c85a23ed69779eb3f6fe0eb17..09e0ec913ce1f3c1dd7fc7158ee9c4a9f635bfae 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 8cda2113c16b50ebe2c631d05118eb74f86ce8a5..62612c22430461aaf48e5c480ff1dd9fd2ce4810 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 b935a374299d5bc97d9305446d1bfeb1e76a89d2..1ccdaa85d4b3c7878ffd615de642b4c4eaca286e 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 96d2fc4a6a0d6fdee5cd2a7d6c279f383601f2d4..9150d12e36fbabcfe951c94dd7357592d0517075 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