From a1ab4c454949c0e02647b935e40912b19835620b Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 6 Jun 2024 09:04:21 +0100 Subject: [PATCH 01/21] Add fibre info to ratpac data structure This will allow us to pass on fiber info like angle and length to processors. --- src/core/include/RAT/Gsim.hh | 3 + src/core/src/Gsim.cc | 3 + src/ds/CMakeLists.txt | 1 + src/ds/include/RAT/DS/LinkDef.hh | 1 + src/ds/include/RAT/DS/NestedTubeInfo.hh | 92 +++++++++ src/ds/include/RAT/DS/Run.hh | 18 ++ src/geo/CMakeLists.txt | 3 + .../RAT/GeoNestedSolidArrayFactoryBase.hh | 22 ++ .../include/RAT/GeoNestedTubeArrayFactory.hh | 16 ++ .../include/RAT/GeoNestedTubeConstruction.hh | 63 ++++++ src/geo/src/GeoNestedSolidArrayFactoryBase.cc | 194 ++++++++++++++++++ src/geo/src/GeoNestedTubeArrayFactory.cc | 73 +++++++ src/geo/src/GeoNestedTubeConstruction.cc | 174 ++++++++++++++++ 13 files changed, 663 insertions(+) create mode 100644 src/ds/include/RAT/DS/NestedTubeInfo.hh create mode 100644 src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh create mode 100644 src/geo/include/RAT/GeoNestedTubeArrayFactory.hh create mode 100644 src/geo/include/RAT/GeoNestedTubeConstruction.hh create mode 100644 src/geo/src/GeoNestedSolidArrayFactoryBase.cc create mode 100644 src/geo/src/GeoNestedTubeArrayFactory.cc create mode 100644 src/geo/src/GeoNestedTubeConstruction.cc diff --git a/src/core/include/RAT/Gsim.hh b/src/core/include/RAT/Gsim.hh index ff04143d..23365866 100644 --- a/src/core/include/RAT/Gsim.hh +++ b/src/core/include/RAT/Gsim.hh @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -86,6 +87,8 @@ class Gsim : public Producer, G4UserRunAction, G4UserEventAction, G4UserTracking std::vector fPMTTime; //< PMT transit time/delay calculator (indexed by modeltype) std::vector fPMTCharge; //< PMT single-pe charge calculator (indexed by modeltype) + RAT::DS::NestedTubeInfo *fNestedTubeInfo; + RAT::DS::Run *run; int runID; TTimeStamp utc; diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index 8df286fc..c2d1baed 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -30,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -179,6 +180,7 @@ void Gsim::BeginOfRunAction(const G4Run * /*aRun*/) { run = DS::RunStore::GetRun(runID); fPMTInfo = run->GetPMTInfo(); + fNestedTubeInfo = run->GetNestedTubeInfo(); for (size_t i = 0; i < fPMTTime.size(); i++) { delete fPMTTime[i]; @@ -421,6 +423,7 @@ void Gsim::MakeRun(int _runID) { run->SetType((unsigned)lrun->GetI("runtype")); run->SetStartTime(utc); run->SetPMTInfo(&PMTFactoryBase::GetPMTInfo()); + run->SetNestedTubeInfo(&GeoNestedSolidArrayFactoryBase::GetNestedTubeInfo()); DS::RunStore::AddNewRun(run); } diff --git a/src/ds/CMakeLists.txt b/src/ds/CMakeLists.txt index b4e2f254..afa0b86c 100644 --- a/src/ds/CMakeLists.txt +++ b/src/ds/CMakeLists.txt @@ -37,6 +37,7 @@ root_generate_dictionary(G__RATDict RAT/DS/RunStore.hh RAT/DS/Run.hh RAT/DS/PMTInfo.hh + RAT/DS/NestedTubeInfo.hh RAT/DS/MCTrack.hh RAT/DS/MCTrackStep.hh RAT/DS/Calib.hh diff --git a/src/ds/include/RAT/DS/LinkDef.hh b/src/ds/include/RAT/DS/LinkDef.hh index b48566dc..adea4fb9 100644 --- a/src/ds/include/RAT/DS/LinkDef.hh +++ b/src/ds/include/RAT/DS/LinkDef.hh @@ -2,6 +2,7 @@ #pragma link C++ class RAT::DS::Root + ; #pragma link C++ class RAT::DS::PMTInfo + ; +#pragma link C++ class RAT::DS::NestedTubeInfo + ; #pragma link C++ class RAT::DS::MC + ; #pragma link C++ class RAT::DS::MCParticle + ; diff --git a/src/ds/include/RAT/DS/NestedTubeInfo.hh b/src/ds/include/RAT/DS/NestedTubeInfo.hh new file mode 100644 index 00000000..820255dc --- /dev/null +++ b/src/ds/include/RAT/DS/NestedTubeInfo.hh @@ -0,0 +1,92 @@ +/** + * @class DS::NestedTubeInfo + * Data Structure: Fiber properties + * + * Information about nested tubes (fibers), including positions, rotations, and lengths + */ + +#ifndef __RAT_DS_NestedTubeInfo__ +#define __RAT_DS_NestedTubeInfo__ + +#include +#include +#include + +namespace RAT { +namespace DS { + +class NestedTubeInfo : public TObject { + public: + NestedTubeInfo() : TObject() {}; + virtual ~NestedTubeInfo() = default; + + virtual void AddNestedTube(const G4ThreeVector& _pos, const G4ThreeVector& _dir, + const double _length, + const double _core_r, + const double _inner_r, + const double _outer_r, + const std::string _core_material, + const std::string _inner_material, + const std::string _outer_material) { + pos.push_back(_pos); + dir.push_back(_dir); + length.push_back(_length); + core_r.push_back(_core_r); + inner_r.push_back(_inner_r); + outer_r.push_back(_outer_r); + core_material.push_back(_core_material); + inner_material.push_back(_inner_material); + outer_material.push_back(_outer_material); + } + + virtual void AddNestedTube(const G4ThreeVector& _pos, const G4ThreeVector& _dir, const int _length) { + AddNestedTube(_pos, _dir, _length, 0.47, 0.485, 0.5, "", "", ""); + } + + virtual Int_t GetNestedTubeCount() const { return pos.size(); } + + virtual G4ThreeVector GetPosition(int id) const { return pos.at(id); } + virtual void SetPosition(int id, const G4ThreeVector& _pos) { pos.at(id) = _pos; } + + virtual G4ThreeVector GetDirection(int id) const { return dir.at(id); } + virtual void SetDirection(int id, const G4ThreeVector& _dir) { dir.at(id) = _dir; } + + virtual int GetLength(int id) const { return length.at(id); } + virtual void SetLength(int id, int _length) { length.at(id) = _length; } + + virtual double GetCoreR(int id) const { return core_r.at(id); } + virtual void SetCoreR(int id, double _core_r) { core_r.at(id) = _core_r; } + + virtual double GetInnerR(int id) const { return inner_r.at(id); } + virtual void SetInnerR(int id, double _inner_r) { inner_r.at(id) = _inner_r; } + + virtual double GetOuterR(int id) const { return outer_r.at(id); } + virtual void SetOuterR(int id, double _outer_r) { outer_r.at(id) = _outer_r; } + + virtual std::string GetCoreMaterial(int id) const { return core_material.at(id); } + virtual void SetCoreMaterial(int id, double _core_material) { core_material.at(id) = _core_material; } + + virtual std::string GetInnerMaterial(int id) const { return inner_material.at(id); } + virtual void SetInnerMaterial(int id, double _inner_material) { inner_material.at(id) = _inner_material; } + + virtual std::string GetOuterMaterial(int id) const { return outer_material.at(id); } + virtual void SetOuterMaterial(int id, double _outer_material) { outer_material.at(id) = _outer_material; } + + ClassDef(NestedTubeInfo, 2); + + protected: + std::vector pos; + std::vector dir; + std::vector length; + std::vector core_r; + std::vector inner_r; + std::vector outer_r; + std::vector core_material; + std::vector inner_material; + std::vector outer_material; +}; + +} // namespace DS +} // namespace RAT + +#endif diff --git a/src/ds/include/RAT/DS/Run.hh b/src/ds/include/RAT/DS/Run.hh index ab0e725a..036de3e7 100644 --- a/src/ds/include/RAT/DS/Run.hh +++ b/src/ds/include/RAT/DS/Run.hh @@ -12,6 +12,7 @@ #include #include +#include #include namespace RAT { @@ -50,6 +51,22 @@ class Run : public TObject { virtual bool ExistPMTInfo() { return !pmtinfo.empty(); } virtual void PrunePMTInfo() { pmtinfo.resize(0); } + /** Nested tube information */ + virtual NestedTubeInfo *GetNestedTubeInfo() { + if (nestedtubeinfo.empty()) { + nestedtubeinfo.resize(1); + } + return &nestedtubeinfo[0]; + } + virtual void SetNestedTubeInfo(const NestedTubeInfo *_nestedtubeinfo) { + if (nestedtubeinfo.empty()) { + nestedtubeinfo.resize(1); + } + nestedtubeinfo[0] = *_nestedtubeinfo; + } + virtual bool ExistNestedTubeInfo() { return !nestedtubeinfo.empty(); } + virtual void PruneNestedTubeInfo() { nestedtubeinfo.resize(0); } + ClassDef(Run, 2); protected: @@ -57,6 +74,7 @@ class Run : public TObject { ULong64_t type; TTimeStamp startTime; std::vector pmtinfo; + std::vector nestedtubeinfo; }; } // namespace DS diff --git a/src/geo/CMakeLists.txt b/src/geo/CMakeLists.txt index fb4da37d..65aec4d5 100644 --- a/src/geo/CMakeLists.txt +++ b/src/geo/CMakeLists.txt @@ -33,6 +33,9 @@ add_library(geo OBJECT src/GeoRevolutionChimneyFactory.cc src/GeoRevolutionFactory.cc src/GeoSolidArrayFactoryBase.cc + src/GeoNestedSolidArrayFactoryBase.cc + src/GeoNestedTubeArrayFactory.cc + src/GeoNestedTubeConstruction.cc src/GeoSolidFactory.cc src/GeoSphereFactory.cc src/GeoSurfaceFactory.cc diff --git a/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh new file mode 100644 index 00000000..785f4f8c --- /dev/null +++ b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh @@ -0,0 +1,22 @@ +#ifndef __RAT_GeoNestedSolidArrayFactoryBase__ +#define __RAT_GeoNestedSolidArrayFactoryBase__ + +#include +#include + +namespace RAT { +class GeoNestedSolidArrayFactoryBase : public GeoFactory { + public: + GeoNestedSolidArrayFactoryBase(const std::string &name) : GeoFactory(name){}; + static const DS::NestedTubeInfo &GetNestedTubeInfo() { return nestedtubeinfo; } + + protected: + virtual G4VPhysicalVolume *Construct(DBLinkPtr table); + + static DS::NestedTubeInfo nestedtubeinfo; +}; + +} // namespace RAT + +#endif + diff --git a/src/geo/include/RAT/GeoNestedTubeArrayFactory.hh b/src/geo/include/RAT/GeoNestedTubeArrayFactory.hh new file mode 100644 index 00000000..ad33b9de --- /dev/null +++ b/src/geo/include/RAT/GeoNestedTubeArrayFactory.hh @@ -0,0 +1,16 @@ +#ifndef __RAT_GeoNestedTubeArrayFactory__ +#define __RAT_GeoNestedTubeArrayFactory__ + +#include + +namespace RAT { +class GeoNestedTubeArrayFactory : public GeoNestedSolidArrayFactoryBase { + public: + GeoNestedTubeArrayFactory() : GeoNestedSolidArrayFactoryBase("nestedtubearray"){}; + using GeoNestedSolidArrayFactoryBase::Construct; + virtual G4VPhysicalVolume *Construct(RAT::DBLinkPtr table); +}; + +} // namespace RAT + +#endif diff --git a/src/geo/include/RAT/GeoNestedTubeConstruction.hh b/src/geo/include/RAT/GeoNestedTubeConstruction.hh new file mode 100644 index 00000000..cc0b2ff4 --- /dev/null +++ b/src/geo/include/RAT/GeoNestedTubeConstruction.hh @@ -0,0 +1,63 @@ +#ifndef __RAT_GeoNestedTubeConstruction__ +#define __RAT_GeoNestedTubeConstruction__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +struct GeoNestedTubeConstructionParams { + GeoNestedTubeConstructionParams() { + invisible = false; + }; + + bool invisible; + + double outer_r; + double inner_r; + double core_r; + double Dz; // half length + + G4Material *outer; + G4Material *inner; + G4Material *core; + + // G4OpticalSurface *outer_inner; + G4OpticalSurface *inner_core; + +}; + +class GeoNestedTubeConstruction { + public: + GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID); + virtual ~GeoNestedTubeConstruction() {} + + virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID); + virtual G4VSolid *BuildSolid(const std::string &prefix); + virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, + G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, + int copyNo); + + protected: + // physical volumes + G4PVPlacement *inner_phys; + G4PVPlacement *core_phys; + + G4LogicalVolume *log_tube; + GeoNestedTubeConstructionParams fParams; + + DBLinkPtr myTable; +}; + +} // namespace RAT + +#endif \ No newline at end of file diff --git a/src/geo/src/GeoNestedSolidArrayFactoryBase.cc b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc new file mode 100644 index 00000000..95c56f4c --- /dev/null +++ b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc @@ -0,0 +1,194 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +DS::NestedTubeInfo GeoNestedSolidArrayFactoryBase::nestedtubeinfo; + +G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { + std::string volume_name = table->GetIndex(); + std::string volume_name_one = volume_name + "one"; + std::string volume_name_two = volume_name + "two"; + std::string mother_name = table->GetS("mother"); + + G4LogicalVolume *mother = FindMother(mother_name); + if (mother == 0) { + Log::Die("GeoNestedTube: Unable to find mother volume " + mother_name + " for " + volume_name); + } + // build logical volume of nested tube. Contains volumes for inner and core. + // GeoNestedTubeConstruction *construction = new GeoNestedTubeConstruction(table, mother); + // G4LogicalVolume *log_tube = construction->BuildVolume(volume_name); + + // Read Solid positions + // TODO: The programmers guide gives an alternative for GetDArray that will be faster when you have large arrays -> our arrays will be on order 6k so we should switch + std::string pos_table_name = table->GetS("pos_table"); + DBLinkPtr lpos_table = DB::Get()->GetLink(pos_table_name); + const std::vector &pos_x = lpos_table->GetDArray("x"); + const std::vector &pos_y = lpos_table->GetDArray("y"); + const std::vector &pos_z = lpos_table->GetDArray("z"); + + // read max number of solids to use + int max_solids = pos_x.size(); // default to read all + try { + max_solids = table->GetI("max_num"); + } catch (DBNotFoundError &e) { + } + if (max_solids > (int)pos_x.size()) max_solids = pos_x.size(); + + // read starting number of solids to use + int start_solid_num = 0; // default to read all + try { + start_solid_num = table->GetI("start_num"); + } catch (DBNotFoundError &e) { + } + if (start_solid_num < 0) start_solid_num = 0; + + // Try to see if a sub type has been specified + int sub_type = -1; // default to read all + try { + sub_type = table->GetI("sub_type"); + } catch (DBNotFoundError &e) { + } + + std::vector sub_type_array; + for (int i = 0; i < max_solids; i++) sub_type_array.push_back(-1); + if (sub_type > -1) { + sub_type_array = lpos_table->GetIArray("sub_type"); + } + + // direction of individual solids. Default is that +z is orientation pointing + // direction optional, default is no rotation + bool rot_manual = false; + std::vector rot_x, rot_y, rot_z; + try { + std::string rotate_str = table->GetS("rotate_solids"); + if (rotate_str == "manual") rot_manual = true; + } catch (DBNotFoundError &e) { + } + if (rot_manual) { + rot_x = lpos_table->GetDArray("rot_x"); + rot_y = lpos_table->GetDArray("rot_y"); + rot_z = lpos_table->GetDArray("rot_z"); + } + + // Orientation of Solids + bool orient_manual = false; + try { + std::string orient_str = table->GetS("orientation"); + if (orient_str == "manual") + orient_manual = true; + else if (orient_str == "point") + orient_manual = false; + else + Log::Die("GeoBuilder error: Unknown solid orientation " + orient_str); + } catch (DBNotFoundError &e) { + } + + std::vector dir_x, dir_y, dir_z; + std::vector orient_point_array; + G4ThreeVector orient_point; + if (orient_manual) { + dir_x = lpos_table->GetDArray("dir_x"); + dir_y = lpos_table->GetDArray("dir_y"); + dir_z = lpos_table->GetDArray("dir_z"); + } else { + orient_point_array = table->GetDArray("orient_point"); + if (orient_point_array.size() != 3) Log::Die("GeoBuilder error: orient_point must have 3 values"); + orient_point.set(orient_point_array[0], orient_point_array[1], orient_point_array[2]); + } + + // Optionally can rescale Solid radius from mother volume center for + // case where Solids have spherical layout symmetry + + bool rescale_radius = false; + double new_radius = 1.0; + try { + new_radius = table->GetD("rescale_radius"); + rescale_radius = true; + } catch (DBNotFoundError &e) { + } + + // get pointer to physical mother volume + // ##outer_tank## + G4VPhysicalVolume *phys_mother = GeoFactory::FindPhysMother(mother_name); + + // create physical volumes for each fibre by placing logiSolid in mother volume + for (int solidID = start_solid_num; solidID < max_solids; solidID++) { + if ((sub_type == -1) || (sub_type == sub_type_array[solidID])) { + // construct + GeoNestedTubeConstruction *construction = new GeoNestedTubeConstruction(table, lpos_table, mother, solidID); + G4LogicalVolume *log_tube = construction->BuildVolume(volume_name, solidID); + // name + std::string tubename = volume_name + "_" + ::to_string(solidID); + + // position + G4ThreeVector tubepos(pos_x[solidID], pos_y[solidID], pos_z[solidID]); + + // direction + G4ThreeVector soliddir; + if (orient_manual) + soliddir.set(dir_x[solidID], dir_y[solidID], dir_z[solidID]); + else + soliddir = orient_point - tubepos; + soliddir = soliddir.unit(); + + // rescale + if (rescale_radius) tubepos.setMag(new_radius); + + // rotation required to point in direction of soliddir + double angle_y = (-1.0) * atan2(soliddir.x(), soliddir.z()); + double angle_x = atan2(soliddir.y(), sqrt(soliddir.x() * soliddir.x() + soliddir.z() * soliddir.z())); + double angle_z = atan2(-1 * soliddir.y() * soliddir.z(), soliddir.x()); + + G4RotationMatrix *tuberot = new G4RotationMatrix(); + + tuberot->rotateY(angle_y); + tuberot->rotateX(angle_x); + tuberot->rotateZ(angle_z); + + if (rot_manual) { + tuberot->rotateZ(rot_z[solidID] * CLHEP::deg); + tuberot->rotateY(rot_y[solidID] * CLHEP::deg); + tuberot->rotateX(rot_x[solidID] * CLHEP::deg); + } + // **************************************************************** + // * Use the constructor that specifies the PHYSICAL mother, since + // * each Solid occurs only once in one physical volume. This saves + // * the GeometryManager some work. -GHS. + // **************************************************************** + + // Write the real fiber positions and directions. + // This goes into the DS by way of Gsim + double length = lpos_table->GetDArray("Dz")[solidID]*2; + double core_r = table->GetD("core_r"); + double inner_r = table->GetD("inner_r"); + double outer_r = table->GetD("outer_r"); + std::string core_material = table->GetS("material_core"); + std::string inner_material = table->GetS("material_inner"); + std::string outer_material = table->GetS("material_outer"); + nestedtubeinfo.AddNestedTube(tubepos, soliddir, length, + core_r, inner_r, outer_r, + core_material, inner_material, outer_material); + + // instance of physical volume for fibre inside mother volume + construction->PlaceNestedTube(tuberot, tubepos, tubename, log_tube, phys_mother, false, solidID); + + } // end loop over solidID + } + return 0; +} + +} // namespace RAT diff --git a/src/geo/src/GeoNestedTubeArrayFactory.cc b/src/geo/src/GeoNestedTubeArrayFactory.cc new file mode 100644 index 00000000..4fe4dace --- /dev/null +++ b/src/geo/src/GeoNestedTubeArrayFactory.cc @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#include +#include + +namespace RAT { + +G4VPhysicalVolume *GeoNestedTubeArrayFactory::Construct(DBLinkPtr table) { + std::string volume_name = table->GetIndex(); + + info << "GeoNestedTubeArrayFactory: Constructing volume " + volume_name << newline; + + // Optional parameters + G4double r_min = 0.0; + try { + r_min = table->GetD("r_min") * CLHEP::mm; + } catch (DBNotFoundError &e) { + }; + G4double phi_start = 0.0; + try { + phi_start = table->GetD("phi_start") * CLHEP::deg; + } catch (DBNotFoundError &e) { + }; + G4double phi_delta = CLHEP::twopi; + try { + phi_delta = table->GetD("phi_delta") * CLHEP::deg; + } catch (DBNotFoundError &e) { + }; + + // can cut out a spherical region from all the solids of + // radius sphere_cut_r. + // requires that rescale_r be std::set. + G4double s_cut_r = -1.0; + try { + s_cut_r = table->GetD("sphere_cut_r") * CLHEP::mm; + } catch (DBNotFoundError &e) { + }; + + // can rescale Solid radius from mother volume center for + // case where Solids have spherical layout symmetry + G4double rescale_r = -1.0; + try { + rescale_r = table->GetD("rescale_radius") * CLHEP::mm; + } catch (DBNotFoundError &e) { + }; + + int preflip = 0; + try { + preflip = table->GetI("preflip"); + } catch (DBNotFoundError &e) { + }; + + // End optional parameters + + if ((s_cut_r > 0) && (rescale_r > 0)) { + G4VSolid *sphere_cutter = new G4Orb("temp_sphere", s_cut_r); // This is the cut out piece + + G4RotationMatrix *sphererot = new G4RotationMatrix(); + + G4ThreeVector spherepos(0.0, 0.0, -1 * rescale_r); + } + + if (preflip) { + G4RotationMatrix *fliprot = new G4RotationMatrix(G4ThreeVector(1, 0, 0), CLHEP::pi); + } + + //TODO: all the above parameters are ignored and only the table parameters are used. For wider usability the above parameters should be included in the array construction + return GeoNestedSolidArrayFactoryBase::Construct(table); +} + +} // namespace RAT diff --git a/src/geo/src/GeoNestedTubeConstruction.cc b/src/geo/src/GeoNestedTubeConstruction.cc new file mode 100644 index 00000000..1aff8bee --- /dev/null +++ b/src/geo/src/GeoNestedTubeConstruction.cc @@ -0,0 +1,174 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +GeoNestedTubeConstruction::GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID) { + inner_phys = 0; + core_phys = 0; + + log_tube = 0; + + myTable = table; + + // Setup NestedTube parameters + fParams.outer_r = table->GetD("outer_r"); + fParams.inner_r = table->GetD("inner_r"); + fParams.core_r = table->GetD("core_r"); + + // Materials + fParams.outer = G4Material::GetMaterial(table->GetS("material_outer")); + fParams.inner = G4Material::GetMaterial(table->GetS("material_inner")); + fParams.core = G4Material::GetMaterial(table->GetS("material_core")); + DB *db = DB::Get(); + std::string core_surface_name = table->GetS("material_core"); + if (Materials::optical_surface.count(core_surface_name) == 0) + Log::Die("GeoSolidFactory: Surface " + core_surface_name + " does not exist"); + fParams.inner_core = Materials::optical_surface[core_surface_name]; + fParams.Dz = postable->GetDArray("Dz")[ID]; + + std::string tube_name = table->GetS("index"); + Log::Assert(fParams.outer_r > 0, "GeoNestedTubeConstruction: " + tube_name + " outer radius must be positive"); + Log::Assert(fParams.inner_r > 0, "GeoNestedTubeConstruction: " + tube_name + " inner radius must be positive"); + Log::Assert(fParams.core_r > 0, "GeoNestedTubeConstruction: " + tube_name + " core radius must be positive"); + Log::Assert(fParams.outer_r > fParams.inner_r, + "GeoNestedTubeConstruction: " + tube_name + " outer radius is smaller than inner radius"); + Log::Assert(fParams.inner_r > fParams.core_r, + "GeoNestedTubeConstruction: " + tube_name + " inner radius is smaller than core radius"); + Log::Assert(fParams.outer, "GeoNestedTubeConstruction: " + tube_name + " has an invalid outer material"); + Log::Assert(fParams.inner, "GeoNestedTubeConstruction: " + tube_name + " has an invalid inner material"); + Log::Assert(fParams.core, "GeoNestedTubeConstruction: " + tube_name + " has an invalid core material"); + Log::Assert(fParams.inner_core, "GeoNestedTubeConstruction: " + tube_name + " has an invalid core surface material"); +} + +G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefix, int ID) { + if (log_tube) { + return log_tube; + } + + // fibre outer + G4Tubs *outer_solid = (G4Tubs *)BuildSolid(prefix + "_outer_solid"); + + // fibre inner + G4Tubs *inner_solid = + new G4Tubs(prefix + "_inner_solid", 0.0, fParams.inner_r, fParams.Dz, 0.0, CLHEP::twopi); + + // fibre core + G4Tubs *core_solid = + new G4Tubs(prefix + "_core_solid", 0.0, fParams.core_r, fParams.Dz, 0.0, CLHEP::twopi); + + // ------------ Logical Volumes ------------- + G4LogicalVolume *outer_log, *inner_log, *core_log; + + outer_log = new G4LogicalVolume(outer_solid, fParams.outer, prefix + "_outer_logic"); + inner_log = new G4LogicalVolume(inner_solid, fParams.inner, prefix + "_inner_logic"); + core_log = new G4LogicalVolume(core_solid, fParams.core, prefix + "_core_logic"); + + // ------------ Physical Volumes ------------- + G4ThreeVector noTranslation(0., 0., 0.); + + // Place the core solids in the inner solid to produce the physical volumes + inner_phys = new G4PVPlacement(0, // no rotation + noTranslation, // place inner tube concentric to outer + inner_log, // the logical volume + prefix + "_inner_phys", // a name for this physical volume + outer_log, // the mother volume + false, // no boolean ops + 0); // copy number + + core_phys = new G4PVPlacement(0, // no rotation + noTranslation, // place inner tube concentric to outer + core_log, // the logical volume + prefix + "_" + std::to_string(ID) + "_core", // a name for this physical volume + inner_log, // the mother volume + false, // no boolean ops + 0); // copy number + + // ------------ Vis Attributes ------------- + G4VisAttributes *vis = new G4VisAttributes(); + G4VisAttributes *vis_inner = new G4VisAttributes(); + G4VisAttributes *vis_core = new G4VisAttributes(); + try { + const std::vector &color = myTable->GetDArray("color"); + if (color.size() == 3) { // RGB + vis->SetColour(G4Colour(color[0], color[1], color[2])); + vis_inner->SetColour(G4Colour(color[2], color[0], color[1])); + vis_core->SetColour(G4Colour(color[1], color[2], color[0])); + } + else if (color.size() == 4) { // RGBA + vis->SetColour(G4Colour(color[0], color[1], color[2], color[3])); + vis_inner->SetColour(G4Colour(color[2], color[0], color[1], color[3])); + vis_core->SetColour(G4Colour(color[1], color[2], color[0], color[3])); + } + else + warn << "GeoNestedTubeConstruction error: " << myTable->GetName() << "[" << myTable->GetIndex() + << "].color must have 3 or 4 components" << newline; + } catch (DBNotFoundError &e) { + }; + try { + std::string drawstyle = myTable->GetS("drawstyle"); + if (drawstyle == "wireframe") { + vis->SetForceWireframe(true); + vis_inner->SetForceWireframe(true); + vis_core->SetForceWireframe(true); + } + else if (drawstyle == "solid") { + vis->SetForceSolid(true); + vis_inner->SetForceSolid(true); + vis_core->SetForceSolid(true); + } + else + warn << "GeoNestedTubeConstruction error: " << myTable->GetName() << "[" << myTable->GetIndex() + << "].drawstyle must be either \"wireframe\" or \"solid\"."; + } catch (DBNotFoundError &e) { + }; + + // Check for invisible flag last + try { + int invisible = myTable->GetI("invisible"); + if (invisible) + { + outer_log->SetVisAttributes(G4VisAttributes::GetInvisible()); + inner_log->SetVisAttributes(G4VisAttributes::GetInvisible()); + core_log->SetVisAttributes(G4VisAttributes::GetInvisible()); + } + } catch (DBNotFoundError &e) { + }; + + outer_log->SetVisAttributes(vis); + inner_log->SetVisAttributes(vis_inner); + core_log->SetVisAttributes(vis_core); + + log_tube = outer_log; + + return log_tube; +} + +G4VSolid *GeoNestedTubeConstruction::BuildSolid(const std::string &name) { + G4Tubs *outer = new G4Tubs(name, 0.0, fParams.outer_r, fParams.Dz, 0.0, CLHEP::twopi); + return outer; +} + +G4PVPlacement *GeoNestedTubeConstruction::PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, + G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, + bool booleanSolid, int copyNo) { + G4PVPlacement *outer_phys = new G4PVPlacement(tuberot, tubepos, name, logi_tube, mother_phys, booleanSolid, copyNo); + + // core surface + new G4LogicalBorderSurface(name + "_core_logsurf1", core_phys, inner_phys, fParams.inner_core); + + // build the inner surface + // new G4LogicalBorderSurface(name + "_inner_logsurf1", outer_phys, inner_phys, fParams.outer_inner); + + return outer_phys; +} + +} // namespace RAT \ No newline at end of file From 3717658c1069f65fcc945792f3478f61008fc878 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Mon, 17 Jun 2024 14:27:15 +0100 Subject: [PATCH 02/21] Add nested tube arrays to docs --- doc/users_guide/data_structure.rst | 2 +- doc/users_guide/database.rst | 4 ++-- doc/users_guide/geometry.rst | 23 +++++++++++++++++++++++ doc/users_guide/processes.rst | 2 +- 4 files changed, 27 insertions(+), 4 deletions(-) diff --git a/doc/users_guide/data_structure.rst b/doc/users_guide/data_structure.rst index eeaab6fe..2513e407 100644 --- a/doc/users_guide/data_structure.rst +++ b/doc/users_guide/data_structure.rst @@ -13,7 +13,7 @@ The Ratpac data structure is defined in the `src/ds` directory and lists everything under the `RAT::DS` namespace. An instance of the data structure is defined in the `RAT::DS::Root` object. The data structure is tree-like, with each instance of `RAT::DS::Root` containing a list of `RAT::DS::MC` objects which contain -the Monte Carl truth infomration, and a list of `RAT::DS::EV` objects which contain +the Monte Carl truth information, and a list of `RAT::DS::EV` objects which contain the reconstructed event information (usually after going through a processor that simulates the detector response). diff --git a/doc/users_guide/database.rst b/doc/users_guide/database.rst index d6c187ae..c0aa7dcb 100644 --- a/doc/users_guide/database.rst +++ b/doc/users_guide/database.rst @@ -63,7 +63,7 @@ Normally, tables and fields are all you have to think about, but RATDB also addresses an additional complication: overriding constants. It is a common use case to have default values of constants, values which are only valid in certain time intervals and can change (like the optical properties of the -scintilator), and user-specified values which are intended to override +scintillator), and user-specified values which are intended to override everything. RATDB handles this by internally grouping tables into three ''planes''. The @@ -79,7 +79,7 @@ When an item is requested, RATDB will attempt to locate it in the user plane first, then the time plane, and finally the default plane. Note that this is all handled in the background. You simply request the index_of_refraction field in the MEDIA[acrylic] table, and RATDB figures out the appropriate plane -from which to retreive the data. +from which to retrieve the data. How do I load data into RATDB? `````````````````````````````` diff --git a/doc/users_guide/geometry.rst b/doc/users_guide/geometry.rst index 8be611c0..0be64013 100644 --- a/doc/users_guide/geometry.rst +++ b/doc/users_guide/geometry.rst @@ -80,6 +80,7 @@ Allowed types: * tubearray - Array of tubes * lgarray - Array of tubes where one end has the PMT face cut out * pmtarray - Array of PMTs + * nestedtubearray - Array of three nested tubes. Useful to simulate optical fibers * waterboxarray - Array of standard cubitainer water boxes * extpolyarray - Array of extruded polygonal solids * bubble - Collection of bubbles @@ -147,6 +148,28 @@ PMTArray Fields: ``rescale_radius`` ``float`` (optional) Assumes all PMTs are spherically arranged around the center of the mother volume and rescales their positions to a particular radius. By default, no rescaling is done. ====================== ========================== =================== +NestedTubeArray Fields: + +====================== ========================== =================== +**Field** **Type** **Description** +====================== ========================== =================== +``pos_table`` ``string`` Specifies the table containing position (and direction) arrays specifying how to place PMTs +``core_r`` ``float`` The radius of the core tube (mm) +``inner_r`` ``float`` The radius of the inner tube (mm) +``outer_r`` ``float`` The radius of the outer tube (mm) +``material_core`` ``string`` The material of the core tube +``material_inner`` ``string`` The material of the inner tube +``material_outer`` ``string`` The material of the outer tube +``Dz`` ``float`` Half-height of tube (mm) +``phi_start`` ``float`` (optional) Angle (deg) where tube segment starts. Default is 0.0 +``phi_delta`` ``float`` (optional) Angle span (deg) of tube segment. Default is 360.0 +``start_idx`` ``int`` (optional) Index to start building nested tubes in the ``NESTEDTUBEINFO`` table specified (inclusive, defaults to 0) +``end_idx`` ``int`` (optional) Index to stop building nested tubes in the ``NESTEDTUBEINFO`` table specified (inclusive, defaults to length-1) +``orientation`` ``string`` Method of determining nested tube direction. "point" will aim all nested tubes at a point in space. "manual" requires that the position table also contain dir_x, dir_y, and dir_z fields which define the direction vector for each PMT. +``orient_point`` ``float[3]`` (optional) Point (mm) in mother volume to aim all tubes toward. +``rescale_radius`` ``float`` (optional) Assumes all tubes are spherically arranged around the center of the mother volume and rescales their positions to a particular radius. By default, no rescaling is done. +====================== ========================== =================== + Creating a parameterized geometry ````````````````````````````````` Using a ``DetectorFactory`` one can build a DB defined geometry on the fly diff --git a/doc/users_guide/processes.rst b/doc/users_guide/processes.rst index 0faa103e..d3491651 100644 --- a/doc/users_guide/processes.rst +++ b/doc/users_guide/processes.rst @@ -128,7 +128,7 @@ below. Wavelength Shifting ``````````````````` There are a few ways of doing bulk wavelength shifting in RAT. The default -behavior is for GLG4Scint to handle opticalphotons as well as charged +behavior is for GLG4Scint to handle optical photons as well as charged particles. Alternatively, you can also let GLG4Scint handle the primary scintillation, then use Geant4's G4OpWLS process or the custom BNLOpWLSModel to do the reemission. From ba998df12de211bb8a050a394d0b1dab2a141b93 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Mon, 17 Jun 2024 15:54:31 +0100 Subject: [PATCH 03/21] Add nested tube array factory to geoBuilder Also added a validation macro and geometry for the nested tubes. --- macros/vis_nestedTubeArrays.mac | 53 +++++++++ ratdb/Validation/NestedTubeArray.geo | 170 +++++++++++++++++++++++++++ src/geo/src/GeoBuilder.cc | 2 + 3 files changed, 225 insertions(+) create mode 100644 macros/vis_nestedTubeArrays.mac create mode 100644 ratdb/Validation/NestedTubeArray.geo diff --git a/macros/vis_nestedTubeArrays.mac b/macros/vis_nestedTubeArrays.mac new file mode 100644 index 00000000..6e54e047 --- /dev/null +++ b/macros/vis_nestedTubeArrays.mac @@ -0,0 +1,53 @@ +/glg4debug/glg4param omit_muon_processes 0.0 +/glg4debug/glg4param omit_hadronic_processes 0.0 + + +#set the detector parameters +/rat/db/set DETECTOR experiment "Validation" +/rat/db/set DETECTOR geo_file "Validation/NestedTubeArray.geo" + +# Colors +/rat/db/set GEO[world] invisible 1 +/rat/db/set GEO[outer_vessel] invisible 1 +/rat/db/set GEO[outer_tank] invisible 1 +#/rat/db/set GEO[outer_vessel] color [0.1,0.1,0.6,0.95] +#/rat/db/set GEO[outer_tank] color [0.42,0.47,0.57,0.3] +#/rat/db/set GEO[fiber_1_inner] color [0.0,0.8,0.0,0.10] + +/run/initialize + +/tracking/storeTrajectory 1 + +##### Visualization ########################## + +/vis/open OGLSQt +/vis/scene/create +/vis/scene/add/trajectories rich smooth +/tracking/storeTrajectory 1 +/tracking/FillPointCont 1 +/vis/scene/add/volume +/vis/scene/add/hits +/vis/sceneHandler/attach +/vis/viewer/set/upVector 0.0 0.0 1.0 +/vis/viewer/set/viewpointThetaPhi -90 135 +/vis/viewer/set/style s +/vis/viewer/flush + +/rat/proc count +/rat/procset update 10 + + +#add ntuple tracking +/rat/proc outntuple +/rat/procset include_tracking 1 +/rat/procset include_mcparticles 1 +/rat/procset include_pmthits 1 +/rat/procset include_untriggered_events 1 +/rat/physics/setOpWLS g4 +##### GENERATORS ################# +/generator/add combo pbomb:point:poisson +/generator/vtx/set 10 430 # 10000 photons, 430nm +/generator/pos/set 0.0 0.0 0.0 + +##### RUN ########### +/run/beamOn 1 \ No newline at end of file diff --git a/ratdb/Validation/NestedTubeArray.geo b/ratdb/Validation/NestedTubeArray.geo new file mode 100644 index 00000000..533ea83e --- /dev/null +++ b/ratdb/Validation/NestedTubeArray.geo @@ -0,0 +1,170 @@ +{ + name: "GEO", + index: "world", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "", + type: "box", + size: [20000.0,20000.0,20000.0], + material: "mirror", +} + +{ + name: "GEO", + index: "outer_vessel", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "world", + type: "box", + size: [10.0,10.0,10.0], + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "aluminum", + color: [0.02,0.2,0.2,0.03], +} + +{ + name: "GEO", + index: "outer_tank", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "outer_vessel", + type: "box", + size: [9.0,9.0,9.0], + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "mirror", + #color: [0.02,0.2,0.2,0.1], +} + +// here we place nested tubes manually +{ + name: "GEO", + index: "fiber_0_outer", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "outer_tank", + type: "tube", + r_max: 0.5, + size_z: 49.5, + position: [-5.0, 0.0, -5.0], + rotation: [-90.0, 0.0, 0.0], + material: "aluminum", + color: [0.0,0.8,0.0,0.01], +} + +{ + name: "GEO", + index: "fiber_0_inner", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "fiber_0_outer", + type: "tube", + r_max: 0.485, + size_z: 49.5, + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "glass", + color: [0.0,0.8,0.0,0.05], +} + +{ + name: "GEO", + index: "fiber_0_core", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "fiber_0_inner", + type: "tube", + r_max: 0.47, + size_z: 49.5, + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "mirror", + color: [0.0,0.8,0.0,0.1], +} +{ + name: "GEO", + index: "fiber_1_outer", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "outer_tank", + type: "tube", + r_max: 0.5, + size_z: 29.5, + position: [-5.0, 0.0, 5.0], + rotation: [-90.0, 0.0, 0.0], + material: "aluminum", + color: [0.0,0.8,0.0,0.01], +} + +{ + name: "GEO", + index: "fiber_1_inner", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "fiber_1_outer", + type: "tube", + r_max: 0.485, + size_z: 29.5, + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "glass", + color: [0.0,0.8,0.0,0.05], +} + +{ + name: "GEO", + index: "fiber_1_core", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "fiber_1_inner", + type: "tube", + r_max: 0.47, + size_z: 29.5, + position: [0.0, 0.0, 0.0], + rotation: [0.0, 0.0, 0.0], + material: "mirror", + color: [0.0,0.8,0.0,0.1], +} +// manual nested tubes ends here + +// this table is the information +// about where to place the nested tubes +// in the array +{ +name: "cable_pos", +valid_begin: [0], +valid_end: [0], +x: [5, 5] +y: [0, 0], +z: [5, -5], +dir_x: [0, 0], +dir_y: [1, 1], +dir_z: [0, 0], +Dz: [49.5, 29.5] +} + +// this is the definition of the +// nested tube 'properties' +// note that we reference the pos_table +// by name, also the pos_table +// doesn't have to be in the .geo file +{ +name: "GEO", +index: "fibers", +enable: 1, +valid_begin: [0, 0], +valid_end: [0, 0], +mother: "outer_tank", +type: "nestedtubearray", +core_r: 0.47, +inner_r: 0.485, +outer_r: 0.5, +pos_table: "cable_pos", +orientation: "manual", +material_outer: "aluminum", +material_inner: "glass", +material_core: "mirror", +#drawstyle: "solid", +color: [0.8,0.0,0.0,0.8] +} diff --git a/src/geo/src/GeoBuilder.cc b/src/geo/src/GeoBuilder.cc index d28c0e00..6818d017 100644 --- a/src/geo/src/GeoBuilder.cc +++ b/src/geo/src/GeoBuilder.cc @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -44,6 +45,7 @@ GeoBuilder::GeoBuilder() { new GeoReflectorFactory(); new GeoReflectorWaveguideFactory(); new PMTArrayFactory(); + new GeoNestedTubeArrayFactory(); new PMTCoverageFactory(); new GeoWaterBoxArrayFactory(); new GeoBubbleFactory(); From de896b8cf6ca748e106d0e81baf62abdab8de0d1 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Mon, 24 Jun 2024 11:24:30 +0100 Subject: [PATCH 04/21] Add clang-formatting to new files --- src/core/include/RAT/Gsim.hh | 2 +- src/core/src/Gsim.cc | 2 +- src/ds/include/RAT/DS/NestedTubeInfo.hh | 17 +++--- src/ds/include/RAT/DS/Run.hh | 2 +- .../RAT/GeoNestedSolidArrayFactoryBase.hh | 5 +- .../include/RAT/GeoNestedTubeConstruction.hh | 11 ++-- src/geo/src/GeoBuilder.cc | 2 +- src/geo/src/GeoNestedSolidArrayFactoryBase.cc | 22 +++---- src/geo/src/GeoNestedTubeArrayFactory.cc | 3 +- src/geo/src/GeoNestedTubeConstruction.cc | 61 +++++++++---------- 10 files changed, 58 insertions(+), 69 deletions(-) diff --git a/src/core/include/RAT/Gsim.hh b/src/core/include/RAT/Gsim.hh index 23365866..9a23c172 100644 --- a/src/core/include/RAT/Gsim.hh +++ b/src/core/include/RAT/Gsim.hh @@ -9,8 +9,8 @@ #include #include #include -#include #include +#include #include #include #include diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index c2d1baed..188a024a 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -30,7 +31,6 @@ #include #include #include -#include #include #include #include diff --git a/src/ds/include/RAT/DS/NestedTubeInfo.hh b/src/ds/include/RAT/DS/NestedTubeInfo.hh index 820255dc..24e372e0 100644 --- a/src/ds/include/RAT/DS/NestedTubeInfo.hh +++ b/src/ds/include/RAT/DS/NestedTubeInfo.hh @@ -2,13 +2,14 @@ * @class DS::NestedTubeInfo * Data Structure: Fiber properties * - * Information about nested tubes (fibers), including positions, rotations, and lengths + * Information about nested tubes (fibers), including positions, rotations, and lengths */ #ifndef __RAT_DS_NestedTubeInfo__ #define __RAT_DS_NestedTubeInfo__ #include + #include #include @@ -20,14 +21,10 @@ class NestedTubeInfo : public TObject { NestedTubeInfo() : TObject() {}; virtual ~NestedTubeInfo() = default; - virtual void AddNestedTube(const G4ThreeVector& _pos, const G4ThreeVector& _dir, - const double _length, - const double _core_r, - const double _inner_r, - const double _outer_r, - const std::string _core_material, - const std::string _inner_material, - const std::string _outer_material) { + virtual void AddNestedTube(const G4ThreeVector& _pos, const G4ThreeVector& _dir, const double _length, + const double _core_r, const double _inner_r, const double _outer_r, + const std::string _core_material, const std::string _inner_material, + const std::string _outer_material) { pos.push_back(_pos); dir.push_back(_dir); length.push_back(_length); @@ -87,6 +84,6 @@ class NestedTubeInfo : public TObject { }; } // namespace DS -} // namespace RAT +} // namespace RAT #endif diff --git a/src/ds/include/RAT/DS/Run.hh b/src/ds/include/RAT/DS/Run.hh index 036de3e7..3792b7ca 100644 --- a/src/ds/include/RAT/DS/Run.hh +++ b/src/ds/include/RAT/DS/Run.hh @@ -11,8 +11,8 @@ #include #include -#include #include +#include #include namespace RAT { diff --git a/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh index 785f4f8c..babe0be7 100644 --- a/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh +++ b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh @@ -1,13 +1,13 @@ #ifndef __RAT_GeoNestedSolidArrayFactoryBase__ #define __RAT_GeoNestedSolidArrayFactoryBase__ -#include #include +#include namespace RAT { class GeoNestedSolidArrayFactoryBase : public GeoFactory { public: - GeoNestedSolidArrayFactoryBase(const std::string &name) : GeoFactory(name){}; + GeoNestedSolidArrayFactoryBase(const std::string &name) : GeoFactory(name) {}; static const DS::NestedTubeInfo &GetNestedTubeInfo() { return nestedtubeinfo; } protected: @@ -19,4 +19,3 @@ class GeoNestedSolidArrayFactoryBase : public GeoFactory { } // namespace RAT #endif - diff --git a/src/geo/include/RAT/GeoNestedTubeConstruction.hh b/src/geo/include/RAT/GeoNestedTubeConstruction.hh index cc0b2ff4..a1b5401d 100644 --- a/src/geo/include/RAT/GeoNestedTubeConstruction.hh +++ b/src/geo/include/RAT/GeoNestedTubeConstruction.hh @@ -16,16 +16,14 @@ namespace RAT { struct GeoNestedTubeConstructionParams { - GeoNestedTubeConstructionParams() { - invisible = false; - }; + GeoNestedTubeConstructionParams() { invisible = false; }; bool invisible; double outer_r; double inner_r; double core_r; - double Dz; // half length + double Dz; // half length G4Material *outer; G4Material *inner; @@ -33,7 +31,6 @@ struct GeoNestedTubeConstructionParams { // G4OpticalSurface *outer_inner; G4OpticalSurface *inner_core; - }; class GeoNestedTubeConstruction { @@ -44,8 +41,8 @@ class GeoNestedTubeConstruction { virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID); virtual G4VSolid *BuildSolid(const std::string &prefix); virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, - G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, - int copyNo); + G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, + int copyNo); protected: // physical volumes diff --git a/src/geo/src/GeoBuilder.cc b/src/geo/src/GeoBuilder.cc index 6818d017..4e3e2f17 100644 --- a/src/geo/src/GeoBuilder.cc +++ b/src/geo/src/GeoBuilder.cc @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -29,7 +30,6 @@ #include #include #include -#include #include #include #include diff --git a/src/geo/src/GeoNestedSolidArrayFactoryBase.cc b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc index 95c56f4c..0f633f75 100644 --- a/src/geo/src/GeoNestedSolidArrayFactoryBase.cc +++ b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc @@ -6,12 +6,12 @@ #include #include #include +#include +#include #include #include #include #include -#include -#include #include namespace RAT { @@ -33,7 +33,8 @@ G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { // G4LogicalVolume *log_tube = construction->BuildVolume(volume_name); // Read Solid positions - // TODO: The programmers guide gives an alternative for GetDArray that will be faster when you have large arrays -> our arrays will be on order 6k so we should switch + // TODO: The programmers guide gives an alternative for GetDArray that will be faster when you have large arrays -> + // our arrays will be on order 6k so we should switch std::string pos_table_name = table->GetS("pos_table"); DBLinkPtr lpos_table = DB::Get()->GetLink(pos_table_name); const std::vector &pos_x = lpos_table->GetDArray("x"); @@ -112,7 +113,7 @@ G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { // Optionally can rescale Solid radius from mother volume center for // case where Solids have spherical layout symmetry - + bool rescale_radius = false; double new_radius = 1.0; try { @@ -169,26 +170,25 @@ G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { // * each Solid occurs only once in one physical volume. This saves // * the GeometryManager some work. -GHS. // **************************************************************** - + // Write the real fiber positions and directions. // This goes into the DS by way of Gsim - double length = lpos_table->GetDArray("Dz")[solidID]*2; + double length = lpos_table->GetDArray("Dz")[solidID] * 2; double core_r = table->GetD("core_r"); double inner_r = table->GetD("inner_r"); double outer_r = table->GetD("outer_r"); std::string core_material = table->GetS("material_core"); std::string inner_material = table->GetS("material_inner"); std::string outer_material = table->GetS("material_outer"); - nestedtubeinfo.AddNestedTube(tubepos, soliddir, length, - core_r, inner_r, outer_r, - core_material, inner_material, outer_material); + nestedtubeinfo.AddNestedTube(tubepos, soliddir, length, core_r, inner_r, outer_r, core_material, inner_material, + outer_material); // instance of physical volume for fibre inside mother volume construction->PlaceNestedTube(tuberot, tubepos, tubename, log_tube, phys_mother, false, solidID); - + } // end loop over solidID } return 0; } -} // namespace RAT +} // namespace RAT diff --git a/src/geo/src/GeoNestedTubeArrayFactory.cc b/src/geo/src/GeoNestedTubeArrayFactory.cc index 4fe4dace..46ee759a 100644 --- a/src/geo/src/GeoNestedTubeArrayFactory.cc +++ b/src/geo/src/GeoNestedTubeArrayFactory.cc @@ -66,7 +66,8 @@ G4VPhysicalVolume *GeoNestedTubeArrayFactory::Construct(DBLinkPtr table) { G4RotationMatrix *fliprot = new G4RotationMatrix(G4ThreeVector(1, 0, 0), CLHEP::pi); } - //TODO: all the above parameters are ignored and only the table parameters are used. For wider usability the above parameters should be included in the array construction + // TODO: all the above parameters are ignored and only the table parameters are used. For wider usability the above + // parameters should be included in the array construction return GeoNestedSolidArrayFactoryBase::Construct(table); } diff --git a/src/geo/src/GeoNestedTubeConstruction.cc b/src/geo/src/GeoNestedTubeConstruction.cc index 1aff8bee..a88142d4 100644 --- a/src/geo/src/GeoNestedTubeConstruction.cc +++ b/src/geo/src/GeoNestedTubeConstruction.cc @@ -1,17 +1,18 @@ -#include #include #include -#include #include +#include +#include +#include #include #include -#include #include #include namespace RAT { -GeoNestedTubeConstruction::GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID) { +GeoNestedTubeConstruction::GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, + int ID) { inner_phys = 0; core_phys = 0; @@ -54,16 +55,14 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi return log_tube; } - // fibre outer + // fibre outer G4Tubs *outer_solid = (G4Tubs *)BuildSolid(prefix + "_outer_solid"); - // fibre inner - G4Tubs *inner_solid = - new G4Tubs(prefix + "_inner_solid", 0.0, fParams.inner_r, fParams.Dz, 0.0, CLHEP::twopi); + // fibre inner + G4Tubs *inner_solid = new G4Tubs(prefix + "_inner_solid", 0.0, fParams.inner_r, fParams.Dz, 0.0, CLHEP::twopi); // fibre core - G4Tubs *core_solid = - new G4Tubs(prefix + "_core_solid", 0.0, fParams.core_r, fParams.Dz, 0.0, CLHEP::twopi); + G4Tubs *core_solid = new G4Tubs(prefix + "_core_solid", 0.0, fParams.core_r, fParams.Dz, 0.0, CLHEP::twopi); // ------------ Logical Volumes ------------- G4LogicalVolume *outer_log, *inner_log, *core_log; @@ -76,7 +75,7 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi G4ThreeVector noTranslation(0., 0., 0.); // Place the core solids in the inner solid to produce the physical volumes - inner_phys = new G4PVPlacement(0, // no rotation + inner_phys = new G4PVPlacement(0, // no rotation noTranslation, // place inner tube concentric to outer inner_log, // the logical volume prefix + "_inner_phys", // a name for this physical volume @@ -84,13 +83,13 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi false, // no boolean ops 0); // copy number - core_phys = new G4PVPlacement(0, // no rotation - noTranslation, // place inner tube concentric to outer - core_log, // the logical volume - prefix + "_" + std::to_string(ID) + "_core", // a name for this physical volume - inner_log, // the mother volume - false, // no boolean ops - 0); // copy number + core_phys = new G4PVPlacement(0, // no rotation + noTranslation, // place inner tube concentric to outer + core_log, // the logical volume + prefix + "_" + std::to_string(ID) + "_core", // a name for this physical volume + inner_log, // the mother volume + false, // no boolean ops + 0); // copy number // ------------ Vis Attributes ------------- G4VisAttributes *vis = new G4VisAttributes(); @@ -98,17 +97,15 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi G4VisAttributes *vis_core = new G4VisAttributes(); try { const std::vector &color = myTable->GetDArray("color"); - if (color.size() == 3) { // RGB + if (color.size() == 3) { // RGB vis->SetColour(G4Colour(color[0], color[1], color[2])); vis_inner->SetColour(G4Colour(color[2], color[0], color[1])); vis_core->SetColour(G4Colour(color[1], color[2], color[0])); - } - else if (color.size() == 4) { // RGBA + } else if (color.size() == 4) { // RGBA vis->SetColour(G4Colour(color[0], color[1], color[2], color[3])); vis_inner->SetColour(G4Colour(color[2], color[0], color[1], color[3])); vis_core->SetColour(G4Colour(color[1], color[2], color[0], color[3])); - } - else + } else warn << "GeoNestedTubeConstruction error: " << myTable->GetName() << "[" << myTable->GetIndex() << "].color must have 3 or 4 components" << newline; } catch (DBNotFoundError &e) { @@ -119,13 +116,11 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi vis->SetForceWireframe(true); vis_inner->SetForceWireframe(true); vis_core->SetForceWireframe(true); - } - else if (drawstyle == "solid") { + } else if (drawstyle == "solid") { vis->SetForceSolid(true); vis_inner->SetForceSolid(true); vis_core->SetForceSolid(true); - } - else + } else warn << "GeoNestedTubeConstruction error: " << myTable->GetName() << "[" << myTable->GetIndex() << "].drawstyle must be either \"wireframe\" or \"solid\"."; } catch (DBNotFoundError &e) { @@ -134,15 +129,14 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi // Check for invisible flag last try { int invisible = myTable->GetI("invisible"); - if (invisible) - { + if (invisible) { outer_log->SetVisAttributes(G4VisAttributes::GetInvisible()); inner_log->SetVisAttributes(G4VisAttributes::GetInvisible()); core_log->SetVisAttributes(G4VisAttributes::GetInvisible()); } } catch (DBNotFoundError &e) { }; - + outer_log->SetVisAttributes(vis); inner_log->SetVisAttributes(vis_inner); core_log->SetVisAttributes(vis_core); @@ -157,9 +151,10 @@ G4VSolid *GeoNestedTubeConstruction::BuildSolid(const std::string &name) { return outer; } -G4PVPlacement *GeoNestedTubeConstruction::PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, - G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, - bool booleanSolid, int copyNo) { +G4PVPlacement *GeoNestedTubeConstruction::PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, + const std::string &name, G4LogicalVolume *logi_tube, + G4VPhysicalVolume *mother_phys, bool booleanSolid, + int copyNo) { G4PVPlacement *outer_phys = new G4PVPlacement(tuberot, tubepos, name, logi_tube, mother_phys, booleanSolid, copyNo); // core surface From e7e5500fe6325877d71c4b154504cba59736ca38 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Mon, 24 Jun 2024 11:29:54 +0100 Subject: [PATCH 05/21] Correct documentation Changed "optical photon" back to "opticalphoton" to reflect the name of the G4 class it is referring to. --- doc/users_guide/processes.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/users_guide/processes.rst b/doc/users_guide/processes.rst index d3491651..0faa103e 100644 --- a/doc/users_guide/processes.rst +++ b/doc/users_guide/processes.rst @@ -128,7 +128,7 @@ below. Wavelength Shifting ``````````````````` There are a few ways of doing bulk wavelength shifting in RAT. The default -behavior is for GLG4Scint to handle optical photons as well as charged +behavior is for GLG4Scint to handle opticalphotons as well as charged particles. Alternatively, you can also let GLG4Scint handle the primary scintillation, then use Geant4's G4OpWLS process or the custom BNLOpWLSModel to do the reemission. From 5ff4aea9800bd00cb5fd2a0f190c2d686689eae6 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Mon, 24 Jun 2024 15:09:41 +0100 Subject: [PATCH 06/21] Correct clang-format Formatting was flagged, despite using clang-format on all edited files (I made sure the .clang-format file was being used). I've changed the flagged files to have the same format as similar existing files. --- src/ds/include/RAT/DS/NestedTubeInfo.hh | 4 ++-- src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ds/include/RAT/DS/NestedTubeInfo.hh b/src/ds/include/RAT/DS/NestedTubeInfo.hh index 24e372e0..e2ee6b61 100644 --- a/src/ds/include/RAT/DS/NestedTubeInfo.hh +++ b/src/ds/include/RAT/DS/NestedTubeInfo.hh @@ -18,8 +18,8 @@ namespace DS { class NestedTubeInfo : public TObject { public: - NestedTubeInfo() : TObject() {}; - virtual ~NestedTubeInfo() = default; + NestedTubeInfo() : TObject() {} + virtual ~NestedTubeInfo() {} virtual void AddNestedTube(const G4ThreeVector& _pos, const G4ThreeVector& _dir, const double _length, const double _core_r, const double _inner_r, const double _outer_r, diff --git a/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh index babe0be7..c7cec1fe 100644 --- a/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh +++ b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh @@ -7,7 +7,7 @@ namespace RAT { class GeoNestedSolidArrayFactoryBase : public GeoFactory { public: - GeoNestedSolidArrayFactoryBase(const std::string &name) : GeoFactory(name) {}; + GeoNestedSolidArrayFactoryBase(const std::string &name) : GeoFactory(name){}; static const DS::NestedTubeInfo &GetNestedTubeInfo() { return nestedtubeinfo; } protected: From 4ef9a593dd2f3b46881eabead90c45fb9852e929 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Wed, 3 Jul 2024 17:07:05 +0100 Subject: [PATCH 07/21] Add sensitive detector capability for nested tubes One can now set the cores of nested tubes as sensitive detectors. Useful if using as optical fibres. --- src/geo/CMakeLists.txt | 2 + .../include/RAT/GeoFiberSensitiveDetector.hh | 54 +++++ .../RAT/GeoFiberSensitiveDetectorHit.hh | 63 +++++ .../include/RAT/GeoNestedTubeConstruction.hh | 88 +++---- src/geo/src/GeoFiberSensitiveDetector.cc | 219 ++++++++++++++++++ src/geo/src/GeoFiberSensitiveDetectorHit.cc | 58 +++++ src/geo/src/GeoNestedSolidArrayFactoryBase.cc | 8 +- src/geo/src/GeoNestedTubeConstruction.cc | 23 +- 8 files changed, 463 insertions(+), 52 deletions(-) create mode 100644 src/geo/include/RAT/GeoFiberSensitiveDetector.hh create mode 100644 src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh create mode 100644 src/geo/src/GeoFiberSensitiveDetector.cc create mode 100644 src/geo/src/GeoFiberSensitiveDetectorHit.cc diff --git a/src/geo/CMakeLists.txt b/src/geo/CMakeLists.txt index 65aec4d5..3ff72c9f 100644 --- a/src/geo/CMakeLists.txt +++ b/src/geo/CMakeLists.txt @@ -18,6 +18,8 @@ add_library(geo OBJECT src/GeoCalibrationStickFactory.cc src/GeoCherenkovSourceFactory.cc src/GeoConvexLensFactory.cc + src/GeoFiberSensitiveDetector.cc + src/GeoFiberSensitiveDetectorHit.cc src/GeoCutTubeFactory.cc src/GeoEosFactory.cc src/GeoFactory.cc diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetector.hh b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh new file mode 100644 index 00000000..b4066d2d --- /dev/null +++ b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh @@ -0,0 +1,54 @@ +#ifndef __RAT_GeoFiberSensitiveDetector__ +#define __RAT_GeoFiberSensitiveDetector__ + +#include +#include + +class G4Step; +class G4HCofThisEvent; +class G4TouchableHistory; + +namespace RAT { + +class GeoFiberSensitiveDetector : public G4VSensitiveDetector +{ + public: + GeoFiberSensitiveDetector(G4String name); + virtual ~GeoFiberSensitiveDetector(); + + virtual void Initialize(G4HCofThisEvent *HCE); + virtual G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist); + virtual void EndOfEvent(G4HCofThisEvent *HCE); + + // Data members which are publicly accessible and can be + // written out to the RAT event tree + + std::vector _hit_x; + /** hit x-coordinate */ + std::vector _hit_y; + /** hit y-coordinate */ + std::vector _hit_z; + /** hit z-coordinate */ + std::vector _hit_E; + /** hit energy deposition */ + std::vector _hit_time; + /** global time of the hit */ + std::vector _hit_uid; + /** unique identifier code of the hit */ + std::vector _hit_pdg; + /** pdg of particle that left the hit */ + std::vector _hit_volume; + /** name of volume of hit */ + + private: + int fLastEventID; + int fLastTrackID; + + GeoFiberSensitiveDetectorHitsCollection *_hitsCollection; + G4int HCID; + G4HCofThisEvent *_HCE; +}; + +} // namespace RAT + +#endif diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh new file mode 100644 index 00000000..78478da1 --- /dev/null +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -0,0 +1,63 @@ +#ifndef __RAT_GeoFiberSensitiveDetectorHit__ +#define __RAT_GeoFiberSensitiveDetectorHit__ + +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +class GeoFiberSensitiveDetectorHit : public G4VHit { + public: + GeoFiberSensitiveDetectorHit(G4int i, G4double t); + virtual ~GeoFiberSensitiveDetectorHit(); + GeoFiberSensitiveDetectorHit(const GeoFiberSensitiveDetectorHit &right); + const GeoFiberSensitiveDetectorHit &operator=(const GeoFiberSensitiveDetectorHit &right); + int operator==(const GeoFiberSensitiveDetectorHit &right) const; + + inline void *operator new(size_t); + inline void operator delete(void *aHit); + + void Draw(); + void Print(); + + private: + G4int id; + G4double time; + G4ThreeVector pos; + G4RotationMatrix rot; + const G4LogicalVolume *pLogV; + + public: + inline G4int GetID() const { return id; } + inline G4double GetTime() const { return time; } + inline void SetTime(G4double val) { time = val; } + inline void SetPos(G4ThreeVector xyz) { pos = xyz; } + inline G4ThreeVector GetPos() const { return pos; } + inline void SetRot(G4RotationMatrix rmat) { rot = rmat; } + inline G4RotationMatrix GetRot() const { return rot; } + inline void SetLogV(G4LogicalVolume *val) { pLogV = val; } + inline const G4LogicalVolume *GetLogV() const { return pLogV; } +}; + +typedef G4THitsCollection GeoFiberSensitiveDetectorHitsCollection; + +extern G4Allocator GeoFiberSensitiveDetectorHitAllocator; + +inline void *GeoFiberSensitiveDetectorHit::operator new(size_t) { + void *aHit; + aHit = (void *)GeoFiberSensitiveDetectorHitAllocator.MallocSingle(); + return aHit; +} + +inline void GeoFiberSensitiveDetectorHit::operator delete(void *aHit) { + GeoFiberSensitiveDetectorHitAllocator.FreeSingle((GeoFiberSensitiveDetectorHit *)aHit); +} + +} // namespace RAT + +#endif diff --git a/src/geo/include/RAT/GeoNestedTubeConstruction.hh b/src/geo/include/RAT/GeoNestedTubeConstruction.hh index a1b5401d..50112f85 100644 --- a/src/geo/include/RAT/GeoNestedTubeConstruction.hh +++ b/src/geo/include/RAT/GeoNestedTubeConstruction.hh @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -15,46 +16,47 @@ namespace RAT { -struct GeoNestedTubeConstructionParams { - GeoNestedTubeConstructionParams() { invisible = false; }; - - bool invisible; - - double outer_r; - double inner_r; - double core_r; - double Dz; // half length - - G4Material *outer; - G4Material *inner; - G4Material *core; - - // G4OpticalSurface *outer_inner; - G4OpticalSurface *inner_core; -}; - -class GeoNestedTubeConstruction { - public: - GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID); - virtual ~GeoNestedTubeConstruction() {} - - virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID); - virtual G4VSolid *BuildSolid(const std::string &prefix); - virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, - G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, - int copyNo); - - protected: - // physical volumes - G4PVPlacement *inner_phys; - G4PVPlacement *core_phys; - - G4LogicalVolume *log_tube; - GeoNestedTubeConstructionParams fParams; - - DBLinkPtr myTable; -}; - -} // namespace RAT - -#endif \ No newline at end of file + struct GeoNestedTubeConstructionParams { + GeoNestedTubeConstructionParams() { invisible = false; }; + + bool invisible; + + double outer_r; + double inner_r; + double core_r; + double Dz; // half length + + G4Material *outer; + G4Material *inner; + G4Material *core; + + // G4OpticalSurface *outer_inner; + G4OpticalSurface *inner_core; + }; + + class GeoNestedTubeConstruction { + public: + GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID); + virtual ~GeoNestedTubeConstruction() {} + + virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID, DBLinkPtr table); + virtual G4VSolid *BuildSolid(const std::string &prefix); + virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, + G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, + int copyNo); + + protected: + // physical volumes + G4PVPlacement *inner_phys; + G4PVPlacement *core_phys; + + G4LogicalVolume *log_tube; + GeoNestedTubeConstructionParams fParams; + + DBLinkPtr myTable; + }; + +} // namespace RAT + +#endif + \ No newline at end of file diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc new file mode 100644 index 00000000..ed1d5e8c --- /dev/null +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -0,0 +1,219 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +GeoFiberSensitiveDetector::GeoFiberSensitiveDetector(G4String name) : G4VSensitiveDetector(name) { + G4String HCname; + collectionName.insert(HCname = "CustomSenDet"); + HCID = -1; +} + +GeoFiberSensitiveDetector::~GeoFiberSensitiveDetector() { ; } + +void GeoFiberSensitiveDetector::Initialize(G4HCofThisEvent *HCE) { + debug << "GeoFiberSensitiveDetector::Initialize start." << newline; + _hitsCollection = new GeoFiberSensitiveDetectorHitsCollection(SensitiveDetectorName, collectionName[0]); + + debug << "GeoFiberSensitiveDetector::Initialize hit collection address is " << _hitsCollection << newline; + if (HCID < 0) { + HCID = G4SDManager::GetSDMpointer()->GetCollectionID(_hitsCollection); + } + debug << "GeoFiberSensitiveDetector::Initialize hit collection ID = " << HCID << newline; + HCE->AddHitsCollection(HCID, _hitsCollection); + + // store pointer to hit collection + _HCE = HCE; + + // + // Initialize the data members used to store information for the RAT Event + // + + // Empty the hit information std::vectors + _hit_x.clear(); + _hit_y.clear(); + _hit_z.clear(); + _hit_E.clear(); + _hit_time.clear(); + _hit_uid.clear(); + _hit_pdg.clear(); + _hit_volume.clear(); + fLastTrackID = fLastEventID = -1; + + debug << "GeoFiberSensitiveDetector::Initialize end." << newline; +} + +G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory * /*ROhist*/) { + // We ONLY want to store optical photons + if (aStep->GetTrack()->GetDefinition()->GetParticleName() == "opticalphoton") { + + debug << "GeoFiberSensitiveDetector::ProcessHits start." << newline; + debug << "GeoFiberSensitiveDetector::ProcessHits getting energy deposited." << newline; + + G4double edep = aStep->GetTotalEnergyDeposit(); + G4double dl = aStep->GetStepLength(); + + debug << " Energy deposited: " << edep << newline; + + // if(edep==0.) return true; + + G4Track *aTrack = aStep->GetTrack(); + G4ParticleDefinition *part = aTrack->GetDefinition(); + int pdg = part->GetPDGEncoding(); + + debug << " track information: " << newline << " G4Track Pointer: " << aTrack << newline + << " Particle Definition Pointer: " << part << newline << " Particle PDG Encoding: " << pdg + << newline; + + debug << "GeoFiberSensitiveDetector::ProcessHits getting global time." << newline; + + G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); + G4Material *m = preStepPoint->GetMaterial(); + G4String mname = m->GetName(); + G4double d = m->GetDensity(); + G4String f = m->GetChemicalFormula(); + + G4TouchableHistory *theTouchable = (G4TouchableHistory *)(preStepPoint->GetTouchable()); + // Get the copy number of this element. This may not be unique, + // if the mother is also a replica and occurs more than once. + G4int copyNo = theTouchable->GetVolume()->GetCopyNo(); + // Get the mother copy number, and offset it by 1000. This will help + // form the basis of the chamber ID. + G4int motherCopyNo = (theTouchable->GetVolume(1)->GetCopyNo()) + 1000; + G4int uid = motherCopyNo + copyNo; // unique identifier of chamber + + // construct a unique identifier for this copy + G4int ivol = 0; + G4int idOffset = 1; + uid = 0; + + debug << "History level: " << theTouchable->GetHistoryDepth() << newline; + + while (ivol < theTouchable->GetHistoryDepth()) { + debug << " * volume layer level = " << ivol << newline; + uid += theTouchable->GetVolume(ivol)->GetCopyNo() * idOffset; + idOffset *= 100; + ivol++; + } + + G4ThreeVector worldPos = preStepPoint->GetPosition(); + + debug << "Hit material name " << mname << newline; + debug << "density " << d << newline; + debug << "formula " << f << newline; + debug << "edep " << G4BestUnit(edep, "Energy") << newline; + debug << "dl " << G4BestUnit(dl, "Length") << newline; + debug << "pid " << pdg << newline; + debug << "Position " << G4BestUnit(worldPos.x(), "Length") << " " << G4BestUnit(worldPos.y(), "Length") << " " + << G4BestUnit(worldPos.z(), "Length") << " " << newline; + + debug << " " << newline; + + G4double hitTime = preStepPoint->GetGlobalTime(); + + debug << "GeoFiberSensitiveDetector::ProcessHits checking for an existing hit " + "in this element." + << newline; + // check if this finger already has a hit + G4int ix = -1; + + debug << "GeoFiberSensitiveDetector::ProcessHits hit collection address is " << _hitsCollection << newline; + G4cerr << "GeoFiberSensitiveDetector: Hit ID = " << uid << " and position: " << worldPos << newline; + + int eventID = G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent(); + int trackID = aStep->GetTrack()->GetTrackID(); + + // get volume name of hit as well + // note that this volume name may not be unique + G4String vol = theTouchable->GetVolume()->GetName(); + + if (fLastEventID != eventID || fLastTrackID != trackID) { + // Fill the hit information + _hit_x.push_back(worldPos.x()); + _hit_y.push_back(worldPos.y()); + _hit_z.push_back(worldPos.z()); + _hit_E.push_back(edep); + _hit_time.push_back(hitTime); + _hit_uid.push_back(uid); + _hit_pdg.push_back(pdg); + _hit_volume.push_back(vol); + fLastEventID = eventID; + fLastTrackID = trackID; + } + + /*G4int hitfreq = G4int(db["veto_hit_frequency"]);*/ + + if (NULL == _hitsCollection) { + debug << "GeoFiberSensitiveDetector::ProcessHits hit collection null. " + "Reloading from HCofEThisEvent." + << newline; + if (_HCE) { + _hitsCollection = (GeoFiberSensitiveDetectorHitsCollection *)(_HCE->GetHC(HCID)); + debug << "GeoFiberSensitiveDetector::ProcessHits * hit collection address is " << _hitsCollection << newline; + } else { + debug << "GeoFiberSensitiveDetector::ProcessHits (E) HCofEThisEvent " + "pointer is NULL!" + << newline; + } + } + + if (_hitsCollection) { + for (G4int i = 0; i < _hitsCollection->entries(); i++) { + // RAT::debug << " * GeoFiberSensitiveDetector::ProcessHits checking hit " + // << i + 1 + // << " of " + // << _hitsCollection->entries() + // << newline; + debug << " * this hit ID is " << (*_hitsCollection)[i]->GetID() << newline; + if ((*_hitsCollection)[i]->GetID() == uid) { + ix = i; + break; + } + } + + // if it has, then take the earlier time + if (ix >= 0) { + debug << "GeoFiberSensitiveDetector::ProcessHits use existing earlier time " + "for hit." + << newline; + if ((*_hitsCollection)[ix]->GetTime() > hitTime) { + (*_hitsCollection)[ix]->SetTime(hitTime); + } + } else + // if not, create a new hit and std::set it to the collection + { + debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; + GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime); + G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); + aHit->SetLogV(thePhysical->GetLogicalVolume()); + G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); + aTrans.Invert(); + aHit->SetRot(aTrans.NetRotation()); + aHit->SetPos(aTrans.NetTranslation()); + _hitsCollection->insert(aHit); + aHit->Print(); + aHit->Draw(); + debug << " * Drawing Hit " << uid << newline; + } + } + debug << "GeoFiberSensitiveDetector::ProcessHits end." << newline; + return true; + } + return true; +} + +void GeoFiberSensitiveDetector::EndOfEvent(G4HCofThisEvent * /*HCE*/) { ; } + +} // namespace RAT diff --git a/src/geo/src/GeoFiberSensitiveDetectorHit.cc b/src/geo/src/GeoFiberSensitiveDetectorHit.cc new file mode 100644 index 00000000..b78ca7fd --- /dev/null +++ b/src/geo/src/GeoFiberSensitiveDetectorHit.cc @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +G4Allocator GeoFiberSensitiveDetectorHitAllocator; + +GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(G4int i, G4double t) { + id = i; + time = t; + pLogV = 0; +} +GeoFiberSensitiveDetectorHit::~GeoFiberSensitiveDetectorHit() { ; } + +GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(const GeoFiberSensitiveDetectorHit &right) : G4VHit() { + id = right.id; + time = right.time; + pos = right.pos; + rot = right.rot; + pLogV = right.pLogV; +} + +const GeoFiberSensitiveDetectorHit &GeoFiberSensitiveDetectorHit::operator=(const GeoFiberSensitiveDetectorHit &right) { + id = right.id; + time = right.time; + pos = right.pos; + rot = right.rot; + pLogV = right.pLogV; + return *this; +} + +int GeoFiberSensitiveDetectorHit::operator==(const GeoFiberSensitiveDetectorHit & /*right*/) const { return 0; } + +void GeoFiberSensitiveDetectorHit::Draw() { + G4VVisManager *pVVisManager = G4VVisManager::GetConcreteInstance(); + if (pVVisManager) { + G4Transform3D trans(rot.inverse(), pos); + G4VisAttributes attribs; + const G4VisAttributes *pVA = pLogV->GetVisAttributes(); + if (pVA) attribs = *pVA; + G4Colour colour(0., 1., 1.); + attribs.SetColour(colour); + attribs.SetForceSolid(true); + pVVisManager->Draw(*pLogV, attribs, trans); + } +} + +void GeoFiberSensitiveDetectorHit::Print() { + debug << " GeoFiberSensitiveDetector[" << id << "] " << time / CLHEP::ns << " (nsec)" << newline; +} + +} // namespace RAT diff --git a/src/geo/src/GeoNestedSolidArrayFactoryBase.cc b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc index 0f633f75..09a4a9d7 100644 --- a/src/geo/src/GeoNestedSolidArrayFactoryBase.cc +++ b/src/geo/src/GeoNestedSolidArrayFactoryBase.cc @@ -131,7 +131,7 @@ G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { if ((sub_type == -1) || (sub_type == sub_type_array[solidID])) { // construct GeoNestedTubeConstruction *construction = new GeoNestedTubeConstruction(table, lpos_table, mother, solidID); - G4LogicalVolume *log_tube = construction->BuildVolume(volume_name, solidID); + G4LogicalVolume *log_tube = construction->BuildVolume(volume_name, solidID, table); // name std::string tubename = volume_name + "_" + ::to_string(solidID); @@ -177,9 +177,9 @@ G4VPhysicalVolume *GeoNestedSolidArrayFactoryBase::Construct(DBLinkPtr table) { double core_r = table->GetD("core_r"); double inner_r = table->GetD("inner_r"); double outer_r = table->GetD("outer_r"); - std::string core_material = table->GetS("material_core"); - std::string inner_material = table->GetS("material_inner"); - std::string outer_material = table->GetS("material_outer"); + std::string core_material = table->GetS("core_material"); + std::string inner_material = table->GetS("inner_material"); + std::string outer_material = table->GetS("outer_material"); nestedtubeinfo.AddNestedTube(tubepos, soliddir, length, core_r, inner_r, outer_r, core_material, inner_material, outer_material); diff --git a/src/geo/src/GeoNestedTubeConstruction.cc b/src/geo/src/GeoNestedTubeConstruction.cc index a88142d4..31bab4eb 100644 --- a/src/geo/src/GeoNestedTubeConstruction.cc +++ b/src/geo/src/GeoNestedTubeConstruction.cc @@ -2,9 +2,12 @@ #include #include #include +#include #include #include #include +#include +#include #include #include #include @@ -26,11 +29,11 @@ GeoNestedTubeConstruction::GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr fParams.core_r = table->GetD("core_r"); // Materials - fParams.outer = G4Material::GetMaterial(table->GetS("material_outer")); - fParams.inner = G4Material::GetMaterial(table->GetS("material_inner")); - fParams.core = G4Material::GetMaterial(table->GetS("material_core")); + fParams.outer = G4Material::GetMaterial(table->GetS("outer_material")); + fParams.inner = G4Material::GetMaterial(table->GetS("inner_material")); + fParams.core = G4Material::GetMaterial(table->GetS("core_material")); DB *db = DB::Get(); - std::string core_surface_name = table->GetS("material_core"); + std::string core_surface_name = table->GetS("core_material"); if (Materials::optical_surface.count(core_surface_name) == 0) Log::Die("GeoSolidFactory: Surface " + core_surface_name + " does not exist"); fParams.inner_core = Materials::optical_surface[core_surface_name]; @@ -50,7 +53,7 @@ GeoNestedTubeConstruction::GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr Log::Assert(fParams.inner_core, "GeoNestedTubeConstruction: " + tube_name + " has an invalid core surface material"); } -G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefix, int ID) { +G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefix, int ID, DBLinkPtr table) { if (log_tube) { return log_tube; } @@ -70,6 +73,16 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi outer_log = new G4LogicalVolume(outer_solid, fParams.outer, prefix + "_outer_logic"); inner_log = new G4LogicalVolume(inner_solid, fParams.inner, prefix + "_inner_logic"); core_log = new G4LogicalVolume(core_solid, fParams.core, prefix + "_core_logic"); + // set as sensitive detector if applicable + try { + std::string sensitive_detector = table->GetS("sensitive_detector"); + GeoFiberSensitiveDetector* sensitive = new GeoFiberSensitiveDetector(sensitive_detector+std::to_string(ID)); + G4SDManager *sdman = G4SDManager::GetSDMpointer(); + sdman->AddNewDetector(sensitive); + core_log->SetSensitiveDetector(sensitive); + } + catch (DBNotFoundError &e) { + } // ------------ Physical Volumes ------------- G4ThreeVector noTranslation(0., 0., 0.); From 03191e36ac7a9a463eae6caf500c88b43f06e08e Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 4 Jul 2024 16:53:33 +0100 Subject: [PATCH 08/21] Add to Rat DS so fiber hit info can be saved --- src/core/include/RAT/Gsim.hh | 2 + src/core/src/Gsim.cc | 36 ++++++++++++++++ src/ds/CMakeLists.txt | 2 + src/ds/include/RAT/DS/LinkDef.hh | 2 + src/ds/include/RAT/DS/MC.hh | 11 +++++ src/ds/include/RAT/DS/MCNestedTube.hh | 51 ++++++++++++++++++++++ src/ds/include/RAT/DS/MCNestedTubeHit.hh | 54 ++++++++++++++++++++++++ src/ds/include/RAT/DS/MCPhoton.hh | 2 +- src/ds/src/Root.cc | 3 ++ src/geo/src/GeoFiberSensitiveDetector.cc | 4 +- 10 files changed, 164 insertions(+), 3 deletions(-) create mode 100644 src/ds/include/RAT/DS/MCNestedTube.hh create mode 100644 src/ds/include/RAT/DS/MCNestedTubeHit.hh diff --git a/src/core/include/RAT/Gsim.hh b/src/core/include/RAT/Gsim.hh index 9a23c172..a5100327 100644 --- a/src/core/include/RAT/Gsim.hh +++ b/src/core/include/RAT/Gsim.hh @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -74,6 +75,7 @@ class Gsim : public Producer, G4UserRunAction, G4UserEventAction, G4UserTracking void Init(); // the real constructor void AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventInfo *exinfo = NULL, std::string process = "unknown"); + void AddMCNestedTubeHit(DS::MCNestedTube *rat_mcnt, const GeoFiberSensitiveDetectorHit *hit); /* Storing optical creation track ID and step */ void PhotonRecurse(std::vector &PhotonIDs, int trackID, int &parentID, int &firstCreatedID); diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index 188a024a..da757bbb 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -22,6 +22,8 @@ #include #include #include +#include +#include #include #include #include @@ -567,6 +569,25 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) { } } mc->SetNumPE(numPE); + + /** hits from fibers */ + G4HCofThisEvent *HC = g4ev->GetHCofThisEvent(); + for (int hc = 0; hc < HC->GetNumberOfCollections(); hc++) { + GeoFiberSensitiveDetectorHitsCollection *hit_collection = (GeoFiberSensitiveDetectorHitsCollection*)HC->GetHC(hc); + if (hit_collection->GetName() != "FiberSenDet") + continue; + DS::MCNestedTube *rat_mcnt = mc->AddNewMCNestedTube(); + G4String det_name = hit_collection->GetSDname(); + std::string fibre_id_str = det_name.erase(0, 6).data(); + int fibre_id = std::stoi(fibre_id_str); + rat_mcnt->SetID(fibre_id); + // only process fibers + // info << hit_collection->GetSDname() << newline; + for (int hit = 0; hit < hit_collection->GetSize(); hit++) { + GeoFiberSensitiveDetectorHit *my_hit = (GeoFiberSensitiveDetectorHit*)hit_collection->GetHit(hit); + AddMCNestedTubeHit(rat_mcnt, my_hit); + } + } } void Gsim::AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventInfo * /*exinfo*/, std::string process) { @@ -594,6 +615,21 @@ void Gsim::AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventI rat_mcphoton->SetCreatorProcess(process); } +void Gsim::AddMCNestedTubeHit(DS::MCNestedTube *rat_mcnt, const GeoFiberSensitiveDetectorHit *hit) { + DS::MCNestedTubeHit *rat_mchit = rat_mcnt->AddNewMCNestedTubeHit(); + // Only real photons are added in Gsim, noise and afterpulsing handled in processors + + double x, y, z; + G4ThreeVector pos_vec = hit->GetPos(); + x = pos_vec.x(); + y = pos_vec.y(); + z = pos_vec.z(); + rat_mchit->SetPosition(TVector3(x, y, z)); + + rat_mchit->SetHitID(hit->GetID()); + rat_mchit->SetHitTime(hit->GetTime()); +} + void Gsim::SetStoreParticleTraj(const G4String &particleName, const bool &gDoStore) { if (gDoStore) { fStoreParticleTraj.insert(particleName); diff --git a/src/ds/CMakeLists.txt b/src/ds/CMakeLists.txt index afa0b86c..26e5521e 100644 --- a/src/ds/CMakeLists.txt +++ b/src/ds/CMakeLists.txt @@ -32,6 +32,8 @@ root_generate_dictionary(G__RATDict RAT/DS/MCParticle.hh RAT/DS/MCPhoton.hh RAT/DS/MCPMT.hh + RAT/DS/MCNestedTube.hh + RAT/DS/MCNestedTubeHit.hh RAT/DS/MCSummary.hh RAT/DS/PMT.hh RAT/DS/RunStore.hh diff --git a/src/ds/include/RAT/DS/LinkDef.hh b/src/ds/include/RAT/DS/LinkDef.hh index adea4fb9..17821562 100644 --- a/src/ds/include/RAT/DS/LinkDef.hh +++ b/src/ds/include/RAT/DS/LinkDef.hh @@ -9,7 +9,9 @@ #pragma link C++ class RAT::DS::MCTrack + ; #pragma link C++ class RAT::DS::MCTrackStep + ; #pragma link C++ class RAT::DS::MCPMT + ; +#pragma link C++ class RAT::DS::MCNestedTube + ; #pragma link C++ class RAT::DS::MCPhoton + ; +#pragma link C++ class RAT::DS::MCNestedTubeHit + ; #pragma link C++ class RAT::DS::MCSummary + ; #pragma link C++ class RAT::DS::Calib + ; #pragma link C++ class RAT::DS::FitResult + ; diff --git a/src/ds/include/RAT/DS/MC.hh b/src/ds/include/RAT/DS/MC.hh index a894b52b..b40b69ac 100644 --- a/src/ds/include/RAT/DS/MC.hh +++ b/src/ds/include/RAT/DS/MC.hh @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -102,6 +103,15 @@ class MC : public TObject { }; virtual void PrunePMT() { pmt.resize(0); } + /** List of optical fibers which captured at least one photon */ + virtual MCNestedTube *GetMCNestedTube(int i) { return &nt[i]; } + virtual int GetMCNestedTubeCount() const { return nt.size(); } + virtual MCNestedTube *AddNewMCNestedTube() { + nt.resize(nt.size() + 1); + return &nt.back(); + }; + virtual void PruneNestedTube() { nt.resize(0); } + /** Total number of photoelectrons generated in this event */ virtual int GetNumPE() const { return numPE; } virtual void SetNumPE(int _numPE) { numPE = _numPE; } @@ -132,6 +142,7 @@ class MC : public TObject { std::vector parent; std::vector track; std::vector pmt; + std::vector nt; }; } // namespace DS diff --git a/src/ds/include/RAT/DS/MCNestedTube.hh b/src/ds/include/RAT/DS/MCNestedTube.hh new file mode 100644 index 00000000..d738b0dd --- /dev/null +++ b/src/ds/include/RAT/DS/MCNestedTube.hh @@ -0,0 +1,51 @@ +/** + * @class DS::MCNestedTube + * Data Structure: Hit NestedTube in Monte Carlo + * + * @author Wilf Shorrock + * + * This class represents a NestedTube in which at least one photon was captured + */ + +#ifndef __RAT_DS_MCNestedTube__ +#define __RAT_DS_MCNestedTube__ + +#include +#include +#include + +namespace RAT { +namespace DS { + +class MCNestedTube : public TObject { + public: + MCNestedTube() : TObject() {} + virtual ~MCNestedTube() {} + + /** ID number */ + virtual Int_t GetID() const { return id; }; + virtual void SetID(Int_t _id) { id = _id; }; + + /** List of photons captured in this NestedTube. */ + MCNestedTubeHit *GetMCNestedTubeHit(Int_t i) { return &photon[i]; } + Int_t GetMCNestedTubeHitCount() const { return photon.size(); } + MCNestedTubeHit *AddNewMCNestedTubeHit() { + photon.resize(photon.size() + 1); + return &photon.back(); + } + void RemoveMCNestedTubeHit(Int_t i) { photon.erase(photon.begin() + i); } + void PruneMCNestedTubeHit() { photon.resize(0); } + void SortMCNestedTubeHits() { std::sort(photon.begin(), photon.end()); } + + ClassDef(MCNestedTube, 3); + + protected: + Int_t id; + Int_t type; + std::vector photon; +}; + +} // namespace DS +} // namespace RAT + +#endif diff --git a/src/ds/include/RAT/DS/MCNestedTubeHit.hh b/src/ds/include/RAT/DS/MCNestedTubeHit.hh new file mode 100644 index 00000000..039a774d --- /dev/null +++ b/src/ds/include/RAT/DS/MCNestedTubeHit.hh @@ -0,0 +1,54 @@ +/** + * @class DS::MCNestedTubeHit + * Data Structure: Photon captured in a nested tube (optical fiber). + * + * @author Wilf Shorrock + * + * This class represents a single photon that is captured within the core of a + * nested tube object. + */ + +#ifndef __RAT_DS_MCNestedTubeHit__ +#define __RAT_DS_MCNestedTubeHit__ + +#include +#include + +namespace RAT { +namespace DS { + +class MCNestedTubeHit : public TObject { + public: + MCNestedTubeHit() : TObject() {} + virtual ~MCNestedTubeHit() {} + + /** Time of photon hit on nested tube relative to event start time (ns). */ + virtual Double_t GetHitTime() const { return hitTime; } + virtual void SetHitTime(Double_t _hitTime) { hitTime = _hitTime; } + + /** Location of photon hit (mm). */ + virtual TVector3 GetPosition() const { return pos; } + virtual void SetPosition(const TVector3 &_pos) { pos = _pos; } + + /** Hit ID of photon */ + virtual void SetHitID(Int_t _hitID) { hitID = _hitID; } + virtual Int_t GetHitID() const { return hitID; } + + /** Operator overload **/ + bool operator<(const MCNestedTubeHit &mcp) const { return (hitTime < mcp.hitTime); } + bool operator>(const MCNestedTubeHit &mcp) const { return (hitTime > mcp.hitTime); } + + ClassDef(MCNestedTubeHit, 4); + + protected: + Double_t hitTime; + TVector3 pos; + + Int_t hitID; + Int_t fiberID; +}; + +} // namespace DS +} // namespace RAT + +#endif diff --git a/src/ds/include/RAT/DS/MCPhoton.hh b/src/ds/include/RAT/DS/MCPhoton.hh index 1793dad4..5ea725d1 100644 --- a/src/ds/include/RAT/DS/MCPhoton.hh +++ b/src/ds/include/RAT/DS/MCPhoton.hh @@ -8,7 +8,7 @@ * photocathode of the PMT. The time jitter and delay in transit to the * anode are not included here, but the distribution of charge is, * which is slightly incongruous. - * + * Note that we require that the photon generates a photoelectron. * Absorbed photons are not included here. */ diff --git a/src/ds/src/Root.cc b/src/ds/src/Root.cc index 3cb9708b..66ad8b44 100644 --- a/src/ds/src/Root.cc +++ b/src/ds/src/Root.cc @@ -4,8 +4,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -15,6 +17,7 @@ ClassImp(RAT::DS::Root); ClassImp(RAT::DS::MC); ClassImp(RAT::DS::MCParticle); ClassImp(RAT::DS::MCPMT); +ClassImp(RAT::DS::MCNestedTube); ClassImp(RAT::DS::Calib); ClassImp(RAT::DS::EV); ClassImp(RAT::DS::PMT); diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index ed1d5e8c..1fbaedb5 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -17,7 +17,7 @@ namespace RAT { GeoFiberSensitiveDetector::GeoFiberSensitiveDetector(G4String name) : G4VSensitiveDetector(name) { G4String HCname; - collectionName.insert(HCname = "CustomSenDet"); + collectionName.insert(HCname = "FiberSenDet"); HCID = -1; } @@ -130,7 +130,7 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory G4int ix = -1; debug << "GeoFiberSensitiveDetector::ProcessHits hit collection address is " << _hitsCollection << newline; - G4cerr << "GeoFiberSensitiveDetector: Hit ID = " << uid << " and position: " << worldPos << newline; + // G4cerr << "GeoFiberSensitiveDetector: Hit ID = " << uid << " and position: " << worldPos << newline; int eventID = G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent(); int trackID = aStep->GetTrack()->GetTrackID(); From 228c8bc96cfbcb10a4a017d6e476d879545b5801 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 4 Jul 2024 16:58:24 +0100 Subject: [PATCH 09/21] Update docs for nested tube arrays --- doc/users_guide/geometry.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/users_guide/geometry.rst b/doc/users_guide/geometry.rst index 0be64013..1ac9b799 100644 --- a/doc/users_guide/geometry.rst +++ b/doc/users_guide/geometry.rst @@ -66,7 +66,7 @@ the volume. The common fields shared by all tables: ``mother`` ``string`` Name of the mother volume. The mother volume should fully contain this volume. The world volume has the mother "". ``enable`` ``int`` (optional) If set to zero, this volume is skipped and not constructed. ``type`` ``string`` Shape of this volume, see below for list. -``sensitive_detector`` ``string`` (optional) Name of sensitive detector if this volume should register hits. Limited to ''/mydet/pmt/inner'' and ''/mydet/veto/genericchamber'' +``sensitive_detector`` ``string`` (optional) Name of sensitive detector if this volume should register hits. Limited to ''/mydet/pmt/inner'', ''/mydet/fibers'' and ''/mydet/veto/genericchamber'' ====================== ====================== =================== Allowed types: @@ -168,6 +168,7 @@ NestedTubeArray Fields: ``orientation`` ``string`` Method of determining nested tube direction. "point" will aim all nested tubes at a point in space. "manual" requires that the position table also contain dir_x, dir_y, and dir_z fields which define the direction vector for each PMT. ``orient_point`` ``float[3]`` (optional) Point (mm) in mother volume to aim all tubes toward. ``rescale_radius`` ``float`` (optional) Assumes all tubes are spherically arranged around the center of the mother volume and rescales their positions to a particular radius. By default, no rescaling is done. +``sensitive_detector`` ``string`` (optional) Name of sensitive detector if this volume should register hits. Limited to ''/mydet/fibers''.' ====================== ========================== =================== Creating a parameterized geometry From a39064c2f57d5eb3be96aed800bbda5d132e330f Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Fri, 5 Jul 2024 10:18:54 +0100 Subject: [PATCH 10/21] Update nested tube documentation --- doc/users_guide/geometry.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/users_guide/geometry.rst b/doc/users_guide/geometry.rst index 1ac9b799..ef9ba731 100644 --- a/doc/users_guide/geometry.rst +++ b/doc/users_guide/geometry.rst @@ -157,9 +157,9 @@ NestedTubeArray Fields: ``core_r`` ``float`` The radius of the core tube (mm) ``inner_r`` ``float`` The radius of the inner tube (mm) ``outer_r`` ``float`` The radius of the outer tube (mm) -``material_core`` ``string`` The material of the core tube -``material_inner`` ``string`` The material of the inner tube -``material_outer`` ``string`` The material of the outer tube +``core_material`` ``string`` The material of the core tube +``inner_material`` ``string`` The material of the inner tube +``outer_material`` ``string`` The material of the outer tube ``Dz`` ``float`` Half-height of tube (mm) ``phi_start`` ``float`` (optional) Angle (deg) where tube segment starts. Default is 0.0 ``phi_delta`` ``float`` (optional) Angle span (deg) of tube segment. Default is 360.0 From ab2567dc986caa372a8c2ad57ab07217fcf85018 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Fri, 5 Jul 2024 13:50:03 +0100 Subject: [PATCH 11/21] Change fiber SD to record all photons --- src/geo/src/GeoFiberSensitiveDetector.cc | 52 +++++++----------------- 1 file changed, 14 insertions(+), 38 deletions(-) diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index 1fbaedb5..05730e48 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -139,7 +139,8 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory // note that this volume name may not be unique G4String vol = theTouchable->GetVolume()->GetName(); - if (fLastEventID != eventID || fLastTrackID != trackID) { + // if (fLastEventID != eventID || fLastTrackID != trackID) { + if (fLastTrackID != trackID) { // Fill the hit information _hit_x.push_back(worldPos.x()); _hit_y.push_back(worldPos.y()); @@ -170,43 +171,18 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory } if (_hitsCollection) { - for (G4int i = 0; i < _hitsCollection->entries(); i++) { - // RAT::debug << " * GeoFiberSensitiveDetector::ProcessHits checking hit " - // << i + 1 - // << " of " - // << _hitsCollection->entries() - // << newline; - debug << " * this hit ID is " << (*_hitsCollection)[i]->GetID() << newline; - if ((*_hitsCollection)[i]->GetID() == uid) { - ix = i; - break; - } - } - - // if it has, then take the earlier time - if (ix >= 0) { - debug << "GeoFiberSensitiveDetector::ProcessHits use existing earlier time " - "for hit." - << newline; - if ((*_hitsCollection)[ix]->GetTime() > hitTime) { - (*_hitsCollection)[ix]->SetTime(hitTime); - } - } else - // if not, create a new hit and std::set it to the collection - { - debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; - GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime); - G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); - aHit->SetLogV(thePhysical->GetLogicalVolume()); - G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); - aTrans.Invert(); - aHit->SetRot(aTrans.NetRotation()); - aHit->SetPos(aTrans.NetTranslation()); - _hitsCollection->insert(aHit); - aHit->Print(); - aHit->Draw(); - debug << " * Drawing Hit " << uid << newline; - } + debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; + GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime); + G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); + aHit->SetLogV(thePhysical->GetLogicalVolume()); + G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); + aTrans.Invert(); + aHit->SetRot(aTrans.NetRotation()); + aHit->SetPos(aTrans.NetTranslation()); + _hitsCollection->insert(aHit); + aHit->Print(); + aHit->Draw(); + debug << " * Drawing Hit " << uid << newline; } debug << "GeoFiberSensitiveDetector::ProcessHits end." << newline; return true; From f9d9d3668aa70aa8eae0cb42d51870ba040f0272 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 11 Jul 2024 09:03:20 +0700 Subject: [PATCH 12/21] Remove drawing of fiber hits --- src/geo/src/GeoFiberSensitiveDetector.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index 05730e48..6b26e28c 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -180,8 +180,8 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory aHit->SetRot(aTrans.NetRotation()); aHit->SetPos(aTrans.NetTranslation()); _hitsCollection->insert(aHit); - aHit->Print(); - aHit->Draw(); + // aHit->Print(); + // aHit->Draw(); debug << " * Drawing Hit " << uid << newline; } debug << "GeoFiberSensitiveDetector::ProcessHits end." << newline; From ddf4449736df9b1ac1c75d7f3d869e7dd762a726 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Fri, 12 Jul 2024 04:25:18 +0200 Subject: [PATCH 13/21] Improve efficiency of nested tube arrays Previously, the nested tube array code would save the information for all nested tubes in the array to the rat DS. This could potentially consume a lot of memory. The behaviour has now been changed so that information is only stored for nested tubes with sensitive detector hits. --- src/core/src/Gsim.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index da757bbb..13586c6b 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -574,7 +574,7 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) { G4HCofThisEvent *HC = g4ev->GetHCofThisEvent(); for (int hc = 0; hc < HC->GetNumberOfCollections(); hc++) { GeoFiberSensitiveDetectorHitsCollection *hit_collection = (GeoFiberSensitiveDetectorHitsCollection*)HC->GetHC(hc); - if (hit_collection->GetName() != "FiberSenDet") + if (hit_collection->GetName() != "FiberSenDet" || hit_collection->GetSize() == 0) continue; DS::MCNestedTube *rat_mcnt = mc->AddNewMCNestedTube(); G4String det_name = hit_collection->GetSDname(); From 340bd41b8ef542e3cd76a66657ed97b816a769dc Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 25 Jul 2024 09:50:02 +0100 Subject: [PATCH 14/21] Record attenuation hits only in fiber sensitive detectors --- src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh | 3 +++ src/geo/src/GeoFiberSensitiveDetector.cc | 8 +++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh index 78478da1..5b21fac9 100644 --- a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -31,6 +31,7 @@ class GeoFiberSensitiveDetectorHit : public G4VHit { G4ThreeVector pos; G4RotationMatrix rot; const G4LogicalVolume *pLogV; + std::string proc; public: inline G4int GetID() const { return id; } @@ -42,6 +43,8 @@ class GeoFiberSensitiveDetectorHit : public G4VHit { inline G4RotationMatrix GetRot() const { return rot; } inline void SetLogV(G4LogicalVolume *val) { pLogV = val; } inline const G4LogicalVolume *GetLogV() const { return pLogV; } + inline void SetProcess(std::string *val) { proc = val; } + inline const G4LogicalVolume *GetProcess() const { return proc; } }; typedef G4THitsCollection GeoFiberSensitiveDetectorHitsCollection; diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index 6b26e28c..bef91593 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -59,6 +59,13 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory // We ONLY want to store optical photons if (aStep->GetTrack()->GetDefinition()->GetParticleName() == "opticalphoton") { + // only store information from steps where the photon is attenuated + G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); + G4StepPoint *postStepPoint = aStep->GetPostStepPoint(); + std::string proc = postStepPoint->GetProcessDefinedStep()->GetProcessName(); + if (proc != "Attenuation") + return true; + debug << "GeoFiberSensitiveDetector::ProcessHits start." << newline; debug << "GeoFiberSensitiveDetector::ProcessHits getting energy deposited." << newline; @@ -79,7 +86,6 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory debug << "GeoFiberSensitiveDetector::ProcessHits getting global time." << newline; - G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); G4Material *m = preStepPoint->GetMaterial(); G4String mname = m->GetName(); G4double d = m->GetDensity(); From 2eac0805a4c90bbeaca30140a7b0c607759ab702 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 25 Jul 2024 10:03:29 +0100 Subject: [PATCH 15/21] Remove faulty functions from header --- src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh index 5b21fac9..2e686b7a 100644 --- a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -43,8 +43,6 @@ class GeoFiberSensitiveDetectorHit : public G4VHit { inline G4RotationMatrix GetRot() const { return rot; } inline void SetLogV(G4LogicalVolume *val) { pLogV = val; } inline const G4LogicalVolume *GetLogV() const { return pLogV; } - inline void SetProcess(std::string *val) { proc = val; } - inline const G4LogicalVolume *GetProcess() const { return proc; } }; typedef G4THitsCollection GeoFiberSensitiveDetectorHitsCollection; From 53731794bb44854d0bb423bfcbd4a8432981bed0 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Tue, 10 Dec 2024 10:35:02 +0000 Subject: [PATCH 16/21] Fix visual display issue The position of the logical volume for each nested tube was set to the position of detected photons in the sensitive detector, which looked strange in the visualiser. Fixed the issue by adding a new hit_position variable to the hit collection and using that instead of the logical volume position. --- src/core/src/Gsim.cc | 2 +- src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh | 3 +++ src/geo/src/GeoFiberSensitiveDetector.cc | 3 ++- src/geo/src/GeoFiberSensitiveDetectorHit.cc | 5 ++++- 4 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index 13586c6b..fbee5284 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -620,7 +620,7 @@ void Gsim::AddMCNestedTubeHit(DS::MCNestedTube *rat_mcnt, const GeoFiberSensitiv // Only real photons are added in Gsim, noise and afterpulsing handled in processors double x, y, z; - G4ThreeVector pos_vec = hit->GetPos(); + G4ThreeVector pos_vec = hit->GetHitPos(); x = pos_vec.x(); y = pos_vec.y(); z = pos_vec.z(); diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh index 2e686b7a..f51268af 100644 --- a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -29,6 +29,7 @@ class GeoFiberSensitiveDetectorHit : public G4VHit { G4int id; G4double time; G4ThreeVector pos; + G4ThreeVector hit_pos; G4RotationMatrix rot; const G4LogicalVolume *pLogV; std::string proc; @@ -39,6 +40,8 @@ class GeoFiberSensitiveDetectorHit : public G4VHit { inline void SetTime(G4double val) { time = val; } inline void SetPos(G4ThreeVector xyz) { pos = xyz; } inline G4ThreeVector GetPos() const { return pos; } + inline void SetHitPos(G4ThreeVector xyz) { hit_pos = xyz; } + inline G4ThreeVector GetHitPos() const { return hit_pos; } inline void SetRot(G4RotationMatrix rmat) { rot = rmat; } inline G4RotationMatrix GetRot() const { return rot; } inline void SetLogV(G4LogicalVolume *val) { pLogV = val; } diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index bef91593..730393be 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -178,7 +178,8 @@ G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory if (_hitsCollection) { debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; - GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime); + G4ThreeVector pos = {worldPos.x(), worldPos.y(), worldPos.z()}; + GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime, pos); G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); aHit->SetLogV(thePhysical->GetLogicalVolume()); G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); diff --git a/src/geo/src/GeoFiberSensitiveDetectorHit.cc b/src/geo/src/GeoFiberSensitiveDetectorHit.cc index b78ca7fd..210b93d3 100644 --- a/src/geo/src/GeoFiberSensitiveDetectorHit.cc +++ b/src/geo/src/GeoFiberSensitiveDetectorHit.cc @@ -11,9 +11,10 @@ namespace RAT { G4Allocator GeoFiberSensitiveDetectorHitAllocator; -GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(G4int i, G4double t) { +GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(G4int i, G4double t, G4ThreeVector p) { id = i; time = t; + hit_pos = p; pLogV = 0; } GeoFiberSensitiveDetectorHit::~GeoFiberSensitiveDetectorHit() { ; } @@ -22,6 +23,7 @@ GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(const GeoFiberSensiti id = right.id; time = right.time; pos = right.pos; + hit_pos = right.hit_pos; rot = right.rot; pLogV = right.pLogV; } @@ -30,6 +32,7 @@ const GeoFiberSensitiveDetectorHit &GeoFiberSensitiveDetectorHit::operator=(cons id = right.id; time = right.time; pos = right.pos; + hit_pos = right.hit_pos; rot = right.rot; pLogV = right.pLogV; return *this; From 7aa8bd4a17389efdfbc3cbc7ddac3725277e1fc2 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Tue, 10 Dec 2024 10:42:24 +0000 Subject: [PATCH 17/21] Update header Forgot to update the header during last bugfix --- src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh index f51268af..0c9a344d 100644 --- a/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -13,7 +13,7 @@ namespace RAT { class GeoFiberSensitiveDetectorHit : public G4VHit { public: - GeoFiberSensitiveDetectorHit(G4int i, G4double t); + GeoFiberSensitiveDetectorHit(G4int i, G4double t, G4ThreeVector p); virtual ~GeoFiberSensitiveDetectorHit(); GeoFiberSensitiveDetectorHit(const GeoFiberSensitiveDetectorHit &right); const GeoFiberSensitiveDetectorHit &operator=(const GeoFiberSensitiveDetectorHit &right); From 02395dec8fc24fd591c9e56335b6802c158e8314 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Tue, 10 Dec 2024 11:24:37 +0000 Subject: [PATCH 18/21] Implement clang formatting --- src/core/src/Gsim.cc | 11 +- .../include/RAT/GeoFiberSensitiveDetector.hh | 3 +- src/geo/src/GeoFiberSensitiveDetector.cc | 269 +++++++++--------- src/geo/src/GeoFiberSensitiveDetectorHit.cc | 1 + 4 files changed, 140 insertions(+), 144 deletions(-) diff --git a/src/core/src/Gsim.cc b/src/core/src/Gsim.cc index fbee5284..fd79d352 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -21,9 +21,9 @@ #include #include #include -#include #include #include +#include #include #include #include @@ -569,13 +569,12 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) { } } mc->SetNumPE(numPE); - + /** hits from fibers */ G4HCofThisEvent *HC = g4ev->GetHCofThisEvent(); for (int hc = 0; hc < HC->GetNumberOfCollections(); hc++) { - GeoFiberSensitiveDetectorHitsCollection *hit_collection = (GeoFiberSensitiveDetectorHitsCollection*)HC->GetHC(hc); - if (hit_collection->GetName() != "FiberSenDet" || hit_collection->GetSize() == 0) - continue; + GeoFiberSensitiveDetectorHitsCollection *hit_collection = (GeoFiberSensitiveDetectorHitsCollection *)HC->GetHC(hc); + if (hit_collection->GetName() != "FiberSenDet" || hit_collection->GetSize() == 0) continue; DS::MCNestedTube *rat_mcnt = mc->AddNewMCNestedTube(); G4String det_name = hit_collection->GetSDname(); std::string fibre_id_str = det_name.erase(0, 6).data(); @@ -584,7 +583,7 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) { // only process fibers // info << hit_collection->GetSDname() << newline; for (int hit = 0; hit < hit_collection->GetSize(); hit++) { - GeoFiberSensitiveDetectorHit *my_hit = (GeoFiberSensitiveDetectorHit*)hit_collection->GetHit(hit); + GeoFiberSensitiveDetectorHit *my_hit = (GeoFiberSensitiveDetectorHit *)hit_collection->GetHit(hit); AddMCNestedTubeHit(rat_mcnt, my_hit); } } diff --git a/src/geo/include/RAT/GeoFiberSensitiveDetector.hh b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh index b4066d2d..dc45791c 100644 --- a/src/geo/include/RAT/GeoFiberSensitiveDetector.hh +++ b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh @@ -10,8 +10,7 @@ class G4TouchableHistory; namespace RAT { -class GeoFiberSensitiveDetector : public G4VSensitiveDetector -{ +class GeoFiberSensitiveDetector : public G4VSensitiveDetector { public: GeoFiberSensitiveDetector(G4String name); virtual ~GeoFiberSensitiveDetector(); diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc index 730393be..7a8717af 100644 --- a/src/geo/src/GeoFiberSensitiveDetector.cc +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -58,142 +58,139 @@ void GeoFiberSensitiveDetector::Initialize(G4HCofThisEvent *HCE) { G4bool GeoFiberSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory * /*ROhist*/) { // We ONLY want to store optical photons if (aStep->GetTrack()->GetDefinition()->GetParticleName() == "opticalphoton") { - - // only store information from steps where the photon is attenuated - G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); - G4StepPoint *postStepPoint = aStep->GetPostStepPoint(); - std::string proc = postStepPoint->GetProcessDefinedStep()->GetProcessName(); - if (proc != "Attenuation") - return true; - - debug << "GeoFiberSensitiveDetector::ProcessHits start." << newline; - debug << "GeoFiberSensitiveDetector::ProcessHits getting energy deposited." << newline; - - G4double edep = aStep->GetTotalEnergyDeposit(); - G4double dl = aStep->GetStepLength(); - - debug << " Energy deposited: " << edep << newline; - - // if(edep==0.) return true; - - G4Track *aTrack = aStep->GetTrack(); - G4ParticleDefinition *part = aTrack->GetDefinition(); - int pdg = part->GetPDGEncoding(); - - debug << " track information: " << newline << " G4Track Pointer: " << aTrack << newline - << " Particle Definition Pointer: " << part << newline << " Particle PDG Encoding: " << pdg - << newline; - - debug << "GeoFiberSensitiveDetector::ProcessHits getting global time." << newline; - - G4Material *m = preStepPoint->GetMaterial(); - G4String mname = m->GetName(); - G4double d = m->GetDensity(); - G4String f = m->GetChemicalFormula(); - - G4TouchableHistory *theTouchable = (G4TouchableHistory *)(preStepPoint->GetTouchable()); - // Get the copy number of this element. This may not be unique, - // if the mother is also a replica and occurs more than once. - G4int copyNo = theTouchable->GetVolume()->GetCopyNo(); - // Get the mother copy number, and offset it by 1000. This will help - // form the basis of the chamber ID. - G4int motherCopyNo = (theTouchable->GetVolume(1)->GetCopyNo()) + 1000; - G4int uid = motherCopyNo + copyNo; // unique identifier of chamber - - // construct a unique identifier for this copy - G4int ivol = 0; - G4int idOffset = 1; - uid = 0; - - debug << "History level: " << theTouchable->GetHistoryDepth() << newline; - - while (ivol < theTouchable->GetHistoryDepth()) { - debug << " * volume layer level = " << ivol << newline; - uid += theTouchable->GetVolume(ivol)->GetCopyNo() * idOffset; - idOffset *= 100; - ivol++; - } - - G4ThreeVector worldPos = preStepPoint->GetPosition(); - - debug << "Hit material name " << mname << newline; - debug << "density " << d << newline; - debug << "formula " << f << newline; - debug << "edep " << G4BestUnit(edep, "Energy") << newline; - debug << "dl " << G4BestUnit(dl, "Length") << newline; - debug << "pid " << pdg << newline; - debug << "Position " << G4BestUnit(worldPos.x(), "Length") << " " << G4BestUnit(worldPos.y(), "Length") << " " - << G4BestUnit(worldPos.z(), "Length") << " " << newline; - - debug << " " << newline; - - G4double hitTime = preStepPoint->GetGlobalTime(); - - debug << "GeoFiberSensitiveDetector::ProcessHits checking for an existing hit " - "in this element." - << newline; - // check if this finger already has a hit - G4int ix = -1; - - debug << "GeoFiberSensitiveDetector::ProcessHits hit collection address is " << _hitsCollection << newline; - // G4cerr << "GeoFiberSensitiveDetector: Hit ID = " << uid << " and position: " << worldPos << newline; - - int eventID = G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent(); - int trackID = aStep->GetTrack()->GetTrackID(); - - // get volume name of hit as well - // note that this volume name may not be unique - G4String vol = theTouchable->GetVolume()->GetName(); - - // if (fLastEventID != eventID || fLastTrackID != trackID) { - if (fLastTrackID != trackID) { - // Fill the hit information - _hit_x.push_back(worldPos.x()); - _hit_y.push_back(worldPos.y()); - _hit_z.push_back(worldPos.z()); - _hit_E.push_back(edep); - _hit_time.push_back(hitTime); - _hit_uid.push_back(uid); - _hit_pdg.push_back(pdg); - _hit_volume.push_back(vol); - fLastEventID = eventID; - fLastTrackID = trackID; - } - - /*G4int hitfreq = G4int(db["veto_hit_frequency"]);*/ - - if (NULL == _hitsCollection) { - debug << "GeoFiberSensitiveDetector::ProcessHits hit collection null. " - "Reloading from HCofEThisEvent." - << newline; - if (_HCE) { - _hitsCollection = (GeoFiberSensitiveDetectorHitsCollection *)(_HCE->GetHC(HCID)); - debug << "GeoFiberSensitiveDetector::ProcessHits * hit collection address is " << _hitsCollection << newline; - } else { - debug << "GeoFiberSensitiveDetector::ProcessHits (E) HCofEThisEvent " - "pointer is NULL!" - << newline; - } - } - - if (_hitsCollection) { - debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; - G4ThreeVector pos = {worldPos.x(), worldPos.y(), worldPos.z()}; - GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime, pos); - G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); - aHit->SetLogV(thePhysical->GetLogicalVolume()); - G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); - aTrans.Invert(); - aHit->SetRot(aTrans.NetRotation()); - aHit->SetPos(aTrans.NetTranslation()); - _hitsCollection->insert(aHit); - // aHit->Print(); - // aHit->Draw(); - debug << " * Drawing Hit " << uid << newline; - } - debug << "GeoFiberSensitiveDetector::ProcessHits end." << newline; - return true; - } + // only store information from steps where the photon is attenuated + G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); + G4StepPoint *postStepPoint = aStep->GetPostStepPoint(); + std::string proc = postStepPoint->GetProcessDefinedStep()->GetProcessName(); + if (proc != "Attenuation") return true; + + debug << "GeoFiberSensitiveDetector::ProcessHits start." << newline; + debug << "GeoFiberSensitiveDetector::ProcessHits getting energy deposited." << newline; + + G4double edep = aStep->GetTotalEnergyDeposit(); + G4double dl = aStep->GetStepLength(); + + debug << " Energy deposited: " << edep << newline; + + // if(edep==0.) return true; + + G4Track *aTrack = aStep->GetTrack(); + G4ParticleDefinition *part = aTrack->GetDefinition(); + int pdg = part->GetPDGEncoding(); + + debug << " track information: " << newline << " G4Track Pointer: " << aTrack << newline + << " Particle Definition Pointer: " << part << newline << " Particle PDG Encoding: " << pdg << newline; + + debug << "GeoFiberSensitiveDetector::ProcessHits getting global time." << newline; + + G4Material *m = preStepPoint->GetMaterial(); + G4String mname = m->GetName(); + G4double d = m->GetDensity(); + G4String f = m->GetChemicalFormula(); + + G4TouchableHistory *theTouchable = (G4TouchableHistory *)(preStepPoint->GetTouchable()); + // Get the copy number of this element. This may not be unique, + // if the mother is also a replica and occurs more than once. + G4int copyNo = theTouchable->GetVolume()->GetCopyNo(); + // Get the mother copy number, and offset it by 1000. This will help + // form the basis of the chamber ID. + G4int motherCopyNo = (theTouchable->GetVolume(1)->GetCopyNo()) + 1000; + G4int uid = motherCopyNo + copyNo; // unique identifier of chamber + + // construct a unique identifier for this copy + G4int ivol = 0; + G4int idOffset = 1; + uid = 0; + + debug << "History level: " << theTouchable->GetHistoryDepth() << newline; + + while (ivol < theTouchable->GetHistoryDepth()) { + debug << " * volume layer level = " << ivol << newline; + uid += theTouchable->GetVolume(ivol)->GetCopyNo() * idOffset; + idOffset *= 100; + ivol++; + } + + G4ThreeVector worldPos = preStepPoint->GetPosition(); + + debug << "Hit material name " << mname << newline; + debug << "density " << d << newline; + debug << "formula " << f << newline; + debug << "edep " << G4BestUnit(edep, "Energy") << newline; + debug << "dl " << G4BestUnit(dl, "Length") << newline; + debug << "pid " << pdg << newline; + debug << "Position " << G4BestUnit(worldPos.x(), "Length") << " " << G4BestUnit(worldPos.y(), "Length") << " " + << G4BestUnit(worldPos.z(), "Length") << " " << newline; + + debug << " " << newline; + + G4double hitTime = preStepPoint->GetGlobalTime(); + + debug << "GeoFiberSensitiveDetector::ProcessHits checking for an existing hit " + "in this element." + << newline; + // check if this finger already has a hit + G4int ix = -1; + + debug << "GeoFiberSensitiveDetector::ProcessHits hit collection address is " << _hitsCollection << newline; + // G4cerr << "GeoFiberSensitiveDetector: Hit ID = " << uid << " and position: " << worldPos << newline; + + int eventID = G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent(); + int trackID = aStep->GetTrack()->GetTrackID(); + + // get volume name of hit as well + // note that this volume name may not be unique + G4String vol = theTouchable->GetVolume()->GetName(); + + // if (fLastEventID != eventID || fLastTrackID != trackID) { + if (fLastTrackID != trackID) { + // Fill the hit information + _hit_x.push_back(worldPos.x()); + _hit_y.push_back(worldPos.y()); + _hit_z.push_back(worldPos.z()); + _hit_E.push_back(edep); + _hit_time.push_back(hitTime); + _hit_uid.push_back(uid); + _hit_pdg.push_back(pdg); + _hit_volume.push_back(vol); + fLastEventID = eventID; + fLastTrackID = trackID; + } + + /*G4int hitfreq = G4int(db["veto_hit_frequency"]);*/ + + if (NULL == _hitsCollection) { + debug << "GeoFiberSensitiveDetector::ProcessHits hit collection null. " + "Reloading from HCofEThisEvent." + << newline; + if (_HCE) { + _hitsCollection = (GeoFiberSensitiveDetectorHitsCollection *)(_HCE->GetHC(HCID)); + debug << "GeoFiberSensitiveDetector::ProcessHits * hit collection address is " << _hitsCollection << newline; + } else { + debug << "GeoFiberSensitiveDetector::ProcessHits (E) HCofEThisEvent " + "pointer is NULL!" + << newline; + } + } + + if (_hitsCollection) { + debug << "GeoFiberSensitiveDetector::ProcessHits creating a new hit." << newline; + G4ThreeVector pos = {worldPos.x(), worldPos.y(), worldPos.z()}; + GeoFiberSensitiveDetectorHit *aHit = new GeoFiberSensitiveDetectorHit(uid, hitTime, pos); + G4VPhysicalVolume *thePhysical = theTouchable->GetVolume(); + aHit->SetLogV(thePhysical->GetLogicalVolume()); + G4AffineTransform aTrans = theTouchable->GetHistory()->GetTopTransform(); + aTrans.Invert(); + aHit->SetRot(aTrans.NetRotation()); + aHit->SetPos(aTrans.NetTranslation()); + _hitsCollection->insert(aHit); + // aHit->Print(); + // aHit->Draw(); + debug << " * Drawing Hit " << uid << newline; + } + debug << "GeoFiberSensitiveDetector::ProcessHits end." << newline; + return true; + } return true; } diff --git a/src/geo/src/GeoFiberSensitiveDetectorHit.cc b/src/geo/src/GeoFiberSensitiveDetectorHit.cc index 210b93d3..a150a0fa 100644 --- a/src/geo/src/GeoFiberSensitiveDetectorHit.cc +++ b/src/geo/src/GeoFiberSensitiveDetectorHit.cc @@ -1,4 +1,5 @@ #include + #include #include #include From 71dc358f27867b6f573cda32ae54e932ddef1894 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 12 Dec 2024 16:53:41 +0000 Subject: [PATCH 19/21] Add example WLS materials and use in validation --- macros/vis_nestedTubeArrays.mac | 10 +- ratdb/Validation/NestedTubeArray.geo | 24 ++--- ratdb/WLSExample.ratdb | 146 +++++++++++++++++++++++++++ 3 files changed, 164 insertions(+), 16 deletions(-) create mode 100644 ratdb/WLSExample.ratdb diff --git a/macros/vis_nestedTubeArrays.mac b/macros/vis_nestedTubeArrays.mac index 6e54e047..045fade7 100644 --- a/macros/vis_nestedTubeArrays.mac +++ b/macros/vis_nestedTubeArrays.mac @@ -29,7 +29,9 @@ /vis/scene/add/hits /vis/sceneHandler/attach /vis/viewer/set/upVector 0.0 0.0 1.0 -/vis/viewer/set/viewpointThetaPhi -90 135 +/vis/viewer/set/viewpointThetaPhi -90 135 +/vis/modeling/trajectories/create/drawByCharge +/vis/modeling/trajectories/drawByCharge-0/setRGBA 0 0.236 0.240 0.002 1 /vis/viewer/set/style s /vis/viewer/flush @@ -39,15 +41,15 @@ #add ntuple tracking /rat/proc outntuple -/rat/procset include_tracking 1 +/rat/procset include_tracking 0 /rat/procset include_mcparticles 1 /rat/procset include_pmthits 1 /rat/procset include_untriggered_events 1 /rat/physics/setOpWLS g4 ##### GENERATORS ################# /generator/add combo pbomb:point:poisson -/generator/vtx/set 10 430 # 10000 photons, 430nm +/generator/vtx/set 1000 350 # 1000 photons, 350nm /generator/pos/set 0.0 0.0 0.0 ##### RUN ########### -/run/beamOn 1 \ No newline at end of file +/run/beamOn 1 diff --git a/ratdb/Validation/NestedTubeArray.geo b/ratdb/Validation/NestedTubeArray.geo index 533ea83e..0f6a090f 100644 --- a/ratdb/Validation/NestedTubeArray.geo +++ b/ratdb/Validation/NestedTubeArray.geo @@ -6,7 +6,7 @@ mother: "", type: "box", size: [20000.0,20000.0,20000.0], - material: "mirror", + material: "air", } { @@ -33,7 +33,7 @@ size: [9.0,9.0,9.0], position: [0.0, 0.0, 0.0], rotation: [0.0, 0.0, 0.0], - material: "mirror", + material: "air", #color: [0.02,0.2,0.2,0.1], } @@ -49,7 +49,7 @@ size_z: 49.5, position: [-5.0, 0.0, -5.0], rotation: [-90.0, 0.0, 0.0], - material: "aluminum", + material: "Fpolyethylene", color: [0.0,0.8,0.0,0.01], } @@ -64,7 +64,7 @@ size_z: 49.5, position: [0.0, 0.0, 0.0], rotation: [0.0, 0.0, 0.0], - material: "glass", + material: "PMMA", color: [0.0,0.8,0.0,0.05], } @@ -79,7 +79,7 @@ size_z: 49.5, position: [0.0, 0.0, 0.0], rotation: [0.0, 0.0, 0.0], - material: "mirror", + material: "WLSExample", color: [0.0,0.8,0.0,0.1], } { @@ -93,7 +93,7 @@ size_z: 29.5, position: [-5.0, 0.0, 5.0], rotation: [-90.0, 0.0, 0.0], - material: "aluminum", + material: "Fpolyethylene", color: [0.0,0.8,0.0,0.01], } @@ -108,7 +108,7 @@ size_z: 29.5, position: [0.0, 0.0, 0.0], rotation: [0.0, 0.0, 0.0], - material: "glass", + material: "PMMA", color: [0.0,0.8,0.0,0.05], } @@ -123,7 +123,7 @@ size_z: 29.5, position: [0.0, 0.0, 0.0], rotation: [0.0, 0.0, 0.0], - material: "mirror", + material: "WLSExample", color: [0.0,0.8,0.0,0.1], } // manual nested tubes ends here @@ -162,9 +162,9 @@ inner_r: 0.485, outer_r: 0.5, pos_table: "cable_pos", orientation: "manual", -material_outer: "aluminum", -material_inner: "glass", -material_core: "mirror", +outer_material: "Fpolyethylene", +inner_material: "PMMA", +core_material: "WLSExample", #drawstyle: "solid", -color: [0.8,0.0,0.0,0.8] +color: [0.0,0.8,0.0,0.2] } diff --git a/ratdb/WLSExample.ratdb b/ratdb/WLSExample.ratdb new file mode 100644 index 00000000..b4d089e4 --- /dev/null +++ b/ratdb/WLSExample.ratdb @@ -0,0 +1,146 @@ +//------------------------------------------------// +//Wavelength Shifting (WLS) Material Example +//------------------------------------------------// +// From https://en.wikipedia.org/wiki/Polystyrene +// Composition (C8H8)n +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf +// Density 1.05 g / cm3 +{ +name: "MATERIAL", +index: "WLSExample", +valid_begin : [0, 0], +valid_end : [0, 0], +density: 1.050, +nelements: 2, +nmaterials: 0, +elements: ["Hydrogen", "Carbon"], +elemprop: [0.0774, 0.9226], +index_of_refraction: 1.59 +} +{ +name: "OPTICS", +index: "WLSExample", +valid_begin : [0, 0], +valid_end : [0, 0], +surface: 1, +finish: "polished", +model: "unified", +polish: 1.0, +LIGHT_YIELD: 0.0, // Does not scintillate! +//------------------------------------------------// +// Refractive Index +//------------------------------------------------// +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf, page 3 +RINDEX_option: "wavelength", +RINDEX_value1: [ 250.0, 800.0 ], +RINDEX_value2: [ 1.590, 1.590 ], +//------------------------------------------------// +// Re-emission spectrum +//------------------------------------------------// +SCINTILLATION_option: "dy_dwavelength", +SCINTILLATION_value1: [ 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 800.0 ], +SCINTILLATION_value2: [ 0.000, 0.000, 0.000, 0.000, 1.000, 1.000, 0.000, 0.000, 0.000 ], +SCINTILLATION_WLS_option: "dy_dwavelength", +SCINTILLATION_WLS_value1: [ 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 800.0 ], +SCINTILLATION_WLS_value2: [ 0.000, 0.000, 0.000, 0.000, 1.000, 1.000, 0.000, 0.000, 0.000 ], +//------------------------------------------------// +// Quantum efficiency of the WLS dye a function of absorbed wavelength +//------------------------------------------------// +REEMISSION_PROB_option: "wavelength", +REEMISSION_PROB_value1: [ 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 800.0 ], +REEMISSION_PROB_value2: [ 0.900, 0.900, 0.900, 0.900, 0.900, 0.900, 0.900, 0.900, 0.900 ], +//------------------------------------------------// +// Time constants for WLS re-emission +//------------------------------------------------// +// Decay time assumed to be 10.0 ns +// Rise time is assumed to be 0.1 ns +REEMITWAVEFORM_value1: [-10.0], +REEMITWAVEFORM_value2: [1.0], +SCINTWAVEFORM_value1: [-10.0], +SCINTWAVEFORM_value2: [1.0], +SCINT_RISE_TIME: 0.1, +//------------------------------------------------// +// ABSORPTION LENGTH +//------------------------------------------------// +ABSLENGTH_option: "wavelength", +ABSLENGTH_value1: [ 250.0, 300.0, 350.0, 400.0, 450.00, 500.00, 550.00, 600.00, 800.00 ], +ABSLENGTH_value2: [ 0.000, 0.500, 0.500, 0.500, 4000.0, 4000.0, 4000.0, 4000.0, 4000.0 ], +PROPERTY_LIST: ["LIGHT_YIELD", "RINDEX", "SCINTILLATION", "SCINTILLATION_WLS", "ABSLENGTH", "SCINT_RISE_TIME", "SCINTWAVEFORM", "REEMISSION_PROB", "REEMITWAVEFORM"] +} +//////////////////////////////////////////////////// +//------------------------------------------------// +//Fiber Cladding Materials: +//------------------------------------------------// +//////////////////////////////////////////////////// +//------------------------------------------------// +//Inner cladding of fiber: PMMA +//------------------------------------------------// +// From https://en.wikipedia.org/wiki/Poly(methyl_methacrylate) +// Composition (C5O2H8)n +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf, page 3 +// Density 1.19 g / cm3 +{ +name: "MATERIAL", +index: "PMMA", +valid_begin : [0, 0], +valid_end : [0, 0], +density: 1.19, +nelements: 3, +nmaterials: 0, +elements: ["Hydrogen", "Carbon", "Oxygen"], +elemprop: [0.080538, 0.599848, 0.319614], +index_of_refraction: 1.49 +} +{ +name: "OPTICS", +index: "PMMA", +valid_begin : [0, 0], +valid_end : [0, 0], +//------------------------------------------------// +// REFRACTIVE INDEX +//------------------------------------------------// +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf, page 3 +RINDEX_option: "wavelength", +RINDEX_value1: [250.0, 800.0], +RINDEX_value2: [1.49, 1.49], +PROPERTY_LIST: ["RINDEX"], +} +//------------------------------------------------// +//Outer cladding of fiber: Fluorinated polyethylene +//------------------------------------------------// +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf, page 3 +// Density 1.43 g / cm3 +{ +name: "MATERIAL", +index: "Fpolyethylene", +valid_begin : [0, 0], +valid_end : [0, 0], +density: 1.430, +nelements: 3, +nmaterials: 0, +elements: ["Hydrogen", "Carbon", "Fluorine"], +// H, C and F atoms.. let's do 8:3:3 (random guess) +elemprop: [0.08, 0.36, 0.56], +index_of_refraction: 1.42 +} +{ +name: "OPTICS", +index: "Fpolyethylene", +valid_begin : [0, 0], +valid_end : [0, 0], +//------------------------------------------------// +// REFRACTIVE INDEX +//------------------------------------------------// +// From +// https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf, page 3 +RINDEX_option: "wavelength", +RINDEX_value1: [250.0, 800.0], +RINDEX_value2: [1.42, 1.42], +PROPERTY_LIST: ["RINDEX"], +} + From a06b9d6d0fdd3cae61acb0f470eac59ef9047930 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Fri, 13 Dec 2024 17:39:23 +0000 Subject: [PATCH 20/21] Add nested tube hit info to ntuple output --- macros/vis_nestedTubeArrays.mac | 2 ++ ratdb/Validation/NestedTubeArray.geo | 3 +- src/ds/include/RAT/DS/LinkDef.hh | 4 +++ src/io/include/RAT/OutNtupleProc.hh | 18 +++++++++- src/io/src/OutNtupleProc.cc | 53 ++++++++++++++++++++++++++++ 5 files changed, 78 insertions(+), 2 deletions(-) diff --git a/macros/vis_nestedTubeArrays.mac b/macros/vis_nestedTubeArrays.mac index 045fade7..7e714aeb 100644 --- a/macros/vis_nestedTubeArrays.mac +++ b/macros/vis_nestedTubeArrays.mac @@ -44,7 +44,9 @@ /rat/procset include_tracking 0 /rat/procset include_mcparticles 1 /rat/procset include_pmthits 1 +/rat/procset include_nestedtubehits 1 /rat/procset include_untriggered_events 1 +#/rat/proc outroot /rat/physics/setOpWLS g4 ##### GENERATORS ################# /generator/add combo pbomb:point:poisson diff --git a/ratdb/Validation/NestedTubeArray.geo b/ratdb/Validation/NestedTubeArray.geo index 0f6a090f..6e21cbbb 100644 --- a/ratdb/Validation/NestedTubeArray.geo +++ b/ratdb/Validation/NestedTubeArray.geo @@ -166,5 +166,6 @@ outer_material: "Fpolyethylene", inner_material: "PMMA", core_material: "WLSExample", #drawstyle: "solid", -color: [0.0,0.8,0.0,0.2] +color: [0.0,0.8,0.0,0.2], +sensitive_detector: "/mydet/fibers" } diff --git a/src/ds/include/RAT/DS/LinkDef.hh b/src/ds/include/RAT/DS/LinkDef.hh index 17821562..054cf376 100644 --- a/src/ds/include/RAT/DS/LinkDef.hh +++ b/src/ds/include/RAT/DS/LinkDef.hh @@ -68,6 +68,8 @@ #pragma link C++ class vector < RAT::DS::MCTrackStep>; #pragma link C++ class vector < RAT::DS::MCPMT>; #pragma link C++ class vector < RAT::DS::MCPhoton>; +#pragma link C++ class vector < RAT::DS::MCNestedTube>; +#pragma link C++ class vector < RAT::DS::MCNestedTubeHit>; #pragma link C++ class vector < RAT::DS::Calib>; #pragma link C++ class vector < RAT::DS::EV>; #pragma link C++ class vector < RAT::DS::PMT>; @@ -89,6 +91,8 @@ #pragma link C++ class vector < RAT::DS::MCTrack *>; #pragma link C++ class vector < RAT::DS::MCPMT *>; #pragma link C++ class vector < RAT::DS::MCPhoton *>; +#pragma link C++ class vector < RAT::DS::MCNestedTube *>; +#pragma link C++ class vector < RAT::DS::MCNestedTubeHit *>; #pragma link C++ class vector < RAT::DS::MCTrackStep *>; #endif diff --git a/src/io/include/RAT/OutNtupleProc.hh b/src/io/include/RAT/OutNtupleProc.hh index 33499982..d6a80495 100644 --- a/src/io/include/RAT/OutNtupleProc.hh +++ b/src/io/include/RAT/OutNtupleProc.hh @@ -51,6 +51,7 @@ class OutNtupleProc : public Processor { bool pmthits; bool untriggered; bool mchits; + bool nthits; }; NtupleOptions options; @@ -74,6 +75,13 @@ class OutNtupleProc : public Processor { std::vector pmtU; std::vector pmtV; std::vector pmtW; + std::vector ntId; + std::vector ntX; + std::vector ntY; + std::vector ntZ; + std::vector ntU; + std::vector ntV; + std::vector ntW; // Data Branches Int_t mcpdg; double mcx, mcy, mcz; @@ -95,6 +103,14 @@ class OutNtupleProc : public Processor { int mcpecount; std::vector mcpmtid; std::vector mcpmtnpe; + // MCNestedTube + int mcnNTs; + int mcnNThits; + std::vector mcNTid; + std::vector mcNThittime; + std::vector mcNThitx; + std::vector mcNThity; + std::vector mcNThitz; // MCPE std::vector mcpehittime; std::vector mcpefrontendtime; @@ -125,11 +141,11 @@ class OutNtupleProc : public Processor { std::vector hitPMTDigitizedTime; std::vector hitPMTDigitizedCharge; std::vector hitPMTNCrossings; - // Tracking std::map processCodeMap; std::vector processCodeIndex; std::vector processName; + // Tracking std::vector trackPDG; std::vector> trackPosX; std::vector> trackPosY; diff --git a/src/io/src/OutNtupleProc.cc b/src/io/src/OutNtupleProc.cc index 0710ee37..ab89de04 100644 --- a/src/io/src/OutNtupleProc.cc +++ b/src/io/src/OutNtupleProc.cc @@ -48,12 +48,14 @@ OutNtupleProc::OutNtupleProc() : Processor("outntuple") { options.pmthits = table->GetZ("include_pmthits"); options.untriggered = table->GetZ("include_untriggered_events"); options.mchits = table->GetZ("include_mchits"); + options.nthits = table->GetZ("include_nestedtubehits"); } catch (DBNotFoundError &e) { options.tracking = false; options.mcparticles = false; options.pmthits = true; options.untriggered = false; options.mchits = true; + options.nthits = false; } } @@ -129,6 +131,22 @@ bool OutNtupleProc::OpenFile(std::string filename) { outputTree->Branch("hitPMTDigitizedCharge", &hitPMTDigitizedCharge); outputTree->Branch("hitPMTNCrossings", &hitPMTNCrossings); } + if (options.nthits) { + outputTree->Branch("mcnNTs", &mcnNTs); + outputTree->Branch("mcnNThits", &mcnNThits); + outputTree->Branch("mcNTid", &mcNTid); + outputTree->Branch("mcNThittime", &mcNThittime); + outputTree->Branch("mcNThitx", &mcNThitx); + outputTree->Branch("mcNThity", &mcNThity); + outputTree->Branch("mcNThitz", &mcNThitz); + metaTree->Branch("ntId", &ntId); + metaTree->Branch("ntX", &ntX); + metaTree->Branch("ntY", &ntY); + metaTree->Branch("ntZ", &ntZ); + metaTree->Branch("ntU", &ntU); + metaTree->Branch("ntV", &ntV); + metaTree->Branch("ntW", &ntW); + } if (options.mchits) { // Save full MC PMT hit information outputTree->Branch("mcPMTID", &mcpmtid); @@ -171,6 +189,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) { } runBranch = DS::RunStore::GetRun(ds); DS::PMTInfo *pmtinfo = runBranch->GetPMTInfo(); + DS::NestedTubeInfo *ntinfo = runBranch->GetNestedTubeInfo(); ULong64_t stonano = 1000000000; dsentries++; // Clear the previous vectors @@ -329,6 +348,23 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) { } } } + if (options.nthits) { + mcnNTs = mc->GetMCNestedTubeCount(); + mcnNThits = 0; + for (int iNT = 0; iNT < mc->GetMCNestedTubeCount(); iNT++) { + DS::MCNestedTube *mcnt = mc->GetMCNestedTube(iNT); + mcnNThits += mcnt->GetMCNestedTubeHitCount(); + mcNTid.push_back(mcnt->GetID()); + G4ThreeVector position = ntinfo->GetPosition(mcnt->GetID()); + for (int ih = 0; ih < mcnt->GetMCNestedTubeHitCount(); ih++) { + RAT::DS::MCNestedTubeHit *mcph = mcnt->GetMCNestedTubeHit(ih); + mcNThittime.push_back(mcph->GetHitTime()); + mcNThitx.push_back(position.x()); + mcNThity.push_back(position.y()); + mcNThitz.push_back(position.z()); + } + } + } // EV Branches for (subev = 0; subev < ds->GetEVCount(); subev++) { @@ -463,6 +499,20 @@ OutNtupleProc::~OutNtupleProc() { pmtV.push_back(direction.Y()); pmtW.push_back(direction.Z()); } + if (options.nthits) { + DS::NestedTubeInfo *ntinfo = runBranch->GetNestedTubeInfo(); + for (int id = 0; id < ntinfo->GetNestedTubeCount(); id++) { + G4ThreeVector position = ntinfo->GetPosition(id); + G4ThreeVector direction = ntinfo->GetDirection(id); + ntId.push_back(id); + ntX.push_back(position.x()); + ntY.push_back(position.y()); + ntZ.push_back(position.z()); + ntU.push_back(direction.x()); + ntV.push_back(direction.y()); + ntW.push_back(direction.z()); + } + } runId = runBranch->GetID(); runType = runBranch->GetType(); // Converting to unix time @@ -526,6 +576,9 @@ void OutNtupleProc::SetI(std::string param, int value) { if (param == "include_pmthits") { options.pmthits = value ? true : false; } + if (param == "include_nestedtubehits") { + options.nthits = value ? true : false; + } if (param == "include_untriggered_events") { options.untriggered = value ? true : false; } From 4c60a0615685d2236e279cd63cf4d0e494bdfc35 Mon Sep 17 00:00:00 2001 From: Wilf Shorrock Date: Thu, 9 Jan 2025 15:54:19 +0000 Subject: [PATCH 21/21] Fix formatting --- src/core/include/RAT/Gsim.hh | 4 +- src/ds/include/RAT/DS/MC.hh | 2 +- src/ds/include/RAT/DS/MCNestedTubeHit.hh | 2 +- src/ds/include/RAT/DS/MCPhoton.hh | 2 +- src/ds/include/RAT/DS/Run.hh | 2 +- src/ds/src/Root.cc | 4 +- .../include/RAT/GeoNestedTubeConstruction.hh | 69 +++++++++---------- src/geo/src/GeoNestedTubeConstruction.cc | 11 ++- 8 files changed, 47 insertions(+), 49 deletions(-) diff --git a/src/core/include/RAT/Gsim.hh b/src/core/include/RAT/Gsim.hh index 3f9ee826..2b681d66 100644 --- a/src/core/include/RAT/Gsim.hh +++ b/src/core/include/RAT/Gsim.hh @@ -9,14 +9,14 @@ #include #include #include -#include #include +#include #include -#include #include #include #include #include +#include #include #include diff --git a/src/ds/include/RAT/DS/MC.hh b/src/ds/include/RAT/DS/MC.hh index b40b69ac..c4030b54 100644 --- a/src/ds/include/RAT/DS/MC.hh +++ b/src/ds/include/RAT/DS/MC.hh @@ -16,8 +16,8 @@ #include #include -#include #include +#include #include #include #include diff --git a/src/ds/include/RAT/DS/MCNestedTubeHit.hh b/src/ds/include/RAT/DS/MCNestedTubeHit.hh index 039a774d..e48f196d 100644 --- a/src/ds/include/RAT/DS/MCNestedTubeHit.hh +++ b/src/ds/include/RAT/DS/MCNestedTubeHit.hh @@ -4,7 +4,7 @@ * * @author Wilf Shorrock * - * This class represents a single photon that is captured within the core of a + * This class represents a single photon that is captured within the core of a * nested tube object. */ diff --git a/src/ds/include/RAT/DS/MCPhoton.hh b/src/ds/include/RAT/DS/MCPhoton.hh index 5ea725d1..7d8e83b6 100644 --- a/src/ds/include/RAT/DS/MCPhoton.hh +++ b/src/ds/include/RAT/DS/MCPhoton.hh @@ -8,7 +8,7 @@ * photocathode of the PMT. The time jitter and delay in transit to the * anode are not included here, but the distribution of charge is, * which is slightly incongruous. - + * Note that we require that the photon generates a photoelectron. * Absorbed photons are not included here. */ diff --git a/src/ds/include/RAT/DS/Run.hh b/src/ds/include/RAT/DS/Run.hh index 54bae730..b8d09d45 100644 --- a/src/ds/include/RAT/DS/Run.hh +++ b/src/ds/include/RAT/DS/Run.hh @@ -11,8 +11,8 @@ #include #include -#include #include +#include #include #include diff --git a/src/ds/src/Root.cc b/src/ds/src/Root.cc index 66ad8b44..2a182461 100644 --- a/src/ds/src/Root.cc +++ b/src/ds/src/Root.cc @@ -3,11 +3,11 @@ #include #include #include -#include #include +#include +#include #include #include -#include #include #include #include diff --git a/src/geo/include/RAT/GeoNestedTubeConstruction.hh b/src/geo/include/RAT/GeoNestedTubeConstruction.hh index 50112f85..4544b5a7 100644 --- a/src/geo/include/RAT/GeoNestedTubeConstruction.hh +++ b/src/geo/include/RAT/GeoNestedTubeConstruction.hh @@ -7,56 +7,55 @@ #include #include #include -#include #include #include #include +#include #include #include namespace RAT { - struct GeoNestedTubeConstructionParams { - GeoNestedTubeConstructionParams() { invisible = false; }; +struct GeoNestedTubeConstructionParams { + GeoNestedTubeConstructionParams() { invisible = false; }; + + bool invisible; + + double outer_r; + double inner_r; + double core_r; + double Dz; // half length - bool invisible; + G4Material *outer; + G4Material *inner; + G4Material *core; - double outer_r; - double inner_r; - double core_r; - double Dz; // half length + // G4OpticalSurface *outer_inner; + G4OpticalSurface *inner_core; +}; - G4Material *outer; - G4Material *inner; - G4Material *core; +class GeoNestedTubeConstruction { + public: + GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID); + virtual ~GeoNestedTubeConstruction() {} - // G4OpticalSurface *outer_inner; - G4OpticalSurface *inner_core; - }; + virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID, DBLinkPtr table); + virtual G4VSolid *BuildSolid(const std::string &prefix); + virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, + G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, + int copyNo); - class GeoNestedTubeConstruction { - public: - GeoNestedTubeConstruction(DBLinkPtr table, DBLinkPtr postable, G4LogicalVolume *mother, int ID); - virtual ~GeoNestedTubeConstruction() {} + protected: + // physical volumes + G4PVPlacement *inner_phys; + G4PVPlacement *core_phys; - virtual G4LogicalVolume *BuildVolume(const std::string &prefix, int ID, DBLinkPtr table); - virtual G4VSolid *BuildSolid(const std::string &prefix); - virtual G4PVPlacement *PlaceNestedTube(G4RotationMatrix *tuberot, G4ThreeVector tubepos, const std::string &name, - G4LogicalVolume *logi_tube, G4VPhysicalVolume *mother_phys, bool booleanSolid, - int copyNo); + G4LogicalVolume *log_tube; + GeoNestedTubeConstructionParams fParams; - protected: - // physical volumes - G4PVPlacement *inner_phys; - G4PVPlacement *core_phys; + DBLinkPtr myTable; +}; - G4LogicalVolume *log_tube; - GeoNestedTubeConstructionParams fParams; +} // namespace RAT - DBLinkPtr myTable; - }; - -} // namespace RAT - #endif - \ No newline at end of file diff --git a/src/geo/src/GeoNestedTubeConstruction.cc b/src/geo/src/GeoNestedTubeConstruction.cc index 31bab4eb..a7df66a1 100644 --- a/src/geo/src/GeoNestedTubeConstruction.cc +++ b/src/geo/src/GeoNestedTubeConstruction.cc @@ -1,13 +1,13 @@ #include #include #include -#include #include +#include #include #include -#include -#include #include +#include +#include #include #include #include @@ -76,12 +76,11 @@ G4LogicalVolume *GeoNestedTubeConstruction::BuildVolume(const std::string &prefi // set as sensitive detector if applicable try { std::string sensitive_detector = table->GetS("sensitive_detector"); - GeoFiberSensitiveDetector* sensitive = new GeoFiberSensitiveDetector(sensitive_detector+std::to_string(ID)); + GeoFiberSensitiveDetector *sensitive = new GeoFiberSensitiveDetector(sensitive_detector + std::to_string(ID)); G4SDManager *sdman = G4SDManager::GetSDMpointer(); sdman->AddNewDetector(sensitive); core_log->SetSensitiveDetector(sensitive); - } - catch (DBNotFoundError &e) { + } catch (DBNotFoundError &e) { } // ------------ Physical Volumes -------------