/**
 * @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)
{

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



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;

}