LBNFTargetModule.cc

Go to the documentation of this file.

// Code that defines the volumes and placements of the target module (e.g. SAT)
// which includes some engineering details. John Back, March 2016.

#include "LBNEVolumePlacements.hh"

#include "G4Material.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4RotationMatrix.hh"
#include "G4Sphere.hh"
#include "G4ThreeVector.hh"
#include "G4Tubs.hh"
#include "G4VPhysicalVolume.hh"
#include "G4UnionSolid.hh"

#include "LBNERunManager.hh"
#include "LBNEPrimaryGeneratorAction.hh"

#include <iostream>
#include <sstream>

void LBNEVolumePlacements::InitTargetModule() {

    std::cout << "Initialising the target module parameters for type = " << fTargetModuleType << std::endl;
    std::cout << "Target radius = " << fTargetRadius << std::endl;
    std::cout << "Target total original length = " << fTargetLength << std::endl;

    // Change the target object length if we have set the required number of interaction lengths
    bool useLambda(false);

    if (fTargetNLambda > -1.0 && fTargetDensity > 0.0) {

        std::cout << "Updating target object length using interaction length info:" << std::endl;
        fTargetLength = fTargetNLambda*fTargetSpecificLambda/fTargetDensity;
        std::cout << "Specific lambda = " << fTargetSpecificLambda*CLHEP::cm2/CLHEP::g
                  << ", rho = " << fTargetDensity*CLHEP::cm3/CLHEP::g << ", NLambda = " 
                  << fTargetNLambda << " => target L = " << fTargetLength << std::endl;

        useLambda = true;

    }

    fNumTargetObjects = 1;

    if (fTargetModuleType == LBNEVolumePlacements::SAT) {

        std::cout << "Checking length for SAT" << std::endl;

        double tolerance = 0.005*CLHEP::mm; // 5 micron tolerance
        double diameter = fTargetRadius*2.0;
        double NumNominal = fTargetLength/(diameter + tolerance);

        // If we have specified the number of interaction lengths, then
        // find the required number of spheres using the full calculation
        // given in LBNE-doc-9547 (Quynh) assuming a Gaussian beam
        if (useLambda) {NumNominal = this->FindNTargetSpheres(1);}
        std::cout << "Nominal number of spheres = " << NumNominal << std::endl;

        // Round to nearest integer
        fNumTargetObjects = (int) (NumNominal + 0.5);
        std::cout << "Number of spheres rounded to "<< fNumTargetObjects <<std::endl;

        // Update the total "useful" length of the SAT, i.e. the sum of the length of the spheres
        fTargetLength = fNumTargetObjects*(diameter + tolerance);
        std::cout << "SAT total length adjusted to " << fTargetLength << std::endl; 

    }

    // The diameter thickness of the supporting rods, which also leaves a gap for He gas flow
    // in the remaining azimuthal areas. For the cylinder target, we have just the gap, no rods
    fTargetSupportD = 2.5*CLHEP::mm;
    std::cout << "Target module support radial thickness = " << fTargetSupportD << std::endl;

    // Thickness of the titanium tube casings
    fTargetCaseT = 1.0*CLHEP::mm;
    std::cout << "Target module Ti casing thickness = " << fTargetCaseT << std::endl;

    // The He gas gap between the two inner and outer Ti cylindrical casings
    fTargetHeGap = 4.0*CLHEP::mm;
    std::cout << "Target module He gap between casings = " << fTargetHeGap << std::endl;

    // Inner radius of the outer casing
    fTargetOutCaseInnerR = fTargetRadius + fTargetSupportD + fTargetCaseT + fTargetHeGap;
    std::cout << "Target module inner radius of outer Ti casing = " << fTargetOutCaseInnerR << std::endl;
                        
    // The length of the inner canister holding the spheres
    fTargetInCaseUpL = 50.0*CLHEP::mm; // Could also use fTargetRadius if we want to scale the size to rSphere?
    fTargetInCaseDnL = 5.0*CLHEP::mm;
    fTargetInCaseL   = fTargetLength + fTargetInCaseUpL + fTargetInCaseDnL;

    // The length of the outer target module canister
    fTargetCaseDiffL = 5.0*CLHEP::mm;
    fTargetOutCaseL  = fTargetInCaseL + fTargetCaseDiffL;

    // Total length of the target module, including the end bulb of the outer canister
    fTargetModuleTotL = fTargetOutCaseL + fTargetOutCaseInnerR + fTargetCaseT;
    std::cout << "Target module total length = " << fTargetModuleTotL << std::endl;

    // Also update the target length outside the horn if we have specified the fraction.
    // By default, this fraction is not set (=-1), with fTargetLengthOutsideHorn = 450.0 mm.
    // The length includes the outer canister end regions
    if (fTargetFracOutHornL > -1.0) {

        fTargetLengthOutsideHorn = fTargetFracOutHornL*fTargetModuleTotL;
    }

    std::cout << "TargetFracOutHornL = " << fTargetFracOutHornL
              << ", TargetLengthOutsideHorn = " << fTargetLengthOutsideHorn << std::endl;
    std::cout << "fTargetLengthOutsideExtra = " << fTargetLengthOutsideExtra << std::endl;

}

