/**
 * @file SolidUtils.cc
 * @author: Ibrahin Pinera
 * @date 2017 SoLid - University of Antwerp
 */


#include "SolidUtils.hh"
#include "SolidMassComputation.hh"

#include "G4Material.hh"
#include "G4VPhysicalVolume.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4TransportationManager.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Trd.hh"
#include "G4PVPlacement.hh"
#include "G4IntersectionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4SystemOfUnits.hh"


// Method to check if given volume exist in the geometry
G4bool FindVolume(G4String VolName)
{
  G4VPhysicalVolume *tempPV      = NULL;
  G4PhysicalVolumeStore *PVStore = 0;
  G4String theRequiredVolumeName = VolName;
  PVStore      = G4PhysicalVolumeStore::GetInstance();
  G4int      i = 0;
  G4bool found = false;
  while (!found && i<G4int(PVStore->size())) {
    tempPV = (*PVStore)[i];
    found  = tempPV->GetName() == theRequiredVolumeName;
    if (!found)
      {i++;}
  }
  // found = true then the volume exists else it doesnt.
  if(found == true)
    {
      return true;
    }
  else
    {
      return false;
    }

}

// Method to get the process Id
G4int GetProcessID(G4String procname)
{
  if (procname == "Generator") return 0;
  //neutron/anti-neutron
  else if (procname == "nCapture") return 1;
  else if (procname == "nFission") return 2;
  else if (procname == "neutronInelastic") return 3;
  else if (procname == "anti_neutronInelastic") return 4;  

  //e+/e-
  else if (procname == "eIoni") return 10;
  else if (procname == "eBrem") return 11;
  else if (procname == "annihil") return 12;
  else if (procname == "electronNuclear") return 13;
  else if (procname == "positronNuclear") return 14;
  
  //phot
  else if (procname == "conv") return 20;
  else if (procname == "phot") return 21;
  else if (procname == "compt") return 22;

  //muons  
  else if (procname == "muIoni") return 30;
  else if (procname == "muPairProd") return 31;
  else if (procname == "muonNuclear") return 32;
  else if (procname == "muBrems") return 33;
  else if (procname == "muMinusCaptureAtRest") return 34;
  
  //hadrons 
  else if (procname == "hIoni") return 40;
  else if (procname == "ionIoni") return 41;
  else if (procname == "hadElastic") return 42;
  else if (procname == "hBertiniCaptureAtRest") return 43;
  else if (procname == "hFritiofCaptureAtRest") return 44;
  else if (procname == "hBertiniCaptureAtRest") return 45;  
  else if (procname == "hBrems") return 46;
  else if (procname == "hPairProd") return 47;
  
  //inelastic
  else if (procname == "protonInelastic") return 50;
  else if (procname == "anti_protonInelastic") return 51;
  
  else if (procname == "pi-Inelastic") return 52;
  else if (procname == "pi+Inelastic") return 53;

  else if (procname == "kaon0SInelastic") return 54;
  else if (procname == "kaon0LInelastic") return 55;

  else if (procname == "kaon+Inelastic") return 56;
  else if (procname == "kaon-Inelastic") return 57;

  else if (procname == "tInelastic") return 58; 
  else if (procname == "anti_tritonInelastic") return 59;
  
  else if (procname == "dInelastic") return 60;
  else if (procname == "anti_deuteronInelastic") return 61;
  
  else if (procname == "sigma+Inelastic") return 62;
  else if (procname == "anti_sigma+Inelastic") return 63;

  else if (procname == "sigma-Inelastic") return 64;
  else if (procname == "anti_sigma-Inelastic") return 65;

  else if (procname == "xi-Inelastic") return 66;
  else if (procname == "anti_xi-Inelastic") return 67;

  else if (procname == "xi0Inelastic") return 68;
  else if (procname == "anti_xi0Inelastic") return 69;

  else if (procname == "omega-Inelastic") return 70;
  else if (procname == "anti_omega-Inelastic") return 71;  

  else if (procname == "lambdaInelastic") return 72;
  else if (procname == "anti-lambdaInelastic") return 73;  

  else if (procname == "alphaInelastic") return 74;
  else if (procname == "anti_alphaInelastic") return 75;
  
  else if (procname == "He3Inelastic") return 76;
  else if (procname == "anti_He3Inelastic") return 77;

  else if (procname == "ionInelastic") return 78;

  //others
  else if (procname == "CoulombScat") return 80;
  else if (procname == "msc") return 81;
  
  //radioactive 
  else if (procname == "RadioactiveDecay") return 90;
  else if (procname == "Decay") return 91;

  else if (procname == "Transportation") return 100;
  else if (procname == "unknown") return 900;
  else return 1000;
				 
}



