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..ef9ba731 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: @@ -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,29 @@ 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) +``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 +``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. +``sensitive_detector`` ``string`` (optional) Name of sensitive detector if this volume should register hits. Limited to ''/mydet/fibers''.' +====================== ========================== =================== + Creating a parameterized geometry ````````````````````````````````` Using a ``DetectorFactory`` one can build a DB defined geometry on the fly diff --git a/macros/vis_nestedTubeArrays.mac b/macros/vis_nestedTubeArrays.mac new file mode 100644 index 00000000..7e714aeb --- /dev/null +++ b/macros/vis_nestedTubeArrays.mac @@ -0,0 +1,57 @@ +/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/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 + +/rat/proc count +/rat/procset update 10 + + +#add ntuple tracking +/rat/proc outntuple +/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 +/generator/vtx/set 1000 350 # 1000 photons, 350nm +/generator/pos/set 0.0 0.0 0.0 + +##### RUN ########### +/run/beamOn 1 diff --git a/ratdb/Validation/NestedTubeArray.geo b/ratdb/Validation/NestedTubeArray.geo new file mode 100644 index 00000000..6e21cbbb --- /dev/null +++ b/ratdb/Validation/NestedTubeArray.geo @@ -0,0 +1,171 @@ +{ + name: "GEO", + index: "world", + valid_begin: [0, 0], + valid_end: [0, 0], + mother: "", + type: "box", + size: [20000.0,20000.0,20000.0], + material: "air", +} + +{ + 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: "air", + #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: "Fpolyethylene", + 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: "PMMA", + 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: "WLSExample", + 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: "Fpolyethylene", + 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: "PMMA", + 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: "WLSExample", + 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", +outer_material: "Fpolyethylene", +inner_material: "PMMA", +core_material: "WLSExample", +#drawstyle: "solid", +color: [0.0,0.8,0.0,0.2], +sensitive_detector: "/mydet/fibers" +} 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"], +} + diff --git a/src/core/include/RAT/Gsim.hh b/src/core/include/RAT/Gsim.hh index e2085528..2b681d66 100644 --- a/src/core/include/RAT/Gsim.hh +++ b/src/core/include/RAT/Gsim.hh @@ -10,11 +10,13 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include @@ -74,6 +76,8 @@ 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); void SetOpticalPhotonIDs(std::string particle_type, int trackID, int parentID); @@ -85,6 +89,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 409c87a2..de1d7d1d 100644 --- a/src/core/src/Gsim.cc +++ b/src/core/src/Gsim.cc @@ -21,6 +21,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -179,6 +182,7 @@ void Gsim::BeginOfRunAction(const G4Run * /*aRun*/) { run = DS::RunStore::GetRun(runID); fPMTInfo = run->GetPMTInfo(); + fNestedTubeInfo = run->GetNestedTubeInfo(); GLG4VEventAction::GetTheHitPMTCollection()->SetChannelStatus(run->GetChannelStatus()); for (size_t i = 0; i < fPMTTime.size(); i++) { @@ -421,6 +425,7 @@ void Gsim::MakeRun(int _runID) { run->SetID(_runID); run->SetType((unsigned)lrun->GetI("runtype")); run->SetStartTime(utc); + run->SetNestedTubeInfo(&GeoNestedSolidArrayFactoryBase::GetNestedTubeInfo()); const DS::PMTInfo *pmtinfo = &PMTFactoryBase::GetPMTInfo(); run->SetPMTInfo(pmtinfo); DS::ChannelStatus ch_status; @@ -567,6 +572,24 @@ 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; + 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) { @@ -596,6 +619,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->GetHitPos(); + 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 946e7009..dc3ea64a 100644 --- a/src/ds/CMakeLists.txt +++ b/src/ds/CMakeLists.txt @@ -31,11 +31,14 @@ 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 RAT/DS/Run.hh RAT/DS/PMTInfo.hh + RAT/DS/NestedTubeInfo.hh RAT/DS/ChannelStatus.hh RAT/DS/MCTrack.hh RAT/DS/MCTrackStep.hh diff --git a/src/ds/include/RAT/DS/LinkDef.hh b/src/ds/include/RAT/DS/LinkDef.hh index eebcaa76..f01fec7d 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::ChannelStatus + ; #pragma link C++ class RAT::DS::MC + ; @@ -9,7 +10,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 + ; @@ -68,6 +71,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>; @@ -93,6 +98,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/ds/include/RAT/DS/MC.hh b/src/ds/include/RAT/DS/MC.hh index a894b52b..c4030b54 100644 --- a/src/ds/include/RAT/DS/MC.hh +++ b/src/ds/include/RAT/DS/MC.hh @@ -16,6 +16,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..e48f196d --- /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..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/NestedTubeInfo.hh b/src/ds/include/RAT/DS/NestedTubeInfo.hh new file mode 100644 index 00000000..e2ee6b61 --- /dev/null +++ b/src/ds/include/RAT/DS/NestedTubeInfo.hh @@ -0,0 +1,89 @@ +/** + * @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() {} + + 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 170d6430..b8d09d45 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 #include @@ -51,6 +52,24 @@ 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); + /** Channel status */ virtual ChannelStatus const *GetChannelStatus() const { return &ch_status; } virtual void SetChannelStatus(const ChannelStatus &_ch_status) { ch_status = _ch_status; } @@ -61,6 +80,7 @@ class Run : public TObject { Int_t id; ULong64_t type; TTimeStamp startTime; + std::vector nestedtubeinfo; std::vector pmtinfo; // ah.. why is this a vector? ChannelStatus ch_status; }; diff --git a/src/ds/src/Root.cc b/src/ds/src/Root.cc index 3cb9708b..2a182461 100644 --- a/src/ds/src/Root.cc +++ b/src/ds/src/Root.cc @@ -3,6 +3,8 @@ #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/CMakeLists.txt b/src/geo/CMakeLists.txt index fb4da37d..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 @@ -33,6 +35,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/GeoFiberSensitiveDetector.hh b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh new file mode 100644 index 00000000..dc45791c --- /dev/null +++ b/src/geo/include/RAT/GeoFiberSensitiveDetector.hh @@ -0,0 +1,53 @@ +#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..0c9a344d --- /dev/null +++ b/src/geo/include/RAT/GeoFiberSensitiveDetectorHit.hh @@ -0,0 +1,67 @@ +#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, G4ThreeVector p); + 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; + G4ThreeVector hit_pos; + G4RotationMatrix rot; + const G4LogicalVolume *pLogV; + std::string proc; + + 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 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; } + 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/GeoNestedSolidArrayFactoryBase.hh b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh new file mode 100644 index 00000000..c7cec1fe --- /dev/null +++ b/src/geo/include/RAT/GeoNestedSolidArrayFactoryBase.hh @@ -0,0 +1,21 @@ +#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..4544b5a7 --- /dev/null +++ b/src/geo/include/RAT/GeoNestedTubeConstruction.hh @@ -0,0 +1,61 @@ +#ifndef __RAT_GeoNestedTubeConstruction__ +#define __RAT_GeoNestedTubeConstruction__ + +#include +#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, 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 diff --git a/src/geo/src/GeoBuilder.cc b/src/geo/src/GeoBuilder.cc index d28c0e00..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 @@ -44,6 +45,7 @@ GeoBuilder::GeoBuilder() { new GeoReflectorFactory(); new GeoReflectorWaveguideFactory(); new PMTArrayFactory(); + new GeoNestedTubeArrayFactory(); new PMTCoverageFactory(); new GeoWaterBoxArrayFactory(); new GeoBubbleFactory(); diff --git a/src/geo/src/GeoFiberSensitiveDetector.cc b/src/geo/src/GeoFiberSensitiveDetector.cc new file mode 100644 index 00000000..7a8717af --- /dev/null +++ b/src/geo/src/GeoFiberSensitiveDetector.cc @@ -0,0 +1,199 @@ +#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 = "FiberSenDet"); + 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") { + // 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; +} + +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..a150a0fa --- /dev/null +++ b/src/geo/src/GeoFiberSensitiveDetectorHit.cc @@ -0,0 +1,62 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +G4Allocator GeoFiberSensitiveDetectorHitAllocator; + +GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(G4int i, G4double t, G4ThreeVector p) { + id = i; + time = t; + hit_pos = p; + pLogV = 0; +} +GeoFiberSensitiveDetectorHit::~GeoFiberSensitiveDetectorHit() { ; } + +GeoFiberSensitiveDetectorHit::GeoFiberSensitiveDetectorHit(const GeoFiberSensitiveDetectorHit &right) : G4VHit() { + id = right.id; + time = right.time; + pos = right.pos; + hit_pos = right.hit_pos; + rot = right.rot; + pLogV = right.pLogV; +} + +const GeoFiberSensitiveDetectorHit &GeoFiberSensitiveDetectorHit::operator=(const GeoFiberSensitiveDetectorHit &right) { + id = right.id; + time = right.time; + pos = right.pos; + hit_pos = right.hit_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 new file mode 100644 index 00000000..09a4a9d7 --- /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, table); + // 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("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); + + // 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..46ee759a --- /dev/null +++ b/src/geo/src/GeoNestedTubeArrayFactory.cc @@ -0,0 +1,74 @@ +#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..a7df66a1 --- /dev/null +++ b/src/geo/src/GeoNestedTubeConstruction.cc @@ -0,0 +1,181 @@ +#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) { + 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("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("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]; + 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, DBLinkPtr table) { + 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"); + // 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.); + + // 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 diff --git a/src/io/include/RAT/OutNtupleProc.hh b/src/io/include/RAT/OutNtupleProc.hh index 88c725c8..04a80c85 100644 --- a/src/io/include/RAT/OutNtupleProc.hh +++ b/src/io/include/RAT/OutNtupleProc.hh @@ -59,6 +59,7 @@ class OutNtupleProc : public Processor { bool digitizerfits; bool untriggered; bool mchits; + bool nthits; }; NtupleOptions options; @@ -88,6 +89,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; u_int32_t digitizerWindowSize; Double_t digitizerSampleRate; Double_t digitizerDynamicRange; @@ -121,6 +129,14 @@ class OutNtupleProc : public Processor { std::vector mcpmtid; std::vector mcpmtnpe; std::vector mcpmtcharge; + // 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; @@ -150,6 +166,9 @@ class OutNtupleProc : public Processor { std::vector hitPMTID; std::vector hitPMTTime; std::vector hitPMTCharge; + std::vector hitPMTDigitizedTime; + std::vector hitPMTDigitizedCharge; + std::vector hitPMTNCrossings; // Store PMT information from digitized waveform int digitNhits; std::vector digitPeak; @@ -173,6 +192,7 @@ class OutNtupleProc : public Processor { std::vector volumeCodeIndex; std::vector volumeName; + // 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 82666a85..4acc0416 100644 --- a/src/io/src/OutNtupleProc.cc +++ b/src/io/src/OutNtupleProc.cc @@ -54,12 +54,14 @@ OutNtupleProc::OutNtupleProc() : Processor("outntuple") { options.digitizerfits = table->GetZ("include_digitizerfits"); 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; } if (options.digitizerfits) { waveform_fitters = table->GetSArray("waveform_fitters"); @@ -161,6 +163,22 @@ bool OutNtupleProc::OpenFile(std::string filename) { outputTree->Branch(TString("fit_charge_" + fitter_name), &fitCharge[fitter_name]); } } + 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); @@ -215,6 +233,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) { } runBranch = DS::RunStore::GetRun(ds); DS::PMTInfo *pmtinfo = runBranch->GetPMTInfo(); + DS::NestedTubeInfo *ntinfo = runBranch->GetNestedTubeInfo(); const DS::ChannelStatus *channel_status = runBranch->GetChannelStatus(); ULong64_t stonano = 1000000000; dsentries++; @@ -391,6 +410,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++) { @@ -629,6 +665,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 @@ -666,6 +716,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; }