void LBNEVolumePlacements::InitTarget2Module() {

    std::cout << "Initialising the 2nd target module parameters for type = " << fTarget2ModuleType << std::endl;
    std::cout << "Target2 radius = " << fTarget2Radius << std::endl;
    std::cout << "Target2 total original length = " << fTarget2Length << std::endl;

    // Change the target object length if we have set the required number of interaction lengths
    bool useLambda(false);

    if (fTarget2NLambda > -1.0 && fTarget2Density > 0.0) {

        std::cout << "Updating target 2 object length using interaction length info:" << std::endl;
        fTarget2Length = fTarget2NLambda*fTarget2SpecificLambda/fTarget2Density;
        std::cout << "Specific lambda = " << fTarget2SpecificLambda*CLHEP::cm2/CLHEP::g
                  << ", rho = " << fTarget2Density*CLHEP::cm3/CLHEP::g << ", NLambda = " 
                  << fTarget2NLambda << " => target2 L = " << fTarget2Length << std::endl;

        useLambda = true;

    }

    fNumTarget2Objects = 1;

    if (fTarget2ModuleType == LBNEVolumePlacements::SAT) {

        std::cout << "Checking length for 2nd SAT" << std::endl;

        double tolerance = 0.005*CLHEP::mm; // 5 micron tolerance
        double diameter = fTarget2Radius*2.0;
        double NumNominal = fTarget2Length/(diameter + tolerance);

        // If we have specified the number of interaction lengths, then
        // find the required number of spheres using the full calculation
        // given in LBNE-doc-9547 (Quynh) assuming a Gaussian beam
        if (useLambda) {NumNominal = this->FindNTargetSpheres(2);}
        std::cout << "Nominal number of spheres = " << NumNominal << std::endl;

        // Round to nearest integer
        fNumTarget2Objects = (int) (NumNominal + 0.5);
        std::cout << "Number of spheres rounded to "<< fNumTarget2Objects <<std::endl;

        // Update the total "useful" length of the SAT, i.e. the sum of the length of the spheres
        fTarget2Length = fNumTarget2Objects*(diameter + tolerance);
        std::cout << "SAT2 total length adjusted to " << fTarget2Length << std::endl; 

    }

    // Inner radius of the outer casing
    fTarget2OutCaseInnerR = fTarget2Radius + fTargetSupportD + fTargetCaseT + fTargetHeGap;
    std::cout << "Target module 2 inner radius of outer Ti casing = " << fTarget2OutCaseInnerR << std::endl;
                        
    // The length of the inner canister holding the spheres
    fTarget2InCaseL    = fTarget2Length + fTargetInCaseUpL + fTargetInCaseDnL;

    // The length of the outer target module canister
    fTarget2OutCaseL   = fTarget2InCaseL + fTargetCaseDiffL;

    // Total length of the target module, including the end bulb of the outer canister
    fTarget2ModuleTotL = fTarget2OutCaseL + fTarget2OutCaseInnerR + fTargetCaseT;
    std::cout << "Target module 2 total length = " << fTarget2ModuleTotL << std::endl;

    // Mirror image 180 degree rotation matrix: rotate around the y axis
    fMirrorRotation = new G4RotationMatrix();
    fMirrorRotation->rotateY(180.0*CLHEP::degree);

    // z displacement for the 2nd target module: depends on length of 1st module
    fZModuleShift = fTargetModuleTotL + fTarget2OutCaseInnerR + fTargetCaseT + 1.0*CLHEP::mm;

}