// Method to check if given point is within the given volume
//G4bool IsPointInsideVolume(G4ThreeVector& point, G4String VolName)
//{
//  G4ThreeVector null(0.,0.,0.);
//  G4ThreeVector *ptr;
//  ptr = &null;

//  // Check position is within VolName, if so true, else false
//  G4VPhysicalVolume *theVolume;
//  G4Navigator *gNavigator =G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
//  theVolume=gNavigator->LocateGlobalPointAndSetup(point,ptr,true);
//  if(!theVolume) return(false);
//  G4String theVolName = theVolume->GetName();
//  if(theVolName == VolName)
//    {
//      return(true);
//    }
//  else
//    return(false);
//}

// Method to get the volume name to which a given point belongs
G4String IsPointInsideVolume(G4ThreeVector& point)
{
  G4ThreeVector null(0.,0.,0.);
  G4ThreeVector *ptr;
  ptr = &null;

  // Check position is within VolName, if so true, else false
  G4VPhysicalVolume *theVolume;
  G4Navigator *gNavigator =G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
  theVolume=gNavigator->LocateGlobalPointAndSetup(point,ptr,true);
  if(!theVolume){
      G4cout << "Particle is not in any volume" << G4endl;
      return " ";
  }else{
      return theVolume->GetName();
  }
}


// Method to sample the IBD volume selection
//G4int SelectIBDVol(std::vector<double> choice_weight)
//{
//    G4int num_choices = choice_weight.size();
//    G4double sum_of_weight = 0.;

//    for(int i=0; i<num_choices; i++) {
//        sum_of_weight += choice_weight[i];
//    }
//    G4double rnd = G4UniformRand()*sum_of_weight;
////    int rnd = random(sum_of_weight);
//    for(int i=0; i<num_choices; i++) {
//        if(rnd < choice_weight[i])
//            return i;
//        rnd -= choice_weight[i];
//    }

//    G4cout << "--> no choice selected !!!" << G4endl;
//    return -1;
//}


//G4VPhysicalVolume* GetPhysVolume(G4String VolName)
//{
//  G4VPhysicalVolume *tempPV      = NULL;
//  G4PhysicalVolumeStore *PVStore = 0;
//  G4String theRequiredVolumeName = VolName;
//  PVStore      = G4PhysicalVolumeStore::GetInstance();
//  G4int      i = 0;
//  G4bool found = false;
//  while (!found && i<G4int(PVStore->size())) {
//    tempPV = (*PVStore)[i];
//    found  = tempPV->GetName() == theRequiredVolumeName;
//    if (!found)
//      {i++;}
//  }
//  // found = true then the volume exists else it doesnt.
//  if(found == true)
//    {
//      return tempPV;
//    }
//  else
//    {
//      G4cout << " **** Error: Volume does not exist **** " << G4endl;
//      return NULL;
//    }
//}



//====================================================
// List of all logVolumes in th geometry
// Mass computation of each selected volume
//====================================================
void ListOfLogVolumes()
{
    G4cout << G4endl << "**************Volumes in the geometry****************" << G4endl;

    SolidMassComputation * MassComputation = new SolidMassComputation(1);

    G4VPhysicalVolume *tempPV      = NULL;
    G4PhysicalVolumeStore *PVStore = 0;
    PVStore = G4PhysicalVolumeStore::GetInstance();

    for (G4int i = 0 ; i < G4int(PVStore->size()) ; i++) {
        tempPV = (*PVStore)[i];
        G4LogicalVolume * logToCheck = tempPV->GetLogicalVolume();
        if (logToCheck == NULL) {
            G4cerr << "ERROR !! : " << tempPV->GetName() << " doesn't exist. Mass check aborted." << G4endl;
        }
        else {
            G4cout << tempPV->GetName() << G4endl;
            MassComputation->ComputeMassForLogical(logToCheck);
        }
    }
    G4cout << G4endl << "**************END Volumes in the geometry****************" << G4endl << G4endl;

}


