/*
 * Fire some particles at a block of water.
 */
#include <numeric>

#if defined(THREADED)
#include "G4MTRunManager.hh"
#else
#include "G4RunManager.hh"
#endif
#include "G4VUserDetectorConstruction.hh"
#include "G4NistManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4Box.hh"
#include "G4PVPlacement.hh"
#include "G4PhantomParameterisation.hh"
#include "G4PVParameterised.hh"
#include "G4PhysListFactory.hh"
#include "G4VUserActionInitialization.hh"
#include "G4UserRunAction.hh"
#include "G4Run.hh"
#include "G4VUserPrimaryGeneratorAction.hh"
#include "G4GeneralParticleSource.hh"
#include "G4UImanager.hh"

class Detector : public G4VUserDetectorConstruction
{
   // 1x1x1mm voxels in a 10x10x10cm mesh.
   const G4double VOXEL_X = 1.0 * mm;
   const G4double VOXEL_Y = 1.0 * mm;
   const G4double VOXEL_Z = 1.0 * mm;
   const int    N_VOXEL_X = 100;
   const int    N_VOXEL_Y = 100;
   const int    N_VOXEL_Z = 100;
   const int    N_VOXELS  = N_VOXEL_X * N_VOXEL_Y * N_VOXEL_Z;
   const G4double TARGET_X = N_VOXEL_X * VOXEL_X;
   const G4double TARGET_Y = N_VOXEL_Y * VOXEL_Y;
   const G4double TARGET_Z = N_VOXEL_Z * VOXEL_Z;

public:
   Detector ( void )
   {}
   virtual ~Detector ( void )
   {}

   virtual G4VPhysicalVolume* Construct( void )
   {
      // Get nist material manager and the usual materials.
      G4NistManager* nist = G4NistManager::Instance();

      // Really don't need the special densities because these are
      // close to the defaults.
      G4Material*   air = nist->FindOrBuildMaterial( "G4_AIR" );
      G4Material* water = nist->FindOrBuildMaterial( "G4_WATER" );

      G4Material* background = air;
      G4Material* target     = water;

      // The bounding box of the simulation.
      G4double bbox_x = 2 * m;
      G4double bbox_y = 2 * m;
      G4double bbox_z = 2 * m;
      G4Box* bbox_geom =
         new G4Box( "World",
                    0.5 * bbox_x, 0.5 * bbox_y, 0.5 * bbox_z );
      G4LogicalVolume* bbox_log =
         new G4LogicalVolume( bbox_geom,  // Geometry of bbox
                              background, // Material of bbox
                              "World" );  // Name of bbox
      G4VPhysicalVolume* bbox_phys =
         new G4PVPlacement( 0,               // no rotation
                            G4ThreeVector(), // at (0,0,0)
                            bbox_log,        // shape and material
                            "World",         // Name
                            0,               // Mother (none)
                            false,           // no boolean operation
                            0,               // copy number
                            false );         // no overlap checking


      // And the target (detector, phantom).
      G4Box* target_geom =
         new G4Box( "Target",
                    0.5 * TARGET_X, 0.5 * TARGET_Y, 0.5 * TARGET_Z );
      G4LogicalVolume* target_log =
         new G4LogicalVolume( target_geom, // Geometry of target
                              target,      // Material of target [unused?]
                              "Target" );  // Name of target
      // I guess its parent manages its lifetime.
      G4VPhysicalVolume* target_phys =
         new G4PVPlacement( 0,               // no rotation
                            // starts at x=0 with centered yz voxel row
                            G4ThreeVector(0.5*TARGET_X,0.*mm,0.*mm),
                            target_log,      // shape and material
                            "Target",        // Name
                            bbox_log,        // Mother (bounding box)
                            false,           // no boolean operation
                            0,               // copy number
                            false );         // no overlap checking

      // Try to turn the Target in to a bunch of voxels.

      G4Box* voxel_solid =
         new G4Box( "Voxel",
                    0.5 * VOXEL_X, 0.5 * VOXEL_Y, 0.5 * VOXEL_Z );
      // We need this in the detector sensitivity method.
      voxel_log =
         new G4LogicalVolume( voxel_solid, target, "Voxel", 0, 0, 0 );
      // This is the nub: Create a parameterization object with the
      // operative layout of voxels.
      G4PhantomParameterisation* params =
         new G4PhantomParameterisation;
      params->SetVoxelDimensions( 0.5 * VOXEL_X, 0.5 * VOXEL_Y, 0.5 * VOXEL_Z );
      params->SetNoVoxel( N_VOXEL_X, N_VOXEL_Y, N_VOXEL_Z );
      // It's all water.
      std::vector<G4Material*> materials{water};
      params->SetMaterials( materials );
      params->BuildContainerSolid( target_phys );
      params->CheckVoxelsFillContainer( target_geom->GetXHalfLength(),
                                        target_geom->GetYHalfLength(),
                                        target_geom->GetZHalfLength() );

      params->SetSkipEqualMaterials( false );

      G4PVParameterised* phantom_phys =
         new G4PVParameterised( "phantom", voxel_log, target_log,
                                kXAxis, N_VOXELS, params );

      // And I guess this is the most important step to use
      // G4RegularNavigation.
      phantom_phys->SetRegularStructureId( 1 );

      return bbox_phys;
   }
private:
   G4LogicalVolume* voxel_log;
};

class PrimaryAction : public G4VUserPrimaryGeneratorAction
{
public:
   PrimaryAction ( void )
      : G4VUserPrimaryGeneratorAction(),
        gps_( nullptr )
   {
      gps_ = new G4GeneralParticleSource();
   }
   virtual ~PrimaryAction ( void )
   {
      delete gps_;
   }
   virtual void GeneratePrimaries ( G4Event* event )
   {
      gps_->GeneratePrimaryVertex( event );
   }
private:
   G4GeneralParticleSource* gps_;
};

// Maybe this is starting to become slightly more clear.
class Run : public G4Run
{
public:
   Run ( void )
   {}
   ~Run ( void )
   {}
};

class RunAction : public G4UserRunAction
{
public:
   RunAction ()
   {}
   virtual ~RunAction ( void )
   {}
   virtual G4Run* GenerateRun ( void )
   {
      G4Run* run = new Run;
      return run;
   }
};

class ActionInitialization : public G4VUserActionInitialization
{
public:
   ActionInitialization ( void )
      : G4VUserActionInitialization()
   {}
   ~ActionInitialization ( void )
   {}
   virtual void BuildForMaster ( void ) const
   {
      SetUserAction( new RunAction );
   }
   virtual void Build ( void ) const
   {
      SetUserAction( new PrimaryAction );
      SetUserAction( new RunAction );
   }
};

int main ( int argc, char* argv[] )
{
   G4String run_file( "electrons.txt" );
   if ( argc > 1 ) {
      run_file = argv[1];
   }

   G4cout << "Run file: \""  << run_file << "\"" << G4endl;

#if defined(THREADED)
   G4MTRunManager* run_manager = new G4MTRunManager;
#else
   G4RunManager* run_manager = new G4RunManager;
#endif

   G4PhysListFactory factory;
   auto physics_list = factory.GetReferencePhysList( "QGSP_BIC_EMV" );
   run_manager->SetUserInitialization( physics_list );
   run_manager->SetUserInitialization( new Detector );
   run_manager->SetUserInitialization( new ActionInitialization );

   // It appears that there's nothing you can't do with the macro
   // file, at least as far as running the problem.
   G4UImanager* ui = G4UImanager::GetUIpointer();

   ui->ApplyCommand( "/control/execute " + run_file );

   delete run_manager;

   return 0;
}