LBNEDownstrVolPlacements.cc

Go to the documentation of this file.

#include "LBNEDetectorConstruction.hh"

#include "G4UIdirectory.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithABool.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithoutParameter.hh"
#include "G4UnitsTable.hh"

#include "G4Material.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Polycone.hh"
#include "G4Trap.hh"
#include "G4Cons.hh"
#include "G4Torus.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4SubtractionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4VisAttributes.hh"
#include "globals.hh"
#include "G4Transform3D.hh"
#include "G4RotationMatrix.hh"
#include "G4PVReplica.hh"
#include "G4AssemblyVolume.hh"
#include "LBNEMagneticField.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4PVPlacement.hh"
#include "G4RegionStore.hh"
#include "G4SolidStore.hh"
#include "G4GeometryManager.hh"
#include "G4FieldManager.hh"

#include "LBNERunManager.hh"

#include "G4VisExtent.hh"
#include "LBNEPlacementMessenger.hh"
#include "LBNESurveyor.hh"

//---  Second file for geometry: the first one (LBNEVolumePlacements) was getting to long --// 

void LBNEVolumePlacements::PlaceFinalDownstrTarget(G4PVPlacement *mother) {
    
    if (fUse1p2MW) {
      if (!fUsePseudoNova) PlaceFinalDownstrTarget1p2MW(mother);
      return;
    }
    
    Create("Horn1TargetDownstrHeContainer");
    G4PVPlacement *vMTop = PlaceFinal(std::string("Horn1TargetDownstrHeContainer"), mother);
    LBNEVolumePlacementData *plDowstrHe = Create("Horn1TargetDownstrHelium");
    G4PVPlacement *vMHe = PlaceFinal(std::string("Horn1TargetDownstrHelium"), vMTop);
    // Deal with alignment ring located in the dowstream part of the target. We assume here that 
    // they have been already defined while dealing with the Upstream target portion 
    // First alignment ring, locate flush with the end plate (within 1 mm ) , left and right      
    std::map<G4String, LBNEVolumePlacementData>::iterator itTmpRLeft = fSubVolumes.find(G4String("TargetAlignmentRingLeft"));
    LBNEVolumePlacementData &infoTmpRLeft = itTmpRLeft->second;
    std::map<G4String, LBNEVolumePlacementData>::iterator itTmpRRight = fSubVolumes.find(G4String("TargetAlignmentRingRight"));
    LBNEVolumePlacementData &infoTmpRRight = itTmpRRight->second;
    std::map<G4String, LBNEVolumePlacementData>::iterator itM1 = fSubVolumes.find(G4String("TargetUpstrDownstrHelium"));
    LBNEVolumePlacementData *plM1 = &itM1 ->second;
    std::map<G4String, LBNEVolumePlacementData>::iterator itMother = fSubVolumes.find(G4String("Horn1TargetDownstrHelium"));
    LBNEVolumePlacementData *plMother = &itMother ->second;
    G4ThreeVector posTmp;
    posTmp[0] = 0.; // The alignment rings are always centered.. 
    posTmp[1] = 0.; // We start upstream to keep the separation among rings. 
    posTmp[2] = -1.0*plM1->fParams[2]/2.0 + infoTmpRLeft.fParams[2]/2. + 1.0*CLHEP::mm; // 1 mm spacing Left and right have the same thickness. 
    int copyNumber = 0;
    int copyNumberHere = 0;
    bool transitionDone = false;
    while (true) {
      if (!transitionDone) {
         posTmp[2] += fTargetAlignRingSpacing;      
         copyNumber++;
         if (posTmp[2] > plM1->fParams[2]/2.0) {
           transitionDone = true;
           const double offset = posTmp[2] - plM1->fParams[2]/2.0;
           posTmp[2] =  -1.0*plMother->fParams[2]/2.0 + infoTmpRLeft.fParams[2]/2. + offset ;
         }
      } else {
        std::ostringstream cNumStrStr; cNumStrStr << "_P" << copyNumber;
        if (copyNumber == (fMaxNumAlignRings -1)) {
         // retract the last one a bit by 5 mm to place the return cooling pipe. 
         posTmp[2] -= 10.0*CLHEP::mm;
        }                                 
        new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, infoTmpRLeft.fCurrent, G4String("TargetAlignmentRingLeft")+cNumStrStr.str(), 
                                          vMHe->GetLogicalVolume(), false, copyNumberHere, fCheckVolumeOverLapWC);
        new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, infoTmpRRight.fCurrent, G4String("TargetAlignmentRingRight")+ cNumStrStr.str(), 
                                          vMHe->GetLogicalVolume(), false, copyNumberHere, fCheckVolumeOverLapWC);
        posTmp[2] += fTargetAlignRingSpacing;      
        copyNumber++; copyNumberHere++; 
        if ( copyNumber == fMaxNumAlignRings) break;
      }
    }
    // Deal with the first, split, target segment.  
    LBNEVolumePlacementData *infoTargSegFirst = 0;
    if (fTargetFinLengthSplitDwnstr > 0.3*CLHEP::mm) { 
      infoTargSegFirst = Create("Horn1TargetSegmentFirst");
      G4PVPlacement *vTargSegFirst = PlaceFinal("Horn1TargetSegmentFirst", vMHe);
      Create("Horn1TargetFinVertFirst");
      LBNEVolumePlacementData * plCoolingTubeFirst = Create("Horn1TargetCoolingTubeFirst");
      Create("Horn1TargetCoolingTubeFirstWater");
      posTmp[0] = 0.; posTmp[1] = fTargetFinHeight/2.; posTmp[2] = 0.;                    
      G4PVPlacement *vTubeUp = new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plCoolingTubeFirst->fCurrent, 
                                    G4String("Horn1TargetSegmentCoolingTubeFirst_PTop"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
      posTmp[1] = -1.0*fTargetFinHeight/2.;                       
      G4PVPlacement *vTubeDown =  new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plCoolingTubeFirst->fCurrent, 
                                    G4String("Horn1TargetSegmentCoolingTubeFirst_PBottom"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
      PlaceFinal("Horn1TargetCoolingTubeFirstWater", vTubeDown);                                          
      PlaceFinal("Horn1TargetCoolingTubeFirstWater", vTubeUp);                                    
      PlaceFinal("Horn1TargetFinVertFirst", vTargSegFirst);
//
// Oct. 29 2013 : add the 4 corners as well (overzealous, too much complexity already.. but makes vis. complete.. 
// 
      LBNEVolumePlacementData *plTargetFinCorDFUpLeft = 
        Create(G4String("TargetFinVertCornerDownstrFirstUpLeft"));
      LBNEVolumePlacementData *plTargetFinCorDFUpRight = 
        Create(G4String("TargetFinVertCornerDownstrFirstUpRight"));
      LBNEVolumePlacementData *plTargetFinCorDFDwnLeft = 
       Create(G4String("TargetFinVertCornerDownstrFirstDwnLeft"));
      LBNEVolumePlacementData *plTargetFinCorDFDwnRight = 
      Create(G4String("TargetFinVertCornerDownstrFirstDwnRight"));
   
      posTmp[0] = -fTargetFinWidth/4.; 
      posTmp[1] = fTargetFinHeight/2. -  fTargetCTubeOuterRadius/2. + 0.125*CLHEP::mm;
      posTmp[2] = 0.;    
      new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plTargetFinCorDFUpLeft->fCurrent, 
                                          G4String("Horn1TargetFinVertCornerDownstrFirstUpLeft_P"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
      posTmp[0] = fTargetFinWidth/4.; posTmp[1] = fTargetFinHeight/2. -  fTargetCTubeOuterRadius/2. + 0.125*CLHEP::mm;
      new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plTargetFinCorDFUpRight->fCurrent, 
                                          G4String("Horn1TargetFinVertCornerDownstrFirstUpRight_P"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
                                          
      posTmp[0] = -fTargetFinWidth/4.; posTmp[1] = -1.0*(fTargetFinHeight/2. -  fTargetCTubeOuterRadius/2. + 0.125*CLHEP::mm);
      new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plTargetFinCorDFDwnLeft->fCurrent, 
                                          G4String("Horn1TargetFinVertCornerDownstrFirstDwnLeft_P"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
      posTmp[0] = fTargetFinWidth/4.; posTmp[1] = -1.0*(fTargetFinHeight/2. -  fTargetCTubeOuterRadius/2. + 0.125*CLHEP::mm);
      new G4PVPlacement((G4RotationMatrix *) 0, 
                                    posTmp, plTargetFinCorDFDwnRight->fCurrent, 
                                          G4String("Horn1TargetFinVertCornerDownstrFirstDwnRight_P"), 
                                          vTargSegFirst->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);

     }
    // Now place the previously defined standard target segment. Note: they already contain their cooling and 
    // real target. 
//    LBNEVolumePlacementData  *plHorn1TargetSegment = Create("Horn1TargetSegment");
    std::map<G4String, LBNEVolumePlacementData>::iterator itTargSeg = fSubVolumes.find(G4String("TargetUpstrDownstrSegment"));
    LBNEVolumePlacementData  *plTargetUpstrDownstrSegment= &itTargSeg->second;    
    
    double zCoordTmp = (infoTargSegFirst == 0) ? 
                       (-1.0*plDowstrHe->fParams[2]/2. + plTargetUpstrDownstrSegment->fParams[2]/2. + 0.005*CLHEP::mm) :
                        ( infoTargSegFirst->fPosition[2] + infoTargSegFirst->fParams[2]/2. + 
                       plTargetUpstrDownstrSegment->fParams[2]/2. + 0.002*CLHEP::mm);          
    for (int iSeg=0; iSeg != fTargetNumFinsInHorn; iSeg++) { // Place with no misalignment
      posTmp[0] = 0.; posTmp[1] = 0.;
      posTmp[2] =  zCoordTmp;
//      std::cerr << " In Horn1, Positioning target segment " << iSeg  << " at Z = zCoord " << zCoordTmp << std::endl;
      std::ostringstream cNumStrStr; cNumStrStr << "_P" << iSeg;
         new G4PVPlacement((G4RotationMatrix *) 0, 
                                  posTmp, plTargetUpstrDownstrSegment->fCurrent, 
                                    G4String("Horn1TargetSegment")+cNumStrStr.str(), 
                                        vMHe->GetLogicalVolume(), false, iSeg+100, fCheckVolumeOverLapWC);
      zCoordTmp += (fTargetFinLength + fTargetFinSpacingLength);
    }
    // Now the end 
    LBNEVolumePlacementData *plCoolingTubeReturn = Create("Horn1TargetCoolingTubeDwnRetTit");
    posTmp[0] = 0.; posTmp[1] = 0.; 
    posTmp[2] = plDowstrHe->fParams[2]/2 - 
               fTargetCTubeReturnDownstrThickWater - 1.5*fTargetCTubeReturnDownstrThickTitanium/2. - 2.0*CLHEP::mm;
               //  Hopefully some roo to spare..                          
    new G4PVPlacement(&plCoolingTubeReturn->fRotation, 
                                    posTmp, plCoolingTubeReturn->fCurrent, 
                                    G4String("Horn1TargetCoolingTubeDwnRetTit_PFront"), 
                                          vMHe->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
    posTmp[2] = plDowstrHe->fParams[2]/2 -  0.5*fTargetCTubeReturnDownstrThickTitanium/2. - 1.3*CLHEP::mm;                        
    new G4PVPlacement(&plCoolingTubeReturn->fRotation, 
                                  posTmp, plCoolingTubeReturn->fCurrent, 
                                  G4String("Horn1TargetCoolingTubeDwnRetTit_PBack"), 
                                        vMHe->GetLogicalVolume(), false, 0, fCheckVolumeOverLapWC);
    
    Create("Horn1TargetCoolingTubeDwnRetWater");
    PlaceFinal("Horn1TargetCoolingTubeDwnRetWater", vMHe);

    Create("Horn1TargetDownstrHeContainerCap");
    PlaceFinal("Horn1TargetDownstrHeContainerCap", vMHe);                                                                         
}
//
// Horn1 calls and utilities. 
//

void LBNEVolumePlacements::RescaleHorn1Lengthwise() {

    if (fHorn1LongRescaleCnt != 0) {
      G4Exception("LBNEVolumePlacements::RescaleHorn1Lengthwise", "", 
                FatalErrorInArgument, " This method can only be called once per job ");
    }
    fHorn1LongRescaleCnt++;
    fHorn1Length  *= fHorn1LongRescale; 
    fHorn1IOTransLength *= fHorn1LongRescale;
     
   for (size_t k=0; k != fHorn1UpstrLengths.size(); ++k) 
      fHorn1UpstrLengths[k] *=  fHorn1LongRescale; 
   for (size_t k=0; k != fHorn1UpstrZPositions.size(); ++k) 
      fHorn1UpstrZPositions[k] *=  fHorn1LongRescale; 
   
  for (size_t k=0; k != fHorn1UpstrOuterIOTransLengths.size(); ++k)
      fHorn1UpstrOuterIOTransLengths[k] *= fHorn1LongRescale; 
      
   for (size_t k=0; k != fHorn1UpstrOuterIOTransPositions.size(); ++k)
      fHorn1UpstrOuterIOTransPositions[k] *= fHorn1LongRescale; 
      
  for (size_t k=0; k != fTargetHorn1Lengths.size(); ++k) 
      fTargetHorn1Lengths[k] *= fHorn1LongRescale;
   for (size_t k=0; k != fTargetHorn1ZPositions.size(); ++k) 
      fTargetHorn1ZPositions[k] *= fHorn1LongRescale;  // This will be approximate...
      
  fHorn1TopUpstrLength *= fHorn1LongRescale;
  fHorn1TopDownstrLength *= fHorn1LongRescale;
  
  // Now we deal with the inner conductor radial equations.
  for (std::vector<LBNEHornRadialEquation>::iterator it=fHorn1Equations.begin();
        it != fHorn1Equations.end(); it++) it->SetLongRescaleFactor(fHorn1LongRescale);
        
  fHorn1NeckLength *= fHorn1LongRescale;
  fHorn1NeckZPosition *= fHorn1LongRescale;
  
//   for (size_t k=0; k != fHorn1SectionLengths.size(); ++k) 
//       fHorn1SectionLengths[k] *= fHorn1LongRescale;  
//  fHorn1OuterTubeInnerLength  *= fHorn1LongRescale; // ??????????????///
  fHorn1OuterConnectorLength *= fHorn1LongRescale;
  fHorn1InnerConnectorLength *= fHorn1LongRescale;
//  fHorn1InnerTopUpstrLength *= fHorn1LongRescale;  // ??????????? Check that all the rescale is consistently... 
//  fHorn1InnerTopDownstrLength  *= fHorn1LongRescale;
  
//
// Aug 18 2014., for simplified geometry.   
  for(std::vector<double>::iterator it = fMotherHorn1AllLengths.begin(); it != fMotherHorn1AllLengths.end(); it++) 
      (*it) *= fHorn1LongRescale;
  
}

void LBNEVolumePlacements::RescaleHorn1Radially() { // Note: not all variables are rescale here, 
    if (fHorn1RadialRescaleCnt != 0) {
      G4Exception("LBNEVolumePlacements::RescaleHorn1Radially", "", 
                FatalErrorInArgument, " This method can only be called once per job ");
    }
    fHorn1RadialRescaleCnt++;
  // Some of this work will be done via the LBNEHornRadialEquation
  fHorn1IOTransInnerRad *= fHorn1RadialRescale;
  fHorn1IOTransOuterRad *= fHorn1RadialRescale;
  
  for (size_t k=0; k != fHorn1UpstrOuterIOTransInnerRads.size(); ++k)
      fHorn1UpstrOuterIOTransInnerRads[k] *= fHorn1RadialRescale; 
  for (size_t k=0; k != fHorn1UpstrOuterIOTransThicks.size(); ++k)
      fHorn1UpstrOuterIOTransThicks[k] *= fHorn1RadialRescale; 
      
  for (size_t k=0; k != fTargetHorn1InnerRadsUpstr.size(); ++k) 
      fTargetHorn1InnerRadsUpstr[k] *= fHorn1RadialRescale;
  for (size_t k=0; k != fTargetHorn1InnerRadsDownstr.size(); ++k) 
      fTargetHorn1InnerRadsDownstr[k] *= fHorn1RadialRescale;
  for (size_t k=0; k != fTargetHorn1TransThick.size(); ++k) 
      fTargetHorn1TransThick[k] *= fHorn1RadialRescale;

  for (size_t k=0; k != fHorn1UpstrInnerRadsUpstr.size(); ++k) 
      fHorn1UpstrInnerRadsUpstr[k] *= fHorn1RadialRescale;
  for (size_t k=0; k != fHorn1UpstrInnerRadsDownstr.size(); ++k) 
      fHorn1UpstrInnerRadsDownstr[k] *= fHorn1RadialRescale;
   
  for (size_t k=0; k != fHorn1UpstrInnerRadsOuterUpstr.size(); ++k)
    fHorn1UpstrInnerRadsOuterUpstr[k] *= fHorn1RadialRescale;
  for (size_t k=0; k != fHorn1UpstrInnerRadsOuterDownstr.size(); ++k)
    fHorn1UpstrInnerRadsOuterDownstr[k] *= fHorn1RadialRescale;
  
//  fHorn1TopUpstrInnerRad  *= fHorn1RadialRescale;  Recomputed in CreateHorn1TopLevelUpstr

  fHorn1OuterTubeInnerRad  *= fHorn1RadialRescale;; 
  fHorn1OuterTubeOuterRad  *= fHorn1RadialRescale;; 
  fHorn1TopUpstrOuterRad  *= fHorn1RadialRescale;
  
  fHorn1TopDownstrOuterRad  *= fHorn1RadialRescale; 
  // Now we deal with the inner conductor radial equations.
  for (std::vector<LBNEHornRadialEquation>::iterator it=fHorn1Equations.begin();
        it != fHorn1Equations.end(); it++) it->SetRadialRescaleFactor(fHorn1RadialRescale);

//
// Aug 18 2014., for simplified geometry.   
  for(std::vector<double>::iterator it = fMotherHorn1AllRads.begin(); it != fMotherHorn1AllRads.end(); it++) 
      (*it) *= fHorn1RadialRescale;
}  

void LBNEVolumePlacements::DeclareHorn1Dims() {
  
  const double in = 2.54*CLHEP::cm;
  fHorn1IOTransLength = 3.0*CLHEP::cm + 3.316*in + 0.005*CLHEP::mm; 
                          // Drawing 8875.112 -MD-363097 The 3 cm is the MCZERO offset, per verbal discussion 
                              // with J. Hylen. The 5 microns if to avoid G4 volume overlaps.  
                              
//  fHorn1RadialSafetyMargin = 2.9*CLHEP::mm; // per agreement between Jim H. and Alberto M., Aug. 22 2013. 
  fHorn1RadialSafetyMargin = 2.5*CLHEP::mm; // per agreement between Jim H. and Alberto M., Oct 4 2013 
  
  fHorn1IOTransInnerRad = 2.520*in/2. - fHorn1RadialSafetyMargin/2. ; // last term is the 
  fHorn1IOTransOuterRad = 16.250*in/2.;
  
  fHorn1UpstrInnerRadsUpstr.resize(4);
  fHorn1UpstrInnerRadsDownstr.resize(4);
  fHorn1UpstrInnerRadsOuterUpstr.resize(4);
  fHorn1UpstrInnerRadsOuterDownstr.resize(4);
  fHorn1UpstrLengths.resize(4);
  fHorn1UpstrZPositions.resize(4);
  
  fHorn1UpstrInnerRadsUpstr[0] = 1.572*in; 
  fHorn1UpstrInnerRadsOuterUpstr[0] = fHorn1UpstrInnerRadsUpstr[0] + 0.32*in; 
  fHorn1UpstrInnerRadsDownstr[0] = 1.41*in; 
  fHorn1UpstrInnerRadsOuterDownstr[0] = fHorn1UpstrInnerRadsDownstr[0] + 0.32*in; 
  fHorn1UpstrLengths[0] = 0.508*in - 0.100*CLHEP::mm;
  fHorn1UpstrZPositions[0] = fHorn1UpstrLengths[0]/2. + 0.025; // With respect to the beginning of mother volume.  
  
  fHorn1UpstrInnerRadsUpstr[1] = fHorn1UpstrInnerRadsDownstr[0]; 
  fHorn1UpstrInnerRadsOuterUpstr[1] = fHorn1UpstrInnerRadsUpstr[1] + 0.32*in; 
  fHorn1UpstrInnerRadsDownstr[1] = 1.288*in; 
  fHorn1UpstrInnerRadsOuterDownstr[1] = fHorn1UpstrInnerRadsDownstr[1] + 0.3*in; 
  fHorn1UpstrLengths[1] = 0.639*in  - 0.100*CLHEP::mm;
  fHorn1UpstrZPositions[1] = fHorn1UpstrZPositions[0] + 0.025*CLHEP::mm + fHorn1UpstrLengths[0]/2 + fHorn1UpstrLengths[1]/2.;
  
  fHorn1UpstrInnerRadsUpstr[2] = fHorn1UpstrInnerRadsDownstr[1]; 
  fHorn1UpstrInnerRadsOuterUpstr[2] = fHorn1UpstrInnerRadsUpstr[2] + 0.28*in; 
  fHorn1UpstrInnerRadsDownstr[2] = 1.268*in; 
  fHorn1UpstrInnerRadsOuterDownstr[2] = fHorn1UpstrInnerRadsDownstr[2] + 0.118*in; 
  fHorn1UpstrLengths[2] = 2.835*in  - 0.100*CLHEP::mm; // Reduce a bit, too tight..
  fHorn1UpstrZPositions[2] = fHorn1UpstrZPositions[1] + 0.025*CLHEP::mm + fHorn1UpstrLengths[1]/2 + fHorn1UpstrLengths[2]/2.;
  
  fHorn1UpstrInnerRadsUpstr[3] = fHorn1UpstrInnerRadsDownstr[2]; 
  fHorn1UpstrInnerRadsOuterUpstr[3] = fHorn1UpstrInnerRadsUpstr[2] + 0.118*in; 
  fHorn1UpstrInnerRadsDownstr[3] =  fHorn1UpstrInnerRadsUpstr[3]; 
  fHorn1UpstrInnerRadsOuterDownstr[3] = fHorn1UpstrInnerRadsDownstr[3] + 0.20*in; 
  fHorn1UpstrLengths[3] = 0.479*in - 0.100*CLHEP::mm; // Nov 19 2013 : change from 0.40 inches to 0.479, per more 
                                               // accurate reading of drawing 363097
  fHorn1UpstrZPositions[3] = fHorn1UpstrZPositions[2] + 0.025*CLHEP::mm + fHorn1UpstrLengths[2]/2 + fHorn1UpstrLengths[3]/2.;
 //
 // These are elements, approximated as tubes, of the Inner Outer transition piece of Horn1 
 // 
  fHorn1UpstrOuterIOTransInnerRads.resize(4);
  fHorn1UpstrOuterIOTransThicks.resize(4);
  fHorn1UpstrOuterIOTransLengths.resize(4);
  fHorn1UpstrOuterIOTransPositions.resize(4);
  
  fHorn1UpstrOuterIOTransInnerRads[0] = 5.625*in;
  fHorn1UpstrOuterIOTransThicks[0] = 0.30*in;
  fHorn1UpstrOuterIOTransLengths[0] = 0.421*in;
  fHorn1UpstrOuterIOTransPositions[0] = fHorn1UpstrOuterIOTransLengths[0]/2. + 0.005*CLHEP::mm;

  fHorn1UpstrOuterIOTransInnerRads[1] = 5.763*in;
  fHorn1UpstrOuterIOTransThicks[1] = 0.32*in;
  fHorn1UpstrOuterIOTransLengths[1] = 0.421*in;
  fHorn1UpstrOuterIOTransPositions[1] = 
  fHorn1UpstrOuterIOTransLengths[0]+ fHorn1UpstrOuterIOTransLengths[1]/2. + 0.015*CLHEP::mm;
  
  fHorn1UpstrOuterIOTransInnerRads[2] = 5.763*in;
  fHorn1UpstrOuterIOTransThicks[2] = 0.4*in;
  fHorn1UpstrOuterIOTransLengths[2] = 0.75*in;
  fHorn1UpstrOuterIOTransPositions[2] = fHorn1UpstrOuterIOTransLengths[0] +
  fHorn1UpstrOuterIOTransLengths[1] + fHorn1UpstrOuterIOTransLengths[2]/2. + 0.050*CLHEP::mm;

  fHorn1UpstrOuterIOTransInnerRads[3] = 5.763*in;
  fHorn1UpstrOuterIOTransThicks[3] = 2.26*in;
  fHorn1UpstrOuterIOTransLengths[3] = fHorn1IOTransLength - 0.250*CLHEP::mm - 
    fHorn1UpstrOuterIOTransLengths[0] - fHorn1UpstrOuterIOTransLengths[1] - fHorn1UpstrOuterIOTransLengths[2];
    
  fHorn1UpstrOuterIOTransPositions[3] = fHorn1IOTransLength - fHorn1UpstrOuterIOTransLengths[3]/2. - 0.01*CLHEP::mm;;

// Outer dimension of the big tube that contains almost everything.. 
// Except the big connedtors rings downsream (max dim are 23.5 in 
     
  fHorn1OuterTubeOuterRad =  13.750*in/2.0;
  fHorn1OuterTubeInnerRad =  11.750*in/2.0; // checked with drawing 364094
 // 

  //
 // Now the top level sections. Only two of them, alignable. 
   fHorn1TopUpstrOuterRad = 16.250*in/2. +  5.0*CLHEP::mm; // safety 5 mm, big radius..// Drawing 8875.112-MD 363097
   fHorn1NeckLength = 1.568*in; // Drawing 8875.112-MD  363105
   fHorn1NeckZPosition = 30.315*in + fHorn1NeckLength/2.;  // Drawing 8875.112-MD  363105
      // In Z Drawing coordinates. 
      
   fHorn1TopDownstrOuterRad =  fHorn1TopUpstrOuterRad;  // Initialization finally fixed on March 18 2014... 
                                                        // Should have been done 9 months ago!   

   fHorn1TopDownstrLength  = 133.752*in; // We will have to subtract the upstream length, once we know it 
   // This include the downstream bulkhead. 
  // Other quantities defined in CreateHorn1TopLevelUpstr..
 // The first one goes up the downstream end of the IO transition peace 
 // and up the downstream end of the target. If the target is completly out of the horn, 
 // it is the part before the neck. (Drawing 8875.112-MD-363104)
 //  This is a tube, set by the inner conductor
 // radius at the end of the target.  
 // 
 // Load the equations Drawing 8875.112-MD 363104, 363105
//
  fHorn1Equations.clear();
  LBNEHornRadialEquation e1(1.975805, -0.055858, -0.078740); fHorn1Equations.push_back(e1);
  LBNEHornRadialEquation e2(1.818869, -0.055858, 0.); fHorn1Equations.push_back(e2);
  LBNEHornRadialEquation eBlank(0., 0.,  0.); fHorn1Equations.push_back(eBlank); // equation 3 not found on drawing. 
  LBNEHornRadialEquation e4(-5.619190, 0.180183, 0.); fHorn1Equations.push_back(e4); 
  LBNEHornRadialEquation e5(-5.462253, 0.180183, -0.078740); fHorn1Equations.push_back(e5); 
  LBNEHornRadialEquation e6(1.97805, -0.055858, 0.); fHorn1Equations.push_back(e6); 
  fHorn1Equations.push_back(eBlank); // equation 7 not found on drawing as well 
  LBNEHornRadialEquation e8(-5.462253, 0.180183, 0.); fHorn1Equations.push_back(e8); 
 
 // The 2nd top lelvel section is the remaining part of that part, up to the end of the Horn1Hall Z location. 
 // This is just a tube. 
 //
  fHorn1InnerCondMat = G4String("Aluminum");  
  fHorn1AllCondMat = G4String("Aluminum");  
 // Then, conical sections, all in a tube that has an inner radius of the neck, our radius 
 // 
 // To compute all that, we need to locate the target first.. 
 // So, done at construction time. 
 //
  
}
void LBNEVolumePlacements::setMotherVolumeForHorn1() {
  const double in = 2.54*CLHEP::cm;
  fMotherHorn1AllRads.clear();
  fMotherHorn1AllLengths.clear();
  const double epsilR = 0.050*CLHEP::mm;
  double zCurrent = 0.;
  fMotherHorn1AllRads.push_back(fHorn1IOTransInnerRad - epsilR);
  fMotherHorn1AllLengths.push_back(zCurrent);
  fMotherHorn1AllRads.push_back(fHorn1IOTransInnerRad - epsilR);
//  zCurrent += fHorn1IOTransLength; // Oversize ? ! Drawing 363092, 363097 says (2.436" +  3.2852"  =  5.752 ) 
  zCurrent += 5.752*in; ; // Correct...
  fMotherHorn1AllLengths.push_back(zCurrent);
  // Up to the neck, equation 
  const int numSectEq1 = 5;
  const double lengthEq1 = 12.896*in/numSectEq1; // Drawin 363104 
  double zz1 = 3.316*in;
  for(int iSect1 = 0; iSect1 != numSectEq1; iSect1++) {
    zz1 += lengthEq1;
    const double rr1 = fHorn1Equations[0].GetVal(zz1) - epsilR - 0.300*CLHEP::mm -0.050*CLHEP::mm*iSect1;// needs a bit more room here, these subsection are too straight..
    zCurrent += lengthEq1;
    fMotherHorn1AllRads.push_back(rr1); 
    fMotherHorn1AllLengths.push_back(zCurrent);
//    std::cerr << " Polycone Mother, subsection 1 , zz1 " << zz1 << " zCurrent " << zCurrent << " rr1 " << rr1 << std::endl;
  }
  const int numSectEq2 = 8;
  const double lengthEq2 = (4.9799*in + 9.2262*in)/numSectEq2; // Drawing 363104, 363105 
  for(int iSect2 = 0; iSect2 != numSectEq2; iSect2++) {
    zz1 += lengthEq2;
    const double rr2 = fHorn1Equations[1].GetVal(zz1)  - epsilR - 0.8*CLHEP::mm - 0.1*CLHEP::mm*iSect2;// empirical, but it will need to be improved on for 
                                                          // misalignment studies. 
    zCurrent += lengthEq2;
    fMotherHorn1AllRads.push_back(rr2); 
    fMotherHorn1AllLengths.push_back(zCurrent);
//    std::cerr << " Polycone Mother, subsection 2 , zz1 " << zz1 << " zCurrent " << zCurrent << " rr2 " << rr2 << std::endl;
  }
  // The neck + all the other sections, set to the neck. 
  const double neckRadius = 0.709*in/2.; // Drawing 363105 
  fMotherHorn1AllRads.push_back(neckRadius - epsilR);
  fMotherHorn1AllLengths.push_back(zCurrent);
  zCurrent += 1.568*in; zz1 += 1.568*in;
  fMotherHorn1AllRads.push_back(neckRadius - epsilR);
  fMotherHorn1AllLengths.push_back(zCurrent);
//  std::cerr << " ... Just downstream of the neck, zz1 = " << zz1 << " zCurrent " << zCurrent << std::endl;
  // Downstream of the neck.. Equation 4 
  const int numSectEq4 = 4;
  const double lengthEq4 = (36.2709 - 31.8827)*in/numSectEq4; // Drawing 363105 
  for(int iSect4 = 0; iSect4 != numSectEq4; iSect4++) {
    zz1 += lengthEq4;
    double rr4 = fHorn1Equations[3].GetVal(zz1)  - epsilR - 0.2*CLHEP::mm - 0.1*CLHEP::mm*iSect4;// empirical, but it will need to be improved on for 
                                                          // misalignment studies. 
    if (iSect4 == 0) rr4 = (rr4 + neckRadius)/2.0; // take the average... Not really important, a realistic plug should nevr fit tightly. 
    zCurrent += lengthEq4;
    fMotherHorn1AllRads.push_back(rr4); 
    fMotherHorn1AllLengths.push_back(zCurrent);
//    std::cerr << " Polycone Mother, subsection 4 , zz1 " << zz1 << " zCurrent " << zCurrent << " rr4 " << rr4 << std::endl;
  }
  // Downstream of the neck.. Equation 5 Up to the downstream flange, drawing 363108, 363096
  const int numSectEq5 = 20;
  const double lengthEq5 = (117.126 - 36.27097)*in/numSectEq5; // Drawing 363105, 363108 
  double rLastParabolic = 0.;
  for(int iSect5 = 0; iSect5 != numSectEq5; iSect5++) {
    zz1 += lengthEq5;
    const double rr5 = fHorn1Equations[4].GetVal(zz1)  - epsilR - 1.0*CLHEP::mm - 0.1*CLHEP::mm*iSect5;// empirical, but it will need to be improved on for 
                                                          // misalignment studies. 
    zCurrent += lengthEq5;
    fMotherHorn1AllRads.push_back(rr5); 
    fMotherHorn1AllLengths.push_back(zCurrent);
//    std::cerr << " Polycone Mother, subsection 5 , zz1 " << zz1 << " zCurrent " << zCurrent << " rr5 " << rr5 << std::endl;
    rLastParabolic = rr5;
  }  
  zCurrent += 14.244*in + 3.0*CLHEP::cm; // xx cm safe for misalignment Drawing 363093 and 363096  
  fMotherHorn1AllRads.push_back(rLastParabolic - epsilR);
  fMotherHorn1AllLengths.push_back(zCurrent);
  //
  // Now the outer...
  // 
  const double rOut = 16.25*in + 15.0*in + 2.5*CLHEP::cm;
  fMotherHorn1AllRads.push_back(rOut);
  fMotherHorn1AllLengths.push_back(zCurrent);
  fMotherHorn1AllRads.push_back(rOut);
  fMotherHorn1AllLengths.push_back(0.);
  //
  // dump this curve, so we can draw it.
//  std::ofstream fileOutTmp("./RZCurve_1.txt");
//  fileOutTmp << " Z R " << std::endl;
//  for (size_t ii=0 ; ii != fMotherHorn1AllLengths.size(); ii++) {
//     std::cerr << " Radius At at  Z(StartHorn1) "<< fMotherHorn1AllLengths[ii]  
//             << "  is " << fMotherHorn1AllRads[ii] << std::endl;
//     fileOutTmp <<  " " << fMotherHorn1AllLengths[ii] << " " << fMotherHorn1AllRads[ii] << std::endl;     
//  }
//  fileOutTmp.close();
//  std::cerr << " And quit !!! ... " << std::endl; exit(2);
}

LBNEVolumePlacementData* LBNEVolumePlacements::CreateHorn1TopLevelUpstr() { 

     // A bit special, case, we have to compute the inner conductor radius..
     // at the length of the target.
     
  const double in = 2.54*CLHEP::cm;   
  fHorn1TopUpstrLength  =  fTargetHeContTubeLengthInHorn -  // Not rescaled!!! 
                           fHorn1LongRescale*(3.0*CLHEP::cm - 3.316*in + 7.0*CLHEP::mm); 
                           // The last factor is a tweak, coming from the various safety marging and/or arbitrariness 
                           // in defining the geometry.  
  // The transition Inner to Outer is now included in the Horn1TopLevelUpstr volume!!
  fHorn1TopUpstrLength += fHorn1IOTransLength;
  //
  // The above equation to determine fHorn1TopUpstrLength looks wrong based on Geantino analysis
  // Why not simply:
  
  fHorn1TopUpstrLength  =  fTargetHeContTubeLengthInHorn + 7.0*CLHEP::mm;
  //
  // Wrong if the entire target is out of Horn1.. 
  //
  if (fUsePseudoNova) fHorn1TopUpstrLength = fHorn1IOTransLength + + 3.36*in + 60.0*CLHEP::mm + 3.0*CLHEP::mm; // The 60 mm 
  // is very arbitrary.  It provide enough space for the weld of the I/O to the most upstream inner conductor 
  // piece.  Note that, if the target is completly outside Horn1, the raison d'etre 
  // of the split Upst/Downstream part of the Horn1 goes away.. But the structure of the code 
  // was there, try to coopt it.. 
  //
  
  
  
//  std::cerr << "  CreateHorn1TopLevelUpstr : fHorn1IOTransLength " << fHorn1IOTransLength
//             <<  " fTargetHeContTubeLengthInHorn " << fTargetHeContTubeLengthInHorn 
//           <<  " fHorn1TopUpstrLength " << fHorn1TopUpstrLength << " and quit " << std::endl;
//           exit(2);              
  double zEnd = fHorn1TopUpstrLength - fHorn1LongRescale*3.0*CLHEP::cm;  // Crappy convention again.. 
  if (zEnd > fHorn1LongRescale*30.3150*in) { // Z (drawing coordinate system) 
     std::ostringstream messStrStr; 
    messStrStr << " Helium Tube target hitting the neck! zEnd " << zEnd << "\n";
    G4String messStr(messStrStr.str());
    G4Exception("LBNEVolumePlacements::CreateHorn1TopLevelUpstr", " ", 
                   FatalErrorInArgument, messStr.c_str()); 
  }
  fHorn1TopUpstrInnerRad = fHorn1Equations[0].GetVal(zEnd) - 0.005*CLHEP::mm; // recaled radially before construction, if need be.
  return this->Create("Horn1TopLevelUpstr");
}
// Checking the possible radial overlap at the downstream end of the target. 
// The size and position of Horn1TopLevelUpstr and Horn1TargetDownstrHeContainer are fixed... 

void LBNEVolumePlacements::CheckHorn1InnerConductAndTargetRadii() {

  if (fRemoveTargetAlltogether || fUseTargetModule) return;
  const double in = 2.54*CLHEP::cm;   
  const LBNEVolumePlacementData *plSmallRing;
  const LBNEVolumePlacementData *plLargeRing;
  plSmallRing  = this->Find(G4String("Check"), 
                       G4String("TargetNoSplitHeContainer"), 
                       G4String("LBNEVolumePlacements::CheckHorn1InnerConductAndTargetRadii"));

  plLargeRing = this->Find(G4String("Check"), 
                        G4String("Horn1PolyM1"), 
                        G4String("LBNEVolumePlacements::CheckHorn1InnerConductAndTargetRadii"));

  const double xOffS = plSmallRing->fPosition[0] + 
                          (plSmallRing->fParams[2]/2.)*plSmallRing->fRotation.xz();     // Assumed misaligned once (withe respect to perfectly 
                                                                                       // aligned (HORN1. )                                         
  const double yOffS = plSmallRing->fPosition[1] + 
                          (plSmallRing->fParams[2]/2.)*plSmallRing->fRotation.yz();
                          
                          // We take the full length, the clash is at the downstream end. 
                          
  double xOffL = plLargeRing->fPosition[0] + 
                          (plLargeRing->fParams[2]/2.)*plLargeRing->fRotation.xz();
                                                                                                                            
  double yOffL = plLargeRing->fPosition[1] + 
                          (plLargeRing->fParams[2]/2.)*plLargeRing->fRotation.yz();
                          // Same, at the downstream end. But this is only approximate. 
  double maxLength = -100000.0*CLHEP::m;
  for (std::vector<double>::const_iterator il = fMotherHorn1AllLengths.begin(); il != fMotherHorn1AllLengths.end(); il++) 
      if (*il > maxLength ) maxLength = *il;
  xOffL = plLargeRing->fPosition[0] + maxLength*plLargeRing->fRotation.xz()/2.;
                                                                                                                         
  yOffL = plLargeRing->fPosition[1] + maxLength*plLargeRing->fRotation.yz()/2.;
  const double xOff = xOffS - xOffL;// now translated to ref. frame of the large ring. 
  const double yOff = yOffS - yOffL;
  
  const double rSmall = plSmallRing->fParams[1]; // Outer radius
  // Compute the large radius based on equation 0, at the relevant Z 
  const double zMaxDC  = fHorn1TopUpstrLength - 3.0*CLHEP::cm;
  const double rMinEqn1 = fHorn1Equations[0].GetVal(zMaxDC); // Equation 1 or 0
  const double rMinEqn2 = fHorn1Equations[1].GetVal(zMaxDC); // Equation 1 or 0
  const double ratio0vs1 = std::min(1.0, (zMaxDC/(21.0888*in*fHorn1LongRescale)));
  const double rLarge = rMinEqn2*ratio0vs1 + (1.0-ratio0vs1)*rMinEqn1; 
              // Equation 1 and 0, weighted by Z ;  // Inner Radius
//  std::cerr << " CheckHorn1InnerConductAndTargetRadii, zCoordEnd " 
//            << zCoordEnd << " rLarge " << rLarge << " rSmall " << rSmall << std::endl;
//  std::cerr << " .. xOffS " <<        xOffS << "  yOffS  " <<  yOffS 
//            << " xOffL " << xOffL  << " yOffL " << yOffL << std::endl;
  G4String messHeader(" The downstream end target is touching (or nearly touching)  the Horn1 inner conductor \n");
  if ((std::abs(xOff) < 0.002*CLHEP::mm) && (std::abs(yOff) < 0.002*CLHEP::mm)) {
    if (rSmall + fHorn1RadialSafetyMargin < rLarge) return;
    
    std::ostringstream messStrStr; 
    messStrStr << messHeader << " Helium Tube radius " << rSmall
                                  << " Inner cond. radius " << rLarge << " \n";
    messStrStr << messHeader << " Suggestion: Reduce fTargetLengthIntoHorn, currently set to  " << 
                                fTargetLengthIntoHorn <<  " mm \n";
    G4String messStr(messStrStr.str());
    G4Exception("LBNEVolumePlacements::CheckHorn1InnerConductAndTargetRadii", " ", 
                   FatalErrorInArgument, messStr.c_str()); 
  } else { // Misalignment
//     const double phic = std::atan2(yOff, xOff); // extremum in phis . We don't care where it is in phi.. 
//     double rMaxDev = -1.0;
//     for (int iQ=0; iQ!=4; iQ++) { //. but which quadrant? (Should be able to figure out this more simply..) 
//        const double xCrit = rSmall*std::sin(phic);                                        
//        const double yCrit = rSmall*std::cos(phic);
//      const double rr = std::sqrt(xCrit*xCrit + yCrit*yCrit);
//      if (rr > rMaxDev) rMaxDev = rr;
//    }
    const double  rMaxDev = std::sqrt(xOff*xOff + yOff*yOff);
//    std::cerr << " Misalignments detected in CheckHorn1InnerConductAndTargetRadii, phic " << phic 
//              << " rMaxDev  " << rMaxDev << std::endl;
    if ((rMaxDev + rSmall) < (rLarge - fHorn1RadialSafetyMargin)) return; // O.K> 
    std::ostringstream messStrStr; 
    messStrStr << messHeader << " Helium Tube radius " << rSmall
                                  << " Inner cond. radius " << rLarge << " \n";
    messStrStr << messHeader << " Tranvser Offset between the two rings  " << xOff << " / " << yOff << "\n";
    G4String messStr(messStrStr.str());
    G4Exception("LBNEVolumePlacements::CheckHorn1InnerConductAndTargetRadii", " ", 
                   FatalErrorInArgument, messStr.c_str()); 
  }  
}

//
// Subdivide a conical section such that the radial equations are correct within a 
// a tolerance of 5 microns. We assume that the number of subdivisions
// estimated at the upstream end of the section  will 
// be adequate over the entire section. 
//
int LBNEVolumePlacements::GetNumberOfInnerHornSubSections(size_t eqn, double z1, double z2, int nMax) const {

   int numSubSect = 1;
   double zEnd = z2;
   const double zLengthTotal = z2 - z1;
   double zLength = zLengthTotal;
   while (true) { 
     const double r0 = fHorn1Equations[eqn].GetVal(z1);
     const double r1 = fHorn1Equations[eqn].GetVal(zEnd);
     const double rMid = fHorn1Equations[eqn].GetVal((z1 + zEnd)/2.);
//     std::cerr << " GetNumberOfInnerHornSubSections z1 " << z1 << " zEnd " << zEnd << 
//                     " r0 " << r0 << " r1 " << r1 << " rMid " << rMid << std::endl;
     if (std::abs(rMid - 0.5*(r0+r1)) < 0.050*CLHEP::mm) break; // 50 microns good enough for Gov. work. 
                                                         // meachanical tolerance is one 5 mills, (127 microns)
     zLength /=2.;
     zEnd = z1 + zLength;
     numSubSect++;
     if (numSubSect > nMax) {
       G4Exception("LBNEVolumePlacements::GetNumberOfInnerHornSubSections", " ", FatalErrorInArgument, " Crazy subdivision! ");
       break;
     }
   }
   return numSubSect;
}
//
// 2nd argument is in drawing coordinate system, the upstream coordinate of the beginning of the ring
//  Third argument, in the reference frame of the mother volume, the long. center position 
//
void LBNEVolumePlacements::Horn1InstallSpiderHanger(const G4String &nameStrH, 
                                                    double zLocTweakedDC, double zPosMotherVolume, 
                                                   G4PVPlacement *vMother ) {

  const double in = 2.54*CLHEP::cm;
  G4String nameStr(nameStrH);nameStr += G4String("Ring");
  const double length = 0.750*in*fHorn1LongRescale;
  const int eqnNum = (zLocTweakedDC < (41.076*fHorn1LongRescale*in)) ? 5 : 7;
  const double zSignLength = (eqnNum == 5) ? -1.0 : 1.0; // to avoid collision with the inner conductor outer radius 
                                                         // at the upstream or downstream end 
  double rTmp1 = fHorn1Equations[eqnNum].GetVal(zLocTweakedDC + length*zSignLength) 
                          + 0.0015*CLHEP::mm + fWaterLayerThickInHorns;
  if (fHorn1RadiusBigEnough) {
     const double rr = fHorn1Equations[eqnNum].GetVal(zLocTweakedDC + length*zSignLength);
     if (rr < fHorn1RMinAllBG ) rTmp1 = fHorn1RMinAllBG + 0.0015*CLHEP::mm + fWaterLayerThickInHorns;
  }                       
  const double rTmp2 = rTmp1 + 0.24*in; // Deduced from 363104 and equation 6
  G4Tubs *aTubs = new G4Tubs(nameStr, rTmp1, rTmp2, 
                                length/2.   , 0., 360.0*CLHEP::degree);
  G4LogicalVolume *pCurrent = new G4LogicalVolume(aTubs, G4Material::GetMaterial(fHorn1InnerCondMat), nameStr);
  G4ThreeVector posTmp; posTmp[0] = 0.; posTmp[1] = 0.;
  posTmp[2] = zPosMotherVolume;                            
  new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrent, nameStr + std::string("_P"), 
                     vMother->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
  
  //The connecting piece ring to the hangers.. There are three of them, at 120 degrees from each other. 
  
  G4String nameStr2(nameStrH); nameStr2 += G4String("Riser");
  const double heightRiser = 0.333*in - 0.020*CLHEP::mm;
  const double widthH = 1.5*in; // See drawing 8875.112-MD 363115
  const double thickH = 0.184*2*in; 
  G4Box *aBoxRiser = new G4Box(nameStr2, widthH/2., heightRiser/2.0, thickH/2.0);  
  G4LogicalVolume *pCurrentRiser = 
    new G4LogicalVolume(aBoxRiser, G4Material::GetMaterial(fHorn1InnerCondMat), nameStr2);
    
  G4String nameStr3(nameStrH); nameStr3 += G4String("Hanger");
  const double heightH = fHorn1OuterTubeInnerRad - rTmp2 - 1.0*CLHEP::mm - heightRiser;
  const double widthH2 = 1.0*in; // 363115 Note: we collapsed both hanger along the horizontal, transverse 
                                // direction. 
  const double thickH2 = 0.031*in; 
  G4Box *aBoxHanger = new G4Box(nameStr3, widthH2/2., heightH/2.0, thickH2/2.0);  
  G4LogicalVolume *pCurrentHanger = 
    new G4LogicalVolume(aBoxHanger, G4Material::GetMaterial(fHorn1InnerCondMat), nameStr3);
  
  for (int iRot=0; iRot != 3; iRot++) {
    std::ostringstream rnStrStr; rnStrStr << "_" << (iRot+1);
    G4String rnStr(rnStrStr.str());
  // Same Z position as above... 
    G4RotationMatrix * rMatPr = 0;
     if (iRot == 1) {
      rMatPr = new G4RotationMatrix; 
      rMatPr->rotateZ(2.0*M_PI/3.);
    } else if (iRot == 2) {
      rMatPr = new G4RotationMatrix; 
      rMatPr->rotateZ(-2.0*M_PI/3.);
    }
    
    const double dHRiser = rTmp2 + 0.010*CLHEP::mm + heightRiser/2.;
    posTmp[0] = 0.;  posTmp[1] = dHRiser;
    if (iRot != 0) {
      posTmp[0] = dHRiser*rMatPr->xy();
      posTmp[1] = dHRiser*rMatPr->yy();
    }
    if (iRot == 0) { 
       new G4PVPlacement(rMatPr, posTmp, pCurrentRiser, nameStr2 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
    } else {
       new G4PVPlacement(G4Transform3D(*rMatPr, posTmp), pCurrentRiser, nameStr2 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
    }
// Now the hanger it self 
    
    const double dHHanger = rTmp2 + 0.010*CLHEP::mm + 0.5*CLHEP::mm + heightRiser + heightH/2.;
    posTmp[0] = 0.;  posTmp[1] = dHHanger;
    if (iRot != 0) {
      posTmp[0] = dHHanger*rMatPr->xy();
      posTmp[1] = dHHanger*rMatPr->yy();
    }
  // Same Z position as above... 
    if (iRot == 0) { 
       new G4PVPlacement(rMatPr, posTmp, pCurrentHanger, nameStr3 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);    
     } else {
       new G4PVPlacement(G4Transform3D(*rMatPr, posTmp), pCurrentHanger, nameStr3 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
     } 
 } // on the 120 degree symmetry point.
}
//
// The set of equations that determine the radii of the inner conductor sections
//
//
// Horn2 code starts here... 
void LBNEVolumePlacements::DeclareHorn2Dims() {
  
  const double in = 2.54*CLHEP::cm;
  fHorn2RadialRescale = 1.0;
  fHorn2RadialRescaleCst = 0.0;
  fHorn2LongRescale = 1.0;
  fHorn2RadialRescaleCnt = 0;
  fHorn2LongRescaleCnt = 0;                             
  
  fHorn2PartsLengths.resize(7);
  fHorn2PartsRadii.resize(fHorn2PartsLengths.size());
  
  fHorn2OuterTubeOuterRadMax =  37.0*in/2.0; // Drawing 8875.112-MD 363385  Need an other drawing. 
  fHorn2OuterTubeOuterRad =  31.30*in/2.0; // Drawing 8875.112-MD 363385  Need an other drawing. 
  fHorn2OuterTubeInnerRad =  29.25*in/2.0; 
  fHorn2LongPosition = 6600.0*CLHEP::mm; // Drawing received fropm Alberto M., who received from Rich Stefanek, we think. 
                               // From MCZero (our Zero in g4lbne ref. frame ), to the z=0. of drawing 8875.112Md 363383 
                               
  for(size_t k=0; k!= fHorn2PartsRadii.size(); ++k) 
    fHorn2PartsRadii[k] = fHorn2OuterTubeOuterRad + 1.0*in;
  
  fHorn2OffsetIOTr1 = 1.889*in; // Offset between the entrance of Horn2TopLevel - margin and Z=0, drawing... 
  fHorn2PartsLengths[0] = fHorn2OffsetIOTr1 + 5.861*in; // Drawing 8875.112-MD 363385
  fHorn2LengthNominal = 142.91*in; // Drawing 8875.112-MD 363382, Revision E
  
  fHorn2PartsRadii[0] = (33.75 + 1.0)*in/2. ; // Info only, take the outer Radius 

  fHorn2PartsLengths[1] = (24.061)*in; // Drawing 8875.112-MD 363386
//  fHorn2PartsRadii[1] = (19.385 + 1.0)*in ; // Info only, take the outer Radius 
  
  fHorn2PartsLengths[2] = (20.092)*in; // Drawing 8875.112-MD 363387
//  fHorn2PartsRadii[2] = (10.70 + 1.0)*in ; // 
  
  fHorn2PartsLengths[3] = (20.061)*in; // Drawing 8875.112-MD 363388
//   fHorn2PartsRadii[3] = (13.293 + 1.0)*in ; // info only. 

  fHorn2PartsLengths[4] = (24.030)*in; // Drawing 8875.112-MD 363389
//   fHorn2PartsRadii[4] = (17.755 + 1.0)*in ; // info only. 

  fHorn2PartsLengths[5] = (24.325)*in; // Drawing 8875.112-MD 363390
//   fHorn2PartsRadii[4] = (21.341 + 1.0)*in ; // info only. 

  fHorn2PartsLengths[6] = (21.054)*in; // Drawing 8875.112-MD 363390
  fHorn2PartsRadii[6] = fHorn2OuterTubeOuterRadMax+ 1.0*in ; // info only. 
// add 1/2 inch margin for misalignment.
 
  fHorn2LengthMargin = 0.5*in;
   
  fHorn2Length = fHorn2LengthMargin; 
  for (size_t k=0; k!= fHorn2PartsLengths.size(); ++k) fHorn2Length +=fHorn2PartsLengths[k];
  fHorn2Length += fHorn2LengthMargin; 
//
//  std::cout << " Check of Horn2 lengths, Sum of indivuals parts + margin " << fHorn2Length 
//            << " nominal " << fHorn2LengthNominal << std::endl;
   
  int numSctTr = 3 + 2; // Number of tubes, otherwise. Approximate.. 
  //  Take 10 for the transition bulge, 2 for the outer flat part.
  const double radIOTrOut = 2.184*in; 
  const double radIOTrIn = 1.869*in; 
  const double radIOTrAv = 0.5*(radIOTrOut + radIOTrIn);
  const double deltaRadIOTr = radIOTrOut - radIOTrIn;
  const double thetaIOTRStart = std::abs(fHorn2OffsetIOTr1 - radIOTrOut)/radIOTrOut; // approximate.. 
  const double surfXSect = (M_PI - thetaIOTRStart)*(radIOTrOut*radIOTrOut - radIOTrIn*radIOTrIn);
  const double surfXSectPart = surfXSect/(numSctTr-2);
  const double bigRadIOTr = 12.756*in;
  // Valid only for 
  fHorn2UpstrOuterIOTransLengths.resize(numSctTr);
  fHorn2UpstrOuterIOTransRadsOne.resize(numSctTr);
  fHorn2UpstrOuterIOTransRadsTwo.resize(numSctTr);
  fHorn2UpstrOuterIOTransPositions.resize(numSctTr);
  
  fHorn2UpstrOuterIOTransLengths[1] = deltaRadIOTr - 0.010*CLHEP::mm; 
  fHorn2UpstrOuterIOTransPositions[1] = deltaRadIOTr/2.; // with respect to the start of the mother volume   
  const double heightTmp1 = surfXSectPart/fHorn2UpstrOuterIOTransLengths[1];
  fHorn2UpstrOuterIOTransRadsOne[1] = bigRadIOTr - heightTmp1/2.;
  fHorn2UpstrOuterIOTransRadsTwo[1] = bigRadIOTr + heightTmp1/2.;
  const double zMaxElemZeroUpstr = deltaRadIOTr + 0.010*CLHEP::mm;
  const double zMaxElemZeroDownstr = fHorn2OffsetIOTr1 - 0.010*CLHEP::mm;
  fHorn2UpstrOuterIOTransLengths[0] = zMaxElemZeroDownstr - zMaxElemZeroUpstr;
  fHorn2UpstrOuterIOTransPositions[0] = (zMaxElemZeroDownstr + zMaxElemZeroUpstr)/2.;
  fHorn2UpstrOuterIOTransLengths[2] = fHorn2UpstrOuterIOTransLengths[0];
  fHorn2UpstrOuterIOTransPositions[2] = fHorn2UpstrOuterIOTransPositions[0];
  const double heightTmp0 = surfXSectPart/fHorn2UpstrOuterIOTransLengths[0];
  fHorn2UpstrOuterIOTransRadsOne[0] = bigRadIOTr - heightTmp0/2. - radIOTrAv/sqrt(2.) ; 
  fHorn2UpstrOuterIOTransRadsTwo[0] = bigRadIOTr + heightTmp0/2. - radIOTrAv/sqrt(2.) ;
  fHorn2UpstrOuterIOTransRadsOne[2] = bigRadIOTr - heightTmp0/2. + radIOTrAv/sqrt(2.) ; 
  fHorn2UpstrOuterIOTransRadsTwo[2] = bigRadIOTr + heightTmp0/2. + radIOTrAv/sqrt(2.) ;
  // Top part. 
  fHorn2UpstrOuterIOTransLengths[3] = 0.45*in - 0.010*CLHEP::mm;
  fHorn2UpstrOuterIOTransPositions[3] = radIOTrOut + fHorn2UpstrOuterIOTransLengths[3]/2. + 0.005*CLHEP::mm;
  fHorn2UpstrOuterIOTransRadsOne[3] = bigRadIOTr + radIOTrIn ;
  fHorn2UpstrOuterIOTransRadsTwo[3] = bigRadIOTr + radIOTrOut;
  // Fange part 
  fHorn2UpstrOuterIOTransLengths[4] = 1.30*in - 0.010*CLHEP::mm;
  fHorn2UpstrOuterIOTransPositions[4] = radIOTrOut +  0.0075*CLHEP::mm + 
                                       fHorn2UpstrOuterIOTransLengths[3] + fHorn2UpstrOuterIOTransLengths[4]/2.;
  fHorn2UpstrOuterIOTransRadsOne[4] = bigRadIOTr + fHorn2UpstrOuterIOTransLengths[3]/2. + 0.005*CLHEP::mm ;
  fHorn2UpstrOuterIOTransRadsTwo[4] = 33.75*in/2;


  fHorn2InnerIOTransLength = 5.8610*in;
  
 // 
 // Load the equations Drawing 8875.112-MD 36385 - 363092
//
  fHorn2Equations.clear();
  LBNEHornRadialEquation e1(114.73006, -2.91414, 0.); fHorn2Equations.push_back(e1);
  LBNEHornRadialEquation e2(22.68402, -0.54203, 0., false); fHorn2Equations.push_back(e2);
  LBNEHornRadialEquation e3(-10.63139, 0.30058, 0., false); fHorn2Equations.push_back(e3); // linear. 
  LBNEHornRadialEquation e4(-56.92263, 1.44583, 0.); fHorn2Equations.push_back(e4); 
  LBNEHornRadialEquation e5(0.12835, 0.0932, 0., false); fHorn2Equations.push_back(e5); 
  LBNEHornRadialEquation e6(1.93227, 0.07398, 0., false); fHorn2Equations.push_back(e6); // 363390
  LBNEHornRadialEquation e7(114.73007, -2.91414, -0.11811);  fHorn2Equations.push_back(e7); // equation 7 found 0n 363386
                                                                        // equation 7 on 363385 differs sliggtly.. in the noise... 
  LBNEHornRadialEquation e8(-56.92263, 1.44583, -0.11811); fHorn2Equations.push_back(e8); 
  LBNEHornRadialEquation e9(0.01064, 0.0932, 0., false); fHorn2Equations.push_back(e9); 
  LBNEHornRadialEquation e10(1.81456, 0.07398, 0., false); fHorn2Equations.push_back(e10); // 363389
  LBNEHornRadialEquation eIO(10.9108, -0.1780, 0., false); fHorn2Equations.push_back(eIO); // 363385
 // The 2nd top lelvel section is the remaining part of that part, up to the end of the Horn1Hall Z location. 
 // This is just a tube. 
 //
    
 // Then, conical sections, all in a tube that has an inner radius of the neck, our radius 
 // 
 // To compute all that, we need to locate the target first.. 
 // So, done at construction time. 
 //
  
}
void LBNEVolumePlacements::RescaleHorn2Lengthwise() {
  if (fHorn2LongRescaleCnt != 0) {
   G4Exception("LBNEVolumePlacements::RescaleHorn2Lengthwise", "", 
                FatalErrorInArgument, " This method can only be called once per job ");
  }
  fHorn2LongRescaleCnt++;
  for (size_t i=0; i!= fHorn2PartsLengths.size(); i++) 
       fHorn2PartsLengths[i] *= fHorn2LongRescale;       
  fHorn2Length *=  fHorn2LongRescale;
  for (size_t i=0; i!= fHorn2UpstrOuterIOTransLengths.size(); i++) 
      fHorn2UpstrOuterIOTransLengths[i] *= fHorn2LongRescale;
  for (size_t i=0; i!= fHorn2UpstrOuterIOTransPositions.size(); i++) 
      fHorn2UpstrOuterIOTransPositions[i] *= fHorn2LongRescale;
  fHorn2InnerIOTransLength *= fHorn2LongRescale;
  fHorn2LengthNominal *=  fHorn2LongRescale;      
  // Now we deal with the inner conductor radial equations.
  for (std::vector<LBNEHornRadialEquation>::iterator it=fHorn2Equations.begin();
        it != fHorn2Equations.end(); it++) it->SetLongRescaleFactor(fHorn2LongRescale);
        
  fHorn2NeckLength *= fHorn2LongRescale;
  fHorn2NeckZPosition *= fHorn2LongRescale;
  // We also have to move the decay pipe snout, as Horn2 is longer.. 
  fDecayPipeWindowZLocation = fHorn2LongPosition + (fHorn2LongRescale-1.0)*fHorn2Length + 201.484*25.4*CLHEP::mm; 

  
}
  
void LBNEVolumePlacements::RescaleHorn2Radially() {
// If the user really want it call twice, so be it.  Amit, Paul, March 27 2015... 
//  if (fHorn2RadialRescaleCnt != 0) {
//   G4Exception("LBNEVolumePlacements::RescaleHorn2Radially", "", 
//                FatalErrorInArgument, " This method can only be called once per job ");
//  }
  fHorn2RadialRescaleCnt++;
//    std::cerr << " LBNEVolumePlacements::RescaleHorn2Radially, fHorn2RadialRescale = " << fHorn2RadialRescale 
//              << "fHorn2RadialRescaleCst " << fHorn2RadialRescaleCst << std::endl;
    for (size_t i=0; i!= fHorn2PartsRadii.size(); i++) {
       fHorn2PartsRadii[i] *= fHorn2RadialRescale; 
    }      
    fHorn2OuterTubeOuterRad *= fHorn2RadialRescale;
    fHorn2OuterTubeOuterRadMax *= fHorn2RadialRescale;
    fHorn2OuterTubeInnerRad *= fHorn2RadialRescale;
    for (size_t i=0; i!= fHorn2UpstrOuterIOTransRadsOne.size(); i++) {  
      fHorn2UpstrOuterIOTransRadsOne[i] *= fHorn2RadialRescale;
    }
    for (size_t i=0; i!= fHorn2UpstrOuterIOTransRadsTwo.size(); i++) { 
      fHorn2UpstrOuterIOTransRadsTwo[i] *= fHorn2RadialRescale;
    }
    
  // Now we deal with the inner conductor radial equations.



  for (std::vector<LBNEHornRadialEquation>::iterator it=fHorn2Equations.begin();
        it != fHorn2Equations.end(); it++) it->SetRadialRescaleFactor(fHorn2RadialRescale);

}

void LBNEVolumePlacements::ShiftHorn2Radially() {

//    std::cerr << " LBNEVolumePlacements::RescaleHorn2Radially, fHorn2RadialRescale = " << fHorn2RadialRescale 
//              << "fHorn2RadialRescaleCst " << fHorn2RadialRescaleCst << std::endl;
    for (size_t i=0; i!= fHorn2PartsRadii.size(); i++) {
       fHorn2PartsRadii[i] += fHorn2RadialRescaleCst;
//       std::cerr << " LBNEVolumePlacements::RescaleHorn2Radially Old Rad " << oldRad 
//                 << " New " << fHorn2PartsRadii[i] << std::endl;
    }      
    fHorn2OuterTubeOuterRad += fHorn2RadialRescaleCst;
    fHorn2OuterTubeOuterRadMax += fHorn2RadialRescaleCst;
    fHorn2OuterTubeInnerRad += fHorn2RadialRescaleCst;
    for (size_t i=0; i!= fHorn2UpstrOuterIOTransRadsOne.size(); i++) {  
      fHorn2UpstrOuterIOTransRadsOne[i] += fHorn2RadialRescaleCst;
    }
    for (size_t i=0; i!= fHorn2UpstrOuterIOTransRadsTwo.size(); i++) { 
      fHorn2UpstrOuterIOTransRadsTwo[i] += fHorn2RadialRescaleCst;
    }
    
  // Now we deal with the inner conductor radial equations.

  for (std::vector<LBNEHornRadialEquation>::iterator it=fHorn2Equations.begin();
        it != fHorn2Equations.end(); it++) it->SetRadialRescaleFactorCst(fHorn2RadialRescaleCst);
}

void LBNEVolumePlacements::PlaceFinalHorn2(G4PVPlacement *vH2Hall) {

  const double in = 2.54*CLHEP::cm;
  LBNEVolumePlacementData *plInfoH2Top = this->Create("Horn2TopLevel");
  G4PVPlacement *vH2 = PlaceFinal("Horn2TopLevel", vH2Hall);
  std::string gasMaterial = fFillHornsWithArgon ? std::string("Argon") : std::string("Air");
/* The Plug code */
  if(fConstructPlug){
    G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.; 
    G4Tubs* plugTube = new G4Tubs("Plug", fPlugInnerRadius, fPlugOuterRadius, fPlugLength*.5, 0., 360.*CLHEP::deg);
    G4LogicalVolume *plugl = new G4LogicalVolume(plugTube, G4Material::GetMaterial(fPlugMaterial.c_str()), "Plug");
    posTmp[2] = (plInfoH2Top->fParams[2]-fPlugLength)*.5-fPlugZPosition;
    new G4PVPlacement((G4RotationMatrix *) 0, posTmp, plugl, "Plug_P",
                      vH2->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
  }
  /* end of The Plug code */
  double zPosPart = -1.0*plInfoH2Top->fParams[2]/2 + fHorn2LengthMargin;
  
  
  const double zShiftDrawingIOTrans = fHorn2PartsLengths[0]/2 - fHorn2InnerIOTransLength;
  fHorn2DeltaZEntranceToZOrigin = zShiftDrawingIOTrans + fHorn2LengthMargin;
  
  double zStartDrawing =  0.; // definition of z=0 
//  std::cerr << " Horn2 Placements, zPostPart " << zPosPart << " Length of Horn2 Top level " << plInfoH2Top->fParams[2] 
//            << " zShiftDrawingIOTrans " << zShiftDrawingIOTrans << " 1rst part length " << fHorn2PartsLengths[0] << std::endl;
  double H2InnerRadius=fPlugOuterRadius;
  if(!fConstructPlug)H2InnerRadius=0.;  
  fHorn2ZEqnChanges.clear();
  for (size_t kPart=0; kPart!= fHorn2PartsLengths.size(); kPart++) {
    std::ostringstream nStrStr; nStrStr << "Horn2Part" << kPart;
    G4String nStr(nStrStr.str());
    G4Tubs* tubsPart = new G4Tubs(nStr, H2InnerRadius, fHorn2PartsRadii[kPart], 
                                  (fHorn2PartsLengths[kPart]/2 - 0.002*CLHEP::mm ), 0., 360.0*CLHEP::degree);
    G4LogicalVolume *tubsL = new G4LogicalVolume(tubsPart, G4Material::GetMaterial(gasMaterial.c_str()), nStr);
    G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.; 
    zPosPart += fHorn2PartsLengths[kPart]/2;
    posTmp[2] = zPosPart;
    std::cerr << " At Horn2 big part " << kPart << "Radius  " << fHorn2PartsRadii[kPart] << std::endl;
    G4PVPlacement *vPart = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsL, nStr + std::string("_P"), 
                     vH2->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
    zPosPart += fHorn2PartsLengths[kPart]/2;
    // Sub volumes for each part. 
    switch (kPart) {
      case 0: // The Inner Outer Transition piece 
      {
        //The inner part. 
        const double zEndDrawing = fHorn2InnerIOTransLength;
        fHorn2ZEqnChanges.push_back(zEndDrawing); // To be used in the magnetic field class.
        int numSubSect = GetNumberOfInnerHorn2SubSections(6, zStartDrawing, 
                                                      zEndDrawing, 10); // These Z position are from the start of the inner conductor.   
        const double deltaZ = (zEndDrawing - zStartDrawing)/numSubSect;
//        std::cerr << " Number of subsection for the IO transition of Horn2 " << numSubSect 
//                  << " deltaz " << deltaZ << std::endl;
       for (int iSub = 0; iSub != numSubSect; iSub++) {                                       
         const double zzBegin = zStartDrawing + iSub*deltaZ;
         const double zzEnd = zzBegin + deltaZ;
         std::ostringstream nameStrStr; nameStrStr << "Horn2InnerPartIOTransSub" << iSub;
         G4String nameStr(nameStrStr.str());
         fHorn2IC.push_back(nameStr);
         const double rMin1 = fHorn2Equations[6].GetVal(zzBegin); // Equation 1
         const double rMin2 = fHorn2Equations[6].GetVal(zzEnd);
         const double rMax1 = fHorn2Equations[10].GetVal(zzBegin) + fWaterLayerThickInHorns + 0.0025; 
         const double rMax2 = fHorn2Equations[10].GetVal(zzEnd) + fWaterLayerThickInHorns + 0.0025;     
         G4Cons *aCons = new G4Cons(nameStr, rMin1, rMax1,rMin2, rMax2,
                                      (deltaZ - 0.005*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
       


         G4LogicalVolume *pCurrent = new G4LogicalVolume(aCons, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
         posTmp[0] = 0.; posTmp[1] = 0.;
         posTmp[2] = zzBegin  + zShiftDrawingIOTrans + deltaZ/2.;                             
         G4PVPlacement *vSub = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrent, nameStr + std::string("_P"), 
                        vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
                        
         if (fWaterLayerThickInHorns > 0.002*CLHEP::mm) {
           nameStrStr.str(""); nameStrStr.clear(); nameStrStr << "Horn2InnerPartIOTransSub" << iSub << "Water";
           nameStr = nameStrStr.str();
           G4Cons *aConsW2 = new G4Cons(nameStr, rMax1 - fWaterLayerThickInHorns, rMax1-0.001*CLHEP::mm,
                                         rMax2 - fWaterLayerThickInHorns, rMax2-0.001*CLHEP::mm,
                                      (deltaZ - 0.0075*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
           G4LogicalVolume *pCurrentW2 = new G4LogicalVolume(aConsW2, G4Material::GetMaterial(std::string("Water")), nameStr);
           posTmp[0] = 0.; posTmp[1] = 0.; posTmp[2] =0.;                             
           new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrentW2, nameStr + std::string("_P"), 
                        vSub->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);     
         }
       } // on the number of sub section of the Inner Outer transition 
       // Outer parts 
       for (size_t kOut=0; kOut != fHorn2UpstrOuterIOTransLengths.size(); kOut++) {
          std::ostringstream nStrStrTR3; nStrStrTR3 << "Horn2IOTransOuterPart" << kOut;
          G4String nStrTR3(nStrStrTR3.str());
          G4Tubs* tubsPartTR3 = new G4Tubs(nStr, fHorn2UpstrOuterIOTransRadsOne[kOut], 
                                              fHorn2UpstrOuterIOTransRadsTwo[kOut], 
                                              fHorn2UpstrOuterIOTransLengths[kOut]/2., 0., 360.0*CLHEP::degree);
          G4LogicalVolume *tubsLTR3 = new G4LogicalVolume(tubsPartTR3, G4Material::GetMaterial(fHorn2AllCondMat), nStr);
          posTmp[0] =0.; posTmp[1] =0.; 
          posTmp[2] = - fHorn2PartsLengths[kPart]/2 + fHorn2UpstrOuterIOTransPositions[kOut];
          new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLTR3, nStrTR3 + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
       }
       // The weld at the end. Shift it a bit upstream such that it does not interfere 
       {
          G4String nStrTW4("Horn2IOTransWeld");
          fHorn2IC.push_back(nStrTW4);
          const double weldLength = 0.386*in;
          const double zzEnd = fHorn2LongRescale*5.8610*in - weldLength - 0.050*CLHEP::mm;
          const double radius = fHorn2Equations[10].GetVal(zzEnd) + fWaterLayerThickInHorns + 0.0025; 
          const double thick = 0.0775*in; 
          G4Tubs* tubsPartTW4 = new G4Tubs(nStrTW4, radius, (radius+thick), weldLength/2., 0., 360.0*CLHEP::degree);
          G4LogicalVolume *tubsLTW4 = new G4LogicalVolume(tubsPartTW4, G4Material::GetMaterial(fHorn2InnerCondMat), nStr);
          posTmp[0] = 0.; posTmp[1] =0.; 
          posTmp[2] =  fHorn2PartsLengths[kPart]/2 - weldLength/2. - 0.025*CLHEP::mm;
          new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLTW4, nStrTW4 + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
        }
        zStartDrawing += fHorn2InnerIOTransLength; // for the next section 
      } // End of Transition piece 
      break;
      case 2: // The neck.  Must subdivide further, radial equation keeps changing.
       // There will be 5 parts .. Drawing 8875.112-MD 363387
       
      {
        const double lengthNeck = fHorn2LongRescale*2.411*in - 0.020*CLHEP::mm; // approximate. Equation not available in usual form 
        fHorn2NeckLength = lengthNeck; // to compute the magnetic field 
        const double zPosNeck = fHorn2LongRescale*(39.8193 - 29.800)*in; // with respect to the start of the mother volume 
        fHorn2NeckZPosition = zPosNeck;
        // We start by the neck 
        {
          G4String nStrPN2("Horn2Part2Neck");
          const double radiusInner =  fHorn2RadialRescale*3.071*in/2. + fHorn2RadialRescaleCst;
          const double radiusOuter = fHorn2RadialRescale*3.465*in/2. + fWaterLayerThickInHorns + 0.0025 + fHorn2RadialRescaleCst; 
          fHorn2NeckOuterRadius = fHorn2RadialRescale*3.465*in/2.+fHorn2RadialRescaleCst;
          fHorn2NeckInnerRadius = radiusInner;
          G4Tubs* tubsPartPN2 = new G4Tubs(nStrPN2, radiusInner, radiusOuter, lengthNeck/2., 0., 360.0*CLHEP::degree );
          G4LogicalVolume *tubsLPN2 = new G4LogicalVolume(tubsPartPN2, G4Material::GetMaterial(fHorn2InnerCondMat), nStr);
          posTmp[0] = 0.; posTmp[1] =0.; 
          posTmp[2] =  -1.0*fHorn2PartsLengths[kPart]/2  + zPosNeck;
          fHorn2IC.push_back(nStrPN2);
          G4PVPlacement *vNeck = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLPN2, nStrPN2 + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
         if (fWaterLayerThickInHorns > 0.002*CLHEP::mm) {
           G4String nameStrW5(nStr); nameStrW5 += std::string("Water");
           G4Tubs* tubsPartW5 = new G4Tubs(nameStrW5, radiusOuter -0.001*CLHEP::mm - fWaterLayerThickInHorns , 
                                          radiusOuter - 0.001*CLHEP::mm , 
                                          lengthNeck/2. - 0.002*CLHEP::mm, 0., 360.0*CLHEP::degree);
           G4LogicalVolume *pCurrentW5 = new G4LogicalVolume(tubsPartW5, G4Material::GetMaterial(std::string("Water")), nameStrW5);
           posTmp[0] = 0.; posTmp[1] = 0.; posTmp[2] =0.;                             
           new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrentW5, nameStrW5 + std::string("_P"), 
                        vNeck->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
         }
        //    std::cerr<< "Rescaled Horn2 Neck Radius is, outer " << fHorn2NeckOuterRadius << "Inner "<<fHorn2NeckInnerRadius<<std::endl;
             
        } // Neck part per se.  
        {
          double zStartDrawingsSubP2[] = {29.800, 36.8888, 41.0248, 43.3660}; // Note the third number is compromise, corresponding 
                                                                              // to the length of the neck going ~2 " to 2.41 "    
          double zEndDrawingsSubP2[] = {36.8888, 38.6138, 43.3660, 49.820};
          size_t eqnInner[] = {6, 6, 7, 7};
          size_t eqnOuter[] = {0, 1, 2 ,3}; 
         // Now loop on these 4 subparts.  (this gona a bit tedious...) 
          for (size_t kSub=0; kSub!=4; kSub++) {
            const double zStartDrawingPP2 = fHorn2LongRescale*zStartDrawingsSubP2[kSub]*in;
            const double zEndDrawingPP2 =  fHorn2LongRescale*zEndDrawingsSubP2[kSub]*in;
            fHorn2ZEqnChanges.push_back(zStartDrawingPP2);
            // Need to insert the neck one.. to get consistent set to compute the magnetic field. 
            if (kSub == 1) fHorn2ZEqnChanges.push_back(zEndDrawingPP2);
            const double zShiftTmp = zStartDrawingPP2 - fHorn2LongRescale*zStartDrawingsSubP2[0]*in;
            int numSubSect = GetNumberOfInnerHorn2SubSections(eqnInner[kSub], zStartDrawingPP2, zEndDrawingPP2, 10);   
            const double deltaZ = (zEndDrawingPP2 - zStartDrawingPP2)/numSubSect;
//            std::cerr << " Number of subsection for the middle part of Horn2 index " << kPart 
//                << " Upstream of the neck,  is " << numSubSect 
//                  << " deltaz " << deltaZ << std::endl;
            for (int iSub = 0; iSub != numSubSect; iSub++) {                                          
              const double zzBegin = zStartDrawingPP2 + iSub*deltaZ;
              const double zzEnd = zzBegin + deltaZ;
              std::ostringstream nameStrStr; 
              nameStrStr << "Horn2InnerPart" << kPart << "NeckkSub" << kSub << "iSub" << iSub;
              G4String nameStr(nameStrStr.str());
              fHorn2IC.push_back(nameStr);
              const double rMin1 = fHorn2Equations[eqnInner[kSub]].GetVal(zzBegin); // Equation 1
              const double rMin2 = fHorn2Equations[eqnInner[kSub]].GetVal(zzEnd);
              const double rMax1 = fHorn2Equations[eqnOuter[kSub]].GetVal(zzBegin) + fWaterLayerThickInHorns + 0.0025; 
              const double rMax2 = fHorn2Equations[eqnOuter[kSub]].GetVal(zzEnd) + fWaterLayerThickInHorns + 0.0025;     
//              std::cerr << " Pre/Post Neck, zzbegin " << zzBegin << " radii " << 
//                       rMin1 << " / " << rMax1 << " / " << rMin2 <<" / " << rMax2 << std::endl;
              G4Cons *aCons = new G4Cons(nameStr, rMin1, rMax1,rMin2, rMax2,
                                      (deltaZ - 0.005*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
                      
              G4LogicalVolume *pCurrent = new G4LogicalVolume(aCons, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
              posTmp[0] = 0.; posTmp[1] = 0.;
              posTmp[2] = -fHorn2PartsLengths[kPart]/2 + iSub*deltaZ + deltaZ/2. + zShiftTmp;                         
              G4PVPlacement *vSub = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrent, nameStr + std::string("_P"), 
                        vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
                        
             if (fWaterLayerThickInHorns > 0.002*CLHEP::mm) {
               G4String nameStrW6(nameStr); nameStrW6 += std::string("Water");
               G4Cons *aConsW6 = new G4Cons(nameStrW6, rMax1 - fWaterLayerThickInHorns, rMax1-0.001*CLHEP::mm,
                                         rMax2 - fWaterLayerThickInHorns, rMax2-0.001*CLHEP::mm,
                                      (deltaZ - 0.0075*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
               G4LogicalVolume *pCurrentW6 = new G4LogicalVolume(aConsW6, G4Material::GetMaterial(std::string("Water")), nameStrW6);
               posTmp[0] = 0.; posTmp[1] = 0.; posTmp[2] =0.;                         
               new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrentW6, nameStrW6 + std::string("_P"), 
                        vSub->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);     
             }
        }
       } // Parts surrounding the neck. 
       // Weld Only the downstream one. The upstream part is comprise in the previous volume 
//       if (std::abs(fHorn2RadialRescaleCst) < 1.0e-10*CLHEP::mm) {
       {
          std::ostringstream nameStrStr; nameStrStr << "Horn2NeckPart" << kPart << "WeldDownstr";
          G4String nameStr(nameStrStr.str());
          fHorn2IC.push_back(nameStr);
          const double weldLength = 2.0 * 0.386*in; // double the length to include the part on the next volume 
          const double zzEndForWeld = fHorn2LongRescale*zEndDrawingsSubP2[3]*in + 0.50*CLHEP::mm; 
          const double radius = fHorn2Equations[3].GetVal(zzEndForWeld) + fWaterLayerThickInHorns + 0.025; // Using equation e8, drawing 363387
//        std::cerr << " Debugging bad Radial Cst.. zzEndForWeld " <<  zzEndForWeld << " radius " << radius 
//                  << " kPart " << kPart << " name " << nameStr << std::endl;
          const double thick = 0.0775*in; 
          G4Tubs* tubsPartN2P2 = new G4Tubs(nameStr, radius, (radius+thick), weldLength/2., 0., 360.0*CLHEP::degree);
          G4LogicalVolume *tubsLN2P2 = new G4LogicalVolume(tubsPartN2P2, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
          posTmp[0] = 0.; posTmp[1] =0.; 
          posTmp[2] =  fHorn2PartsLengths[kPart]/2 - weldLength - 0.025*CLHEP::mm;
          new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLN2P2, nameStr + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
        }
        zStartDrawing =  fHorn2LongRescale*49.820*in;
      } 
      break;
     } // end of case 2 
     case 6: // the downstream end piece. simple tubes. 
      {
        {
         G4String nameStr("Horn2DownstrPart6");
         const double length = fHorn2LongRescale*(21.054 - 1.25)*in - 0.020*CLHEP::mm;
         const double radInner = fHorn2RadialRescale*21.103*in/2. + fHorn2RadialRescaleCst;
         const double radOuter = fHorn2RadialRescale*21.853*in/2. + fHorn2RadialRescaleCst;
         G4Tubs* tubsDP6 = new G4Tubs(nameStr, radInner, radOuter, length/2., 0., 360.0*CLHEP::degree);
         G4LogicalVolume *tubsLDP6 = new G4LogicalVolume(tubsDP6, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
         posTmp[0] = 0.; posTmp[1] =0.; 
         posTmp[2] =  -fHorn2PartsLengths[kPart]/2 + length/2. + 0.0125*CLHEP::mm;
         fHorn2IC.push_back(nameStr);
         new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLDP6, nameStr + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
        }
         // the flange/connector 
          
         {
         G4String nameStr("Horn2DownstrPart6Flange");
         const double length = fHorn2LongRescale*(1.25)*in - 0.020*CLHEP::mm;
         const double radInner = fHorn2RadialRescale*21.103*in/2. + fHorn2RadialRescaleCst;
         const double radOuter = fHorn2RadialRescale*25.040*in/2. + fHorn2RadialRescaleCst;
         G4Tubs* tubsFL6 = new G4Tubs(nameStr, radInner, radOuter, length/2., 0., 360.0*CLHEP::degree);
         G4LogicalVolume *tubsLFL6 = new G4LogicalVolume(tubsFL6, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
         posTmp[0] = 0.; posTmp[1] =0.; 
         posTmp[2] =  fHorn2PartsLengths[kPart]/2 - length/2. - 0.0125*CLHEP::mm;
         fHorn2IC.push_back(nameStr);
         new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLFL6, nameStr + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
                     
        }
        break;     
      }
      default: 
      { // One set of conical sections. + weld at the end
        const double zEndDrawing = zStartDrawing + fHorn2PartsLengths[kPart];
        size_t eqnInner = 6; // kPart =1 
        size_t eqnOuter = 0;
        if (kPart == 3 ) {
          eqnOuter = 3;
          fHorn2ZEqnChanges.push_back(zStartDrawing);
          eqnInner = 7;
        } else if (kPart == 4) {
          eqnOuter = 4;
          fHorn2ZEqnChanges.push_back(zStartDrawing);
          eqnInner = 8;
        } else if (kPart == 5) {
          eqnOuter = 5;
          fHorn2ZEqnChanges.push_back(zStartDrawing);
          eqnInner = 9;
          fHorn2ZEqnChanges.push_back(zEndDrawing);
        }
        int numSubSect = GetNumberOfInnerHorn2SubSections(eqnOuter, zStartDrawing, 
                                                      zEndDrawing, 10);   
        const double deltaZ = (zEndDrawing - zStartDrawing)/numSubSect;
//        std::cerr << " Number of subsection for the generic part of Horn2 index " << kPart << " is " << numSubSect 
//                  << " deltaz " << deltaZ << std::endl;
        double zzEndForWeld = 0.;         
       for (int iSub = 0; iSub != numSubSect; iSub++) {                                       
         const double zzBegin = zStartDrawing + iSub*deltaZ;
         const double zzEnd = zzBegin + deltaZ;
         zzEndForWeld = zzEnd;
         std::ostringstream nameStrStr; nameStrStr << "Horn2InnerPart" << kPart << "Sub" << iSub;
         G4String nameStr(nameStrStr.str());
         fHorn2IC.push_back(nameStr);
         const double rMin1 = fHorn2Equations[eqnInner].GetVal(zzBegin); // Equation 1
         const double rMin2 = fHorn2Equations[eqnInner].GetVal(zzEnd);
         const double rMax1 = fHorn2Equations[eqnOuter].GetVal(zzBegin) + fWaterLayerThickInHorns + 0.0025; 
         const double rMax2 = fHorn2Equations[eqnOuter].GetVal(zzEnd) + fWaterLayerThickInHorns + 0.0025;     
         G4Cons *aCons = new G4Cons(nameStr, rMin1, rMax1,rMin2, rMax2,
                                      (deltaZ - 0.005*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
                                                              
         G4LogicalVolume *pCurrent = new G4LogicalVolume(aCons, 
                                    G4Material::GetMaterial(std::string(fHorn2InnerCondMat)), nameStr);
         posTmp[0] = 0.; posTmp[1] = 0.;
         posTmp[2] = -fHorn2PartsLengths[kPart]/2 + iSub*deltaZ + deltaZ/2.;                          
         G4PVPlacement *vSub = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrent, nameStr + std::string("_P"), 
                        vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
                        
         if (fWaterLayerThickInHorns > 0.002*CLHEP::mm) {
           G4String nameStrW(nameStr); nameStrW += std::string("Water");
           G4Cons *aConsW7 = new G4Cons(nameStrW, rMax1 - fWaterLayerThickInHorns, rMax1-0.001*CLHEP::mm,
                                         rMax2 - fWaterLayerThickInHorns, rMax2-0.001*CLHEP::mm,
                                      (deltaZ - 0.0075*CLHEP::mm)/2., 0., 360.0*CLHEP::degree);
                                      
           G4LogicalVolume *pCurrentW7 = new G4LogicalVolume(aConsW7, G4Material::GetMaterial(std::string("Water")), nameStrW);
           posTmp[0] = 0.; posTmp[1] = 0.; posTmp[2] =0.;                             
           new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrentW7, nameStrW + std::string("_P"), 
                        vSub->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);     
         }
        
      }
        // The weld at the end. Shift it a bit upstream such that it does not interfere 
       {
          std::ostringstream nameStrStr; nameStrStr << "Horn2InnerPart" << kPart << "Weld";
          G4String nameStr(nameStrStr.str());
          fHorn2IC.push_back(nameStr);
          const double weldLength = 2.0 * 0.386*in; // double to include the next one.. Perhaps a bit oversized. 
          if (kPart == 1) zzEndForWeld -= weldLength - 0.50*CLHEP::mm; // oversizing a bit by moving while computing the radius
          else zzEndForWeld += 0.50*CLHEP::mm; 
          const double radius = fHorn2Equations[eqnOuter].GetVal(zzEndForWeld) + fWaterLayerThickInHorns + 0.0025; 
          const double thick = 0.0775*in; 
//        std::cerr << " Debugging bad Radial Cst.. zzEndForWeld " <<  zzEndForWeld << " radius " << radius 
//                  << " kPart " << kPart << " name " << nameStr << std::endl;
          G4Tubs* tubsPartWe7 = new G4Tubs(nameStr, radius, (radius+thick), weldLength/2., 0., 360.0*CLHEP::degree);
          G4LogicalVolume *tubsLWe7 = new G4LogicalVolume(tubsPartWe7, G4Material::GetMaterial(fHorn2InnerCondMat), nameStr);
          posTmp[0] = 0.; posTmp[1] =0.; 
          posTmp[2] =  fHorn2PartsLengths[kPart]/2 - weldLength/2. - 0.025*CLHEP::mm;
          new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLWe7, nameStr + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
        }
     
      zStartDrawing += fHorn2PartsLengths[kPart]; // for the next section 
     }
    } //end of special cases. 
    // Outer tube (not for k=0 ) 
    if (kPart != 0) {
      std::ostringstream nStrStrOU; nStrStrOU << "Horn2OuterConductor" << kPart;
      G4String nStrOU(nStrStrOU.str());
      const double leffOU = (fHorn2PartsLengths[kPart] - 0.040*CLHEP::mm);
      G4Tubs* tubsPartOU = new G4Tubs(nStrOU, fHorn2OuterTubeInnerRad, fHorn2OuterTubeOuterRad,
         leffOU/2.0, 0., 360.0*CLHEP::degree);
      G4LogicalVolume *tubsLOU = new G4LogicalVolume(tubsPartOU, G4Material::GetMaterial(fHorn2AllCondMat), nStrOU);
      posTmp[0] =0.; posTmp[1] =0.; posTmp[2] =0.; 
      new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLOU, nStrOU + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
    } else {
    
        // April 2 2014: per careful observations by Seongtae Park 
      // Talk delivered on March 27 2014. 
      // The I/O transition outer conductor does not completly covers the range. 
       std::ostringstream nStrStrOU2; nStrStrOU2 << "Horn2OuterConductor" << kPart << "Extra";
      G4String nStrOU2(nStrStr.str());
      const double leffOU2 = fHorn2LongRescale*96.5*CLHEP::mm;
      G4Tubs* tubsPartOU2 = new G4Tubs(nStrOU2, fHorn2OuterTubeInnerRad, fHorn2OuterTubeOuterRad,
         leffOU2/2.0, 0., 360.0*CLHEP::degree);
      G4LogicalVolume *tubsLOU2 = new G4LogicalVolume(tubsPartOU2, G4Material::GetMaterial(fHorn2AllCondMat), nStrOU2);
      posTmp[0] =0.; posTmp[1] =0.;
      posTmp[2] = fHorn2PartsLengths[0]/2.0 - leffOU2/2.0 - 0.010*CLHEP::mm;
      new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsLOU2, nStrOU2 + std::string("_P"), 
                     vPart->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC); 
 
    }
    // Eventual Hangers. only one set 
    if (kPart == 3) this->Horn2InstallSpiderHanger(nStr, vPart);
    
    // More 
  } // Loop over the parts, daughters of top level. 
  // Final flange downstream flange Approximate Drawing 363382-E1  Length is still a guess.. 
  G4String nStr("Horn2DownstrOuterFlange");
  G4Tubs* tubsOF1 = new G4Tubs(nStr, fHorn2OuterTubeOuterRad + 1.0*CLHEP::mm, 
                                fHorn2OuterTubeOuterRadMax-1.0*CLHEP::mm, 0.50*in, 0., 360.0*CLHEP::degree);
  G4LogicalVolume *tubsL = new G4LogicalVolume(tubsOF1, G4Material::GetMaterial(fHorn2AllCondMat), nStr);
  G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.;
  posTmp[2] = plInfoH2Top->fParams[2]/2 - fHorn2LongRescale*0.6*in; // approximate location
  new G4PVPlacement((G4RotationMatrix *) 0, posTmp, tubsL, nStr + std::string("_P"), 
                     vH2->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);   
}

void LBNEVolumePlacements::PlaceFinalDecayPipeSnout(G4PVPlacement* tunnel) {

  // ref : Drawing number 2251-000-ME-487107 and phone converstion with Glenn Waver (MD/AD) 
  const double in = 2.54*CLHEP::cm;
  G4PVPlacement* vTop = this->PlaceFinal(G4String("DecayPipeSnout"), tunnel);
  const LBNEVolumePlacementData *plInfoTop = this->Find(G4String("Snout"), G4String("DecayPipeSnout"),
                                     G4String("LBNEVolumePlacements::PlaceFinalDecayPipeSnout"));
  //
  // The window is tilted, true vertical. Pick the rotation. Create a rotated buffer volume   
  //  
  LBNERunManager *pRunManager=static_cast<LBNERunManager*> (LBNERunManager::GetRunManager());
  const LBNEDetectorConstruction *pDet = 
            static_cast<const LBNEDetectorConstruction*> (pRunManager->GetUserDetectorConstruction());
  G4double angle = pDet->GetBeamlineAngle();
  G4RotationMatrix *beamAngleRot = new G4RotationMatrix;
  beamAngleRot->rotateX(-angle);
//  std::cerr << " Window angle, yz term " << beamAngleRot->yz() << " yy " 
//            << beamAngleRot->yy() << " zz " << beamAngleRot->zz() << std::endl; 
  double boxX = 50.0*in + 0.1*CLHEP::mm; 
  double boxY = 48.0*in + 0.1*CLHEP::mm; 
  double boxZ = 24.0*in + 0.1*CLHEP::mm;  
  G4String boxCStr("DecayPipeUpstrWindowCont");
  G4Box* aBoxC = new G4Box(boxCStr, boxX/2., boxY/2., boxZ/2.);
  G4LogicalVolume *contL = new G4LogicalVolume(aBoxC, G4Material::GetMaterial("DecayPipeGas"), boxCStr);
  G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.;
  posTmp[2] = -plInfoTop->fParams[2]/2 + boxZ/2. + 0.5*CLHEP::cm; // approximate location, but accurate enough 
  posTmp[2] += 25.0*in*std::abs( beamAngleRot->yz()); // inclined! Shift a bit upstream 
//  std::cerr << " Position of the rotated container  " << posTmp[2] << std::endl;
  G4PVPlacement *vCont = new G4PVPlacement(beamAngleRot, posTmp, contL, boxCStr + std::string("_P"), 
                     vTop->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
// if air, we stop here, the aluminum window and the window itself don't exist. 
// (3.5 million bugs saving, I am told..) 
                      
   if (fDecayPipeGas == G4String("Air")) return;
   boxX -= 0.1*CLHEP::mm;
   boxY -= 0.1*CLHEP::mm;
   boxZ -= 0.1*CLHEP::mm;
   G4String boxFStr("DecayPipeUpstrWindowFrame");
   G4Box* aBoxF = new G4Box(boxFStr, boxX/2., boxY/2., boxZ/2.);
   G4LogicalVolume *frameL = new G4LogicalVolume(aBoxF, G4Material::GetMaterial("Aluminum"), boxFStr);
   posTmp[0] =0.; posTmp[1] =0.; posTmp[2] =0.;
   G4PVPlacement *vFrame = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, 
                          frameL, boxFStr + std::string("_P"), 
                     vCont->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
                     
   G4String holeFStr("DecayPipeUpstrWindowHole");
   double holeRadius= fDecayPipeWindowRadiusAlum + 0.01*CLHEP::mm;
   G4Tubs* holeR = new G4Tubs(holeFStr, 0., holeRadius, boxZ/2., 0., 360.0*CLHEP::degree);
                                // The upstream part is Air.. 
   G4LogicalVolume *holeL = new G4LogicalVolume(holeR, G4Material::GetMaterial("Air"),holeFStr);
   G4PVPlacement *vHole = new G4PVPlacement((G4RotationMatrix *) 0, posTmp, 
                          holeL, holeFStr + std::string("_P"), 
                     vFrame->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
//
   G4String windAlumStr("DecayPipeUpstrWindowAlum");
   G4Tubs* windAlum = new G4Tubs(windAlumStr, fDecayPipeWindowRadiusBeryl, 
                             fDecayPipeWindowRadiusAlum , fDecayPipeWindowThickAlum/2., 0., 360.0*CLHEP::degree);
                                // The upstream part is Air.. 
   G4LogicalVolume *windAlumL = new G4LogicalVolume(windAlum, G4Material::GetMaterial("Aluminum"), windAlumStr);
   posTmp[2] = -boxZ/2. + 17.1*in + 3.0*CLHEP::mm; // first number picked up from drawing. We place the Alum a bit behind..  
   new G4PVPlacement((G4RotationMatrix *) 0, posTmp, 
                          windAlumL, windAlumStr + std::string("_P"), 
                     vHole->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
                     
   G4String windBerylStr("DecayPipeUpstrWindowBeryl");
   G4Tubs* windBeryl = new G4Tubs(windBerylStr, 0., 
                             fDecayPipeWindowRadiusBeryl , fDecayPipeWindowThickBeryl/2., 0., 360.0*CLHEP::degree);
                                // The upstream part is Air.. 
   G4LogicalVolume *windBerylL = new G4LogicalVolume(windBeryl, G4Material::GetMaterial("Beryllium"), windBerylStr);
   posTmp[2] = -boxZ/2. + 17.1*in + 1.0*CLHEP::mm; // first number picked up from drawing.
   new G4PVPlacement((G4RotationMatrix *) 0, posTmp, 
                          windBerylL, windBerylStr + std::string("_P"), 
                     vHole->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);
} 
// Clone for Horn1.. Could avoid a bit of code bloat, but that's not much... 
int LBNEVolumePlacements::GetNumberOfInnerHorn2SubSections(size_t eqn, double z1, double z2, int nMax) const {

   int numSubSect = 1;
   double zEnd = z2;
   const double zLengthTotal = z2 - z1;
   double zLength = zLengthTotal;
   while (true) { 
     const double r0 = fHorn2Equations[eqn].GetVal(z1);
     const double r1 = fHorn2Equations[eqn].GetVal(zEnd);
     const double rMid = fHorn2Equations[eqn].GetVal((z1 + zEnd)/2.);
//     std::cerr << " GetNumberOfInnerHorn2SubSections z1 " << z1 << " zEnd " << zEnd << 
//                     " r0 " << r0 << " r1 " << r1 << " rMid " << rMid << std::endl;
     if (std::abs(rMid - 0.5*(r0+r1)) < 0.050*CLHEP::mm) break; // 50 microns good enough for Gov. work. 
                                                         // meachanical tolerance is one 5 mills, (127 microns)
     zLength /=2.;
     zEnd = z1 + zLength;
     numSubSect++;
     if (numSubSect > nMax) {
       G4Exception("LBNEVolumePlacements::PlaceFinalHorn2", " ", FatalErrorInArgument, " Crazy subdivision! ");
       break;
     }
   }
   return numSubSect;
}
// Clone from Horn1 .. But simpler, less arguments because there is one set of such contraptions. 

void LBNEVolumePlacements::Horn2InstallSpiderHanger(const G4String &nameStrH, G4PVPlacement *vMother ) {

  const double in = 2.54*CLHEP::cm;
  G4String nameStr(nameStrH);nameStr += G4String("Ring");
  // Drawing 8875.112-MD 363388 and 363405
  const double length = 0.750*in*fHorn2LongRescale;
  const int eqnNum = 3;
  const double zLocDrawing = fHorn2LongRescale*60.4681*in;
  const double rTmp1 = fHorn2Equations[eqnNum].GetVal(zLocDrawing) + 0.015*CLHEP::mm + fWaterLayerThickInHorns;
//  const double rTmp2 = fHorn2RadialRescale*11.625*in + fWaterLayerThickInHorns + 0.015*CLHEP::mm; // 
// 11.625 inches is the outer diameter of this ring.  Not it's radius. Arghh! Feb 4 2014, P. Lebrun
  const double rTmp2 = fHorn2RadialRescale*11.625*in/2.0 + fWaterLayerThickInHorns + 0.015*CLHEP::mm + fHorn2RadialRescaleCst; // 
// Re instate the error to study it.. Feb 19 2014 
//  const double rTmp2 = fHorn2RadialRescale*11.625*in + fWaterLayerThickInHorns + 0.015*CLHEP::mm; // 
//  std::cerr <<  " LBNEVolumePlacements::Horn2InstallSpiderHanger, radii of the rings.. " << rTmp1 << "  " << rTmp2 << std::endl;
//  std::cerr <<  " .... Volume name ...  " << nameStr << " ... Aud quit !!! " << std::endl; exit(2);
  G4Tubs *aTubs = new G4Tubs(nameStr, rTmp1, rTmp2, length/2.   , 0., 360.0*CLHEP::degree);
  G4LogicalVolume *pCurrent = new G4LogicalVolume(aTubs, G4Material::GetMaterial(fHorn2AllCondMat), nameStr);
  G4ThreeVector posTmp; posTmp[0] = 0.; posTmp[1] = 0.;
  posTmp[2] = fHorn2LongRescale*in*(-0.0180); // computed from drawing 8875.112MD-363388                           
  new G4PVPlacement((G4RotationMatrix *) 0, posTmp, pCurrent, nameStr + std::string("_P"), 
                     vMother->GetLogicalVolume(), false, 1, fCheckVolumeOverLapWC);    
  
  //The connecting piece ring to the hangers.. There are three of them, at 120 degrees from each other. 
  
  G4String nameStr2(nameStrH); nameStr2 += G4String("Riser");
  const double heightRiser = 0.333*in - 0.020*CLHEP::mm;
  const double widthH = 1.5*in; // See drawing 8875.112-MD 363115
  const double thickH = 0.184*2*in; 
  G4Box *aBoxRiser = new G4Box(nameStr2, widthH/2., heightRiser/2.0, thickH/2.0);  
  G4LogicalVolume *pCurrentRiser = 
    new G4LogicalVolume(aBoxRiser, G4Material::GetMaterial(fHorn2AllCondMat), nameStr2);
    
  G4String nameStr3(nameStrH); nameStr3 += G4String("Hanger");
  const double heightH = fHorn2OuterTubeInnerRad - rTmp2 - 1.0*CLHEP::mm - heightRiser;
  const double widthH2 = 1.0*in; // 363115 Note: we collapsed both hanger along the horizontal, transverse 
                                // direction. 
  const double thickH2 = 0.031*in; 
  G4Box *aBoxHanger = new G4Box(nameStr3, widthH2/2., heightH/2.0, thickH2/2.0);  
  G4LogicalVolume *pCurrentHanger = 
    new G4LogicalVolume(aBoxHanger, G4Material::GetMaterial(fHorn2AllCondMat), nameStr3);
  
  for (int iRot=0; iRot != 3; iRot++) {
    std::ostringstream rnStrStr; rnStrStr << "_" << (iRot+1);
    G4String rnStr(rnStrStr.str());
  // Same Z position as above... 
    G4RotationMatrix * rMatPr = 0;
     if (iRot == 1) {
      rMatPr = new G4RotationMatrix; 
      rMatPr->rotateZ(2.0*M_PI/3.);
    } else if (iRot == 2) {
      rMatPr = new G4RotationMatrix; 
      rMatPr->rotateZ(-2.0*M_PI/3.);
    }
    
    const double dHRiser = rTmp2 + 0.010*CLHEP::mm + heightRiser/2.;
    posTmp[0] = 0.;  posTmp[1] = dHRiser;
    if (iRot != 0) {
      posTmp[0] = dHRiser*rMatPr->xy();
      posTmp[1] = dHRiser*rMatPr->yy();
    }
    if (iRot == 0) { 
       new G4PVPlacement(rMatPr, posTmp, pCurrentRiser, nameStr2 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
    } else {
       new G4PVPlacement(G4Transform3D(*rMatPr, posTmp), pCurrentRiser, nameStr2 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
    }
// Now the hanger it self 
    
    const double dHHanger = rTmp2 + 0.010*CLHEP::mm + 0.5*CLHEP::mm + heightRiser + heightH/2.;
    posTmp[0] = 0.;  posTmp[1] = dHHanger;
    if (iRot != 0) {
      posTmp[0] = dHHanger*rMatPr->xy();
      posTmp[1] = dHHanger*rMatPr->yy();
    }
  // Same Z position as above... 
    if (iRot == 0) { 
       new G4PVPlacement(rMatPr, posTmp, pCurrentHanger, nameStr3 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);    
     } else {
       new G4PVPlacement(G4Transform3D(*rMatPr, posTmp), pCurrentHanger, nameStr3 + std::string("_P") + rnStr, 
                     vMother->GetLogicalVolume(), false, iRot, fCheckVolumeOverLapWC);
     } 
 } // on the 120 degree symmetry point.
}
void LBNEVolumePlacements::ExtendChaseLengthForHorn2(){
  //
  // Pior to construct the long chase that contains both the target and the horns, we check the length of this mother volume
  // is sufficient. If not, we do not extend the tunnel, but reduce the decay pipe length. Paul L. Lebrun, Feb 2015. 
  // Based on discussion with Alberto M., Feb 6 2015. 
  //
  const double aH2AllLength = fHorn2Length + 4.0*fHorn2LengthMargin + 2.5*CLHEP::cm; // Last is the Amit's tracking plane
  // We assume here that the decay pipe snout will not be placed.. 
  const double aEndH2 = aH2AllLength + fHorn2LongPosition - fHorn2OffsetIOTr1; // See Placement of Horn2Hall
  if (aEndH2 > fDecayPipeLongPosition ) {
    const double delta = aEndH2 - fDecayPipeLongPosition + 1.0*CLHEP::cm; // Last if for good measure..
    fDecayPipeLongPosition += delta;
    fDecayPipeLength -= delta;
    std::cerr << " LBNEVolumePlacements::ExtendChaseLengthForHorn2, moving beginning of decay pipe downstream by " 
              << delta/CLHEP::m  << " meters " << std::endl;
    std::cerr << " New Decay length = " <<  fDecayPipeLength/CLHEP::m << " Starting at " 
              <<  fDecayPipeLongPosition/CLHEP::m << " meters"  << std::endl;
  }
  
}

double LBNEHornRadialEquation::inchDef = 2.54*CLHEP::cm;

LBNEHornRadialEquation::LBNEHornRadialEquation(double rc, double zc, double rOffTmp, bool isParabolic):
parabolic(isParabolic),
rCoeff(rc),
zCoeff(zc),
rOff(rOffTmp),
rResc(1.0),
rRescCst(0.),
zResc(1.0)
{ ; } 

double LBNEHornRadialEquation::GetVal(double z) const {
  // z = 0 above is by arbitrary choice the Z coordinate of the starting point of the Horn1TopLevelUpstr logical volume
  // We shift the drawing coordinate system
 // const double zD = z + 3.2852*inchDef*zResc; // Only valid for Horn1 !!!!! 
// By definition, drawing 8875.000-ME-363028 (Z=0 is ACTRNT)
 
 const double zR = z/inchDef; // back in inches. 
 const double argR = rCoeff + zR*zCoeff/zResc;
 if (argR < 0.) {
    std::ostringstream mStrStr; mStrStr << " Negative argument, z = " << z 
                                        << " argR " << argR << " zResc " << zResc << std::endl; 
    G4String mStr(mStrStr.str());
    G4Exception("LBNEHornRadialEquation::GetVal", " ", FatalErrorInArgument, mStr.c_str());
  } 
 const double radius = parabolic ? (std::sqrt(argR) + rOff) : argR;
// return radius*zResc*inchDef; v3r0p10
 return (radius*rResc*inchDef + rRescCst); // Dec 16: add a constant term for this linear transform 
}

void LBNEHornRadialEquation::test1() const {

  //Test at the end of first section of Horn1 (8875.112-MD 363104) 
  const double argZ = (3.2645 + 17.876 - 0.031);
  const double rTest = this->GetVal(argZ*25.4*CLHEP::mm*zResc);
  std::cerr << " LBNEHornRadialEquation::test1, argZ " << argZ << "  rTest (mm) " << rTest <<  std::endl;
  std::cerr << " inchDef " << inchDef << " zCoeff " << zCoeff << " rCoeff " << rCoeff << " rOff " << rOff << std::endl;
  const double delta = 2.0*rTest - 1.6326*25.4*CLHEP::mm*rResc; // Actually, the drawing says 1.6 !!! mistake of 700 microns. 
  std::cerr << " delta (mm) " << delta << " zResc " << zResc << " rResc " << rResc <<  std::endl;
  if (std::abs(delta) > 0.127*CLHEP::mm) { // 5 mill tolerance 
    G4Exception("LBNEHornRadialEquation::test1", " ", FatalErrorInArgument,
                " Horn1 Equation 0 inconsistent with drawing 363104"); 
  } 
}