BeamabortCreator.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#include <beast/beamabort/geometry/BeamabortCreator.h>
#include <beast/beamabort/simulation/SensitiveDetector.h>
 
#include <geometry/Materials.h>
#include <geometry/CreatorFactory.h>
#include <framework/gearbox/GearDir.h>
#include <framework/logging/Logger.h>
 
#include <cmath>
#include <boost/format.hpp>
#include <boost/foreach.hpp>
#include <boost/algorithm/string.hpp>
 
#include <G4LogicalVolume.hh>
#include <G4PVPlacement.hh>
 
//Shapes
#include <G4Box.hh>
#include "G4SubtractionSolid.hh"
#include <G4UserLimits.hh>
 
//Visualization
#include "G4Colour.hh"
#include <G4VisAttributes.hh>
 
using namespace std;
using namespace boost;
 
namespace Belle2 {
  using namespace geometry;
 
  namespace beamabort {
 
    // Register the creator
    geometry::CreatorFactory <BeamabortCreator> BeamabortFactory("BEAMABORTCreator");
 
    BeamabortCreator::BeamabortCreator()
    {}
 
    BeamabortCreator::~BeamabortCreator()
    {}
 
    void BeamabortCreator::createGeometry(G4LogicalVolume& topVolume, GeometryTypes)
    {
      //get parameters from .xml file
      std::string prep = "Active.";
 
      B2INFO("BeamabortCreator phase3");
      //Visualization Attributes
      G4VisAttributes* orange = new G4VisAttributes(G4Colour(1, 2, 0));
      orange->SetForceAuxEdgeVisible(true);
      G4VisAttributes* magenta = new G4VisAttributes(G4Colour(1, 0, 1));
      magenta->SetForceAuxEdgeVisible(true);
      //lets get the stepsize parameter with a default value of 5 µm
      double stepSize = m_config.getParameter("stepSize", 5 * CLHEP::um);
 
      //Lets loop over all the Active nodes
      int phase = m_config.getParameter(prep + "phase");
      //Positioned Diamonds
      std::vector<double> x_pos;
      std::vector<double> y_pos;
      std::vector<double> z_pos;
      std::vector<double> phi;
      std::vector<double> thetaX;
      std::vector<double> thetaY;
      std::vector<double> thetaZ;
      std::vector<double> svdAngle;
      std::vector<double> r;
      std::vector<double> deltaX;
      unsigned int dimz = 0;
 
      for (double z : m_config.getParArray(prep + "z", {0})) {
        z *= CLHEP::cm;
        z_pos.push_back(z);
        dimz++;
      }
 
      for (double ThetaZ : m_config.getParArray(prep + "ThetaZ", {0})) {
        thetaZ.push_back(ThetaZ);
      }
 
      if (thetaZ.size() != dimz) { B2ERROR("Diamond data not consistent (i.e. not same number of all position parmeters)"); return;}
 
      if (phase == 2 || phase == 3) {
 
        for (double Phi : m_config.getParArray(prep + "Phi", {0})) {
          phi.push_back(Phi);
        }
 
        for (double r_dia : m_config.getParArray(prep + "r_dia", {0})) {
          r_dia *= CLHEP::cm;
          r.push_back(r_dia);
        }
        for (double dX : m_config.getParArray(prep + "deltaX", {0})) {
          dX *= CLHEP::cm;
          deltaX.push_back(dX);
        }
        for (double addAngle : m_config.getParArray(prep + "addAngle", {0})) {
          svdAngle.push_back(addAngle);
        }
        if (phi.size() != dimz || r.size() != dimz || deltaX.size() != dimz || svdAngle.size() != dimz) { B2ERROR("Diamond data not consistent (i.e. not same number of all position parmeters)"); return;}
      }
      if (phase == 1) {
        for (double x : m_config.getParArray(prep + "x", {0})) {
          x *= CLHEP::cm;
          x_pos.push_back(x);
        }
        for (double y : m_config.getParArray(prep + "y", {0})) {
          y *= CLHEP::cm;
          y_pos.push_back(y);
        }
 
        for (double ThetaX : m_config.getParArray(prep + "ThetaX", {0})) {
          thetaX.push_back(ThetaX);
        }
        for (double ThetaY : m_config.getParArray(prep + "ThetaY", {0})) {
          thetaY.push_back(ThetaY);
        }
        if (x_pos.size() != dimz || y_pos.size() != dimz || thetaX.size() != dimz || thetaY.size() != dimz) { B2ERROR("Diamond data not consistent (i.e. not same number of all position parmeters)"); return;}
      }
 
      //create beamabort package
      G4double dx_opa = 12. / 2. * CLHEP::mm;
      G4double dy_opa = 18. / 2. * CLHEP::mm;
      G4double dz_opa = 3.1 / 2. * CLHEP::mm;
      G4VSolid* s_pa = new G4Box("s_opa", dx_opa, dy_opa, dz_opa);
      G4double dx_ipa = 6. / 2. * CLHEP::mm;
      G4double dy_ipa = 6. / 2. * CLHEP::mm;
      G4double dz_ipa = 5.6 / 2. * CLHEP::mm;
      G4VSolid* s_ipa = new G4Box("s_ipa", dx_ipa, dy_ipa, dz_ipa);
      s_pa = new G4SubtractionSolid("s_pa", s_pa, s_ipa, 0, G4ThreeVector(0., 4.0, 0.));
      G4LogicalVolume* l_pa = new G4LogicalVolume(s_pa, G4Material::GetMaterial("Al6061"), "l_pa");
      l_pa->SetVisAttributes(magenta);
      G4Transform3D transform;
      for (unsigned int i = 0; i < dimz; i++) {
        B2INFO("DIA-" << i << "RotateZ3D phi: " << phi[i]);
 
        if (phase == 1) {
          transform = G4Translate3D(x_pos[i], y_pos[i], z_pos[i]) * G4RotateX3D(thetaX[i]) *
                      G4RotateY3D(thetaY[i]) * G4RotateZ3D(thetaZ[i]);
        } else if (phase == 2 || phase == 3) {
          transform = G4Translate3D(deltaX[i], 0, 0) * G4RotateZ3D(phi[i]) *
                      G4Translate3D(r[i], 0, z_pos[i]) * G4RotateX3D(-M_PI / 2 - thetaZ[i]) *
                      G4RotateZ3D(svdAngle[i]) * G4RotateY3D(M_PI / 2);
        }
        new G4PVPlacement(transform, l_pa, TString::Format("p_dia_pa_%d", i).Data(), &topVolume, false, 0);
        B2INFO("DIA-" << i << " placed at: " << transform.getTranslation() << " mm ");
      }
 
      //create beamabort volumes
      G4double dx_ba = 4.5 / 2. * CLHEP::mm;
      G4double dy_ba = 4.5 / 2. * CLHEP::mm;
      G4double dz_ba = 0.5 / 2. * CLHEP::mm;
      G4Box* s_BEAMABORT = new G4Box("s_BEAMABORT", dx_ba, dy_ba, dz_ba);
      G4LogicalVolume* l_BEAMABORT = new G4LogicalVolume(s_BEAMABORT, geometry::Materials::get("Diamond"), "l_BEAMABORT");
      l_BEAMABORT->SetVisAttributes(orange);
 
      //Lets limit the Geant4 stepsize inside the volume
      l_BEAMABORT->SetUserLimits(new G4UserLimits(stepSize));
      l_BEAMABORT->SetVisAttributes(orange);
 
      // make the volume sensitive if the beam background study flag is enabled
      if (m_config.getParameter("BeamBackgroundStudy")) l_BEAMABORT->SetSensitiveDetector(new SensitiveDetector());
 
      for (unsigned int i = 0; i < dimz; i++) {
 
        if (phase == 1)
          transform = G4Translate3D(x_pos[i], y_pos[i],
                                    z_pos[i]) * G4RotateX3D(thetaX[i]) * G4RotateY3D(thetaY[i]) * G4RotateZ3D(thetaZ[i]);
        if (phase == 2 || phase == 3)
          transform = G4Translate3D(deltaX[i], 0, 0) * G4RotateZ3D(phi[i]) *
                      G4Translate3D(r[i], 0, z_pos[i]) * G4RotateX3D(-M_PI / 2 - thetaZ[i]) *
                      G4RotateZ3D(svdAngle[i]) * G4RotateY3D(M_PI / 2) * G4Translate3D(0, 4.2, 0);
 
        new G4PVPlacement(transform, l_BEAMABORT, TString::Format("p_dia_%d", i).Data(), &topVolume, false, i);
 
        B2INFO("DIA-sensitive volume-" << i << " placed at: " << transform.getTranslation() << " mm "
               << "  at phi angle = " << phi[i] << "   at theta angle = "
               << thetaZ[i]);
        B2INFO("DIA-sensitive volume-" << i << " G4RotateZ3D of: phi= " << phi[i] << "  G4RotateX3D = "
               << (-M_PI / 2 - thetaZ[i]));
      }
    }
  }
}