From 88a2bc9961d5840dcf871f8ec3a88e3eb21a45f6 Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Tue, 19 Mar 2013 15:36:12 +0000
Subject: [PATCH] * Fixing bug in initial condition * Fixing
GeometricalEfficiency.C macro to fit new TInitialCondition * Adding support
for TInitialCondition in Sharc
---
Inputs/DetectorConfiguration/Sharc.detector | 28 ++++-----
Inputs/EventGenerator/proton.source | 6 +-
NPAnalysis/macros/GeometricalEfficiency.C | 30 +++++++---
.../InitialConditions/TInitialConditions.cxx | 28 +++++++++
NPLib/InitialConditions/TInitialConditions.h | 59 +++++++++++--------
.../TInteractionCoordinates.h | 10 ++--
NPSimulation/include/Sharc.hh | 1 -
NPSimulation/src/ParticleStack.cc | 1 +
NPSimulation/src/Sharc.cc | 20 +++++++
NPSimulation/src/SharcScorers.cc | 31 +++++++++-
NPSimulation/vis.mac | 2 +-
11 files changed, 159 insertions(+), 57 deletions(-)
diff --git a/Inputs/DetectorConfiguration/Sharc.detector b/Inputs/DetectorConfiguration/Sharc.detector
index 9c8290941..2ede44d27 100644
--- a/Inputs/DetectorConfiguration/Sharc.detector
+++ b/Inputs/DetectorConfiguration/Sharc.detector
@@ -15,29 +15,29 @@ Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream CD
SharcQQQ
- Z= -100
+ Z= -66
R= 0
Phi= 0
ThicknessDector= 100
SharcQQQ
- Z= -100
+ Z= -66
R= 0
Phi= 90
ThicknessDector= 100
SharcQQQ
- Z= -100
+ Z= -66
R= 0
Phi= 180
ThicknessDector= 100
SharcQQQ
- Z= -100
+ Z= -66
R= 0
Phi= 270
ThicknessDector= 100
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream Box
SharcBOX
- Z= -40
+ Z= -30.775
ThicknessDector1= 100
ThicknessDector2= 100
ThicknessDector3= 100
@@ -49,7 +49,7 @@ Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Down Stream Box
SharcBOX
- Z= 40
+ Z= 34.3
ThicknessDector1= 100
ThicknessDector2= 100
ThicknessDector3= 100
@@ -60,23 +60,23 @@ Sharc
ThicknessPAD4= 1000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Downstream CD
- SharcQQQ
- Z= 100
+ %SharcQQQ
+ Z= 66
R= 0
Phi= 0
ThicknessDector= 100
- SharcQQQ
- Z= 100
+ %SharcQQQ
+ Z= 66
R= 0
Phi= 90
ThicknessDector= 100
- SharcQQQ
- Z= 100
+ % SharcQQQ
+ Z= 66
R= 0
Phi= 180
ThicknessDector= 100
- SharcQQQ
- Z= 100
+ %SharcQQQ
+ Z= 66
R= 0
Phi= 270
ThicknessDector= 100
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index ca319dcce..b6a19c47f 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 300
- EnergyHigh= 300
+ EnergyLow= 5
+ EnergyHigh= 5
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 30
+ HalfOpenAngleMax= 180
x0= 0
y0= 0
z0= 0
diff --git a/NPAnalysis/macros/GeometricalEfficiency.C b/NPAnalysis/macros/GeometricalEfficiency.C
index 8c8023369..34d552501 100644
--- a/NPAnalysis/macros/GeometricalEfficiency.C
+++ b/NPAnalysis/macros/GeometricalEfficiency.C
@@ -27,21 +27,23 @@
*****************************************************************************/
#include <iostream>
+#include <cmath>
#include "TROOT.h"
#include "TSystem.h"
#include "TFile.h"
+#include "TF1.h"
#include "TString.h"
#include "TTree.h"
#include "TBranch.h"
#include "TH1F.h"
+#include "TCanvas.h"
#include "TInitialConditions.h"
#include "TInteractionCoordinates.h"
using namespace std ;
-
void GeometricalEfficiency(const char * fname = "myResult")
{
// Open output ROOT file from NPTool simulation run
@@ -56,11 +58,11 @@ void GeometricalEfficiency(const char * fname = "myResult")
// TInitialConditions branch
TInitialConditions *initCond = 0;
tree->SetBranchAddress("InitialConditions", &initCond);
- tree->SetBranchStatus("InitialConditions", 1);
+ tree->SetBranchStatus("InitialConditions", true);
// TInteractionCoordinates branch
TInteractionCoordinates *interCoord = 0;
tree->SetBranchAddress("InteractionCoordinates", &interCoord);
- tree->SetBranchStatus("InteractionCoordinates", 1);
+ tree->SetBranchStatus("InteractionCoordinates", true);
// Prepare histograms
TH1F *hDetecTheta = new TH1F("hDetecTheta", "DetecTheta", 180, 0, 180);
@@ -73,13 +75,13 @@ void GeometricalEfficiency(const char * fname = "myResult")
//if (i%1000 == 0) cout << "Entry " << i << endl;
tree->GetEntry(i);
// Fill histos
- hEmittTheta->Fill(initCond->GetICEmittedAngleThetaLabWorldFrame(0));
+ hEmittTheta->Fill(initCond->GetThetaLab_WorldFrame(0));
+
if (interCoord->GetDetectedMultiplicity() > 0)
hDetecTheta->Fill(interCoord->GetDetectedAngleTheta(0));
}
- TCanvas *c1 = new TCanvas("c1", "c1");
- c1->Draw();
+ TCanvas* c1 = new TCanvas("c1", "c1");
// Compute relative efficiency in %
TH1F *efficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
// efficiency->SetTitle("Efficiency GASPARD (Spheric version);#Theta [deg];#epsilon [%]");
@@ -89,9 +91,21 @@ void GeometricalEfficiency(const char * fname = "myResult")
efficiency->Draw();
- TCanvas *c2 = new TCanvas("c2", "c2");
- c2->Draw();
+ TCanvas* c2 = new TCanvas("c2", "c2");
hEmittTheta->Draw();
hDetecTheta->SetLineColor(kRed);
hDetecTheta->Draw("same");
+
+ TCanvas* c3 = new TCanvas("c3", "c3");
+ TH1F* SolidA = new TH1F(*hDetecTheta);
+ SolidA->Sumw2();
+ TF1* C = new TF1("C",Form("%i /(4*%f)",nentries,M_PI),0,180);
+ SolidA->Divide(C,1);
+ SolidA->Rebin(2);
+ SolidA->Draw();
+ TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *2*%f/180.",M_PI,M_PI,M_PI),0,180);
+ f->Draw("SAME");
+ SolidA->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ SolidA->GetYaxis()->SetTitle("d#Omega (sr) ");
+ c3->Update();
}
diff --git a/NPLib/InitialConditions/TInitialConditions.cxx b/NPLib/InitialConditions/TInitialConditions.cxx
index 926ca611a..8baf85729 100644
--- a/NPLib/InitialConditions/TInitialConditions.cxx
+++ b/NPLib/InitialConditions/TInitialConditions.cxx
@@ -98,3 +98,31 @@ void TInitialConditions::Dump() const{
}
+
+
+TVector3 TInitialConditions::GetBeamDirection() const{
+ return TVector3( sin(fIC_Incident_Emittance_ThetaX),
+ sin(fIC_Incident_Emittance_PhiY),
+ cos(fIC_Incident_Emittance_ThetaX) + cos(fIC_Incident_Emittance_PhiY));
+}
+
+TVector3 TInitialConditions::GetParticleDirection (const int &i) const {
+ return TVector3( fIC_Momentum_Direction_X[i],
+ fIC_Momentum_Direction_Y[i],
+ fIC_Momentum_Direction_Z[i]);
+}
+
+
+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 3e2efd4d9..2e4eb9cd0 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -30,11 +30,19 @@
// STL Header
#include <vector>
#include <string>
-
+#include <cmath>
using namespace std ;
// Root Header
#include "TObject.h"
+#include "TVector3.h"
+
+// NPTOOL headers
+#include "../include/NPGlobalSystemOfUnits.h"
+#include "../include/NPPhysicalConstants.h"
+#ifdef NP_SYSTEM_OF_UNITS_H
+using namespace NPUNITS;
+#endif
class TInitialConditions : public TObject{
private:
@@ -69,29 +77,29 @@ public:
///////////////////// GetTERS ////////////////////////
// Incident beam parameter
- void SetIncidentParticleName (string Incident_Particle_Name) {fIC_Incident_Particle_Name = Incident_Particle_Name;}
- void SetIncidentInitialKineticEnergy (double Incident_Initial_Kinetic_Energy)
+ void SetIncidentParticleName (const string &Incident_Particle_Name) {fIC_Incident_Particle_Name = Incident_Particle_Name;}
+ void SetIncidentInitialKineticEnergy (const double &Incident_Initial_Kinetic_Energy)
{fIC_Incident_Initial_Kinetic_Energy = Incident_Initial_Kinetic_Energy;}
- void SetIncidentFinalKineticEnergy (double Incident_Final_Kinetic_Energy)
+ void SetIncidentFinalKineticEnergy (const double &Incident_Final_Kinetic_Energy)
{fIC_Incident_Final_Kinetic_Energy = Incident_Final_Kinetic_Energy;}
- void SetIncidentEmittanceTheta (double Incident_Emittance_Theta) {fIC_Incident_Emittance_Theta = Incident_Emittance_Theta;}
- void SetIncidentEmittancePhi (double Incident_Emittance_Phi) {fIC_Incident_Emittance_Phi = Incident_Emittance_Phi;}
- void SetIncidentEmittanceThetaX (double Incident_Emittance_ThetaX) {fIC_Incident_Emittance_ThetaX = Incident_Emittance_ThetaX;}
- void SetIncidentEmittancePhiY (double Incident_Emittance_PhiY) {fIC_Incident_Emittance_PhiY = Incident_Emittance_PhiY;}
+ void SetIncidentEmittanceTheta (const double &Incident_Emittance_Theta) {fIC_Incident_Emittance_Theta = Incident_Emittance_Theta;}
+ void SetIncidentEmittancePhi (const double &Incident_Emittance_Phi) {fIC_Incident_Emittance_Phi = Incident_Emittance_Phi;}
+ void SetIncidentEmittanceThetaX (const double &Incident_Emittance_ThetaX) {fIC_Incident_Emittance_ThetaX = Incident_Emittance_ThetaX;}
+ void SetIncidentEmittancePhiY (const double &Incident_Emittance_PhiY) {fIC_Incident_Emittance_PhiY = Incident_Emittance_PhiY;}
// Beam status at the initial interaction point
- void SetIncidentPositionX (double Incident_Position_X) {fIC_Incident_Position_X = Incident_Position_X;}
- void SetIncidentPositionY (double Incident_Position_Y) {fIC_Incident_Position_Y = Incident_Position_Y;}
- void SetIncidentPositionZ (double Incident_Position_Z) {fIC_Incident_Position_Z = Incident_Position_Z;}
+ void SetIncidentPositionX (const double &Incident_Position_X) {fIC_Incident_Position_X = Incident_Position_X;}
+ void SetIncidentPositionY (const double &Incident_Position_Y) {fIC_Incident_Position_Y = Incident_Position_Y;}
+ void SetIncidentPositionZ (const double &Incident_Position_Z) {fIC_Incident_Position_Z = Incident_Position_Z;}
// emmitted particle
- void SetParticleName (string Particle_Name) {fIC_Particle_Name.push_back(Particle_Name);}
- void SetThetaCM (double ThetaCM) {fIC_ThetaCM.push_back(ThetaCM);}
- void SetKineticEnergy (double Kinetic_Energy) {fIC_Kinetic_Energy.push_back(Kinetic_Energy);}
- void SetMomentumDirectionX (double Momentum_Direction_X) {fIC_Momentum_Direction_X.push_back(Momentum_Direction_X);}
- void SetMomentumDirectionY (double Momentum_Direction_Y) {fIC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}
- void SetMomentumDirectionZ (double Momentum_Direction_Z) {fIC_Momentum_Direction_Z.push_back(Momentum_Direction_Z);}
+ void SetParticleName (const string &Particle_Name) {fIC_Particle_Name.push_back(Particle_Name);}
+ void SetThetaCM (const double &ThetaCM) {fIC_ThetaCM.push_back(ThetaCM);}
+ void SetKineticEnergy (const double &Kinetic_Energy) {fIC_Kinetic_Energy.push_back(Kinetic_Energy);}
+ void SetMomentumDirectionX (const double &Momentum_Direction_X) {fIC_Momentum_Direction_X.push_back(Momentum_Direction_X);}
+ void SetMomentumDirectionY (const double &Momentum_Direction_Y) {fIC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}
+ void SetMomentumDirectionZ (const double &Momentum_Direction_Z) {fIC_Momentum_Direction_Z.push_back(Momentum_Direction_Z);}
///////////////////// GETTERS ////////////////////////
// Incident beam parameter
@@ -107,12 +115,17 @@ public:
double GetIncidentPositionZ () const {return fIC_Incident_Position_Z ;}
// emmitted particle
- string GetParticleName (int i) const {return fIC_Particle_Name[i];}
- double GetThetaCM (int i) const {return fIC_ThetaCM[i];}
- double GetKineticEnergy (int i) const {return fIC_Kinetic_Energy[i];}
- double GetMomentumDirectionX (int i) const {return fIC_Momentum_Direction_X[i];}
- double GetMomentumDirectionY (int i) const {return fIC_Momentum_Direction_Y[i];}
- double GetMomentumDirectionZ (int i) const {return fIC_Momentum_Direction_Z[i];}
+ string GetParticleName (const int &i) const {return fIC_Particle_Name[i];}
+ double GetThetaCM (const int &i) const {return fIC_ThetaCM[i];}
+ double GetKineticEnergy (const int &i) const {return fIC_Kinetic_Energy[i];}
+ double GetMomentumDirectionX (const int &i) const {return fIC_Momentum_Direction_X[i];}
+ double GetMomentumDirectionY (const int &i) const {return fIC_Momentum_Direction_Y[i];}
+ double GetMomentumDirectionZ (const int &i) const {return fIC_Momentum_Direction_Z[i];}
+
+ TVector3 GetBeamDirection () const ;
+ TVector3 GetParticleDirection (const int &i) const ;
+ double GetThetaLab_WorldFrame (const int &i) const ;
+ double GetThetaLab_IncidentFrame (const int &i) const ;
unsigned int GetEmittedMult() const {return fIC_Particle_Name.size();}
diff --git a/NPLib/InteractionCoordinates/TInteractionCoordinates.h b/NPLib/InteractionCoordinates/TInteractionCoordinates.h
index 051fa8c17..ec2d303e3 100644
--- a/NPLib/InteractionCoordinates/TInteractionCoordinates.h
+++ b/NPLib/InteractionCoordinates/TInteractionCoordinates.h
@@ -70,12 +70,12 @@ public:
Int_t GetDetectedMultiplicity() {return fDetected_Position_X.size();}
// Incident particle properties (before interactions in the target)
// Vertex of interaction
- Double_t GetDetectedPositionX(Int_t i) {return fDetected_Position_X.at(i);}
- Double_t GetDetectedPositionY(Int_t i) {return fDetected_Position_Y.at(i);}
- Double_t GetDetectedPositionZ(Int_t i) {return fDetected_Position_Z.at(i);}
+ Double_t GetDetectedPositionX(Int_t i) {return fDetected_Position_X[i];}
+ Double_t GetDetectedPositionY(Int_t i) {return fDetected_Position_Y[i];}
+ Double_t GetDetectedPositionZ(Int_t i) {return fDetected_Position_Z[i];}
// Incident particle angles
- Double_t GetDetectedAngleTheta(Int_t i) {return fDetected_Angle_Theta.at(i);}
- Double_t GetDetectedAnglePhi(Int_t i) {return fDetected_Angle_Phi.at(i);}
+ Double_t GetDetectedAngleTheta(Int_t i) {return fDetected_Angle_Theta[i];}
+ Double_t GetDetectedAnglePhi(Int_t i) {return fDetected_Angle_Phi[i];}
ClassDef(TInteractionCoordinates, 1) // InteractionCoordinates structure
};
diff --git a/NPSimulation/include/Sharc.hh b/NPSimulation/include/Sharc.hh
index 5dbf7f4dc..8a8fafe3b 100644
--- a/NPSimulation/include/Sharc.hh
+++ b/NPSimulation/include/Sharc.hh
@@ -194,7 +194,6 @@ private:
G4MultiFunctionalDetector* m_BOXScorer ;
G4MultiFunctionalDetector* m_PADScorer ;
G4MultiFunctionalDetector* m_QQQScorer ;
-
private:
// Initialize material used in detector definition
diff --git a/NPSimulation/src/ParticleStack.cc b/NPSimulation/src/ParticleStack.cc
index 5c4eb7617..61ec26be6 100644
--- a/NPSimulation/src/ParticleStack.cc
+++ b/NPSimulation/src/ParticleStack.cc
@@ -211,6 +211,7 @@ string ParticleStack::ChangeNameToG4Standard(string OriginalName){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ParticleStack::ShootAllParticle(G4Event* anEvent){
unsigned int size = m_ParticleStack.size();
+ m_InitialConditions->Clear();
for(unsigned int i = 0 ; i < size ; i++){
diff --git a/NPSimulation/src/Sharc.cc b/NPSimulation/src/Sharc.cc
index 3a6c6e9bd..c1d34f433 100644
--- a/NPSimulation/src/Sharc.cc
+++ b/NPSimulation/src/Sharc.cc
@@ -661,6 +661,15 @@ void Sharc::ReadSensitive(const G4Event* event){
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
@@ -687,6 +696,7 @@ void Sharc::ReadSensitive(const G4Event* event){
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
@@ -722,6 +732,14 @@ void Sharc::ReadSensitive(const G4Event* event){
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
@@ -762,6 +780,8 @@ void Sharc::InitializeScorers(){
QQQ_Wafer_NumberOf_AnnularStrip,
EnergyThreshold);
+
+
//and register it to the multifunctionnal detector
m_BOXScorer->RegisterPrimitive(BOXScorer);
m_PADScorer->RegisterPrimitive(PADScorer);
diff --git a/NPSimulation/src/SharcScorers.cc b/NPSimulation/src/SharcScorers.cc
index 51c08da2f..8cda2113c 100644
--- a/NPSimulation/src/SharcScorers.cc
+++ b/NPSimulation/src/SharcScorers.cc
@@ -48,7 +48,7 @@ PS_Silicon_Rectangle::~PS_Silicon_Rectangle(){
G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
// contain Energy Time, DetNbr, StripFront and StripBack
- G4double* EnergyAndTime = new G4double[5];
+ G4double* EnergyAndTime = new G4double[10];
EnergyAndTime[0] = aStep->GetTotalEnergyDeposit();
// If the energy deposit is below the threshold, the deposit is ignored
@@ -61,6 +61,14 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
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_StripLengthNumber = (int)((m_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength ) + 1 ;
@@ -70,6 +78,7 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
EnergyAndTime[3] = m_StripLengthNumber;
EnergyAndTime[4] = m_StripWidthNumber;
+
//Rare case where particle is close to edge of silicon plan
if (m_StripLengthNumber == m_NumberOfStripLength+1) m_StripLengthNumber = m_StripLengthNumber;
if (m_StripWidthNumber == m_NumberOfStripWidth+1) m_StripWidthNumber = m_StripWidthNumber;
@@ -94,6 +103,11 @@ 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();
}
@@ -139,7 +153,7 @@ PS_Silicon_Annular::~PS_Silicon_Annular(){
G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){
// contain Energy Time, DetNbr, StripFront and StripBack
- G4double* EnergyAndTime = new G4double[5];
+ G4double* EnergyAndTime = new G4double[10];
EnergyAndTime[0] = aStep->GetTotalEnergyDeposit();
// If the energy deposit is below the threshold, the deposit is ignored
@@ -152,6 +166,14 @@ G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){
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 ;
@@ -185,6 +207,11 @@ 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();
}
diff --git a/NPSimulation/vis.mac b/NPSimulation/vis.mac
index d0df3eacf..96d2fc4a6 100644
--- a/NPSimulation/vis.mac
+++ b/NPSimulation/vis.mac
@@ -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