// check if the volume is insed the full detector volume
// (everithing inside the Modules array)
bool IsVolInsideDetector(G4String volName)
{

    if ( IsVolInsideCubeMod(volName)  // cubes
         || IsVolInsideFibre(volName) // fibres
         /// other parts attached to fibres
         || volName=="physFibreGrease"
         || volName=="physFibreMirror"
         || volName=="physFibreMPPC"
         /// modules
         || volName=="physDetector"
         || volName=="physModule"
         || volName=="physModuleSupport"
         || volName=="physTopPost"
         || volName=="physVerPost"
         || volName=="physBottomPostX"
         || volName=="physBottomPostZin"
         || volName=="physBottomPostZout"
         || volName=="physSupportFoot"
         /// Al Frames
         || volName=="physAlFrame"
         || volName=="physFrameAlBars"
         || volName=="physFrameAlCovers"
         /// Full Plane
         || volName=="physPlane"
         || volName=="physPlanePEReflector"
         || volName=="physPlaneCoverFront"
         || volName=="physPlaneCoverBack"
         /// Al electronic box
         || volName=="physElectronicBox"
         || volName=="physEBSupport"
         || volName=="physEBSideVertExt"
         || volName=="physEBSideVertInt"
         || volName=="physEBHorSideBox"
         || volName=="physEBCover"
         )

        return true;
    else
        return false;

}


// check if the volume is insed the DetMod (PVT + LiF:ZnS + Tyvek + air gaps)
bool IsVolInsideCubeMod(G4String volName)
{

    if ( volName=="physDetMod"
         || volName=="physPVT"
         || volName=="physLi"
         || volName=="physBacking"
         || volName=="physTyvek"
         )

        return true;
    else
        return false;

}


// check if the volume is insed the fibres
bool IsVolInsideFibre(G4String volName)
{

    if ( volName=="physFibreCore"
         || volName=="physFibreClad"
         || volName=="physFibreClad2"
//         || volName=="physFibreGrease"
//         || volName=="physFibreMirror"
//         || volName=="physFibreMPPC"
         )

        return true;
    else
        return false;

}



