From fe69bdaf77d7881e549b89ce983a701d8f944433 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?GIRARD=20ALCINDOR=20Val=C3=A9rian?=
<girard-alcindor@ijclab.in2p3.fr>
Date: Tue, 8 Oct 2024 10:23:06 +0200
Subject: [PATCH] Modifiying Exogam simu
---
NPSimulation/Detectors/Exogam/Exogam.cc | 148 ++++++++++++++++++++----
NPSimulation/Detectors/Exogam/Exogam.hh | 13 ++-
2 files changed, 135 insertions(+), 26 deletions(-)
diff --git a/NPSimulation/Detectors/Exogam/Exogam.cc b/NPSimulation/Detectors/Exogam/Exogam.cc
index 3c001d15a..209391827 100644
--- a/NPSimulation/Detectors/Exogam/Exogam.cc
+++ b/NPSimulation/Detectors/Exogam/Exogam.cc
@@ -70,7 +70,8 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Exogam_NS {
// Energy and time Resolution
- const double EnergyThreshold = 10 * keV;
+ // const double EnergyThreshold = 10 * keV;
+ const double EnergyThreshold = -100 * keV;
// const double ResoTime = 4.5*ns ; //not used
const double ResoEnergy = 2. * keV;
} // namespace Exogam_NS
@@ -79,7 +80,7 @@ namespace Exogam_NS {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Exogam Specific Method
Exogam::Exogam() {
- m_Event = new TExogamData();
+ m_Event = new TExogamCalData();
m_ExogamScorer = 0;
InitializeMaterials();
@@ -127,12 +128,65 @@ void Exogam::BuildClover(int i_clo, G4LogicalVolume* world) {
G4LogicalVolume* logicSupClover = new G4LogicalVolume(solidSupClover, m_Vacuum, "SupClover");
Offset = dzEnv; //-distCollimatorToGeCan;
+ // Offset = 0;
- G4RotationMatrix rm;
- rm.rotateX(m_ThetaX[i_clo] / rad).rotateY(m_ThetaY[i_clo] / rad).rotateZ(m_ThetaZ[i_clo] / rad);
-
- new G4PVPlacement(G4Transform3D(rm, G4ThreeVector(m_X[i_clo] * mm, m_Y[i_clo] * mm, m_Z[i_clo] * mm + Offset)),
- logicSupClover, "Clover", world, false, i_clo + 1, false); // this void overlaps the whole setup
+ if (m_X[i_clo] >= 0) {
+ G4RotationMatrix rm;
+ rm.rotateX(m_ThetaX[i_clo] / rad).rotateY(m_ThetaY[i_clo] / rad).rotateZ(m_ThetaZ[i_clo] / rad);
+ new G4PVPlacement(G4Transform3D(rm, G4ThreeVector(m_X[i_clo] * mm, m_Y[i_clo] * mm, m_Z[i_clo] * mm + Offset)),
+ logicSupClover, "Clover", world, false, i_clo + 1, false); // this void overlaps the whole setup
+ }
+ else if (m_R[i_clo] >= 0) {
+
+ 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);
+ G4double Theta = m_Theta[i_clo];
+ G4double Phi = m_Phi[i_clo];
+
+ // (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_clo]+Offset) * sin(Theta / rad) * cos(Phi / rad);
+ G4double wY = (m_R[i_clo]+Offset) * sin(Theta / rad) * sin(Phi / rad);
+ G4double wZ = (m_R[i_clo]+Offset) * 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_clo], MMu);
+ // MMrot->rotate(m_beta_v[i_clo], MMv);
+ // MMrot->rotate(m_beta_w[i_clo], MMw);
+ // translation to place Telescope
+ double Length = 0 * cm;
+ MMpos = MMw * Length * 0.5 + CT;
+ std::cout << i_clo << std::endl;
+ std::string clover_name = "Clover_" + std::to_string(i_clo + 1);
+ new G4PVPlacement(G4Transform3D(*MMrot, MMpos), logicSupClover, clover_name, world, false, i_clo + 1,
+ false); // this void overlaps the whole setup
+ }
// The Cryostat
////////////////
@@ -157,7 +211,7 @@ void Exogam::BuildClover(int i_clo, G4LogicalVolume* world) {
rOuterCan[1] = 6.2 * cm;
rOuterCan[2] = 6.2 * cm;
- rInnerCan[0] = rInnerCan[1] = rInnerCan[2] = 0. * cm;
+ rInnerCan[0] = rInnerCan[1] = rInnerCan[2] = 0.1 * cm;
G4Polyhedra* solidCloverCan =
new G4Polyhedra("CloverCan", PhiStartCan, PhiTotCan, 4, 3, zPlaneCan, rInnerCan, rOuterCan);
@@ -414,15 +468,24 @@ void Exogam::BuildClover(int i_clo, G4LogicalVolume* world) {
// Now the individual diode is built; create logical volumes for each of
// the four individual diodes A, B, C and D:
- G4LogicalVolume* logicGeA = new G4LogicalVolume(solidGe, m_Germanium, "GeA");
- G4LogicalVolume* logicGeB = new G4LogicalVolume(solidGe, m_Germanium, "GeB");
- G4LogicalVolume* logicGeC = new G4LogicalVolume(solidGe, m_Germanium, "GeC");
- G4LogicalVolume* logicGeD = new G4LogicalVolume(solidGe, m_Germanium, "GeD");
+ std::string GeA_name = "GeA_" + std::to_string(i_clo + 1);
+ G4LogicalVolume* logicGeA = new G4LogicalVolume(solidGe, m_Germanium, GeA_name);
+ std::string GeB_name = "GeB_" + std::to_string(i_clo + 1);
+ G4LogicalVolume* logicGeB = new G4LogicalVolume(solidGe, m_Germanium, GeB_name);
+ std::string GeC_name = "GeC_" + std::to_string(i_clo + 1);
+ G4LogicalVolume* logicGeC = new G4LogicalVolume(solidGe, m_Germanium, GeC_name);
+ std::string GeD_name = "GeD_" + std::to_string(i_clo + 1);
+ G4LogicalVolume* logicGeD = new G4LogicalVolume(solidGe, m_Germanium, GeD_name);
+ std::cout << "SetSensitiveDetector: " << i_clo << std::endl;
logicGeA->SetSensitiveDetector(m_ExogamScorer);
+ std::cout << "GeA SensitiveDetector: " << logicGeA->GetSensitiveDetector()->GetName() << std::endl;
logicGeB->SetSensitiveDetector(m_ExogamScorer);
+ std::cout << "GeB SensitiveDetector: " << logicGeB->GetSensitiveDetector()->GetName() << std::endl;
logicGeC->SetSensitiveDetector(m_ExogamScorer);
+ std::cout << "GeC SensitiveDetector: " << logicGeC->GetSensitiveDetector()->GetName() << std::endl;
logicGeD->SetSensitiveDetector(m_ExogamScorer);
+ std::cout << "GeD SensitiveDetector: " << logicGeD->GetSensitiveDetector()->GetName() << std::endl;
// positioning the tapered partial diodes (A to D)
// into the real vacuum of the can
@@ -443,10 +506,15 @@ void Exogam::BuildClover(int i_clo, G4LogicalVolume* world) {
G4ThreeVector posDumVacD(-xDumVac, yDumVac, zDumVac);
G4Transform3D positionVacD(rm90, posDumVacD);
- new G4PVPlacement(0, positionVacA, logicGeA, "GeA", logicVac, false, 1, true); // There is an overlap with vacumm Vac
- new G4PVPlacement(positionVacB, logicGeB, "GeB", logicVac, false, 2, true);
- new G4PVPlacement(positionVacC, logicGeC, "GeC", logicVac, false, 3, true);
- new G4PVPlacement(positionVacD, logicGeD, "GeD", logicVac, false, 4, true);
+ std::string GeA_name2 = "GeA_" + std::to_string(i_clo + 1);
+ new G4PVPlacement(0, positionVacA, logicGeA, GeA_name, logicVac, false, 1,
+ true); // There is an overlap with vacumm Vac
+ std::string GeB_name2 = "GeB_" + std::to_string(i_clo + 1);
+ new G4PVPlacement(positionVacB, logicGeB, GeB_name, logicVac, false, 2, true);
+ std::string GeC_name2 = "GeC_" + std::to_string(i_clo + 1);
+ new G4PVPlacement(positionVacC, logicGeC, GeC_name, logicVac, false, 3, true);
+ std::string GeD_name2 = "GeD_" + std::to_string(i_clo + 1);
+ new G4PVPlacement(positionVacD, logicGeD, GeD_name, logicVac, false, 4, true);
//
// some material between the diodes to reproduce the experimental addback factor ...
@@ -530,6 +598,24 @@ void Exogam::AddDetector(double X, double Y, double Z, double ThetaX, double The
m_ThetaX.push_back(ThetaX);
m_ThetaY.push_back(ThetaY);
m_ThetaZ.push_back(ThetaZ);
+
+ m_R.push_back(-1000);
+ m_Theta.push_back(-1000);
+ m_Phi.push_back(-1000);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Exogam::AddDetector(double R, double Theta, double Phi) {
+ m_R.push_back(R);
+ m_Theta.push_back(Theta);
+ m_Phi.push_back(Phi);
+
+ m_X.push_back(-1000);
+ m_Y.push_back(-1000);
+ m_Z.push_back(-1000);
+ m_ThetaX.push_back(-1000);
+ m_ThetaY.push_back(-1000);
+ m_ThetaZ.push_back(-1000);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -546,6 +632,7 @@ void Exogam::ReadConfiguration(NPL::InputParser parser) {
cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> coord = {"X", "Y", "Z", "ThetaX", "ThetaY", "ThetaZ"};
+ vector<string> sphe = {"R", "Theta", "Phi"};
for (unsigned int i = 0; i < blocks.size(); i++) {
if (blocks[i]->HasTokenList(coord)) {
@@ -559,6 +646,14 @@ void Exogam::ReadConfiguration(NPL::InputParser parser) {
double ThetaZ = blocks[i]->GetDouble("ThetaZ", "deg");
AddDetector(X, Y, Z, ThetaX, ThetaY, ThetaZ);
}
+ else if (blocks[i]->HasTokenList(sphe)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Exogam " << i + 1 << endl;
+ double R = blocks[i]->GetDouble("R", "mm");
+ double Theta = blocks[i]->GetDouble("Theta", "deg");
+ double Phi = blocks[i]->GetDouble("Phi", "deg");
+ AddDetector(R, Theta, Phi);
+ }
else {
cout << "ERROR: check your input file formatting " << endl;
exit(1);
@@ -596,7 +691,7 @@ void Exogam::InitializeRootOutput() {
RootOutput* pAnalysis = RootOutput::getInstance();
TTree* pTree = pAnalysis->GetTree();
if (!pTree->FindBranch("Exogam")) {
- pTree->Branch("Exogam", "TExogamData", &m_Event);
+ pTree->Branch("Exogam", "TExogamCalData", &m_Event);
}
pTree->SetBranchAddress("Exogam", &m_Event);
}
@@ -612,17 +707,24 @@ void Exogam::ReadSensitive(const G4Event*) {
CalorimeterScorers::PS_Calorimeter* Scorer = (CalorimeterScorers::PS_Calorimeter*)m_ExogamScorer->GetPrimitive(0);
unsigned int size = Scorer->GetMult();
+ // if (size > 0)
+ // std::cout << "/////////: " << size << " "<< Scorer->GetName()<< std::endl;
for (unsigned int i = 0; i < size; i++) {
- double Energy = RandGauss::shoot(Scorer->GetEnergy(i), Exogam_NS::ResoEnergy);
+ // double Energy = RandGauss::shoot(Scorer->GetEnergy(i) * MeV, Exogam_NS::ResoEnergy);
+ double Energy = Scorer->GetEnergy(i) * MeV;
+
if (Energy > Exogam_NS::EnergyThreshold) {
double Time = Scorer->GetTime(i);
int CristalNbr = Scorer->GetLevel(i)[0];
int CloverNbr = Scorer->GetLevel(i)[1];
- //FIXME
- //m_Event->SetCristal(CristalNbr);
- //m_Event->SetClover(CloverNbr);
- //m_Event->SetEnergy(Energy);
- //m_Event->SetTime(Time);
+ // FIXME
+ // m_Event->SetCristal(CristalNbr);
+ // m_Event->SetClover(CloverNbr);
+ // m_Event->SetEnergy(Energy);
+ // m_Event->SetTime(Time);
+ int CristalID = CristalNbr + (CloverNbr - 1) * 4;
+
+ m_Event->SetExo(CristalID, Energy, -1000, Time, -1000, -1000, -1000, -1000, -1000, -1000, -1000);
}
}
}
diff --git a/NPSimulation/Detectors/Exogam/Exogam.hh b/NPSimulation/Detectors/Exogam/Exogam.hh
index fd2141a00..5f3921cba 100644
--- a/NPSimulation/Detectors/Exogam/Exogam.hh
+++ b/NPSimulation/Detectors/Exogam/Exogam.hh
@@ -34,7 +34,7 @@ using namespace std;
// NPTool header
#include "NPSVDetector.hh"
-#include "TExogamData.h"
+#include "TExogamCalData.h"
#include "NPInputParser.h"
class Exogam : public NPS::VDetector{
@@ -49,8 +49,10 @@ class Exogam : public NPS::VDetector{
/////// Specific Function of this Class ///////////
////////////////////////////////////////////////////
public:
- // Spherical
+ // Cartesian
void AddDetector(double X,double Y, double Z, double ThetaX, double ThetaY, double ThetaZ);
+ // Spherical
+ void AddDetector(double R,double Theta, double Phi);
G4int InitializeMaterials();
void BuildClover(int i_clo, G4LogicalVolume* world);
@@ -91,7 +93,7 @@ class Exogam : public NPS::VDetector{
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
private:
- TExogamData* m_Event ;
+ TExogamCalData* m_Event ;
////////////////////////////////////////////////////
///////////////Private intern Data//////////////////
@@ -105,6 +107,11 @@ class Exogam : public NPS::VDetector{
vector<double> m_ThetaX; //rotation angles to X, Y, Z axis
vector<double> m_ThetaY;
vector<double> m_ThetaZ;
+
+ // Detector Coordinate in spherical
+ vector<double> m_R;
+ vector<double> m_Theta;
+ vector<double> m_Phi;
private:
//materials
--
GitLab