G4double LBNEVolumePlacements::FindNTargetSpheres(int number) const {

    // Find the length of the target (in mm) that gives the required number of interaction lengths
    // for the specified radius, assuming a Gaussian beam distribution based on the calculation
    // given in LBNE-doc-9547 (Quynh).
    // Performs the double integral of the effective interaction length for sigma = radius/3.
    // Note that the beam parameters are defined after the geometry, so the user needs to make
    // sure the beam sigmas are set to radius/3...
 
    double targetRadius = fTargetRadius;
    double targetNLambda = fTargetNLambda;
    double targetSpecificLambda = fTargetSpecificLambda;
    double targetDensity = fTargetDensity;

    if (number == 2) {
        targetRadius = fTarget2Radius;
        targetNLambda = fTarget2NLambda;
        targetSpecificLambda = fTarget2SpecificLambda;
        targetDensity = fTarget2Density;
    }

    std::cout << "Finding number of target spheres for NLambda = " << targetNLambda << std::endl;

    double xSigma = targetRadius/3.0;
    double ySigma = targetRadius/3.0;
    double x0(0.0), y0(0.0);
    
    /*
    // This part is what we could do if the generator action is available during geometry construction:
    LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
    if (theRunManager) {
        
        const LBNEPrimaryGeneratorAction* PGA = static_cast<const LBNEPrimaryGeneratorAction*> (theRunManager->GetUserPrimaryGeneratorAction());
        
        if (PGA) {
        
            xSigma = PGA->GetBeamSigmaX();
            ySigma = PGA->GetBeamSigmaY();
            x0 = PGA->GetBeamOffsetX();
            y0 = PGA->GetBeamOffsetY();
    
        } else {

            std::cout << "Primary generator pointer not found" << std::endl;
        }

    }
    */
    
    std::cout << "Assuming beam sigma_x = " << xSigma << ", sigma_y = " << ySigma << std::endl;
    std::cout << "Assuming beam offset_x = " << x0 << ", offset_y = " << y0 << std::endl;

    int N(1000);
    int N1(N+1);
    double RSq = targetRadius*targetRadius;

    double xMin(-targetRadius);
    double xMax(targetRadius);
    double dx = (xMax - xMin)/(N*1.0);
    double dy(dx);
    double yMin(-targetRadius);

    double invXSigmaSq(0.0);
    if (xSigma > 0.0) {invXSigmaSq = 1.0/(xSigma*xSigma);}
    double invYSigmaSq(0.0);
    if (ySigma > 0.0) {invYSigmaSq = 1.0/(ySigma*ySigma);}

    double norm(1.0);
    if (xSigma > 0.0 && ySigma > 0.0) {
        norm = 1.0/(M_PI*xSigma*ySigma);
    }

    double result(0.0);

    int i(0), j(0);
    for (i = 0; i < N1; i++) {

        double x = dx*i + xMin - x0;
        double xSq = x*x;

        for (j = 0; j < N1; j++) {

            double y = dy*j + yMin - y0;
            double ySq = y*y;

            double rSq = xSq + ySq;

            if (rSq < RSq) {

                double expPower = -0.5*(xSq*invXSigmaSq + ySq*invYSigmaSq);
                double expTerm(0.0);
                if (fabs(expPower) < 30.0) {expTerm = exp(expPower);}

                result += expTerm*sqrt(RSq - rSq);

            }

        } // j

    } // i

    // Normalisation
    result *= norm*dx*dy;
    std::cout << "Sphere double integral = " << result << std::endl;

    G4double nSpheres = targetNLambda*targetSpecificLambda/(targetDensity*result);
    return nSpheres;

}

void LBNEVolumePlacements::PlaceTargetModule(int number) {

    std::cout << "Calling PlaceTargetModule" << number << std::endl;
    // InitTargetModule should already be called by LBNEVolumePlacements::setEntireTargetDims()

    this->PlaceTargetOuterCan(number);
    this->PlaceTargetInnerCan(number);
    this->PlaceTargetSupport(number);
    this->PlaceTargetObject(number);

}