int GetVolID(G4String volName, G4int LiID)
{

    std::map<G4String, int> volNamesID = {
        {	"physWorld"                 ,	1	}	,
        {	"physLab"                   ,	2	}	,
        /// DetMod (PVT + Tyvek + ZnS)
        {	"physDetMod"               	,	10	}	,
        {	"physPVT"               	,	15	}	,
        {	"physLi"                    ,	20	}	,
        {	"physBacking"               ,	30	}	,
        {	"physTyvek"             	,	40	}	,
        /// Fibers
        {	"physFibreCore"             ,	50	}	,
        {	"physFibreClad"             ,	51	}	,
        {	"physFibreClad2"            ,	52	}	,
        {	"physFibreGrease"           ,	53	}	,
        {	"physFibreMirror"           ,	54	}	,
        {	"physFibreMPPC"             ,	55	}	,
        /// Al Frames
        {	"physAlFrame"               ,	61	}	,
        {	"physFrameAlBars"           ,	61	}	,
        {	"physFrameAlCovers"         ,	61	}	,
        /// Full Plane
        {	"physPlane"                 ,	60	}	,
        {	"physPlanePEReflector"      ,	62	}	,
        {	"physPlaneCoverFront"       ,	63	}	,
        {	"physPlaneCoverBack"        ,	63	}	,
        /// Al electronic box
        {	"physElectronicBox"     	,	65	}	,
        {	"physEBSupport"             ,	65	}	,
        {	"physEBSideVertExt"         ,	65	}	,
        {	"physEBSideVertInt"     	,	65	}	,
        {	"physEBHorSideBox"          ,	65	}	,
        {	"physEBCover"           	,	65	}	,
        /// module
        {	"physModule"                ,	70	}	,
        {	"physModuleSupport"         ,	71	}	,
        {	"physTopPost"               ,	72	}	,
        {	"physVerPost"               ,	72	}	,
        {	"physBottomPostX"           ,	72	}	,
        {	"physBottomPostZin"         ,	72	}	,
        {	"physBottomPostZout"        ,	72	}	,
        {	"physSupportFoot"           ,	72	}	,
        /// detector virtual volume
        {	"physDetector"              ,	75	}	,
        /// PE neutron reflectors
        {	"physNeutronReflectorBack"	,	80	}	,
        {	"physNeutronReflectorFront"	,	81	}	,
        /// Support rail system
        {	"physSupportRails"          ,	85	}	,
        {	"physSupportRailsZ"         ,	86	}	,
        {	"physSupportRailsX"         ,	86	}	,
        {	"physSupportCover"          ,	86	}	,
        {	"physSupportTube"           ,	86	}	,
        {	"physSupportTubeBase"       ,	86	}	,
        /// Container
        {	"physContainer"                     ,	90	}	,
        {	"physContainerISOFrame"             ,	91	}	,
        {	"physISOFrameBackPost"              ,	91	}	,
        {	"physISOFrameFrontPost"             ,	91	}	,
        {	"physISOFrameTopBackPost"           ,	91	}	,
        {	"physISOFrameBottomFrontBackPost"   ,	91	}	,
        {	"physISOFrameBackWall"              ,	91	}	,
        {	"physISOFrameTopFrontPost"          ,	91	}	,
        {	"physISOFrameFrontWall"             ,	91	}	,
        {	"physISOFrameSideStruct"            ,	91	}	,
        {	"physISOFrameSideWall"              ,	91	}	,
        {	"physISOFrameFloorBar"              ,	91	}	,
        {	"physISOFrameFloorPanel"            ,	91	}	,
        {	"physISOFrameRoofPanel"             ,	91	}	,
        {	"physContainerInsulationAlcover"    ,	92	}	,
        {	"physContainerInsulation"           ,	93	}	,
        {	"physContainerAlFoil"               ,	94	}	,
        {	"physContainerTopPanel"             ,	95	}	,
        {	"physContainerBottomPanel"          ,	96	}	,
        {	"physContainerFoot"                 ,	97	}	,
        /// Shielding
        {	"physShielding"                 ,	100	}	,
        {	"physShieldingExternalSupport"  ,	101	}	,
        {	"physShieldingInternalSupport"  ,	102	}	,
        {	"physShieldingWaterTank"        ,	103	}	,
        {	"physShieldingWater"            ,	104	}	,
        {	"physShieldingTopPlane"         ,	105	}	,
        {	"physShieldingBottomPlane"      ,	106	}	,
        {	"physShieldingBottomCd"         ,	107	}	,
        {	"physShieldingTopCd"            ,	108	}	,
        {	"physShieldingBackCd"           ,	109	}	,
        {	"physShieldingTopSupport"       ,	110	}	,
        {	"physShieldingSteelGutter"      ,	111	}	,
        /// BR2
        {	"physBR2CB"                 ,	200	}	,
        {	"physBR2CBDome"             ,	201	}	,
        {	"physBR2ExternalWall"       ,	202	}	,
        {	"physBR2InternalWall"     	,	203	}	,
        {	"physBR2ConcreteWall"       ,	204	}	,
        {	"physBR2Level2Floor"        ,	205	}	,
        {	"physBR2Level3Floor"        ,	206	}	,
        {	"physBR2Level3FloorEnforced",	206	}	,
        {	"physBR2Level4Floor"        ,	207	}	,
        {	"physBR2Level5Floor"        ,	208	}	,
        {	"physBR2Level6Floor"        ,	209	}	,
        {	"physBR2Level7Floor"        ,	210	}	,
        {	"physBR2Level8Floor"        ,	211	}	,
        {	"physBR2PillarC9"        	,	212	}	,
        {	"physBR2ReactorCore"       	,	213	}	,
        {	"physBR2WaterPool"          ,	214	}	,
        {	"physBR2LeadWallR1"         ,	215	}	,
        {	"physBR2LeadWallT2"         ,	216	}	,
        {	"physBR2LeadWallT3"         ,	217	}	,
        {	"physBR2AlWall"             ,	218	}	,
        {	"physBR2ParaffinBlock"      ,	219	}	,
        /// CROSS structures
        {	"physCrossPiece0"   ,	230	}	,
        {	"physCrossPiece1"   ,	230	}	,
        {	"physCrossPiece2"   ,	230	}	,
        {	"physCrossPiece3"   ,	230	}	,
        {	"physCrossPiece4"   ,	230	}	,
        {	"physCrossPiece5"   ,	230	}	,
        {	"physCrossPiece6"   ,	230	}	,
        {	"physCrossPiece7"   ,	230	}	,
        {	"physCrossPiece8"   ,	230	}	,
        {	"physCrossPiece9"   ,	230	}	,
        {	"physCrossPiece10"  ,	230	}	,
        {	"physCrossGutter"   ,	230	}	,
        /// CROSS Carrier structures
        {	"physCrossCarrier"                      ,	231	}	,
        {	"physSideTube1CrossCarrier"             ,	231	}	,
        {	"physSideTube2CrossCarrier"             ,	231	}	,
        {	"physCentralRackCrossCarrier"           ,	231	}	,
        {	"physSupportTube1CrossCarrier"          ,	231	}	,
        {	"physSupportTube2CrossCarrier"          ,	231	}	,
        {	"physSupportTube3CrossCarrier"          ,	231	}	,
        {	"physSupportTube4CrossCarrier"          ,	231	}	,
        {	"physSupportBarCrossCarrier"            ,	231	}	,
        {	"physSupportRackCrossCarrier"           ,	231	}	,
        {	"physCentralFinalSupportCrossCarrier"   ,	231	}	,
        {	"physSideSupportCrossCarrier"           ,	231	}	,
        {	"physSide2SupportCrossCarrier"          ,	231	}	,
        {	"physCrossPiece4"                       ,	231	}	,
        {	"physCrossPiece5"                       ,	231	}	,
        {	"physCrossPiece6"                       ,	231	}	,
        {	"physCrossCarrierHead"                  ,	232	}	,
        {	"physCrossCarrierHeadSource"            ,	232	}	,
        {	"physCrossSourceCapsule"                ,	233	}	,
    };

    if( !volNamesID[volName] )
        return -999999;
    else
        return volNamesID[volName] + LiID;

}


