uLib/src/HEP/Geant/Solid.cpp

/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
  All rights reserved

  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >

  ------------------------------------------------------------------
  This library is free software;  you  can  redistribute  it  and/or
  modify it  under the  terms  of  the  GNU  Lesser  General  Public
  License as published  by  the  Free  Software  Foundation;  either
  version 3.0 of the License, or (at your option) any later version.

  This library is  distributed in  the hope that it will  be useful,
  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of  the GNU Lesser General  Public
  License along with this library.

//////////////////////////////////////////////////////////////////////////////*/

// G4 Solid //
#include <CLHEP/Units/SystemOfUnits.h>
#include <Geant4/G4LogicalVolume.hh>
#include <Geant4/G4Material.hh>
#include <Geant4/G4NistManager.hh>

// Tessellated solid //
#include <Geant4/G4TessellatedSolid.hh>
#include <Geant4/G4ThreeVector.hh>
#include <Geant4/G4TriangularFacet.hh>
#include <Geant4/G4Box.hh>
#include <Geant4/G4PVPlacement.hh>


#include "Math/Dense.h"
#include "Math/Transform.h"

#include "Solid.h"

namespace uLib {
namespace Geant {

class DetectorsSolidImpl {
public:
  static G4ThreeVector getG4Vector3f(const Vector3f &vector) {
    return G4ThreeVector(vector(0), vector(1), vector(2));
  }

  
};

Solid::Solid()
    : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), 
      m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {}

Solid::Solid(const char *name)
    : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), 
      m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {}

Solid::~Solid() {
  if (m_Position) delete m_Position;
  if (m_Rotation) delete m_Rotation;
}

void Solid::SetNistMaterial(const char *name) {
  G4NistManager *nist = G4NistManager::Instance();
  G4Material *mat = nist->FindOrBuildMaterial(name);
  if (mat) SetMaterial(mat);
}

void Solid::SetMaterial(G4Material *material) {
  if (material) {
    m_Material = material;
    if (m_Logical) {
      m_Logical->SetMaterial(material);
    } else if (GetG4Solid()) {
      m_Logical = new G4LogicalVolume(GetG4Solid(), m_Material, GetName());
    }
  }
}

void Solid::SetTransform(Matrix4f transform) {
    uLib::AffineTransform t;
  t.SetMatrix(transform);

    // 2. Extract position and rotation for Geant4
    Vector3f pos = t.GetPosition();
    if (!m_Position) m_Position = new G4ThreeVector();
    *m_Position = G4ThreeVector(pos(0), pos(1), pos(2));

    // Create a G4 rotation matrix from the 4x4 matrix
    Matrix3f m = t.GetRotation();
    if (!m_Rotation) m_Rotation = new G4RotationMatrix();
    m_Rotation->set(G4ThreeVector(m(0,0), m(1,0), m(2,0)), 
                    G4ThreeVector(m(0,1), m(1,1), m(2,1)), 
                    G4ThreeVector(m(0,2), m(1,2), m(2,2))); 

    // 3. If object is already placed, update its transformation
    if (m_Physical) {
        m_Physical->SetTranslation(*m_Position);
        m_Physical->SetRotation(m_Rotation);
    }
}

void Solid::SetParent(Solid *parent) {
  if (!m_Logical) {
    std::cerr << "logical volume not created for solid " << GetName() << std::endl;
    return;
  }

  if(m_Physical) {
    std::cerr << "physical volume already created for solid " << GetName() << std::endl;
    return;
  }

  G4LogicalVolume* parentLogical = nullptr;
  if (parent) {
    parentLogical = parent->GetLogical();
    if (!parentLogical) {
      std::cerr << "parent logical volume not created for solid " << parent->GetName() << std::endl;
      return;
    }
  }

  // G4PVPlacement
  m_Physical = new G4PVPlacement(
      m_Rotation,            // Rotation
      *m_Position,           // Position (translation) inside the parent
      m_Logical,             // The logical volume of this solid (the child)
      GetName(),             // Name of the physical volume
      parentLogical,         // The logical volume of the parent (nullptr if it's the World volume)
      false,                 // Boolean operations (usually false)
      0,                     // Copy number
      true                   // Check overlaps (useful to enable in debug phase)
  );
}







TessellatedSolid::TessellatedSolid() 
    : BaseClass("unnamed_tessellated"), m_Solid(new G4TessellatedSolid("unnamed_tessellated")) {}

TessellatedSolid::TessellatedSolid(const char *name)
    : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) {
    }

void TessellatedSolid::SetMesh(TriangleMesh &mesh) {
  this->m_Mesh = mesh;
  G4TessellatedSolid *ts = this->m_Solid;
  for (int i = 0; i < mesh.Triangles().size(); ++i) {
    const Vector3i &trg = mesh.Triangles().at(i);
    G4TriangularFacet *facet = new G4TriangularFacet(
        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(0))),
        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(1))),
        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(2))), ABSOLUTE);
    ts->AddFacet((G4VFacet *)facet);
  }
  if (this->m_Logical) {
    this->m_Logical->SetSolid(ts);
  }
}

void TessellatedSolid::Update() {
}





BoxSolid::BoxSolid(const char *name) : 
    BaseClass(name), 
    m_ContainerBox(new ContainerBox()), 
    m_Solid(new G4Box(name, 1, 1, 1))
  {}

BoxSolid::BoxSolid(const char *name, ContainerBox *box) : BaseClass(name) {
  m_Solid = new G4Box(name, 1, 1, 1);
  m_ContainerBox = box;
  Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update);
  if (m_Logical) {
    m_Logical->SetSolid(m_Solid);
  }
  Update();
}

void BoxSolid::Update() {
  if (m_ContainerBox) {
    Vector3f size = m_ContainerBox->GetSize();
    m_Solid->SetXHalfLength(size(0) * 0.5);
    m_Solid->SetYHalfLength(size(1) * 0.5);
    m_Solid->SetZHalfLength(size(2) * 0.5);
    
    // Geant4 placement is relative to center. uLib Box is anchored at corner.
    // 1. Get position and rotation (clean, without scale)
    Vector3f pos = m_ContainerBox->GetPosition();
    Matrix3f rot = m_ContainerBox->GetRotation();

    // 2. Center = Corner + Rotation * (Half-Size)
    // We must rotate the offset vector because uLib box can be rotated.
    Vector3f center = pos + rot * (size * 0.5);

    uLib::AffineTransform t;
    t.SetPosition(center);
    t.SetRotation(rot);
    
    this->SetTransform(t.GetMatrix());
  }



}


} // namespace Geant
} // namespace uLib