void LBNEVolumePlacements::PlaceTargetOuterCan(int number) {

    // Similar to TargetNoSplitM1 (and also TargetNoSplitHeContainer)
    std::ostringstream cNumStrStr; cNumStrStr << number;
    G4String numString = cNumStrStr.str();
   
    // The target hall mother volume
    G4String outerCanName("Target");
    outerCanName += numString; outerCanName += "OuterCan";

    const LBNEVolumePlacementData *mother = Find(outerCanName.data(), "TargetHallAndHorn1", "");
    if (!mother) {return;}
    double hallLength = mother->fParams[2];

    // The horn1 volume
    const LBNEVolumePlacementData* horn1 = Find(outerCanName.data(), "Horn1PolyM1", "");
    if (!horn1) {return;}

    // The outer canister volume
    LBNEVolumePlacementData* canInfo = this->CreateTargetVol(outerCanName, number);
    
    G4ThreeVector position(0.0, 0.0, 0.0);

    // Correction for the placement of the upstream end section of the horn
    G4double horn1Length = horn1->fParams[2];

    // This positions the downstream end of the target module at the upstream edge of the first horn
    // which may require an additional ad-hoc z-shift correction
    double targetModuleTotL = fTargetModuleTotL;
    if (number == 2) {targetModuleTotL = fTarget2ModuleTotL;}

    position[2] = 0.5*hallLength - horn1Length - 0.5*targetModuleTotL;
    // Take into account how much of the target needs to be inside the horn
    G4double targetInHorn = targetModuleTotL - fTargetLengthOutsideHorn;
    position[2] += targetInHorn;

    // Specify the rotation matrix for what will become the mother volume of the target module.
    // All other internal volumes will respect this rotation
    G4RotationMatrix* rotMatrix = 0;
    if (number == 2) {
        position[2] += fZModuleShift;
        rotMatrix = fMirrorRotation;
    }

    // If using the conceptual design horn set-up LBNFConceptDesignHorns::PlaceFinalLBNFConceptHornA()
    // then we need to shift the target module further upstream relative to the mother volume Horn1PolyM1,
    // owing to the larger gap downstream between the target end and the remaining hall/horn1 space.
    // This positions the start of the target module very close to z = 0 m
    if (fUseLBNFOptimConceptDesignHornA) {position[2] -= 31.0*CLHEP::cm;}

    new G4PVPlacement(rotMatrix, position, canInfo->fCurrent,
                      (outerCanName + "_P"), mother->fCurrent, false, 0, fCheckVolumeOverLapWC);

    // The helium gas inside the outer canister
    G4String heName("Target");
    heName += numString; heName += "OuterHeGas";
    LBNEVolumePlacementData* heInfo = this->CreateTargetVol(heName, number);

    // Simply place this helium region inside its mother volume, TargetOuterCan
    new G4PVPlacement((G4RotationMatrix*) 0, G4ThreeVector(0.0, 0.0, 0.0), heInfo->fCurrent,
                      (heName + "_P"), canInfo->fCurrent, false, 0, fCheckVolumeOverLapWC);
    

}

void LBNEVolumePlacements::PlaceTargetInnerCan(int number) {

    std::ostringstream cNumStrStr; cNumStrStr << number;
    G4String numString = cNumStrStr.str();

    G4String innerCanName("Target");
    innerCanName += numString; innerCanName += "InnerCan";

    G4String outerHeName("Target");
    outerHeName += numString; outerHeName += "OuterHeGas";

    const LBNEVolumePlacementData* mother = Find(innerCanName, outerHeName, "");
    if (!mother) {return;}

    // The inner canister volume
    LBNEVolumePlacementData* canInfo = this->CreateTargetVol(innerCanName, number);

    double zShift = -0.5*fTargetCaseDiffL;
    G4ThreeVector position(0.0, 0.0, zShift);
    
    new G4PVPlacement((G4RotationMatrix*) 0, position, canInfo->fCurrent,
                      (innerCanName + "_P"), mother->fCurrent, false, 0, fCheckVolumeOverLapWC);

    // The inner canister helium volume
    G4String innerHeName("Target");
    innerHeName += numString; innerHeName += "InnerHeGas";
    LBNEVolumePlacementData* heInfo = this->CreateTargetVol(innerHeName, number);

    new G4PVPlacement((G4RotationMatrix*) 0, G4ThreeVector(0.0, 0.0, 0.0), heInfo->fCurrent,
                      (innerHeName + "_P"), canInfo->fCurrent, false, 0, fCheckVolumeOverLapWC);

}