double CubesDistance(int cube1, int cube2)
{
    int cube1z = cube1/10000;
    int cube1x = (cube1%10000)/100;
    int cube1y = cube1%100;
    int cube2z = cube2/10000;
    int cube2x = (cube2%10000)/100;
    int cube2y = cube2%100;
    int distx = abs(cube2x - cube1x);
    int disty = abs(cube2y - cube1y);
    int distz = abs(cube2z - cube1z);

    return sqrt(distx*distx + disty*disty + distz*distz);

}


// Extruded piece (used for Modules support in Phase1)
// type parameter indicates the number of repetition in X of the basic unit
G4LogicalVolume* AlExtrusionConstruction(G4int type, G4double length, G4double width, G4VisAttributes* vis)
{

    G4Material* mate;

    G4double holeThickness = 2.*mm;


    // construct log volume box
    G4Box* volBox = new G4Box("volBox", 0.5*width*type, 0.5*width, 0.5*length);

    mate = G4Material::GetMaterial("Air");
    G4LogicalVolume* logVolPiece = new G4LogicalVolume(volBox, mate, "logVolPiece", 0, 0, 0);
    logVolPiece->SetVisAttributes(G4VisAttributes::Invisible);


    // construct piece box
    G4Box* pieceBox = new G4Box("pieceBox", 0.5*width, 0.5*width, 0.5*length);

    // construct hole in the center
    G4double tubeDiam = 7.*mm;
    G4Tubs* tube = new G4Tubs("tube", 0, 0.5*tubeDiam, 0.5*length+1.*mm, 0, 360.*deg);
    G4SubtractionSolid* piece = new G4SubtractionSolid("piece", pieceBox, tube, 0, G4ThreeVector(0,0,0));

    // construct square holes in the corners
    G4double hole1 = 12.*mm;
    G4double pos = 0.5*width - holeThickness - 0.5*hole1;
    G4Box* hole1Box = new G4Box("hole1Box", 0.5*hole1, 0.5*hole1, 0.5*length+1.*mm);
    piece = new G4SubtractionSolid("piece", piece, hole1Box, 0, G4ThreeVector( pos,  pos, 0));
    piece = new G4SubtractionSolid("piece", piece, hole1Box, 0, G4ThreeVector( pos, -pos, 0));
    piece = new G4SubtractionSolid("piece", piece, hole1Box, 0, G4ThreeVector(-pos,  pos, 0));
    piece = new G4SubtractionSolid("piece", piece, hole1Box, 0, G4ThreeVector(-pos, -pos, 0));

    // construct holes in the sides
    // as intersection of triangle & rectangle
    // and the result in union with another small rectangle in the bottom
/*
           /\
          /	 \
         /	  \
        /	   \		________		 _______				 _______
       /		\	∩  |   		|  =	/       \	U		=	/	    \
      /			 \	   |		|		|       |				|	    |
     /	Hole21	  \	   |Hole22	|		|Hole2  |	  Hole23	|Hole3  |
    /______________\   |________|		|_______|	   ____		|_______|
                                                       |___|	  |___|
*/

    G4double hole21Base = 32.*mm;
    G4double hole21Height = 16.*mm;
    G4double hole22Width = 13.*mm;
    G4double hole22Height = 12.*mm;
    G4double hole23Width = 7.*mm;
    G4double hole23Height = 4.5*mm;
    G4Trd* hole21Box = new G4Trd("hole21Box", 0.5*hole21Base, 0, 0.5*length+1.*mm, 0.5*length+1.*mm, 0.5*hole21Height);
    G4Box* hole22Box = new G4Box("hole22Box", 0.5*hole22Width, 0.5*hole22Height, 0.5*length+1.*mm);
    G4Box* hole23Box = new G4Box("hole23Box", 0.5*hole23Width, 0.5*hole23Height+1.*mm, 0.5*length+1.*mm);
    G4RotationMatrix* rot = new G4RotationMatrix();
    rot->rotateX(90*deg);
    G4IntersectionSolid* hole2 = new G4IntersectionSolid("hole2", hole22Box, hole21Box, rot, G4ThreeVector(0, 0.5*(hole21Height - hole22Height), 0));
    G4UnionSolid* hole3 = new G4UnionSolid("hole3", hole2, hole23Box, 0, G4ThreeVector(0, -0.5*(hole22Height + hole23Height) + 0.5*mm, 0));
    pos = 0.5*width - 0.5*(hole22Height + hole23Height) - 2.*mm;
    piece = new G4SubtractionSolid("piece", piece, hole3, 0, G4ThreeVector(0, -pos, 0));
    rot = new G4RotationMatrix();
    rot->rotateZ(180*deg);
    piece = new G4SubtractionSolid("piece", piece, hole3, rot, G4ThreeVector(0, pos, 0));
    rot = new G4RotationMatrix();
    rot->rotateZ(90*deg);
    piece = new G4SubtractionSolid("piece", piece, hole3, rot, G4ThreeVector(-pos, 0, 0));
    rot = new G4RotationMatrix();
    rot->rotateZ(270*deg);
    piece = new G4SubtractionSolid("piece", piece, hole3, rot, G4ThreeVector(pos, 0, 0));

    mate = G4Material::GetMaterial("Aluminium_6061");
    G4LogicalVolume* logPiece = new G4LogicalVolume(piece, mate, "logPiece", 0, 0, 0);
    logPiece->SetVisAttributes(vis);

    for(int i=0; i<type; i++)
    {
        new G4PVPlacement(0, G4ThreeVector(0.5*width*(1 - type + 2.*i), 0., 0.),
                          logPiece, "physPiece", logVolPiece, false, 0);
    }

    return logVolPiece;

}


