diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx
index f8a1933c8f7c8712f0598ce13ce3dfc5042719de..0d3385c78ab1f0a26f92b00863484e2bbf2a670e 100644
--- a/NPLib/Core/NPDetectorManager.cxx
+++ b/NPLib/Core/NPDetectorManager.cxx
@@ -268,7 +268,7 @@ void NPL::DetectorManager::BuildPhysicalEvent(){
(it->second->*m_ClearEventPhysicsPtr)();
(it->second->*m_BuildPhysicalPtr)();
(it->second->*m_FillSpectra)();
- if(m_CheckSpectra)
+ if(m_CheckSpectra) {
(it->second->*m_CheckSpectra)();
}
}
diff --git a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h
index e63650cc71b1fbe501339488bfa1790c5463c32b..ae07f7afbbb134c6f5df82bd5d8bcdf1690e3697 100644
--- a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h
+++ b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h
@@ -57,18 +57,18 @@ class TComptonTelescopePhysics : public TObject, public NPL::VDetector
// DSSD
int EventMultiplicity;
vector<int> EventType;
- //vector<int> TowerNumber;
+ vector<int> TowerNumber;//
vector<int> DetectorNumber;
vector<int> Strip_Front;
vector<int> Strip_Back;
- //vector<double> Strip_E;
+ vector<double> Strip_E;//
vector<double> Strip_T;
vector<double> Front_Energy;
vector<double> Back_Energy;
- //vector<double> Half_Energy;
+ vector<double> Half_Energy;//
vector<double> Front_Time;
vector<double> Back_Time;
- //vector<bool> Same_FBTime;
+ vector<bool> Same_FBTime;//
// Calorimeter
vector<double> Calor_E;
vector<double> Calor_T;
@@ -240,8 +240,8 @@ class TComptonTelescopePhysics : public TObject, public NPL::VDetector
Int_t m_CounterEvt[50]; //!
Int_t m_CounterHit[50]; //!
// physical events
- vector<int> TowerNumber;
- vector<double> Half_Energy;
+ //vector<int> TowerNumber;
+ //vector<double> Half_Energy;
//vector<bool> Same_FBTime;
};
diff --git a/NPSimulation/Detectors/Strasse/CADMesh.hh b/NPSimulation/Detectors/Strasse/CADMesh.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1f2db2a16b14c3e2f1794fcdbbb502bdfa39fd1c
--- /dev/null
+++ b/NPSimulation/Detectors/Strasse/CADMesh.hh
@@ -0,0 +1,2704 @@
+// The MIT License (MIT)
+//
+// Copyright (c) 2011-2020 Christopher M. Poole <mail@christopherpoole.net>
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+// CADMESH - Load CAD files into Geant4 quickly and easily.
+//
+// Read all about it at: https://github.com/christopherpoole/CADMesh
+//
+// Basic usage:
+//
+// #include "CADMesh.hh" // This file.
+// ...
+// auto mesh = CADMesh::TessellatedMesh::FromSTL("mesh.stl");
+// G4VSolid* solid = mesh->GetSolid();
+// ...
+//
+
+// !! THIS FILE HAS BEEN AUTOMATICALLY GENERATED. DON'T EDIT IT DIRECTLY. !!
+
+#pragma once
+
+class BuiltInReader;
+class CADMeshTemplate;
+class Mesh;
+class STLReader;
+class ASSIMPReader;
+class PLYReader;
+class FileTypes;
+class OBJReader;
+class Exceptions;
+class TetrahedralMesh;
+class Reader;
+class Lexer;
+class LexerMacros;
+class TessellatedMesh;
+
+#include "G4String.hh"
+
+#include <algorithm>
+#include <map>
+
+namespace CADMesh {
+
+namespace File {
+
+enum Type {
+ Unknown,
+ PLY,
+ STL,
+ DAE,
+ OBJ,
+ TET,
+ OFF,
+};
+
+static std::map<Type, G4String> Extension = {
+ {Unknown, "unknown"}, {PLY, "ply"}, {STL, "stl"}, {DAE, "dae"},
+ {OBJ, "obj"}, {TET, "tet"}, {OFF, "off"}};
+
+static std::map<Type, G4String> TypeString = {
+ {Unknown, "UNKNOWN"}, {PLY, "PLY"}, {STL, "STL"}, {DAE, "DAE"},
+ {OBJ, "OBJ"}, {TET, "TET"}, {OFF, "OFF"}};
+
+static std::map<Type, G4String> TypeName = {
+ {Unknown, "Unknown File Format"}, {PLY, "Stanford Triangle Format (PLY)"},
+ {STL, "Stereolithography (STL)"}, {DAE, "COLLADA (DAE)"},
+ {OBJ, "Wavefront (OBJ)"}, {TET, "TetGet (TET)"},
+ {OFF, "Object File Format (OFF)"}};
+
+Type TypeFromExtension(G4String extension);
+Type TypeFromName(G4String name);
+}
+}
+
+#include "G4ThreeVector.hh"
+#include "G4TriangularFacet.hh"
+
+#include <memory>
+#include <vector>
+
+namespace CADMesh {
+
+typedef std::vector<G4ThreeVector> Points;
+typedef std::vector<G4TriangularFacet *> Triangles;
+
+class Mesh {
+public:
+ Mesh(Points points, Triangles triangles, G4String name = "");
+
+ static std::shared_ptr<Mesh> New(Points points, Triangles triangles,
+ G4String name = "");
+
+ static std::shared_ptr<Mesh> New(Triangles triangles, G4String name = "");
+
+ static std::shared_ptr<Mesh> New(std::shared_ptr<Mesh> mesh,
+ G4String name = "");
+
+public:
+ G4String GetName();
+ Points GetPoints();
+ Triangles GetTriangles();
+
+ G4bool IsValidForNavigation();
+
+private:
+ G4String name_ = "";
+
+ Points points_;
+ Triangles triangles_;
+};
+
+typedef std::vector<std::shared_ptr<Mesh>> Meshes;
+}
+
+namespace CADMesh {
+
+namespace File {
+
+class Reader {
+public:
+ Reader(G4String reader_name);
+ ~Reader();
+
+ virtual G4bool Read(G4String filepath) = 0;
+ virtual G4bool CanRead(Type file_type) = 0;
+
+public:
+ G4String GetName();
+
+ std::shared_ptr<Mesh> GetMesh();
+ std::shared_ptr<Mesh> GetMesh(size_t index);
+ std::shared_ptr<Mesh> GetMesh(G4String name, G4bool exact = true);
+
+ size_t GetNumberOfMeshes();
+
+ Meshes GetMeshes();
+
+protected:
+ size_t AddMesh(std::shared_ptr<Mesh> mesh);
+ void SetMeshes(Meshes meshs);
+
+private:
+ Meshes meshes_;
+
+ G4String name_ = "";
+};
+}
+}
+
+#include <iostream>
+#include <string>
+
+namespace CADMesh {
+
+namespace File {
+
+struct Token {
+ std::string name;
+
+ bool operator==(Token other) { return name == other.name; };
+ bool operator!=(Token other) { return name != other.name; };
+};
+
+static Token ErrorToken{"ErrorToken"};
+static Token EnfOfFileToken{"EndOfFileToken"};
+
+static Token ParentToken{"Parent"};
+static Token WordToken{"Word"};
+static Token NumberToken{"Number"};
+static Token ThreeVectorToken{"ThreeVector"};
+static Token VertexToken{"Vertex"};
+static Token VerticesToken{"Vertices"};
+static Token FacetToken{"Facet"};
+static Token LineToken{"Line"};
+static Token NormalToken{"Normal"};
+static Token TextureToken{"Texture"};
+static Token SolidToken{"Solid"};
+static Token ObjectToken{"Object"};
+static Token CommentToken{"Comment"};
+static Token BlankLineToken{"BlankLine"};
+
+struct Item {
+ Token token;
+
+ size_t position;
+ size_t line;
+
+ std::string value;
+ std::string error;
+
+ Item *parent;
+
+ std::vector<Item> children;
+};
+
+typedef std::vector<Item> Items;
+typedef Items::iterator ItemsIterator;
+
+class Lexer;
+
+struct State {
+ virtual State *operator()(Lexer *) const = 0;
+};
+
+struct __FinalState : public State {
+ State *operator()(Lexer *) const { return nullptr; }
+};
+
+class Lexer {
+public:
+ Lexer(std::string filepath, State *initial_state = nullptr);
+
+public:
+ std::string String();
+
+ void Run(State *initial_state, size_t lines = 0);
+ Items GetItems();
+
+ void Backup();
+ void BackupTo(int position);
+
+ std::string Next();
+ std::string Peek();
+
+ void Skip();
+
+ Item *ThisIsA(Token token, std::string error = "");
+ Item *StartOfA(Token token, std::string error = "");
+ Item *EndOfA(Token token, std::string error = "");
+ Item *MaybeEndOfA(Token token, std::string error = "");
+
+ bool OneOf(std::string possibles);
+ bool ManyOf(std::string possibles);
+ bool Until(std::string match);
+ bool MatchExactly(std::string match);
+
+ bool OneDigit();
+ bool ManyDigits();
+
+ bool OneLetter();
+ bool ManyLetters();
+
+ bool ManyCharacters();
+
+ bool Integer();
+ bool Float();
+ bool Number();
+
+ bool SkipWhiteSpace();
+ bool SkipLineBreak();
+ bool SkipLineBreaks();
+ bool SkipLine();
+
+ State *Error(std::string message);
+ State *LastError();
+
+ bool TestState(State *state);
+
+ bool IsDryRun();
+
+ void PrintMessage(std::string name, std::string message);
+ void PrintItem(Item item);
+
+ size_t LineNumber();
+
+private:
+ State *state_;
+
+ Item *parent_item_ = nullptr;
+ Items items_;
+
+ std::string input_;
+
+ size_t position_ = 0;
+ size_t start_ = 0;
+ size_t width_ = 1;
+ size_t line_ = 1;
+ size_t end_line_ = 0;
+
+ bool dry_run_ = false;
+
+ int depth_ = 0;
+
+ std::string last_error_ = "";
+};
+}
+}
+
+#ifdef USE_CADMESH_ASSIMP_READER
+
+#include "assimp/Importer.hpp"
+#include "assimp/postprocess.h"
+#include "assimp/scene.h"
+
+namespace CADMesh {
+
+namespace File {
+
+class ASSIMPReader : public File::Reader {
+public:
+ ASSIMPReader();
+ ~ASSIMPReader();
+
+public:
+ G4bool Read(G4String filepath);
+ G4bool CanRead(Type file_type);
+
+private:
+ Assimp::Importer *importer_;
+};
+
+std::shared_ptr<ASSIMPReader> ASSIMP();
+}
+}
+#endif
+
+namespace CADMesh {
+
+namespace File {
+
+class BuiltInReader : public Reader {
+public:
+ BuiltInReader() : Reader("BuiltInReader"){};
+
+public:
+ G4bool Read(G4String filepath);
+ G4bool CanRead(File::Type file_type);
+};
+
+std::shared_ptr<BuiltInReader> BuiltIn();
+}
+}
+#ifndef CADMESH_DEFAULT_READER
+#define CADMESH_DEFAULT_READER BuiltIn
+#endif
+
+#include "G4AssemblyVolume.hh"
+#include "G4LogicalVolume.hh"
+#include "G4Material.hh"
+#include "G4TessellatedSolid.hh"
+#include "G4Tet.hh"
+#include "G4UIcommand.hh"
+
+namespace CADMesh {
+
+template <typename T> class CADMeshTemplate {
+public:
+ CADMeshTemplate();
+
+ CADMeshTemplate(G4String file_name);
+
+ CADMeshTemplate(G4String file_name, File::Type file_type);
+
+ CADMeshTemplate(std::shared_ptr<File::Reader> reader);
+
+ CADMeshTemplate(G4String file_name, std::shared_ptr<File::Reader> reader);
+
+ CADMeshTemplate(G4String file_name, File::Type file_type,
+ std::shared_ptr<File::Reader> reader);
+
+ static std::shared_ptr<T> From(G4String file_name);
+
+ static std::shared_ptr<T> From(G4String file_name,
+ std::shared_ptr<File::Reader> reader);
+
+ static std::shared_ptr<T> FromPLY(G4String file_name);
+
+ static std::shared_ptr<T> FromPLY(G4String file_name,
+ std::shared_ptr<File::Reader> reader);
+
+ static std::shared_ptr<T> FromSTL(G4String file_name);
+
+ static std::shared_ptr<T> FromSTL(G4String file_name,
+ std::shared_ptr<File::Reader> reader);
+
+ static std::shared_ptr<T> FromOBJ(G4String file_name);
+
+ static std::shared_ptr<T> FromOBJ(G4String file_name,
+ std::shared_ptr<File::Reader> reader);
+
+ ~CADMeshTemplate();
+
+public:
+ virtual G4VSolid *GetSolid() = 0;
+ virtual G4VSolid *GetSolid(G4int index) = 0;
+ virtual G4VSolid *GetSolid(G4String name, G4bool exact = true) = 0;
+
+ virtual std::vector<G4VSolid *> GetSolids() = 0;
+
+ virtual G4AssemblyVolume *GetAssembly() = 0;
+
+ bool IsValidForNavigation();
+
+public:
+ G4String GetFileName();
+
+ File::Type GetFileType();
+
+ void SetVerbose(G4int verbose);
+ G4int GetVerbose();
+
+ void SetScale(G4double scale);
+ G4double GetScale();
+
+ void SetOffset(G4double x, G4double y, G4double z);
+ void SetOffset(G4ThreeVector offset);
+ G4ThreeVector GetOffset();
+
+protected:
+ G4String file_name_;
+ File::Type file_type_;
+ G4int verbose_;
+ G4double scale_;
+ G4ThreeVector offset_;
+
+ G4AssemblyVolume *assembly_ = nullptr;
+
+ std::shared_ptr<File::Reader> reader_ = nullptr;
+};
+}
+
+#include "globals.hh"
+
+namespace CADMesh {
+
+namespace Exceptions {
+
+void FileNotFound(G4String origin, G4String filepath);
+
+void LexerError(G4String origin, G4String message);
+
+void ParserError(G4String origin, G4String message);
+
+void ReaderCantReadError(G4String origin, File::Type file_type,
+ G4String filepath);
+
+void MeshNotFound(G4String origin, size_t index);
+void MeshNotFound(G4String origin, G4String name);
+}
+}
+
+namespace CADMesh {
+
+class TessellatedMesh : public CADMeshTemplate<TessellatedMesh> {
+ using CADMeshTemplate::CADMeshTemplate;
+
+public:
+ G4VSolid *GetSolid();
+ G4VSolid *GetSolid(G4int index);
+ G4VSolid *GetSolid(G4String name, G4bool exact = true);
+
+ std::vector<G4VSolid *> GetSolids();
+
+ G4TessellatedSolid *GetTessellatedSolid();
+ G4TessellatedSolid *GetTessellatedSolid(G4int index);
+ G4TessellatedSolid *GetTessellatedSolid(G4String name, G4bool exact = true);
+ G4TessellatedSolid *GetTessellatedSolid(std::shared_ptr<Mesh> mesh);
+
+ G4AssemblyVolume *GetAssembly();
+
+public:
+ void SetReverse(G4bool reverse) { this->reverse_ = reverse; };
+
+ G4bool GetReverse() { return this->reverse_; };
+
+private:
+ G4bool reverse_;
+};
+}
+
+#ifdef USE_CADMESH_TETGEN
+
+#include "tetgen.h"
+
+namespace CADMesh {
+
+class TetrahedralMesh : public CADMeshTemplate<TetrahedralMesh> {
+public:
+ using CADMeshTemplate::CADMeshTemplate;
+
+ TetrahedralMesh();
+ ~TetrahedralMesh();
+
+public:
+ G4VSolid *GetSolid();
+ G4VSolid *GetSolid(G4int index);
+ G4VSolid *GetSolid(G4String name, G4bool exact = true);
+
+ std::vector<G4VSolid *> GetSolids();
+
+ G4AssemblyVolume *GetAssembly();
+
+public:
+ void SetMaterial(G4Material *material) { this->material_ = material; };
+
+ G4Material *GetMaterial() { return this->material_; };
+
+ void SetQuality(G4double quality) { this->quality_ = quality; };
+
+ G4double GetQuality() { return this->quality_; };
+
+ std::shared_ptr<tetgenio> GetTetgenInput() { return in_; };
+
+ std::shared_ptr<tetgenio> GetTetgenOutput() { return out_; };
+
+private:
+ G4ThreeVector GetTetPoint(G4int index_offset);
+
+private:
+ std::shared_ptr<tetgenio> in_ = nullptr;
+ std::shared_ptr<tetgenio> out_ = nullptr;
+
+ G4double quality_;
+
+ G4Material *material_;
+};
+}
+#endif
+
+namespace CADMesh {
+
+namespace File {
+
+inline Type TypeFromExtension(G4String extension) {
+ std::for_each(extension.begin(), extension.end(),
+ [](char &e) { e = ::tolower(e); });
+
+ for (auto e = Extension.begin(); e != Extension.end(); ++e) {
+ if (e->second == extension) {
+ return e->first;
+ }
+ }
+
+ return Unknown;
+}
+
+inline Type TypeFromName(G4String name) {
+ auto extension = name.substr(name.find_last_of(".") + 1);
+
+ return TypeFromExtension(extension);
+}
+}
+}
+
+namespace CADMesh {
+
+inline Mesh::Mesh(Points points, Triangles triangles, G4String name)
+ : name_(name), points_(points), triangles_(triangles) {}
+
+inline std::shared_ptr<Mesh> Mesh::New(Points points, Triangles triangles,
+ G4String name) {
+ return std::make_shared<Mesh>(points, triangles, name);
+}
+
+inline std::shared_ptr<Mesh> Mesh::New(Triangles triangles, G4String name) {
+ Points points;
+
+ return New(points, triangles, name);
+}
+
+inline std::shared_ptr<Mesh> Mesh::New(std::shared_ptr<Mesh> mesh,
+ G4String name) {
+ return New(mesh->GetPoints(), mesh->GetTriangles(), name);
+}
+
+inline G4String Mesh::GetName() { return name_; }
+
+inline Points Mesh::GetPoints() { return points_; }
+
+inline Triangles Mesh::GetTriangles() { return triangles_; }
+
+inline G4bool Mesh::IsValidForNavigation() {
+ std::map<G4ThreeVector, size_t> point_index;
+
+ size_t index = 0;
+ for (auto point : points_) {
+ point_index[point] = index;
+ index++;
+ }
+
+ typedef std::pair<size_t, size_t> Edge;
+ std::map<Edge, G4int> edge_use_count;
+
+ for (size_t i = 0; i < triangles_.size(); i++) {
+ auto triangle = triangles_[i];
+
+ size_t a = point_index[triangle->GetVertex(0)];
+ size_t b = point_index[triangle->GetVertex(1)];
+ size_t c = point_index[triangle->GetVertex(2)];
+
+ if (a < b) {
+ edge_use_count[Edge(a, b)] += 1;
+ }
+
+ else if (a > b) {
+ edge_use_count[Edge(b, a)] += 1;
+ }
+
+ if (b < c) {
+ edge_use_count[Edge(b, c)] += 1;
+ }
+
+ else if (b > c) {
+ edge_use_count[Edge(c, b)] += 1;
+ }
+
+ if (c < a) {
+ edge_use_count[Edge(c, a)] += 1;
+ }
+
+ else if (c > a) {
+ edge_use_count[Edge(a, c)] += 1;
+ }
+ }
+
+ for (auto count : edge_use_count) {
+ if (count.second != 2) {
+ return false;
+ }
+ }
+
+ return true;
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+inline Reader::Reader(G4String reader_name) : name_(reader_name) {}
+
+inline Reader::~Reader() {}
+
+inline G4String Reader::GetName() { return name_; }
+
+inline std::shared_ptr<Mesh> Reader::GetMesh() {
+ if (meshes_.size() > 0) {
+ return meshes_[0];
+ }
+
+ return nullptr;
+}
+
+inline std::shared_ptr<Mesh> Reader::GetMesh(size_t index) {
+ if (index < meshes_.size()) {
+ return meshes_[index];
+ }
+
+ Exceptions::MeshNotFound("Reader::GetMesh", index);
+
+ return nullptr;
+}
+
+inline std::shared_ptr<Mesh> Reader::GetMesh(G4String name, G4bool exact) {
+ for (auto mesh : meshes_) {
+ if (exact) {
+ if (mesh->GetName() == name)
+ return mesh;
+ }
+
+ else {
+ if (mesh->GetName().find(name) != std::string::npos)
+ return mesh;
+ }
+ }
+
+ Exceptions::MeshNotFound("Reader::GetMesh", name);
+
+ return nullptr;
+}
+
+inline Meshes Reader::GetMeshes() { return meshes_; }
+
+inline size_t Reader::GetNumberOfMeshes() { return meshes_.size(); }
+
+inline size_t Reader::AddMesh(std::shared_ptr<Mesh> mesh) {
+ meshes_.push_back(mesh);
+
+ return meshes_.size();
+}
+
+inline void Reader::SetMeshes(Meshes meshes) { meshes_ = meshes; }
+}
+}
+
+#include <fstream>
+#include <sstream>
+#include <streambuf>
+
+namespace CADMesh {
+
+namespace File {
+
+inline Lexer::Lexer(std::string filepath, State *initial_state) {
+ std::ifstream file(filepath);
+ input_ = std::string((std::istreambuf_iterator<char>(file)),
+ std::istreambuf_iterator<char>());
+
+ if (initial_state) {
+ Run(initial_state);
+ }
+}
+
+inline std::string Lexer::String() {
+ return input_.substr(start_, position_ - start_);
+}
+
+inline void Lexer::Run(State *initial_state, size_t lines) {
+ parent_item_ = new Item{ParentToken, position_, line_, "", "",
+ nullptr, std::vector<Item>()};
+
+ state_ = initial_state;
+
+ end_line_ = 0;
+
+ if (lines > 0)
+ end_line_ = line_ + lines;
+
+ while (state_) {
+ state_ = (*state_)(this);
+ }
+}
+
+inline Items Lexer::GetItems() { return parent_item_->children; }
+
+inline void Lexer::Backup() {
+ position_ -= width_;
+
+ auto next = input_.substr(position_, 1);
+
+ if (next == "\n") {
+ line_--;
+ }
+}
+
+inline void Lexer::BackupTo(int position) {
+ auto s = input_.substr(position, position_ - position);
+ line_ -= std::count(s.begin(), s.end(), '\n');
+
+ position_ = position;
+}
+
+inline std::string Lexer::Next() {
+ if (position_ >= input_.length()) {
+ return "";
+ }
+
+ auto next = input_.substr(position_, 1);
+
+ width_ = 1;
+ position_ += width_;
+
+ if (next == "\n")
+ line_++;
+
+ return next;
+}
+
+inline std::string Lexer::Peek() {
+ auto next = Next();
+
+ if (next != "")
+ Backup();
+
+ return next;
+}
+
+inline void Lexer::Skip() { start_ = position_; }
+
+inline Item *Lexer::ThisIsA(Token token, std::string error) {
+ if (dry_run_)
+ return nullptr;
+
+ auto item =
+ Item{token, position_, line_, String(), error, parent_item_, Items()};
+ Skip();
+
+ if (parent_item_) {
+ PrintItem(item);
+
+ parent_item_->children.push_back(item);
+ return &(parent_item_->children.back());
+ }
+
+ else {
+ depth_++;
+ PrintItem(item);
+
+ items_.push_back(item);
+ return &(items_.back());
+ }
+}
+
+inline Item *Lexer::StartOfA(Token token, std::string error) {
+ if (dry_run_)
+ return nullptr;
+
+ parent_item_ = ThisIsA(token, error);
+
+ depth_++;
+
+ return parent_item_;
+}
+
+inline Item *Lexer::EndOfA(Token token, std::string /*error*/) {
+ if (dry_run_)
+ return nullptr;
+
+ depth_--;
+
+ PrintItem(*parent_item_);
+
+ if (parent_item_->token.name != token.name) {
+ Exceptions::LexerError("Lexer::EndOfA",
+ "Trying to end a '" + parent_item_->token.name +
+ "' with a '" + token.name + "' token.");
+ }
+
+ if (parent_item_) {
+ parent_item_ = parent_item_->parent;
+ }
+
+ return nullptr;
+}
+
+inline Item *Lexer::MaybeEndOfA(Token token, std::string error) {
+ if (parent_item_->token.name == token.name) {
+ return EndOfA(token, error);
+ }
+
+ else {
+ return nullptr;
+ }
+}
+
+inline bool Lexer::OneOf(std::string possibles) {
+ auto peek = Peek();
+
+ size_t position = possibles.find(Peek());
+
+ if (position != std::string::npos) {
+ auto next = Next();
+ return next != "";
+ }
+
+ return false;
+}
+
+inline bool Lexer::ManyOf(std::string possibles) {
+ bool has = false;
+
+ while (OneOf(possibles)) {
+ has = true;
+ }
+
+ return has;
+}
+
+inline bool Lexer::Until(std::string match) {
+ while (!OneOf(match)) {
+ if (Next() == "")
+ return false;
+ }
+
+ return true;
+}
+
+inline bool Lexer::MatchExactly(std::string match) {
+ auto start_position = position_;
+
+ for (auto m : match) {
+ if (!OneOf(std::string(1, m))) {
+ BackupTo(start_position);
+ return false;
+ }
+ }
+
+ return true;
+}
+
+inline bool Lexer::OneDigit() { return OneOf("0123456789"); }
+
+inline bool Lexer::ManyDigits() { return ManyOf("0123456789"); }
+
+inline bool Lexer::OneLetter() {
+ return OneOf("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ");
+}
+
+inline bool Lexer::ManyLetters() {
+ return ManyOf("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ");
+}
+
+inline bool Lexer::ManyCharacters() {
+ return ManyOf("!\"#$%&\\\'()*+,-./"
+ "0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[|]^_`"
+ "abcdefghijklmnopqrstuvwxyz{}~");
+}
+
+inline bool Lexer::Integer() {
+ auto start_position = position_;
+
+ OneOf("+-");
+
+ if (!ManyDigits()) {
+ BackupTo(start_position);
+ return false;
+ }
+
+ return true;
+}
+
+inline bool Lexer::Float() {
+ auto start_position = position_;
+
+ OneOf("+-");
+
+ bool has_integer = Integer();
+
+ if (!OneOf(".")) {
+ BackupTo(start_position);
+ return false;
+ }
+
+ bool has_decimal = ManyDigits();
+
+ if (!has_integer || !has_decimal) {
+ BackupTo(start_position);
+ return false;
+ }
+
+ return true;
+}
+
+inline bool Lexer::Number() {
+ if (!Float()) {
+ if (!Integer()) {
+ return false;
+ }
+ }
+
+ auto start_position = position_;
+
+ if (OneOf("eE")) {
+ if (!Float()) {
+ if (!Integer()) {
+ BackupTo(start_position);
+ }
+ }
+ }
+
+ return true;
+}
+
+inline bool Lexer::SkipWhiteSpace() {
+ if (!ManyOf(" \t\r")) {
+ Skip();
+ return false;
+ }
+
+ Skip();
+ return true;
+}
+
+inline bool Lexer::SkipLineBreak() {
+ if (!OneOf("\n")) {
+ return false;
+ }
+
+ Skip();
+ return true;
+}
+
+inline bool Lexer::SkipLineBreaks() {
+ if (!ManyOf("\n")) {
+ return false;
+ }
+
+ Skip();
+ return true;
+}
+
+inline bool Lexer::SkipLine() {
+ if (!Until("\n")) {
+ return false;
+ }
+
+ Skip();
+ return true;
+}
+
+inline State *Lexer::Error(std::string message) {
+
+ std::stringstream error;
+ error << "Error around line " << line_ << ": " << message << std::endl;
+
+ last_error_ = error.str();
+
+ if (dry_run_)
+ return nullptr;
+
+ Item item{ErrorToken, position_, line_, "",
+ error.str(), parent_item_, Items()};
+ items_.push_back(item);
+
+ Exceptions::LexerError("Lexer", error.str());
+
+ return nullptr;
+}
+
+inline State *Lexer::LastError() {
+ if (last_error_ == "") {
+ Exceptions::LexerError("Lexer", "Something went wrong.");
+ }
+
+ else {
+ Exceptions::LexerError("Lexer", last_error_);
+ }
+
+ return nullptr;
+}
+
+inline bool Lexer::TestState(State *state) {
+ if (dry_run_)
+ return false;
+
+ if (end_line_ > 0 && line_ > end_line_)
+ return false;
+
+ auto start_position = position_;
+
+ dry_run_ = true;
+
+ bool state_transition = ((*state)(this) != nullptr);
+
+ dry_run_ = false;
+
+ BackupTo(start_position);
+
+ return state_transition;
+}
+
+inline bool Lexer::IsDryRun() { return dry_run_; }
+
+inline size_t Lexer::LineNumber() { return line_; }
+
+#ifdef CADMESH_LEXER_VERBOSE
+inline void Lexer::PrintMessage(std::string name, std::string message) {
+ std::cout << "Lexer::" << name << " : " << message << std::endl;
+}
+#else
+inline void Lexer::PrintMessage(std::string, std::string) {}
+#endif
+
+#ifdef CADMESH_LEXER_VERBOSE
+inline void Lexer::PrintItem(Item item) {
+ auto depth = std::max(0, depth_) * 2;
+ std::cout << std::string(depth, ' ') << item.token.name << ": " << item.value
+ << std::endl;
+}
+#else
+inline void Lexer::PrintItem(Item) {}
+#endif
+}
+}
+
+namespace CADMesh {
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate() : CADMeshTemplate<T>("") {}
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate(G4String file_name)
+ : CADMeshTemplate<T>(file_name, File::Unknown) {}
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate(G4String file_name, File::Type file_type)
+ : CADMeshTemplate<T>(file_name, file_type, File::CADMESH_DEFAULT_READER()) {
+}
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate(std::shared_ptr<File::Reader> reader)
+ : CADMeshTemplate<T>("", reader) {}
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate(G4String file_name,
+ std::shared_ptr<File::Reader> reader)
+ : CADMeshTemplate<T>(file_name, File::Unknown, reader) {}
+
+template <typename T>
+CADMeshTemplate<T>::CADMeshTemplate(G4String file_name, File::Type file_type,
+ std::shared_ptr<File::Reader> reader) {
+ if (!reader->CanRead(file_type)) {
+ Exceptions::ReaderCantReadError(reader->GetName(), file_type, file_name);
+ }
+
+ file_name_ = file_name;
+ file_type_ = file_type;
+
+ scale_ = 1.0;
+
+ offset_ = G4ThreeVector();
+
+ verbose_ = 0;
+
+ reader_ = reader;
+
+ reader_->Read(file_name_);
+}
+
+template <typename T>
+std::shared_ptr<T> CADMeshTemplate<T>::From(G4String file_name) {
+ return std::make_shared<T>(file_name);
+}
+
+template <typename T>
+std::shared_ptr<T>
+CADMeshTemplate<T>::From(G4String file_name,
+ std::shared_ptr<File::Reader> reader) {
+ return std::make_shared<T>(file_name, reader);
+}
+
+template <typename T>
+std::shared_ptr<T> CADMeshTemplate<T>::FromPLY(G4String file_name) {
+ return std::make_shared<T>(file_name, File::PLY);
+}
+
+template <typename T>
+std::shared_ptr<T>
+CADMeshTemplate<T>::FromPLY(G4String file_name,
+ std::shared_ptr<File::Reader> reader) {
+ return std::make_shared<T>(file_name, File::PLY, reader);
+}
+
+template <typename T>
+std::shared_ptr<T> CADMeshTemplate<T>::FromSTL(G4String file_name) {
+ return std::make_shared<T>(file_name, File::STL);
+}
+
+template <typename T>
+std::shared_ptr<T>
+CADMeshTemplate<T>::FromSTL(G4String file_name,
+ std::shared_ptr<File::Reader> reader) {
+ return std::make_shared<T>(file_name, File::STL, reader);
+}
+
+template <typename T>
+std::shared_ptr<T> CADMeshTemplate<T>::FromOBJ(G4String file_name) {
+ return std::make_shared<T>(file_name, File::OBJ);
+}
+
+template <typename T>
+std::shared_ptr<T>
+CADMeshTemplate<T>::FromOBJ(G4String file_name,
+ std::shared_ptr<File::Reader> reader) {
+ return std::make_shared<T>(file_name, File::OBJ, reader);
+}
+
+template <typename T> CADMeshTemplate<T>::~CADMeshTemplate() {}
+
+template <typename T> bool CADMeshTemplate<T>::IsValidForNavigation() {
+ return reader_->GetMesh()->IsValidForNavigation();
+}
+
+template <typename T> G4String CADMeshTemplate<T>::GetFileName() {
+ return file_name_;
+}
+
+template <typename T> File::Type CADMeshTemplate<T>::GetFileType() {
+ return file_type_;
+}
+
+template <typename T> void CADMeshTemplate<T>::SetVerbose(G4int verbose) {
+ verbose_ = verbose;
+}
+
+template <typename T> G4int CADMeshTemplate<T>::GetVerbose() {
+ return verbose_;
+}
+
+template <typename T> void CADMeshTemplate<T>::SetScale(G4double scale) {
+ scale_ = scale;
+}
+
+template <typename T> G4double CADMeshTemplate<T>::GetScale() { return scale_; }
+
+template <typename T>
+void CADMeshTemplate<T>::SetOffset(G4double x, G4double y, G4double z) {
+ SetOffset(G4ThreeVector(x, y, z));
+}
+
+template <typename T> void CADMeshTemplate<T>::SetOffset(G4ThreeVector offset) {
+ offset_ = offset;
+}
+
+template <typename T> G4ThreeVector CADMeshTemplate<T>::GetOffset() {
+ return offset_;
+}
+}
+
+namespace CADMesh {
+
+namespace Exceptions {
+
+inline void FileNotFound(G4String origin, G4String filepath) {
+ G4Exception(
+ ("CADMesh in " + origin).c_str(), "FileNotFound", FatalException,
+ ("\nThe file: \n\t" + filepath + "\ncould not be found.").c_str());
+}
+
+inline void LexerError(G4String origin, G4String message) {
+ G4Exception(
+ ("CADMesh in " + origin).c_str(), "LexerError", FatalException,
+ ("\nThe CAD file appears to contain incorrect syntax:\n\t" + message)
+ .c_str());
+}
+
+inline void ParserError(G4String origin, G4String message) {
+ G4Exception(("CADMesh in " + origin).c_str(), "ParserError", FatalException,
+ ("\nThe CAD file appears to contain invalid data:\n\t" + message)
+ .c_str());
+}
+
+inline void ReaderCantReadError(G4String origin, File::Type file_type,
+ G4String filepath) {
+ G4Exception(
+ ("CADMesh in " + origin).c_str(), "ReaderCantReadError", FatalException,
+ (G4String("\nThe the reader can't read files of type '") +
+ File::TypeString[file_type] +
+ ".'\n\tSpecified the incorrect file type (?) for file:\n\t\t" + filepath)
+ .c_str());
+}
+
+inline void MeshNotFound(G4String origin, size_t index) {
+ std::stringstream message;
+ message << "\nThe mesh with index '" << index << "' could not be found.";
+
+ G4Exception(("CADMesh in " + origin).c_str(), "MeshNotFound", FatalException,
+ message.str().c_str());
+}
+
+inline void MeshNotFound(G4String origin, G4String name) {
+ G4Exception(
+ ("CADMesh in " + origin).c_str(), "MeshNotFound", FatalException,
+ ("\nThe mesh with name '" + name + "' could not be found.").c_str());
+}
+}
+}
+
+#include "Randomize.hh"
+
+namespace CADMesh {
+
+inline G4VSolid *TessellatedMesh::GetSolid() {
+ return (G4VSolid *)GetTessellatedSolid();
+}
+
+inline G4VSolid *TessellatedMesh::GetSolid(G4int index) {
+ return (G4VSolid *)GetTessellatedSolid(index);
+}
+
+inline G4VSolid *TessellatedMesh::GetSolid(G4String name, G4bool exact) {
+ return (G4VSolid *)GetTessellatedSolid(name, exact);
+}
+
+inline std::vector<G4VSolid *> TessellatedMesh::GetSolids() {
+ std::vector<G4VSolid *> solids;
+
+ for (auto mesh : reader_->GetMeshes()) {
+ solids.push_back(GetTessellatedSolid(mesh));
+ }
+
+ return solids;
+}
+
+inline G4AssemblyVolume *TessellatedMesh::GetAssembly() {
+ if (assembly_) {
+ return assembly_;
+ }
+
+ for (auto mesh : reader_->GetMeshes()) {
+ auto solid = GetTessellatedSolid(mesh);
+
+ G4Material *material = nullptr;
+
+ auto logical =
+ new G4LogicalVolume(solid, material, mesh->GetName() + "_logical");
+
+ G4ThreeVector position = G4ThreeVector();
+ G4RotationMatrix *rotation = new G4RotationMatrix();
+
+ assembly_->AddPlacedVolume(logical, position, rotation);
+ }
+
+ return assembly_;
+}
+
+inline G4TessellatedSolid *TessellatedMesh::GetTessellatedSolid() {
+ return GetTessellatedSolid(0);
+}
+
+inline G4TessellatedSolid *TessellatedMesh::GetTessellatedSolid(G4int index) {
+ return GetTessellatedSolid(reader_->GetMesh(index));
+}
+
+inline G4TessellatedSolid *TessellatedMesh::GetTessellatedSolid(G4String name,
+ G4bool exact) {
+ return GetTessellatedSolid(reader_->GetMesh(name, exact));
+}
+
+inline G4TessellatedSolid *
+TessellatedMesh::GetTessellatedSolid(std::shared_ptr<Mesh> mesh) {
+ auto volume_solid = new G4TessellatedSolid(mesh->GetName());
+
+ for (auto triangle : mesh->GetTriangles()) {
+ auto a = triangle->GetVertex(0) * scale_ + offset_;
+ auto b = triangle->GetVertex(1) * scale_ + offset_;
+ auto c = triangle->GetVertex(2) * scale_ + offset_;
+
+ auto t = new G4TriangularFacet(a, b, c, ABSOLUTE);
+
+ if (reverse_) {
+ volume_solid->AddFacet((G4VFacet *)t->GetFlippedFacet());
+ }
+
+ else {
+ volume_solid->AddFacet((G4VFacet *)t);
+ }
+ }
+
+ volume_solid->SetSolidClosed(true);
+ if (volume_solid->GetNumberOfFacets() == 0) {
+ G4Exception("TessellatedMesh::GetTessellatedSolid",
+ "The loaded mesh has 0 faces.", FatalException,
+ "The file may be empty.");
+ return nullptr;
+ }
+ return volume_solid;
+}
+}
+
+#ifdef USE_CADMESH_TETGEN
+
+namespace CADMesh {
+
+inline TetrahedralMesh::TetrahedralMesh() {}
+
+inline TetrahedralMesh::~TetrahedralMesh() {}
+
+inline G4VSolid *TetrahedralMesh::GetSolid() { return GetSolid(0); }
+
+inline G4VSolid *TetrahedralMesh::GetSolid(G4int /*index*/) { return nullptr; }
+
+inline G4VSolid *TetrahedralMesh::GetSolid(G4String /*name*/,
+ G4bool /*exact*/) {
+
+ return nullptr;
+}
+
+inline std::vector<G4VSolid *> TetrahedralMesh::GetSolids() {
+
+ return std::vector<G4VSolid *>();
+}
+
+inline G4AssemblyVolume *TetrahedralMesh::GetAssembly() {
+ if (assembly_) {
+ return assembly_;
+ }
+
+ assembly_ = new G4AssemblyVolume();
+
+ in_ = std::make_shared<tetgenio>();
+ out_ = std::make_shared<tetgenio>();
+
+ char *fn = (char *)file_name_.c_str();
+
+ G4bool do_tet = true;
+
+ if (file_type_ == File::STL) {
+ in_->load_stl(fn);
+ }
+
+ else if (file_type_ == File::PLY) {
+ in_->load_ply(fn);
+ }
+
+ else if (file_type_ == File::TET) {
+ out_->load_tetmesh(fn, 0);
+ do_tet = false;
+ }
+
+ else if (file_type_ == File::OFF) {
+ out_->load_off(fn);
+ do_tet = false;
+ }
+
+ if (do_tet) {
+ tetgenbehavior *behavior = new tetgenbehavior();
+ behavior->nobisect = 1;
+ behavior->plc = 1;
+ behavior->quality = quality_;
+
+ tetrahedralize(behavior, in_.get(), out_.get());
+ }
+
+ G4RotationMatrix *element_rotation = new G4RotationMatrix();
+ G4ThreeVector element_position = G4ThreeVector();
+ G4Transform3D assembly_transform = G4Translate3D();
+
+ for (int i = 0; i < out_->numberoftetrahedra; i++) {
+ int index_offset = i * 4;
+
+ G4ThreeVector p1 = GetTetPoint(index_offset);
+ G4ThreeVector p2 = GetTetPoint(index_offset + 1);
+ G4ThreeVector p3 = GetTetPoint(index_offset + 2);
+ G4ThreeVector p4 = GetTetPoint(index_offset + 3);
+
+ G4String tet_name =
+ file_name_ + G4String("_tet_") + G4UIcommand::ConvertToString(i);
+
+ auto tet_solid =
+ new G4Tet(tet_name + G4String("_solid"), p1, p2, p3, p4, 0);
+
+ auto tet_logical = new G4LogicalVolume(
+ tet_solid, material_, tet_name + G4String("_logical"), 0, 0, 0);
+
+ assembly_->AddPlacedVolume(tet_logical, element_position, element_rotation);
+ }
+
+ return assembly_;
+}
+
+inline G4ThreeVector TetrahedralMesh::GetTetPoint(G4int index_offset) {
+ return G4ThreeVector(
+ out_->pointlist[out_->tetrahedronlist[index_offset] * 3] * scale_ -
+ offset_.x(),
+ out_->pointlist[out_->tetrahedronlist[index_offset] * 3 + 1] * scale_ -
+ offset_.y(),
+ out_->pointlist[out_->tetrahedronlist[index_offset] * 3 + 2] * scale_ -
+ offset_.z());
+}
+}
+#endif
+
+#define CADMeshLexerToken(name) \
+ static Token name##Token { #name }
+
+#define CADMeshLexerStateDefinition(name) \
+ struct name##State : public State { \
+ State *operator()(Lexer *lexer) const; \
+ }
+
+#define CADMeshLexerState(name) name##State::operator()(Lexer *lexer) const
+
+#define ThisIsA(name) lexer->ThisIsA(name##Token)
+#define StartOfA(name) lexer->StartOfA(name##Token)
+#define EndOfA(name) lexer->EndOfA(name##Token)
+#define MaybeEndOfA(name) lexer->MaybeEndOfA(name##Token)
+
+#define Skip() lexer->Skip()
+
+#define Next() lexer->Peek()
+#define PrintNext() std::cout << lexer->Peek() << std::endl;
+
+#define OneOf(possibles) lexer->OneOf(possibles)
+#define NotOneOf(possibles) !OneOf(possibles)
+
+#define ManyOf(possibles) lexer->ManyOf(possibles)
+#define NotManyOf(possibles) !ManyOf(possibles)
+
+#define Until(match) lexer->Until(match)
+
+#define RestOfLine() \
+ if (Until("\n\r")) \
+ lexer->Backup()
+
+#define MatchExactly(match) lexer->MatchExactly(match)
+#define DoesNotMatchExactly(match) !MatchExactly(match)
+
+#define OneDigit() lexer->OneDigit()
+#define NotOneDigit() !OneDigit()
+
+#define ManyDigits() lexer->ManyDigits()
+#define NotManyDigits() !ManyDigits()
+
+#define OneLetter() lexer->OneLetter()
+#define NotOneLetter() !OneLetter()
+
+#define ManyLetters() lexer->ManyLetters()
+#define NotManyLetters() !ManyLetters()
+
+#define ManyCharacters() lexer->ManyCharacters()
+#define NotManyCharacters() !ManyCharacters()
+
+#define Integer() lexer->Integer()
+#define NotAnInteger() !Integer()
+
+#define Float() lexer->Float()
+#define NotAFloat() !Float()
+
+#define Number() lexer->Number()
+#define NotANumber() !Number()
+
+#define SkipWhiteSpace() lexer->SkipWhiteSpace()
+#define DidNotSkipWhiteSpace() !SkipWhiteSpace()
+
+#define SkipLineBreak() lexer->SkipLineBreak()
+#define DidNotSkipLineBreak() !SkipLineBreak()
+
+#define SkipLineBreaks() lexer->SkipLineBreaks()
+#define DidNotSkipLineBreaks() !SkipLineBreaks()
+
+#define SkipLine() lexer->SkipLine()
+#define DidNotSkipLine() !SkipLine()
+
+#define AtEndOfLine() Next() == "\n" || Next() == "\r"
+
+#define Error(message) \
+ { \
+ lexer->Error(message); \
+ return nullptr; \
+ }
+
+#define NextState(next) return new next##State()
+#define TestState(next) lexer->TestState(new next##State())
+#define TryState(next) \
+ if (TestState(next)) \
+ NextState(next)
+#define FinalState() return new __FinalState();
+
+#define RunLexer(filepath, start) Lexer(filepath, new start##State).GetItems()
+
+namespace CADMesh {
+
+namespace File {
+
+class STLReader : public Reader {
+public:
+ STLReader() : Reader("STLReader"){};
+
+ G4bool Read(G4String filepath);
+ G4bool CanRead(Type file_type);
+
+protected:
+ CADMeshLexerStateDefinition(StartSolid);
+ CADMeshLexerStateDefinition(EndSolid);
+
+ CADMeshLexerStateDefinition(StartFacet);
+ CADMeshLexerStateDefinition(EndFacet);
+
+ CADMeshLexerStateDefinition(StartVertices);
+ CADMeshLexerStateDefinition(EndVertices);
+
+ CADMeshLexerStateDefinition(Vertex);
+
+ CADMeshLexerStateDefinition(ThreeVector);
+
+ std::shared_ptr<Mesh> ParseMesh(Items items);
+ G4TriangularFacet *ParseFacet(Items items);
+ G4TriangularFacet *ParseVertices(Items items);
+ G4ThreeVector ParseThreeVector(Items items);
+};
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+class OBJReader : public Reader {
+public:
+ OBJReader() : Reader("OBJReader"){};
+
+ G4bool Read(G4String filepath);
+ G4bool CanRead(Type file_type);
+
+protected:
+ CADMeshLexerStateDefinition(StartSolid);
+ CADMeshLexerStateDefinition(EndSolid);
+
+ CADMeshLexerStateDefinition(Ignore);
+ CADMeshLexerStateDefinition(Vertex);
+ CADMeshLexerStateDefinition(Facet);
+ CADMeshLexerStateDefinition(Object);
+
+ std::shared_ptr<Mesh> ParseMesh(Items items);
+ G4ThreeVector ParseVertex(Items items);
+ G4TriangularFacet *ParseFacet(Items items, G4bool quad);
+
+private:
+ Points vertices_;
+};
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+CADMeshLexerToken(Header);
+CADMeshLexerToken(Element);
+CADMeshLexerToken(Property);
+
+class PLYReader : public Reader {
+public:
+ PLYReader() : Reader("PLYReader"){};
+
+ G4bool Read(G4String filepath);
+ G4bool CanRead(Type file_type);
+
+protected:
+ CADMeshLexerStateDefinition(StartHeader);
+ CADMeshLexerStateDefinition(EndHeader);
+
+ CADMeshLexerStateDefinition(Element);
+ CADMeshLexerStateDefinition(Property);
+ CADMeshLexerStateDefinition(Ignore);
+
+ CADMeshLexerStateDefinition(Vertex);
+ CADMeshLexerStateDefinition(Facet);
+
+ void ParseHeader(Items items);
+
+ std::shared_ptr<Mesh> ParseMesh(Items vertex_items, Items face_items);
+ G4ThreeVector ParseVertex(Items items);
+ G4TriangularFacet *ParseFacet(Items items, Points vertices);
+
+ size_t vertex_count_ = 0;
+ size_t facet_count_ = 0;
+
+ size_t x_index_ = 0;
+ size_t y_index_ = 0;
+ size_t z_index_ = 0;
+
+ size_t facet_index_ = 0;
+};
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+inline State *STLReader::CADMeshLexerState(StartSolid) {
+ if (DoesNotMatchExactly("solid"))
+ Error("STL files start with 'solid'. Make sure you are using an ASCII STL "
+ "file.");
+
+ SkipWhiteSpace();
+
+ RestOfLine();
+
+ StartOfA(Solid);
+
+ SkipLineBreak();
+ NextState(StartFacet);
+}
+
+inline State *STLReader::CADMeshLexerState(EndSolid) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("endsolid"))
+ Error("STL files end with 'endsolid'.");
+
+ Skip();
+ EndOfA(Solid);
+
+ FinalState();
+}
+
+inline State *STLReader::CADMeshLexerState(StartFacet) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("facet normal"))
+ Error("Facets are indicated by the tag 'facet normal'.");
+
+ SkipWhiteSpace();
+
+ StartOfA(Facet);
+
+ SkipLine();
+
+ NextState(StartVertices);
+}
+
+inline State *STLReader::CADMeshLexerState(EndFacet) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("endfacet"))
+ Error("The end of a facets is indicated by the tag 'endfacet'.");
+
+ SkipWhiteSpace();
+ SkipLineBreak();
+
+ EndOfA(Facet);
+
+ TryState(StartFacet);
+
+ NextState(EndSolid);
+}
+
+inline State *STLReader::CADMeshLexerState(StartVertices) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("outer loop"))
+ Error("The start of the vertices is indicated by the tag 'outer loop'.");
+
+ SkipWhiteSpace();
+ SkipLineBreak();
+
+ StartOfA(Vertices);
+
+ NextState(Vertex);
+}
+
+inline State *STLReader::CADMeshLexerState(EndVertices) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("endloop"))
+ Error("The end of the vertices is indicated by the tag 'endloop'.");
+
+ SkipWhiteSpace();
+ SkipLineBreak();
+
+ EndOfA(Vertices);
+
+ NextState(EndFacet);
+}
+
+inline State *STLReader::CADMeshLexerState(Vertex) {
+ SkipWhiteSpace();
+ SkipLineBreaks();
+ SkipWhiteSpace();
+
+ if (DoesNotMatchExactly("vertex"))
+ Error("A vertex is indicated by the tag 'vertex'.");
+
+ SkipWhiteSpace();
+
+ NextState(ThreeVector);
+}
+
+inline State *STLReader::CADMeshLexerState(ThreeVector) {
+ SkipWhiteSpace();
+
+ StartOfA(ThreeVector);
+
+ if (NotANumber())
+ Error("First number in three vector not found.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Second number in three vector not found.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Third number in three vector not found.");
+
+ ThisIsA(Number);
+ EndOfA(ThreeVector);
+
+ SkipWhiteSpace();
+
+ if (DidNotSkipLineBreak())
+ Error("Expecting a new line at the end of a three vector.");
+
+ TryState(StartVertices);
+
+ TryState(Vertex);
+
+ NextState(EndVertices);
+}
+
+inline G4bool STLReader::Read(G4String filepath) {
+ auto items = RunLexer(filepath, StartSolid);
+
+ if (items.size() == 0) {
+ Exceptions::ParserError("STLReader::Read",
+ "The STL file appears to be empty.");
+ }
+
+ for (auto item : items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The mesh appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("STLReader::Read", error.str());
+ }
+
+ auto mesh = ParseMesh(item.children);
+
+ auto name = G4String(item.value);
+ AddMesh(Mesh::New(mesh, name));
+ }
+
+ return true;
+}
+
+inline G4bool STLReader::CanRead(Type file_type) { return (file_type == STL); }
+
+inline std::shared_ptr<Mesh> STLReader::ParseMesh(Items items) {
+ Triangles triangles;
+
+ for (auto item : items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The facet appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("STLReader::Mesh", error.str());
+ }
+
+ triangles.push_back(ParseFacet(item.children));
+ }
+
+ return Mesh::New(triangles);
+}
+
+inline G4TriangularFacet *STLReader::ParseFacet(Items items) {
+ Triangles triangles;
+
+ for (auto item : items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The vertex appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("STLReader::ParseFacet", error.str());
+ }
+
+ triangles.push_back(ParseVertices(item.children));
+ }
+
+ if (triangles.size() != 1) {
+ std::stringstream error;
+ error << "STL files expect exactly 1 triangle per facet.";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("STLReader::ParseFacet", error.str());
+ }
+
+ return triangles[0];
+}
+
+inline G4TriangularFacet *STLReader::ParseVertices(Items items) {
+ std::vector<G4ThreeVector> vertices;
+
+ for (auto item : items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The vertex appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("STLReader::ParseVertices", error.str());
+ }
+
+ vertices.push_back(ParseThreeVector(item.children));
+ }
+
+ if (vertices.size() != 3) {
+ std::stringstream error;
+ error << "STL files expect exactly 3 vertices for a triangular facet. ";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("STLReader::ParseVertices", error.str());
+ }
+
+ return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
+}
+
+inline G4ThreeVector STLReader::ParseThreeVector(Items items) {
+ std::vector<double> numbers;
+
+ for (auto item : items) {
+ numbers.push_back((double)atof(item.value.c_str()));
+ }
+
+ if (numbers.size() != 3) {
+ std::stringstream error;
+ error << "Three vectors in STL files require exactly 3 numbers";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("STLReader::ParseThreeVector", error.str());
+ }
+
+ return G4ThreeVector(numbers[0], numbers[1], numbers[2]);
+}
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+inline State *OBJReader::CADMeshLexerState(StartSolid) {
+ StartOfA(Solid);
+
+ TryState(Object);
+ TryState(Vertex);
+ TryState(Ignore);
+
+ Error("Invalid element tag.");
+}
+
+inline State *OBJReader::CADMeshLexerState(EndSolid) {
+ if (Next() != "")
+ lexer->LastError();
+
+ EndOfA(Solid);
+ FinalState();
+}
+
+inline State *OBJReader::CADMeshLexerState(Ignore) {
+ if (DidNotSkipLine())
+ NextState(EndSolid);
+
+ TryState(Object);
+ TryState(Vertex);
+ TryState(Facet);
+ TryState(Ignore);
+
+ NextState(EndSolid);
+}
+
+inline State *OBJReader::CADMeshLexerState(Vertex) {
+ SkipLineBreaks();
+
+ if (DoesNotMatchExactly("v "))
+ Error("A vertex is indicated by the tag 'v'.");
+
+ SkipWhiteSpace();
+ StartOfA(Vertex);
+
+ if (NotANumber())
+ Error("First number in three vector not found.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Second number in three vector not found.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Third number in three vector not found.");
+
+ ThisIsA(Number);
+
+ EndOfA(Vertex);
+
+ SkipLine();
+
+ TryState(Vertex);
+ TryState(Object);
+ TryState(Facet);
+ TryState(Ignore);
+
+ NextState(EndSolid);
+}
+
+inline State *OBJReader::CADMeshLexerState(Facet) {
+ SkipLineBreaks();
+
+ if (DoesNotMatchExactly("f "))
+ Error("A facet is indicated by the tag 'f'.");
+
+ SkipWhiteSpace();
+ StartOfA(Facet);
+
+ if (NotANumber())
+ Error("First number in three vector not found.");
+
+ ThisIsA(Number);
+
+ OneOf("/");
+ Number();
+ OneOf("/");
+ Number();
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Second number in three vector not found.");
+
+ ThisIsA(Number);
+
+ OneOf("/");
+ Number();
+ OneOf("/");
+ Number();
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Third number in three vector not found.");
+
+ ThisIsA(Number);
+
+ OneOf("/");
+ Number();
+ OneOf("/");
+ Number();
+ SkipWhiteSpace();
+
+ if (Number())
+ ThisIsA(Number);
+
+ EndOfA(Facet);
+
+ SkipLine();
+
+ TryState(Facet);
+ TryState(Vertex);
+ TryState(Object);
+ TryState(Ignore);
+
+ NextState(EndSolid);
+}
+
+inline State *OBJReader::CADMeshLexerState(Object) {
+ SkipLineBreaks();
+
+ if (DoesNotMatchExactly("o "))
+ Error("An object is indicated by the tag 'o'.");
+
+ EndOfA(Solid);
+
+ SkipWhiteSpace();
+ StartOfA(Solid);
+
+ ManyCharacters();
+ ThisIsA(Word);
+ SkipWhiteSpace();
+
+ TryState(Vertex);
+ TryState(Facet);
+ TryState(Object);
+ TryState(Ignore);
+
+ NextState(EndSolid);
+}
+
+inline G4bool OBJReader::Read(G4String filepath) {
+ auto items = RunLexer(filepath, StartSolid);
+
+ if (items.size() == 0) {
+ Exceptions::ParserError("OBJReader::Read",
+ "The OBJ file appears to be empty.");
+ }
+
+ for (auto item : items) {
+ if (item.children.size() == 0) {
+ continue;
+ }
+
+ auto mesh = ParseMesh(item.children);
+
+ if (mesh->GetTriangles().size() == 0) {
+ continue;
+ }
+
+ G4String name;
+
+ if (item.children[0].token == WordToken) {
+ name = item.children[0].value;
+ }
+
+ AddMesh(Mesh::New(mesh, name));
+ }
+
+ return true;
+}
+
+inline G4bool OBJReader::CanRead(Type file_type) { return (file_type == OBJ); }
+
+inline std::shared_ptr<Mesh> OBJReader::ParseMesh(Items items) {
+ Triangles facets;
+
+ for (auto item : items) {
+ if (item.token != VertexToken) {
+ continue;
+ }
+
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The vertex appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("OBJReader::Mesh", error.str());
+ }
+
+ vertices_.push_back(ParseVertex(item.children));
+ }
+
+ for (auto item : items) {
+ if (item.token != FacetToken) {
+ continue;
+ }
+
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The facet appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("OBJReader::Mesh", error.str());
+ }
+
+ facets.push_back(ParseFacet(item.children, false));
+
+ if (item.children.size() == 4) {
+ facets.push_back(ParseFacet(item.children, true));
+ }
+ }
+
+ return Mesh::New(facets);
+}
+
+inline G4ThreeVector OBJReader::ParseVertex(Items items) {
+ std::vector<double> numbers;
+
+ for (auto item : items) {
+ numbers.push_back((double)atof(item.value.c_str()));
+ }
+
+ if (numbers.size() != 3) {
+ std::stringstream error;
+ error << "Three vectors in OBJ files require exactly 3 numbers";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("OBJReader::ParseThreeVector", error.str());
+ }
+
+ return G4ThreeVector(numbers[0], numbers[1], numbers[2]);
+}
+
+inline G4TriangularFacet *OBJReader::ParseFacet(Items items, G4bool quad) {
+ std::vector<int> indices;
+
+ for (auto item : items) {
+ indices.push_back((int)atoi(item.value.c_str()));
+ }
+
+ if (indices.size() < 3) {
+ std::stringstream error;
+ error << "Facets in OBJ files require at least 3 indicies";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("OBJReader::ParseFacet", error.str());
+ }
+
+ if (quad && indices.size() != 4) {
+ std::stringstream error;
+ error << "Trying to triangle-ify a facet that isn't a quad";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("OBJReader::ParseFacet", error.str());
+ }
+
+ if (quad) {
+ return new G4TriangularFacet(vertices_[indices[0] - 1],
+ vertices_[indices[2] - 1],
+ vertices_[indices[3] - 1], ABSOLUTE);
+ }
+
+ else {
+ return new G4TriangularFacet(vertices_[indices[0] - 1],
+ vertices_[indices[1] - 1],
+ vertices_[indices[2] - 1], ABSOLUTE);
+ }
+}
+}
+}
+
+namespace CADMesh {
+
+namespace File {
+
+inline State *PLYReader::CADMeshLexerState(StartHeader) {
+ if (DoesNotMatchExactly("ply"))
+ Error("PLY files start with 'ply'.");
+
+ StartOfA(Header);
+
+ SkipLine();
+
+ TryState(Element);
+ TryState(Ignore);
+
+ Error("Invalid header tag.");
+}
+
+inline State *PLYReader::CADMeshLexerState(EndHeader) {
+ if (DoesNotMatchExactly("end_header"))
+ Error("PLY file headers end with 'end_header'.");
+
+ MaybeEndOfA(Element);
+
+ EndOfA(Header);
+ FinalState();
+}
+
+inline State *PLYReader::CADMeshLexerState(Element) {
+ if (DoesNotMatchExactly("element "))
+ Error("An element is indicated by the tag 'element'.");
+
+ MaybeEndOfA(Element);
+
+ SkipWhiteSpace();
+ StartOfA(Element);
+
+ if (NotManyCharacters())
+ Error("Element type not found.");
+
+ ThisIsA(Word);
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Element count not found.");
+
+ ThisIsA(Number);
+ SkipLine();
+
+ TryState(Property);
+ TryState(Ignore);
+
+ NextState(EndHeader);
+}
+
+inline State *PLYReader::CADMeshLexerState(Property) {
+ if (DoesNotMatchExactly("property "))
+ Error("An property is indicated by the tag 'property'.");
+
+ SkipWhiteSpace();
+ StartOfA(Property);
+
+ if (NotManyCharacters())
+ Error("Element type not found.");
+
+ ThisIsA(Word);
+ SkipWhiteSpace();
+
+ RestOfLine();
+ ThisIsA(Word);
+
+ EndOfA(Property);
+
+ SkipLineBreak();
+
+ TryState(Property);
+ TryState(Element);
+ TryState(EndHeader);
+ TryState(Ignore);
+
+ NextState(EndHeader);
+}
+
+inline State *PLYReader::CADMeshLexerState(Ignore) {
+ if (DidNotSkipLine())
+ NextState(EndHeader);
+
+ TryState(Element);
+ TryState(Property);
+ TryState(EndHeader);
+
+ NextState(Ignore);
+}
+
+inline State *PLYReader::CADMeshLexerState(Vertex) {
+ SkipLineBreaks();
+ SkipWhiteSpace();
+ SkipLineBreaks();
+
+ StartOfA(Vertex);
+
+ size_t i = 0;
+
+ while (i < 32) {
+ if (AtEndOfLine())
+ break;
+
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Expecting only numbers in the vertex specification.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ i++;
+ }
+
+ EndOfA(Vertex);
+
+ SkipLine();
+
+ TryState(Vertex);
+
+ FinalState();
+}
+
+inline State *PLYReader::CADMeshLexerState(Facet) {
+ SkipLineBreaks();
+ SkipWhiteSpace();
+ SkipLineBreaks();
+
+ StartOfA(Facet);
+
+ size_t i = 0;
+
+ while (i < 32) {
+ if (AtEndOfLine())
+ break;
+
+ SkipWhiteSpace();
+
+ if (NotANumber())
+ Error("Expecting only numbers in the facet specification.");
+
+ ThisIsA(Number);
+ SkipWhiteSpace();
+
+ i++;
+ }
+
+ EndOfA(Facet);
+
+ SkipLine();
+
+ TryState(Facet);
+
+ FinalState();
+}
+
+inline G4bool PLYReader::Read(G4String filepath) {
+ auto lexer = Lexer(filepath, new StartHeaderState);
+ auto items = lexer.GetItems();
+
+ if (items.size() == 0) {
+ std::stringstream error;
+ error << "The header appears to be empty.";
+
+ Exceptions::ParserError("PLYReader::Read", error.str());
+ }
+
+ ParseHeader(items);
+
+ lexer.Run(new VertexState, vertex_count_);
+ auto vertex_items = lexer.GetItems();
+
+ if (vertex_items.size() == 0) {
+ Exceptions::ParserError("PLYReader::Read",
+ "The PLY file appears to have no vertices.");
+ }
+
+ if (vertex_items.size() != vertex_count_) {
+ Exceptions::ParserError("PLYReader::Read",
+ "The PLY file appears to be missing vertices.");
+ }
+
+ lexer.Run(new FacetState, facet_count_);
+ auto face_items = lexer.GetItems();
+
+ if (face_items.size() == 0) {
+ Exceptions::ParserError("PLYReader::Read",
+ "The PLY file appears to have no facets.");
+ }
+
+ if (face_items.size() != facet_count_) {
+ Exceptions::ParserError("PLYReader::Read",
+ "The PLY file appears to be missing facets");
+ }
+
+ auto mesh = ParseMesh(vertex_items, face_items);
+ AddMesh(Mesh::New(mesh));
+
+ return true;
+}
+
+inline G4bool PLYReader::CanRead(Type file_type) { return (file_type == PLY); }
+
+inline void PLYReader::ParseHeader(Items items) {
+ if (items.size() != 1) {
+ std::stringstream error;
+ error << "The header appears to be invalid or missing."
+ << "Error around line 1.";
+
+ Exceptions::ParserError("PLYReader::ParseHeader", error.str());
+ }
+
+ for (auto item : items[0].children) {
+ if (item.token == ElementToken) {
+ if (item.children.size() < 2) {
+ std::stringstream error;
+ error << "Invalid element information in header. Expecting 'vertex' or "
+ "'face' and a number."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("PLYReader::ParseHeader", error.str());
+ }
+
+ if (item.children[0].token == WordToken &&
+ item.children[1].token == NumberToken) {
+ if (item.children[0].value == "vertex") {
+ vertex_count_ = atoi(item.children[1].value.c_str());
+
+ for (size_t i = 2; i < item.children.size(); i++) {
+ auto property = item.children[i];
+
+ if (property.children.size() > 1) {
+ if (property.children[1].token == WordToken) {
+ if (property.children[1].value == "x") {
+ x_index_ = i - 2;
+ }
+
+ if (property.children[1].value == "y") {
+ y_index_ = i - 2;
+ }
+
+ if (property.children[1].value == "z") {
+ z_index_ = i - 2;
+ }
+ }
+ }
+ }
+ }
+
+ else if (item.children[0].value == "face") {
+ facet_count_ = atoi(item.children[1].value.c_str());
+
+ for (size_t i = 2; i < item.children.size(); i++) {
+ auto property = item.children[i];
+
+ if (property.children.size() > 1) {
+ if (property.children[1].token == WordToken) {
+ if (property.children[1].value == "uchar int vertex_indices") {
+ facet_index_ = i - 2;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ if (vertex_count_ == 0) {
+ std::stringstream error;
+ error << "The number of vertices was not found in the header.";
+
+ Exceptions::ParserError("PLYReader::ParseHeader", error.str());
+ }
+
+ if (facet_count_ == 0) {
+ std::stringstream error;
+ error << "The number of faces was not found in the header.";
+
+ Exceptions::ParserError("PLYReader::ParseHeader", error.str());
+ }
+
+ if ((x_index_ == y_index_) || (y_index_ == z_index_) ||
+ (x_index_ == z_index_)) {
+ std::stringstream error;
+ error << "The vertex x, y, z indices were not found in the header.";
+
+ Exceptions::ParserError("PLYReader::ParseHeader", error.str());
+ }
+}
+
+inline std::shared_ptr<Mesh> PLYReader::ParseMesh(Items vertex_items,
+ Items face_items) {
+ Points vertices;
+ Triangles facets;
+
+ for (auto item : vertex_items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The vertex appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("PLYReader::ParseMesh", error.str());
+ }
+
+ if (item.token == VertexToken) {
+ vertices.push_back(ParseVertex(item.children));
+ }
+ }
+
+ for (auto item : face_items) {
+ if (item.children.size() == 0) {
+ std::stringstream error;
+ error << "The facet appears to be empty."
+ << "Error around line " << item.line << ".";
+
+ Exceptions::ParserError("PLYReader::Mesh", error.str());
+ }
+
+ if (item.token == FacetToken) {
+ facets.push_back(ParseFacet(item.children, vertices));
+ }
+ }
+
+ return Mesh::New(facets);
+}
+
+inline G4ThreeVector PLYReader::ParseVertex(Items items) {
+ std::vector<double> numbers;
+
+ for (auto item : items) {
+ numbers.push_back((double)atof(item.value.c_str()));
+ }
+
+ if (numbers.size() < 3) {
+ std::stringstream error;
+ error << "Vertices in PLY files require atleast 3 numbers. ";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("PLYReader::ParseVertex", error.str());
+ }
+
+ return G4ThreeVector(numbers[x_index_], numbers[y_index_], numbers[z_index_]);
+}
+
+inline G4TriangularFacet *PLYReader::ParseFacet(Items items, Points vertices) {
+ std::vector<int> indices;
+
+ for (auto item : items) {
+ indices.push_back((int)atoi(item.value.c_str()));
+ }
+
+ if (indices.size() < 4) {
+ std::stringstream error;
+ error << "Facets in PLY files require 3 indicies";
+
+ if (items.size() != 0) {
+ error << "Error around line " << items[0].line << ".";
+ }
+
+ Exceptions::ParserError("PLYReader::ParseFacet", error.str());
+ }
+
+ return new G4TriangularFacet(vertices[indices[1 + facet_index_]],
+ vertices[indices[2 + facet_index_]],
+ vertices[indices[3 + facet_index_]], ABSOLUTE);
+}
+}
+}
+
+#ifdef USE_CADMESH_ASSIMP_READER
+
+namespace CADMesh {
+
+namespace File {
+
+inline ASSIMPReader::ASSIMPReader() : Reader("ASSIMPReader") {
+ importer_ = new Assimp::Importer();
+}
+
+inline ASSIMPReader::~ASSIMPReader() { delete importer_; }
+
+inline G4bool ASSIMPReader::Read(G4String filepath) {
+ auto scene = importer_->ReadFile(filepath.c_str(),
+ aiProcess_Triangulate |
+ aiProcess_JoinIdenticalVertices |
+ aiProcess_CalcTangentSpace);
+
+ if (!scene) {
+ Exceptions::FileNotFound("ASSIMP::Read", filepath);
+ return false;
+ }
+
+ for (unsigned int index = 0; index < scene->mNumMeshes; index++) {
+ aiMesh *mesh = scene->mMeshes[index];
+ auto name = mesh->mName.C_Str();
+
+ Triangles triangles;
+
+ for (unsigned int i = 0; i < mesh->mNumFaces; i++) {
+ const aiFace &face = mesh->mFaces[i];
+
+ G4ThreeVector a(mesh->mVertices[face.mIndices[0]].x,
+ mesh->mVertices[face.mIndices[0]].y,
+ mesh->mVertices[face.mIndices[0]].z);
+
+ G4ThreeVector b(mesh->mVertices[face.mIndices[1]].x,
+ mesh->mVertices[face.mIndices[1]].y,
+ mesh->mVertices[face.mIndices[1]].z);
+
+ G4ThreeVector c(mesh->mVertices[face.mIndices[2]].x,
+ mesh->mVertices[face.mIndices[2]].y,
+ mesh->mVertices[face.mIndices[2]].z);
+
+ triangles.push_back(new G4TriangularFacet(a, b, c, ABSOLUTE));
+ }
+
+ AddMesh(Mesh::New(triangles, name));
+ }
+
+ return true;
+}
+
+inline G4bool ASSIMPReader::CanRead(Type /*file_type*/) { return true; }
+
+std::shared_ptr<ASSIMPReader> ASSIMP() {
+ return std::make_shared<ASSIMPReader>();
+}
+}
+}
+#endif
+
+namespace CADMesh {
+
+namespace File {
+
+inline G4bool BuiltInReader::Read(G4String filepath) {
+ File::Reader *reader = nullptr;
+
+ auto type = TypeFromName(filepath);
+
+ if (type == STL) {
+ reader = new File::STLReader();
+ }
+
+ else if (type == OBJ) {
+ reader = new File::OBJReader();
+ }
+
+ else if (type == PLY) {
+ reader = new File::PLYReader();
+ }
+
+ else {
+ Exceptions::ReaderCantReadError("BuildInReader::Read", type, filepath);
+ }
+
+ if (!reader->Read(filepath)) {
+ return false;
+ }
+
+ SetMeshes(reader->GetMeshes());
+ return true;
+}
+
+inline G4bool BuiltInReader::CanRead(Type type) {
+ return type == STL || type == OBJ || type == PLY;
+}
+
+inline std::shared_ptr<BuiltInReader> BuiltIn() {
+ return std::make_shared<BuiltInReader>();
+}
+}
+}
diff --git a/NPSimulation/Detectors/Strasse/Strasse.cc b/NPSimulation/Detectors/Strasse/Strasse.cc
index 47e27e931ffd16034842e0504ec7d246cda232c3..0ee0ab3daf49125fde7313d92332e070399c7085 100644
--- a/NPSimulation/Detectors/Strasse/Strasse.cc
+++ b/NPSimulation/Detectors/Strasse/Strasse.cc
@@ -55,6 +55,9 @@
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
+// CAD Mesh
+#include "CADMesh.hh"
+
using namespace std;
using namespace CLHEP;
@@ -85,6 +88,8 @@ namespace Strasse_NS{
double Inner_PCB_DownstreamWidth=2*mm;
double Inner_PCB_MidWidth=2*mm;
double Inner_PCB_Thickness=3*mm;
+ double Inner_PCB_Ledge = 1*mm ;
+ double Inner_PCB_Step = 2*mm ;
double Inner_Wafer_TransverseStrips= 128;
double Inner_Wafer_LongitudinalStrips= 128;
@@ -108,6 +113,8 @@ namespace Strasse_NS{
double Outer_PCB_DownstreamWidth=2*mm;
double Outer_PCB_MidWidth=2*mm;
double Outer_PCB_Thickness=3*mm;
+ double Outer_PCB_Ledge = 1*mm ;
+ double Outer_PCB_Step = 2*mm ;
double Outer_Wafer_TransverseStrips= 128;
double Outer_Wafer_LongitudinalStrips= 128;
@@ -140,8 +147,14 @@ Strasse::Strasse(){
m_InnerDetector=0;
m_OuterDetector=0;
m_Chamber=0;
- m_Frame=0;
+ m_Blades=0;
+ m_Stars=0;
+ m_Base=0;
m_Electronic=0;
+ found_chamber = false;
+ found_stars = false;
+ found_blades = false;
+ found_base = false;
// Dark Grey
SiliconVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
// Green
@@ -149,11 +162,13 @@ Strasse::Strasse(){
// Gold Yellow
PADVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.2)) ;
// Light Grey
- FrameVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
+ StarsVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
// Space transparent
ChamberVisAtt = new G4VisAttributes(G4Colour(0.3, 0.4, 0.5,0.2)) ;
// Light Blue
- GuardRingVisAtt = new G4VisAttributes(G4Colour(0, 0, 0,0.5)) ;
+ GuardRingVisAtt = new G4VisAttributes(G4Colour(0.85, 0.85, 0.85,0.5)) ;
+ // Light Blue
+ BladeVisAtt = new G4VisAttributes(G4Colour(1, 0.65, 0.0,0.7)) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -199,12 +214,15 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
+Inner_PCB_PortWidth;
vector<G4TwoVector> PCBCrossSection;
- double l1 = Inner_PCB_Thickness*0.5/tan(Inner_PCB_BevelAngle);
+
+ double l1;
+ if(Inner_PCB_BevelAngle==90) l1 = 0;
+ else l1 = Inner_PCB_Thickness*0.5/tan(Inner_PCB_BevelAngle);
PCBCrossSection.push_back(G4TwoVector(Inner_PCB_Width*0.5-l1,-Inner_PCB_Thickness*0.5));
PCBCrossSection.push_back(G4TwoVector(Inner_PCB_Width*0.5,Inner_PCB_Thickness*0.5));
- PCBCrossSection.push_back(G4TwoVector(-Inner_PCB_Width*0.5-l1,Inner_PCB_Thickness*0.5));
- PCBCrossSection.push_back(G4TwoVector(-Inner_PCB_Width*0.5,-Inner_PCB_Thickness*0.5));
+ PCBCrossSection.push_back(G4TwoVector(-Inner_PCB_Width*0.5,Inner_PCB_Thickness*0.5));
+ PCBCrossSection.push_back(G4TwoVector(-Inner_PCB_Width*0.5+l1,-Inner_PCB_Thickness*0.5));
G4ExtrudedSolid* PCBFull =
new G4ExtrudedSolid("PCBFull",
@@ -219,7 +237,7 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
m_InnerDetector->SetVisAttributes(G4VisAttributes::Invisible);
///////////////////////////////////////////////////////////////////////////
- // Build the external PCB layer
+ // Build the external PCB frame
// Calculate the hole shift within the PCB
double Width_Shift= -0.5*Inner_PCB_Width + 0.5*Inner_Wafer_Width // Flush to border
+Inner_PCB_PortWidth; // add the port side shift
@@ -256,33 +274,23 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
false,0);
///////////////////////////////////////////////////////////////////////////
- // Build the internal PCB layer
- double offsetPCB2 = -0.3;
- double Inner_PCB2_Thickness = Inner_PCB_Thickness + offsetPCB2;
+ // Build the PCB Ledge on wich Silicon is glued
+ double Inner_PCB2_Thickness = Inner_PCB_Step; //Step size
+ double offsetPCB2 = Inner_PCB2_Thickness - Inner_PCB_Thickness;
- double Inner_PCB2_StarboardWidth = Inner_PCB_StarboardWidth;
- double Inner_PCB2_PortWidth = Inner_PCB_PortWidth;
- double Inner_PCB2_UpstreamWidth = Inner_PCB_UpstreamWidth;
double Inner_PCB2_MidWidth = Inner_PCB_MidWidth;
- double Inner_PCB2_DownstreamWidth =Inner_PCB_DownstreamWidth;
// perpendicular to beam axis
- double Inner_PCB2_Width= Inner_Wafer_Width-Inner_PCB_PortWidth*2
- +Inner_PCB2_StarboardWidth
- +Inner_PCB2_PortWidth;
+ double Inner_PCB2_Width= Inner_Wafer_Width;
vector<G4TwoVector> PCB2CrossSection;
- double l21 = Inner_PCB2_Thickness*0.5/tan(Inner_PCB_BevelAngle);
-
- PCB2CrossSection.push_back(G4TwoVector(Inner_PCB2_Width*0.5-l21,-Inner_PCB2_Thickness*0.5));
+ PCB2CrossSection.push_back(G4TwoVector(Inner_PCB2_Width*0.5,-Inner_PCB2_Thickness*0.5));
PCB2CrossSection.push_back(G4TwoVector(Inner_PCB2_Width*0.5,Inner_PCB2_Thickness*0.5));
- PCB2CrossSection.push_back(G4TwoVector(-Inner_PCB2_Width*0.5-l21,Inner_PCB2_Thickness*0.5));
+ PCB2CrossSection.push_back(G4TwoVector(-Inner_PCB2_Width*0.5,Inner_PCB2_Thickness*0.5));
PCB2CrossSection.push_back(G4TwoVector(-Inner_PCB2_Width*0.5,-Inner_PCB2_Thickness*0.5));
- double Inner_PCB2_Length= 2*(Inner_Wafer_Length-Inner_PCB_PortWidth)
- +Inner_PCB2_UpstreamWidth
- +Inner_PCB2_MidWidth
- +Inner_PCB2_DownstreamWidth;
+ //double Inner_PCB2_Length= Inner_PCB_Length;
+ double Inner_PCB2_Length= 2*Inner_Wafer_Length + Inner_PCB_MidWidth;
G4ExtrudedSolid* PCB2Full =
new G4ExtrudedSolid("PCB2Full",
@@ -291,36 +299,54 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
G4TwoVector(0,0),1,// offset, scale
G4TwoVector(0,0),1);// offset, scale
- G4ThreeVector HoleShift21 = G4ThreeVector(0, 0, 0);
+
+ double Length_Shift21 = -0.5*Inner_PCB_Length // Flush to border
+ + 0.5*(Inner_PCB_UpstreamWidth+Inner_PCB_DownstreamWidth) // add Upstream side shift
+ + 0.5*Inner_Wafer_Length;
+ double Length_Shift22 = Length_Shift21 // overlap detector 1
+ + Inner_Wafer_Length // at opposing edge
+ + Inner_PCB_MidWidth; // after mid width
+
+ G4ThreeVector HoleShift21 = G4ThreeVector(0, 0, Length_Shift21);
+ G4ThreeVector HoleShift22 = G4ThreeVector(0, 0, Length_Shift22);
+
G4Box* HoleShape2 = new G4Box("HoleShape2",
- Inner_PCB2_Width*0.5 -1,
+ Inner_Wafer_Width*0.5 - Inner_PCB_Ledge,
Inner_PCB2_Thickness,
- Inner_PCB2_Length*0.5-1);
+ Inner_Wafer_Length*0.5 - Inner_PCB_Ledge);
// Substracting the hole Shape from the Stock PCB
G4SubtractionSolid* PCB2_1 = new G4SubtractionSolid("PCB2_1", PCB2Full, HoleShape2,
new G4RotationMatrix,HoleShift21);
-
+ G4SubtractionSolid* PCB2_2 = new G4SubtractionSolid("PCB2_2", PCB2_1, HoleShape2,
+ new G4RotationMatrix,HoleShift22);
+
G4ThreeVector HoleCenterBar = G4ThreeVector(0, 0, 0);
G4Box* HoleShapeCenterBar = new G4Box("HoleShapeCenterBar",
- Inner_PCB2_Width*0.5,
+ Inner_PCB2_Width*0.5+0.1,
Inner_PCB2_Thickness,
- 0.5*mm);
+ Inner_PCB2_MidWidth*0.5);
- // Substracting the hole Shape from the Stock PCB
- G4SubtractionSolid* PCB2_2 = new G4SubtractionSolid("PCB2_2", PCB2_1, HoleShapeCenterBar,
+ G4SubtractionSolid* PCB2_3 = new G4SubtractionSolid("PCB2_3", PCB2_2, HoleShapeCenterBar,
new G4RotationMatrix,HoleCenterBar);
// Sub Volume PCB
G4LogicalVolume* logicPCB2 =
- new G4LogicalVolume(PCB2_2,m_MaterialPCB,"logicPCB2", 0, 0, 0);
- logicPCB2->SetVisAttributes(PCBVisAtt);
+ new G4LogicalVolume(PCB2_3,m_MaterialPCB,"logicPCB2", 0, 0, 0);
+ logicPCB2->SetVisAttributes(PADVisAtt);
+
+ // Offset along beam axis between PCB middle and (2*Wafer+MiddleBar) volume center
+ double CentralZOffset = - Inner_PCB_Length*0.5
+ + Inner_PCB_UpstreamWidth
+ + Inner_Wafer_Length
+ + Inner_PCB_MidWidth*0.5;
new G4PVPlacement(new G4RotationMatrix(0,0,0),
- G4ThreeVector(0,-0.5*offsetPCB2,0),
+ G4ThreeVector(0,-0.5*offsetPCB2,CentralZOffset),
logicPCB2,"Strasse_Inner_PCB2",m_InnerDetector,
false,0);
+
///////////////////////////////////////////////////////////////////////////
// Build the Wafer
// Sub volume Wafer
@@ -337,19 +363,23 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
new G4LogicalVolume(WaferShape,m_MaterialSilicon,"logicWafer2", 0, 0, 0);
logicWafer2->SetVisAttributes(GuardRingVisAtt);
+ // Shift along Y to flush the wafer to the pcb ledge on one side
+ G4ThreeVector WaferShiftY = G4ThreeVector(0,-0.5*Inner_Wafer_Thickness
+ -Inner_Wafer_AlThickness
+ -0.5*Inner_PCB_Thickness
+ -offsetPCB2-0.05,0);
+
+ G4ThreeVector WaferShiftZ = G4ThreeVector(0,0,-Inner_PCB_UpstreamWidth-Inner_PCB_MidWidth);
+
new G4PVPlacement(new G4RotationMatrix(0,0,0),
- G4ThreeVector(0,0.5*Inner_Wafer_Thickness
- +Inner_Wafer_AlThickness
- -0.5*Inner_PCB_Thickness,0) // flush the wafer to the pcb on one side
- +HoleShift1, // Shift wafer in the hole
+ WaferShiftY+WaferShiftZ // Shift along Y
+ +HoleShift21, // Shift along Z to putwafer in the 1st hole
logicWafer1,"Strasse_Inner_Wafer1",m_InnerDetector,
false,0);
new G4PVPlacement(new G4RotationMatrix(0,0,0),
- G4ThreeVector(0,0.5*Inner_Wafer_Thickness
- +Inner_Wafer_AlThickness
- -0.5*Inner_PCB_Thickness,0)// flush the wafer to the pcb on one side
- +HoleShift2, // Shift wafer in the hole
+ WaferShiftY+WaferShiftZ// Shift along Y
+ +HoleShift22, // Shift along Z to put wafer in the 2nd hole
logicWafer2,"Strasse_Inner_Wafer2",m_InnerDetector,
false,0);
@@ -378,7 +408,9 @@ G4LogicalVolume* Strasse::BuildInnerDetector(){
G4ThreeVector(0,0,-0.5*(Inner_Wafer_PADExternal-Inner_Wafer_PADInternal)), // assymetric pading for bounding
logicActiveWafer2,"Strasse_Inner_ActiveWafer2",logicWafer2,
false,2);
+
}
+
return m_InnerDetector;
}
@@ -399,12 +431,14 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){
vector<G4TwoVector> PCBCrossSection;
- double l1 = Outer_PCB_Thickness*0.5/tan(Outer_PCB_BevelAngle);
+ double l1;
+ if(Outer_PCB_BevelAngle==90) l1 = 0;
+ else l1 = Outer_PCB_Thickness*0.5/tan(Outer_PCB_BevelAngle);
PCBCrossSection.push_back(G4TwoVector(Outer_PCB_Width*0.5-l1,-Outer_PCB_Thickness*0.5));
PCBCrossSection.push_back(G4TwoVector(Outer_PCB_Width*0.5,Outer_PCB_Thickness*0.5));
- PCBCrossSection.push_back(G4TwoVector(-Outer_PCB_Width*0.5-l1,Outer_PCB_Thickness*0.5));
- PCBCrossSection.push_back(G4TwoVector(-Outer_PCB_Width*0.5,-Outer_PCB_Thickness*0.5));
+ PCBCrossSection.push_back(G4TwoVector(-Outer_PCB_Width*0.5,Outer_PCB_Thickness*0.5));
+ PCBCrossSection.push_back(G4TwoVector(-Outer_PCB_Width*0.5+l1,-Outer_PCB_Thickness*0.5));
G4ExtrudedSolid* PCBFull =
new G4ExtrudedSolid("PCBFull",
@@ -456,32 +490,21 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){
///////////////////////////////////////////////////////////////////////////
// Build the internal PCB layer
- double offsetPCB2 = -0.3;
- double Outer_PCB2_Thickness = Outer_PCB_Thickness + offsetPCB2;
+ double Outer_PCB2_Thickness = Outer_PCB_Step;//Step size
+ double offsetPCB2 = Outer_PCB2_Thickness - Outer_PCB_Thickness;
- double Outer_PCB2_StarboardWidth = Outer_PCB_StarboardWidth;
- double Outer_PCB2_PortWidth = Outer_PCB_PortWidth;
- double Outer_PCB2_UpstreamWidth =Outer_PCB_UpstreamWidth;
- double Outer_PCB2_MidWidth =Outer_PCB_MidWidth;
- double Outer_PCB2_DownstreamWidth =Outer_PCB_DownstreamWidth;
+ double Outer_PCB2_MidWidth = Outer_PCB_MidWidth;
// perpendicular to beam axis
- double Outer_PCB2_Width= Outer_Wafer_Width-Outer_PCB_PortWidth*2
- +Outer_PCB2_StarboardWidth
- +Outer_PCB2_PortWidth;
+ double Outer_PCB2_Width= Outer_Wafer_Width;
vector<G4TwoVector> PCB2CrossSection;
- double l21 = Outer_PCB2_Thickness*0.5/tan(Outer_PCB_BevelAngle);
-
- PCB2CrossSection.push_back(G4TwoVector(Outer_PCB2_Width*0.5-l21,-Outer_PCB2_Thickness*0.5));
+ PCB2CrossSection.push_back(G4TwoVector(Outer_PCB2_Width*0.5,-Outer_PCB2_Thickness*0.5));
PCB2CrossSection.push_back(G4TwoVector(Outer_PCB2_Width*0.5,Outer_PCB2_Thickness*0.5));
- PCB2CrossSection.push_back(G4TwoVector(-Outer_PCB2_Width*0.5-l21,Outer_PCB2_Thickness*0.5));
+ PCB2CrossSection.push_back(G4TwoVector(-Outer_PCB2_Width*0.5,Outer_PCB2_Thickness*0.5));
PCB2CrossSection.push_back(G4TwoVector(-Outer_PCB2_Width*0.5,-Outer_PCB2_Thickness*0.5));
- double Outer_PCB2_Length= 2*(Outer_Wafer_Length-Outer_PCB_PortWidth)
- +Outer_PCB2_UpstreamWidth
- +Outer_PCB2_MidWidth
- +Outer_PCB2_DownstreamWidth;
+ double Outer_PCB2_Length= 2*Outer_Wafer_Length + Outer_PCB_MidWidth;
G4ExtrudedSolid* PCB2Full =
new G4ExtrudedSolid("PCB2Full",
@@ -490,33 +513,51 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){
G4TwoVector(0,0),1,// offset, scale
G4TwoVector(0,0),1);// offset, scale
- G4ThreeVector HoleShift21 = G4ThreeVector(0, 0, 0);
+ // Offset along beam axis between PCB middle and (2*Wafer+MiddleBar) volume center
+ double CentralZOffset = - Outer_PCB_Length*0.5
+ + Outer_PCB_UpstreamWidth
+ + Outer_Wafer_Length
+ + Outer_PCB_MidWidth*0.5;
+
+ double Length_Shift21 = -0.5*Outer_PCB_Length // Flush to border
+ + 0.5*(Outer_PCB_UpstreamWidth+Outer_PCB_DownstreamWidth) // add Upstream side shift
+ + 0.5*Outer_Wafer_Length;
+ double Length_Shift22 = Length_Shift21 // overlap detector 1
+ + Outer_Wafer_Length // at opposing edge
+ + Outer_PCB_MidWidth; // after mid width
+
+ G4ThreeVector HoleShift21 = G4ThreeVector(0, 0, Length_Shift21);
+ G4ThreeVector HoleShift22 = G4ThreeVector(0, 0, Length_Shift22);
+
G4Box* HoleShape2 = new G4Box("HoleShape2",
- Outer_PCB2_Width*0.5 -1,
+ Outer_Wafer_Width*0.5 - Outer_PCB_Ledge,
Outer_PCB2_Thickness,
- Outer_PCB2_Length*0.5-1);
+ Outer_Wafer_Length*0.5 - Outer_PCB_Ledge);
// Substracting the hole Shape from the Stock PCB
G4SubtractionSolid* PCB2_1 = new G4SubtractionSolid("PCB2_1", PCB2Full, HoleShape2,
new G4RotationMatrix,HoleShift21);
+ G4SubtractionSolid* PCB2_2 = new G4SubtractionSolid("PCB2_2", PCB2_1, HoleShape2,
+ new G4RotationMatrix,HoleShift22);
G4ThreeVector HoleCenterBar = G4ThreeVector(0, 0, 0);
G4Box* HoleShapeCenterBar = new G4Box("HoleShapeCenterBar",
- Outer_PCB2_Width*0.5,
+ Outer_PCB2_Width*0.5+0.1,
Outer_PCB2_Thickness,
- 0.5*mm);
+ Outer_PCB2_MidWidth*0.5);
// Substracting the hole Shape from the Stock PCB
- G4SubtractionSolid* PCB2_2 = new G4SubtractionSolid("PCB2_2", PCB2_1, HoleShapeCenterBar,
+ G4SubtractionSolid* PCB2_3 = new G4SubtractionSolid("PCB2_3", PCB2_2, HoleShapeCenterBar,
new G4RotationMatrix,HoleCenterBar);
// Sub Volume PCB
G4LogicalVolume* logicPCB2 =
- new G4LogicalVolume(PCB2_2,m_MaterialPCB,"logicPCB2", 0, 0, 0);
- logicPCB2->SetVisAttributes(PCBVisAtt);
+ new G4LogicalVolume(PCB2_3,m_MaterialPCB,"logicPCB2", 0, 0, 0);
+ logicPCB2->SetVisAttributes(PADVisAtt);
new G4PVPlacement(new G4RotationMatrix(0,0,0),
- G4ThreeVector(0,-0.5*offsetPCB2,0),
+ //G4ThreeVector(0,-0.5*offsetPCB2,Outer_PCB_UpstreamWidth-Outer_PCB_MidWidth),
+ G4ThreeVector(0,-0.5*offsetPCB2,CentralZOffset),
logicPCB2,"Strasse_Outer_PCB2",m_OuterDetector,
false,0);
@@ -577,6 +618,7 @@ G4LogicalVolume* Strasse::BuildOuterDetector(){
G4ThreeVector(0,0,-0.5*(Outer_Wafer_PADExternal-Outer_Wafer_PADInternal)), // assymetric pading for bounding
logicActiveWafer2,"Strasse_Outer_ActiveWafer2",logicWafer2,
false,2);
+
}
return m_OuterDetector;
}
@@ -655,6 +697,78 @@ G4LogicalVolume* Strasse::BuildChamber(){
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Strasse::BuildChamberFromCAD(string path){
+ if(!m_Chamber){
+ auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
+ mesh->SetScale(mm);
+ //mesh->SetOffset(offset);
+ mesh->SetReverse(false);
+
+ auto cad_solid = mesh->GetSolid();
+ m_Chamber = new G4LogicalVolume(cad_solid,
+ m_MaterialAl,
+ "Strasse_Chamber", 0, 0, 0);
+
+ m_Chamber->SetVisAttributes(ChamberVisAtt);
+ }
+ return m_Chamber;
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Strasse::BuildBlades(string path){
+ if(!m_Blades){
+ auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
+ mesh->SetScale(mm);
+ //mesh->SetOffset(offset);
+ mesh->SetReverse(false);
+
+ auto cad_solid = mesh->GetSolid();
+ m_Blades = new G4LogicalVolume(cad_solid,
+ m_MaterialCu,
+ "Strasse_Blades", 0, 0, 0);
+
+ m_Blades->SetVisAttributes(BladeVisAtt);
+ }
+ return m_Blades;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Strasse::BuildStars(string path){
+ if(!m_Stars){
+ auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
+ mesh->SetScale(mm);
+ //mesh->SetOffset(offset);
+ mesh->SetReverse(false);
+
+ auto cad_solid = mesh->GetSolid();
+ m_Stars = new G4LogicalVolume(cad_solid,
+ m_MaterialAl,
+ "Strasse_Stars", 0, 0, 0);
+
+ m_Stars->SetVisAttributes(StarsVisAtt);
+ }
+ return m_Stars;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Strasse::BuildBase(string path){
+ if(!m_Base){
+ auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
+ mesh->SetScale(mm);
+ //mesh->SetOffset(offset);
+ mesh->SetReverse(false);
+
+ auto cad_solid = mesh->GetSolid();
+ m_Base = new G4LogicalVolume(cad_solid,
+ m_MaterialAl,
+ "Strasse_Base", 0, 0, 0);
+
+ m_Base->SetVisAttributes(StarsVisAtt);
+ }
+ return m_Base;
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -688,6 +802,8 @@ void Strasse::ReadConfiguration(NPL::InputParser parser){
"Inner_PCB_DownstreamWidth",
"Inner_PCB_MidWidth",
"Inner_PCB_Thickness",
+ "Inner_PCB_Ledge",
+ "Inner_PCB_Step",
"Inner_Wafer_TransverseStrips",
"Inner_Wafer_LongitudinalStrips",
"Outer_Wafer_Length",
@@ -704,6 +820,8 @@ void Strasse::ReadConfiguration(NPL::InputParser parser){
"Outer_PCB_DownstreamWidth",
"Outer_PCB_MidWidth",
"Outer_PCB_Thickness",
+ "Outer_PCB_Ledge",
+ "Outer_PCB_Step",
"Outer_Wafer_TransverseStrips",
"Outer_Wafer_LongitudinalStrips",
"Chamber_Thickness",
@@ -734,6 +852,8 @@ void Strasse::ReadConfiguration(NPL::InputParser parser){
Inner_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_DownstreamWidth","mm");
Inner_PCB_MidWidth = blocks_info[0]->GetDouble("Inner_PCB_MidWidth","mm");
Inner_PCB_Thickness = blocks_info[0]->GetDouble("Inner_PCB_Thickness","mm");
+ Inner_PCB_Ledge = blocks_info[0]->GetDouble("Inner_PCB_Ledge","mm");
+ Inner_PCB_Step = blocks_info[0]->GetDouble("Inner_PCB_Step","mm");
Outer_Wafer_Length = blocks_info[0]->GetDouble("Outer_Wafer_Length","mm");
Outer_Wafer_Width = blocks_info[0]->GetDouble("Outer_Wafer_Width","mm");
Outer_Wafer_Thickness = blocks_info[0]->GetDouble("Outer_Wafer_Thickness","mm");
@@ -750,6 +870,8 @@ void Strasse::ReadConfiguration(NPL::InputParser parser){
Outer_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_DownstreamWidth","mm");
Outer_PCB_MidWidth = blocks_info[0]->GetDouble("Outer_PCB_MidWidth","mm");
Outer_PCB_Thickness = blocks_info[0]->GetDouble("Outer_PCB_Thickness","mm");
+ Outer_PCB_Ledge = blocks_info[0]->GetDouble("Outer_PCB_Ledge","mm");
+ Outer_PCB_Step = blocks_info[0]->GetDouble("Outer_PCB_Step","mm");
Chamber_Thickness= blocks_info[0]->GetDouble("Chamber_Thickness","mm");
Chamber_Cylinder_Length= blocks_info[0]->GetDouble("Chamber_Cylinder_Length","mm");
Chamber_Radius= blocks_info[0]->GetDouble("Chamber_Radius","mm");
@@ -837,6 +959,32 @@ void Strasse::ReadConfiguration(NPL::InputParser parser){
}
+ // Inactive material inside chamber imported form CAD drawings
+ vector<NPL::InputBlock*> blocks_material = parser.GetAllBlocksWithTokenAndValue("Strasse","InactiveMaterial");
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks_material.size() << " inactive material found " << endl;
+
+ for(unsigned int i = 0 ; i < blocks_material.size() ; i++){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Strasse Inactive material from CAD " << i+1 << endl;
+
+ if(blocks_material[i]->HasToken("Chamber")){
+ ChamberPath= blocks_material[i]->GetString("Chamber");
+ found_chamber = true;
+ }
+ if(blocks_material[i]->HasToken("Stars")){
+ StarsPath= blocks_material[i]->GetString("Stars");
+ found_stars = true;
+ }
+ if(blocks_material[i]->HasToken("Blades")){
+ BladesPath= blocks_material[i]->GetString("Blades");
+ found_blades = true;
+ }
+ if(blocks_material[i]->HasToken("Base")){
+ BasePath= blocks_material[i]->GetString("Base");
+ found_base = true;
+ }
+ }
}
@@ -848,7 +996,7 @@ void Strasse::ConstructDetector(G4LogicalVolume* world){
// Inner Barrel
for (unsigned short i = 0 ; i < m_Inner_R.size() ; i++) {
- G4ThreeVector Det_pos = G4ThreeVector(m_Inner_Shift[i],m_Inner_R[i]+0.5*Inner_PCB_Thickness,m_Inner_Z[i]) ;
+ G4ThreeVector Det_pos = G4ThreeVector(m_Inner_Shift[i],m_Inner_R[i]+0.5*Inner_PCB_Thickness+0.001,m_Inner_Z[i]) ; //0.001 offset just to avoid overlap with frame
Det_pos.rotate(-m_Inner_Phi[i],G4ThreeVector(0,0,1));
G4RotationMatrix* Rot = new G4RotationMatrix(0*deg,0*deg,m_Inner_Phi[i]);
@@ -859,7 +1007,7 @@ void Strasse::ConstructDetector(G4LogicalVolume* world){
// Outer Barrel
for (unsigned short i = 0 ; i < m_Outer_R.size() ; i++) {
- G4ThreeVector Det_pos = G4ThreeVector(m_Outer_Shift[i],m_Outer_R[i]+0.5*Inner_PCB_Thickness,m_Outer_Z[i]) ;
+ G4ThreeVector Det_pos = G4ThreeVector(m_Outer_Shift[i],m_Outer_R[i]+0.5*Inner_PCB_Thickness+0.001,m_Outer_Z[i]) ;//0.001 offset just to avoid overlap with frame
Det_pos.rotate(-m_Outer_Phi[i],G4ThreeVector(0,0,1));
G4RotationMatrix* Rot = new G4RotationMatrix(0*deg,0*deg,m_Outer_Phi[i]);
@@ -869,6 +1017,7 @@ void Strasse::ConstructDetector(G4LogicalVolume* world){
}
// Chamber
+ /*
for (unsigned short i = 0 ; i < m_Chamber_Z.size() ; i++) {
G4ThreeVector Det_pos = G4ThreeVector(0,0,-m_Chamber_Z[i]) ;
G4RotationMatrix* Rot = new G4RotationMatrix();
@@ -877,6 +1026,35 @@ void Strasse::ConstructDetector(G4LogicalVolume* world){
BuildChamber(),
"Strasse",world,false,i+1);
}
+ */
+
+
+ //G4ThreeVector Det_pos = G4ThreeVector(0,0,+11.5) ;
+ G4ThreeVector Det_pos = G4ThreeVector(0,0,0) ;
+ G4RotationMatrix* Rot = new G4RotationMatrix();
+ Rot->rotateY(270.*deg);
+ Rot->rotateX(0.*deg);
+
+ if(found_chamber){
+ new G4PVPlacement(Rot, Det_pos, BuildChamberFromCAD(ChamberPath),
+ "Strasse_Chamber",world, false, 0);
+ }
+
+ if(found_blades){
+ new G4PVPlacement(Rot, Det_pos, BuildBlades(BladesPath),
+ "Strasse_Blades",world, false, 0);
+ }
+ if(found_stars){
+ G4ThreeVector Det_pos2 = G4ThreeVector(0,0,0) ;
+ new G4PVPlacement(Rot, Det_pos2, BuildStars(StarsPath),
+ "Strasse_Stars",world, false, 0);
+ }
+ if(found_base){
+ G4ThreeVector Det_pos3 = G4ThreeVector(0,0,0) ;
+ new G4PVPlacement(Rot, Det_pos3, BuildBase(BasePath),
+ "Strasse_Base",world, false, 0);
+ }
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
@@ -1087,6 +1265,7 @@ void Strasse::InitializeMaterial(){
m_MaterialSilicon = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
m_MaterialAl = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
m_MaterialPCB = MaterialManager::getInstance()->GetMaterialFromLibrary("PCB");
+ m_MaterialCu = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
m_MaterialVacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
}
diff --git a/NPSimulation/Detectors/Strasse/Strasse.hh b/NPSimulation/Detectors/Strasse/Strasse.hh
index f2012f6ee502b60b8c250f8b7a2c858450866c8a..0250f105e330526b09a2647d71867cc38f665c5c 100644
--- a/NPSimulation/Detectors/Strasse/Strasse.hh
+++ b/NPSimulation/Detectors/Strasse/Strasse.hh
@@ -57,16 +57,29 @@ class Strasse : public NPS::VDetector{
G4LogicalVolume* BuildInnerDetector();
G4LogicalVolume* BuildOuterDetector();
G4LogicalVolume* BuildElectronic();
- G4LogicalVolume* BuildFrame();
G4LogicalVolume* BuildChamber();
+ G4LogicalVolume* BuildChamberFromCAD(string path);
+ G4LogicalVolume* BuildStars(string path);
+ G4LogicalVolume* BuildBlades(string path);
+ G4LogicalVolume* BuildBase(string path);
private:
G4LogicalVolume* m_InnerDetector;
G4LogicalVolume* m_OuterDetector;
G4LogicalVolume* m_Electronic;
- G4LogicalVolume* m_Frame;
+ G4LogicalVolume* m_Stars;
G4LogicalVolume* m_Chamber;
-
+ G4LogicalVolume* m_Blades;
+ G4LogicalVolume* m_Base;
+
+ string ChamberPath;
+ string BasePath;
+ string StarsPath;
+ string BladesPath;
+ bool found_chamber;
+ bool found_blades;
+ bool found_stars;
+ bool found_base;
private:
// Initialize material used in detector definition
@@ -78,6 +91,7 @@ class Strasse : public NPS::VDetector{
G4Material* m_MaterialAl ;
G4Material* m_MaterialVacuum ;
G4Material* m_MaterialPCB ;
+ G4Material* m_MaterialCu ;
// calculated dimension
double m_Active_InnerWafer_Width;
@@ -141,8 +155,6 @@ class Strasse : public NPS::VDetector{
vector<double> m_Chamber_Z;
- // Visualisation Attribute
- //G4VisAttributes* m_VisTrap;
// Needed for dynamic loading of the library
public:
@@ -153,9 +165,10 @@ class Strasse : public NPS::VDetector{
G4VisAttributes* SiliconVisAtt ;
G4VisAttributes* PCBVisAtt;
G4VisAttributes* PADVisAtt ;
- G4VisAttributes* FrameVisAtt ;
+ G4VisAttributes* StarsVisAtt ;
G4VisAttributes* ChamberVisAtt ;
G4VisAttributes* GuardRingVisAtt ;
+ G4VisAttributes* BladeVisAtt ;
};
diff --git a/NPSimulation/Scorers/InteractionScorers.cc b/NPSimulation/Scorers/InteractionScorers.cc
index 1ab9bed3eee5aae3836793c0927582ab22cdbb4e..00f040f1ecfc58a711fa0c45cdd064ca7ca31954 100644
--- a/NPSimulation/Scorers/InteractionScorers.cc
+++ b/NPSimulation/Scorers/InteractionScorers.cc
@@ -66,7 +66,7 @@ G4bool PS_Interactions::ProcessHits(G4Step* aStep, G4TouchableHistory*){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_Interactions::Initialize(G4HCofThisEvent*){
// Clear is called by EventAction
-
+ clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -76,7 +76,6 @@ void PS_Interactions::EndOfEvent(G4HCofThisEvent*){
for(unsigned int i = 0 ; i < size ; i++)
m_InterractionCoordinates->SetInteraction(m_DataVector[i]->GetEnergy(),m_DataVector[i]->GetTime(),m_DataVector[i]->GetPositionX(),m_DataVector[i]->GetPositionY(),m_DataVector[i]->GetPositionZ(),m_DataVector[i]->GetTheta()/deg,m_DataVector[i]->GetPhi()/deg);
- clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/Projects/ComptonTelescope/online/src/DecodeD.cpp b/Projects/ComptonTelescope/online/src/DecodeD.cpp
index 3ef51c8149014f24bed107e60cabd8fb34956b30..3a09dd750da7436982bdd3c75c199e770ad85892 100644
--- a/Projects/ComptonTelescope/online/src/DecodeD.cpp
+++ b/Projects/ComptonTelescope/online/src/DecodeD.cpp
@@ -100,7 +100,7 @@ __EVENTTYPE__* DecodeD::getEvent()
void DecodeD::decodeEvent()
{
- this -> Clear();
+ /*this -> */Clear();
switch (datatype) {
case D_ROOT:
if (cursor < length) {
diff --git a/Projects/ComptonTelescope/online/src/online_coinc.cpp b/Projects/ComptonTelescope/online/src/online_coinc.cpp
index a8ad64d2676ed419b98101a0e2e2cd48aae51346..3b6cedac406382013936dbc79bb7685b8364181c 100644
--- a/Projects/ComptonTelescope/online/src/online_coinc.cpp
+++ b/Projects/ComptonTelescope/online/src/online_coinc.cpp
@@ -16,6 +16,10 @@
#include "DecodeD.h"
#include "DecodeT.h"
+#define __RUN__ 3
+
+#define __1DET__
+
#define __TEST_ZONE__
#undef __TEST_ZONE__
@@ -105,15 +109,47 @@ int main(int argc, char** argv)
auto ccamPhys = (TComptonTelescopePhysics*) m_NPDetectorManager->GetDetector("ComptonTelescope");
ccamPhys->SetRawDataPointer(ccamData);
+#ifdef __1DET__
+ ifstream is;
+ is.open("/disk/proto-data/data/20210510_Bi207/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-05-10_07_41_50.raw");
+ is.seekg(0, ios::end);
+ int len = is.tellg();
+ is.seekg(0, ios::beg);
+ char* buff = new char[len];
+ is.read(buff, len);
+ is.close();
+ DecodeR* DR = new DecodeR(false, buff);
+ while(DR -> getCursor() < len) {
+ DR -> decodeBlobMFM();
+ m_NPDetectorManager->ClearEventPhysics();
+ m_NPDetectorManager->ClearEventData();
+ ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), 64);
+ m_NPDetectorManager->BuildPhysicalEvent();
+ m_OutputTree->Fill();
+ m_NPDetectorManager->CheckSpectraServer();
+ }
+ delete DR;
+ delete [] buff;
+ m_NPDetectorManager -> WriteSpectra();
+#else
+
// read data file/flux and fill ccamData object
std::cout << "Reading data\n";
DecodeR* DR = new DecodeR(false); // Instantiates DecodeR object reading calorimeter data flux
DecodeT* DT = new DecodeT(false); // Instantiates DecodeT object reading trigger data flux
DecodeD* DD = new DecodeD(false); // Instantiates DecodeD object reading DSSSD(s) data flux
// newframe_t* event;
- //DD -> setTree("/disk/proto-data/data/20210304_run2/bb7_3309-7_cs137_20210304_14h35_conv.root");
+ #if __RUN__ == 0
+ DD -> setTree("../data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root");
+ #elif __RUN__ == 1
+ DD -> setTree("/disk/proto-data/data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root");
+ #elif __RUN__ == 2
+ DD -> setTree("/disk/proto-data/data/20210304_run2/bb7_3309-7_cs137_20210304_14h35_conv.root");
+ #elif __RUN__ == 3
DD -> setTree("/disk/proto-data/data/20210305_run3/bb7_3309-7_cs137_20210305_14h53_conv.root");
- //DD -> setTree("../data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root");
+ #elif __RUN__ == 4
+ DD -> setTree("/disk/proto-data/data/20210407_run4/bb7_3309-7_cs137_20210407_14h53_conv.root");
+ #endif
int dlen = DD -> getLength();
int i = 0;// ROSMAP files loop counter
@@ -126,8 +162,16 @@ int main(int argc, char** argv)
ifstream iros, itrig;
cout << "Loading data files " << std::flush;
+ #if __RUN__ == 0
+ itrig.open("../data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary);
+ #elif __RUN__ == 1
+ itrig.open("/disk/proto-data/data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary);
+ #elif __RUN__ == 2
+ #elif __RUN__ == 3
itrig.open("/disk/proto-data/data/20210305_run3/mfm_trigger_20210305_run3.raw", ios::binary);
- //itrig.open("../data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary);
+ #elif __RUN__ == 4
+ itrig.open("/disk/proto-data/data/20210407_run4/mfm_trigger_20210407_run4.raw", ios::binary);
+ #endif
itrig.seekg(0, ios::end);
int tlen = itrig.tellg();
itrig.seekg(0, ios::beg);
@@ -136,8 +180,16 @@ int main(int argc, char** argv)
itrig.close();
cout << "... " << std::flush;
+ #if __RUN__ == 0
+ iros.open("../data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary);
+ #elif __RUN__ == 1
+ iros.open("/disk/proto-data/data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary);
+ #elif __RUN__ == 2
+ #elif __RUN__ == 3
iros.open("/disk/proto-data/data/20210305_run3/mfm_rosmap_20210305_run3.raw", ios::binary);
- //iros.open("../data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary);
+ #elif __RUN__ == 4
+ iros.open("/disk/proto-data/data/20210407_run4/mfm_rosmap_20210407_run4.raw", ios::binary);
+ #endif
iros.seekg(0, ios::end);
int rlen = iros.tellg();
iros.seekg(0, ios::beg);
@@ -157,6 +209,22 @@ int main(int argc, char** argv)
resetCount = DT->getResetCount() - resetCount;
resetCount = -resetCount; // T B C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
cout << "Found reset count: " << resetCount << endl;
+
+ if (resetCount == 0) {
+ string ans;
+ cout << "reset count is 0. Continue ? y/[n]";
+ cin >> ans;
+ if (ans != "y" and ans != "Y") {
+ DT -> setRaw(tbuff);
+ DT -> decodeBlobMFM();
+ resetCount = DT -> getResetCount();
+ while (not(DT->hasTrigged(2))) {
+ DT -> decodeBlobMFM();
+ }
+ resetCount = DT->getResetCount() - resetCount;
+ cout << "Found reset count: " << resetCount << endl;
+ }
+ }
int cr, cd, tr, td;
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
@@ -261,7 +329,7 @@ int main(int argc, char** argv)
DR -> decodeBlobMFM();
//#ifdef __TEST_ZONE__
- cout << "Entering test zone." << endl;
+// cout << "Entering test zone." << endl;
//#else
DD -> rewind();
@@ -271,7 +339,7 @@ int main(int argc, char** argv)
tr = DR -> getTime();
td = DD -> getTime();
- while(DR -> getCursor() < rlen and DD -> getCursor() < dlen)
+ while(DR -> getCursor() < rlen and DD -> getCursor() < dlen and c < 10000)
{
if (cr == cd) {
#ifndef __TEST_ZONE__
@@ -283,7 +351,6 @@ int main(int argc, char** argv)
//DR -> Dump();
//DD -> Dump();
c++;
- cout << cc << " " << c /*<< "(" << cr << ", " << cd << ") : " << tr << " " << td*/ << "\n";
// Clear raw and physics data
m_NPDetectorManager->ClearEventPhysics();
m_NPDetectorManager->ClearEventData();
@@ -307,6 +374,7 @@ int main(int argc, char** argv)
m_OutputTree->Fill();
cc++;
#endif
+ cout << cc << " " << c /*<< "(" << cr << ", " << cd << ") : " << tr << " " << td*/ << " |\t";
}
// check spectra
@@ -353,7 +421,7 @@ int main(int argc, char** argv)
deltaT->Write();
}
fout->Close();
-
+#endif
// Essential
#if __cplusplus > 199711L && NPMULTITHREADING
m_NPDetectorManager->StopThread();
diff --git a/Projects/Strasse/geometry/strasse_CAD.detector b/Projects/Strasse/geometry/strasse_CAD.detector
new file mode 100644
index 0000000000000000000000000000000000000000..1acb1346fd57bcc8452680ae3d0d678122e51f2d
--- /dev/null
+++ b/Projects/Strasse/geometry/strasse_CAD.detector
@@ -0,0 +1,101 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 150 mm
+ ANGLE= 0 deg
+ RADIUS= 15 mm
+ MATERIAL= LH2
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+ NbSlices= 10
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Strasse Info
+ Inner_Wafer_Length= 124 mm
+ Inner_Wafer_Width= 32 mm
+ Inner_Wafer_Thickness= 200 micrometer
+ Inner_Wafer_AlThickness= 0.4 micrometer
+ Inner_Wafer_PADExternal= 0 mm
+ Inner_Wafer_PADInternal= 0 mm
+ Inner_Wafer_GuardRing= 1.0 mm
+ Inner_Wafer_TransverseStrips= 610
+ Inner_Wafer_LongitudinalStrips= 150
+ Inner_PCB_PortWidth= 1 mm
+ Inner_PCB_StarboardWidth= 1 mm
+ Inner_PCB_BevelAngle= 90 deg
+ Inner_PCB_UpstreamWidth= 12 mm
+ Inner_PCB_DownstreamWidth= 38 mm
+ Inner_PCB_MidWidth= 1 mm
+ Inner_PCB_Thickness= 2.4 mm
+ Inner_PCB_Ledge= 1 mm
+ Inner_PCB_Step= 2 mm
+ Outer_Wafer_Length= 123 mm
+ Outer_Wafer_Width= 64.6 mm
+ Outer_Wafer_Thickness= 300 micrometer
+ Outer_Wafer_AlThickness= 0.4 micrometer
+ Outer_Wafer_PADExternal= 0 mm
+ Outer_Wafer_PADInternal= 0 mm
+ Outer_Wafer_GuardRing= 1.0 mm
+ Outer_PCB_PortWidth= 1 mm
+ Outer_PCB_StarboardWidth= 1 mm
+ Outer_PCB_BevelAngle= 90 deg
+ Outer_PCB_UpstreamWidth= 40 mm
+ Outer_PCB_DownstreamWidth= 12 mm
+ Outer_PCB_MidWidth= 1 mm
+ Outer_PCB_Thickness= 2.4 mm
+ Outer_PCB_Ledge= 1 mm
+ Outer_PCB_Step= 2 mm
+ Outer_Wafer_TransverseStrips= 605
+ Outer_Wafer_LongitudinalStrips= 313
+ % unused if using CAD file (.stl) chamber
+ Chamber_Thickness= 3 mm
+ Chamber_Cylinder_Length= 360 mm
+ Chamber_Radius= 180 mm
+ Chamber_ExitTube_Radius= 79.5 mm
+ Chamber_ExitTube_Length= 100 mm
+ Chamber_Flange_Inner_Radius= 50 mm
+ Chamber_Sphere_Radius= 220 mm
+ Chamber_Sphere_Shift= 60 mm
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias InnerPhi
+ Action= Copy
+% Value= 116
+ Value= -4 56 116 176 236 296
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Inner
+ Radius= 27.6 mm
+ Z= 76.5 mm
+ Phi= @InnerPhi deg
+ Shift= 3 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias OuterPhi
+ Action= Copy
+% Value= 0
+ Value= 0 60 120 180 240 300
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Outer
+ Radius= 58.5 mm
+ Z= 76.5 mm
+ Phi= @OuterPhi deg
+ Shift= 0 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse InactiveMaterial
+ Chamber= ./geometry/STRASSE_Chamber.stl
+ Stars= ./geometry/STRASSE_StarSupports.stl
+ Base= ./geometry/STRASSE_Base.stl
+ Blades= ./geometry/STRASSE_Blades.stl
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+%Catana CSV
+% Path= geometry/Catana.csv
+% Pos= 0 0 100 mm
+% Rshift= 100 micrometer
+
+