void LBNEVolumePlacements::PlaceTargetSupport(int number) {

    // We do not need the internal support rods for the cylinder, but we still keep
    // the fTargetSupportD gap to allow for the He gas flow
    int targetModuleType = fTargetModuleType;
    if (number == 2) {targetModuleType = fTarget2ModuleType;}
    if (targetModuleType == LBNEVolumePlacements::Cylinder) {return;}

    std::ostringstream cNumStrStr; cNumStrStr << number;
    G4String numString = cNumStrStr.str();

    G4String supportName("Target");
    supportName += numString; supportName += "Support";

    G4String innerHeName("Target");
    innerHeName += numString; innerHeName += "InnerHeGas";

    // The target support rods are placed inside the inner He gas region
    const LBNEVolumePlacementData* mother = Find(supportName, innerHeName, "");
    if (!mother) {return;}

    // The support rods; place 3 of these separated by 120 degrees in the x-y plane
    LBNEVolumePlacementData* rodInfo = this->CreateTargetVol(supportName, number);

    G4double x0 = 0.0;
    G4double y0 = fTargetRadius + 0.5*fTargetSupportD;
    if (number == 2) {y0 = fTarget2Radius + 0.5*fTargetSupportD;}

    int i(0);
    for (i = 0; i < 3; i++) {

        G4double phi = -M_PI*120.0*i/180.0;
        G4double cPhi = cos(phi);
        G4double sPhi = sin(phi);
        G4double x = x0*cPhi - y0*sPhi;
        G4double y = y0*cPhi + x0*sPhi;
        G4double z = 0.0;

        //std::cout<<"Target support rod: phi = "<<phi*180.0/M_PI<<", x = "<<x<<", y = "<<y<<std::endl;
        G4ThreeVector position(x, y, z);
        
        std::ostringstream sNumStrStr; sNumStrStr << i; 
        new G4PVPlacement((G4RotationMatrix*) 0, position, rodInfo->fCurrent,
                          (supportName + "_P" + sNumStrStr.str()),
                          mother->fCurrent, false, 0, fCheckVolumeOverLapWC);

    }

}

void LBNEVolumePlacements::PlaceTargetObject(int number) {

    // Retrieve the mother volume for the spheres, which is the inner He gas regions
    std::ostringstream cNumStrStr; cNumStrStr << number;
    G4String numString = cNumStrStr.str();

    G4String label("Target"); label += numString;
    int targetModuleType = fTargetModuleType;
    if (number == 2) {targetModuleType = fTarget2ModuleType;}

    if (targetModuleType == LBNEVolumePlacements::Cylinder) {
        label += "Cylinder";
    } else {
        label += "Sphere";
    }

    G4String innerHeName("Target");
    innerHeName += numString; innerHeName += "InnerHeGas";
    const LBNEVolumePlacementData *mother = Find(label, innerHeName, "");
    if (!mother) {return;}

    LBNEVolumePlacementData* placeInfo = this->CreateTargetVol(label, number);
    
    G4ThreeVector position(0.0, 0.0, 0.0);

    double targetRadius = fTargetRadius;
    double targetInCaseL = fTargetInCaseL;
    int numTargetObjects = fNumTargetObjects;
    if (number == 2) {
        targetRadius = fTarget2Radius;
        targetInCaseL = fTarget2InCaseL;
        numTargetObjects = fNumTarget2Objects;
    }
    
    // Spherical array target set-up
    if (targetModuleType == LBNEVolumePlacements::SAT) {

        // The start of the first sphere
        double z0 = -0.5*targetInCaseL + fTargetInCaseUpL + targetRadius;
        position[2] = z0;
    
        double tolerance = 0.005*CLHEP::mm;
    
        int iSph(0);
        // Loop over all spheres
        for (iSph = 0; iSph < numTargetObjects; iSph++) {
        
            position[2] = z0 + iSph*(targetRadius*2.0 + tolerance);
        
            std::ostringstream sNumStrStr; sNumStrStr << iSph;
            new G4PVPlacement((G4RotationMatrix*) 0, position, placeInfo->fCurrent, 
                              (label + "_P" + sNumStrStr.str()),
                              mother->fCurrent, false, iSph, fCheckVolumeOverLapWC);
        
        }

    } else if (targetModuleType == LBNEVolumePlacements::Cylinder) {

        // Just place the cylinder in the middle of the inner He gas volume
        new G4PVPlacement((G4RotationMatrix*) 0, position, placeInfo->fCurrent,
                          (label + "_P"), mother->fCurrent, false, 0, fCheckVolumeOverLapWC);

    }

}