// alternative for the extruded piece above
// constructed just as hollow bar
// type parameter indicates the number of repetition in X of the basic unit
/*
 *      _________
 *     |  _____  |
 *     | |     | |
 *     | |     | |
 *     | |_____| |
 *     |_________|
*/
G4LogicalVolume* AlternativeExtrusionConstruction(G4int type, G4double length, G4double width, G4VisAttributes* vis)
{

    G4Material* mate;

    G4double thickness = 4.25*mm;

    // construct log volume box
    G4Box* volBox = new G4Box("volBox", 0.5*width*type, 0.5*width, 0.5*length);
    mate = G4Material::GetMaterial("Air");
    G4LogicalVolume* logVolPiece = new G4LogicalVolume(volBox, mate, "logVolPiece", 0, 0, 0);
    logVolPiece->SetVisAttributes(G4VisAttributes::Invisible);

    // construct piece box
    G4Box* pieceOuterBox = new G4Box("pieceOuterBox", 0.5*width, 0.5*width, 0.5*length);
    G4Box* pieceInnerBox = new G4Box("pieceInnerBox", 0.5*width - thickness, 0.5*width - thickness, 0.5*length+0.5*mm);
    // make hole in the center
    G4SubtractionSolid* piece = new G4SubtractionSolid("piece", pieceOuterBox, pieceInnerBox, 0, G4ThreeVector(0, 0, 0));

    mate = G4Material::GetMaterial("Aluminium_6061");
    G4LogicalVolume* logPiece = new G4LogicalVolume(piece, mate, "logPiece", 0, 0, 0);
    logPiece->SetVisAttributes(vis);

    for(int i=0; i<type; i++)
    {
        new G4PVPlacement(0, G4ThreeVector(0.5*width*(1 - type + 2.*i), 0., 0.),
                          logPiece, "physPiece", logVolPiece, false, 0);
    }

    return logVolPiece;

}