summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/G4chroma.cc230
-rw-r--r--src/G4chroma.hh49
-rw-r--r--src/root.C147
3 files changed, 426 insertions, 0 deletions
diff --git a/src/G4chroma.cc b/src/G4chroma.cc
new file mode 100644
index 0000000..1aba133
--- /dev/null
+++ b/src/G4chroma.cc
@@ -0,0 +1,230 @@
+#include "G4chroma.hh"
+#include <geant4/G4OpticalPhysics.hh>
+#include <geant4/G4EmPenelopePhysics.hh>
+
+ChromaPhysicsList::ChromaPhysicsList(): G4VModularPhysicsList()
+{
+ // default cut value (1.0mm)
+ defaultCutValue = 1.0*mm;
+
+ // General Physics
+ RegisterPhysics( new G4EmPenelopePhysics(0) );
+ // Optical Physics
+ G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
+ RegisterPhysics( opticalPhysics );
+ opticalPhysics->SetTrackSecondariesFirst(false);
+}
+
+ChromaPhysicsList::~ChromaPhysicsList()
+{
+}
+
+void ChromaPhysicsList::SetCuts(){
+ // " G4VUserPhysicsList::SetCutsWithDefault" method sets
+ // the default cut value for all particle types
+ SetCutsWithDefault();
+}
+
+
+PhotonTrackingAction::PhotonTrackingAction()
+{
+}
+
+PhotonTrackingAction::~PhotonTrackingAction()
+{
+}
+
+int PhotonTrackingAction::GetNumPhotons() const
+{
+ return pos.size();
+}
+
+void PhotonTrackingAction::Clear()
+{
+ pos.clear();
+ dir.clear();
+ pol.clear();
+ wavelength.clear();
+ t0.clear();
+}
+
+void PhotonTrackingAction::GetX(double *x) const
+{
+ for (unsigned i=0; i < pos.size(); i++) x[i] = pos[i].x();
+}
+
+void PhotonTrackingAction::GetY(double *y) const
+{
+ for (unsigned i=0; i < pos.size(); i++) y[i] = pos[i].y();
+}
+
+void PhotonTrackingAction::GetZ(double *z) const
+{
+ for (unsigned i=0; i < pos.size(); i++) z[i] = pos[i].z();
+}
+
+void PhotonTrackingAction::GetDirX(double *dir_x) const
+{
+ for (unsigned i=0; i < dir.size(); i++) dir_x[i] = dir[i].x();
+}
+
+void PhotonTrackingAction::GetDirY(double *dir_y) const
+{
+ for (unsigned i=0; i < dir.size(); i++) dir_y[i] = dir[i].y();
+}
+
+void PhotonTrackingAction::GetDirZ(double *dir_z) const
+{
+ for (unsigned i=0; i < dir.size(); i++) dir_z[i] = dir[i].z();
+}
+
+void PhotonTrackingAction::GetPolX(double *pol_x) const
+{
+ for (unsigned i=0; i < pol.size(); i++) pol_x[i] = pol[i].x();
+}
+
+void PhotonTrackingAction::GetPolY(double *pol_y) const
+{
+ for (unsigned i=0; i < pol.size(); i++) pol_y[i] = pol[i].y();
+}
+
+void PhotonTrackingAction::GetPolZ(double *pol_z) const
+{
+ for (unsigned i=0; i < pol.size(); i++) pol_z[i] = pol[i].z();
+}
+
+void PhotonTrackingAction::GetWavelength(double *wl) const
+{
+ for (unsigned i=0; i < wavelength.size(); i++) wl[i] = wavelength[i];
+}
+
+void PhotonTrackingAction::GetT0(double *t) const
+{
+ for (unsigned i=0; i < t0.size(); i++) t[i] = t0[i];
+}
+
+void PhotonTrackingAction::PreUserTrackingAction(const G4Track *track)
+{
+ G4ParticleDefinition *particle = track->GetDefinition();
+ if (particle->GetParticleName() == "opticalphoton") {
+ pos.push_back(track->GetPosition()/m);
+ dir.push_back(track->GetMomentumDirection());
+ pol.push_back(track->GetPolarization());
+ wavelength.push_back( (h_Planck * c_light / track->GetKineticEnergy()) / nanometer );
+ t0.push_back(track->GetGlobalTime() / s);
+ const_cast<G4Track *>(track)->SetTrackStatus(fStopAndKill);
+ }
+}
+
+#include <boost/python.hpp>
+#include <pyublas/numpy.hpp>
+
+using namespace boost::python;
+
+pyublas::numpy_vector<double> PTA_GetX(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetX(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetY(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetY(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetZ(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetZ(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetDirX(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetDirX(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetDirY(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetDirY(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetDirZ(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetDirZ(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetPolX(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetPolX(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetPolY(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetPolY(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetPolZ(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetPolZ(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetWave(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetWavelength(&r[0]);
+ return r;
+}
+
+pyublas::numpy_vector<double> PTA_GetT0(const PhotonTrackingAction *pta)
+{
+ pyublas::numpy_vector<double> r(pta->GetNumPhotons());
+ pta->GetT0(&r[0]);
+ return r;
+}
+
+
+void export_Chroma()
+{
+ class_<ChromaPhysicsList, ChromaPhysicsList*, bases<G4VModularPhysicsList>, boost::noncopyable > ("ChromaPhysicsList", "EM+Optics physics list")
+ .def(init<>())
+ ;
+
+ class_<PhotonTrackingAction, PhotonTrackingAction*, bases<G4UserTrackingAction>,
+ boost::noncopyable > ("PhotonTrackingAction", "Tracking action that saves photons")
+ .def(init<>())
+ .def("GetNumPhotons", &PhotonTrackingAction::GetNumPhotons)
+ .def("Clear", &PhotonTrackingAction::Clear)
+ .def("GetX", PTA_GetX)
+ .def("GetY", PTA_GetY)
+ .def("GetZ", PTA_GetZ)
+ .def("GetDirX", PTA_GetDirX)
+ .def("GetDirY", PTA_GetDirY)
+ .def("GetDirZ", PTA_GetDirZ)
+ .def("GetPolX", PTA_GetPolX)
+ .def("GetPolY", PTA_GetPolY)
+ .def("GetPolZ", PTA_GetPolZ)
+ .def("GetWavelength", PTA_GetWave)
+ .def("GetT0", PTA_GetT0)
+ ;
+}
+
+BOOST_PYTHON_MODULE(G4chroma)
+{
+ export_Chroma();
+}
diff --git a/src/G4chroma.hh b/src/G4chroma.hh
new file mode 100644
index 0000000..4f085aa
--- /dev/null
+++ b/src/G4chroma.hh
@@ -0,0 +1,49 @@
+#ifndef __G4chroma_hh__
+#define __G4chroma_hh__
+
+#include <geant4/G4VModularPhysicsList.hh>
+class ChromaPhysicsList: public G4VModularPhysicsList
+{
+public:
+ ChromaPhysicsList();
+ virtual ~ChromaPhysicsList();
+ virtual void SetCuts();
+};
+
+#include <geant4/G4UserTrackingAction.hh>
+#include <vector>
+#include <geant4/G4ThreeVector.hh>
+
+class PhotonTrackingAction : public G4UserTrackingAction
+{
+public:
+ PhotonTrackingAction();
+ virtual ~PhotonTrackingAction();
+
+ int GetNumPhotons() const;
+ void Clear();
+
+ void GetX(double *x) const;
+ void GetY(double *y) const;
+ void GetZ(double *z) const;
+ void GetDirX(double *dir_x) const;
+ void GetDirY(double *dir_y) const;
+ void GetDirZ(double *dir_z) const;
+ void GetPolX(double *pol_x) const;
+ void GetPolY(double *pol_y) const;
+ void GetPolZ(double *pol_z) const;
+
+ void GetWavelength(double *wl) const;
+ void GetT0(double *t) const;
+
+ virtual void PreUserTrackingAction(const G4Track *);
+
+protected:
+ std::vector<G4ThreeVector> pos;
+ std::vector<G4ThreeVector> dir;
+ std::vector<G4ThreeVector> pol;
+ std::vector<double> wavelength;
+ std::vector<double> t0;
+};
+
+#endif
diff --git a/src/root.C b/src/root.C
new file mode 100644
index 0000000..2b39b1b
--- /dev/null
+++ b/src/root.C
@@ -0,0 +1,147 @@
+#include <TVector3.h>
+#include <vector>
+#include <TTree.h>
+#include <string>
+
+struct Vertex {
+ std::string particle_name;
+ TVector3 pos;
+ TVector3 dir;
+ TVector3 pol;
+ double ke;
+ double t0;
+
+ ClassDef(Vertex, 1);
+};
+
+struct Photon {
+ double t;
+ TVector3 pos;
+ TVector3 dir;
+ TVector3 pol;
+ double wavelength; // nm
+ unsigned int flag;
+ int last_hit_triangle;
+
+ ClassDef(Photon, 1);
+};
+
+struct Channel {
+ Channel() : id(-1), t(-1e9), q(-1e9) { };
+ int id;
+ double t;
+ double q;
+ unsigned int flag;
+
+ ClassDef(Channel, 1);
+};
+
+struct Event {
+ int id;
+ unsigned int nhit;
+ unsigned int nchannels;
+
+ Vertex primary_vertex;
+
+ std::vector<Vertex> vertices;
+ std::vector<Photon> photons_beg;
+ std::vector<Photon> photons_end;
+ std::vector<Channel> channels;
+
+ ClassDef(Event, 1);
+};
+
+void fill_channels(Event *ev, unsigned int nhit, unsigned int *ids, float *t,
+ float *q, unsigned int *flags, unsigned int nchannels)
+{
+ ev->nhit = 0;
+ ev->nchannels = nchannels;
+ ev->channels.resize(0);
+
+ Channel ch;
+ unsigned int id;
+ for (unsigned int i=0; i < nhit; i++) {
+ ev->nhit++;
+ id = ids[i];
+ ch.id = id;
+ ch.t = t[id];
+ ch.q = q[id];
+ ch.flag = flags[id];
+ ev->channels.push_back(ch);
+ }
+}
+
+void get_channels(Event *ev, int *hit, float *t, float *q, unsigned int *flags)
+{
+ for (unsigned int i=0; i < ev->nchannels; i++) {
+ hit[i] = 0;
+ t[i] = -1e9f;
+ q[i] = -1e9f;
+ flags[i] = 0;
+ }
+
+ unsigned int id;
+ for (unsigned int i=0; i < ev->channels.size(); i++) {
+ id = ev->channels[i].id;
+
+ if (id < ev->nchannels) {
+ hit[id] = 1;
+ t[id] = ev->channels[i].t;
+ q[id] = ev->channels[i].q;
+ flags[id] = ev->channels[i].flag;
+ }
+ }
+}
+
+void get_photons(const std::vector<Photon> &photons, float *pos, float *dir,
+ float *pol, float *wavelengths, float *t,
+ int *last_hit_triangles, unsigned int *flags)
+{
+ for (unsigned int i=0; i < photons.size(); i++) {
+ const Photon &photon = photons[i];
+ pos[3*i] = photon.pos.X();
+ pos[3*i+1] = photon.pos.Y();
+ pos[3*i+2] = photon.pos.Z();
+
+ dir[3*i] = photon.dir.X();
+ dir[3*i+1] = photon.dir.Y();
+ dir[3*i+2] = photon.dir.Z();
+
+ pol[3*i] = photon.pol.X();
+ pol[3*i+1] = photon.pol.Y();
+ pol[3*i+2] = photon.pol.Z();
+
+ wavelengths[i] = photon.wavelength;
+ t[i] = photon.t;
+ flags[i] = photon.flag;
+ last_hit_triangles[i] = photon.last_hit_triangle;
+ }
+}
+
+void fill_photons(std::vector<Photon> &photons,
+ unsigned int nphotons, float *pos, float *dir,
+ float *pol, float *wavelengths, float *t,
+ int *last_hit_triangles, unsigned int *flags)
+{
+ photons.resize(nphotons);
+
+ for (unsigned int i=0; i < nphotons; i++) {
+ Photon &photon = photons[i];
+ photon.t = t[i];
+ photon.pos.SetXYZ(pos[3*i], pos[3*i + 1], pos[3*i + 2]);
+ photon.dir.SetXYZ(dir[3*i], dir[3*i + 1], dir[3*i + 2]);
+ photon.pol.SetXYZ(pol[3*i], pol[3*i + 1], pol[3*i + 2]);
+ photon.wavelength = wavelengths[i];
+ photon.last_hit_triangle = last_hit_triangles[i];
+ photon.flag = flags[i];
+
+ }
+}
+
+#ifdef __MAKECINT__
+#pragma link C++ class vector<Vertex>;
+#pragma link C++ class vector<Photon>;
+#pragma link C++ class vector<Channel>;
+#endif
+
+
generated by cgit (git 2.34.1) at 2025-05-07 14:59:30 +0000