From 2e4e7115d659f3f2cfaf47cc9ee695339a75bbee Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Mon, 15 Dec 2014 15:58:53 +0000
Subject: [PATCH] * Updating Paris: - Now follow the new logique - Less
file - Updated geometry
---
NPLib/Paris/TParisData.cxx | 82 +-
NPLib/Paris/TParisData.h | 349 +++--
NPLib/Paris/TParisPhysics.cxx | 55 +-
NPLib/Paris/TParisPhysics.h | 22 +-
NPSimulation/Paris/Paris.cc | 649 +++++++--
NPSimulation/Paris/Paris.hh | 184 ++-
NPSimulation/Paris/ParisCluster.cc | 1195 -----------------
NPSimulation/Paris/ParisCluster.hh | 169 ---
NPSimulation/Paris/ParisModule.cc | 64 -
NPSimulation/Paris/ParisModule.hh | 97 --
NPSimulation/Paris/ParisPhoswich.cc | 1068 ---------------
NPSimulation/Paris/ParisPhoswich.hh | 173 ---
NPSimulation/Paris/ParisScorers.cc | 516 -------
NPSimulation/Paris/ParisScorers.hh | 199 ---
...eStripScorers.hh => CalorimeterScorers.hh} | 76 +-
NPSimulation/include/SiliconScorers.hh | 48 +-
...eStripScorers.cc => CalorimeterScorers.cc} | 96 +-
NPSimulation/src/SiliconScorers.cc | 109 +-
18 files changed, 1127 insertions(+), 4024 deletions(-)
delete mode 100644 NPSimulation/Paris/ParisCluster.cc
delete mode 100644 NPSimulation/Paris/ParisCluster.hh
delete mode 100644 NPSimulation/Paris/ParisModule.cc
delete mode 100644 NPSimulation/Paris/ParisModule.hh
delete mode 100644 NPSimulation/Paris/ParisPhoswich.cc
delete mode 100644 NPSimulation/Paris/ParisPhoswich.hh
delete mode 100644 NPSimulation/Paris/ParisScorers.cc
delete mode 100644 NPSimulation/Paris/ParisScorers.hh
rename NPSimulation/include/{ResistiveStripScorers.hh => CalorimeterScorers.hh} (60%)
rename NPSimulation/src/{ResistiveStripScorers.cc => CalorimeterScorers.cc} (50%)
diff --git a/NPLib/Paris/TParisData.cxx b/NPLib/Paris/TParisData.cxx
index 116b693ce..3a18f0220 100644
--- a/NPLib/Paris/TParisData.cxx
+++ b/NPLib/Paris/TParisData.cxx
@@ -23,56 +23,36 @@
#include "TParisData.h"
ClassImp(TParisData)
-TParisData::TParisData()
-{
- // Default constructor
- Clear();
+////////////////////////////////////////////////////////////////////////////////
+TParisData::TParisData(){
}
-
-
-TParisData::~TParisData()
-{
+TParisData::~TParisData(){
}
-
-
-void TParisData::Clear()
-{
-
- /*
- fParis_Energy.clear();
- fParis_Number.clear();
- fParis_Time.clear();
- */
-
-
- fPARIS_LaBr3Stage_E_DetectorNbr.clear();
- fPARIS_LaBr3Stage_E_CrystalNbr.clear();
+////////////////////////////////////////////////////////////////////////////////
+void TParisData::Clear(){
+ fPARIS_LaBr3Stage_E_ClusterNbr.clear();
+ fPARIS_LaBr3Stage_E_PhoswichNbr.clear();
fPARIS_LaBr3Stage_E_Energy.clear();
- fPARIS_LaBr3Stage_Eff_phpeak.clear();
// Time
- fPARIS_LaBr3Stage_T_DetectorNbr.clear();
- fPARIS_LaBr3Stage_T_CrystalNbr.clear();
+ fPARIS_LaBr3Stage_T_ClusterNbr.clear();
+ fPARIS_LaBr3Stage_T_PhoswichNbr.clear();
fPARIS_LaBr3Stage_T_Time.clear();
-
- // Second Stage CsI
- // CsI
+ // Second Stage NaI
+ // NaI
// Energy
- fPARIS_CsIStage_E_DetectorNbr.clear();
- fPARIS_CsIStage_E_CrystalNbr.clear();
- fPARIS_CsIStage_E_Energy.clear();
- fPARIS_CsIStage_Eff_phpeak.clear();
+ fPARIS_NaIStage_E_ClusterNbr.clear();
+ fPARIS_NaIStage_E_PhoswichNbr.clear();
+ fPARIS_NaIStage_E_Energy.clear();
// Time
- fPARIS_CsIStage_T_DetectorNbr.clear();
- fPARIS_CsIStage_T_CrystalNbr.clear();
- fPARIS_CsIStage_T_Time.clear();
-
+ fPARIS_NaIStage_T_ClusterNbr.clear();
+ fPARIS_NaIStage_T_PhoswichNbr.clear();
+ fPARIS_NaIStage_T_Time.clear();
}
-
-
+////////////////////////////////////////////////////////////////////////////////
void TParisData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
@@ -82,23 +62,21 @@ void TParisData::Dump() const
}
*/
- for (UShort_t i = 0; i < fPARIS_LaBr3Stage_E_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fPARIS_LaBr3Stage_E_DetectorNbr[i] << " Crystal: " << fPARIS_LaBr3Stage_E_CrystalNbr[i] << " Energy: " << fPARIS_LaBr3Stage_E_Energy[i] << endl;
+ for (UShort_t i = 0; i < fPARIS_LaBr3Stage_E_ClusterNbr.size(); i++)
+ cout << "DetNbr: " << fPARIS_LaBr3Stage_E_ClusterNbr[i] << " Phoswich: " << fPARIS_LaBr3Stage_E_PhoswichNbr[i] << " Energy: " << fPARIS_LaBr3Stage_E_Energy[i] << endl;
// (X,T)
- cout << "PARIS_LaBr3Stage_T_Mult = " << fPARIS_LaBr3Stage_T_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fPARIS_LaBr3Stage_T_DetectorNbr.size(); i++)
- cout << "DetNbr: " << fPARIS_LaBr3Stage_T_DetectorNbr[i] << " Crystal: " << fPARIS_LaBr3Stage_T_CrystalNbr[i] << " Time: " << fPARIS_LaBr3Stage_T_Time[i] << endl;
+ cout << "PARIS_LaBr3Stage_T_Mult = " << fPARIS_LaBr3Stage_T_ClusterNbr.size() << endl;
+ for (UShort_t i = 0; i < fPARIS_LaBr3Stage_T_ClusterNbr.size(); i++)
+ cout << "DetNbr: " << fPARIS_LaBr3Stage_T_ClusterNbr[i] << " Phoswich: " << fPARIS_LaBr3Stage_T_PhoswichNbr[i] << " Time: " << fPARIS_LaBr3Stage_T_Time[i] << endl;
// Second Stage
// Energy
- cout << "PARIS_CsIStage_E_Mult = " << fPARIS_CsIStage_E_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fPARIS_CsIStage_E_DetectorNbr.size(); i++)
- cout << "Det: " << fPARIS_CsIStage_E_DetectorNbr[i] << " Crystal: " << fPARIS_CsIStage_E_CrystalNbr[i] << " Energy: " << fPARIS_CsIStage_E_Energy[i] << endl;
+ cout << "PARIS_NaIStage_E_Mult = " << fPARIS_NaIStage_E_ClusterNbr.size() << endl;
+ for (UShort_t i = 0; i < fPARIS_NaIStage_E_ClusterNbr.size(); i++)
+ cout << "Det: " << fPARIS_NaIStage_E_ClusterNbr[i] << " Phoswich: " << fPARIS_NaIStage_E_PhoswichNbr[i] << " Energy: " << fPARIS_NaIStage_E_Energy[i] << endl;
// Time
- cout << "PARIS_CsIStage_T_Mult = " << fPARIS_CsIStage_T_DetectorNbr.size() << endl;
- for (UShort_t i = 0; i < fPARIS_CsIStage_T_DetectorNbr.size(); i++)
- cout << "Det: " << fPARIS_CsIStage_T_DetectorNbr[i] << " Crystal: " << fPARIS_CsIStage_T_CrystalNbr[i] << " Time: " << fPARIS_CsIStage_T_Time[i] << endl;
-
-
-
+ cout << "PARIS_NaIStage_T_Mult = " << fPARIS_NaIStage_T_ClusterNbr.size() << endl;
+ for (UShort_t i = 0; i < fPARIS_NaIStage_T_ClusterNbr.size(); i++)
+ cout << "Det: " << fPARIS_NaIStage_T_ClusterNbr[i] << " Phoswich: " << fPARIS_NaIStage_T_PhoswichNbr[i] << " Time: " << fPARIS_NaIStage_T_Time[i] << endl;
}
+
diff --git a/NPLib/Paris/TParisData.h b/NPLib/Paris/TParisData.h
index 28aa4f4cd..0e893d90b 100644
--- a/NPLib/Paris/TParisData.h
+++ b/NPLib/Paris/TParisData.h
@@ -29,200 +29,169 @@ using namespace std ;
class TParisData : public TObject {
- protected:
- // First Stage LaBr
- // LaBr3
- // Energy
- vector<UShort_t> fPARIS_LaBr3Stage_E_DetectorNbr;
- vector<UShort_t> fPARIS_LaBr3Stage_E_CrystalNbr;
- vector<Double_t> fPARIS_LaBr3Stage_E_Energy;
- vector<Double_t> fPARIS_LaBr3Stage_Eff_phpeak;
- // Time
- vector<UShort_t> fPARIS_LaBr3Stage_T_DetectorNbr;
- vector<UShort_t> fPARIS_LaBr3Stage_T_CrystalNbr;
- vector<Double_t> fPARIS_LaBr3Stage_T_Time;
-
-
- // Second Stage CsI
- // CsI
- // Energy
- vector<UShort_t> fPARIS_CsIStage_E_DetectorNbr;
- vector<UShort_t> fPARIS_CsIStage_E_CrystalNbr;
- vector<Double_t> fPARIS_CsIStage_E_Energy;
- vector<Double_t> fPARIS_CsIStage_Eff_phpeak;
- // Time
- vector<UShort_t> fPARIS_CsIStage_T_DetectorNbr;
- vector<UShort_t> fPARIS_CsIStage_T_CrystalNbr;
- vector<Double_t> fPARIS_CsIStage_T_Time;
-
- /*
- private:
- vector<double> fParis_Energy;
- vector<double> fParis_Time;
- vector<short> fParis_Number;
- */
-
- public:
- TParisData();
- virtual ~TParisData();
-
- void Clear();
- void Clear(const Option_t*) {};
- void Dump() const;
-
- ///////////////////// GETTERS ////////////////////////
-
- //
- // First stage
- //
- // (E)
- UShort_t GetPARISLaBr3StageEMult() {
- return fPARIS_LaBr3Stage_E_DetectorNbr.size(); // TODO: Maybe change to CrystalNbr for multiplicity
- }
- UShort_t GetPARISLaBr3StageEDetectorNbr(Int_t i) {
- return fPARIS_LaBr3Stage_E_DetectorNbr.at(i);
- }
- UShort_t GetPARISLaBr3StageECrystalNbr(Int_t i) {
- return fPARIS_LaBr3Stage_E_CrystalNbr.at(i);
- }
- Double_t GetPARISLaBr3StageEEnergy(Int_t i) {
+ protected:
+ // First Stage LaBr3
+ // Energy
+ vector<UShort_t> fPARIS_LaBr3Stage_E_ClusterNbr;
+ vector<UShort_t> fPARIS_LaBr3Stage_E_PhoswichNbr;
+ vector<Double_t> fPARIS_LaBr3Stage_E_Energy;
+ // Time
+ vector<UShort_t> fPARIS_LaBr3Stage_T_ClusterNbr;
+ vector<UShort_t> fPARIS_LaBr3Stage_T_PhoswichNbr;
+ vector<Double_t> fPARIS_LaBr3Stage_T_Time;
+
+ // Second Stage NaI
+ // Energy
+ vector<UShort_t> fPARIS_NaIStage_E_ClusterNbr;
+ vector<UShort_t> fPARIS_NaIStage_E_PhoswichNbr;
+ vector<Double_t> fPARIS_NaIStage_E_Energy;
+ // Time
+ vector<UShort_t> fPARIS_NaIStage_T_ClusterNbr;
+ vector<UShort_t> fPARIS_NaIStage_T_PhoswichNbr;
+ vector<Double_t> fPARIS_NaIStage_T_Time;
+
+ public:
+ TParisData();
+ virtual ~TParisData();
+
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+ ///////////////////// GETTERS ////////////////////////
+ // First stage
+ // (E)
+ UShort_t GetPARISLaBr3StageEMult() {
+ return fPARIS_LaBr3Stage_E_ClusterNbr.size();
+ }
+ UShort_t GetPARISLaBr3StageEClusterNbr(Int_t i) {
+ return fPARIS_LaBr3Stage_E_ClusterNbr.at(i);
+ }
+ UShort_t GetPARISLaBr3StageEPhoswichNbr(Int_t i) {
+ return fPARIS_LaBr3Stage_E_PhoswichNbr.at(i);
+ }
+ Double_t GetPARISLaBr3StageEEnergy(Int_t i) {
return fPARIS_LaBr3Stage_E_Energy.at(i);
- }
- Double_t GetPARISLaBr3StageEffphpeak(Int_t i) {
- return fPARIS_LaBr3Stage_Eff_phpeak.at(i);
- }
-
-
- // (T)
- UShort_t GetPARISLaBr3StageTMult() {
- return fPARIS_LaBr3Stage_E_DetectorNbr.size();
- }
- UShort_t GetPARISLaBr3StageTDetectorNbr(Int_t i) {
- return fPARIS_LaBr3Stage_T_DetectorNbr.at(i);
- }
- UShort_t GetPARISLaBr3StageTCrystalNbr(Int_t i) {
- return fPARIS_LaBr3Stage_T_CrystalNbr.at(i);
- }
- Double_t GetPARISLaBr3StageTTime(Int_t i) {
+ }
+
+
+ // (T)
+ UShort_t GetPARISLaBr3StageTMult() {
+ return fPARIS_LaBr3Stage_E_ClusterNbr.size();
+ }
+ UShort_t GetPARISLaBr3StageTClusterNbr(Int_t i) {
+ return fPARIS_LaBr3Stage_T_ClusterNbr.at(i);
+ }
+ UShort_t GetPARISLaBr3StageTPhoswichNbr(Int_t i) {
+ return fPARIS_LaBr3Stage_T_PhoswichNbr.at(i);
+ }
+ Double_t GetPARISLaBr3StageTTime(Int_t i) {
return fPARIS_LaBr3Stage_T_Time.at(i);
- }
-
- //
- // Second stage (CsI
- //
- // (E)
- UShort_t GetPARISCsIStageEMult() {
- return fPARIS_CsIStage_E_DetectorNbr.size();
- }
- UShort_t GetPARISCsIStageEDetectorNbr(Int_t i) {
- return fPARIS_CsIStage_E_DetectorNbr.at(i);
- }
- UShort_t GetPARISCsIStageECrystalNbr(Int_t i) {
- return fPARIS_CsIStage_E_CrystalNbr.at(i);
- }
- Double_t GetPARISCsIStageEEnergy(Int_t i) {
- return fPARIS_CsIStage_E_Energy.at(i);
- }
- Double_t GetPARISCsIStageEffphpeak(Int_t i) {
- return fPARIS_CsIStage_Eff_phpeak.at(i);
- }
-
- // (T)
- UShort_t GetPARISCsIStageTMult() {
- return fPARIS_CsIStage_E_DetectorNbr.size();
- }
- UShort_t GetPARISCsIStageTDetectorNbr(Int_t i) {
- return fPARIS_CsIStage_T_DetectorNbr.at(i);
- }
- UShort_t GetPARISCsIStageTCrystalNbr(Int_t i) {
- return fPARIS_CsIStage_T_CrystalNbr.at(i);
- }
- Double_t GetPARISCsIStageTTime(Int_t i) {
- return fPARIS_CsIStage_T_Time.at(i);
- }
-
- /*
- // (E)
- //double GetEnergy(int i) {return fParis_Energy[i];}
- // (T)
- //double GetTime(int i) {return fParis_Time[i];}
- // (N)
- int GetParisNumber(int i) {return fParis_Number[i];}
- double GetEnergySize() {return fParis_Energy.size();}
- // (T)
- double GetTimeSize() {return fParis_Time.size();}
- // (N)
- int GetParisNumberSize() {return fParis_Number.size();}
- */
-
- ///////////////////// SETTERS ////////////////////////
-
- //
- // First stage
- //
- // (E)
-
- void SetPARISLaBr3StageEDetectorNbr(UShort_t DetNbr) {
- fPARIS_LaBr3Stage_E_DetectorNbr.push_back(DetNbr);
- }
- void SetPARISLaBr3StageECrystalNbr(UShort_t CrystalNbr) {
- fPARIS_LaBr3Stage_E_CrystalNbr.push_back(CrystalNbr);
- }
- void SetPARISLaBr3StageEEnergy(Double_t Energy) {
+ }
+
+ // Second stage NaI
+ // (E)
+ UShort_t GetPARISNaIStageEMult() {
+ return fPARIS_NaIStage_E_ClusterNbr.size();
+ }
+ UShort_t GetPARISNaIStageEClusterNbr(Int_t i) {
+ return fPARIS_NaIStage_E_ClusterNbr.at(i);
+ }
+ UShort_t GetPARISNaIStageEPhoswichNbr(Int_t i) {
+ return fPARIS_NaIStage_E_PhoswichNbr.at(i);
+ }
+ Double_t GetPARISNaIStageEEnergy(Int_t i) {
+ return fPARIS_NaIStage_E_Energy.at(i);
+ }
+
+ // (T)
+ UShort_t GetPARISNaIStageTMult() {
+ return fPARIS_NaIStage_E_ClusterNbr.size();
+ }
+ UShort_t GetPARISNaIStageTClusterNbr(Int_t i) {
+ return fPARIS_NaIStage_T_ClusterNbr.at(i);
+ }
+ UShort_t GetPARISNaIStageTPhoswichNbr(Int_t i) {
+ return fPARIS_NaIStage_T_PhoswichNbr.at(i);
+ }
+ Double_t GetPARISNaIStageTTime(Int_t i) {
+ return fPARIS_NaIStage_T_Time.at(i);
+ }
+
+ ///////////////////// SETTERS ////////////////////////
+ // First stage
+ // (E)
+ void SetParisLaBr3E(UShort_t DetNbr,UShort_t PhoswichNbr, Double_t E){
+ fPARIS_LaBr3Stage_E_ClusterNbr.push_back(DetNbr);
+ fPARIS_LaBr3Stage_E_PhoswichNbr.push_back(PhoswichNbr);
+ fPARIS_LaBr3Stage_E_Energy.push_back(E);
+ }
+
+ void SetPARISLaBr3StageEClusterNbr(UShort_t DetNbr) {
+ fPARIS_LaBr3Stage_E_ClusterNbr.push_back(DetNbr);
+ }
+ void SetPARISLaBr3StageEPhoswichNbr(UShort_t PhoswichNbr) {
+ fPARIS_LaBr3Stage_E_PhoswichNbr.push_back(PhoswichNbr);
+ }
+ void SetPARISLaBr3StageEEnergy(Double_t Energy) {
fPARIS_LaBr3Stage_E_Energy.push_back(Energy);
- }
- void SetPARISLaBr3StageEffphpeak(Double_t Energy) {
- fPARIS_LaBr3Stage_Eff_phpeak.push_back(Energy);
- }
-
-
- // (T)
- void SetPARISLaBr3StageTDetectorNbr(UShort_t DetNbr) {
- fPARIS_LaBr3Stage_T_DetectorNbr.push_back(DetNbr);
- }
- void SetPARISLaBr3StageTCrystalNbr(UShort_t CrystalNbr) {
- fPARIS_LaBr3Stage_T_CrystalNbr.push_back(CrystalNbr);
- }
- void SetPARISLaBr3StageTTime(Double_t Time) {
+ }
+
+ // (T)
+ void SetParisLaBr3T(UShort_t DetNbr,UShort_t PhoswichNbr, Double_t T){
+ fPARIS_LaBr3Stage_T_ClusterNbr.push_back(DetNbr);
+ fPARIS_LaBr3Stage_T_PhoswichNbr.push_back(PhoswichNbr);
+ fPARIS_LaBr3Stage_T_Time.push_back(T);
+ }
+
+
+
+ void SetPARISLaBr3StageTClusterNbr(UShort_t DetNbr) {
+ fPARIS_LaBr3Stage_T_ClusterNbr.push_back(DetNbr);
+ }
+ void SetPARISLaBr3StageTPhoswichNbr(UShort_t PhoswichNbr) {
+ fPARIS_LaBr3Stage_T_PhoswichNbr.push_back(PhoswichNbr);
+ }
+ void SetPARISLaBr3StageTTime(Double_t Time) {
fPARIS_LaBr3Stage_T_Time.push_back(Time);
- }
-
- //
- // Second stage (CsI
- //
- // (E)
- void SetPARISCsIStageEDetectorNbr(UShort_t DetNbr) {
- fPARIS_CsIStage_E_DetectorNbr.push_back(DetNbr);
- }
- void SetPARISCsIStageECrystalNbr(UShort_t CrystalNbr) {
- fPARIS_CsIStage_E_CrystalNbr.push_back(CrystalNbr);
- }
- void SetPARISCsIStageEEnergy(Double_t Energy) {
- fPARIS_CsIStage_E_Energy.push_back(Energy);
- }
- void SetPARISCsIStageEffphpeak(Double_t Energy) {
- fPARIS_CsIStage_Eff_phpeak.push_back(Energy);
- }
-
- // (T)
- void SetPARISCsIStageTDetectorNbr(UShort_t DetNbr) {
- fPARIS_CsIStage_T_DetectorNbr.push_back(DetNbr);
- }
- void SetPARISCsIStageTCrystalNbr(UShort_t CrystalNbr) {
- fPARIS_CsIStage_T_CrystalNbr.push_back(CrystalNbr);
- }
- void SetPARISCsIStageTTime(Double_t Time) {
- fPARIS_CsIStage_T_Time.push_back(Time);
- }
-
-
- /*
- // (E)
- void SetEnergy(double E) {fParis_Energy.push_back(E);}
- void SetTime(double T) {fParis_Time.push_back(T);}
- void SetParisNumber(int N) {fParis_Number.push_back(N);}
- */
- ClassDef(TParisData,1) // ParisData structure
+ }
+
+ // Second stage NaI
+ // (E)
+ void SetParisNaIE(UShort_t DetNbr,UShort_t PhoswichNbr, Double_t E){
+ fPARIS_NaIStage_E_ClusterNbr.push_back(DetNbr);
+ fPARIS_NaIStage_E_PhoswichNbr.push_back(PhoswichNbr);
+ fPARIS_NaIStage_E_Energy.push_back(E);
+ }
+
+ void SetPARISNaIStageEClusterNbr(UShort_t DetNbr) {
+ fPARIS_NaIStage_E_ClusterNbr.push_back(DetNbr);
+ }
+ void SetPARISNaIStageEPhoswichNbr(UShort_t PhoswichNbr) {
+ fPARIS_NaIStage_E_PhoswichNbr.push_back(PhoswichNbr);
+ }
+ void SetPARISNaIStageEEnergy(Double_t Energy) {
+ fPARIS_NaIStage_E_Energy.push_back(Energy);
+ }
+
+ // (T)
+ void SetParisNaIT(UShort_t DetNbr,UShort_t PhoswichNbr, Double_t T){
+ fPARIS_NaIStage_T_ClusterNbr.push_back(DetNbr);
+ fPARIS_NaIStage_T_PhoswichNbr.push_back(PhoswichNbr);
+ fPARIS_NaIStage_T_Time.push_back(T);
+ }
+
+ void SetPARISNaIStageTClusterNbr(UShort_t DetNbr) {
+ fPARIS_NaIStage_T_ClusterNbr.push_back(DetNbr);
+ }
+ void SetPARISNaIStageTPhoswichNbr(UShort_t PhoswichNbr) {
+ fPARIS_NaIStage_T_PhoswichNbr.push_back(PhoswichNbr);
+ }
+ void SetPARISNaIStageTTime(Double_t Time) {
+ fPARIS_NaIStage_T_Time.push_back(Time);
+ }
+
+ ClassDef(TParisData,1) // ParisData structure
};
#endif
diff --git a/NPLib/Paris/TParisPhysics.cxx b/NPLib/Paris/TParisPhysics.cxx
index f7d56f835..724855ab9 100644
--- a/NPLib/Paris/TParisPhysics.cxx
+++ b/NPLib/Paris/TParisPhysics.cxx
@@ -63,12 +63,12 @@ void TParisPhysics::BuildSimplePhysicalEvent(){
void TParisPhysics::BuildPhysicalEvent(){
int multLaBrE = m_EventData->GetPARISLaBr3StageEMult();
- int multCsIE = m_EventData->GetPARISCsIStageEMult();
+ int multNaIE = m_EventData->GetPARISNaIStageEMult();
//cout << "multLaBr= " << multLaBrE << endl;
- //cout << "multCsI= " << multCsIE << endl;
+ //cout << "multNaI= " << multNaIE << endl;
- ParisEventMult=multLaBrE+multCsIE;
+ ParisEventMult=multLaBrE+multNaIE;
// get energy from strips and store it
if(ParisEventMult>=1)
{
@@ -98,18 +98,18 @@ void TParisPhysics::BuildPhysicalEvent(){
ParisInTotalEnergy.push_back(EnergyTot);
}
- if(multCsIE>=1){
+ if(multNaIE>=1){
double EnergySecond;
double EnergyTotSecond;
- for(int j=0;j<multCsIE;j++)
+ for(int j=0;j<multNaIE;j++)
{
- EnergySecond = m_EventData->GetPARISCsIStageEEnergy(j);
- ParisCsI_E.push_back(EnergySecond);
+ EnergySecond = m_EventData->GetPARISNaIStageEEnergy(j);
+ ParisNaI_E.push_back(EnergySecond);
EnergyTotSecond +=EnergySecond;
EnergyTot += EnergySecond;
- //cout << "Energy CsI=" << EnergySecond << endl;
- //cout << "Energytot CsI=" << EnergyTot << endl;
+ //cout << "Energy NaI=" << EnergySecond << endl;
+ //cout << "Energytot NaI=" << EnergyTot << endl;
}
// Fill total energy in outter shell
@@ -138,13 +138,13 @@ void TParisPhysics::Clear(){
//FirstStage_X.clear();
//FirstStage_Y.clear();
- // CsI
- ParisCsI_E.clear();
+ // NaI
+ ParisNaI_E.clear();
//SecondStage_T.clear();
//SecondStage_N.clear();
/*
- // CsI
+ // NaI
ThirdStage_E.clear();
ThirdStage_T.clear();
ThirdStage_N.clear();
@@ -158,11 +158,6 @@ void TParisPhysics::ReadConfiguration(string Path) {
string LineBuffer ;
string DataBuffer ;
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
double Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Az, Bz, Cz, Dz;
TVector3 A, B, C, D;
double Theta = 0, Phi = 0, R = 0, beta_u = 0 , beta_v = 0 , beta_w = 0;
@@ -214,7 +209,7 @@ void TParisPhysics::ReadConfiguration(string Path) {
}
// Position method
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ else if (DataBuffer=="A=") {
check_A = true;
ConfigFile >> DataBuffer ;
Ax = atof(DataBuffer.c_str()) ;
@@ -229,7 +224,7 @@ void TParisPhysics::ReadConfiguration(string Path) {
A = TVector3(Ax, Ay, Az);
cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
}
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
+ else if (DataBuffer=="B=") {
check_B = true;
ConfigFile >> DataBuffer ;
Bx = atof(DataBuffer.c_str()) ;
@@ -244,7 +239,7 @@ void TParisPhysics::ReadConfiguration(string Path) {
B = TVector3(Bx, By, Bz);
cout << "X128 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
}
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
+ else if (DataBuffer=="C=") {
check_C = true;
ConfigFile >> DataBuffer ;
Cx = atof(DataBuffer.c_str()) ;
@@ -259,7 +254,7 @@ void TParisPhysics::ReadConfiguration(string Path) {
C = TVector3(Cx, Cy, Cz);
cout << "X1 Y128 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
}
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
+ else if (DataBuffer=="D=") {
check_D = true;
ConfigFile >> DataBuffer ;
Dx = atof(DataBuffer.c_str()) ;
@@ -276,28 +271,28 @@ void TParisPhysics::ReadConfiguration(string Path) {
} // End Position Method
// Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
+ else if (DataBuffer== "Theta="){
check_Theta = true;
ConfigFile >> DataBuffer ;
Theta = atof(DataBuffer.c_str()) ;
Theta = Theta ;
cout << "Theta: " << Theta << endl;
}
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
+ else if (DataBuffer == "Phi=") {
check_Phi = true;
ConfigFile >> DataBuffer ;
Phi = atof(DataBuffer.c_str()) ;
Phi = Phi ;
cout << "Phi: " << Phi << endl;
}
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ else if (DataBuffer =="R="){
check_R = true;
ConfigFile >> DataBuffer ;
R = atof(DataBuffer.c_str()) ;
R = R ;
cout << "R: " << R << endl;
}
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
+ else if (DataBuffer=="beta=") {
check_beta = true;
ConfigFile >> DataBuffer ;
beta_u = atof(DataBuffer.c_str()) ;
@@ -518,7 +513,7 @@ void TParisPhysics::ReadCalibrationFile(string Path){
/* int Calibration_Si_E_Order;
int Calibration_Si_T_Order;
int Calibration_SiLi_E_Order;
- int Calibration_CsI_E_Order;
+ int Calibration_NaI_E_Order;
*/
// Calibration_Si_X_E[DetectorNumber][StripNumber][Order of Coeff]
vector< vector< vector< double > > > Calibration_Si_X_E ;
@@ -530,13 +525,13 @@ void TParisPhysics::ReadCalibrationFile(string Path){
vector< vector< vector< double > > > Calibration_SiLi_E ;
// Calibration_SiLi_E[DetectorNumber][CrystalNumber][Order of Coeff]
- vector< vector< vector< double > > > Calibration_CsI_E ;
+ vector< vector< vector< double > > > Calibration_NaI_E ;
if (Path == "Simulation") { // Simulation case: data already calibrated
/* Calibration_Si_E_Order = 1;
Calibration_Si_T_Order = 1;
Calibration_SiLi_E_Order = 1;
- Calibration_CsI_E_Order = 1;
+ Calibration_NaI_E_Order = 1;
*/
vector<double> Coef;
// Order 0 Order 1
@@ -551,7 +546,7 @@ void TParisPhysics::ReadCalibrationFile(string Path){
Calibration_Si_Y_T.resize( m_NumberOfModule , StripLine) ;
Calibration_SiLi_E.resize( m_NumberOfModule , StripLine) ;
- Calibration_CsI_E .resize( m_NumberOfModule , StripLine) ;
+ Calibration_NaI_E .resize( m_NumberOfModule , StripLine) ;
}
}
@@ -652,7 +647,7 @@ void TParisPhysics::AddModuleSquare(double theta,
TVector3 U ;
// Vector V on Module Face (parallele to X Strip)
TVector3 V ;
- // Vector W normal to Module Face (pointing CsI)
+ // Vector W normal to Module Face (pointing NaI)
TVector3 W ;
// Vector position of Module Face center
TVector3 C ;
diff --git a/NPLib/Paris/TParisPhysics.h b/NPLib/Paris/TParisPhysics.h
index c43841ba8..602eedbf3 100644
--- a/NPLib/Paris/TParisPhysics.h
+++ b/NPLib/Paris/TParisPhysics.h
@@ -44,7 +44,7 @@ class TParisPhysics : public TObject, public NPA::VDetector{
public:
void Clear();
- void Clear(const Option_t*) {};
+ void Clear(const Option_t*) {};
public:
/////////////////////////////////////
@@ -150,29 +150,11 @@ class TParisPhysics : public TObject, public NPA::VDetector{
// Provide Physical Multiplicity
Int_t ParisEventMult;
- // Provide a Classification of Event
- //vector<int> EventType;
-
- // Telescope
- //vector<int> ModuleNumber;
-
// FirstStage
vector<double> ParisLaBr_E;
- //vector<double> FirstStage_T;
- //vector<int> FirstStage_X;
- //vector<int> FirstStage_Y;
// SecondStage
- vector<double> ParisCsI_E;
- //vector<double> SecondStage_T;
- //vector<int> SecondStage_N;
-
- /*
- // ThirdStage
- vector<double> ThirdStage_E;
- vector<double> ThirdStage_T;
- vector<int> ThirdStage_N;
- */
+ vector<double> ParisNaI_E;
// Physical Value
vector<double> ParisTotalEnergy;
diff --git a/NPSimulation/Paris/Paris.cc b/NPSimulation/Paris/Paris.cc
index 450b582de..6b61ffcd8 100644
--- a/NPSimulation/Paris/Paris.cc
+++ b/NPSimulation/Paris/Paris.cc
@@ -6,147 +6,606 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
- * Creation Date : 04/12/09 *
+ * Creation Date : November 2012 *
* Last update : *
*---------------------------------------------------------------------------*
- * Decription: This class manages different shapes of module for the Paris *
- * detector. It allows to have Paris geometries with an *
- * heterogeneous set of modules *
+ * Decription: *
+ * This class describe the Paris Silicon array *
+ * *
*---------------------------------------------------------------------------*
* Comment: *
* *
- * *
*****************************************************************************/
// C++ headers
-#include <fstream>
#include <sstream>
-#include <string>
#include <cmath>
+#include <limits>
+using namespace std;
-// NPTool headers
+//Geant4
+#include "G4VSolid.hh"
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4UnionSolid.hh"
+#include "G4SubtractionSolid.hh"
+#include "G4SDManager.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4Colour.hh"
+
+// NPS
#include "Paris.hh"
-#include "ParisPhoswich.hh"
-#include "ParisCluster.hh"
+using namespace PARIS;
-using namespace std;
+#include "CalorimeterScorers.hh"
+#include "MaterialManager.hh"
+// NPL
+#include "NPOptionManager.h"
+#include "RootOutput.h"
+
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+using namespace CLHEP;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Paris Specific Method
+Paris::Paris(){
+ m_Event = new TParisData();
+
+ // Orange
+ m_NaIVisAtt = new G4VisAttributes(G4Colour(1, 0.5, 0)) ;
+
+ // Blue
+ m_LaBr3VisAtt = new G4VisAttributes(G4Colour(0, 0.5, 1));
+
+ // Dark Grey
+ m_PMTVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1));
+
+ // Grey wireframe
+ m_PhoswichCasingVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5,0.2));
+ // White wireframe
+ m_ClusterCasingVisAtt = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7));
-Paris::Paris()
-{
+ m_LogicalPhoswich = 0;
+ m_LogicalCluster = 0;
+ m_NaIScorer = 0 ;
+ m_LaBr3Scorer = 0 ;
}
+Paris::~Paris(){
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Paris::AddCluster(G4ThreeVector Pos1, G4ThreeVector Pos2, G4ThreeVector Pos3, G4ThreeVector Pos4){
+ G4ThreeVector Pos=(Pos1+Pos2+Pos3+Pos4)/4.;
+ G4ThreeVector u = Pos1-Pos2;
+ G4ThreeVector v = Pos1-Pos4;
+ u = u.unit(); v = v.unit();
+ G4ThreeVector w = Pos.unit();
+ Pos = Pos + w*Length*0.5;
+
+ m_Type.push_back(1);
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(new G4RotationMatrix(u,v,w));
+}
-Paris::~Paris()
-{
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Paris::AddCluster(G4ThreeVector Pos, double beta_u, double beta_v, double beta_w){
+ double Theta = Pos.theta();
+ double Phi = Pos.phi();
+
+ // vector parallel to one axis of silicon plane
+ G4double ii = cos(Theta / rad) * cos(Phi / rad);
+ G4double jj = cos(Theta / rad) * sin(Phi / rad);
+ G4double kk = -sin(Theta / rad);
+ G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
+
+ G4ThreeVector w = Pos.unit();
+ G4ThreeVector u = w.cross(Y);
+ G4ThreeVector v = w.cross(u);
+ v = v.unit();
+ u = u.unit();
+
+ G4RotationMatrix* r = new G4RotationMatrix(u,v,w);
+ r->rotate(beta_u,u);
+ r->rotate(beta_v,v);
+ r->rotate(beta_w,w);
+
+ m_Type.push_back(1);
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(r);
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Paris::AddPhoswich(G4ThreeVector Pos1, G4ThreeVector Pos2, G4ThreeVector Pos3, G4ThreeVector Pos4){
+ G4ThreeVector Pos=(Pos1+Pos2+Pos3+Pos4)/4.;
+ G4ThreeVector u = Pos1-Pos2;
+ G4ThreeVector v = Pos1-Pos4;
+ u = u.unit(); v = v.unit();
+ G4ThreeVector w = Pos.unit();
+ Pos = Pos + w*Length*0.5;
+
+ m_Type.push_back(0);
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(new G4RotationMatrix(u,v,w));
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Paris::AddPhoswich(G4ThreeVector Pos, double beta_u, double beta_v, double beta_w){
+ double Theta = Pos.theta();
+ double Phi = Pos.phi();
+
+ // vector parallel to one axis of silicon plane
+ G4double ii = cos(Theta / rad) * cos(Phi / rad);
+ G4double jj = cos(Theta / rad) * sin(Phi / rad);
+ G4double kk = -sin(Theta / rad);
+ G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
+
+ G4ThreeVector w = Pos.unit();
+ G4ThreeVector u = w.cross(Y);
+ G4ThreeVector v = w.cross(u);
+ v = v.unit();
+ u = u.unit();
+
+ G4RotationMatrix* r = new G4RotationMatrix(u,v,w);
+ r->rotate(beta_u,u);
+ r->rotate(beta_v,v);
+ r->rotate(beta_w,w);
+
+ m_Type.push_back(0);
+ m_Pos.push_back(Pos);
+ m_Rot.push_back(r);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of VDetector class
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void Paris::ReadConfiguration(string Path)
-{
- // open configuration file
- ifstream ConfigFile;
- ConfigFile.open(Path.c_str());
-
- bool bParisPhoswich = false;
- bool bParisCluster = false;
-
- string LineBuffer;
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
- if (LineBuffer.compare(0, 14, "ParisPhoswich") == 0 && bParisPhoswich == false) {
- bParisPhoswich = true;
-
- // instantiate a new "detector" corresponding to the Square elements
- ParisModule* myDetector = new ParisPhoswich();
-
- // read part of the configuration file corresponding to square elements
- ConfigFile.close();
- myDetector->ReadConfiguration(Path);
- ConfigFile.open(Path.c_str());
-
- // ms_InterCoord comes from VDetector
- myDetector->SetInterCoordPointer(ms_InterCoord);
-
- // store ParisPhoswich "detector"
- m_Modules.push_back(myDetector);
+void Paris::ReadConfiguration(string Path){
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer, DataBuffer;
+
+ // A,B,C,D are the four corner of the detector
+
+ G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
+ G4ThreeVector A , B , C , D ;
+ G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
+ bool Type = true;
+
+ bool ReadingStatus = false;
+
+ bool check_A = false;
+ bool check_C = false;
+ bool check_B = false;
+ bool check_D = false;
+
+ bool check_Theta = false;
+ bool check_Phi = false;
+ bool check_R = false;
+
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+ if (LineBuffer.compare(0, 12, "ParisCluster") == 0) {
+ G4cout << "///" << G4endl ;
+ G4cout << "Cluster element found: " << G4endl ;
+ ReadingStatus = true ;
+ Type = true;
+ }
+
+ else if (LineBuffer.compare(0, 13, "ParisPhoswich") == 0) {
+ G4cout << "///" << G4endl ;
+ G4cout << "Phoswich element found: " << G4endl ;
+ ReadingStatus = true ;
+ Type = false;
+ }
+
+
+ while (ReadingStatus) {
+ ConfigFile >> DataBuffer;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
+
+ // Position method
+ else if (DataBuffer == "A=") {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax * mm ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay * mm ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az * mm ;
+
+ A = G4ThreeVector(Ax, Ay, Az);
+ G4cout << "Corner A position : " << A << G4endl;
+ }
+ else if (DataBuffer == "B=") {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx * mm ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By * mm ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz * mm ;
+
+ B = G4ThreeVector(Bx, By, Bz);
+ G4cout << "Corner B position : " << B << G4endl;
}
- else if (LineBuffer.compare(0, 12, "ParisCluster") == 0 && bParisCluster == false) {
- bParisCluster = true;
- // instantiate a new "detector" corresponding to the Trapezoid elements
- ParisModule* myDetector = new ParisCluster();
+ else if (DataBuffer == "C=") {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx * mm ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy * mm ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz * mm ;
+
+ C = G4ThreeVector(Cx, Cy, Cz);
+ G4cout << "Corner C position : " << C << G4endl;
+ }
+ else if (DataBuffer == "D=") {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx * mm ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy * mm ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz * mm ;
+
+ D = G4ThreeVector(Dx, Dy, Dz);
+ G4cout << "Corner D position : " << D << G4endl;
+ }
- // read part of the configuration file corresponding to trapezoid elements
- ConfigFile.close();
- myDetector->ReadConfiguration(Path);
- ConfigFile.open(Path.c_str());
+ // Angle method
+ else if (DataBuffer=="Theta=") {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta * deg;
+ G4cout << "Theta: " << Theta / deg << G4endl;
+ }
+ else if (DataBuffer=="Phi=") {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi * deg;
+ G4cout << "Phi: " << Phi / deg << G4endl;
+ }
+ else if (DataBuffer=="R=") {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R * mm;
+ G4cout << "R: " << R / mm << G4endl;
+ }
+ else if (DataBuffer=="Beta=") {
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w * deg ;
+ G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
+ }
+
+ else G4cout << "WARNING: Wrong Token, ParisCluster: Cluster Element not added" << G4endl;
+
+ // Add The previously define telescope
+ // With position method
+ if ((check_A && check_B && check_C && check_D) &&
+ !(check_Theta && check_Phi && check_R)) {
+ ReadingStatus = false;
+ check_A = false;
+ check_C = false;
+ check_B = false;
+ check_D = false;
+
+ if(Type)
+ AddCluster(A, B, C, D);
+ else
+ AddPhoswich(A, B, C, D);
+ }
- // ms_InterCoord comes from VDetector
- myDetector->SetInterCoordPointer(ms_InterCoord);
+ // With angle method
+ if ((check_Theta && check_Phi && check_R ) &&
+ !(check_A && check_B && check_C && check_D)) {
+ ReadingStatus = false;
+ check_Theta = false;
+ check_Phi = false;
+ check_R = false;
- // store ParisCluster "detector"
- m_Modules.push_back(myDetector);
+ G4ThreeVector Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
+
+ if(Type)
+ AddCluster(Pos, beta_u, beta_v, beta_w);
+ else
+ AddPhoswich(Pos, beta_u, beta_v, beta_w);
}
+ }
+ }
+}
- }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Construct detector and inialise sensitive part.
+// Called After DetecorConstruction::AddDetector Method
+void Paris::ConstructDetector(G4LogicalVolume* world){
+ unsigned int mysize = m_Pos.size();
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ if(m_Type[i])
+ new G4PVPlacement(G4Transform3D(*m_Rot[i], m_Pos[i]), ConstructCluster(), "ParisCluster", world, false, i+1);
+ else
+ new G4PVPlacement(G4Transform3D(*m_Rot[i], m_Pos[i]), ConstructPhoswich(), "ParisPhoswich", world, false, i+1);
+ }
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Paris::ConstructPhoswich(){
+
+ if(!m_LogicalPhoswich){
+ G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
+ G4Material* Alu = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4Material* NaI = MaterialManager::getInstance()->GetMaterialFromLibrary("NaI");
+ G4Material* LaBr3 = MaterialManager::getInstance()->GetMaterialFromLibrary("LaBr3");
+
+ G4Box* solidParisPhoswich = new G4Box("SolidParisPhoswich", 0.5*PhoswichCasingWidth, 0.5*PhoswichCasingWidth, 0.5*PhoswichCasingLenth+0.5*PMTLength);
+ m_LogicalPhoswich = new G4LogicalVolume(solidParisPhoswich, Vacuum, "LogicParisPhoswich", 0, 0, 0);
+ m_LogicalPhoswich->SetVisAttributes(G4VisAttributes::Invisible);
+
+ // Phoswich construction
+ // Casing
+ G4Box* solidCasing = new G4Box("SolidParisPhoswichCasing", 0.5*PhoswichCasingWidth, 0.5*PhoswichCasingWidth, 0.5*PhoswichCasingLenth);
+ G4LogicalVolume* LogicalCasing = new G4LogicalVolume(solidCasing, Alu, "LogicParisPhoswichCasing", 0, 0, 0);
+ LogicalCasing->SetVisAttributes(m_PhoswichCasingVisAtt);
+
+ new G4PVPlacement(0,
+ G4ThreeVector(0,0,-PMTLength),
+ LogicalCasing,
+ "ParisPhoswich_Casing",
+ m_LogicalPhoswich,
+ false,
+ 0);
+
+ // PMT
+ G4Tubs* solidPMT = new G4Tubs("SolidPMT",0,PMTRadius,0.5*PMTLength,0,360*deg);
+ G4LogicalVolume* LogicalPMT = new G4LogicalVolume(solidPMT, Alu, "LogicPMT" ,0,0,0);
+ new G4PVPlacement(0,
+ G4ThreeVector(0,0,0.5*Length-PMTLength),
+ LogicalPMT,
+ "ParisPhoswich_PMT",
+ m_LogicalPhoswich,
+ false,
+ 0);
+ LogicalPMT->SetVisAttributes(m_PMTVisAtt);
+
+ // LaBr3
+ G4Box* solidLaBr3Stage = new G4Box("solidLaBr3Stage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*LaBr3Thickness);
+ G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, LaBr3, "logicLaBr3Stage", 0, 0, 0);
+ G4ThreeVector positionLaBr3Stage = G4ThreeVector(0, 0, LaBr3Stage_PosZ);
+
+ new G4PVPlacement(0,
+ positionLaBr3Stage,
+ logicLaBr3Stage,
+ "ParisPhoswich_LaBr3Stage",
+ LogicalCasing,
+ false,
+ 0);
+
+ // Set LaBr3 sensible
+ logicLaBr3Stage->SetSensitiveDetector(m_LaBr3Scorer);
+
+ // Visualisation of LaBr3 stage
+ logicLaBr3Stage->SetVisAttributes(m_LaBr3VisAtt);
+
+ // NaI
+ G4ThreeVector positionNaIStage = G4ThreeVector(0, 0, NaIStage_PosZ);
+
+ G4Box* solidNaIStage = new G4Box("solidNaIStage", 0.5*NaIFace, 0.5*NaIFace, 0.5*NaIThickness);
+ G4LogicalVolume* logicNaIStage = new G4LogicalVolume(solidNaIStage, NaI, "logicNaIStage", 0, 0, 0);
+
+ new G4PVPlacement(0,
+ positionNaIStage,
+ logicNaIStage,
+ "ParisPhoswich_NaIStage",
+ LogicalCasing,
+ false,
+ 0);
+
+ // Set NaI sensible
+ logicNaIStage->SetSensitiveDetector(m_NaIScorer);
+
+ // Visualisation of the NaI stage
+ logicNaIStage->SetVisAttributes(m_NaIVisAtt);
+ }
+
+ return m_LogicalPhoswich;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Paris::ConstructCluster(){
+
+ if(!m_LogicalCluster){
+ G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
+ G4Material* Alu = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+
+ // Mother Volume
+ G4Box* solidParisCluster1 = new G4Box("SolidParisCluster1", 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
+ G4Box* solidParisCluster2 = new G4Box("SolidParisCluster2", 0.5*ClusterFrameWidth, 0.5*ClusterFrameWidth, 0.5*ClusterFrameLength);
+
+ // The cluster is made of a metal frame that hold the phoswich via the pmt
+ G4ThreeVector FramePos(0,0,Length*0.5-PMTLength+ClusterFrameLength*0.5+ClusterOffset);
+ G4UnionSolid* solidParisCluster = new G4UnionSolid("solidParisCluster",solidParisCluster1,solidParisCluster2,0,FramePos);
+ m_LogicalCluster = new G4LogicalVolume(solidParisCluster, Vacuum, "LogicSolidParisCluster", 0, 0, 0);
+ m_LogicalCluster->SetVisAttributes(G4VisAttributes::Invisible);
+
+ if(!m_LogicalPhoswich) ConstructPhoswich();
+
+ // The frame is a plate with hole for the PMT
+ G4VSolid* solidClusterFrame = solidParisCluster2;
+ G4ThreeVector Origin(-PhoswichCasingWidth,-PhoswichCasingWidth,0.5*PMTLength);
+ G4Tubs* solidPMT = new G4Tubs("SolidPMT",0,PMTRadius,10*PMTLength,0,360*deg);
+
+ // A cluster is a 3 by 3 aggregat of phoswich
+ unsigned int PhoswichNbr = 1;
+ for(unsigned int i = 0 ; i < 3 ; i++){
+ for(unsigned int j = 0 ; j <3 ; j++){
+ G4ThreeVector Pos = Origin + G4ThreeVector(i*PhoswichCasingWidth,j*PhoswichCasingWidth,0);
+ new G4PVPlacement(0,
+ Pos,
+ m_LogicalPhoswich,
+ "Paris_Phoswich",
+ m_LogicalCluster,
+ false,
+ PhoswichNbr++);
+
+ // make room for the PMT in the cluster frame
+ solidClusterFrame = new G4SubtractionSolid("ClusterFrame",solidClusterFrame,solidPMT,0,Pos);
+ }
+ }
-// Construct detector and initialize sensitive part.
-// Called After DetectorConstruction::AddDetector Method
-void Paris::ConstructDetector(G4LogicalVolume* world)
-{
- // loop on sub-detectors belonging to Paris
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ConstructDetector(world);
-}
+ G4LogicalVolume* LogicalClusterFrame = new G4LogicalVolume(solidClusterFrame, Alu, "LogicSolidClusterFrame", 0, 0, 0);
+ LogicalClusterFrame->SetVisAttributes(m_ClusterCasingVisAtt);
+ new G4PVPlacement(0,
+ FramePos,
+ LogicalClusterFrame,
+ "Paris_ClusterFrame",
+ m_LogicalCluster,
+ false,
+ 0);
+ }
+ return m_LogicalCluster;
-// Connect the ParisData class to the output TTree
-// of the simulation
-void Paris::InitializeRootOutput()
-{
- // loop on sub-detectors belonging to Paris
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeRootOutput();
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void Paris::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ pTree->Branch("Paris", "TParisData", &m_Event) ;
+ pTree->SetBranchAddress("Paris", &m_Event) ;
+}
-// Initialize all scorers necessary for each detector
-void Paris::InitializeScorers()
-{
- // loop on sub-detectors belonging to Paris
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeScorers();
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void Paris::ReadSensitive(const G4Event* event){
+ m_Event->Clear();
+
+ ///////////
+ // LaBr3
+ G4THitsMap<G4double*>* LaBr3HitMap;
+ std::map<G4int, G4double**>::iterator LaBr3_itr;
+
+ G4int LaBr3CollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Paris_LaBr3Scorer/ParisLaBr3");
+ LaBr3HitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(LaBr3CollectionID));
+
+ // Loop on the LaBr3 map
+ for (LaBr3_itr = LaBr3HitMap->GetMap()->begin() ; LaBr3_itr != LaBr3HitMap->GetMap()->end() ; LaBr3_itr++){
+
+ G4double* Info = *(LaBr3_itr->second);
+
+ double Energy = Info[0];
+
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int PhoswichNbr = (int) Info[2];
+ int ClusterNbr = (int) Info[3];
+
+ m_Event->SetParisLaBr3E(ClusterNbr,PhoswichNbr,Energy);
+ m_Event->SetParisLaBr3T(ClusterNbr,PhoswichNbr,Time);
+ }
+ }
+ // clear map for next event
+ LaBr3HitMap->clear();
+
+ ///////////
+ // NaI
+ G4THitsMap<G4double*>* NaIHitMap;
+ std::map<G4int, G4double**>::iterator NaI_itr;
+
+ G4int NaICollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Paris_NaIScorer/ParisNaI");
+ NaIHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(NaICollectionID));
+
+ // Loop on the NaI map
+ for (NaI_itr = NaIHitMap->GetMap()->begin() ; NaI_itr != NaIHitMap->GetMap()->end() ; NaI_itr++){
+
+ G4double* Info = *(NaI_itr->second);
+
+ double Energy = Info[0];
+
+ if(Energy>EnergyThreshold){
+ double Time = Info[1];
+ int PhoswichNbr = (int) Info[2];
+ int ClusterNbr = (int) Info[3];
+
+ m_Event->SetParisNaIE(ClusterNbr,PhoswichNbr,Energy);
+ m_Event->SetParisNaIT(ClusterNbr,PhoswichNbr,Time);
+ }
+ }
+
+ // clear map for next event
+ NaIHitMap->clear();
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Paris::InitializeScorers(){
+ vector<G4int> NestingLevel;
+ NestingLevel.push_back(2);
+ NestingLevel.push_back(3);
+ // LaBr3 Associate Scorer
+ bool already_exist = false;
+ m_LaBr3Scorer = CheckScorer("Paris_LaBr3Scorer",already_exist);
+
+ // if the scorer were created previously nothing else need to be made
+ if(already_exist) return;
+
+ G4VPrimitiveScorer* LaBr3Scorer =
+ new CALORIMETERSCORERS::PS_Calorimeter("ParisLaBr3",NestingLevel);
+ //and register it to the multifunctionnal detector
+ m_LaBr3Scorer->RegisterPrimitive(LaBr3Scorer);
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3Scorer) ;
+
+
+ // NaI Associate Scorer
+ already_exist = false;
+ m_NaIScorer = CheckScorer("Paris_NaIScorer",already_exist);
+
+ // if the scorer were created previously nothing else need to be made
+ if(already_exist) return;
+
+ G4VPrimitiveScorer* NaIScorer =
+ new CALORIMETERSCORERS::PS_Calorimeter("ParisNaI",NestingLevel);
+ //and register it to the multifunctionnal detector
+ m_NaIScorer->RegisterPrimitive(NaIScorer);
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_NaIScorer) ;
-// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAction
-void Paris::ReadSensitive(const G4Event* event)
-{
- // Before looping on each sub-detector, clear the static variable
- // ms_InterCoord
- // This is done on the first element of the m_Modules vector.
- // This should be done here since this variable (of type TIneractionCoordinates)
- // deals with multiplicity of events > 1.
- m_Modules[0]->GetInterCoordPointer()->Clear();
-
- // We do the same for the static variable ms_Event
- m_Modules[0]->GetEventPointer()->Clear();
-
- // loop on sub-detectors belonging to Paris
- int nbDetectors = m_Modules.size();
- for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ReadSensitive(event);
}
diff --git a/NPSimulation/Paris/Paris.hh b/NPSimulation/Paris/Paris.hh
index 3ee1aa4f5..b885e1063 100644
--- a/NPSimulation/Paris/Paris.hh
+++ b/NPSimulation/Paris/Paris.hh
@@ -8,67 +8,169 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
* *
- * Creation Date : 04/12/09 *
+ * Creation Date : November 2012 *
* Last update : *
*---------------------------------------------------------------------------*
- * Decription: This class manages different shapes of module for the Paris *
- * detector. It allows to have Paris geometries with an *
- * heterogeneous set of modules *
+ * Decription: *
+ * This class describe the Paris Silicon detector *
+ * *
*---------------------------------------------------------------------------*
* Comment: *
* *
- * *
*****************************************************************************/
+// C++ header
+#include <string>
+#include <vector>
+// G4 header defining G4 types
+#include "globals.hh"
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4VisAttributes.hh"
+#include "G4MultiFunctionalDetector.hh"
-// C++ headers
-#include <vector>
-
-// NPTool header
+// NPSimulation header
#include "VDetector.hh"
-#include "ParisModule.hh"
+
+// NPLib
+#include "TParisData.h"
using namespace std;
-#include "CLHEP/Units/SystemOfUnits.h"
using namespace CLHEP;
-
-class Paris : public VDetector
-{
- ////////////////////////////////////////////////////
- /////// Default Constructor and Destructor /////////
- ////////////////////////////////////////////////////
+class Paris : public VDetector{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
public:
- Paris();
- virtual ~Paris();
-
- ////////////////////////////////////////////////////
- ///////// Inherite from VDetector class ///////////
- ////////////////////////////////////////////////////
+ Paris() ;
+ ~Paris() ;
+
+ ////////////////////////////////////////////////////
+ //////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
-
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world);
+ // To add a Cluster detector
+ // By corner position
+ void AddCluster(G4ThreeVector Pos1, G4ThreeVector Pos2, G4ThreeVector Pos3, G4ThreeVector Pos4);
+ // By Center position
+ void AddCluster(G4ThreeVector Pos, double beta_u=0, double beta_v=0, double beta_w=0);
+ // To add a Phoswich detector
+ // By corner position
+ void AddPhoswich(G4ThreeVector Pos1,G4ThreeVector Pos2,G4ThreeVector Pos3,G4ThreeVector Pos4);
+ // By Center position
+ void AddPhoswich(G4ThreeVector Pos, double beta_u=0, double beta_v=0, double beta_w=0);
+
+ // Return a logical volume of the detector
+ G4LogicalVolume* ConstructPhoswich();
+ G4LogicalVolume* ConstructCluster();
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput();
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
+private: // Guarranty that each volume is created only once
+ G4LogicalVolume* m_LogicalPhoswich;
+ G4LogicalVolume* m_LogicalCluster;
+ ////////////////////////////////////////////////////
+ ///////// Inherite from VDetector class ///////////
+ ////////////////////////////////////////////////////
public:
- // Initialize all scorers necessary for each detector
- void InitializeScorers();
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path) ;
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world) ;
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event) ;
+
+
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+private:
+ TParisData* m_Event ;
+
+ ////////////////////////////////////////////////////
+ ///////////////// Scorer Related ///////////////////
+ ////////////////////////////////////////////////////
+
+private:
+ // Initialize all Scorer
+ void InitializeScorers() ;
+
+ // Scorer Associate to the Silicon
+ G4MultiFunctionalDetector* m_NaIScorer ;
+ G4MultiFunctionalDetector* m_LaBr3Scorer ;
private:
- vector<ParisModule*> m_Modules;
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private:
+ // True if the detector is a Cluster, false if a single Phoswich
+ vector<bool> m_Type;
+
+ // Detector position
+ vector<G4ThreeVector> m_Pos;
+ vector<G4RotationMatrix*> m_Rot;
+
+private:/// Visualisation Attribute:
+ G4VisAttributes* m_NaIVisAtt ;
+ G4VisAttributes* m_LaBr3VisAtt;
+ G4VisAttributes* m_PhoswichCasingVisAtt ;
+ G4VisAttributes* m_ClusterCasingVisAtt ;
+ G4VisAttributes* m_PMTVisAtt;
};
+
+namespace PARIS{
+
+ // Resolution and Threshold
+ const G4double EnergyThreshold = 100*keV;
+ const G4double ResoFirstStage = 0.0085; // = 3% at .662 MeV of Resolution // Unit is MeV/2.35
+ const G4double ResoSecondStage = 0.0366; // = 13% at .662 MeV of resolution // Unit is MeV/2.35
+ const G4double ResoTime = 0.212765957;// = 500ps // Unit is ns/2.35
+
+ // Phoswich geometry
+ const G4double PhoswichCasingLenth = 226*mm;
+ const G4double PhoswichCasingWidth = 56*mm;
+
+ // LaBr3
+ const G4double LaBr3Face = 50.8*mm;
+ const G4double LaBr3Thickness = 50.8*mm; // for phoswich 2x2x3
+ const G4double InterSpace = 0*mm;
+
+ // NaI
+ const G4double NaIFace = LaBr3Face;
+ const G4double NaIThickness = 152.4*mm; // for phoswich 2x2x6
+
+ // PMT geometry
+ const G4double PMTRadius = 26*mm;
+ const G4double PMTLength = 150*mm; // PMT + connector
+
+ // Starting form the LaBr3 and going to the NaI
+ const G4double CasingThickness = 0.5*(PhoswichCasingWidth-LaBr3Face);
+ const G4double LaBr3Stage_PosZ =CasingThickness+PhoswichCasingLenth* -0.5 + 0.5*LaBr3Thickness;
+ const G4double NaIStage_PosZ = CasingThickness+ PhoswichCasingLenth* -0.5 + LaBr3Thickness + 0.5*NaIThickness;
+
+ // Cluster geometry
+ const G4double ClusterFrameWidth = 188*mm;
+ const G4double ClusterFrameLength = 15*mm;
+ const G4double ClusterOffset = 30*mm;
+ const G4double FaceFront = 3*PhoswichCasingWidth;
+ const G4double FaceBack = FaceFront;
+ const G4double Length = PhoswichCasingLenth+PMTLength;
+
+
+}
+
#endif
diff --git a/NPSimulation/Paris/ParisCluster.cc b/NPSimulation/Paris/ParisCluster.cc
deleted file mode 100644
index 752909d8f..000000000
--- a/NPSimulation/Paris/ParisCluster.cc
+++ /dev/null
@@ -1,1195 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : 20/07/11 *
- *---------------------------------------------------------------------------*
- * Decription: Define a cluster of 9 phoswich modules for the Paris detector *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-// C++ headers
-#include <sstream>
-#include <string>
-#include <cmath>
-
-// G4 Geometry headers
-#include "G4Trd.hh"
-#include "G4Box.hh"
-#include "G4Trap.hh"
-
-// G4 various headers
-#include "G4Material.hh"
-#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
-#include "G4RotationMatrix.hh"
-#include "G4Transform3D.hh"
-#include "G4PVPlacement.hh"
-#include "G4PVDivision.hh"
-
-// G4 sensitive
-#include "G4SDManager.hh"
-#include "G4MultiFunctionalDetector.hh"
-
-// NPTool headers
-#include "ParisCluster.hh"
-#include "ParisScorers.hh"
-#include "ObsoleteGeneralScorers.hh"
-#include "MaterialManager.hh"
-#include "RootOutput.h"
-#include "VDetector.hh"
-// CLHEP
-#include "CLHEP/Random/RandGauss.h"
-
-using namespace std;
-using namespace CLHEP;
-using namespace PARISCLUSTER;
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-ParisCluster::ParisCluster()
-{
- ms_InterCoord = 0;
- m_LaBr3StageScorer=0;
- m_CsIStageScorer=0;
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-ParisCluster::~ParisCluster()
-{
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisCluster::AddModule(G4ThreeVector X1_Y1,
- G4ThreeVector X128_Y1,
- G4ThreeVector X1_Y128,
- G4ThreeVector X128_Y128)
-{
- m_DefinitionType.push_back(true);
-
- m_X1_Y1.push_back(X1_Y1);
- m_X128_Y1.push_back(X128_Y1);
- m_X1_Y128.push_back(X1_Y128);
- m_X128_Y128.push_back(X128_Y128);
-
- m_R.push_back(0);
- m_Theta.push_back(0);
- m_Phi.push_back(0);
- m_beta_u.push_back(0);
- m_beta_v.push_back(0);
- m_beta_w.push_back(0);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisCluster::AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w)
-{
- G4ThreeVector empty = G4ThreeVector(0, 0, 0);
-
- m_DefinitionType.push_back(false);
-
- m_X1_Y1.push_back(empty);
- m_X128_Y1.push_back(empty);
- m_X1_Y128.push_back(empty);
- m_X128_Y128.push_back(empty);
-
- m_R.push_back(R);
- m_Theta.push_back(Theta);
- m_Phi.push_back(Phi);
- m_beta_u.push_back(beta_u);
- m_beta_v.push_back(beta_v);
- m_beta_w.push_back(beta_w);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisCluster::VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world)
-{
- G4double NbrTelescopes = DetecNumber;
- G4String DetectorNumber;
- ostringstream Number;
- Number << NbrTelescopes;
- DetectorNumber = Number.str();
-
- G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
- G4Material* NaI = MaterialManager::getInstance()->GetMaterialFromLibrary("NaI");
- G4Material* LaBr3 = MaterialManager::getInstance()->GetMaterialFromLibrary("LaBr3");
-
- ////////////////////////////////////////////////////////////////
- ////////////// Starting Volume Definition //////////////////////
- ////////////////////////////////////////////////////////////////
- G4String Name = "ParisCluster" + DetectorNumber ;
-
- // Mother Volume
- G4Box* solidParisCluster = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
- G4LogicalVolume* logicParisCluster = new G4LogicalVolume(solidParisCluster, Vacuum, Name, 0, 0, 0);
-
- new G4PVPlacement(G4Transform3D(*MMrot, MMpos) ,
- logicParisCluster ,
- Name ,
- world ,
- false ,
- 0);
-
- G4VisAttributes* ClusterVisAtt = new G4VisAttributes(G4Colour(0.90, 0.90, 0.90));
- ClusterVisAtt->SetForceWireframe(true);
- logicParisCluster->SetVisAttributes(ClusterVisAtt);
- // Phoswich volume construction (empty volume)
-
- G4ThreeVector positionPhoSwStage = G4ThreeVector(LaBr3Face+InterSpace, LaBr3Face+InterSpace, 0.);
- G4Box* solidPhoSwStage = new G4Box(Name+"_PhoSwStage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*Length);
- G4LogicalVolume* logicPhoSwStage = new G4LogicalVolume(solidPhoSwStage, Vacuum, "logicPhoSwStage", 0, 0, 0);
- logicPhoSwStage->SetVisAttributes(G4VisAttributes::Invisible);
-
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false,8);
- // 0);
-
- positionPhoSwStage = G4ThreeVector(0., LaBr3Face+InterSpace, 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false,7);
- // 1);
-
- positionPhoSwStage = G4ThreeVector(-(LaBr3Face+InterSpace), LaBr3Face+InterSpace, 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false,6);
- //2);
-
- positionPhoSwStage = G4ThreeVector(LaBr3Face+InterSpace, 0., 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false,5);
- //3);
-
- positionPhoSwStage = G4ThreeVector(0., 0., 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false, 4);
- //4);
-
- positionPhoSwStage = G4ThreeVector(-(LaBr3Face+InterSpace), 0., 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false, 3);
- //5);
-
- positionPhoSwStage = G4ThreeVector(LaBr3Face+InterSpace, -(LaBr3Face+InterSpace), 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false,2);
- //6);
-
- positionPhoSwStage = G4ThreeVector(0.,-(LaBr3Face+InterSpace), 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false, 1);
- //7);
-
- positionPhoSwStage = G4ThreeVector(-(LaBr3Face+InterSpace),-(LaBr3Face+InterSpace), 0.);
- new G4PVPlacement(0,
- positionPhoSwStage,
- logicPhoSwStage,
- Name+"_PhoSwStage",
- logicParisCluster,
- false, 0);
- //8);
-
-
- // LaBr3
- G4ThreeVector positionLaBr3Stage = G4ThreeVector(0., 0., LaBr3Stage_PosZ);
-
- G4Box* solidLaBr3Stage = new G4Box(Name+"_LaBr3Stage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*LaBr3Thickness);
- G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, LaBr3, "logicLaBr3Stage", 0, 0, 0);
-
- new G4PVPlacement(0,
- positionLaBr3Stage,
- logicLaBr3Stage,
- Name+"_LaBr3Stage",
- logicPhoSwStage,
- false,
- 0);
-
-
-
-
- // Set LaBr3 sensible
- logicLaBr3Stage->SetSensitiveDetector(m_LaBr3StageScorer);
-
- // Visualisation of LaBr3Stage Strip
- G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(0., 0.5, 1.));
- logicLaBr3Stage->SetVisAttributes(LaBr3VisAtt);
-
- // CsI or NaI
- G4ThreeVector positionCsIStage = G4ThreeVector(0.,0., CsIStage_PosZ);
-
- G4Box* solidCsIStage = new G4Box(Name+"_CsIStage", 0.5*CsIFace, 0.5*CsIFace, 0.5*CsIThickness);
- //G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, CsI, "logicCsIStage", 0, 0, 0);
- G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, NaI, "logicCsIStage", 0, 0, 0);
-
- new G4PVPlacement(0,
- positionCsIStage,
- logicCsIStage,
- Name + "_CsIStage",
- logicPhoSwStage,
- false,
- 0);
-
- // Set CsI sensible
- logicCsIStage->SetSensitiveDetector(m_CsIStageScorer);
-
- // Visualisation of CsIStage Strip
- G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1., 0.25, 0.));
- logicCsIStage->SetVisAttributes(CsIVisAtt);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Virtual Method of VDetector class
-
-// Read stream at Configfile to pick-up parameters of detector (Position,...)
-// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void ParisCluster::ReadConfiguration(string Path)
-{
- ifstream ConfigFile;
- ConfigFile.open(Path.c_str());
- string LineBuffer, DataBuffer;
-
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
- G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
- G4ThreeVector A , B , C , D ;
- G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
-
- bool ReadingStatus = false;
-
- bool check_A = false;
- bool check_C = false;
- bool check_B = false;
- bool check_D = false;
-
- bool check_Theta = false;
- bool check_Phi = false;
- bool check_R = false;
-
- bool checkVis = false;
-
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
- if (LineBuffer.compare(0, 12, "ParisCluster") == 0) {
- G4cout << "///" << G4endl ;
- G4cout << "Cluster element found: " << G4endl ;
- ReadingStatus = true ;
- }
-
- while (ReadingStatus) {
- ConfigFile >> DataBuffer;
- // Comment Line
- if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
-
- // Position method
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- Ax = atof(DataBuffer.c_str()) ;
- Ax = Ax * mm ;
- ConfigFile >> DataBuffer ;
- Ay = atof(DataBuffer.c_str()) ;
- Ay = Ay * mm ;
- ConfigFile >> DataBuffer ;
- Az = atof(DataBuffer.c_str()) ;
- Az = Az * mm ;
-
- A = G4ThreeVector(Ax, Ay, Az);
- G4cout << "X1 Y1 corner position : " << A << G4endl;
- }
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- Bx = atof(DataBuffer.c_str()) ;
- Bx = Bx * mm ;
- ConfigFile >> DataBuffer ;
- By = atof(DataBuffer.c_str()) ;
- By = By * mm ;
- ConfigFile >> DataBuffer ;
- Bz = atof(DataBuffer.c_str()) ;
- Bz = Bz * mm ;
-
- B = G4ThreeVector(Bx, By, Bz);
- G4cout << "X128 Y1 corner position : " << B << G4endl;
- }
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- Cx = atof(DataBuffer.c_str()) ;
- Cx = Cx * mm ;
- ConfigFile >> DataBuffer ;
- Cy = atof(DataBuffer.c_str()) ;
- Cy = Cy * mm ;
- ConfigFile >> DataBuffer ;
- Cz = atof(DataBuffer.c_str()) ;
- Cz = Cz * mm ;
-
- C = G4ThreeVector(Cx, Cy, Cz);
- G4cout << "X1 Y128 corner position : " << C << G4endl;
- }
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- Dx = atof(DataBuffer.c_str()) ;
- Dx = Dx * mm ;
- ConfigFile >> DataBuffer ;
- Dy = atof(DataBuffer.c_str()) ;
- Dy = Dy * mm ;
- ConfigFile >> DataBuffer ;
- Dz = atof(DataBuffer.c_str()) ;
- Dz = Dz * mm ;
-
- D = G4ThreeVector(Dx, Dy, Dz);
- G4cout << "X128 Y128 corner position : " << D << G4endl;
- }
-
- // Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- Theta = Theta * deg;
- G4cout << "Theta: " << Theta / deg << G4endl;
- }
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- Phi = Phi * deg;
- G4cout << "Phi: " << Phi / deg << G4endl;
- }
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
- R = R * mm;
- G4cout << "R: " << R / mm << G4endl;
- }
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- beta_u = beta_u * deg ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- beta_v = beta_v * deg ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- beta_w = beta_w * deg ;
- G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
- }
-
- else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
- checkVis = true ;
- ConfigFile >> DataBuffer;
- if (DataBuffer.compare(0, 3, "all") == 0) m_non_sensitive_part_visiualisation = true;
- }
-
- else G4cout << "WARNING: Wrong Token, ParisCluster: Cluster Element not added" << G4endl;
-
- // Add The previously define telescope
- // With position method
- if ((check_A && check_B && check_C && check_D && checkVis) &&
- !(check_Theta && check_Phi && check_R)) {
- ReadingStatus = false;
- check_A = false;
- check_C = false;
- check_B = false;
- check_D = false;
- checkVis = false;
-
- AddModule(A, B, C, D);
- }
-
- // With angle method
- if ((check_Theta && check_Phi && check_R && checkVis) &&
- !(check_A && check_B && check_C && check_D)) {
- ReadingStatus = false;
- check_Theta = false;
- check_Phi = false;
- check_R = false;
- checkVis = false;
-
- AddModule(R, Theta, Phi, beta_u, beta_v, beta_w);
- }
- }
- }
-}
-
-
-
-// Construct detector and inialise sensitive part.
-// Called After DetecorConstruction::AddDetector Method
-void ParisCluster::ConstructDetector(G4LogicalVolume* world)
-{
- G4RotationMatrix* MMrot = NULL ;
- G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
-
- G4int NumberOfTelescope = m_DefinitionType.size() ;
-
- for (G4int i = 0; i < NumberOfTelescope; i++) {
- // By Point
- if (m_DefinitionType[i]) {
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- MMu = m_X128_Y1[i] - m_X1_Y1[i];
- MMu = MMu.unit();
-
- MMv = m_X1_Y128[i] - m_X1_Y1[i];
- MMv = MMv.unit();
-
- MMw = MMu.cross(MMv);
-// if (MMw.z() > 0) MMw = MMv.cross(MMu) ;
- MMw = MMw.unit();
-
- MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4;
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + MMCenter;
- }
-
- // By Angle
- else {
- G4double Theta = m_Theta[i] ;
- G4double Phi = m_Phi[i] ;
-
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- // Phi is angle between X axis and projection in (X,Y) plan
- // Theta is angle between position vector and z axis
- G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
- G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
- G4double wZ = m_R[i] * cos(Theta / rad);
- MMw = G4ThreeVector(wX, wY, wZ);
-
- // vector corresponding to the center of the module
- G4ThreeVector CT = MMw;
-
- // vector parallel to one axis of silicon plane
- G4double ii = cos(Theta / rad) * cos(Phi / rad);
- G4double jj = cos(Theta / rad) * sin(Phi / rad);
- G4double kk = -sin(Theta / rad);
- G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
-
- MMw = MMw.unit();
- MMu = MMw.cross(Y);
- MMv = MMw.cross(MMu);
- MMv = MMv.unit();
- MMu = MMu.unit();
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // Telescope is rotate of Beta angle around MMv axis.
- MMrot->rotate(m_beta_u[i], MMu);
- MMrot->rotate(m_beta_v[i], MMv);
- MMrot->rotate(m_beta_w[i], MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + CT ;
- }
-
- VolumeMaker(i + 1, MMpos, MMrot, world);
- }
-
- delete MMrot ;
-}
-
-
-
-// Connect the ParisData class to the output TTree
-// of the simulation
-void ParisCluster::InitializeRootOutput()
-{
-}
-
-
-
-// Set the TinteractionCoordinates object from VDetector to the present class
-void ParisCluster::SetInterCoordPointer(TInteractionCoordinates* interCoord)
-{
- ms_InterCoord = interCoord;
-}
-
-
-
-// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAvtion
-void ParisCluster::ReadSensitive(const G4Event* event)
-{
- //////////////////////////////////////////////////////////////////////////////////////
- //////////////////////// Used to Read Event Map of detector //////////////////////////
- //////////////////////////////////////////////////////////////////////////////////////
-
- momentum = event->GetPrimaryVertex()->GetPrimary()->GetMomentum();
- G4double EGamma = momentum.getR(); // for photon E=p
- EGamma *= 1;
- //G4double EGammaMin = EGamma-4*ResoFirstStage;
- //G4double EGammaMax = EGamma+4*ResoFirstStage;
-
- // First Stage (LaBr3)
- std::map<G4int, G4int*>::iterator DetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CrystalNumber_itr;
- std::map<G4int, G4double*>::iterator Energy_itr;
- std::map<G4int, G4double*>::iterator Time_itr;
- //std::map<G4int, G4double*>::iterator X_itr;
- //std::map<G4int, G4double*>::iterator Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_X_itr;
- //std::map<G4int, G4double*>::iterator Pos_Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_Z_itr;
- //std::map<G4int, G4double*>::iterator Ang_Theta_itr;
- //std::map<G4int, G4double*>::iterator Ang_Phi_itr;
-
- G4THitsMap<G4int>* DetectorNumberHitMap;
- G4THitsMap<G4int>* CrystalNumberHitMap;
- G4THitsMap<G4double>* EnergyHitMap;
- G4THitsMap<G4double>* TimeHitMap;
- // G4THitsMap<G4double>* XHitMap;
- // G4THitsMap<G4double>* YHitMap;
- //G4THitsMap<G4double>* PosXHitMap;
- //G4THitsMap<G4double>* PosYHitMap;
- //G4THitsMap<G4double>* PosZHitMap;
- //G4THitsMap<G4double>* AngThetaHitMap;
- //G4THitsMap<G4double>* AngPhiHitMap;
-
- // NULL pointer are given to avoid warning at compilation
- // Second Stage (CsI)
- std::map<G4int, G4int*>::iterator CsIDetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CsICrystalNumber_itr; // added by Anna
- std::map<G4int, G4double*>::iterator CsIStageEnergy_itr ;
- G4THitsMap<G4int>* CsIDetectorNumberHitMap = NULL ;
- G4THitsMap<G4int>* CsICrystalNumberHitMap = NULL ; // added by Anna
- G4THitsMap<G4double>* CsIStageEnergyHitMap = NULL ;
-
-
- // Read the Scorer associate to the LaBr
- //Detector Number
- G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/DetectorNumber") ;
- DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ;
- DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ;
-
- //Crystal Number
- G4int CrystDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/CrystalNumber") ;
- CrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CrystDetCollectionID)) ;
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin() ;
-
- //Energy
- G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/StripEnergy") ;
- EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
- Energy_itr = EnergyHitMap->GetMap()->begin() ;
-
- //Time of Flight
- G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/StripTime") ;
- TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ;
- Time_itr = TimeHitMap->GetMap()->begin() ;
-
- /*
- //Strip Number X
- G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/StripIDFront") ;
- XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ;
- X_itr = XHitMap->GetMap()->begin() ;
-
- //Strip Number Y
- G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/StripIDBack");
- YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ;
- Y_itr = YHitMap->GetMap()->begin() ;
-
-
- //Interaction Coordinate X
- G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/InterCoordX") ;
- PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ;
- Pos_X_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Y
- G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/InterCoordY") ;
- PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ;
- Pos_Y_itr = PosYHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Z
- G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/InterCoordZ") ;
- PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ;
- Pos_Z_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Theta
- G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/InterCoordAngTheta") ;
- AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ;
- Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Phi
- G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisCluster/InterCoordAngPhi") ;
- AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ;
- Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ;
- */
-
-
- // Read the Scorer associate to the SecondStage
- //CsI Detector Number
- G4int CsIDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisCluster/CsIDetectorNumber") ;
- CsIDetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsIDetCollectionID)) ;
- CsIDetectorNumber_itr = CsIDetectorNumberHitMap->GetMap()->begin() ;
-
- //CsI Crystal Number // added by Anna
- G4int CsICryCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisCluster/CsICrystalNumber") ; // added by Anna
- CsICrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICryCollectionID)) ; // added by Anna
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin() ; // added by Anna
-
-
-
- //CsI Energy
- G4int CsIStageEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisCluster/CsIStageEnergy") ;
- CsIStageEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIStageEnergyCollectionID)) ;
- CsIStageEnergy_itr = CsIStageEnergyHitMap->GetMap()->begin() ;
-
-
- // Check the size of different map
- G4int sizeN = DetectorNumberHitMap->entries(); // number of hits by trackID=1 (can be the same object hit several time)
- G4int sizeC = CrystalNumberHitMap->entries(); // number of hits in all crystal by trackID=1
- G4int sizeE = EnergyHitMap->entries(); // = number of steps with edep non null
- G4int sizeT = TimeHitMap->entries();
- //G4int sizeX = PosXHitMap->entries();
- //G4int sizeY = PosYHitMap->entries();
- //G4int sizeX = XHitMap->entries();
- //G4int sizeY = YHitMap->entries();
-
- G4int sizeNCsI= CsIDetectorNumberHitMap->entries();
- G4int sizeCCsI= CsICrystalNumberHitMap->entries(); // added by Anna
- G4int sizeECsI= CsIStageEnergyHitMap->entries();
-
- sizeC *= 1; // remove warning at compilation added by Anna
- sizeCCsI *= 1; // remove warning at compilation added by Anna
- sizeECsI *= 1; // remove warning at compilation added by Anna
- /*
- G4cout <<"SizeN=" << sizeN << G4endl;
- G4cout <<"SizeC=" << sizeC << G4endl;
- G4cout <<"SizeE=" << sizeE << G4endl;
- G4cout <<"SizeT=" << sizeT << G4endl;
- G4cout <<"SizeN CsI =" << sizeNCsI << G4endl;
- G4cout <<"SizeC CsI =" << sizeCCsI << G4endl;
- G4cout <<"SizeE CsI =" << sizeECsI << G4endl;
- */
-
- //DetectorNumberHitMap->PrintAllHits();
-
-
- if (sizeE != sizeT) {
-// G4cout << "No match size PARIS Event Map: sE:"
-// << sizeE << " sT:" << sizeT << G4endl ;
-
- // if (sizeE != sizeX) {
- //G4cout << "No match size PARIS Event Map: sE:"
- //<< sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << G4endl ;
- return;
- }
-
-
- //G4cout <<"SizeN=" << sizeN << G4endl;
-
- //if at least 1 cluster is hit:
- if(sizeN>0)
- {
-
- // Deal with trackID=1:
- G4int N_first= *(DetectorNumber_itr->second); // ID of first hit det (det=cluster)
- G4int NTrackID = DetectorNumber_itr->first - N_first; // first trackID dealt with (not always =1)
- G4double E = *(Energy_itr->second);
- G4double T = *(Time_itr->second);
- G4int NCryst= *(CrystalNumber_itr->second); // ID of first hit crystal
-
- NCryst *= 1; //remove warning at compilation // added by Anna
- NTrackID *= 1; //remove warning at compilation
- /*
- G4cout <<"NTrackID=" << NTrackID << G4endl;
- G4cout <<"N_first=" << N_first << G4endl;
- G4cout <<"CrystalNumber_first=" << NCryst << G4endl;
- G4cout <<"Energy first=" << E << G4endl;
- G4cout <<"Time first =" << T << G4endl;
- G4cout<<"*******"<<endl;
-
- //added by Nicolas on the model of gaspard scorers 8/7/11
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin();
-
- for (G4int l = 0 ; l < sizeC ; l++) {
- G4double NCryst = *(CrystalNumber_itr->second);
- G4int CTrackID = CrystalNumber_itr->first - N_first - NCryst;
- G4cout << N_first << " " << NTrackID << " " << CTrackID << " " << NCryst << G4endl;
- if (CTrackID == NTrackID) {
- G4cout << "****************** NCryst " << NCryst << G4endl;
- //ms_Event->SetGPDTrkFirstStageFrontEEnergy(RandGauss::shoot(E, ResoFirstStage));
- //ms_Event->SetGPDTrkFirstStageBackEEnergy(RandGauss::shoot(E, ResoFirstStage));
- }
- CrystalNumber_itr++;
- }
- //rewind:
- for (G4int l = 0 ; l < sizeC ; l++) {
- CrystalNumber_itr--;
- }
- */
-
- //if there is more than 1 interaction in crystals (1 interaction = 1 hit)
- if(sizeN>1)
- {
-
- //At clusters'level
- for (G4int l = 1; l < sizeN ; l++) { // loop on all the other hits
-
- Energy_itr++;
- Time_itr++;
- DetectorNumber_itr++;
- CrystalNumber_itr++;
-
- G4int N= *(DetectorNumber_itr->second); // ID of hit cluster
- G4int Nc= *(CrystalNumber_itr->second); // ID of hit crystal
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << DetectorNumber_itr->first << G4endl;
- //G4cout <<"TrackID=" << NTrackID << G4endl;
- //G4cout <<"CTrackID=" << CTrackID << G4endl;
-
- // At clusters'level:
- if(N==N_first) // if the hit is still in the same cluster
- {
- // At crystals'level:
- if(Nc==NCryst) // if we are still in the same crystal than the first hit crystal
- {
- E += *(Energy_itr->second); // if same detector then sum up the energy
-
- // G4cout <<"Energy=" << E << G4endl;
-
- }else // we fill the tree for the first hit crystal and move to the next hit crystal if any
- {
- if(E!=0)
- {
- //G4cout <<" -1- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Cluster"] + N_first );
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Cluster"] + N_first );
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Anna
- ms_Event->SetPARISLaBr3StageTCrystalNbr(NCryst); // added by Anna
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree with energy deposited in first detector
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
-
- NCryst=Nc; // continue with the next hit crystal
- E=*(Energy_itr->second);
- }
-
-
- }else // if it's a new cluster, one fills the tree with first hit and inspect this new cluster
- {
-
- if(E!=0)
- {
-// G4cout <<" -2- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Cluster"] + N_first );
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Cluster"] + N_first );
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Anna
- ms_Event->SetPARISLaBr3StageTCrystalNbr(NCryst); // added by Anna
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree with energy deposited in first detector
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
- N_first=N; // continue with the new hit cluster
- NCryst=Nc; // and continue with the new hit crystal in this new cluster
- E=*(Energy_itr->second);
-
- }
-
- //G4cout <<"Energy=" << E << G4endl;
- //G4cout <<"Time =" << T << G4endl;
-
-
- // Always fill the tree at the end of the loop:
- if(l==(sizeN-1) && E!=0)
- {
-// G4cout <<"-3- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Cluster"] + N_first);
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Cluster"] + N_first);
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Anna
- ms_Event->SetPARISLaBr3StageTCrystalNbr(NCryst); // added by Anna
-// G4cout<<NTrackID<<" filled at the end "<<NCryst<<endl;
-
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
-
-
- } // end of the l loop
-
- }else
- {
- // Fill the tree if sizeN=1 and E>0
- if(E!=0)
- {
-// G4cout <<" -4- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Cluster"] + N_first);
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Cluster"] + N_first);
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst);
- ms_Event->SetPARISLaBr3StageTCrystalNbr(NCryst); // added by Anna
-// G4cout<<NTrackID<<" filled with sizeN=1 "<<NCryst<<endl;
-
-
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
- }
-
- }
-
-
- ///////// For CsI stage:
- //EGammaMin = EGamma-4*ResoSecondStage;
- //EGammaMax = EGamma+4*ResoSecondStage;
- if(sizeNCsI>0)
- {
- // Deal with trackID=1:
- G4int NCsI_first= *(CsIDetectorNumber_itr->second); // ID of first det hit
- G4int NCsITrackID = CsIDetectorNumber_itr->first - NCsI_first; // first trackID dealt with (not always =1)
- G4double E_CsI=*(CsIStageEnergy_itr->second); // Energy second stage
- G4int NCsICryst= *(CsICrystalNumber_itr->second); //9/6 Nicolas and Anna added
-
-
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
- //G4cout <<"NCsI_first=" << NCsI_first << G4endl;
- //G4cout <<"CsI Energy first=" << E_CsI << G4endl;
-
- //added by Anna and Nicolas on the model of gaspard scorers 8/7/11
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin();
- /*
- for (G4int l = 0 ; l < sizeCCsI ; l++) {
- G4double NCsICryst = *(CsICrystalNumber_itr->second);
- G4int CCsITrackID = CsICrystalNumber_itr->first - NCsI_first - NCsICryst;
- G4cout << NCsI_first << " " << NCsITrackID << " " << CCsITrackID << " " << NCsICryst << G4endl;
- if (CCsITrackID == NCsITrackID) {
- G4cout << "****************** NCsICryst " << NCsICryst << G4endl;
- //ms_Event->SetGPDTrkFirstStageFrontEEnergy(RandGauss::shoot(E, ResoFirstStage));
- //ms_Event->SetGPDTrkFirstStageBackEEnergy(RandGauss::shoot(E, ResoFirstStage));
- }
- CsICrystalNumber_itr++;
- }
- //rewind:
- for (G4int l = 0 ; l < sizeCCsI ; l++) {
- CsICrystalNumber_itr--;
- }
- */
-
- if(sizeNCsI>1)
- {
- for (G4int l = 1; l < sizeNCsI ; l++) { // loop on all the other tracks
-
- CsIStageEnergy_itr++;
- CsIDetectorNumber_itr++;
- CsICrystalNumber_itr++;
-
- G4int NCsI= *(CsIDetectorNumber_itr->second); // ID of hit cluster
- //NCsITrackID = CsIDetectorNumber_itr->first - NCsI; // ID of the track
- G4int NcCsI= *(CsICrystalNumber_itr->second); // ID of hit CsI crystal
- //NcCsITrackID = CsIDetectorNumber_itr->first - NCsI; // ID of the track
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"NCsI=" << NCsI << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << CsIDetectorNumber_itr->first << G4endl;
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
-
-
- // At clusters'level:
- if(NCsI==NCsI_first) // if the hit is still in the same cluster
- {
-
- // At crystals'level:
- if(NcCsI==NCsICryst) // if we are still in the same crystal than the first hit CsI crystal
- {
- E_CsI += *(CsIStageEnergy_itr->second); // if same detector then sum up the energy
-
- }else // we fill the tree for the first detector hit and move to the next detector hit
- {
- if(E_CsI!=0)
- {
-// G4cout <<" -1- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- ms_Event->SetPARISCsIStageTCrystalNbr(NCsICryst);
- // Fill Energy
- // ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- //ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst); // added by Anna
- //cout<<NCsITrackID<<" filled with sizeN=1 "<<NCsICryst<<endl;
- }
-
- NCsICryst=NcCsI; // continue with the next hit CsI crystal
- E_CsI=*(CsIStageEnergy_itr->second);
- }
-
- }else // if it's a new cluster, one fills the tree with first hit and inspect the new cluster
- {
- if(E_CsI!=0)
- {
-// G4cout <<" -2- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Cluster"] + NCsI_first );
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Cluster"] + NCsI_first );
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst); // added by Anna
- ms_Event->SetPARISCsIStageTCrystalNbr(NCsICryst); // added by Anna
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoFirstStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree with energy deposited in first detector
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- }
-
- NCsI_first=NCsI; // continue with the new hit cluster
- NCsICryst=NcCsI; // and continue with the new hit CsI crystal in this new cluster
- E_CsI=*(CsIStageEnergy_itr->second);
-
- }
-
-
- // Always fill the tree at the end of the loop:
-
- if(l==(sizeNCsI-1) && E_CsI!=0)
- {
-// G4cout <<" -3- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- ms_Event->SetPARISCsIStageTCrystalNbr(NCsICryst);
- // Fill Energy
- // ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst); // added by Anna
- G4cout<<NCsITrackID<<" filled with sizeN=1 "<<NCsICryst<<endl;
- }
-
- } // end of l loop
-
- }else
- {
- // Fill the tree if sizeN=1:
- if(E_CsI!=0)
- {
-// G4cout <<" -4- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Cluster"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- ms_Event->SetPARISCsIStageTCrystalNbr(NCsICryst);
- // Fill Energy
- // ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- // Fill Crystal number // added by Anna
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- //cout<<NCsITrackID<<" filled with sizeN=1 "<<NCsICryst<<endl;
-
- }
- }
-
- }
-
-
-
-
- // clear map for next event
- DetectorNumberHitMap -> clear();
- CrystalNumberHitMap -> clear(); // added by anna
- EnergyHitMap -> clear();
- TimeHitMap -> clear();
- //XHitMap -> clear();
- //YHitMap -> clear();
- //PosXHitMap -> clear();
- //PosYHitMap -> clear();
- //PosZHitMap -> clear();
- //AngThetaHitMap -> clear();
- //AngPhiHitMap -> clear();
-
- CsIDetectorNumberHitMap -> clear();
- CsICrystalNumberHitMap -> clear(); // added by Anna
- CsIStageEnergyHitMap -> clear();
- // SecondStageEnergyHitMap -> clear();
- // ThirdStageEnergyHitMap -> clear();
-
-
-}
-
-
-
-void ParisCluster::InitializeScorers()
-{
- bool already_exist = false;
-
- // LaBr3 Associate Scorer
- m_LaBr3StageScorer = VDetector::CheckScorer("LaBr3StageScorerParisCluster",already_exist);
- m_CsIStageScorer = VDetector::CheckScorer("CsIStageScorerParisCluster",already_exist);
-
- if(already_exist) return;
-
- /**/
- // G4VPrimitiveScorer* DetNbr = new OBSOLETEGENERALSCORERS::PSDetectorNumber("DetectorNumber", "ParisCluster", 0);
- G4VPrimitiveScorer* DetNbr = new PARISScorerLaBr3StageDetectorNumber("DetectorNumber", "ParisCluster", 0);
- // G4VPrimitiveScorer* TOF = new OBSOLETEGENERALSCORERS::PSTOF("StripTime","ParisCluster", 0);
- G4VPrimitiveScorer* TOF = new PARISScorerLaBr3StageTOF("StripTime","ParisCluster", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesX = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","ParisCluster", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesY = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","ParisCluster", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesZ = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","ParisCluster", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","ParisCluster", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","ParisCluster", 0);
-
- G4VPrimitiveScorer* Energy = new PARISScorerLaBr3StageEnergy("StripEnergy", "ParisCluster", 0);
- G4VPrimitiveScorer* CrystNbr = new PARISScorerLaBr3StageCrystal("CrystalNumber", "ParisCluster", 0);
-
- // G4VPrimitiveScorer* StripPositionX = new PARIScorerLaBr3StageFrontStripDummyShape("StripIDFront", 0, NumberOfStrips);
- // G4VPrimitiveScorer* StripPositionY = new PARISScorerLaBr3StageBackStripDummyShape("StripIDBack", 0, NumberOfStrips);
-
- //and register it to the multifunctionnal detector
- m_LaBr3StageScorer->RegisterPrimitive(DetNbr);
- m_LaBr3StageScorer->RegisterPrimitive(CrystNbr);
- m_LaBr3StageScorer->RegisterPrimitive(Energy);
- m_LaBr3StageScorer->RegisterPrimitive(TOF);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionX);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesX);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesZ);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
-
-
-
- /**/
-
- // Second stage Associate Scorer
- /**/
- G4VPrimitiveScorer* CsIDetNbr = new PARISScorerCsIStageDetectorNumber("CsIDetectorNumber", "ParisCluster", 0);
- G4VPrimitiveScorer* CsICryNbr = new PARISScorerCsIStageCrystalNumber("CsICrystalNumber", "ParisCluster", 0);
- G4VPrimitiveScorer* CsIStageEnergy = new PARISScorerCsIStageEnergy("CsIStageEnergy", "ParisCluster", 0);
-
- m_CsIStageScorer->RegisterPrimitive(CsIDetNbr);
- m_CsIStageScorer->RegisterPrimitive(CsICryNbr); // Added by anna
- m_CsIStageScorer->RegisterPrimitive(CsIStageEnergy);
- /**/
-
- // Add All Scorer to the Global Scorer Manager
- G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3StageScorer);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIStageScorer);
-}
diff --git a/NPSimulation/Paris/ParisCluster.hh b/NPSimulation/Paris/ParisCluster.hh
deleted file mode 100644
index 6dcb2bddf..000000000
--- a/NPSimulation/Paris/ParisCluster.hh
+++ /dev/null
@@ -1,169 +0,0 @@
-#ifndef ParisCluster_h
-#define ParisCluster_h 1
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: Define a cluster of 9 phoswich modules for the Paris detector *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-
-// C++ headers
-#include <vector>
-
-// NPTool header
-#include "ParisModule.hh"
-#include "TInteractionCoordinates.h"
-
-using namespace std;
-#include "CLHEP/Units/SystemOfUnits.h"
-using namespace CLHEP;
-
-
-
-class ParisCluster : public ParisModule
-{
- ////////////////////////////////////////////////////
- /////// Default Constructor and Destructor /////////
- ////////////////////////////////////////////////////
-public:
- ParisCluster();
- virtual ~ParisCluster();
-
- ////////////////////////////////////////////////////
- //////// Specific Function of this Class ///////////
- ////////////////////////////////////////////////////
-public:
- // By Position Method
- void AddModule(G4ThreeVector TL,
- G4ThreeVector BL,
- G4ThreeVector BR,
- G4ThreeVector CT);
-
- // By Angle Method
- void AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w);
-
- // Effectively construct Volume
- // Avoid to have two time same code for Angle and Point definition
- void VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world);
-
-
- ///////////////////////////////////////////
- //// Inherite from ParisModule class /////
- ///////////////////////////////////////////
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
-
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world);
-
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput();
-
- // Initialize all scorers necessary for the detector
- void InitializeScorers();
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
-
- // Give the static TInteractionCoordinates from VDetector to the classes
- // deriving from ParisModule
- // This is mandatory since the Paris*** does not derive from VDetector
- void SetInterCoordPointer(TInteractionCoordinates* interCoord);
- TInteractionCoordinates* GetInterCoordPointer() {return ms_InterCoord;};
-
-
- ////////////////////////////////////////////////////
- ///////////////Private intern Data//////////////////
- ////////////////////////////////////////////////////
-private:
- // Interaction Coordinates coming from VDetector through the
- // SetInteractionCoordinatesPointer method
- TInteractionCoordinates* ms_InterCoord;
-
- // True if Define by Position, False is Define by angle
- vector<bool> m_DefinitionType ;
-
- // Used for "By Point Definition"
- vector<G4ThreeVector> m_X1_Y1 ; // Top Left Corner Position Vector
- vector<G4ThreeVector> m_X1_Y128 ; // Bottom Left Corner Position Vector
- vector<G4ThreeVector> m_X128_Y1 ; // Bottom Right Corner Position Vector
- vector<G4ThreeVector> m_X128_Y128 ; // Center Corner Position Vector
-
- // Used for "By Angle Definition"
- vector<G4double> m_R ; // |
- vector<G4double> m_Theta ; // > Spherical coordinate of Strips Silicium Plate
- vector<G4double> m_Phi ; // |
-
- vector<G4double> m_beta_u ; // |
- vector<G4double> m_beta_v ; // > Tilt angle of the Telescope
- vector<G4double> m_beta_w ; // |
-
- G4ThreeVector momentum;
-
-};
-
-
-
-namespace PARISCLUSTER
-{
- // Resolution
-// const G4double ResoFirstStage = 0; // = 50 keV of Resolution // Unit is MeV/2.35
- const G4double ResoFirstStage = 0.0085; // = 3% at .662 MeV of Resolution // Unit is MeV/2.35
- const G4double ResoSecondStage = 0.0366; // = 13% at .662 MeV of resolution // Unit is MeV/2.35
- const G4double ResoThirdStage = 0.0; // = 50 keV of resolution // Unit is MeV/2.35
- const G4double ResoTimeGpd = 0.212765957;// = 500ps // Unit is ns/2.35
-
- // Geometry for the mother volume containing the different layers of your dummy shape module
- const G4double FaceFront = 16.9*cm;
- const G4double FaceBack = 16.9*cm;
- //const G4double Length = 21*cm; // for 2x2x2 LaBr + 2x2x6 NaI
- const G4double Length = 23.6*cm; // for 2x2x3 LaBr + 2x2x6 NaI
-
- // Geometry for the phoswitch
- // LaBr3
- const G4double LaBr3Face = 5.08*cm;
- //const G4double LaBr3Thickness = 5.08*cm; // for phoswich 2x2x2
- const G4double LaBr3Thickness = 7.62*cm; // for phoswich 2x2x3
- //const G4double LaBr3Thickness = 15.24*cm; // for long LaBr3 2x2x6
- const G4double InterSpace = 0.3*cm;
-
- // CsI
- const G4double CsIFace = LaBr3Face;
- const G4double CsIThickness = 15.24*cm; // for phoswich 2x2x6
- //const G4double CsIThickness = 0.001*cm; // for long LaBr3
-
- // Starting form the LaBr3 and going to the CsI
- const G4double LaBr3Stage_PosZ = Length* -0.5 + 0.5*LaBr3Thickness;
- const G4double CsIStage_PosZ = Length* -0.5 + LaBr3Thickness + 0.5*CsIThickness;
-}
-
-#endif
diff --git a/NPSimulation/Paris/ParisModule.cc b/NPSimulation/Paris/ParisModule.cc
deleted file mode 100644
index 84ea44e65..000000000
--- a/NPSimulation/Paris/ParisModule.cc
+++ /dev/null
@@ -1,64 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class is an Abstract Base Class (ABC) from which should *
- * derive all different modules from the Paris detector. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-#include "ParisModule.hh"
-#include "RootOutput.h"
-
-
-TParisData *ParisModule::ms_Event = 0;
-
-
-
-ParisModule::ParisModule()
-{
- if (ms_Event == 0) ms_Event = new TParisData();
-
- InitializeRootOutput();
- InitializeIndex();
-}
-
-
-
-ParisModule::~ParisModule()
-{
-}
-
-
-
-void ParisModule::InitializeRootOutput()
-{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- // if the branch does not exist yet, create it
- if (!pTree->GetBranch("PARIS"))
- pTree->Branch("PARIS", "TParisData", &ms_Event);
-
- pTree->SetBranchAddress("PARIS", &ms_Event);
-
-}
-
-
-
-void ParisModule::InitializeIndex()
-{
- m_index["Phoswich"] = 0; // 24 phoswish max
- m_index["Cluster"] = 100; // 18 cluster max
-}
diff --git a/NPSimulation/Paris/ParisModule.hh b/NPSimulation/Paris/ParisModule.hh
deleted file mode 100644
index 89524b999..000000000
--- a/NPSimulation/Paris/ParisModule.hh
+++ /dev/null
@@ -1,97 +0,0 @@
-#ifndef ParisModule_h
-#define ParisModule_h 1
-
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class is an Abstract Base Class (ABC) from which should *
- * derive all different modules from the Paris detector. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-// C++ headers
-#include <string>
-#include <vector>
-
-// G4 headers
-#include "G4LogicalVolume.hh"
-#include "G4Event.hh"
-#include "G4MultiFunctionalDetector.hh"
-
-// NPTool - ROOT headers
-#include "TInteractionCoordinates.h"
-#include "TParisData.h"
-
-using namespace std;
-#include"CLHEP/Units/SystemOfUnits.h"
-using namespace CLHEP;
-
-
-
-class ParisModule
-{
-public:
- ParisModule();
- virtual ~ParisModule();
-
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- virtual void ReadConfiguration(string Path) = 0;
-
- // Construct detector and inialise sensitive part.
- virtual void ConstructDetector(G4LogicalVolume* world) = 0;
-
- // Read sensitive part of a the ParisModule detector and fill the Root tree.
- virtual void ReadSensitive(const G4Event* event) = 0;
-
- // Add Detector branch to the ROOT output tree
- virtual void InitializeRootOutput();
-
- // Initialize all scorers necessary for each detector
- virtual void InitializeScorers() = 0;
-
- // Give the static TInteractionCoordinates from VDetector to the classes
- // deriving from ParisModule
- // This is mandatory since the Paris*** does not derive from VDetector
- virtual void SetInterCoordPointer(TInteractionCoordinates* interCoord) = 0;
- virtual TInteractionCoordinates* GetInterCoordPointer() = 0;
-
- // Initialize the Index map for the different modules of Paris
- void InitializeIndex();
-
-public:
- TParisData* GetEventPointer() {return ms_Event;};
-
-protected:
- // Class to store the result of simulation
- static TParisData* ms_Event;
-
- // Set to true if you want to see Telescope Frame in your visualisation
- bool m_non_sensitive_part_visiualisation;
-
-protected:
- // First stage Associate Scorer
- G4MultiFunctionalDetector* m_LaBr3StageScorer;
-
- // Second stage Associate Scorer
- G4MultiFunctionalDetector* m_CsIStageScorer;
-
-protected:
- map<string, int> m_index;
-};
-
-#endif
diff --git a/NPSimulation/Paris/ParisPhoswich.cc b/NPSimulation/Paris/ParisPhoswich.cc
deleted file mode 100644
index bb43ca833..000000000
--- a/NPSimulation/Paris/ParisPhoswich.cc
+++ /dev/null
@@ -1,1068 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : 20/07/11 *
- *---------------------------------------------------------------------------*
- * Decription: Define a phoswich module for the Paris detector. *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-// C++ headers
-#include <sstream>
-#include <string>
-#include <cmath>
-
-// G4 Geometry headers
-#include "G4Trd.hh"
-#include "G4Box.hh"
-#include "G4Trap.hh"
-
-// G4 various headers
-#include "G4Material.hh"
-#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
-#include "G4RotationMatrix.hh"
-#include "G4Transform3D.hh"
-#include "G4PVPlacement.hh"
-#include "G4PVDivision.hh"
-
-// G4 sensitive
-#include "G4SDManager.hh"
-#include "G4MultiFunctionalDetector.hh"
-
-// NPTool headers
-#include "ParisPhoswich.hh"
-#include "ObsoleteGeneralScorers.hh"
-#include "ParisScorers.hh"
-#include "MaterialManager.hh"
-#include "RootOutput.h"
-#include "VDetector.hh"
-// CLHEP
-#include "CLHEP/Random/RandGauss.h"
-
-using namespace std;
-using namespace CLHEP;
-using namespace PARISPHOSWICH;
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-ParisPhoswich::ParisPhoswich()
-{
- ms_InterCoord = 0;
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-ParisPhoswich::~ParisPhoswich()
-{
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisPhoswich::AddModule(G4ThreeVector X1_Y1,
- G4ThreeVector X128_Y1,
- G4ThreeVector X1_Y128,
- G4ThreeVector X128_Y128)
-{
- m_DefinitionType.push_back(true);
-
- m_X1_Y1.push_back(X1_Y1);
- m_X128_Y1.push_back(X128_Y1);
- m_X1_Y128.push_back(X1_Y128);
- m_X128_Y128.push_back(X128_Y128);
-
- m_R.push_back(0);
- m_Theta.push_back(0);
- m_Phi.push_back(0);
- m_beta_u.push_back(0);
- m_beta_v.push_back(0);
- m_beta_w.push_back(0);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisPhoswich::AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w)
-{
- G4ThreeVector empty = G4ThreeVector(0, 0, 0);
-
- m_DefinitionType.push_back(false);
-
- m_X1_Y1.push_back(empty);
- m_X128_Y1.push_back(empty);
- m_X1_Y128.push_back(empty);
- m_X128_Y128.push_back(empty);
-
- m_R.push_back(R);
- m_Theta.push_back(Theta);
- m_Phi.push_back(Phi);
- m_beta_u.push_back(beta_u);
- m_beta_v.push_back(beta_v);
- m_beta_w.push_back(beta_w);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParisPhoswich::VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world)
-{
- G4double NbrTelescopes = DetecNumber;
- G4String DetectorNumber;
- ostringstream Number;
- Number << NbrTelescopes;
- DetectorNumber = Number.str();
-
- G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
- G4Material* NaI = MaterialManager::getInstance()->GetMaterialFromLibrary("NaI");
- G4Material* LaBr3 = MaterialManager::getInstance()->GetMaterialFromLibrary("LaBr3");
-
- ////////////////////////////////////////////////////////////////
- ////////////// Starting Volume Definition //////////////////////
- ////////////////////////////////////////////////////////////////
- G4String Name = "ParisPhoswich" + DetectorNumber ;
-
- // Mother Volume
- // G4Box* solidParisPhoswich = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
- G4Box* solidParisPhoswich = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length);
- G4LogicalVolume* logicParisPhoswich = new G4LogicalVolume(solidParisPhoswich, Vacuum, Name, 0, 0, 0);
-
- new G4PVPlacement(G4Transform3D(*MMrot, MMpos) ,
- logicParisPhoswich ,
- Name ,
- world ,
- false ,
- 0);
-
- logicParisPhoswich->SetVisAttributes(G4VisAttributes::Invisible);
- if (m_non_sensitive_part_visiualisation) logicParisPhoswich->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
-
- // Phoswich construction
- // LaBr3
- G4ThreeVector positionLaBr3Stage = G4ThreeVector(0, 0, LaBr3Stage_PosZ);
-
- G4Box* solidLaBr3Stage = new G4Box("solidLaBr3Stage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*LaBr3Thickness);
- G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, LaBr3, "logicLaBr3Stage", 0, 0, 0);
- //G4LogicalVolume* logicLaBr3Stage = new G4LogicalVolume(solidLaBr3Stage, Ge, "logicLaBr3Stage", 0, 0, 0);
-
- new G4PVPlacement(0,
- positionLaBr3Stage,
- logicLaBr3Stage,
- Name + "_LaBr3Stage",
- logicParisPhoswich,
- false,
- 0);
-
- // Set LaBr3 sensible
- logicLaBr3Stage->SetSensitiveDetector(m_LaBr3StageScorer);
-
- // Visualisation of LaBr3Stage Strip
- G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));
- logicLaBr3Stage->SetVisAttributes(LaBr3VisAtt);
-
- // CsI or NaI
- G4ThreeVector positionCsIStage = G4ThreeVector(0, 0, CsIStage_PosZ);
-
- G4Box* solidCsIStage = new G4Box("solidCsIStage", 0.5*CsIFace, 0.5*CsIFace, 0.5*CsIThickness);
- //G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, CsI, "logicCsIStage", 0, 0, 0);
- G4LogicalVolume* logicCsIStage = new G4LogicalVolume(solidCsIStage, NaI, "logicCsIStage", 0, 0, 0);
-
- new G4PVPlacement(0,
- positionCsIStage,
- logicCsIStage,
- Name + "_CsIStage",
- logicParisPhoswich,
- false,
- 0);
-
- // Set CsI sensible
- logicCsIStage->SetSensitiveDetector(m_CsIStageScorer);
-
- // Visualisation of CsIStage Strip
- G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1., 0., 0.));
- logicCsIStage->SetVisAttributes(CsIVisAtt);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Virtual Method of VDetector class
-
-// Read stream at Configfile to pick-up parameters of detector (Position,...)
-// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void ParisPhoswich::ReadConfiguration(string Path)
-{
- ifstream ConfigFile;
- ConfigFile.open(Path.c_str());
- string LineBuffer, DataBuffer;
-
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
- G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
- G4ThreeVector A , B , C , D ;
- G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
-
- bool ReadingStatus = false;
-
- bool check_A = false;
- bool check_C = false;
- bool check_B = false;
- bool check_D = false;
-
- bool check_Theta = false;
- bool check_Phi = false;
- bool check_R = false;
-
- bool checkVis = false;
-
- while (!ConfigFile.eof()) {
- getline(ConfigFile, LineBuffer);
- if (LineBuffer.compare(0, 14, "ParisPhoswich") == 0) {
- G4cout << "///" << G4endl ;
- G4cout << "Phoswich element found: " << G4endl ;
- ReadingStatus = true ;
- }
-
- while (ReadingStatus) {
- ConfigFile >> DataBuffer;
- // Comment Line
- if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;}
-
- // Position method
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- Ax = atof(DataBuffer.c_str()) ;
- Ax = Ax * mm ;
- ConfigFile >> DataBuffer ;
- Ay = atof(DataBuffer.c_str()) ;
- Ay = Ay * mm ;
- ConfigFile >> DataBuffer ;
- Az = atof(DataBuffer.c_str()) ;
- Az = Az * mm ;
-
- A = G4ThreeVector(Ax, Ay, Az);
- G4cout << "X1 Y1 corner position : " << A << G4endl;
- }
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- Bx = atof(DataBuffer.c_str()) ;
- Bx = Bx * mm ;
- ConfigFile >> DataBuffer ;
- By = atof(DataBuffer.c_str()) ;
- By = By * mm ;
- ConfigFile >> DataBuffer ;
- Bz = atof(DataBuffer.c_str()) ;
- Bz = Bz * mm ;
-
- B = G4ThreeVector(Bx, By, Bz);
- G4cout << "X128 Y1 corner position : " << B << G4endl;
- }
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- Cx = atof(DataBuffer.c_str()) ;
- Cx = Cx * mm ;
- ConfigFile >> DataBuffer ;
- Cy = atof(DataBuffer.c_str()) ;
- Cy = Cy * mm ;
- ConfigFile >> DataBuffer ;
- Cz = atof(DataBuffer.c_str()) ;
- Cz = Cz * mm ;
-
- C = G4ThreeVector(Cx, Cy, Cz);
- G4cout << "X1 Y128 corner position : " << C << G4endl;
- }
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- Dx = atof(DataBuffer.c_str()) ;
- Dx = Dx * mm ;
- ConfigFile >> DataBuffer ;
- Dy = atof(DataBuffer.c_str()) ;
- Dy = Dy * mm ;
- ConfigFile >> DataBuffer ;
- Dz = atof(DataBuffer.c_str()) ;
- Dz = Dz * mm ;
-
- D = G4ThreeVector(Dx, Dy, Dz);
- G4cout << "X128 Y128 corner position : " << D << G4endl;
- }
-
- // Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- Theta = Theta * deg;
- G4cout << "Theta: " << Theta / deg << G4endl;
- }
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- Phi = Phi * deg;
- G4cout << "Phi: " << Phi / deg << G4endl;
- }
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
- R = R * mm;
- G4cout << "R: " << R / mm << G4endl;
- }
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- beta_u = beta_u * deg ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- beta_v = beta_v * deg ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- beta_w = beta_w * deg ;
- G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
- }
-
- else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
- checkVis = true ;
- ConfigFile >> DataBuffer;
- if (DataBuffer.compare(0, 3, "all") == 0) m_non_sensitive_part_visiualisation = true;
- }
-
- else G4cout << "WARNING: Wrong Token, ParisPhoswich: Phoswich Element not added" << G4endl;
-
- // Add The previously define telescope
- // With position method
- if ((check_A && check_B && check_C && check_D && checkVis) &&
- !(check_Theta && check_Phi && check_R)) {
- ReadingStatus = false;
- check_A = false;
- check_C = false;
- check_B = false;
- check_D = false;
- checkVis = false;
-
- AddModule(A, B, C, D);
- }
-
- // With angle method
- if ((check_Theta && check_Phi && check_R && checkVis) &&
- !(check_A && check_B && check_C && check_D)) {
- ReadingStatus = false;
- check_Theta = false;
- check_Phi = false;
- check_R = false;
- checkVis = false;
-
- AddModule(R, Theta, Phi, beta_u, beta_v, beta_w);
- }
- }
- }
-}
-
-
-
-// Construct detector and inialise sensitive part.
-// Called After DetecorConstruction::AddDetector Method
-void ParisPhoswich::ConstructDetector(G4LogicalVolume* world)
-{
- G4RotationMatrix* MMrot = NULL ;
- G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
-
- G4int NumberOfTelescope = m_DefinitionType.size() ;
-
- for (G4int i = 0; i < NumberOfTelescope; i++) {
- // By Point
- if (m_DefinitionType[i]) {
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- MMu = m_X128_Y1[i] - m_X1_Y1[i];
- MMu = MMu.unit();
-
- MMv = m_X1_Y128[i] - m_X1_Y1[i];
- MMv = MMv.unit();
-
- MMw = MMu.cross(MMv);
-// if (MMw.z() > 0) MMw = MMv.cross(MMu) ;
- MMw = MMw.unit();
-
- MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4;
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + MMCenter;
- }
-
- // By Angle
- else {
- G4double Theta = m_Theta[i] ;
- G4double Phi = m_Phi[i] ;
-
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing ThirdStage
- // Phi is angle between X axis and projection in (X,Y) plan
- // Theta is angle between position vector and z axis
- G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
- G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
- G4double wZ = m_R[i] * cos(Theta / rad);
- MMw = G4ThreeVector(wX, wY, wZ);
-
- // vector corresponding to the center of the module
- G4ThreeVector CT = MMw;
-
- // vector parallel to one axis of silicon plane
- G4double ii = cos(Theta / rad) * cos(Phi / rad);
- G4double jj = cos(Theta / rad) * sin(Phi / rad);
- G4double kk = -sin(Theta / rad);
- G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
-
- MMw = MMw.unit();
- MMu = MMw.cross(Y);
- MMv = MMw.cross(MMu);
- MMv = MMv.unit();
- MMu = MMu.unit();
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // Telescope is rotate of Beta angle around MMv axis.
- MMrot->rotate(m_beta_u[i], MMu);
- MMrot->rotate(m_beta_v[i], MMv);
- MMrot->rotate(m_beta_w[i], MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + CT ;
- }
-
- VolumeMaker(i + 1, MMpos, MMrot, world);
- }
-
- delete MMrot ;
-}
-
-
-
-// Connect the ParisData class to the output TTree
-// of the simulation
-void ParisPhoswich::InitializeRootOutput()
-{
-}
-
-
-
-// Set the TinteractionCoordinates object from VDetector to the present class
-void ParisPhoswich::SetInterCoordPointer(TInteractionCoordinates* interCoord)
-{
- ms_InterCoord = interCoord;
-}
-
-
-
-// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAvtion
-void ParisPhoswich::ReadSensitive(const G4Event* event)
-{
- //////////////////////////////////////////////////////////////////////////////////////
- //////////////////////// Used to Read Event Map of detector //////////////////////////
- //////////////////////////////////////////////////////////////////////////////////////
-
- momentum = event->GetPrimaryVertex()->GetPrimary()->GetMomentum();
- G4double EGamma = momentum.getR(); // for photon E=p
- EGamma *= 1;
- //G4double EGammaMin = EGamma-4*ResoFirstStage;
- //G4double EGammaMax = EGamma+4*ResoFirstStage;
-
-
- // First Stage (LaBr3)
- std::map<G4int, G4int*>::iterator DetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CrystalNumber_itr;
- std::map<G4int, G4double*>::iterator Energy_itr;
- std::map<G4int, G4double*>::iterator Time_itr;
- //std::map<G4int, G4double*>::iterator X_itr;
- //std::map<G4int, G4double*>::iterator Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_X_itr;
- //std::map<G4int, G4double*>::iterator Pos_Y_itr;
- //std::map<G4int, G4double*>::iterator Pos_Z_itr;
- //std::map<G4int, G4double*>::iterator Ang_Theta_itr;
- //std::map<G4int, G4double*>::iterator Ang_Phi_itr;
-
- G4THitsMap<G4int>* DetectorNumberHitMap;
- G4THitsMap<G4int>* CrystalNumberHitMap;
- G4THitsMap<G4double>* EnergyHitMap;
- G4THitsMap<G4double>* TimeHitMap;
- // G4THitsMap<G4double>* XHitMap;
- // G4THitsMap<G4double>* YHitMap;
- //G4THitsMap<G4double>* PosXHitMap;
- //G4THitsMap<G4double>* PosYHitMap;
- //G4THitsMap<G4double>* PosZHitMap;
- //G4THitsMap<G4double>* AngThetaHitMap;
- //G4THitsMap<G4double>* AngPhiHitMap;
-
- // NULL pointer are given to avoid warning at compilation
- // Second Stage (CsI)
- std::map<G4int, G4int*>::iterator CsIDetectorNumber_itr;
- std::map<G4int, G4int*>::iterator CsICrystalNumber_itr; // added by Marc
- std::map<G4int, G4double*>::iterator CsIStageEnergy_itr ;
- G4THitsMap<G4int>* CsIDetectorNumberHitMap = NULL ;
- G4THitsMap<G4int>* CsICrystalNumberHitMap = NULL ; // added by Marc
- G4THitsMap<G4double>* CsIStageEnergyHitMap = NULL ;
-
-
- // Read the Scorer associate to the LaBr
- //Detector Number
- G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/DetectorNumber") ;
- DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ;
- DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ;
-
- //Crystal Number
- G4int CrystDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/CrystalNumber") ;
- CrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CrystDetCollectionID)) ;
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin() ;
-
- //Energy
- G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/StripEnergy") ;
- EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
- Energy_itr = EnergyHitMap->GetMap()->begin() ;
-
- //Time of Flight
- G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/StripTime") ;
- TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ;
- Time_itr = TimeHitMap->GetMap()->begin() ;
-
- /*
- //Strip Number X
- G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/StripIDFront") ;
- XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ;
- X_itr = XHitMap->GetMap()->begin() ;
-
- //Strip Number Y
- G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/StripIDBack");
- YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ;
- Y_itr = YHitMap->GetMap()->begin() ;
-
-
- //Interaction Coordinate X
- G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/InterCoordX") ;
- PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ;
- Pos_X_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Y
- G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/InterCoordY") ;
- PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ;
- Pos_Y_itr = PosYHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Z
- G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/InterCoordZ") ;
- PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ;
- Pos_Z_itr = PosXHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Theta
- G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/InterCoordAngTheta") ;
- AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ;
- Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ;
-
- //Interaction Coordinate Angle Phi
- G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LaBr3StageScorerParisPhoswich/InterCoordAngPhi") ;
- AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ;
- Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ;
- */
-
-
- // Read the Scorer associate to the SecondStage
- //CsI Detector Number
- G4int CsIDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisPhoswich/CsIDetectorNumber") ;
- CsIDetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsIDetCollectionID)) ;
- CsIDetectorNumber_itr = CsIDetectorNumberHitMap->GetMap()->begin() ;
-
- //CsI Crystal Number // added by Marc
- G4int CsICryCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisPhoswich/CsICrystalNumber") ; // added by Marc
- CsICrystalNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICryCollectionID)) ; // added by Anna
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin() ; // added by Anna
-
- //CsI Energy
- G4int CsIStageEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIStageScorerParisPhoswich/CsIStageEnergy") ;
- CsIStageEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIStageEnergyCollectionID)) ;
- CsIStageEnergy_itr = CsIStageEnergyHitMap->GetMap()->begin() ;
-
-
- // Check the size of different map
- G4int sizeN = DetectorNumberHitMap->entries(); // number of objects hit by trackID=1 (can be the same object hit several time)
- G4int sizeC = CrystalNumberHitMap->entries();
- G4int sizeE = EnergyHitMap->entries(); // = number of steps with edep non null
- G4int sizeT = TimeHitMap->entries();
- //G4int sizeX = PosXHitMap->entries();
- //G4int sizeY = PosYHitMap->entries();
- //G4int sizeX = XHitMap->entries();
- //G4int sizeY = YHitMap->entries();
-
- G4int sizeNCsI= CsIDetectorNumberHitMap->entries();
- G4int sizeCCsI= CsICrystalNumberHitMap->entries();
- G4int sizeECsI= CsIStageEnergyHitMap->entries();
-
- sizeC *= 1; // remove warning at compilation added by Marc
- sizeCCsI *= 1; // remove warning at compilation added by Marc
- sizeECsI *= 1; // remove warning at compilation added by Marc
- //G4cout <<"SizeN=" << sizeN << G4endl;
- //G4cout <<"SizeC=" << sizeC << G4endl;
- //G4cout <<"SizeN CsI =" << sizeNCsI << G4endl;
- //G4cout <<"SizeE CsI =" << sizeECsI << G4endl;
-
- //DetectorNumberHitMap->PrintAllHits();
-
-
- if (sizeE != sizeT) {
- G4cout << "No match size PARIS Event Map: sE:"
- << sizeE << " sT:" << sizeT << G4endl ;
-
- // if (sizeE != sizeX) {
- //G4cout << "No match size PARIS Event Map: sE:"
- //<< sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << G4endl ;
- return;
- }
-
-
- //G4cout <<"SizeN=" << sizeN << G4endl;
-
-
- //if at least 1 phoswich is hit:
- if(sizeN>0)
- {
-
- // Deal with trackID=1:
- G4int N_first= *(DetectorNumber_itr->second); // ID of first det hit
- G4int NTrackID = DetectorNumber_itr->first - N_first; // first trackID dealt with (not always =1)
- G4double E = *(Energy_itr->second);
- G4double T = *(Time_itr->second);
- G4int NCryst= *(CrystalNumber_itr->second);
-
- NCryst *= 1; //remove warning at compilation // added by Marc
- /*
- G4cout <<"NTrackID=" << NTrackID << G4endl;
- G4cout <<"N_first=" << N_first << G4endl;
- G4cout <<"CrystalNumber_first=" << NCryst << G4endl;
- G4cout <<"Energy first=" << E << G4endl;
- G4cout <<"Time first =" << T << G4endl;
- G4cout<<"*******"<<endl;
-
- //added by Nicolas on the model of gaspard scorers 8/7/11
- CrystalNumber_itr = CrystalNumberHitMap->GetMap()->begin();
-
- for (G4int l = 0 ; l < sizeC ; l++) {
- G4double NCryst = *(CrystalNumber_itr->second);
- G4int CTrackID = CrystalNumber_itr->first - N_first - NCryst;
- G4cout << N_first << " " << NTrackID << " " << CTrackID << " " << NCryst << G4endl;
- if (CTrackID == NTrackID) {
- G4cout << "****************** NCryst " << NCryst << G4endl;
- //ms_Event->SetGPDTrkFirstStageFrontEEnergy(RandGauss::shoot(E, ResoFirstStage));
- //ms_Event->SetGPDTrkFirstStageBackEEnergy(RandGauss::shoot(E, ResoFirstStage));
- }
- CrystalNumber_itr++;
- }
- //rewind:
- for (G4int l = 0 ; l < sizeC ; l++) {
- CrystalNumber_itr--;
- }
- */
-
- //if there is more than 1 interaction in crystals (1 interaction = 1 hit)
- if(sizeN>1)
- {
-
- //At clusters'level
- for (G4int l = 1; l < sizeN ; l++) { // loop on all the other tracks
-
- Energy_itr++;
- Time_itr++;
- DetectorNumber_itr++;
- CrystalNumber_itr++;
-
- G4int N= *(DetectorNumber_itr->second); // ID of Phoswich hit
- G4int Nc= *(CrystalNumber_itr->second); // ID of hit crystal (always =0 for Phoswich)
-
- NTrackID = DetectorNumber_itr->first - N; // ID of the track
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"N=" << N << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << DetectorNumber_itr->first << G4endl;
- //G4cout <<"NTrackID=" << NTrackID << G4endl;
-
- if(N==N_first)
- {
- // At crystals'level:
- if(Nc==NCryst) // if we are still in the same crystal than the first hit crystal
- {
- E += *(Energy_itr->second);
-
- }else // we fill the tree for the first detector hit and move to the next detector hit
- {
- if(E!=0)
- {
- G4cout <<" -1- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Marc
- // Fill Energy
- // ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
- NCryst=Nc; // continue with the next hit crystal
- E=*(Energy_itr->second);
- }
-
-
- }else // if it's a new phoswich, one fills the tree with first hit and inspect this new cluster
- {
-
- if(E!=0)
- {
- G4cout <<" -2- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Phoswich"] + N_first );
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Phoswich"] + N_first );
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Anna
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree with energy deposited in first detector
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
-
- }
-
- N_first=N;
- NCryst=Nc; // and continue with the new hit crystal in this new cluster
- E=*(Energy_itr->second);
-
- }
-
-
- //G4cout <<"Energy=" << E << G4endl;
- //G4cout <<"Time =" << T << G4endl;
-
- // Always fill the tree at the end of the loop:
- if(l==(sizeN-1) && E!=0)
- {
- G4cout <<"-3- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst); // added by Anna
- G4cout<<NTrackID<<" filled at the end "<<NCryst<<endl;
- // Fill Energy
- // ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
-
- }
-
- }else
- {
- // Fill the tree if sizeN=1:
- if(E!=0)
- {
- G4cout <<" -4- filling tree with energy=" << E << G4endl;
- // Fill detector number
- ms_Event->SetPARISLaBr3StageEDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageTDetectorNbr(m_index["Phoswich"] + N_first);
- ms_Event->SetPARISLaBr3StageECrystalNbr(NCryst);
- G4cout<<NTrackID<<" filled with sizeN=1 "<<NCryst<<endl;
- // Fill Energy
- // ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E=RandGauss::shoot(E, ResoFirstStage);
- ms_Event->SetPARISLaBr3StageEEnergy(E); // Fill the tree
- //if(E>EGammaMin && E<EGammaMax)ms_Event->SetPARISLaBr3StageEffphpeak(EGamma);
-
- // Fill Time
- ms_Event->SetPARISLaBr3StageTTime(RandGauss::shoot(T, ResoTimeGpd));
- }
- }
-
-
-
- }
-
-
- ///////// For CsI stage:
- //EGammaMin = EGamma-4*ResoSecondStage;
- //EGammaMax = EGamma+4*ResoSecondStage;
- if(sizeNCsI>0)
- {
- // Deal with trackID=1:
- G4int NCsI_first= *(CsIDetectorNumber_itr->second); // ID of first det hit
- G4int NCsITrackID = CsIDetectorNumber_itr->first - NCsI_first; // first trackID dealt with (not always =1)
- NCsITrackID *= 1;
- G4double E_CsI=*(CsIStageEnergy_itr->second); // Energy second stage
- G4int NCsICryst= *(CsICrystalNumber_itr->second); // added by Marc
-
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
- //G4cout <<"NCsI_first=" << NCsI_first << G4endl;
- //G4cout <<"CsI Energy first=" << E_CsI << G4endl;
-
- CsICrystalNumber_itr = CsICrystalNumberHitMap->GetMap()->begin(); // added by Marc
- /*
- for (G4int l = 0 ; l < sizeCCsI ; l++) {
- G4double NCsICryst = *(CsICrystalNumber_itr->second);
- G4int CCsITrackID = CsICrystalNumber_itr->first - NCsI_first - NCsICryst;
- G4cout << NCsI_first << " " << NCsITrackID << " " << CCsITrackID << " " << NCsICryst << G4endl;
- if (CCsITrackID == NCsITrackID) {
- G4cout << "****************** NCsICryst " << NCsICryst << G4endl;
- }
- CsICrystalNumber_itr++;
- }
- //rewind:
- for (G4int l = 0 ; l < sizeCCsI ; l++) {
- CsICrystalNumber_itr--;
- }
- */
-
- if(sizeNCsI>1)
- {
- for (G4int l = 1; l < sizeNCsI ; l++) { // loop on all the other tracks
-
- CsIStageEnergy_itr++;
- CsIDetectorNumber_itr++;
- CsICrystalNumber_itr++;
-
-
- G4int NCsI= *(CsIDetectorNumber_itr->second); // ID of Phoswich hit
- G4int NcCsI= *(CsICrystalNumber_itr->second); // ID of crystal hit (always =0 for phoswich)
- NCsITrackID = CsIDetectorNumber_itr->first - NCsI; // ID of the track
-
- //G4cout <<"l=" << l << G4endl;
- //G4cout <<"NCsI=" << NCsI << G4endl;
- //G4cout <<"DetectorNumber_itr->first =" << CsIDetectorNumber_itr->first << G4endl;
- //G4cout <<"NCsITrackID=" << NCsITrackID << G4endl;
-
- // At phoswichs'level:
- if(NCsI==NCsI_first)
- {
- // At crystals'level:
- if(NcCsI==NCsICryst) // if we are still in the same crystal than the first hit CsI crystal
- {
-
- E_CsI += *(CsIStageEnergy_itr->second);
-
- }else // we fill the tree for the first detector hit and move to the next detector hit
- {
- if(E_CsI!=0)
- {
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- //ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
- }
-
- NCsICryst=NcCsI;
- E_CsI=*(CsIStageEnergy_itr->second);
- }
-
- }else // if it's a new cluster, one fills the tree with first hit and inspect the new cluster
- {
- if(E_CsI!=0)
- {
- G4cout <<" -2- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Phoswich"] + NCsI_first );
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Phoswich"] + NCsI_first );
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst); // added by Marc
- // Fill Energy
- //ms_Event->SetPARISLaBr3StageEEnergy(RandGauss::shoot(E, ResoFirstStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoFirstStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree with energy deposited in first detector
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
-
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- }
-
- NCsI_first=NCsI; // continue with the new hit cluster
- NCsICryst=NcCsI; // and continue with the new hit CsI crystal in this new cluster
- E_CsI=*(CsIStageEnergy_itr->second);
-
- }
-
- // Always fill the tree at the end of the loop:
-
- if(l==(sizeNCsI-1) && E_CsI!=0)
- {
- G4cout <<" -3- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- // ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
-
- }
-
- }
-
- }else
- {
- // Fill the tree if sizeN=1:
- if(E_CsI!=0)
- {
- G4cout <<" -4- filling tree with energy in CsI =" << E_CsI << G4endl;
- // Fill detector number
- ms_Event->SetPARISCsIStageEDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageTDetectorNbr(m_index["Phoswich"] + NCsI_first);
- ms_Event->SetPARISCsIStageECrystalNbr(NCsICryst);
- // Fill Energy
- //ms_Event->SetPARISCsIStageEEnergy(RandGauss::shoot(E_CsI, ResoSecondStage));
- E_CsI=RandGauss::shoot(E_CsI, ResoSecondStage);
- ms_Event->SetPARISCsIStageEEnergy(E_CsI); // Fill the tree
- //if(E_CsI>EGammaMin && E_CsI<EGammaMax)ms_Event->SetPARISCsIStageEffphpeak(EGamma);
- // Fill Time
- //ms_Event->SetPARISCsIStageTTime(RandGauss::shoot(T_CsI, ResoTimeGpd));
- }
- }
-
- }
-
-
-
-
- // clear map for next event
- DetectorNumberHitMap -> clear();
- CrystalNumberHitMap -> clear(); // added by Marc
- EnergyHitMap -> clear();
- TimeHitMap -> clear();
- //XHitMap -> clear();
- //YHitMap -> clear();
- //PosXHitMap -> clear();
- //PosYHitMap -> clear();
- //PosZHitMap -> clear();
- //AngThetaHitMap -> clear();
- //AngPhiHitMap -> clear();
-
- CsIDetectorNumberHitMap -> clear();
- CsICrystalNumberHitMap -> clear(); // added by Marc
- CsIStageEnergyHitMap -> clear();
-
-
-
-}
-
-
-
-void ParisPhoswich::InitializeScorers()
-{
- bool already_exist = false;
-
- // LaBr3 Associate Scorer
- m_LaBr3StageScorer = VDetector::CheckScorer("LaBr3StageScorerParisCluster",already_exist);
- m_CsIStageScorer = VDetector::CheckScorer("CsIStageScorerParisCluster",already_exist);
-
- // G4VPrimitiveScorer* DetNbr = new OBSOLETEGENERALSCORERS::PSDetectorNumber("DetectorNumber", "ParisPhoswich", 0);
- G4VPrimitiveScorer* DetNbr = new PARISScorerLaBr3StageDetectorNumber("DetectorNumber", "ParisPhoswich", 0);
- // G4VPrimitiveScorer* TOF = new OBSOLETEGENERALSCORERS::PSTOF("StripTime","ParisPhoswich", 0);
- G4VPrimitiveScorer* TOF = new PARISScorerLaBr3StageTOF("StripTime","ParisPhoswich", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesX = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","ParisPhoswich", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesY = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","ParisPhoswich", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesZ = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","ParisPhoswich", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","ParisPhoswich", 0);
- //G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","ParisPhoswich", 0);
-
- G4VPrimitiveScorer* Energy = new PARISScorerLaBr3StageEnergy("StripEnergy", "ParisPhoswich", 0);
- G4VPrimitiveScorer* CrystNbr = new PARISScorerLaBr3StageCrystal("CrystalNumber", "ParisPhoswich", 0);
-
- // G4VPrimitiveScorer* StripPositionX = new PARIScorerLaBr3StageFrontStripDummyShape("StripIDFront", 0, NumberOfStrips);
- // G4VPrimitiveScorer* StripPositionY = new PARISScorerLaBr3StageBackStripDummyShape("StripIDBack", 0, NumberOfStrips);
-
- //and register it to the multifunctionnal detector
- m_LaBr3StageScorer->RegisterPrimitive(DetNbr);
- m_LaBr3StageScorer->RegisterPrimitive(CrystNbr);
- m_LaBr3StageScorer->RegisterPrimitive(Energy);
- m_LaBr3StageScorer->RegisterPrimitive(TOF);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionX);
- //m_LaBr3StageScorer->RegisterPrimitive(StripPositionY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesX);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesY);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesZ);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
- //m_LaBr3StageScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
-
-
-
- /**/
-
- // Second stage Associate Scorer
- /**/
- G4VPrimitiveScorer* CsIDetNbr = new PARISScorerCsIStageDetectorNumber("CsIDetectorNumber", "ParisPhoswich", 0);
- G4VPrimitiveScorer* CsICryNbr = new PARISScorerCsIStageCrystalNumber("CsICrystalNumber", "ParisPhoswich", 0); // added by Marc
- G4VPrimitiveScorer* CsIStageEnergy = new PARISScorerCsIStageEnergy("CsIStageEnergy", "ParisPhoswich", 0);
-
- m_CsIStageScorer->RegisterPrimitive(CsIDetNbr);
- m_CsIStageScorer->RegisterPrimitive(CsICryNbr); // Added by Marc
- m_CsIStageScorer->RegisterPrimitive(CsIStageEnergy);
- /**/
-
- // Add All Scorer to the Global Scorer Manager
- G4SDManager::GetSDMpointer()->AddNewDetector(m_LaBr3StageScorer);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIStageScorer);
-
-}
diff --git a/NPSimulation/Paris/ParisPhoswich.hh b/NPSimulation/Paris/ParisPhoswich.hh
deleted file mode 100644
index 1d6008141..000000000
--- a/NPSimulation/Paris/ParisPhoswich.hh
+++ /dev/null
@@ -1,173 +0,0 @@
-#ifndef ParisPhoswich_h
-#define ParisPhoswich_h 1
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 04/12/09 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: Define a phoswich module for the Paris detector. *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-
-// C++ headers
-#include <vector>
-
-// NPTool header
-#include "ParisModule.hh"
-#include "TInteractionCoordinates.h"
-
-using namespace std;
-
using namespace CLHEP;
-
-
-
-class ParisPhoswich : public ParisModule
-{
- ////////////////////////////////////////////////////
- /////// Default Constructor and Destructor /////////
- ////////////////////////////////////////////////////
-public:
- ParisPhoswich();
- virtual ~ParisPhoswich();
-
- ////////////////////////////////////////////////////
- //////// Specific Function of this Class ///////////
- ////////////////////////////////////////////////////
-public:
- // By Position Method
- void AddModule(G4ThreeVector TL,
- G4ThreeVector BL,
- G4ThreeVector BR,
- G4ThreeVector CT);
-
- // By Angle Method
- void AddModule(G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w);
-
- // Effectively construct Volume
- // Avoid to have two time same code for Angle and Point definition
- void VolumeMaker(G4int DetecNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- G4LogicalVolume* world);
-
-
- ///////////////////////////////////////////
- //// Inherite from ParisModule class /////
- ///////////////////////////////////////////
-public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
-
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world);
-
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput();
-
- // Initialize all scorers necessary for the detector
- void InitializeScorers();
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
-
- // Give the static TInteractionCoordinates from VDetector to the classes
- // deriving from ParisModule
- // This is mandatory since the Paris*** does not derive from VDetector
- void SetInterCoordPointer(TInteractionCoordinates* interCoord);
- TInteractionCoordinates* GetInterCoordPointer() {return ms_InterCoord;};
-
-
- ////////////////////////////////////////////////////
- ///////////////Private intern Data//////////////////
- ////////////////////////////////////////////////////
-private:
- // Interaction Coordinates coming from VDetector through the
- // SetInteractionCoordinatesPointer method
- TInteractionCoordinates* ms_InterCoord;
-
-
-
- // True if Define by Position, False is Define by angle
- vector<bool> m_DefinitionType ;
-
- // Used for "By Point Definition"
- vector<G4ThreeVector> m_X1_Y1 ; // Top Left Corner Position Vector
- vector<G4ThreeVector> m_X1_Y128 ; // Bottom Left Corner Position Vector
- vector<G4ThreeVector> m_X128_Y1 ; // Bottom Right Corner Position Vector
- vector<G4ThreeVector> m_X128_Y128 ; // Center Corner Position Vector
-
- // Used for "By Angle Definition"
- vector<G4double> m_R ; // |
- vector<G4double> m_Theta ; // > Spherical coordinate of Strips Silicium Plate
- vector<G4double> m_Phi ; // |
-
- vector<G4double> m_beta_u ; // |
- vector<G4double> m_beta_v ; // > Tilt angle of the Telescope
- vector<G4double> m_beta_w ; // |
-
- G4ThreeVector momentum;
-
-};
-
-
-
-namespace PARISPHOSWICH
-{
-
-
-
- // Resolution
-// const G4double ResoFirstStage = 0; // = 50 keV of Resolution // Unit is MeV/2.35
- const G4double ResoFirstStage = 0.0099; // = 3.5% at .662MeV of Resolution // Unit is MeV/2.35
- const G4double ResoSecondStage = 0.0366; // = 13% at .662 MeV of resolution // Unit is MeV/2.35
- const G4double ResoThirdStage = 0.0213; // = 50 keV of resolution // Unit is MeV/2.35
- const G4double ResoTimeGpd = 0.212765957;// = 500ps // Unit is ns/2.35
-
- // Geometry for the mother volume containing the different layers of your dummy shape module
- const G4double FaceFront = 5.7*cm; //6.3*cm;
- const G4double FaceBack = 5.7*cm; //6.3*cm;
- //const G4double Length = 21*cm; // for 2x2x2 LaBr + 2x2x6 NaI
- const G4double Length = 23.6*cm;// for 2x2x3 LaBr + 2x2x6 NaI
-
- // Geometry for the phoswitch
- // LaBr3
- const G4double LaBr3Face = 5.08*cm;
- //const G4double LaBr3Thickness = 5.08*cm; // for phoswich 2x2x2
- const G4double LaBr3Thickness = 7.62*cm; // for phoswich 2x2x3
- //const G4double LaBr3Thickness = 15.24*cm; // for LaBr3 long 2x2x6 inches
-
- // CsI
- const G4double CsIFace = LaBr3Face;
- const G4double CsIThickness = 15.24*cm; // for phoswich
- //const G4double CsIThickness = 0.001*cm; // For LaBr long
-
- // Starting form the LaBr3 and going to the CsI
- const G4double LaBr3Stage_PosZ = Length* -0.5 + 0.5*LaBr3Thickness;
- const G4double CsIStage_PosZ = Length* -0.5 + LaBr3Thickness + 0.5*CsIThickness;
-
-}
-
-#endif
diff --git a/NPSimulation/Paris/ParisScorers.cc b/NPSimulation/Paris/ParisScorers.cc
deleted file mode 100644
index 44a8459b9..000000000
--- a/NPSimulation/Paris/ParisScorers.cc
+++ /dev/null
@@ -1,516 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 30/02/10 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class holds all the scorers needed by the *
- * Paris*** objects. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-// G4 headers
-#include "G4UnitsTable.hh"
-
-// NPTool headers
-#include "ObsoleteGeneralScorers.hh"
-#include "ParisScorers.hh"
-
-#include "ParisCluster.hh"
-#include "ParisPhoswich.hh"
-
-using namespace PARISCLUSTER;
-using namespace PARISPHOSWICH;
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Added by Adrien MATTA:
-// Those Scorer use TrackID as map index. This way ones can rebuild energy deposit,
-// time of flight or position,... particle by particle for each event. Because standard
-// scorer provide by G4 don't work this way but using a global ID for each event you should
-// not use those scorer with some G4 provided ones or being very carefull doing so.
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage Detector Number Scorer (deal with multiple particle hit)
-PARISScorerLaBr3StageDetectorNumber::PARISScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerLaBr3StageDetectorNumber::~PARISScorerLaBr3StageDetectorNumber()
-{
-}
-
-G4bool PARISScorerLaBr3StageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
- return TRUE;
-}
-void PARISScorerLaBr3StageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerLaBr3StageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerLaBr3StageDetectorNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerLaBr3StageDetectorNumber::DrawAll()
-{
-}
-
-void PARISScorerLaBr3StageDetectorNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
- << G4endl;
- }
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage Energy Scorer (deal with multiple particle hit)
-PARISScorerLaBr3StageEnergy::PARISScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerLaBr3StageEnergy::~PARISScorerLaBr3StageEnergy()
-{
-}
-
-G4bool PARISScorerLaBr3StageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
- POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, edep);
- return TRUE;
-}
-
-void PARISScorerLaBr3StageEnergy::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerLaBr3StageEnergy::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerLaBr3StageEnergy::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerLaBr3StageEnergy::DrawAll()
-{
-}
-
-void PARISScorerLaBr3StageEnergy::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
-PARISScorerLaBr3StageCrystal::PARISScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerLaBr3StageCrystal::~PARISScorerLaBr3StageCrystal()
-{
-}
-
-G4bool PARISScorerLaBr3StageCrystal::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
-
- //G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1);
- //Adde by Anna:
- G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); // daughter crystal volume
- //depth:0 is the world, 1 is the mother...
- //AC instead of GetReplicaNumber(1)
- //cout<<"****replica****"<<endl;
- //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
- //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- //EvtMap->add(DetNbr + index, CrystalNbr); // Marc
- //cout<<"index "<<index<<endl;
- //cout<<"DetNbr + index "<<DetNbr + index<<endl;
- EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);// modified by Anna
- return TRUE;
-
-}
-
-void PARISScorerLaBr3StageCrystal::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerLaBr3StageCrystal::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerLaBr3StageCrystal::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerLaBr3StageCrystal::DrawAll()
-{
-}
-
-void PARISScorerLaBr3StageCrystal::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
- << G4endl;
- }
-}
-
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// FirstStage ToF Scorer (deal with multiple particle hit)
-PARISScorerLaBr3StageTOF::PARISScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerLaBr3StageTOF::~PARISScorerLaBr3StageTOF()
-{
-}
-
-G4bool PARISScorerLaBr3StageTOF::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get TOF
- G4double TOF = aStep->GetPreStepPoint()->GetGlobalTime();
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, TOF);
- return TRUE;
-}
-
-void PARISScorerLaBr3StageTOF::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerLaBr3StageTOF::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerLaBr3StageTOF::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerLaBr3StageTOF::DrawAll()
-{
-}
-
-void PARISScorerLaBr3StageTOF::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
-}
-
-
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// SecondStage (CsI) Detector Number Scorer (deal with multiple particle hit)
-PARISScorerCsIStageDetectorNumber::PARISScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerCsIStageDetectorNumber::~PARISScorerCsIStageDetectorNumber()
-{
-}
-
-G4bool PARISScorerCsIStageDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- // int DetNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0); // this line does exactly the same than the line above
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4double edep = aStep->GetTotalEnergyDeposit();
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->set(CrystalNbr + DetNbr + index, DetNbr);
- return TRUE;
-}
-void PARISScorerCsIStageDetectorNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerCsIStageDetectorNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerCsIStageDetectorNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerCsIStageDetectorNumber::DrawAll()
-{
-}
-
-void PARISScorerCsIStageDetectorNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " Detector Number : " << G4BestUnit(*(itr->second), "DetectorNumber")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// CsIStage Energy Scorer (deal with multiple particle hit)
-PARISScorerCsIStageEnergy::PARISScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerCsIStageEnergy::~PARISScorerCsIStageEnergy()
-{
-}
-
-G4bool PARISScorerCsIStageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
- int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
-
- // get energy
- G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
- POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- //if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- EvtMap->add(CrystalNbr + DetNbr + index, edep);
- return TRUE;
-}
-
-void PARISScorerCsIStageEnergy::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerCsIStageEnergy::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerCsIStageEnergy::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerCsIStageEnergy::DrawAll()
-{
-}
-
-void PARISScorerCsIStageEnergy::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " energy deposit: " << G4BestUnit(*(itr->second), "Energy")
- << G4endl;
- }
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// all added by Anna:
-
-// FirstStage CrystalNbr Scorer (deal with multiple particle hit)
-PARISScorerCsIStageCrystalNumber::PARISScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth)
- : G4VPrimitiveScorer(name, depth), HCID(-1)
-{
- m_VolumeName = volumeName;
-}
-
-PARISScorerCsIStageCrystalNumber::~PARISScorerCsIStageCrystalNumber()
-{
-}
-
-G4bool PARISScorerCsIStageCrystalNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
-{
- // get detector number
- int DetNbr = OBSOLETEGENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName);
-
- // get energy
- G4int CrystalNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
- //depth:0 is the world, 1 is the mother...
- //AC instead of GetReplicaNumber(1)
- //cout<<"****replica****"<<endl;
- //cout<<"copy nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1)<<endl;
- //cout<<"replica nbr"<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetReplicaNumber(1)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2)<<endl;
- //cout<<aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0)<<endl;
-
- G4double edep = aStep->GetTotalEnergyDeposit();
- // if (edep < 100*keV) return FALSE;
- if (edep < 0*keV) return FALSE; // = HERE IS THE DIFFERENCE WITH GENERALSCORER
- G4int index = aStep->GetTrack()->GetTrackID();
- //cout<<"index "<<index<<endl;
- //cout<<"DetNbr + index "<<DetNbr + index<<endl;
- EvtMap->set(CrystalNbr + DetNbr + index, CrystalNbr);
- return TRUE;
-
-}
-
-void PARISScorerCsIStageCrystalNumber::Initialize(G4HCofThisEvent* HCE)
-{
- EvtMap = new G4THitsMap<G4int>(GetMultiFunctionalDetector()->GetName(), GetName());
- if (HCID < 0) {
- HCID = GetCollectionID(0);
- }
- HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
-}
-
-void PARISScorerCsIStageCrystalNumber::EndOfEvent(G4HCofThisEvent*)
-{
-}
-
-void PARISScorerCsIStageCrystalNumber::Clear()
-{
- EvtMap->clear();
-}
-
-void PARISScorerCsIStageCrystalNumber::DrawAll()
-{
-}
-
-void PARISScorerCsIStageCrystalNumber::PrintAll()
-{
- G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
- G4cout << " PrimitiveScorer " << GetName() << G4endl;
- G4cout << " Number of entries " << EvtMap->entries() << G4endl;
- std::map<G4int, G4int*>::iterator itr = EvtMap->GetMap()->begin();
- for (; itr != EvtMap->GetMap()->end(); itr++) {
- G4cout << " copy no.: " << itr->first
- << " crystal nbr: " << G4BestUnit(*(itr->second), "Crystal")
- << G4endl;
- }
-}
diff --git a/NPSimulation/Paris/ParisScorers.hh b/NPSimulation/Paris/ParisScorers.hh
deleted file mode 100644
index 1b5fa31f5..000000000
--- a/NPSimulation/Paris/ParisScorers.hh
+++ /dev/null
@@ -1,199 +0,0 @@
-#ifndef PARISScorer_h
-#define PARISScorer_h 1
-/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: M. Labiche contact address: marc.labiche@stfc.ac.uk *
- * *
- * Creation Date : 30/02/10 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: This class holds all the scorers needed by the *
- * GaspardTracker*** objects. *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-
-
-#include "G4VPrimitiveScorer.hh"
-#include "G4THitsMap.hh"
-
-namespace PARISSCORERS
-{
-
-
-class PARISScorerLaBr3StageDetectorNumber : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerLaBr3StageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerLaBr3StageDetectorNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-
-class PARISScorerLaBr3StageEnergy : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerLaBr3StageEnergy(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerLaBr3StageEnergy();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-
-
-class PARISScorerLaBr3StageCrystal : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerLaBr3StageCrystal(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerLaBr3StageCrystal();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-class PARISScorerLaBr3StageTOF : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerLaBr3StageTOF(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerLaBr3StageTOF();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-
-
-class PARISScorerCsIStageDetectorNumber : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerCsIStageDetectorNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerCsIStageDetectorNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-// Added by Anna
-class PARISScorerCsIStageCrystalNumber : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerCsIStageCrystalNumber(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerCsIStageCrystalNumber();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4int>* EvtMap;
-};
-
-
-class PARISScorerCsIStageEnergy : public G4VPrimitiveScorer
-{
-public: // with description
- PARISScorerCsIStageEnergy(G4String name, G4String volumeName, G4int depth = 0);
- virtual ~PARISScorerCsIStageEnergy();
-
-protected: // with description
- virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
-public:
- virtual void Initialize(G4HCofThisEvent*);
- virtual void EndOfEvent(G4HCofThisEvent*);
- virtual void Clear();
- virtual void DrawAll();
- virtual void PrintAll();
-
-private:
- G4String m_VolumeName;
- G4int HCID;
- G4THitsMap<G4double>* EvtMap;
-};
-
-}
-
-using namespace PARISSCORERS;
-#endif
diff --git a/NPSimulation/include/ResistiveStripScorers.hh b/NPSimulation/include/CalorimeterScorers.hh
similarity index 60%
rename from NPSimulation/include/ResistiveStripScorers.hh
rename to NPSimulation/include/CalorimeterScorers.hh
index dc05e1882..a08fd3a4e 100644
--- a/NPSimulation/include/ResistiveStripScorers.hh
+++ b/NPSimulation/include/CalorimeterScorers.hh
@@ -1,5 +1,5 @@
-#ifndef SharcScorers_h
-#define SharcScorers_h 1
+#ifndef CalorimeterScorers_h
+#define CalorimeterScorers_h 1
/*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
@@ -14,7 +14,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * File old the scorer specific to the Sharc Detector *
+ * File old the scorer specific to the Silicon Detector *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -31,47 +31,37 @@
using namespace std;
using namespace CLHEP;
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-namespace SILICONSCORERS{
- class PS_Silicon_Resistive : public G4VPrimitiveScorer{
-
- public: // with description
- PS_Silicon_Resistive(G4String name,
- G4double StripPlaneLength, G4double StripPlaneWidth,
- G4int NumberOfStripWidth,G4int depth=0);
-
- ~PS_Silicon_Resistive();
-
- protected: // with description
- G4bool ProcessHits(G4Step*, G4TouchableHistory*);
-
- public:
- void Initialize(G4HCofThisEvent*);
- void EndOfEvent(G4HCofThisEvent*);
- void clear();
- void DrawAll();
- void PrintAll();
-
- private: // Threshold
- G4double m_TriggerThreshold;
-
- private: // Geometry of the detector
- G4double m_StripPlaneLength;
- G4double m_StripPlaneWidth;
- G4int m_NumberOfStripWidth;
- G4double m_StripPitchWidth;
-
- private: // inherited from G4VPrimitiveScorer
- G4int HCID;
- G4THitsMap<G4double*>* EvtMap;
-
- private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
- G4ThreeVector m_Position ;
- G4int m_DetectorNumber ;
- G4int m_StripWidthNumber ;
- G4int m_Index ;
+namespace CALORIMETERSCORERS {
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ class PS_Calorimeter : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_Calorimeter(G4String name, vector<G4int> NestingLevel,G4int depth=0);
+ ~PS_Calorimeter();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+ private: // Threshold
+ G4double m_TriggerThreshold;
+
+ private: // How much level of volume nesting should be considered
+ // Give the list of the nesting level at which the copy number should be return.
+ // 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
+ vector<G4int> m_NestingLevel;
+ G4int m_Index;
+ private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ G4THitsMap<G4double*>* EvtMap;
};
-
}
+
#endif
diff --git a/NPSimulation/include/SiliconScorers.hh b/NPSimulation/include/SiliconScorers.hh
index 93eae93f1..af9e75869 100644
--- a/NPSimulation/include/SiliconScorers.hh
+++ b/NPSimulation/include/SiliconScorers.hh
@@ -32,8 +32,7 @@ using namespace std;
using namespace CLHEP;
namespace SILICONSCORERS {
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class PS_Silicon_Rectangle : public G4VPrimitiveScorer{
public: // with description
@@ -50,9 +49,6 @@ namespace SILICONSCORERS {
void DrawAll();
void PrintAll();
- private: // Threshold
- G4double m_TriggerThreshold;
-
private: // Geometry of the detector
G4double m_StripPlaneLength;
G4double m_StripPlaneWidth;
@@ -74,6 +70,7 @@ namespace SILICONSCORERS {
};
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class PS_Silicon_Annular : public G4VPrimitiveScorer{
public: // with description
@@ -90,9 +87,6 @@ namespace SILICONSCORERS {
void DrawAll();
void PrintAll();
- private: // Threshold
- G4double m_TriggerThreshold;
-
private: // Geometry of the detector
G4double m_StripPlaneInnerRadius;
G4double m_StripPlaneOuterRadius;
@@ -119,6 +113,44 @@ namespace SILICONSCORERS {
G4int m_Index ;
};
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ class PS_Silicon_Resistive : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_Silicon_Resistive(G4String name,
+ G4double StripPlaneLength, G4double StripPlaneWidth,
+ G4int NumberOfStripWidth,G4int depth=0);
+
+ ~PS_Silicon_Resistive();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+ private: // Geometry of the detector
+ G4double m_StripPlaneLength;
+ G4double m_StripPlaneWidth;
+ G4int m_NumberOfStripWidth;
+ G4double m_StripPitchWidth;
+
+ private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ G4THitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4ThreeVector m_Position ;
+ G4int m_DetectorNumber ;
+ G4int m_StripWidthNumber ;
+ G4int m_Index ;
+
+ };
+
}
#endif
diff --git a/NPSimulation/src/ResistiveStripScorers.cc b/NPSimulation/src/CalorimeterScorers.cc
similarity index 50%
rename from NPSimulation/src/ResistiveStripScorers.cc
rename to NPSimulation/src/CalorimeterScorers.cc
index 66ef34d7e..c0b6d87a0 100644
--- a/NPSimulation/src/ResistiveStripScorers.cc
+++ b/NPSimulation/src/CalorimeterScorers.cc
@@ -1,4 +1,4 @@
- /****************************************************************************
+ /*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
@@ -8,94 +8,66 @@
/*****************************************************************************
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
- * Creation Date : December 2013 *
+ * Creation Date : February 2013 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * Generic scorer for Resistive Strip Silicon detector *
+ * File old the scorer specific to the Sharc Detector *
* *
*---------------------------------------------------------------------------*
* Comment: *
- * Energy Sharing and Balistic deficit effect. *
- * *
+ * This new type of scorer is aim to become the standard for DSSD,SSSD and *
+ * PAD detector (any Silicon Detector) *
*****************************************************************************/
-#include "ResistiveStripScorers.hh"
+#include "CalorimeterScorers.hh"
#include "G4UnitsTable.hh"
-using namespace SILICONSCORERS;
+using namespace CALORIMETERSCORERS ;
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_Silicon_Resistive::PS_Silicon_Resistive(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripWidth,G4int depth)
+PS_Calorimeter::PS_Calorimeter(G4String name, vector<G4int> NestingLevel,G4int depth)
:G4VPrimitiveScorer(name, depth),HCID(-1){
- m_StripPlaneLength = StripPlaneLength;
- m_StripPlaneWidth = StripPlaneWidth;
- m_NumberOfStripWidth = NumberOfStripWidth;
- m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
-
- m_TriggerThreshold = 1*eV;
- m_Position = G4ThreeVector(-1000,-1000,-1000);
- m_DetectorNumber = -1;
- m_StripWidthNumber = -1;
- m_Index = -1 ;
+ m_NestingLevel = NestingLevel;
}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_Silicon_Resistive::~PS_Silicon_Resistive(){
+PS_Calorimeter::~PS_Calorimeter(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4bool PS_Silicon_Resistive::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+G4bool PS_Calorimeter::ProcessHits(G4Step* aStep, G4TouchableHistory*){
- // contain Energy Total, E1, E2, Time, DetNbr, and StripWidth
- G4double* EnergyAndTime = new G4double[10];
-
- EnergyAndTime[2] = aStep->GetPreStepPoint()->GetGlobalTime();
-
- m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
- 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_StripWidthNumber = (int)((m_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ;
- m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ;
-
- // The energy is divided in two depending on the position
- // position along the resistive strip
- double P = (m_Position.x())/(0.5*m_StripPlaneLength);
+ // contain Energy, Time + as many copy number as nested volume
+ unsigned int mysize = m_NestingLevel.size();
+ G4double* Infos = new G4double[2+mysize];
+ Infos[0] = aStep->GetTotalEnergyDeposit();
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
- // Upstream Energy
- EnergyAndTime[0] = aStep->GetTotalEnergyDeposit()*(1+P)*0.5;
-
- // Downstream Energy
- EnergyAndTime[1] = aStep->GetTotalEnergyDeposit()-EnergyAndTime[0];
+ for(G4int i = 0 ; i < mysize ; i++){
+ Infos[i+2] = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i]);
+ }
- EnergyAndTime[3] = m_DetectorNumber;
- EnergyAndTime[4] = m_StripWidthNumber;
-
- //Rare case where particle is close to edge of silicon plan
- if (m_StripWidthNumber > m_NumberOfStripWidth) m_StripWidthNumber = m_NumberOfStripWidth;
-
- m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripWidthNumber * 1e6;
+ m_Index = 0 ;
+ G4int multiplier = 1;
+ for(G4int i = 0 ; i < mysize ; i++){
+ m_Index+= aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i])*multiplier;
+ multiplier*=10;
+ }
// 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];
+ Infos[0]+=dummy[0];
}
- EvtMap->set(m_Index, EnergyAndTime);
+ EvtMap->set(m_Index, Infos);
return TRUE;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_Silicon_Resistive::Initialize(G4HCofThisEvent* HCE){
+void PS_Calorimeter::Initialize(G4HCofThisEvent* HCE){
EvtMap = new G4THitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
if (HCID < 0) {
HCID = GetCollectionID(0);
@@ -104,11 +76,11 @@ void PS_Silicon_Resistive::Initialize(G4HCofThisEvent* HCE){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_Silicon_Resistive::EndOfEvent(G4HCofThisEvent*){
+void PS_Calorimeter::EndOfEvent(G4HCofThisEvent*){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_Silicon_Resistive::clear(){
+void PS_Calorimeter::clear(){
std::map<G4int, G4double**>::iterator MapIterator;
for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
delete *(MapIterator->second);
@@ -118,12 +90,12 @@ void PS_Silicon_Resistive::clear(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_Silicon_Resistive::DrawAll(){
+void PS_Calorimeter::DrawAll(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_Silicon_Resistive::PrintAll(){
+void PS_Calorimeter::PrintAll(){
G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
G4cout << " PrimitiveScorer " << GetName() << G4endl ;
G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
diff --git a/NPSimulation/src/SiliconScorers.cc b/NPSimulation/src/SiliconScorers.cc
index bb08e7766..504861f0c 100644
--- a/NPSimulation/src/SiliconScorers.cc
+++ b/NPSimulation/src/SiliconScorers.cc
@@ -33,7 +33,6 @@ PS_Silicon_Rectangle::PS_Silicon_Rectangle(G4String name, G4double StripPlaneLen
m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength;
m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
- m_TriggerThreshold=1*eV;
m_Position = G4ThreeVector(-1000,-1000,-1000);
m_DetectorNumber = -1;
m_StripLengthNumber = -1;
@@ -142,7 +141,6 @@ PS_Silicon_Annular::PS_Silicon_Annular(G4String name, G4double StripPlaneInnerRa
m_StripPitchSector = (m_PhiStop-m_PhiStart)/m_NumberOfStripSector;
m_StripPitchQuadrant = (m_PhiStop-m_PhiStart)/m_NumberOfStripQuadrant;
- m_TriggerThreshold = 1*eV;
m_Position = G4ThreeVector(-1000,-1000,-1000);
m_uz = G4ThreeVector(0,0,1);
m_DetectorNumber = -1;
@@ -239,3 +237,110 @@ void PS_Silicon_Annular::PrintAll(){
G4cout << " PrimitiveScorer " << GetName() << G4endl ;
G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Resistive::PS_Silicon_Resistive(G4String name, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripWidth,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_StripPlaneLength = StripPlaneLength;
+ m_StripPlaneWidth = StripPlaneWidth;
+ m_NumberOfStripWidth = NumberOfStripWidth;
+ m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
+
+ m_Position = G4ThreeVector(-1000,-1000,-1000);
+ m_DetectorNumber = -1;
+ m_StripWidthNumber = -1;
+ m_Index = -1 ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_Silicon_Resistive::~PS_Silicon_Resistive(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_Silicon_Resistive::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Total, E1, E2, Time, DetNbr, and StripWidth
+ G4double* EnergyAndTime = new G4double[10];
+
+ EnergyAndTime[2] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1);
+ 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_StripWidthNumber = (int)((m_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ;
+ m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ;
+
+ // The energy is divided in two depending on the position
+ // position along the resistive strip
+ double P = (m_Position.x())/(0.5*m_StripPlaneLength);
+
+ // Upstream Energy
+ EnergyAndTime[0] = aStep->GetTotalEnergyDeposit()*(1+P)*0.5;
+
+ // Downstream Energy
+ EnergyAndTime[1] = aStep->GetTotalEnergyDeposit()-EnergyAndTime[0];
+
+ EnergyAndTime[3] = m_DetectorNumber;
+ EnergyAndTime[4] = m_StripWidthNumber;
+
+ //Rare case where particle is close to edge of silicon plan
+ if (m_StripWidthNumber > m_NumberOfStripWidth) m_StripWidthNumber = m_NumberOfStripWidth;
+
+ m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripWidthNumber * 1e6;
+
+ // 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;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Resistive::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new G4THitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Resistive::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Resistive::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Resistive::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_Silicon_Resistive::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+
+
--
GitLab