LBNEVolumePlacementData* LBNEVolumePlacements::CreateTargetVol(const G4String &volName, int number) {

    // Based on LBNEVolumePlacements::Create() but put here in order to reduce
    // the huge size of that function...

    // Check to see if the volume has already been created
    std::map<G4String, LBNEVolumePlacementData>::const_iterator itDupl = fSubVolumes.find(volName);
    if (itDupl != fSubVolumes.end()) {
        std::ostringstream mStrStr;
        mStrStr << " Volume named " << volName << " Already defined. Fatal ";
        G4String mStr(mStrStr.str());
        G4Exception("LBNEVolumePlacements::Create", " ", FatalErrorInArgument, mStr.c_str()); 
    }

    LBNEVolumePlacementData info;
    info.initialize();

    // Check that the volume name at least contains "Target". Otherwise ignore it
    std::ostringstream cNumStrStr; cNumStrStr << number;
    G4String numString = cNumStrStr.str();
    G4String targetLabel("Target"); targetLabel += numString;

    if (!volName.contains("Target")) {return 0;}

    double targetRadius = fTargetRadius;
    double targetLength = fTargetLength;
    double targetOutCaseInnerR = fTargetOutCaseInnerR;
    double targetOutCaseL = fTargetOutCaseL;
    double targetModuleTotL = fTargetModuleTotL;
    double targetInCaseL = fTargetInCaseL;
    if (number == 2) {
        targetRadius = fTarget2Radius;
        targetLength = fTarget2Length;
        targetOutCaseInnerR = fTarget2OutCaseInnerR;
        targetOutCaseL = fTarget2OutCaseL;
        targetModuleTotL = fTarget2ModuleTotL;
        targetInCaseL = fTarget2InCaseL;
    }

    //G4cout<<"CreateTargetVol: r = "<<targetRadius<<", L = "<<targetLength<<", OCIR = "<<targetOutCaseInnerR
    //      <<", OCL = "<<targetOutCaseL<<", ModTotL = "<<targetModuleTotL<<", ICL = "<<targetInCaseL<<G4endl;

    if (volName.contains("Sphere")) { 

        // One of the target spheres
        info.fParams[0] = 0.0;
        info.fParams[1] = targetRadius;
        info.fParams[2] = 0.0; 
        G4Sphere* aSphere = new G4Sphere(volName,
                                         info.fParams[0], //inner radius = 0
                                         info.fParams[1], //outer radius
                                         0.0, 360.0*CLHEP::degree,
                                         0.0, 180.0*CLHEP::degree); //Phi and Theta angles

        info.fCurrent = new G4LogicalVolume(aSphere, G4Material::GetMaterial(targetLabel), volName);
          
    } else if (volName.contains("Cylinder")) { 

        // The target cylinder object
        info.fParams[0] = 0.0;
        info.fParams[1] = targetRadius;
        info.fParams[2] = targetLength; 
        G4Tubs* aTube = new G4Tubs(volName, 0.0, targetRadius,
                                   0.5*targetLength, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        info.fCurrent = new G4LogicalVolume(aTube, G4Material::GetMaterial(targetLabel), volName);
          
    } else if (volName.contains("OuterCan")) {

        // The outer canister, containing the union of a tube with a downstream spherical "bulb"
        double innerR = 0.0;
        double outerR = targetOutCaseInnerR + fTargetCaseT;
        double halfTubeL = 0.5*targetOutCaseL;

        // The cylindrical tube section
        G4String tubeName(volName); tubeName += "Tube";
        G4Tubs* aTube = new G4Tubs(tubeName, innerR, outerR, 
                                   halfTubeL, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        // The end bulb section
        G4String sphereName(volName); sphereName += "Sphere";
        G4Sphere* aSphere = new G4Sphere(sphereName, innerR, outerR,
                                         0.0, 360.0*CLHEP::degree,
                                         0.0, 180.0*CLHEP::degree); //Phi and Theta angles

        // The union of the two volumes
        G4UnionSolid* outerCan = new G4UnionSolid(volName, aTube, aSphere, 0, G4ThreeVector(0.0, 0.0, halfTubeL));

        info.fCurrent = new G4LogicalVolume(outerCan, G4Material::GetMaterial(std::string("Titanium")), volName);

        // Parameters: inner radius, outer radius, full length
        info.fParams[0] = innerR;
        info.fParams[1] = outerR;
        info.fParams[2] = targetModuleTotL;

    } else if (volName.contains("OuterHeGas")) {

        // The outer canister, containing the union of a tube with a downstream spherical "bulb"
        double innerR = 0.0;
        double outerR = targetOutCaseInnerR;
        double halfTubeL = 0.5*targetOutCaseL;

        // The cylindrical tube section
        G4String tubeName(volName); tubeName += "Tube";
        G4Tubs* aTube = new G4Tubs(tubeName, innerR, outerR, 
                                   halfTubeL, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        // The end bulb section
        G4String sphereName(volName); sphereName += "Sphere";
        G4Sphere* aSphere = new G4Sphere(sphereName, innerR, outerR,
                                         0.0, 360.0*CLHEP::degree,
                                         0.0, 180.0*CLHEP::degree); //Phi and Theta angles

        // The union of the two volumes
        G4UnionSolid* outerCan = new G4UnionSolid(volName, aTube, aSphere, 0, G4ThreeVector(0.0, 0.0, halfTubeL));

        info.fCurrent = new G4LogicalVolume(outerCan, G4Material::GetMaterial(std::string("HeliumTarget")), volName);

        // Parameters: inner radius, outer radius, full length
        info.fParams[0] = innerR;
        info.fParams[1] = outerR;
        info.fParams[2] = targetModuleTotL - fTargetCaseT;

    } else if (volName.contains("InnerCan")) {

        // The inner, cylindrical canister of the target module
        double innerR = 0.0;
        double outerR = targetRadius + fTargetSupportD + fTargetCaseT;
        double halfTubeL = 0.5*targetInCaseL;

        // The cylindrical tube section
        G4String tubeName(volName); tubeName += "Tube";
        G4Tubs* aTube = new G4Tubs(tubeName, innerR, outerR, 
                                   halfTubeL, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        info.fCurrent = new G4LogicalVolume(aTube, G4Material::GetMaterial(std::string("Titanium")), volName);

        // Parameters: inner radius, outer radius, full length
        info.fParams[0] = innerR;
        info.fParams[1] = outerR;
        info.fParams[2] = targetInCaseL;

    } else if (volName.contains("InnerHeGas")) {

        // The inner, cylindrical canister of the target module
        double innerR = 0.0;
        double outerR = targetRadius + fTargetSupportD;
        double halfTubeL = 0.5*targetInCaseL;

        // The cylindrical tube section
        G4String tubeName(volName); tubeName += "Tube";
        G4Tubs* aTube = new G4Tubs(tubeName, innerR, outerR, 
                                   halfTubeL, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        info.fCurrent = new G4LogicalVolume(aTube, G4Material::GetMaterial(std::string("HeliumTarget")), volName);

        // Parameters: inner radius, outer radius, full length
        info.fParams[0] = innerR;
        info.fParams[1] = outerR;
        info.fParams[2] = targetInCaseL;

    } else if (volName.contains("Support")) {

        // The target cylinder object
        double radius = 0.5*fTargetSupportD;
        info.fParams[0] = 0.0;
        info.fParams[1] = radius;
        info.fParams[2] = targetInCaseL; 
        G4Tubs* aTube = new G4Tubs(volName, 0.0, radius,
                                   0.5*targetInCaseL, // half-length
                                   0.0, 360.0*CLHEP::degree); // Phi angle extent

        info.fCurrent = new G4LogicalVolume(aTube, G4Material::GetMaterial(std::string("Titanium")), volName); 

    }

    // Keep track of the volume we've made
    fSubVolumes.insert(std::pair<G4String, LBNEVolumePlacementData>(volName, info));

    return &(fSubVolumes.find(volName)->second);

}