LBNEDetectorConstruction.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------// 
// $Id: LBNEDetectorConstruction.cc,v 1.3.2.40 2013/12/24 20:31:46 lebrun Exp $
//---------------------------------------------------------------------------// 

#include <fstream>
#include <vector>
#include <string>

#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 "LBNFMagneticFieldPolyconeHorn.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"

#include "G4RegionStore.hh"
#include "G4SolidStore.hh"
#include "G4GeometryManager.hh"
#include "G4FieldManager.hh"
#include "LBNEVolumePlacements.hh"
#include "LBNEDetectorMessenger.hh"
#include "LBNERunManager.hh"
#include "G4GDMLParser.hh"
#include "G4UIcmdWithADouble.hh"

#include "G4RunManager.hh"

#include "G4VisExtent.hh"

//---------------------------------------------------------------------------// 
// Constructor, Destructor, and Initialization
//---------------------------------------------------------------------------// 

LBNEDetectorConstruction::LBNEDetectorConstruction()
{
  fPlacementHandler = LBNEVolumePlacements::Instance(); // Minimal setup for the Placement algorithm. 
  fDetectorMessenger = new LBNEDetectorMessenger(this);
//
  // Used only in placing the absorber.. 
  fBeamlineAngle = -101*CLHEP::mrad;
  fRotBeamlineAngle.rotateX(-fBeamlineAngle);  

  InitializeMaterials();
  Initialize();
  fHasBeenConstructed = false; 
//  Construct(); Not yet!  Need to read the data card first... 
  fHornCurrent = 200.*CLHEP::ampere*1000; // in kA, defined via Detector GUImessenger if need be
  fSkinDepthInnerRad = 1.0e10*CLHEP::mm; // infinitly long skin depth. (assume in default ) 
  fDeltaEccentricityIO.resize(5); // max number of horns 
  fDeltaEllipticityI.resize(5); // max number of horns 
  fCurrentEqualizerLongAbsLength.resize(5);
  fFileNameFieldMapForCE.resize(5);
  fCurrentEqualizerQuadAmpl.resize(5);
  fCurrentEqualizerOctAmpl.resize(5);
  fCurrentMultiplier.resize(5);
  for (size_t k=0; k!= fDeltaEccentricityIO.size(); k++) fDeltaEccentricityIO[k] = 0.;
  for (size_t k=0; k!= fDeltaEllipticityI.size(); k++) fDeltaEllipticityI[k] = 0.;
  for (size_t k=0; k!= fCurrentEqualizerLongAbsLength.size(); k++) fCurrentEqualizerLongAbsLength[k] = 0.;
  for (size_t k=0; k!= fCurrentEqualizerQuadAmpl.size(); k++) fCurrentEqualizerQuadAmpl[k] = 0.;
  for (size_t k=0; k!= fCurrentEqualizerOctAmpl.size(); k++) fCurrentEqualizerOctAmpl[k] = 0.;
  for (size_t k=0; k!= fCurrentMultiplier.size(); k++) fCurrentMultiplier[k] = 1.;
  for (size_t k=0; k!= fFileNameFieldMapForCE.size(); k++) fFileNameFieldMapForCE[k] = std::string("?");
// 
// for Perfect Focusing setting.. All Pions tracks withh be destroyed and recereated with px = py =pz 
// in Stepping Action, should the Post Step Z coord. be greater than.. 
//
   fZCoordForPerfectFocusing = 3. * 13.5 * 1.e9 * 24.* 3600. * 365.25 * 3.0e8 *CLHEP::m; // There...  
   fConstructSimpAbsorber = false;
   fConstructSculptedAbsorber = true;
   fDisableSpoiler = false;
   fDisableSculptedLayers = false;  
   fDisableMask = false;  
   fExpandAlLayers = false;
   fRemoveLayers = 0;
   //
   // August/September 2016.  Under construction..  
   //
   fConstructLBNFNano = false; // We do not want to depress Jim H., so, it is off.. 
   fInstallLBNFNanoUpstreamMicroStrip = false;
   fInstallLBNFNanoMiddleMicroStrip = false;
   fInstallLBNFNanoDownstreamMicroStrip = false;
   fInstallLBNFNanoMagnetUpstr = false;
   fInstallLBNFNanoMagnetDwnstr = false;
   fInstallLBNFNanoCerenkov = false;
   fInstallLBNFNanoCalorimeter = false;
      
}


LBNEDetectorConstruction::~LBNEDetectorConstruction()
{

  DestroyMaterials();

  //for(size_t i = 0; i< fSubVolumes.size(); i++){
  //  delete fSubVolumes[i];
  //}

  delete fDetectorMessenger;
//  delete fPlacementHandler; A static struct now, no point deleting it 
}

 // Obsolete.
void LBNEDetectorConstruction::InitializeSubVolumes()
{
//  fDecayPipe = new LBNEDecayPipe("DecayPipe");
// fHadronAbsorber = new LBNEHadronAbsorber("HadronAbsorber");
//  fStandardPerson = new LBNEStandardPerson("StandardPerson");
  /*
  fTarget = new LBNETarget();
  fBaffle = new LBNEBaffle();
  fHornAssembly = new LBNEHornAssembly();
  fHadronAbsorber = new LBNEHadronAbsorber();
  fTarget->SetDefaults();
  fBaffle->SetDefaults();
  fHornAssembly->SetDefaults();
  fDecayPipe->SetDefaults();
  fHadronAbsorber->SetDefaults();
  */
 // fSubVolumes.clear();
}

void LBNEDetectorConstruction::Initialize()
{
  // Set standard (and safe) values for class variables
  fSimulationType = "Standard Neutrino Beam"; 
  fConstructTarget = true;
  
}


//-------------------------------------------------------------------------

void LBNEDetectorConstruction::InitializeMaterials() {

  G4Element* elH  = new G4Element("Hydrogen","H" , 1., 1.01*CLHEP::g/CLHEP::mole);
  new G4Element("Helium","He" , 2., 4.003*CLHEP::g/CLHEP::mole);
  elC  = new G4Element("Carbon","C" , 6., 12.01*CLHEP::g/CLHEP::mole);
  elN  = new G4Element("Nitrogen","N" , 7., 14.01*CLHEP::g/CLHEP::mole);
  elO  = new G4Element("Oxygen"  ,"O" , 8., 16.00*CLHEP::g/CLHEP::mole); 
  G4Element* elNa = new G4Element("Sodium"  ,"Na" , 11., 22.99*CLHEP::g/CLHEP::mole); 
  G4Element* elMg = new G4Element("Magnesium"  ,"Mg" , 12., 24.305*CLHEP::g/CLHEP::mole); 
  G4Element* elAl = new G4Element("Aluminum"  ,"Al" , 13., 26.98*CLHEP::g/CLHEP::mole); 
  G4Element* elSi = new G4Element("Silicon"  ,"Si" , 14., 28.09*CLHEP::g/CLHEP::mole); 
  G4Element* elP  = new G4Element("Phosphorus"  ,"P" , 15., 30.974*CLHEP::g/CLHEP::mole); 
  G4Element* elS  = new G4Element("Sulfur"  ,"S" , 16., 32.065*CLHEP::g/CLHEP::mole); 
  G4Element* elK  = new G4Element("Potassium"  ,"K" , 19., 39.1*CLHEP::g/CLHEP::mole); 
  G4Element* elCa = new G4Element("Calcium"  ,"Ca" , 20., 40.09*CLHEP::g/CLHEP::mole); 
  G4Element* elTi = new G4Element("Titanium"  ,"Ti" , 22., 47.867*CLHEP::g/CLHEP::mole); 
  G4Element* elCr = new G4Element("Chromium"  ,"Cr" , 24., 51.9961*CLHEP::g/CLHEP::mole); 
  G4Element* elMn = new G4Element("Manganese"  ,"Mn" , 25., 54.938*CLHEP::g/CLHEP::mole); 
  elFe = new G4Element("Iron"  ,"Fe" , 26., 55.85*CLHEP::g/CLHEP::mole); 
  G4Element* elNi = new G4Element("Nickel"  ,"Ni" , 28., 58.6934*CLHEP::g/CLHEP::mole); 
  G4Element* elCu = new G4Element("Copper"  ,"Cu" , 29., 63.546*CLHEP::g/CLHEP::mole); 
  G4Element* elHg = new G4Element("Mercury"  ,"Hg" , 80., 200.59*CLHEP::g/CLHEP::mole); 
  G4Element* elMo = new G4Element("Molybdenum"  ,"Mo" , 42., 95.96*CLHEP::g/CLHEP::mole); 
  G4Element* elNiob = new G4Element("Niobium"  ,"Nb" , 41., 92.90637*CLHEP::g/CLHEP::mole); 
  G4Element* elCo = new G4Element("Cobalt"  ,"Nb" , 27., 58.9932*CLHEP::g/CLHEP::mole); 
  G4Element* elV = new G4Element("Vanadium"  ,"V" , 23, 50.9415*CLHEP::g/CLHEP::mole); 


  Air = new G4Material("Air"  , 1.290*CLHEP::mg/CLHEP::cm3, 2);
  Air->AddElement(elN, 0.7);
  Air->AddElement(elO, 0.3);
  
  CT852 = new G4Material("CT852", 7.75*CLHEP::g/CLHEP::cm3, 10); 
  CT852->AddElement(elC,  0.001); 
  CT852->AddElement(elSi, 0.008); 
  CT852->AddElement(elMn, 0.008); 
  CT852->AddElement(elCr, 0.13); 
  CT852->AddElement(elS,  0.00025); 
  CT852->AddElement(elP,  0.0003); 
  CT852->AddElement(elTi, 0.002); 
  CT852->AddElement(elCu, 0.003); 
  CT852->AddElement(elNi, 0.006); 
  CT852->AddElement(elFe, 0.84145); 

  Steel316 = new G4Material("Steel316", 8.0*CLHEP::g/CLHEP::cm3, 9); 
  // Reference: Google search, found Anderson Schumake company.  
  Steel316->AddElement(elC,  0.015); 
  Steel316->AddElement(elSi, 0.005); 
  Steel316->AddElement(elMn, 0.008); 
  Steel316->AddElement(elMo, 0.01); 
  Steel316->AddElement(elCr, 0.17); 
  Steel316->AddElement(elS,  0.00015); 
  Steel316->AddElement(elP,  0.0003); 
  Steel316->AddElement(elNi, 0.12); 
  Steel316->AddElement(elFe, 0.6716); 
  
  Titanium = new G4Material("Titanium", 22, 47.867*CLHEP::g/CLHEP::mole, 4.506*CLHEP::g/CLHEP::cm3);
  
  TitaniumG5 = new G4Material("TitaniumG5", 4.42*CLHEP::g/CLHEP::cm3, 3);
  TitaniumG5->AddElement(elAl, 0.06);
  TitaniumG5->AddElement(elV, 0.04);
  TitaniumG5->AddElement(elTi, 0.9);
  
 // ASTM836 steel
  Slab_Stl = new G4Material("Slab_Stl", 7.8416*CLHEP::g/CLHEP::cm3, 6);
  Slab_Stl->AddElement(elC,  0.001);
  Slab_Stl->AddElement(elSi, 0.001);
  Slab_Stl->AddElement(elMn, 0.004);
  Slab_Stl->AddElement(elFe, 0.982);
  Slab_Stl->AddElement(elNi, 0.01);
  Slab_Stl->AddElement(elCu, 0.002);

  // BluBlock Steel
  Blu_Stl = new G4Material("Blu_Stl", 7.25*CLHEP::g/CLHEP::cm3, 1);
  Blu_Stl->AddElement(elFe, 1.0);

  Water = new G4Material("Water", 1.0*CLHEP::g/CLHEP::cm3, 2);
  Water->AddElement(elH, 2);
  Water->AddElement(elO, 1);

  Vacuum = new G4Material("Vacuum", 2.376e-15*CLHEP::g/CLHEP::cm3,1,kStateGas,300.*CLHEP::kelvin,2.0e-7*CLHEP::bar);
  Vacuum->AddMaterial(Air, 1.);

  new G4Material("Helium", 2, 4.0026*CLHEP::g/CLHEP::mole, 2.55*0.1785*CLHEP::kg/CLHEP::m3, kStateGas,
                          300*CLHEP::kelvin, 2.55*CLHEP::atmosphere); // to fill the canister.For the decay pipe, see below. 
                          
  new G4Material("HeliumTarget", 2, 4.0026*CLHEP::g/CLHEP::mole, 1.7436*0.1785*CLHEP::kg/CLHEP::m3, kStateGas,
                          350*CLHEP::kelvin, 1.36*CLHEP::atmosphere); //20 psi.  The factor of 1.7436 assume perfect gas. 
                          // density proportional to temperature and pression.  

  new G4Material("Aluminum", 13, 26.98*CLHEP::g/CLHEP::mole, 2.7*CLHEP::g/CLHEP::cm3);
  new G4Material("Argon", 18, 39.948*CLHEP::g/CLHEP::mole, 1.784*CLHEP::kg/CLHEP::m3, kStateGas,
                          300*CLHEP::kelvin, CLHEP::atmosphere);
  new G4Material("Lead", 82, 207.19*CLHEP::g/CLHEP::mole, 11.35*CLHEP::g/CLHEP::cm3);
  new G4Material("Iron", 26, 55.85*CLHEP::g/CLHEP::mole, 7.86999*CLHEP::g/CLHEP::cm3);
  new G4Material("Uranium", 92, 238.028*CLHEP::g/CLHEP::mole, 19.1*CLHEP::g/CLHEP::cm3);

  Concrete = new G4Material("Concrete", 2.03*CLHEP::g/CLHEP::cm3, 10);
  Concrete->AddElement( elH,  0.01);
  Concrete->AddElement( elO,  0.529);
  Concrete->AddElement( elNa, 0.016);
  Concrete->AddElement( elHg, 0.002);
  Concrete->AddElement( elAl, 0.034);
  Concrete->AddElement( elSi, 0.337);
  Concrete->AddElement( elK,  0.013);
  Concrete->AddElement( elCa, 0.044);
  Concrete->AddElement( elFe, 0.014);
  Concrete->AddElement( elC,  0.001);
  
  G4Material *rockMat = new G4Material( "rockMat", 2.78*CLHEP::g/CLHEP::cm3, 4 ); //CaMg(CO3)2
  rockMat->AddElement( elCa, 1);
  rockMat->AddElement( elMg, 1); 
  rockMat->AddElement( elC,  2); 
  rockMat->AddElement( elO,  6); 

  G4Material *aluminaMat = new G4Material( "Alumina", 3.92*CLHEP::g/CLHEP::cm3, 2 ); //Al2O3 See 
  // http://www.ortechceramics.com/creamic-materials/alumina-ceramics/?gclid=CKnTttqOns0CFYM2aQod540Nhg
  aluminaMat->AddElement( elAl, 2);
  aluminaMat->AddElement( elO,  3);
   
  G4Material *inconelMat = new G4Material( "Inconel718", 8.1933*CLHEP::g/CLHEP::cm3, 8 ); 
  // See http://wwwspecialmetals.com/  We take the middle range.. 
  inconelMat->AddElement( elNi, 0.525);
  inconelMat->AddElement( elCr,  0.19);
  inconelMat->AddElement( elFe,  0.17925);
  inconelMat->AddElement( elNiob,  0.05125);
  inconelMat->AddElement( elMo,  0.0305);
  inconelMat->AddElement( elTi,  0.009);
  inconelMat->AddElement( elAl,  0.005);
  inconelMat->AddElement( elCo,  0.01);
 
   
  graphiteBaffle = new G4Material( "GraphiteBaffle", 1.78*CLHEP::g/CLHEP::cm3, 3 ); //Carbon, air (Nitrogen and oxygen)
  graphiteBaffle->AddElement( elC,  0.99); // 
  graphiteBaffle->AddElement( elN,  0.007); //  
  graphiteBaffle->AddElement( elO,  0.003); // 
 
  G4Material* graphite = new G4Material( "Graphite", 1.78*CLHEP::g/CLHEP::cm3, 3 ); //Carbon, air (Nitrogen and oxygen)
  graphite->AddElement( elC,  0.99); // 
  graphite->AddElement( elN,  0.007); //  
  graphite->AddElement( elO,  0.003); // 
  
  new G4Material("Beryllium", 4, 9.0122*CLHEP::g/CLHEP::mole, 1.85*CLHEP::g/CLHEP::cm3); 
  
  G4Material *Mylar = new G4Material("Mylar", 1.397*CLHEP::g/CLHEP::cm3, 3);
  Mylar->AddElement(elC, 10);
  Mylar->AddElement(elH, 8);
  Mylar->AddElement(elO, 4);

}
//
// Declare the material for the target after the data cards have been read. 
// 
void LBNEDetectorConstruction::InitializeMaterialsPostPreIdle() {
  // In case they are different.. 
  G4String aNameT(fPlacementHandler->GetTargetMaterialName());
  if ((aNameT == G4String("Graphite")) || (aNameT == G4String("graphite")) ||
      (aNameT == G4String("Carbon")) || (aNameT == G4String("carbon"))) { 
  
    Target = new G4Material( "Target", 
             fPlacementHandler->GetTargetDensity(), 3 ); //Carbon, air (Nitrogen and oxigen),
                                                                  // Assume density of POCO ZXF-5Q  
    Target->AddElement( elC,  0.99); 
    Target->AddElement( elN,  0.007); 
    Target->AddElement( elO,  0.003);
  } else if ((aNameT == G4String("Beryllium")) 
               || (aNameT == G4String("beryllium"))) {
     new G4Material("Target", 4, 9.0122*CLHEP::g/CLHEP::mole, 1.85*CLHEP::g/CLHEP::cm3); 
     // Set new specific nuclear interaction length
     fPlacementHandler->SetTargetSpecificLambda(77.8*CLHEP::g/CLHEP::cm2);
     // Let the placement handler know about the new target density
     fPlacementHandler->SetTargetDensity(1.85*CLHEP::g/CLHEP::cm3);

  } else if ((aNameT == G4String("Aluminium")) 
            || (aNameT == G4String("aluminium"))) {   
     new G4Material("Target", 13, 26.98*CLHEP::g/CLHEP::mole, 2.7*CLHEP::g/CLHEP::cm3);
     // Set new specific nuclear interaction length
     fPlacementHandler->SetTargetSpecificLambda(107.2*CLHEP::g/CLHEP::cm2);
     // Let the placement handler know about the new target density
     fPlacementHandler->SetTargetDensity(2.7*CLHEP::g/CLHEP::cm3);

  } else if (aNameT == G4String("lightCarbon")) { 
      Target = new G4Material( "Target", 
             fPlacementHandler->GetTargetDensity()/100., 3 );
      Target->AddElement( elC,  0.99); 
      Target->AddElement( elN,  0.007); 
      Target->AddElement( elO,  0.003);
  } else {
    G4String mess(" Non-standard material for the target: "); 
    mess += aNameT + G4String (" .  \n");
    mess += G4String("... Please upgrade the code after consultation with mechanical engineers\n.");
    G4Exception("LBNEDetectorConstruction::InitializeMaterialsPostPreIdle", 
                " ",  FatalErrorInArgument, mess.c_str());
  }
  
  G4String aNameG(fPlacementHandler->GetDecayPipeGas());
  if ((aNameG == G4String("Air")) || (aNameG == G4String("air"))) {
    G4Material *gas = new G4Material("DecayPipeGas" , 1.290*CLHEP::mg/CLHEP::cm3, 2);
    gas->AddElement(elN, 0.7);
    gas->AddElement(elO, 0.3);
  } else if ((aNameG == G4String("Helium")) || (aNameG == G4String("helium"))) {
      new G4Material("DecayPipeGas", 2, 4.0026*CLHEP::g/CLHEP::mole, 0.1785*CLHEP::kg/CLHEP::m3, kStateGas,
                          300*CLHEP::kelvin, 1.0*CLHEP::atmosphere);
  }  else {
    G4String mess(" Non-standard gas  : "); 
    mess += aNameG + G4String (" in the decay pipe .  \n");
    mess += G4String("... Please upgrade the code after consultation with mechanical engineers\n. ");
    G4Exception("LBNEDetectorConstruction::InitializeMaterialsPostPreIdle", 
                " ",  FatalErrorInArgument, mess.c_str());
  }
  
  if (fConstructLBNFNano) {
    new G4Material("Silicon", 14, 28.0855*CLHEP::g/CLHEP::mole, 2.329*CLHEP::g/CLHEP::cm3); 
    G4Element* elSr = new G4Element("Strontium"  ,"Sr" , 38., 87.62*CLHEP::g/CLHEP::mole); 
    G4Material* myFerrite = new G4Material("FerriteC8", 5.*CLHEP::g/CLHEP::cm3, 3);
    myFerrite->AddElement(elSr, 1);
    myFerrite->AddElement(elO, 19);
    myFerrite->AddElement(elFe, 12);
  }
  
  // For the 1st and possibly 2nd target modules. These are separated from the previous target
  // material definitions so that the other target definitions are left unchanged. This repeats
  // earlier code but allows us to define materials that won't interfere with the other target setups
  
  if (fPlacementHandler->GetUseTargetModule()) {

      G4String matName = fPlacementHandler->GetTargetMaterialName();
      matName.toLower();

      // Define the Target1 material: elements, density, interaction lengths
      // Allow graphite, beryllium, aluminium (low-Z)
      if (matName == G4String("graphite")) {

          // Carbon, air (Nitrogen and oxygen), assume density of POCO ZXF-5Q 
          G4Material* target1 = new G4Material("Target1", 1.78*CLHEP::g/CLHEP::cm3, 3);
          target1->AddElement(elC, 0.99);
          target1->AddElement(elN, 0.007); 
          target1->AddElement(elO, 0.003);
          fPlacementHandler->SetTargetSpecificLambda(85.8*CLHEP::g/CLHEP::cm2);
          fPlacementHandler->SetTargetDensity(1.78*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("beryllium")) {

          new G4Material("Target1", 4, 9.0122*CLHEP::g/CLHEP::mole, 1.85*CLHEP::g/CLHEP::cm3); 
          // Set new specific nuclear interaction length
          fPlacementHandler->SetTargetSpecificLambda(77.8*CLHEP::g/CLHEP::cm2);
          // Let the placement handler know about the new target density
          fPlacementHandler->SetTargetDensity(1.85*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("aluminium")) {

          new G4Material("Target1", 13, 26.98*CLHEP::g/CLHEP::mole, 2.7*CLHEP::g/CLHEP::cm3);
          // Set new specific nuclear interaction length
          fPlacementHandler->SetTargetSpecificLambda(107.2*CLHEP::g/CLHEP::cm2);
          // Let the placement handler know about the new target density
          fPlacementHandler->SetTargetDensity(2.7*CLHEP::g/CLHEP::cm3);

      }

  }

  if (fPlacementHandler->GetUseTarget2Module()) {

      // Define the Target2 material. Allow either low or high Z:
      // graphite, beryllium, aluminium, tungsten, tantalum, ...

      G4String matName = fPlacementHandler->GetTarget2MaterialName();
      matName.toLower();
      if (matName == G4String("graphite")) {

          // Carbon, air (Nitrogen and oxygen), assume density of POCO ZXF-5Q 
          G4Material* target2 = new G4Material("Target2", 1.78*CLHEP::g/CLHEP::cm3, 3);
          target2->AddElement(elC, 0.99);
          target2->AddElement(elN, 0.007); 
          target2->AddElement(elO, 0.003);
          fPlacementHandler->SetTarget2SpecificLambda(85.8*CLHEP::g/CLHEP::cm2);
          fPlacementHandler->SetTarget2Density(1.78*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("beryllium")) {

          new G4Material("Target2", 4, 9.0122*CLHEP::g/CLHEP::mole, 1.85*CLHEP::g/CLHEP::cm3); 
          // Set new specific nuclear interaction length
          fPlacementHandler->SetTarget2SpecificLambda(77.8*CLHEP::g/CLHEP::cm2);
          // Let the placement handler know about the new target density
          fPlacementHandler->SetTarget2Density(1.85*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("aluminium")) {

          new G4Material("Target2", 13, 26.98*CLHEP::g/CLHEP::mole, 2.7*CLHEP::g/CLHEP::cm3);
          // Set new specific nuclear interaction length
          fPlacementHandler->SetTarget2SpecificLambda(107.2*CLHEP::g/CLHEP::cm2);
          // Let the placement handler know about the new target density
          fPlacementHandler->SetTarget2Density(2.7*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("tantalum")) {

          new G4Material("Target2", 73, 180.94788*CLHEP::g/CLHEP::mole, 16.65*CLHEP::g/CLHEP::cm3);
          // Nuclear interaction length
          fPlacementHandler->SetTarget2SpecificLambda(191.0*CLHEP::g/CLHEP::cm2);
          // Target density
          fPlacementHandler->SetTarget2Density(16.65*CLHEP::g/CLHEP::cm3);

      } else if (matName == G4String("tungsten")) {

          new G4Material("Target2", 74, 183.84*CLHEP::g/CLHEP::mole, 19.3*CLHEP::g/CLHEP::cm3);
          // Nuclear interaction length
          fPlacementHandler->SetTarget2SpecificLambda(191.9*CLHEP::g/CLHEP::cm2);
          // Target density
          fPlacementHandler->SetTarget2Density(19.3*CLHEP::g/CLHEP::cm3);

      } 



  }

}
G4VisAttributes* LBNEDetectorConstruction::GetMaterialVisAttrib(G4String mName){
  G4VisAttributes* visAtt;
  if(mName == "Vacuum")  visAtt = new G4VisAttributes(false);
  else if(mName=="Aluminum") visAtt = new G4VisAttributes(G4Color(0.2, 0.8, 1));
  else if(mName=="Air") visAtt = new G4VisAttributes(G4Color(0.6,0.7,0.8));
  else if(mName=="Iron" || mName=="Slab_Stl") visAtt=new G4VisAttributes(G4Color(0.5,0.3,0));
  else if(mName=="Concrete") visAtt = new G4VisAttributes(G4Color(0.75,0.85,0.95));
  else visAtt = new G4VisAttributes(G4Color(1,0,0));
  return visAtt;
}

void LBNEDetectorConstruction::DestroyMaterials()
{
  // Destroy all allocated elements and materials
  size_t i;
  G4MaterialTable* matTable = (G4MaterialTable*)G4Material::GetMaterialTable();
  for(i=0;i<matTable->size();i++)
  { delete (*(matTable))[i]; }
  matTable->clear();
  G4ElementTable* elemTable = (G4ElementTable*)G4Element::GetElementTable();
  for(i=0;i<elemTable->size();i++)
  { delete (*(elemTable))[i]; }
  elemTable->clear();
}

G4VPhysicalVolume* LBNEDetectorConstruction::Construct() {

  if (fHasBeenConstructed) {
     std::cerr << " WARNING: LBNEDetectorConstruction::Construct, already done, skip " << std::endl;
     return fRock;
  }
  
  this->InitializeMaterialsPostPreIdle(); 
  
  std::cout << " LBNEDetectorConstruction::Construct, Start !!! " << std::endl;
  std::cerr << " LBNEDetectorConstruction::Construct, Start !!! " << std::endl;
  
  fRockX = 60.0*CLHEP::m;
  fRockY = 60.0*CLHEP::m;
  fRockLength = fPlacementHandler->GetTotalLengthOfRock() + 4.0*CLHEP::cm;
//  std::cerr << "  ............  fRockLength " << fRockLength << std::endl;
   
      // See LBNEVolumePlacements constructor. 
  G4Box* ROCK_solid = new G4Box("ROCK_solid",fRockX/2, fRockY/2, fRockLength/2);
  G4LogicalVolume *RockLogical = 
            new G4LogicalVolume(ROCK_solid,
                                G4Material::GetMaterial("Concrete"),
                                "RockLogical",0,0,0); 
  fPlacementHandler->Initialize(RockLogical); // sort of a noop for now..                               
  //RockLogical->SetVisAttributes(G4VisAttributes::Invisible);
  fRock = new G4PVPlacement(0,G4ThreeVector(),RockLogical,"ROCK",0,false,0);

  
  // First create the Target Hall, Pipe Hall, and Absorber Hall, and then
  // connect them together.
  fPlacementHandler->SetTotalLengthOfRock(fRockLength);
  LBNEVolumePlacementData *plTunnel =fPlacementHandler->Create(G4String("Tunnel")); 
    std::cerr << " Placement data for volume Tunnel, half length  " 
              << plTunnel->fParams[2]/2. << " At z  " << plTunnel->fPosition[2] << std::endl;
  G4VPhysicalVolume* tunnel = fPlacementHandler->PlaceFinal(G4String("Tunnel"), fRock ); // like Rock, oversized. Air. 

  if (fPlacementHandler->GetUseHorn1Polycone())
           fPlacementHandler->UpdateParamsForHorn1MotherPoly(); // This could change the length of Horn1.. 

    // 
    // Adjust first the mother volumes. 
    //
  const size_t usedNumberOfHornsPoly = static_cast<size_t>(fPlacementHandler->GetNumberOfHornsPolycone());
  const bool useHornsPolycone = usedNumberOfHornsPoly > 0;
  bool implementHorn1AsLBNEHorn = !useHornsPolycone;// New, April 12 2016.
    // The field map are expected to be the same for Horn1 except that the LBNF field map 
    // are more realistic.. Ellicpticity, eccentricity and imprefect current equalizer are in, 
    // optionally 
  if (useHornsPolycone) {  
         for (size_t iH=0; iH != usedNumberOfHornsPoly; iH++) {
            fPlacementHandler->UpdateParamsForHornMotherPolyNum(iH);
         }
  }
  if (!fPlacementHandler->GetUseMarsTargetHorns()) {
  
     fPlacementHandler->ExtendChaseLengthForHorn2(); // Feb 2015: Optimization.. 
     
     LBNEVolumePlacementData *plTHAH = fPlacementHandler->Create(G4String("TargetHallAndHorn1"));
    std::cerr << " Placement data for volume TargetHallAndHorn1, half length  " 
              << plTHAH->fParams[2]/2. << " At z  " << plTHAH->fPosition[2] << std::endl;
    // 
    // Before placing the container volume for the target region + horn1, define these two volumes, 
    // as these two are adjacent. The boundary is "coordinate zero.", respecting older convention. 
    //
    LBNEVolumePlacementData *plDatUTA = fPlacementHandler->Create(G4String("UpstreamTargetAssembly"));
    std::cerr << " Placement data for volume UpstreamTargetAssembly, half length  " 
              << plDatUTA->fParams[2]/2. << " At z  " << plDatUTA->fPosition[2] << std::endl;
    
    // Is this needed ?  I guess so.. but plH1Dat is not used anywhere.. for reference only..  
    fPlacementHandler->Create(G4String("Horn1Hall"));
    
    G4VPhysicalVolume* targethorn1Phys = fPlacementHandler->PlaceFinal(G4String("TargetHallAndHorn1"), tunnel);
    std::cerr << " Detector Construction, TargetHallAndHorn1 is placed " << std::endl;
    fPlacementHandler->SetTargetHorn1HallPhysPtr(targethorn1Phys);
   //
    G4PVPlacement * upstreamTargetAssPhys = 
        fPlacementHandler->PlaceFinal(G4String("UpstreamTargetAssembly"), targethorn1Phys); 
    std::cerr << " Detector Construction, UpstreamTargetAssembly is placed " << std::endl;
        
    //
    // Horn code. Include multiple Polycone horns..
    //
    
    std::vector<G4PVPlacement*> vHornsAfterHornX; // where X is 1 or zero, depending if we define Horn1 
    // field as for LBNF CDR 2015. 
    LBNEVolumePlacementData *plH1PolyDat = fPlacementHandler->Create(G4String("Horn1PolyM1"));
    G4PVPlacement *vHorn1 = fPlacementHandler->PlaceFinal(G4String("Horn1PolyM1"), targethorn1Phys);
    std::cerr << " Horn1PolyM1 Half Length " << plH1PolyDat->fParams[2]/2. << " At position " << plH1PolyDat->fPosition[2] << std::endl;
    fPlacementHandler->SetHorn1HallPhysPtr(vHorn1);
    
    fPlacementHandler->PlaceFinalUpstrTarget((G4PVPlacement*) upstreamTargetAssPhys); 
            // Method is a misnomer, if after Aug 2014, as the target is no longer split up  
    
    if (fPlacementHandler->GetUseLBNFOptimConceptDesignHornA()) {   
      fPlacementHandler->PlaceFinalLBNFConceptHornA();
      fPlacementHandler->PlaceFinalLBNFConceptTgtSupport();
    } else fPlacementHandler->PlaceFinalNoSplitHorn1((G4PVPlacement*) vHorn1, (G4PVPlacement*) targethorn1Phys);
    if (implementHorn1AsLBNEHorn)
       vHornsAfterHornX.push_back(0); // ==> not used for Horn1, used the old (LBNE < Summer 2015) magentic field. 
    else vHornsAfterHornX.push_back(vHorn1);
    LBNEVolumePlacementData *plH2Dat = 0;
    G4PVPlacement *vHorn2  = 0;
    if (fPlacementHandler->GetUseCDR2015Optimized() || (!useHornsPolycone) 
        || (fPlacementHandler->GetUseHorn1Polycone())) { 
         std::cerr << " Old Horn2 installation... " << std::endl;
        plH2Dat = fPlacementHandler->Create(G4String("Horn2Hall"));
    
       vHorn2 = fPlacementHandler->PlaceFinal(G4String("Horn2Hall"), tunnel); 
    
       fPlacementHandler->PlaceFinalHorn2(vHorn2);
    } else {
     // 
//       plH2Dat = fPlacementHandler->Create(G4String("LBNFChaseDwnstrHorn1"));
    
//       G4PVPlacement *vHorns = fPlacementHandler->PlaceFinal(G4String("LBNFChaseDwnstrHorn1"), tunnel);
         std::cerr << " Polycone Horn2/Horn3 (or B & C)  installation...usedNumberOfHornsPoly " << usedNumberOfHornsPoly << std::endl;
       if (usedNumberOfHornsPoly > 1) {  
         for (size_t iH=1; iH != usedNumberOfHornsPoly; iH++) {
           std::cerr << " ... Placing Horn " << (iH+1) << " in the tunnel.. " << std::endl;
           std::ostringstream nameHStr; nameHStr << "LBNFSimpleHorn" << (iH+1) << "Container";
           std::string nameH(nameHStr.str());
           fPlacementHandler->Create(nameH);   
           G4PVPlacement *vHorn = fPlacementHandler->PlaceFinal(nameH, tunnel); 
           vHornsAfterHornX.push_back(vHorn); // to define the magnetic field. 
           std::cerr << " ... Mother volume placed... " << std::endl;
           switch (iH) {
             case 1:
               if (fPlacementHandler->GetUseLBNFOptimConceptDesignHornB()) { 
                 fPlacementHandler->PlaceFinalLBNFConceptHornB();
                 fPlacementHandler->PlaceFinalLBNFConceptStripLinesConnectHornB();
               } else fPlacementHandler->PlaceFinalSimpleHornPolyNumber(iH, vHorn);
               break;
             case 2:
               if (fPlacementHandler->GetUseLBNFOptimConceptDesignHornC()) { 
                 fPlacementHandler->PlaceFinalLBNFConceptHornC();
                 fPlacementHandler->PlaceFinalLBNFConceptStripLinesConnectHornC();
               } else fPlacementHandler->PlaceFinalSimpleHornPolyNumber(iH, vHorn);
               break;
           }
         }   
       }
       if (fPlacementHandler->GetConstructPlug()) {
          G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.;
          const double aRadIn =  fPlacementHandler->GetPlugInnerRadius();
          const double aRadOut =  fPlacementHandler->GetPlugOuterRadius();
          const double aLength =  fPlacementHandler->GetPlugLength();
          const std::string aPlugMaterial = fPlacementHandler->GetPlugMaterial();
          G4Tubs* plugTube = new G4Tubs("PlugOr2ndTgt", aRadIn, aRadOut, aLength*.5, 0., 360.*CLHEP::deg);
          G4LogicalVolume *plugl = new G4LogicalVolume(plugTube, G4Material::GetMaterial(aPlugMaterial.c_str()), "PlugOr2ndTgt");
          posTmp[2] = fPlacementHandler->GetPlugZPosition();
          G4String motherVolForPlug("TargetHallAndHorn1");
          const LBNEVolumePlacementData *plH1Tmp = 
              fPlacementHandler->Find(G4String("ConstructPlug"), motherVolForPlug, 
                                 G4String("LBNEDetectorConstruction::Construct"));
          const double zMinForPlug =  plH1Tmp->fPosition[2] + plH1Tmp->fParams[2]/2.;
          if ((fPlacementHandler->GetPlugZPosition() - aLength/2.- 0.5*CLHEP::mm) <  zMinForPlug) {  
             posTmp[2] -= plH1Tmp->fPosition[2];
             new G4PVPlacement((G4RotationMatrix *) 0, posTmp, plugl, "PlugOr2ndTgt_P",
                      targethorn1Phys->GetLogicalVolume(), false, 1, true);
             std::cerr << " ..... Plug (e.g., 2nd target) installed in Horn1 starting at Z = " << posTmp[2]-aLength/2. << std::endl;
          } else { 
             new G4PVPlacement((G4RotationMatrix *) 0, posTmp, plugl, "PlugOr2ndTgt_P",
                     tunnel->GetLogicalVolume(), false, 1, true);
             std::cerr << " ..... Plug (e.g., 2nd target) installed in Tunnel starting at Z = " << posTmp[2]-aLength/2. << std::endl;
          }
       }
    }
    // We now turn on the magnetic fields 
    if (implementHorn1AsLBNEHorn) { 
       LBNEMagneticFieldHorn *fieldHorn1 = new LBNEMagneticFieldHorn(true);
       fieldHorn1->SetHornCurrent(fHornCurrent);
       fieldHorn1->SetSkinDepthInnerRad(fSkinDepthInnerRad);
    // Not yet implemented via the old (non-Polycone ) geometry.  
//    if (std::abs(fDeltaEccentricityIO[0]) > 1.0e-6) fieldHorn1->SetDeltaEccentricityIO(fDeltaEccentricityIO[0]);
//    if (std::abs(fDeltaEllipticityI[0]) > 1.0e-6) fieldHorn1->SetDeltaEllipticityI(fDeltaEllipticityI[0]);
       G4FieldManager* aFieldMgr = new G4FieldManager(fieldHorn1); //create a local field                
      aFieldMgr->SetDetectorField(fieldHorn1); //set the field 
      aFieldMgr->CreateChordFinder(fieldHorn1); //create the objects which calculate the trajectory
      const LBNEVolumePlacementData *plH1Dat =
                 fPlacementHandler->Find("FieldHorn1", "Horn1PolyM1", "DetectorConstruction");
       plH1Dat->fCurrent->SetFieldManager(aFieldMgr,true); //attach the local field to logical volume
       fieldHorn1->dumpField();
    }
    this->ConstructLBNEHorn1TrackingPlane(tunnel); //CAll the Horn1TrackingPlane Amit Bashyal
    // Horn2 
    if ((!useHornsPolycone) || (fPlacementHandler->GetUseHorn1Polycone()) 
                            || (fPlacementHandler->GetUseCDR2015Optimized())) { 
        // UseHorn1Polycone was a flag that preceeded the use useHornsPolycones (plural) 
        // and the re-vamped UseCDR2015Optimized.. Convoluted..                     
      LBNEMagneticFieldHorn *fieldHorn2 = new LBNEMagneticFieldHorn(false);
      fieldHorn2->SetHornCurrent(fHornCurrent);
      fieldHorn2->SetSkinDepthInnerRad(fSkinDepthInnerRad);
      G4FieldManager* aFieldMgr2 = new G4FieldManager(fieldHorn2); //create a local field                
      aFieldMgr2->SetDetectorField(fieldHorn2); //set the field 
      aFieldMgr2->CreateChordFinder(fieldHorn2); //create the objects which calculate the trajectory
      plH2Dat->fCurrent->SetFieldManager(aFieldMgr2,true); //attach the local field to logical volume
    
      this->ConstructLBNEHorn2TrackingPlane(tunnel);  //Call the Horn2TrackingPlane Amit Bashyal
    } 
    //
    // Convoluted logic here...See above... 
    // 
    // variable numbers of horn1.. Typically, from 1 to 3. Horn2 already installed if UseCDR2015Optimized
    //
    if ((fPlacementHandler->GetUseCDR2015Optimized()) || 
        (useHornsPolycone) || (fPlacementHandler->GetUseHorn1Polycone())) {
      std::cerr << " Defining Polycone Mag field, usedNumberOfHornsPoly " << usedNumberOfHornsPoly << std::endl;
      for (size_t iH=0; iH != usedNumberOfHornsPoly; iH++) {
        if (implementHorn1AsLBNEHorn) continue;
        LBNFMagneticFieldPolyconeHorn *fieldHorn = new LBNFMagneticFieldPolyconeHorn(iH);
        std::ostringstream vStrStr; vStrStr << "LBNFSimpleHorn" << iH+1 << "Container";
        std::string vStr(vStrStr.str());
        if (iH == 0) vStr = std::string("Horn1PolyM1");
        const LBNEVolumePlacementData *plHxPolyDat = 
           fPlacementHandler->Find(G4String("FieldSetting") , vStr.c_str(), "FieldSetting" );
        const double zStart = plHxPolyDat->fPosition[2] - plHxPolyDat->fParams[2]/2.; 
        fieldHorn->SetZShiftDrawingCoordinate(zStart); // only used to produce the field map.
        std::cerr << " Horn1 will start at Z =" << zStart << std::endl;
        fieldHorn->SetEffectiveLength(plHxPolyDat->fParams[2]); // only used to produce the field map.
        fieldHorn->SetHornCurrent(fCurrentMultiplier[iH]*fHornCurrent);
        fieldHorn->SetSkinDepthInnerRad(fSkinDepthInnerRad);
        if (std::abs(fDeltaEccentricityIO[iH]) > 1.0e-20) fieldHorn->SetDeltaEccentricityIO(fDeltaEccentricityIO[iH]);
        if (std::abs(fDeltaEllipticityI[iH]) > 1.0e-20) fieldHorn->SetDeltaEllipticityI(fDeltaEllipticityI[iH]);
        if (std::abs(fCurrentEqualizerQuadAmpl[iH]) > 1.0e-20)
          fieldHorn->SetCurrentEqualizerQuadAmpl(fCurrentEqualizerQuadAmpl[iH]);
        if (std::abs(fCurrentEqualizerOctAmpl[iH]) > 1.0e-20)
          fieldHorn->SetCurrentEqualizerOctAmpl(fCurrentEqualizerOctAmpl[iH]);    
        fieldHorn->SetCurrentEqualizerLongAbsLength(fCurrentEqualizerLongAbsLength[iH]);
        if (fFileNameFieldMapForCE[iH] != std::string("?")) 
          fieldHorn->SetFileNameFieldMapForCE(fFileNameFieldMapForCE[iH]);
        G4FieldManager* aFieldMgr = new G4FieldManager(fieldHorn); //create a local field                
        aFieldMgr->SetDetectorField(fieldHorn); //set the field 
        aFieldMgr->CreateChordFinder(fieldHorn); //create the objects which calculate the trajectory
        //
        // If the Inner conductors are a bit elliptical, then one expect a small dipole field in the so-called field-free region,
        // inside the the voume comprised by the inner-inner conductor shape. 
        //
        if ((iH == 0) && ((std::abs(fDeltaEllipticityI[iH]) > 1.0e-6) || std::abs(fDeltaEccentricityIO[iH] > 1.0e-6)) ) {
           targethorn1Phys->GetLogicalVolume()->SetFieldManager(aFieldMgr,true);
           std::cerr << " Magnetic field attached to volume " << 
            targethorn1Phys->GetLogicalVolume()->GetName() << " Expecting field inside inner inner conductor " << std::endl;
        } else {
          vHornsAfterHornX[iH]->GetLogicalVolume()->SetFieldManager(aFieldMgr,true); //attach the local field to logical volume
          std::cerr << " Magnetic field attached to volume " << 
            vHornsAfterHornX[iH]->GetLogicalVolume()->GetName() << std::endl;
        }
         //
         // Some testing.. Now the horn current is turned on and eventual correction term are uploaded. 
         //
//      if (iH == 1) {
//         fieldHorn->TestDivergence(0);  
//         fieldHorn->TestDivergence(1);  
//         fieldHorn->TestDivergence(2); 
//      } 
//       if (iH == 0) {
//        std::cerr << " Generating field map Horn A " << std::endl;
//        fieldHorn->dumpField();
//       }
      }
    }
   // we forgot the baffle. So we install it now, that's O.K., but only conditonally. 
   // Agust 1 2016 : Laura asked us to remove it, as it confused the study of the neutrino flux with 
   // Conceptual Design HornA. What did shift the baffle upstream by 17.8 cm., compare to the simple optimized horn. 
   //
    if ((fPlacementHandler->GetInstallBaffle()) &&
        (!fPlacementHandler->GetUseSimpleTargetCylinder()) && (!fPlacementHandler->GetUseSimpleTargetBox())) {  
      fPlacementHandler->Create(G4String("Baffle"));
    // This will be a surveyed elements, but let us skip this step for now.    
      fPlacementHandler->PlaceFinal(G4String("Baffle"), upstreamTargetAssPhys);
      fPlacementHandler->Create(G4String("BaffleWindowUpstr"));
    // This will be a surveyed elements, but let us skip this step for now.    
      fPlacementHandler->PlaceFinal(G4String("BaffleWindowUpstr"), upstreamTargetAssPhys);
      fPlacementHandler->Create(G4String("BaffleWindowDownstr"));
    // This will be a surveyed elements, but let us skip this step for now.    
      fPlacementHandler->PlaceFinal(G4String("BaffleWindowDownstr"), upstreamTargetAssPhys);
    }
  //
  // Place the decay Pipe Snout, which contains the window, in case we have Helium gas. 
  //
  if (!fPlacementHandler->GetRemoveDecayPipeSnout()) {
    fPlacementHandler->Create(G4String("DecayPipeSnout")); // Now in Snout region 
    fPlacementHandler->PlaceFinalDecayPipeSnout((G4PVPlacement*) tunnel);
  }
  //
  // Place the LBNFNano Spectrometer 
  // Not yet uploaded in the repository. 
  //
  // if (fConstructLBNFNano) this->DoConstructLBNFNano(); // most of volume parameters are hardtyped... 
  //   
  // Place the decay pipe 
  //   
   fPlacementHandler->Create(G4String("DecayPipeHall"));
   G4PVPlacement *vDecayPipe = fPlacementHandler->PlaceFinal(G4String("DecayPipeHall"), tunnel);
   fPlacementHandler->Create(G4String("DecayPipeConcrete"));
   fPlacementHandler->Create(G4String("DecayPipeOuterWall"));
   fPlacementHandler->Create(G4String("DecayPipeWall"));
   fPlacementHandler->Create(G4String("DecayPipeVolume"));
   //   fPlacementHandler->Create(G4String("DecayPipeUpstrWindow")); // Now in Snout region 
   this->ConstructLBNEDecayPipeTrackingPlane(tunnel); //Amit Bashyal
   
   fPlacementHandler->PlaceFinal(G4String("DecayPipeConcrete"), vDecayPipe);
   fPlacementHandler->PlaceFinal(G4String("DecayPipeOuterWall"), vDecayPipe);
   fPlacementHandler->PlaceFinal(G4String("DecayPipeWall"), vDecayPipe);
   fPlacementHandler->PlaceFinal(G4String("DecayPipeVolume"), vDecayPipe);
    //
    // March 2015: Milind Diwan suggested to install a wire down the beam pipe. 
    //
    if (std::abs(fPlacementHandler->GetCurrentMilindWire()) > 1.0e-6) {
    
      LBNEMagneticFieldDecayPipe *fieldDecayPipe = new LBNEMagneticFieldDecayPipe(false); // we use the horn field
      fieldDecayPipe->SetWireCurrent(fPlacementHandler->GetCurrentMilindWire());
      fieldDecayPipe->SetWireRadius(fPlacementHandler->GetRadiusMilindWire());
      G4FieldManager* aFieldMgr3 = new G4FieldManager(fieldDecayPipe); //create a local field            
      aFieldMgr3->SetDetectorField(fieldDecayPipe); //set the field 
      aFieldMgr3->CreateChordFinder(fieldDecayPipe); //create the objects which calculate the trajectory
      const LBNEVolumePlacementData *plDCV =
                 fPlacementHandler->Find("FieldDecayPipe", "DecayPipeVolume", "DetectorConstruction");
       plDCV->fCurrent->SetFieldManager(aFieldMgr3,true); //attach the local field to logical volume
    }
    
   //Define geometry absorber from /LBNE/det/ConstructSimpAbsorber true/false option
   if (!fConstructLBNFNano) { 
     if (this->fConstructSimpAbsorber){
       this->ConstructLBNEHadronAbsorberSimple(tunnel);
     }
     else if (this->fConstructSculptedAbsorber){
       this->ConstructLBNEHadronAbsorberSculpted(tunnel);
     } else {
       this->ConstructLBNEHadronAbsorber(tunnel);
   
     }
   }
   if (fPlacementHandler->GetDoInstallShield()) { 
     this->ConstructLBNEShieldingHorn1(targethorn1Phys);
     if ((!useHornsPolycone) || (fPlacementHandler->GetUseHorn1Polycone())) {
       this->ConstructLBNEShieldingHorn2(vHorn2);
       this->ConstructLBNEShieldingBetweenHorns(tunnel);
     } else {
       this->ConstructLBNFShielding(tunnel);
     }
   } 
  }
  else {

    DropMarsTargetHorns(tunnel);

  }
  
   fHasBeenConstructed = true;
   
   if(fPlacementHandler->GetWriteGDMLFile()) {
     G4GDMLParser Parser;
     Parser.Write("g4lbnf.gdml",tunnel);
   }

  return fRock;
}
void LBNEDetectorConstruction::ConstructLBNEHadronAbsorber(G4VPhysicalVolume *mother)
{

   G4cout << "Importing hadron absorber gdml file... " << G4endl;
   
   G4String filename(fPlacementHandler->GetAbsorberGDMLFilename());
   std::ifstream gdmlfile(filename.c_str());
   if (!gdmlfile.is_open()) {
     std::string mess(" AbsorberGDML file "); 
     mess += filename + G4String(" could not be found \n");
     G4Exception("LBNEDetectorConstruction::ConstructLBNEHadronAbsorber", " ", 
                 FatalErrorInArgument, mess.c_str());
     return; // perfunctory. 
     
   } else {
     gdmlfile.close();
   }
   G4GDMLParser parser;
   parser.Read( filename );
   //std::cerr << " And stop after parsing the gdml file " << std::endl; exit(2);
   
     G4LogicalVolume *topAbs = parser.GetVolume( "TOP" );
     // We dump the volume hierarchy.  Hoopefully not too deep,
     /*
     for (int i=0; i != topAbs->GetNoDaughters(); ++i) {
        G4VPhysicalVolume *pVol = topAbs->GetDaughter(i);
        G4LogicalVolume *lVol = pVol->GetLogicalVolume();
        std::cerr << " Top level daughter # " << i << " Name " << lVol->GetName() 
                  << " at " << pVol->GetObjectTranslation() << std::endl;
           if (lVol->GetName().find("Airbox") != std::string::npos) {
             G4Box *aBox = static_cast<G4Box*>(lVol->GetSolid());
             std::cerr << " Airbox size " << aBox->GetXHalfLength() 
                       << " / " << aBox->GetYHalfLength() << " / " << aBox->GetZHalfLength() << std::endl; 
           }
       for (int ii=0; ii != lVol->GetNoDaughters(); ++ii) {
          G4VPhysicalVolume *pVol2 = lVol->GetDaughter(ii);
          G4LogicalVolume *lVol2 = pVol2->GetLogicalVolume();
          std::cerr << "  .. 2nd level daughter # " << ii << " Name " << lVol2->GetName()
           << " at " << pVol2->GetObjectTranslation() << std::endl;
           if (lVol2->GetName().find("Airbox") != std::string::npos) {
             G4Box *aBox = static_cast<G4Box*>(lVol2->GetSolid());
             std::cerr << " Airbox size " << aBox->GetXHalfLength() 
                       << " / " << aBox->GetYHalfLength() << " / " << aBox->GetZHalfLength() << std::endl; 
           }
         for (int iii=0; iii != lVol2->GetNoDaughters(); ++iii) {
           G4VPhysicalVolume *pVol3 = lVol2->GetDaughter(iii);
           G4LogicalVolume *lVol3 = pVol3->GetLogicalVolume();
           std::cerr << "  ... 3rd level daughter # " << iii << " Name " << lVol3->GetName() 
            << " at " << pVol3->GetObjectTranslation() << std::endl;
           for (int i4=0; i4 != lVol3->GetNoDaughters(); ++i4) {
            G4VPhysicalVolume *pVol4 = lVol3->GetDaughter(i4);
            G4LogicalVolume *lVol4 = pVol4->GetLogicalVolume();
            std::cerr << "  .... 4rth level daughter # " << i4 << " Name " << lVol4->GetName()
             << " at " << pVol4->GetObjectTranslation()  << std::endl;
            for (int i5=0; i5 != lVol4->GetNoDaughters(); ++i5) {
              G4VPhysicalVolume *pVol5 = lVol4->GetDaughter(i5);
              G4LogicalVolume *lVol5 = pVol5->GetLogicalVolume();
              std::cerr << "  ..... 5rth level daughter # " << i5 << " Name " << lVol5->GetName() 
              << " at " << pVol5->GetObjectTranslation()  << std::endl;
              for (int i6=0; i6 != lVol5->GetNoDaughters(); ++i6) {
                G4VPhysicalVolume *pVol6 = lVol5->GetDaughter(i6);
                G4LogicalVolume *lVol6 = pVol6->GetLogicalVolume();
                std::cerr << "  ...... 6rth level daughter # " << i6 << " Name " << lVol6->GetName() << std::endl;
              }
            }
          }
        }
      }
    }
*/       
    
     const G4Box *topSol = static_cast<const G4Box *>(topAbs->GetSolid());
//     const double marsTopWidth = topSol->GetYHalfLength(); // Used for debugging only, so comment out for now.. 
//     const double marsTopHeight = topSol->GetXHalfLength();
//     const double marsTopLength = topSol->GetZHalfLength();
     std::cerr << " Dimension of top level Hadron absorber MARS container, X " << topSol->GetXHalfLength() << 
         " Y  "  << topSol->GetYHalfLength() <<   " Z  "  << topSol->GetZHalfLength() << std::endl;
     std::cerr << " Number of daughters for TOP " << topAbs->GetNoDaughters() << std::endl;
     double maxHalfHeight = -1.0;
     double maxHalfWidth = -1.0;
     double maxHalfLength = -1.0;
     for (int i=0; i != topAbs->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = topAbs->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
       std::cerr << " Daughther " << lVol->GetName();
       const G4Box *aBox = static_cast<const G4Box *>(lVol->GetSolid());
       G4ThreeVector loc = pVol->GetObjectTranslation();
       std::cerr << " at MARS coordinates " << loc[0] << ", " <<loc[1] << ", " << loc[2] << 
                     " zLength " << 2.0*aBox->GetZHalfLength() << std::endl;
       // Compute the maximum height, width.  Note the confusion about  X and Y X is up, vertical, in MArs  
       if ((std::abs(loc[2]) + aBox->GetZHalfLength()) > maxHalfLength)
           maxHalfLength = std::abs(loc[2]) + aBox->GetZHalfLength();
       if ((std::abs(loc[1]) + aBox->GetYHalfLength()) >  maxHalfWidth)
           maxHalfWidth = std::abs(loc[1]) + aBox->GetYHalfLength(); // Width is along X G4lbne orientation, which Y MARS 
       if ((std::abs(loc[0]) + aBox->GetXHalfLength()) >  maxHalfHeight)
           maxHalfHeight = std::abs(loc[0]) + aBox->GetXHalfLength();
            // Height is along Y G4lbne orientation, which is negative X in MARS 
     }
     maxHalfHeight += 5.0*CLHEP::cm;
     maxHalfWidth += 5.0*CLHEP::cm;
     maxHalfLength += std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) + 5.0*CLHEP::cm;;
     std::cerr << " Container volume for Hadron absorber, 1/2 width " 
               << maxHalfWidth << " 1/2 Height " << maxHalfHeight 
             << " 1/2 length " << maxHalfLength << std::endl;
     G4Box *aHABoxTop = new G4Box(G4String("HadronAbsorberTop"), maxHalfWidth, maxHalfHeight, maxHalfLength);
     G4LogicalVolume *aHATopL = 
        new G4LogicalVolume(aHABoxTop, G4Material::GetMaterial("Air"), G4String("HadronAbsorberTop"));
      const LBNEVolumePlacementData *plDecayPipe = 
         fPlacementHandler->Find(G4String("HadronAbsorber"), G4String("DecayPipeHall"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEHadronAbsorber"));
     const double zzz = maxHalfLength + plDecayPipe->fParams[2]/2 + plDecayPipe->fPosition[2] + 
         std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) +  5.0*CLHEP::cm;
                
     std::cerr << " half length Decay Pipe " << plDecayPipe->fParams[2]/2 
               << " Position in tunnel (center) " << plDecayPipe->fPosition[2] 
               << " half width decay pipe " << plDecayPipe->fParams[0]/2 << " 1/2 height " << plDecayPipe->fParams[1]/2
               << " ZPosAbs " << zzz <<  std::endl;
// One must adjust the height of the Absorber such the the beam line crosses the center of the 3 rth Core block in it's 
// center.. 
//
     const double yyy = 2.3*CLHEP::m; // Adjusted by hand, looking Geantino's Acuracy: for crossing of block 3 and 4 
                               // Average of crossing -4 mm . Min Max average is 0.6 mm.  
     G4ThreeVector posTopHA(0., yyy, zzz);
     
     new G4PVPlacement(&fRotBeamlineAngle, posTopHA, aHATopL, "HadronAbsorberTop", 
                     mother->GetLogicalVolume(), false, 1, true);
     
     G4RotationMatrix *marsRot = new G4RotationMatrix;
     marsRot->rotateZ(-M_PI/2.);
     for (int i=0; i != topAbs->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = topAbs->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
//       const G4Box *aBox = static_cast<const G4Box *>(lVol->GetSolid());
       G4ThreeVector loc = pVol->GetObjectTranslation();
//       const double yyy = loc[0]  - marsTopHeight + maxHalfHeight;  
       double yyyI = loc[0] ;  // Up to a sign!!! 
       const double xxx = loc[1]; // Y G4LBNE = -X Mars. Was centered in MARS, set to 0., o.k.
         // X G4LBNE = Y Mars. !! Not centered in Mars! Perhaps, ned a shift due to the 
        // the different size of the mother volume  
//       const double zzz = loc[2] - marsTopLength + maxHalfLength;
       double zzzI = loc[2] - 27.9*CLHEP::m; // Setting up this last to avoid air to air volume overlap. 
       if (lVol->GetName().find("AH_top") != std::string::npos) {
          yyyI += 0.1*CLHEP::mm; // To avoid clash.. Cosmetic 
       // Note: for release v3r0px up to v3r0p10, this was set to 27.9 m 
          zzzI = loc[2] - 27.0*CLHEP::m; // Also to avoid clash.  Note that this volume is not 
                                 // in the pion beam, it is up ~ 8 m. above the beam line. 
       }
       G4ThreeVector posTmp(xxx, yyyI, zzzI);
//       std::cerr << " Placing volume " << lVol->GetName() << " at " << posTmp << " 1/2 sizes (G4 coord)  " 
//        << aBox->GetYHalfLength() << " , " << aBox->GetXHalfLength() << " , " 
//      <<  aBox->GetZHalfLength() << std::endl;
//       std::cerr << "  .... Extend in X " << posTmp[0] - aBox->GetYHalfLength() 
//               << " to " << posTmp[0] + aBox->GetYHalfLength() << std::endl;
//       std::cerr << "  .... Extend in Y " << posTmp[1] - aBox->GetXHalfLength() 
//               << " to " << posTmp[1] + aBox->GetXHalfLength() << std::endl;
//       std::cerr << "  .... Extend in Z " << posTmp[2] - aBox->GetZHalfLength() 
//               << " to " << posTmp[2] + aBox->GetZHalfLength() << std::endl;
       new G4PVPlacement(marsRot, posTmp, lVol, lVol->GetName() + std::string("_P"), aHATopL, false, 1, true);
     }


      //---Tracking planes--------
       //-- to place stuff in the absorber, we must first get the
       // muon alcove region.

     
       for (int i=0; i != topAbs->GetNoDaughters(); ++i) {
         G4VPhysicalVolume *pVol1 = topAbs->GetDaughter(i);
         G4LogicalVolume *lVol1 = pVol1->GetLogicalVolume();

        
         for (int ii=0; ii != lVol1->GetNoDaughters(); ++ii) {
           G4VPhysicalVolume *pVol2 = lVol1->GetDaughter(ii);
           G4LogicalVolume *lVol2 = pVol2->GetLogicalVolume();
           std::string aVolNameTmp2(lVol2->GetName());
           if (aVolNameTmp2 != std::string("AH_Muon_alkair")) continue;
           std::cerr << " Found AH_Muon_alk again  ... " << std::endl;

           G4Box *detSolid = new G4Box("detSolid",2.5*CLHEP::m,2.5*CLHEP::m,1.0*CLHEP::cm);
      
           TrackingPlaneLogical= new G4LogicalVolume(detSolid , G4Material::GetMaterial("Air"), "detLogical", 0,0,0);

           G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =0.;
           posTmp[2] = -6.0*CLHEP::m; // placed right after absorber

           std::cerr << " Placing test planes in absorber hall" << std::endl;                    
           new G4PVPlacement((G4RotationMatrix *) 0,posTmp, TrackingPlaneLogical, 
                          "trackPln1",lVol2, false, 1, true);          

           G4ThreeVector posTmp2; posTmp2[0] =0.; posTmp2[1] =0.;
           posTmp2[2] = -3.3*CLHEP::m; // placed after blue block     
           new G4PVPlacement((G4RotationMatrix *) 0,posTmp2, TrackingPlaneLogical, 
                          "trackPln2",lVol2, false, 1, true);          
           std::cerr << " Top level daughter # " << i << " Name " << lVol1->GetName();                   
         }
       }


     //
     // We now clear the MARS top absorbers, simplify the geometry search 
     //
     topAbs->ClearDaughters();

     // save logical volume
     HadrBox_log = aHATopL;
}

void LBNEDetectorConstruction::ConstructLBNEHadronAbsorberSimple(G4VPhysicalVolume *mother)
{

   //*************************** Start creating ********************************
  //Define dynamically measure for Steel and Concrete volumes 
  G4double SteelWidth=this->GetDwStrAbsSteelWidth();
  G4double ConcWidth=this->GetDwStrAbsConcreteWidth();
  
   //Container Absorber(Air) (those measures are hard-coded but it would be better to define them dynamically)
   G4double air_box_x = 16.223*CLHEP::m;
   G4double air_box_y = 19.355*CLHEP::m;
   G4double air_box_z = 20*CLHEP::m;


   G4Box *air_Box = new G4Box("AirBox",air_box_x/2,air_box_y/2,air_box_z/2);
   G4ThreeVector air_posBox; air_posBox[0] =0.; air_posBox[1] =0.45*CLHEP::m; air_posBox[2] = 231.0*CLHEP::m;

   G4LogicalVolume *air_logBox = 
     new G4LogicalVolume(air_Box, G4Material::GetMaterial("Air"), G4String("LogicalVAirBox"),0,0,0);
   
   //fRotBeamlineAngle is the rotational matrix that defines the rotation of the beam; if we want to tilt the absorber we use it, otherwise we use "(G4RotationMatrix *)0"
     new G4PVPlacement(&fRotBeamlineAngle,air_posBox,air_logBox, 
       "PhysAirBox",mother->GetLogicalVolume(), false, 1, true);
   
   air_logBox->SetVisAttributes(GetMaterialVisAttrib("Air"));
  
   //Set Attributes of Container(Air)
   /*G4VisAttributes *airBox_Att = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
   airBox_Att->G4VisAttributes::SetForceWireframe(true);
   air_logBox->SetVisAttributes(airBox_Att);*/
  

   //------------------------------ Define Absorber Volumes ------------------------------------------------
   
   //~~~~~ Absorber (concrete volume) ~~~~~~~
   G4double AH_ConcreteX = 10.23*CLHEP::m;
   G4double AH_ConcreteY = 12.125*CLHEP::m;
   G4double AH_ConcreteZ = 6.24*CLHEP::m + SteelWidth + ConcWidth;

   G4Box *ConcreteBox1= new G4Box("ConcreteBox1",AH_ConcreteX/2,AH_ConcreteY/2,AH_ConcreteZ/2);

   G4Box *box1 = new G4Box("box1",3.5/2*CLHEP::m, 4.05/2*CLHEP::m,5*CLHEP::m );
   G4SubtractionSolid *AH_Concrete= 
     new G4SubtractionSolid("AH_Concrete",ConcreteBox1,box1,0,G4ThreeVector(0.,0.1*CLHEP::m,-4.37*CLHEP::m-AH_ConcreteZ/2));
   
   G4ThreeVector posAH_Concrete; posAH_Concrete[0] =0.; posAH_Concrete[1] =0.4*CLHEP::m; posAH_Concrete[2] =AH_ConcreteZ/2-air_box_z/2;
  
   G4ThreeVector posAH_AirBox1; posAH_AirBox1[0] =0.; posAH_AirBox1[1] =0.1*CLHEP::m; posAH_AirBox1[2] = (0.63/2)*CLHEP::m-AH_ConcreteZ/2;
  
   G4LogicalVolume *AH_ConcreteLV= 
     new G4LogicalVolume(AH_Concrete, G4Material::GetMaterial("Concrete"), G4String("LV_AHConcrete"),0,0,0);
   new G4PVPlacement((G4RotationMatrix *) 0,posAH_Concrete,AH_ConcreteLV, 
   "PV_AHConcrete",air_logBox, false, 1, true);
   
   //AH_ConcreteLV->SetVisAttributes(GetMaterialVisAttrib("Concrete"));

   //Set Attributes Absorber(concrete)
   G4VisAttributes *AHConcrete_Att = new G4VisAttributes(G4Colour(G4Colour::Magenta()));
   AHConcrete_Att->G4VisAttributes::SetForceWireframe(true);
   AH_ConcreteLV->SetVisAttributes(AHConcrete_Att);
   //----------------------------end Vol #1-----------------


   //~~~~~ Daughter 1 (steel volume) ~~~~~~~
   
   G4double AH_ShieldStlX = 5.46*CLHEP::m;
   G4double AH_ShieldStlY = 5.59*CLHEP::m;
   G4double AH_ShieldStlZ = 5.41*CLHEP::m+SteelWidth;

   G4Box *AH_ShieldStl= new G4Box("AH_ShieldStl",AH_ShieldStlX/2,AH_ShieldStlY/2,AH_ShieldStlZ/2);
   

   G4ThreeVector posAH_ShieldStl; posAH_ShieldStl[0] =0.; 
   posAH_ShieldStl[1] =0.; posAH_ShieldStl[2] =AH_ShieldStlZ/2-AH_ConcreteZ/2+0.63*CLHEP::m; //0.63m is the width of concrete front hole

   G4LogicalVolume *AH_ShieldStlLV= 
     new G4LogicalVolume(AH_ShieldStl,/*G4Material::GetMaterial("Air")*/G4Material::GetMaterial("Slab_Stl"), G4String("LV_AHShieldStl"),0,0,0);
   new G4PVPlacement((G4RotationMatrix *) 0,posAH_ShieldStl,AH_ShieldStlLV, 
                     "PV_AHShieldStl",AH_ConcreteLV, false, 1, true);
   // AH_ShieldStlLV->SetVisAttributes(GetMaterialVisAttrib("Iron"));

   //Set Atrributes Absorber(steel)
   G4VisAttributes *AHShieldStl_Att = new G4VisAttributes(G4Colour(0.,1.,0.));
   AHShieldStl_Att->G4VisAttributes::SetForceWireframe(true);
   AH_ShieldStlLV->SetVisAttributes(AHShieldStl_Att);
   //----------------------------end Vol #2-----------------
   //--- volume 3 is no longer used.

   //~~~~~ GrandDaughter 1 (Al hole volume) ~~~~~~~

   G4double AH_AlholeX = 1.96*CLHEP::m;
   G4double AH_AlholeY = 2.18*CLHEP::m;
   G4double AH_AlholeZ = 1.83*CLHEP::m;

   G4Box *AH_Alhole= new G4Box("AH_Alhole",AH_AlholeX/2,AH_AlholeY/2,AH_AlholeZ/2);

   //Create Polycone volume
   G4double al = (-1.83/2)*CLHEP::m;
   G4double posZ [8];posZ[0]=al; posZ[1]=al+0.305*CLHEP::m; 
   posZ[2]=al+0.305*CLHEP::m; posZ[3]=al+0.610*CLHEP::m; posZ[4]=al+0.610*CLHEP::m;
   posZ[5]=0.305*CLHEP::m;posZ[6]=0.305*CLHEP::m;
   posZ[7]=0.915*CLHEP::m;
   G4double Irad [8];Irad[0]=0.; Irad[1]=0.; Irad[2]=0.;Irad[3]=0.;Irad[4]=0.;Irad[5]=0.;Irad[6]=0.;Irad[7]=0.;
   G4double Orad [8];Orad[0]=0.365*CLHEP::m; Orad[1]=0.365*CLHEP::m; 
   Orad[2]=0.615*CLHEP::m;Orad[3]=0.615*CLHEP::m;Orad[4]=0.5*CLHEP::m;
   Orad[5]=0.5*CLHEP::m;Orad[6]=0.365*CLHEP::m;Orad[7]=0.365*CLHEP::m;
   // Convert to an array of z, r, wrapping on itself (old signature gives a run time exception..)  
   std::vector<double> zzPolycone(16, 0.);
   std::vector<double> rrPolycone(16, 0.);
   for (size_t kPoly=0; kPoly != 8; kPoly++) {
     zzPolycone[kPoly] = posZ[kPoly];
     zzPolycone[kPoly+8] = posZ[7-kPoly];
     rrPolycone[kPoly] = Irad[kPoly];
     rrPolycone[kPoly+8] = Orad[kPoly];
   }
    
   G4Polycone* pcone = new G4Polycone("pcone", 0.,2*M_PI*CLHEP::rad,8, &rrPolycone[0], &zzPolycone[0]);

   G4ThreeVector posAH_Alhole; posAH_Alhole[0] =0.; posAH_Alhole[1] =0.; posAH_Alhole[2] =AH_AlholeZ/2-AH_ShieldStlZ/2;

   G4LogicalVolume *AH_AlholeLV= 
     new G4LogicalVolume(AH_Alhole, G4Material::GetMaterial("Aluminum"), G4String("LV_AHAlhole"),0,0,0);
   new G4PVPlacement((G4RotationMatrix *) 0,posAH_Alhole,AH_AlholeLV, 
                     "PV_AHAlhole",AH_ShieldStlLV, false, 1, true);
   
    
    G4LogicalVolume *AH_AirholeLV= 
     new G4LogicalVolume(pcone, G4Material::GetMaterial("Air"), G4String("LV_AHAirhole"),0,0,0);
    new G4PVPlacement((G4RotationMatrix *) 0,G4ThreeVector(0,0,0),AH_AirholeLV, 
    "PV_AHAirhole",AH_AlholeLV, false, 1, true);
 

   //Set Atrributes Absorber(Al hole)
   G4VisAttributes *AHAlhole_Att = new G4VisAttributes(G4Colour(0.,192.,255.));
   AHAlhole_Att->G4VisAttributes::SetForceWireframe(true);
   AH_AlholeLV->SetVisAttributes(AHAlhole_Att);
   //----------------------------end Vol #4-----------------


   //~~~~~ GrandDaughter 2 (Al core volume) ~~~~~~~
   G4double AH_AlcoreX = 1.52*CLHEP::m;
   G4double AH_AlcoreY = 1.52*CLHEP::m;
   G4double AH_AlcoreZ = 2.46*CLHEP::m;

   G4Box *AH_Alcore= new G4Box("AH_Alcore",AH_AlcoreX/2,AH_AlcoreY/2,AH_AlcoreZ/2);
   

   G4ThreeVector posAH_Alcore; posAH_Alcore[0] =0.; 
   posAH_Alcore[1] =-0.2*CLHEP::m; posAH_Alcore[2] =AH_AlcoreZ/2-AH_ShieldStlZ/2+AH_AlholeZ;

   G4LogicalVolume *AH_AlcoreLV= 
     new G4LogicalVolume(AH_Alcore,G4Material::GetMaterial("Aluminum"), G4String("LV_AHAlcore"),0,0,0);
   new G4PVPlacement((G4RotationMatrix *) 0,posAH_Alcore,AH_AlcoreLV, 
                     "PV_AHAlcore",AH_ShieldStlLV, false, 1, true);

   //Set Atrributes Absorber(Al hole)
   G4VisAttributes *AHAlcore_Att = new G4VisAttributes(G4Colour(0.,0.,1.));
   AHAlcore_Att->G4VisAttributes::SetForceWireframe(true);
   AH_AlcoreLV->SetVisAttributes(AHAlcore_Att);
   //----------------------------end Vol #5-----------------


 //~~~~~ GrandDaughter 3 (Steel absorber volume) ~~~~~~~
   G4double AH_SteelX = 1.52*CLHEP::m;
   G4double AH_SteelY = 1.52*CLHEP::m;
   G4double AH_SteelZ = 1.22*CLHEP::m;

   G4Box *AH_Steel= new G4Box("Ah_Steel",AH_SteelX/2,AH_SteelY/2,AH_SteelZ/2);

   G4ThreeVector posAH_Steel; posAH_Steel[0] =0.; 
   posAH_Steel[1] =-0.2*CLHEP::m; posAH_Steel[2] =AH_SteelZ/2-AH_ShieldStlZ/2+(AH_AlholeZ+AH_AlcoreZ);

   G4LogicalVolume *AH_SteelLV= 
     new G4LogicalVolume(AH_Steel, G4Material::GetMaterial("Slab_Stl"), G4String("LV_AHSteel"),0,0,0);
   new G4PVPlacement((G4RotationMatrix *) 0,posAH_Steel,AH_SteelLV, 
                     "PV_AHSteel",AH_ShieldStlLV, false, 1, true);

   //Set Atrributes Absorber(steel)
   G4VisAttributes *AHSteel_Att = new G4VisAttributes(G4Colour(1.,0.,0.));
   AHSteel_Att->G4VisAttributes::SetForceWireframe(true);
   AH_SteelLV->SetVisAttributes(AHSteel_Att);
   //----------------------------end Vol #6-----------------


   //~~~~~Tracking Planes ~~~~~~~
   //Position of planes is dynamically coded so that they are 2.7m apart from each other and 4m away from the end of the absorber and surround the Blueblock. 
   //If you change geometry of "simple detector" consider changing those measures (2.7m and 4m). 
    G4Box *detSolid = new G4Box("detSolid",2.5*CLHEP::m,2.5*CLHEP::m,1.0*CLHEP::cm);
      
   TrackingPlaneLogical= new G4LogicalVolume(detSolid , G4Material::GetMaterial("Air"), "detLogical", 0,0,0);

   G4ThreeVector posTmp; posTmp[0] =0.; posTmp[1] =-0.45*CLHEP::m;
   posTmp[2] = ConcWidth+4*CLHEP::m; // placed right after absorber
   std::cerr << " Placing test planes in absorber hall" << std::endl;                    
   new G4PVPlacement((G4RotationMatrix *) 0,posTmp, TrackingPlaneLogical, 
                          "trackPln1",air_logBox, false, 1, true);          

   G4ThreeVector posTmp2; posTmp2[0] =0.; posTmp2[1] =-0.45*CLHEP::m;
   posTmp2[2] =posTmp[2]+2.7*CLHEP::m; // placed after blue block     
   new G4PVPlacement((G4RotationMatrix *) 0,posTmp2, TrackingPlaneLogical, 
   "trackPln2",air_logBox, false, 1, true);

  // G4VPhysicalVolume *trackerPhysBox =
  //  new G4VPhysicalVolume((G4RotationMatrix *)0,posBox,"TrackerBox",logBox, mother->GetLogicalVolume(), false, 1, true);         
   
           //**********************DONE**************************
        
}
//end simpler Absorber


void LBNEDetectorConstruction::ConstructLBNEShieldingHorn1(G4VPhysicalVolume *mother) {

   const double in = 25.4*CLHEP::mm;
//
// Install steel shielding around the Horn1 and the target. 
//
// Based on Docdb 5339, page 30 and 31. 
// See also drawing 2251.000-ME-487105 
// 
// We simply install steel, assumed low carbon...
//
   const double maxR1 =  fPlacementHandler->GetMaxRadiusMotherHorn1();
         fPlacementHandler->Find(G4String("LBNFShielding"), G4String("Horn1PolyM1"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingHorn2"));
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   const LBNEVolumePlacementData *plTop = 
         fPlacementHandler->Find(G4String("ShieldingHorn1"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingHorn1"));
                                 
   const double horn1HeightAtMCZERO = 66.7*in;                           
//
// Bottom 
//
   G4String bName = G4String("ShieldingHorn1Bottom"); 
   const double bWidth = 208*in;
   const double bHeight = 72.0*in;
   const double bLength = plTop->fParams[2] - 2.*bHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *bBox = new G4Box(bName, bWidth/2., bHeight/2., bLength/2.);
   G4LogicalVolume *bLVol = new G4LogicalVolume(bBox, G4Material::GetMaterial(std::string("Slab_Stl")), bName); 
   G4ThreeVector posTmp(0., 0., 0.); 
   posTmp[1] = -horn1HeightAtMCZERO  - bHeight/2.; 
   // MCZERO is G4 coordinate 0.0, the center of tunnel, so the above number needs to be corrected for the 
   // shift in the center of the mother volume  
   const double zShift = -1.0*plTop->fPosition[2];
   const double yCorr = -1.0*zShift*std::sin(fBeamlineAngle);
   posTmp[1]  -= yCorr;
//   std::cerr << "ConstructLBNEShieldingHorn1 .... "<< bName << " zShift bottom plate " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr << std::endl;
//           
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, bLVol, bName + std::string("_P"),
                       mother->GetLogicalVolume(), false, 1, true);
//
// Sides 
//
   G4String sName = G4String("ShieldingHorn1Side"); 
   const double sWidth = 65.0*in; // Reduced from 72 inches down to 65 to give room for Norn1 Mother volume 
                                  // for Optimized design.
   const double sHeight = 270.8*in - 72.0*in - 4.0*CLHEP::cm; // oversize a bit, but no matter.. 
   const double sLength = plTop->fParams[2] - 2.*sHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *sBox = new G4Box(sName, sWidth/2., sHeight/2., sLength/2.);
   G4LogicalVolume *sLVol = new G4LogicalVolume(sBox, G4Material::GetMaterial(std::string("Slab_Stl")), sName); 
   posTmp[0] = -maxR1 - 15.0*CLHEP::cm - sWidth/2.;
   posTmp[1] = sHeight/2. - horn1HeightAtMCZERO + 2.0*CLHEP::cm; 
   // MCZERO is G4 coordinate 0.0, the center of tunnel, so the above number needs to be corrected for the 
   // shift in the center of the mother volume  
   posTmp[1]  -= yCorr;
//   std::cerr << "ConstructLBNEShieldingHorn1 .... "<< bName << " zShift side plate, left  " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr 
//           << std::endl << " .......  X = " << posTmp[0] << " Y " << posTmp[1] << std::endl;
             
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_LeftP"),
                       mother->GetLogicalVolume(), false, 1, true);
   posTmp[0] = maxR1 + 15.0*CLHEP::cm + sWidth/2.;
//   std::cerr << "ConstructLBNEShieldingHorn1 .... "<< bName << " zShift side plate, Right  " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr 
//           << std::endl << " .......  X = " << posTmp[0] << " Y " << posTmp[1] << std::endl;
             
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_RightP"),
                       mother->GetLogicalVolume(), false, 2, true);
//
// Top. Simply put about 2m of steel 
//
   G4String tName = G4String("ShieldingHorn1Top"); 
   const double tWidth = 208*in;
   const double tHeight = 94.0*in;
   const double tLength = plTop->fParams[2] - 2.*tHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *tBox = new G4Box(tName, tWidth/2., tHeight/2., tLength/2.);
   G4LogicalVolume *tLVol = new G4LogicalVolume(tBox, G4Material::GetMaterial(std::string("Slab_Stl")), tName); 
   
   posTmp[0] =0.; 
   posTmp[1] += sHeight/2. + tHeight/2. + 4.0*CLHEP::cm; 
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, tLVol, bName + std::string("_P"),
                       mother->GetLogicalVolume(), false, 1, true);
 
   
}

//Construct the Tracking plane at the end of Horn 1 Amit Bashyal
void LBNEDetectorConstruction::ConstructLBNEHorn1TrackingPlane(G4VPhysicalVolume *mother) {

   const G4String nameM = mother->GetLogicalVolume()->GetName();
   if(nameM != G4String("Tunnel")){
     std::cerr
     <<"Unexpected Mother Volume in LBNEDetectorConstruction::ConstructLBNEHorn1TrackingPlane(G4VPhysicalVolume *mother)"<<std::endl;
      exit(2);
   }
         fPlacementHandler->Find(G4String("Horn1TrackingPlane"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEHorn1TrackingPlane"));
   const LBNEVolumePlacementData *plH1 = 
         fPlacementHandler->Find(G4String("Horn1TrackingPlane"), G4String("TargetHallAndHorn1"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEHorn1TrackingPlane"));
   const double z = plH1->fPosition[2] + plH1->fParams[2]/2.0 + 2.0*CLHEP::cm; //Position set at 2 cm after the end of Horn 1.

   const double widthTop =  80.00*CLHEP::cm; //Set the Dimension such that it just covers the mouth of Horn 1 end.
   const double heightTop = 80.00*CLHEP::cm;
   
//   G4String topNameBetween = G4String("Horn1TrackingPlane");
   G4Box *topBox = new G4Box("Horn1TrackingPlane", widthTop/2., heightTop/2.,0.5*CLHEP::mm); //Thickness of  0.5mm
   TrackingPlaneH1Logical = new G4LogicalVolume(topBox, G4Material::GetMaterial(std::string("Air")),"Horn1TrackingPlaneLogical"); 
   G4ThreeVector posTopLevel (0., 0., z);
std::cerr
      <<"The value of Z in Horn1Plane is:" << z<<"   "<< plH1->fPosition[2]<<"   "<<plH1->fParams[2]<<std::endl;
std::cerr
     <<"Width of Horn1TrackingPlane is:" <<widthTop/2<<" Height of Horn1TrackingPlane is: "<<heightTop/2<<std::endl;
   G4PVPlacement *vTopLevel = new G4PVPlacement((G4RotationMatrix *) 0, posTopLevel, 
                                                TrackingPlaneH1Logical,"Horn1TrackingPlane_P",
                                                mother->GetLogicalVolume(), false, 1, true);
    if(vTopLevel == 0){
      std::cerr << " Can't place Horn1 tracking plane.. Unexpected.. " << std::endl; // This error should never occur.. 
    }
                                                
}



void LBNEDetectorConstruction::ConstructLBNEShieldingHorn2(G4PVPlacement *mother) {

   const double in = 25.4*CLHEP::mm;
//
// Install steel shielding around the Horn2 and the target. 
//
// Based on Docdb 5339, page 30 and 31. 
// See also drawing 2251.000-ME-487105 
// 
// We simply install steel, assumed low carbon...
//
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   const LBNEVolumePlacementData *plTop = 
         fPlacementHandler->Find(G4String("ShieldingHorn2"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingHorn2"));
                                 
   const LBNEVolumePlacementData *plHorn2 = 
         fPlacementHandler->Find(G4String("LBNFShielding"), G4String("Horn2TopLevel"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingHorn2"));
//   std::cerr << " Horn2TopLevel params " << plHorn2->fParams[0] << " / " 
//             << plHorn2->fParams[1] <<        " / " << plHorn2->fParams[2] <<   std::endl;  
//   std::cerr << " Horn2TopLevel position " << plHorn2->fPosition[0] << " / " 
//             << plHorn2->fPosition[1] <<      " / " << plHorn2->fPosition[2] <<   std::endl;  
   const double horn2Height = (37. + 17.3)*in; // ???? Checked final position of the bottom slab 
   // from drawing of page 30 of LBNE Doc 5339-v5 Vertical distance between the top of the bottom shielding 
   // and the beam line, at the geometricla center of Horn2.                             
//
// Bottom 
//
   G4String bName = G4String("ShieldingHorn2Bottom"); 
   const double bWidth = 168*in;
   const double bHeight = 52.0*in;
   const double bLength = plTop->fParams[2] - 2.*bHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *bBox = new G4Box(bName, bWidth/2., bHeight/2., bLength/2.);
   G4LogicalVolume *bLVol = new G4LogicalVolume(bBox, G4Material::GetMaterial(std::string("Slab_Stl")), bName); 
   G4ThreeVector posTmp(0., 0., 0.); 
   posTmp[1] = -horn2Height  - bHeight/2.; 
   // MCZERO is G4 coordinate 0.0, the center of tunnel, so the above number needs to be corrected for the 
   // shift in the center of the mother volume  
//   const double zShift = -1.0*plTop->fPosition[2];
//   const double yCorr = -1.0*zShift*std::sin(fBeamlineAngle); v3r0p8
   const double yCorr = 0; // directly measured of the drawing 
   posTmp[1]  += yCorr;
//  std::cerr << "ConstructLBNEShieldingHorn2 .... "<< bName << " zShift bottom plate " << zShift 
//        << " thus, y position " << posTmp[1] << " coorection " << yCorr <<  " And quit " << std::endl;
//        exit(2);
//           
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, bLVol, bName + std::string("_P"),
                       mother->GetLogicalVolume(), false, 1, true);
//
// Sides 
//
//
// December 2015: Upgrade to v4.10, for which volume overlap problem cause real trouble. 
// So, if we rescale Horn2 radially, we may not have enough room to place shielding on the sides. 
// Note: this is still the old Chase width..
// If LBNF Polycones horns are used, then the chase is extended 

   G4String sName = G4String("ShieldingHorn2Side"); 
   double sWidth = 52.0*in;
   double sHeight = 256.8*in - 52.0*in - 4.0*CLHEP::cm; // oversize a bit, but no matter.. 
   double transvShift = plHorn2->fParams[1] + 1.0*CLHEP::cm;
//   const double aHorn2RadialRescale = fPlacementHandler->GetHorn2RadialRescale();
//   if ((aHorn2RadialRescale - 1.) > 0.) {
//     sWidth -= (aHorn2RadialRescale - 1.)*20*in ; // Approximate radius of the Horn2 
//     std::cerr << " Reducing the width of the side shielding at Horn2 from 1.32 m inches to " 
//               << sWidth/CLHEP::m << " m " << std::endl;                                        
//     sHeight = 256.8*in - sWidth - 4.0*CLHEP::cm; 
//   }
   const double sLength = plTop->fParams[2] - 2.*sHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *sBox = new G4Box(sName, sWidth/2., sHeight/2., sLength/2.);
   G4LogicalVolume *sLVol = new G4LogicalVolume(sBox, G4Material::GetMaterial(std::string("Slab_Stl")), sName); 
   posTmp[0] = -transvShift - sWidth/2.;
   posTmp[1] = sHeight/2. - horn2Height + 2.0*CLHEP::cm; 
   // MCZERO is G4 coordinate 0.0, the center of tunnel, so the above number needs to be corrected for the 
   // shift in the center of the mother volume  
     posTmp[1]  -= yCorr;
//   std::cerr << "ConstructLBNEShieldingHorn2 .... "<< bName << " zShift side plate, left  " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr 
//           << std::endl << " .......  X = " << posTmp[0] << " Y " << posTmp[1] << std::endl;
             
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_LeftP"),
                       mother->GetLogicalVolume(), false, 1, true);
    posTmp[0] = transvShift + sWidth/2.;
//   std::cerr << "ConstructLBNEShieldingHorn2 .... "<< bName << " zShift side plate, Right  " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr 
//           << std::endl << " .......  X = " << posTmp[0] << " Y " << posTmp[1] << std::endl;
             
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_RightP"),
                       mother->GetLogicalVolume(), false, 2, true);
//
// Top. Simply put about 2m of steel 
//
   G4String tName = G4String("ShieldingHorn2Top"); 
   const double tWidth = 168*in;
   const double tHeight = 94.0*in;
   const double tLength = plTop->fParams[2] - 2.*tHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *tBox = new G4Box(tName, tWidth/2., tHeight/2., tLength/2.);
   G4LogicalVolume *tLVol = new G4LogicalVolume(tBox, G4Material::GetMaterial(std::string("Slab_Stl")), tName); 
   
   posTmp[0] =0.; 
   posTmp[1] += sHeight/2. + tHeight/2. + 4.0*CLHEP::cm; 
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, tLVol, bName + std::string("_P"),
                       mother->GetLogicalVolume(), false, 1, true);
 

}

//Call the Function to construct the Horn2 tracking plane Amit Bashyal
void LBNEDetectorConstruction::ConstructLBNEHorn2TrackingPlane(G4VPhysicalVolume *mother) {
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   std::cerr<<"Logical Volumes are as follows: "<<nameM<<std::endl;
   if(nameM != G4String("Tunnel")){
     std::cerr
     <<"Unexpected Mother Volume in LBNEDetectorConstruction::ConstructLBNEHorn2TrackingPlane(G4VPhysicalVolume *mother)"<<std::endl;
      exit(2);
   }
   // Mother Placement not used in this context? Probably O.K.  Paul Lebrun, Fer. 2015 
//    const LBNEVolumePlacementData *plMother = 
//         fPlacementHandler->Find(G4String("Horn2TrackingPlane"), nameM, 
//                               G4String("LBNEDetectorConstruction::ConstructLBNEHorn2TrackingPlane"));
   const LBNEVolumePlacementData *plH2 = 
         fPlacementHandler->Find(G4String("Horn2TrackingPlane"), G4String("Horn2Hall"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEHorn2TrackingPlane"));
   const size_t usedNumberOfHornsPoly = static_cast<size_t>(fPlacementHandler->GetNumberOfHornsPolycone());
   //
   // Re-install Tracking 2, after HornB. It's position is Horn Geometry dependent.
   //
   double z = 1.0e20;
   if (usedNumberOfHornsPoly > 1) {
      std::cerr << " .. Using Conceptual Design 2015, fixed position of the HornBtracking place at z = 7500. mm " << std::endl;
      z = 75000*CLHEP::mm;
   } else {
      const LBNEVolumePlacementData *plH2 = 
         fPlacementHandler->Find(G4String("Horn2TrackingPlane"), G4String("Horn2Hall"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEHorn2TrackingPlane"));
        z = plH2->fPosition[2] + plH2->fParams[2]/2.0 + 2*CLHEP::cm; //2cm after the end of Horn 2
  }
   const double widthTop = 400*CLHEP::cm; //
   const double heightTop =400*CLHEP::cm;
//   G4String topNameBetween = G4String("Horn2TrackingPlane");
   G4Box *topBox = new G4Box("Horn2TrackingPlane", widthTop/2., heightTop/2.,0.5*CLHEP::mm);
   TrackingPlaneH2Logical = new G4LogicalVolume(topBox, G4Material::GetMaterial(std::string("Air")),"Horn2TrackingPlaneLogical"); 
   G4ThreeVector posTopLevel (0., 0., z);
    std::cerr
      <<"The value of Z in Horn2Plane is:" << z<<"   "<< plH2->fPosition[2]<<"   "<<plH2->fParams[2]<<std::endl;
    std::cerr
     <<"Width of Horn2TrackingPlane is:" <<widthTop/2<<" Height of Horn2TrackingPlane is: "<<heightTop/2<<std::endl;
    new G4PVPlacement((G4RotationMatrix *) 0, posTopLevel, 
                                                TrackingPlaneH2Logical,"Horn2TrackingPlane_P",
                                                mother->GetLogicalVolume(), false, 1, true);
}

void LBNEDetectorConstruction::ConstructLBNEShieldingBetweenHorns(G4VPhysicalVolume *mother) {
//
// We have to construct a mother volume located in the tunnel, in between the 2 horns 
// Based on the Horn1 and Honr 2 placements. 
//
   const double in = 25.4*CLHEP::mm;
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   if (nameM != G4String("Tunnel")) {
     std::cerr 
     << " Unexpected Mother volume in LBNEDetectorConstruction::ConstructLBNEShieldingBetweenHorns !" << std::endl;
     exit(2);
   }
   const LBNEVolumePlacementData *plMother = 
         fPlacementHandler->Find(G4String("ShieldingBetweenHorns"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingBetweenHorns"));
   const LBNEVolumePlacementData *plH1 = 
         fPlacementHandler->Find(G4String("ShieldingBetweenHornsHorn"), G4String("TargetHallAndHorn1"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingHorn1"));
   const LBNEVolumePlacementData *plH2 = 
         fPlacementHandler->Find(G4String("ShieldingBetweenHornsHorn"), G4String("Horn2Hall"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEShieldingBetweenHorns"));
   const double zStart = plH1->fPosition[2] + plH1->fParams[2]/2.0 + 5.0*CLHEP::cm;
   const double zEnd = plH2->fPosition[2] - plH2->fParams[2]/2.0 - 5.0*CLHEP::cm;
   const double lengthTop = zEnd - zStart;
//   std::cerr << " Length of the Horn1 to Horn2 corridor " << lengthTop << std::endl;
   const double widthTop = plMother->fParams[0] - 10.0*CLHEP::cm;
   const double heightTop = plMother->fParams[1] - 10.0*CLHEP::cm;
   G4String topNameBetween = G4String("Horn1ToHorn2Corridor");
   G4Box *topBox = new G4Box(topNameBetween, widthTop/2., heightTop/2.,lengthTop/2.);
   G4LogicalVolume *topLevVol = new G4LogicalVolume(topBox, G4Material::GetMaterial(G4String("Air")), topNameBetween); 
   G4ThreeVector posTopLevel (0., 0., 0.5*(zStart + zEnd));
   G4PVPlacement *vTopLevel = new G4PVPlacement((G4RotationMatrix *) 0, posTopLevel, 
                                                topLevVol , topNameBetween + std::string("_P"),
                                                mother->GetLogicalVolume(), false, 1, true);
//                     
// Now we proceed with the shielding blocks. 
//  
   const double beamlineHeight = 0.5*(60.6 + 66.7)*in; // we take the average between Horn1 and Horn2.                           
//
// Bottom 
//
   G4String bName = G4String("ShieldingHornCorridorBottom"); 
   const double bWidth = 168*in; // take the width of Horn2.  
   const double bHeight = 52.0*in;
   const double bLength = lengthTop - 2.*bHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *bBox = new G4Box(bName, bWidth/2., bHeight/2., bLength/2.);
   G4LogicalVolume *bLVol = new G4LogicalVolume(bBox, G4Material::GetMaterial(std::string("Slab_Stl")), bName); 
   G4ThreeVector posTmp(0., 0., 0.); 
   posTmp[1] = -beamlineHeight  - bHeight/2.; 
//   std::cerr << "ConstructLBNEShieldingHorn2 .... "<< bName << " zShift bottom plate " << zShift 
//           << " thus, y position " << posTmp[1] << " coorection " << yCorr << std::endl;
//           
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, bLVol, bName + std::string("_P"),
                       vTopLevel->GetLogicalVolume(), false, 1, true);
//
// Sides 
//
   G4String sName = G4String("ShieldingHornCorridorSide"); 
   const double sWidth = 52.0*in;
   const double sHeight = 256.8*in - 52.0*in - 5.0*CLHEP::cm ; // approximate  
   const double sLength = lengthTop - 2.*sHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *sBox = new G4Box(sName, sWidth/2., sHeight/2., sLength/2.);
   G4LogicalVolume *sLVol = new G4LogicalVolume(sBox, G4Material::GetMaterial(std::string("Slab_Stl")), sName); 
   posTmp[0] = (-32.0 - 26.0)*in;
   posTmp[1] = sHeight/2. - beamlineHeight + 2.0*CLHEP::cm; 
//   std::cerr << "ConstructLBNEShieldingBetweenHorns .... "<< sName  
//         <<  "   X pos. = " << posTmp[0] << " Y " << posTmp[1] << std::endl;
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_LeftP"),
                       vTopLevel->GetLogicalVolume(), false, 1, true);
   posTmp[0] = (+32.0 + 26.0)*in;
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, sLVol, sName + std::string("_RightP"),
                       vTopLevel->GetLogicalVolume(), false, 2, true);
//
// Top. Simply put about 2m of steel 
//
   G4String tName = G4String("ShieldingHornCorridorTop"); 
   const double tWidth = 168*in;
   const double tHeight = 94.0*in;
   const double tLength = lengthTop - 2.*tHeight*std::abs(std::sin(fBeamlineAngle)) - 5.0*CLHEP::cm;
   G4Box *tBox = new G4Box(tName, tWidth/2., tHeight/2., tLength/2.);
   G4LogicalVolume *tLVol = new G4LogicalVolume(tBox, G4Material::GetMaterial(std::string("Slab_Stl")), tName); 
   
   posTmp[0] =0.; 
   posTmp[1] += sHeight/2. + tHeight/2. + 4.0*CLHEP::cm; 
   new G4PVPlacement(&fRotBeamlineAngle, posTmp, tLVol, bName + std::string("_P"),
                       vTopLevel->GetLogicalVolume(), false, 1, true);
}
void LBNEDetectorConstruction::ConstructLBNFShielding(G4VPhysicalVolume *mother) {
//
// We place steel, continuous plates, after Horn1..inclined along the beam line. Unrealsitic, but simple...  
//
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   if (nameM != G4String("Tunnel")) {
     std::cerr 
     << " Unexpected Mother volume in LBNEDetectorConstruction::ConstructLBNEShieldingBetweenHorns !" << std::endl;
     exit(2);
   }
   const LBNEVolumePlacementData *plMother = 
         fPlacementHandler->Find(G4String("LBNFShielding"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNFShielding"));
   const LBNEVolumePlacementData *plH1 = 
         fPlacementHandler->Find(G4String("LBNFShielding"), G4String("TargetHallAndHorn1"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNFShielding"));
   const double zStart = plH1->fPosition[2] + plH1->fParams[2]/2.0 + 5.0*CLHEP::cm;
   const double zEnd = fPlacementHandler->GetDecayPipeLongPosition() - 5.0*CLHEP::cm;
   const double lengthAll = zEnd - zStart;
   const double widthSide = 0.5*(plMother->fParams[0] - 2.0*CLHEP::cm - fPlacementHandler->GetChaseWidthForLBNF());
   const double heightTopBot = 0.5*(plMother->fParams[1] - 2.0*CLHEP::cm - fPlacementHandler->GetChaseWidthForLBNF());
   const double widthTopBot = plMother->fParams[0] - 2.0*CLHEP::cm;
   const double heightSide = (plMother->fParams[1] - 1.0*CLHEP::cm - 2.0*heightTopBot);
   std::cerr << " Shielding parameters, width  " << widthSide << " htB " << heightTopBot 
             << " widthTopBot " << widthTopBot << " heightSide " << heightSide << std::endl;
  // Bottom, or top, since the tunnel is symmetric. But we keep names, duplicate volume, but easier to debug.  
   G4String botName = G4String("LBNFChaseShieldBottom");
   G4Box *botBox = new G4Box(botName, widthTopBot/2., heightTopBot/2., lengthAll/2.);
   G4LogicalVolume *botLevVol = new G4LogicalVolume(botBox, G4Material::GetMaterial(G4String("Slab_Stl")), botName);
   const double botY = -1.0*(plMother->fParams[1]/2. - heightTopBot/2. - 0.5*CLHEP::cm);
   G4ThreeVector botPosV (0., botY, 0.5*(zStart + zEnd));
   new G4PVPlacement((G4RotationMatrix *) 0, botPosV, botLevVol , botName + std::string("_P"),
                                                mother->GetLogicalVolume(), false, 1, true);
   G4String topName = G4String("LBNFChaseShieldTop");
   G4Box *topBox = new G4Box(topName, widthTopBot/2., heightTopBot/2., lengthAll/2.);
   G4LogicalVolume *topLevVol = new G4LogicalVolume(topBox, G4Material::GetMaterial(G4String("Slab_Stl")), topName);
   const double topY = -botY;
   G4ThreeVector topPosV (0., topY, 0.5*(zStart + zEnd));
   new G4PVPlacement((G4RotationMatrix *) 0, topPosV, topLevVol , topName + std::string("_P"),
                                                mother->GetLogicalVolume(), false, 1, true);
// 
// O.K. kind silly, for the sides,  we place here two copy of the same logical volume.                                                                                          
//
   G4String sideName = G4String("LBNFChaseShieldSide");
   G4Box *sideBox = new G4Box(sideName, widthSide/2., heightSide/2., lengthAll/2.);
   G4LogicalVolume *sideVol = new G4LogicalVolume(sideBox, G4Material::GetMaterial(G4String("Slab_Stl")), sideName);
   const double sideX = -1.0 *(plMother->fParams[0]/2. -  widthSide/2. -0.5*CLHEP::cm);
   G4ThreeVector sidePosV (-sideX, 0., 0.5*(zStart + zEnd));
   new G4PVPlacement((G4RotationMatrix *) 0, sidePosV, sideVol , sideName + std::string("_PxN"),
                                                mother->GetLogicalVolume(), false, 1, true);
   G4ThreeVector sidePosV2 (sideX, 0., 0.5*(zStart + zEnd));
   new G4PVPlacement((G4RotationMatrix *) 0, sidePosV2, sideVol , sideName + std::string("_PxP"),
                                                mother->GetLogicalVolume(), false, 2, true);                   
}
void LBNEDetectorConstruction::ConstructLBNEDecayPipeTrackingPlane(G4VPhysicalVolume *mother) 
    {
   const G4String nameM = mother->GetLogicalVolume()->GetName();
   if(nameM != G4String("Tunnel")){
     std::cerr
     <<"Unexpected Mother Volume in LBNEDetectorConstruction::ConstructLBNEDecayPipeTrackingPlane(G4VPhysicalVolume *mother)"<<std::endl;
      exit(2);
   }
         fPlacementHandler->Find(G4String("DPTrackingPlane"), nameM, 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEDecayPipeTrackingPlane"));
   const LBNEVolumePlacementData *plD1 = 
         fPlacementHandler->Find(G4String("DPTrackingPlane"), G4String("DecayPipeHall"), 
                                 G4String("LBNEDetectorConstruction::ConstructLBNEDecayPipeTrackingPlane"));
   const double z = plD1->fPosition[2] - plD1->fParams[2]/2.0 - 0.5*CLHEP::cm;

   const double widthTop =  250.00*CLHEP::cm;
   const double heightTop = 250.00*CLHEP::cm;
   
//   G4String topNameBetween = G4String("Horn1TrackingPlane");
   G4Box *topBox = new G4Box("DPTrackingPlane", widthTop/2., heightTop/2.,0.5*CLHEP::mm);
   TrackingPlaneDPLogical = new G4LogicalVolume(topBox, G4Material::GetMaterial(std::string("Air")),"DecayPipeTrackingPlaneLogical"); 
   G4ThreeVector posTopLevel (0., 0., z);
std::cerr
      <<"The value of Z in DecayPipePlane is:" << z<<"   "<< plD1->fPosition[2]<<"   "<<plD1->fParams[2]<<std::endl;
std::cerr
     <<" Half Width of DecayPipeTrackingPlane is (mm):" <<widthTop/2<<"Half Height of DecayPipeTrackingPlane is (mm): "<<heightTop/2<<std::endl;
   new G4PVPlacement((G4RotationMatrix *) 0, posTopLevel, 
                                                TrackingPlaneDPLogical,"DecayPipeTrackingPlane_P",
                                                mother->GetLogicalVolume(), false, 1, true);
}



void LBNEDetectorConstruction::DropMarsTargetHorns(G4VPhysicalVolume *tunnel) {

  // Place the MARS target and horns
   G4cout << "Importing mars target/horn gdml file... " << G4endl;
   
   G4String filename(fPlacementHandler->GetMarsTargetHornsGDMLFilename());
   std::ifstream gdmlfile(filename.c_str());
   if (!gdmlfile.is_open()) {
      std::string mess(" AbsorberGDML file "); 
      mess += filename + G4String(" could not be found \n");
      G4Exception("LBNEDetectorConstruction::DropMarsTargetHorns", " ", 
                    FatalErrorInArgument, mess.c_str());
      return; // perfunctory. 
      
    } else {
     gdmlfile.close();
    }
     G4GDMLParser parser;
     parser.Read( filename );

     G4LogicalVolume *topMars = parser.GetVolume( "Marstop" );

    
     // We dump the volume hierarchy.  Hoopefully not too deep, 
     for (int i=0; i != topMars->GetNoDaughters(); ++i) {
        G4VPhysicalVolume *pVol = topMars->GetDaughter(i);
        G4LogicalVolume *lVol = pVol->GetLogicalVolume();
        std::cerr << " Top level daughter # " << i << " Name " << lVol->GetName() 
                  << " at " << pVol->GetObjectTranslation() << std::endl;
           if (lVol->GetName().find("Airbox") != std::string::npos) {
             G4Box *aBox = static_cast<G4Box*>(lVol->GetSolid());
             std::cerr << " Airbox size " << aBox->GetXHalfLength() 
                       << " / " << aBox->GetYHalfLength() << " / " << aBox->GetZHalfLength() << std::endl; 
           }
       for (int ii=0; ii != lVol->GetNoDaughters(); ++ii) {
          G4VPhysicalVolume *pVol2 = lVol->GetDaughter(ii);
          G4LogicalVolume *lVol2 = pVol2->GetLogicalVolume();
          std::cerr << "  .. 2nd level daughter # " << ii << " Name " << lVol2->GetName()
           << " at " << pVol2->GetObjectTranslation() << std::endl;
           if (lVol2->GetName().find("Airbox") != std::string::npos) {
             G4Box *aBox = static_cast<G4Box*>(lVol2->GetSolid());
             std::cerr << " Airbox size " << aBox->GetXHalfLength() 
                       << " / " << aBox->GetYHalfLength() << " / " << aBox->GetZHalfLength() << std::endl; 
           }
         for (int iii=0; iii != lVol2->GetNoDaughters(); ++iii) {
           G4VPhysicalVolume *pVol3 = lVol2->GetDaughter(iii);
           G4LogicalVolume *lVol3 = pVol3->GetLogicalVolume();
           std::cerr << "  ... 3rd level daughter # " << iii << " Name " << lVol3->GetName() 
            << " at " << pVol3->GetObjectTranslation() << std::endl;
           for (int i4=0; i4 != lVol3->GetNoDaughters(); ++i4) {
            G4VPhysicalVolume *pVol4 = lVol3->GetDaughter(i4);
            G4LogicalVolume *lVol4 = pVol4->GetLogicalVolume();
            std::cerr << "  .... 4rth level daughter # " << i4 << " Name " << lVol4->GetName()
             << " at " << pVol4->GetObjectTranslation()  << std::endl;
            for (int i5=0; i5 != lVol4->GetNoDaughters(); ++i5) {
              G4VPhysicalVolume *pVol5 = lVol4->GetDaughter(i5);
              G4LogicalVolume *lVol5 = pVol5->GetLogicalVolume();
              std::cerr << "  ..... 5rth level daughter # " << i5 << " Name " << lVol5->GetName() 
              << " at " << pVol5->GetObjectTranslation()  << std::endl;
              for (int i6=0; i6 != lVol5->GetNoDaughters(); ++i6) {
                G4VPhysicalVolume *pVol6 = lVol5->GetDaughter(i6);
                G4LogicalVolume *lVol6 = pVol6->GetLogicalVolume();
                std::cerr << "  ...... 6rth level daughter # " << i6 << " Name " << lVol6->GetName() << std::endl;
              }
            }
          }
        }
      }
    }

     G4LogicalVolume *marsHorn1 = parser.GetVolume( "horn1" );

     const G4Box *horn1Sol = static_cast<const G4Box *>(marsHorn1->GetSolid());
//     const double marsHorn1Width = horn1Sol->GetYHalfLength();
//     const double marsHorn1Height = horn1Sol->GetXHalfLength();
//     const double marsHorn1Length = horn1Sol->GetZHalfLength();
     std::cerr << " Dimension of horn1 level MARS container, X " << horn1Sol->GetXHalfLength() << 
         " Y  "  << horn1Sol->GetYHalfLength() <<   " Z  "  << horn1Sol->GetZHalfLength() << std::endl;
     std::cerr << " Number of daughters for HORN1 " << marsHorn1->GetNoDaughters() << std::endl;
     double maxHalfHeight = -1.0;
     double maxHalfWidth = -1.0;
     double maxHalfLength = -1.0;
     for (int i=0; i != marsHorn1->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = marsHorn1->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
       std::cerr << " Daughther " << lVol->GetName();
       const G4Tubs *aTube = static_cast<const G4Tubs *>(lVol->GetSolid());

       G4ThreeVector loc = pVol->GetObjectTranslation();
       std::cerr << " at MARS coordinates " << loc[0] << ", " <<loc[1] << ", " << loc[2] << 
         " zLength " << 2.0*aTube->GetZHalfLength() << std::endl;
       // Compute the maximum height, width.  Note the confusion about  X and Y X is up, vertical, in MArs  
       if ((std::abs(loc[2]) + aTube->GetZHalfLength()) > maxHalfLength)
         maxHalfLength = std::abs(loc[2]) + aTube->GetZHalfLength();
       if ((std::abs(loc[1]) + aTube->GetOuterRadius()) >  maxHalfWidth)
         maxHalfWidth = std::abs(loc[1]) + aTube->GetOuterRadius(); // Width is along X G4lbne orientation, which Y MARS 
       if ((std::abs(loc[0]) + aTube->GetOuterRadius()) >  maxHalfHeight)
         maxHalfHeight = std::abs(loc[0]) + aTube->GetOuterRadius();
       // Height is along Y G4lbne orientation, which is negative X in MARS 
     }
     maxHalfHeight += 5.0*CLHEP::cm;
     maxHalfWidth += 5.0*CLHEP::cm;
     maxHalfLength += std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) + 5.0*CLHEP::cm;;
     std::cerr << " Container volume for MARS, 1/2 width " 
               << maxHalfWidth << " 1/2 Height " << maxHalfHeight 
             << " 1/2 length " << maxHalfLength << std::endl;

     //std::cerr << " And shorn1 after parsing the gdml file " << std::endl; exit(2);

     G4Box *aMarsBoxHorn1 = new G4Box(G4String("MarsHorn1"), maxHalfWidth, maxHalfHeight, maxHalfLength);
     G4LogicalVolume *aMarsHorn1L = 
        new G4LogicalVolume(aMarsBoxHorn1, G4Material::GetMaterial("Air"), G4String("MarsHorn1"));
     const LBNEVolumePlacementData *plTunnel = 
       fPlacementHandler->Find(G4String("MarsHorn1"), G4String("Tunnel"), 
                                 G4String("LBNEDetectorConstruction::DropMarsTargetHorns"));
     //const double zzz = maxHalfLength + plDecayPipe->fParams[2]/2 + plDecayPipe->fPosition[2] + 
     //    std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) +  5.0*CLHEP::cm;
        
     double zzz = 1.365*CLHEP::m;
        
     std::cerr << " Position in tunnel (center) " << plTunnel->fPosition[2] 
               << " ZPosMars " << zzz <<  std::endl;
       
     // One may have to adjust the height of the mars geometry ... not done yet

     //const double yyy = 2.3*CLHEP::m; 
     const double yyy = 0.; 
     G4ThreeVector posHorn1Mars(0., yyy, zzz);
     
     new G4PVPlacement(0, posHorn1Mars, aMarsHorn1L, "MarsHorn1", 
                     tunnel->GetLogicalVolume(), false, 1, true);

     G4RotationMatrix *marsRot = new G4RotationMatrix;
     marsRot->rotateZ(-M_PI/2.);
     for (int i=0; i != marsHorn1->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = marsHorn1->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
       const G4Tubs *aBox = static_cast<const G4Tubs *>(lVol->GetSolid());
       G4ThreeVector loc = pVol->GetObjectTranslation();
//       const double yyy = loc[0]  - marsHorn1Height + maxHalfHeight;  
       double yyyI = loc[0] ;  // Up to a sign!!! 
       const double xxx = loc[1]; // Y G4LBNE = -X Mars. Was centered in MARS, set to 0., o.k.
         // X G4LBNE = Y Mars. !! Not centered in Mars! Perhaps, ned a shift due to the 
        // the different size of the mother volume  
//       const double zzz = loc[2] - marsHorn1Length + maxHalfLength;
       double zzzI = loc[2]; // Setting up this last to avoid air to air volume overlap. 
       if (lVol->GetName().find("AH_horn1") != std::string::npos) {
          yyyI += 0.1*CLHEP::mm; // To avoid clash.. Cosmetic 
       // Note: for release v3r0px up to v3r0p10, this was set to 27.9 m 
          zzzI = loc[2] - 27.0*CLHEP::m; 
       }
       G4ThreeVector posTmp(xxx, yyyI, zzzI);
       std::cerr << " Placing volume " << lVol->GetName() << " at " << posTmp << " 1/2 length (G4 coord)  " 
        <<  aBox->GetZHalfLength() << std::endl;
       //std::cerr << "  .... Extend in X " << posTmp[0] - aBox->GetYHalfLength() 
       //                << " to " << posTmp[0] + aBox->GetYHalfLength() << std::endl;
       //std::cerr << "  .... Extend in Y " << posTmp[1] - aBox->GetXHalfLength() 
       //                << " to " << posTmp[1] + aBox->GetXHalfLength() << std::endl;
       //std::cerr << "  .... Extend in Z " << posTmp[2] - aBox->GetZHalfLength() 
       //                << " to " << posTmp[2] + aBox->GetZHalfLength() << std::endl;
       new G4PVPlacement(marsRot, posTmp, lVol, lVol->GetName() + std::string("_P"), aMarsHorn1L, false, 1, true);
     }

     //
     // We now clear the MARS horn1, simplify the geometry search 
     //
     //marsHorn1->ClearDaughters();

     
     // Place Horn 2
     G4LogicalVolume *marsHorn2 = parser.GetVolume( "horn2" );

     const G4Box *horn2Sol = static_cast<const G4Box *>(marsHorn2->GetSolid());
//     const double marsHorn2Width = horn2Sol->GetYHalfLength();
//     const double marsHorn2Height = horn2Sol->GetXHalfLength();
//     const double marsHorn2Length = horn2Sol->GetZHalfLength();
     std::cerr << " Dimension of horn2 level MARS container, X " << horn2Sol->GetXHalfLength() << 
         " Y  "  << horn2Sol->GetYHalfLength() <<   " Z  "  << horn2Sol->GetZHalfLength() << std::endl;
     std::cerr << " Number of daughters for HORN2 " << marsHorn2->GetNoDaughters() << std::endl;
     maxHalfHeight = -1.0;
     maxHalfWidth = -1.0;
     maxHalfLength = -1.0;
     for (int i=0; i != marsHorn2->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = marsHorn2->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
       std::cerr << " Daughther " << lVol->GetName();
       const G4Tubs *aTube = static_cast<const G4Tubs *>(lVol->GetSolid());

       G4ThreeVector loc = pVol->GetObjectTranslation();
       std::cerr << " at MARS coordinates " << loc[0] << ", " <<loc[1] << ", " << loc[2] << 
         " zLength " << 2.0*aTube->GetZHalfLength() << std::endl;
       // Compute the maximum height, width.  Note the confusion about  X and Y X is up, vertical, in MArs  
       if ((std::abs(loc[2]) + aTube->GetZHalfLength()) > maxHalfLength)
         maxHalfLength = std::abs(loc[2]) + aTube->GetZHalfLength();
       if ((std::abs(loc[1]) + aTube->GetOuterRadius()) >  maxHalfWidth)
         maxHalfWidth = std::abs(loc[1]) + aTube->GetOuterRadius(); // Width is along X G4lbne orientation, which Y MARS 
       if ((std::abs(loc[0]) + aTube->GetOuterRadius()) >  maxHalfHeight)
         maxHalfHeight = std::abs(loc[0]) + aTube->GetOuterRadius();
       // Height is along Y G4lbne orientation, which is negative X in MARS 
     }
     maxHalfHeight += 5.0*CLHEP::cm;
     maxHalfWidth += 5.0*CLHEP::cm;
     maxHalfLength += std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) + 5.0*CLHEP::cm;
     std::cerr << " Container volume for MARS, 1/2 width " 
               << maxHalfWidth << " 1/2 Height " << maxHalfHeight 
             << " 1/2 length " << maxHalfLength << std::endl;

     //std::cerr << " And shorn2 after parsing the gdml file " << std::endl; exit(2);

     G4Box *aMarsBoxHorn2 = new G4Box(G4String("MarsHorn2"), maxHalfWidth, maxHalfHeight, maxHalfLength);
     G4LogicalVolume *aMarsHorn2L = 
        new G4LogicalVolume(aMarsBoxHorn2, G4Material::GetMaterial("Air"), G4String("MarsHorn2"));
     const LBNEVolumePlacementData *plTunnel2 = 
       fPlacementHandler->Find(G4String("MarsHorn2"), G4String("Tunnel"), 
                                 G4String("LBNEDetectorConstruction::DropMarsTargetHorns"));
     //const double zzz = maxHalfLength + plDecayPipe->fParams[2]/2 + plDecayPipe->fPosition[2] + 
     //    std::abs(maxHalfHeight*std::sin(fBeamlineAngle)) +  5.0*CLHEP::cm;
        
     zzz = 8.365*CLHEP::m;
        
     std::cerr << " Position in tunnel (center) " << plTunnel2->fPosition[2] 
               << " ZPosMars " << zzz <<  std::endl;
       
     // One may have to adjust the height of the mars geometry ... not done yet

     const double yyy2 = 0; 
     G4ThreeVector posHorn2Mars(0., yyy2, zzz);
     
     new G4PVPlacement(0, posHorn2Mars, aMarsHorn2L, "MarsHorn2", 
                     tunnel->GetLogicalVolume(), false, 1, true);

     G4RotationMatrix *marsRot2 = new G4RotationMatrix;
     marsRot2->rotateZ(-M_PI/2.);
     for (int i=0; i != marsHorn2->GetNoDaughters(); ++i) {
       G4VPhysicalVolume *pVol = marsHorn2->GetDaughter(i);
       G4LogicalVolume *lVol = pVol->GetLogicalVolume();
       const G4Tubs *aBox = static_cast<const G4Tubs *>(lVol->GetSolid());
       G4ThreeVector loc = pVol->GetObjectTranslation();
//       const double yyy = loc[0]  - marsHorn2Height + maxHalfHeight;  
       double yyy3 = loc[0] ;  // Up to a sign!!! 
       const double xxx = loc[1]; // Y G4LBNE = -X Mars. Was centered in MARS, set to 0., o.k.
         // X G4LBNE = Y Mars. !! Not centered in Mars! Perhaps, ned a shift due to the 
        // the different size of the mother volume  
//       const double zzz = loc[2] - marsHorn2Length + maxHalfLength;
       double zzz3 = loc[2]; // Setting up this last to avoid air to air volume overlap. 
       if (lVol->GetName().find("AH_horn2") != std::string::npos) {
          yyy3 += 0.1*CLHEP::mm; // To avoid clash.. Cosmetic 
       // Note: for release v3r0px up to v3r0p10, this was set to 27.9 m 
          zzz3 = loc[2] - 27.0*CLHEP::m; 
       }
       G4ThreeVector posTmp(xxx, yyy3, zzz3);
       std::cerr << " Placing volume " << lVol->GetName() << " at " << posTmp << " 1/2 length (G4 coord)  " 
        <<  aBox->GetZHalfLength() << std::endl;
       //std::cerr << "  .... Extend in X " << posTmp[0] - aBox->GetYHalfLength() 
       //                << " to " << posTmp[0] + aBox->GetYHalfLength() << std::endl;
       //std::cerr << "  .... Extend in Y " << posTmp[1] - aBox->GetXHalfLength() 
       //                << " to " << posTmp[1] + aBox->GetXHalfLength() << std::endl;
       //std::cerr << "  .... Extend in Z " << posTmp[2] - aBox->GetZHalfLength() 
       //                << " to " << posTmp[2] + aBox->GetZHalfLength() << std::endl;
       new G4PVPlacement(marsRot2, posTmp, lVol, lVol->GetName() + std::string("_P"), aMarsHorn2L, false, 1, true);
     }

     //
     // We now clear the MARS horn2, simplify the geometry search 
     //
     //marsHorn2->ClearDaughters();

    //
    // Place the decay Pipe Snout, which contains the window, in case we have Helium gas. 
    //
    fPlacementHandler->Create(G4String("DecayPipeSnout")); // Now in Snout region 
    fPlacementHandler->PlaceFinalDecayPipeSnout((G4PVPlacement*) tunnel);

    //   
    // Place the decay pipe 
    //   
    fPlacementHandler->Create(G4String("DecayPipeHall"));
    G4PVPlacement *vDecayPipe = fPlacementHandler->PlaceFinal(G4String("DecayPipeHall"), tunnel);
    fPlacementHandler->Create(G4String("DecayPipeConcrete"));
    fPlacementHandler->Create(G4String("DecayPipeOuterWall"));
    fPlacementHandler->Create(G4String("DecayPipeWall"));
    fPlacementHandler->Create(G4String("DecayPipeVolume"));
    //   fPlacementHandler->Create(G4String("DecayPipeUpstrWindow")); // Now in Snout region 
    
    fPlacementHandler->PlaceFinal(G4String("DecayPipeConcrete"), vDecayPipe);
    fPlacementHandler->PlaceFinal(G4String("DecayPipeOuterWall"), vDecayPipe);
    fPlacementHandler->PlaceFinal(G4String("DecayPipeWall"), vDecayPipe);
    fPlacementHandler->PlaceFinal(G4String("DecayPipeVolume"), vDecayPipe);
    
   //!!: Define geometry absorber from /LBNE/det/ConstructSimpAbsorber true/false option
    if (this->fConstructSimpAbsorber){
     this->ConstructLBNEHadronAbsorberSimple(tunnel);
    }
    else{
     this->ConstructLBNEHadronAbsorber(tunnel);
    }

}


LBNEDetectorMessenger::LBNEDetectorMessenger( LBNEDetectorConstruction* LBNEDet):LBNEDetector(LBNEDet)
{
   LBNEDir = new G4UIdirectory("/LBNE/");
   LBNEDir->SetGuidance("UI commands for detector geometry");
   
   detDir = new G4UIdirectory("/LBNE/det/");
   detDir->SetGuidance("detector control");
   
   
   ConstructTarget = new G4UIcmdWithABool("/LBNE/det/constructTarget",this); 
   ConstructTarget->SetParameterName("constructTarget", false);
   ConstructTarget->SetGuidance("Target construction on/off"); 
   ConstructTarget->SetParameterName("constructTarget",true); 
   ConstructTarget->AvailableForStates(G4State_PreInit,G4State_Idle);
   
   SetBeamlineAngle = new
   G4UIcmdWithADoubleAndUnit("/LBNE/det/setBeamlineAngle",this);
   SetBeamlineAngle->SetParameterName("angle", false);
   SetBeamlineAngle->SetGuidance("Set the angle of the beamline");
   
   UpdateCmd = new G4UIcmdWithoutParameter("/LBNE/det/update",this);
   UpdateCmd->SetGuidance("Update or Construct LBNE geometry. Same difference ");
   UpdateCmd->SetGuidance("This command MUST be applied before \"beamOn\" ");
   UpdateCmd->SetGuidance("if you changed geometrical value(s).");
   UpdateCmd->AvailableForStates(G4State_PreInit);
   
   ConstructCmd = new G4UIcmdWithoutParameter("/LBNE/det/construct",this);
   ConstructCmd->SetGuidance("Construct LBNE geometry. Should be one and only time ");
   ConstructCmd->SetGuidance("This command MUST be applied before \"beamOn\" ");
   ConstructCmd->SetGuidance("if you changed geometrical value(s).");
   ConstructCmd->AvailableForStates(G4State_PreInit);
  
    //Messenger for simple geometry of absorber
    ConstructSimpAbsorber = new G4UIcmdWithABool("/LBNE/det/ConstructSimpAbsorber",this);
    ConstructSimpAbsorber->SetParameterName("Simple Absorber",true);
    ConstructSimpAbsorber->SetGuidance("Set geometry of Hadron Absorber. If true simulate simple version of absorber");
    ConstructSimpAbsorber->SetDefaultValue(false);
    ConstructSimpAbsorber->AvailableForStates(G4State_PreInit,G4State_Idle);

    ConstructSculptedAbsorber = new G4UIcmdWithABool("/LBNE/det/ConstructSculptedAbsorber",this);
    ConstructSculptedAbsorber->SetParameterName("Aluminum, Sculpted Absorber",true);
    ConstructSculptedAbsorber->SetGuidance("Set geometry of Hadron Absorber. If true simulate the 2015 Sculpted absorber");
    ConstructSculptedAbsorber->SetDefaultValue(false);
    ConstructSculptedAbsorber->AvailableForStates(G4State_PreInit,G4State_Idle);

    DisableSpoiler = new G4UIcmdWithABool("/LBNE/det/DisableSpoiler",this);
    DisableSpoiler->SetParameterName("Remove spoiler",true);
    DisableSpoiler->SetGuidance("If true simulate the 2015 Sculpted absorber without the spoiler");
    DisableSpoiler->SetDefaultValue(false);
    DisableSpoiler->AvailableForStates(G4State_PreInit,G4State_Idle);

    DisableSculptedLayers = new G4UIcmdWithABool("/LBNE/det/DisableSculptedLayers",this);
    DisableSculptedLayers->SetParameterName("Remove sculpting",true);
    DisableSculptedLayers->SetGuidance("If true simulate the 2015 Sculpted absorber without the sculpting on the aluminum layers");
    DisableSculptedLayers->SetDefaultValue(false);
    DisableSculptedLayers->AvailableForStates(G4State_PreInit,G4State_Idle);

    DisableMask = new G4UIcmdWithABool("/LBNE/det/DisableMask",this);
    DisableMask->SetParameterName("Remove sculpting",true);
    DisableMask->SetGuidance("If true simulate the 2015 Sculpted absorber without the mask layers between the spoiler and core part of the absorber");
    DisableMask->SetDefaultValue(false);
    DisableMask->AvailableForStates(G4State_PreInit,G4State_Idle);

    ExpandAlLayers = new G4UIcmdWithABool("/LBNE/det/ExpandAlLayers",this);
    ExpandAlLayers->SetParameterName("Expand Sculpted Absorber Al layers to 2.8x2.8 m",true);
    ExpandAlLayers->SetGuidance("If true simulate the 2015 Sculpted absorber with a core region to the size of the decay pipe");
    ExpandAlLayers->SetDefaultValue(false);
    ExpandAlLayers->AvailableForStates(G4State_PreInit,G4State_Idle);

    RemoveLayers = new G4UIcmdWithAnInteger("/LBNE/det/RemoveLayers",this);
    RemoveLayers->SetParameterName("Number of layers to remove from the downstream end of the absorber",true);
    RemoveLayers->SetGuidance("Simulate the 2015 Sculpted absorber with fewer layers");
    RemoveLayers->SetDefaultValue(0);
    RemoveLayers->AvailableForStates(G4State_PreInit,G4State_Idle);

    DwStrAbsSteelWidth = new G4UIcmdWithADoubleAndUnit("/LBNE/det/DwStrAbsSteelWidth",this);
    DwStrAbsSteelWidth->SetParameterName("Steel Downstr Width",true);
    DwStrAbsSteelWidth->SetGuidance("Define the downstream width (m) of the steel volume in the simple geometry of absorber. Default value is ...");
    LBNEDetector->SetDwStrAbsSteelWidth(1.865*CLHEP::m);
    DwStrAbsSteelWidth->SetDefaultUnit("m");
    DwStrAbsSteelWidth->AvailableForStates(G4State_PreInit);

    DwStrAbsConcreteWidth = new G4UIcmdWithADoubleAndUnit("/LBNE/det/DwStrAbsConcreteWidth",this);
    DwStrAbsConcreteWidth->SetParameterName("Concrete Downstr Width",true);
    DwStrAbsConcreteWidth->SetGuidance("Define the downstream width (m) of the concrete volume in the simple geometry of absorber. Default value is ...");
    LBNEDetector->SetDwStrAbsConcreteWidth(1.33*CLHEP::m);
    DwStrAbsConcreteWidth ->SetDefaultUnit("m");
    DwStrAbsConcreteWidth->AvailableForStates(G4State_PreInit);

   new G4UnitDefinition("kiloampere" , "kA", "Electric current", 1000.*CLHEP::ampere);   
   SetHornCurrent = new
     G4UIcmdWithADoubleAndUnit("/LBNE/det/seHornCurrent",this);
   SetHornCurrent->SetParameterName("Horn Current", false);
   SetHornCurrent->SetGuidance(" The current for the horn system. (Horn1 and Horn2 are in series.");
   SetHornCurrent ->SetDefaultValue(200.0); // CDR, March 2012
   SetHornCurrent->SetDefaultUnit("kA");
   SetHornCurrent->SetUnitCandidates("kA");
   SetHornCurrent->AvailableForStates(G4State_PreInit);
   {
   SetSkinDepthInnerRad  = new G4UIcmdWithADoubleAndUnit("/LBNE/det/SetSkinDepthInnerRad",this);
   std::string aGuidance("The skin depth for the inner conductor. \n ");
   aGuidance += std::string(" The default assumes infinite.  Using Zarko' thesis model, \n"); 
   aGuidance += std::string("  MINOS Document 5694-v1, section 5.2.4, p150 - 156 \n");
   SetSkinDepthInnerRad->SetGuidance(aGuidance.c_str());
   SetSkinDepthInnerRad->SetParameterName("SetSkinDepthInnerRad",true);
   SetSkinDepthInnerRad->SetDefaultValue (1.0e10*CLHEP::mm);
   SetSkinDepthInnerRad->AvailableForStates(G4State_PreInit);
   }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "DeltaEccentricityIO"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADoubleAndUnit* aGUICmd  = new G4UIcmdWithADoubleAndUnit(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The Eccentricity between Inner/Outer conductors for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxDeltaEccentricityIO",true);
       aGUICmd->SetDefaultValue (0.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsDeltaEccentricityIO.push_back(aGUICmd);
     }
    }
    {
     //
     
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "DeltaEllipticityI"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADoubleAndUnit* aGUICmd  = new G4UIcmdWithADoubleAndUnit(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The Ellipticity of the  Inner conductor for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxDeltaEllipticityI",true);
       aGUICmd->SetDefaultValue (0.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsDeltaEllipticityI.push_back(aGUICmd);
     }
   }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; 
       aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "CurrentMultipiler"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADouble* aGUICmd  = new G4UIcmdWithADouble(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "A magic current amplifier to increase/decrease the current for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxCurrentMultiplier",true);
       aGUICmd->SetDefaultValue (1.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsCurrentMultiplier.push_back(aGUICmd);
     }
    }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; 
       aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "CurrentEqualizerLongAbsLength"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADoubleAndUnit* aGUICmd  = new G4UIcmdWithADoubleAndUnit(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The charact. length of the decrease of the azimuthal anisotropy \n " 
         << " due imperfect current equal. section  for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxCurrentEqualizerLongAbsLength",true);
       aGUICmd->SetDefaultValue (100.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsCurrentEqualizerLongAbsLength.push_back(aGUICmd);
     }
    }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; 
       aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "CurrentEqualizerQuadAmpl"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADouble* aGUICmd  = new G4UIcmdWithADouble(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The amplitude of the azimuthal, quadrupolar, anisotropy \n "  
         << " due imperfect current equal. section  for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxCurrentEqualizerQuadAmpl",true);
       aGUICmd->SetDefaultValue (0.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsCurrentEqualizerQuadAmpl.push_back(aGUICmd);
     }
    }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; 
       aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "CurrentEqualizerOctAmpl"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithADouble* aGUICmd  = new G4UIcmdWithADouble(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The amplitude of the azimuthal, quadrupolar, anisotropy \n "  
         << " due imperfect current equal. section  for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxCurrentEqualizerOctAmpl",true);
       aGUICmd->SetDefaultValue (0.);
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetHornsCurrentEqualizerOctAmpl.push_back(aGUICmd);
     }
    }
   {
     for (size_t kH=0; kH != 3; kH++) { // We will work with only 3 horns. Per executive decision of Laura F.
       std::ostringstream aNStrStr; 
       aNStrStr << "/LBNE/det/SetHorn" << (kH+1) << "CurrentEqualizerMapFileName"; 
       std::string aNStr(aNStrStr.str());
       G4UIcmdWithAString* aGUICmd  = new G4UIcmdWithAString(aNStr.c_str(),this);
       std::ostringstream aNGStrStr; aNGStrStr 
         << "The file name for the field map created from a priot run, \n "  
         << " related to imperfect current equal. section  for Horn number " << (kH+1) << "\n."; 
       std::string aNGStr(aNGStrStr.str()); 
       aGUICmd->SetGuidance(aNGStr.c_str());
       aGUICmd->SetParameterName("SetHornxCurrentEqualizerFileName",true);
       aGUICmd->SetDefaultValue("?");
       aGUICmd->AvailableForStates(G4State_PreInit);
       SetFileNameFieldMapForCE.push_back(aGUICmd);
     }
    }
   {
   ZCoordForPerfectFocusing  = new G4UIcmdWithADoubleAndUnit("/LBNE/det/ZCoordForPerfectFocusing",this);
   std::string aGuidance("The Z-location at which the perfect focusing process set the track polar angle to zero \n ");
   aGuidance += std::string(" If this Z  is > 0. , we trigger the instantaneous foucing at this Z  \n");
   aGuidance += std::string(" else, Z < 0., we trigger the instantaneous foucing at the radius of Target Canister  \n");
   aGuidance += std::string(" which, for the LBNE 1.2 MW option is 17.6 mm  \n");
   ZCoordForPerfectFocusing->SetGuidance(aGuidance.c_str());
   ZCoordForPerfectFocusing->SetParameterName("ZCoordForPerfectFocusing", true);
   ZCoordForPerfectFocusing->SetDefaultUnit ("m");
   ZCoordForPerfectFocusing->SetUnitCandidates ("mm cm m");
   ZCoordForPerfectFocusing->SetDefaultValue (1.0e103);
   ZCoordForPerfectFocusing->AvailableForStates(G4State_PreInit);
   }
   /*   
   {
   ConstructLBNFNano = new G4UIcmdWithABool("/LBNF/det/ConstructLBNFNano",this); 
   ConstructLBNFNano->SetParameterName("ConstructLBNFNano", false);
   ConstructLBNFNano->SetGuidance("Neutrino Progenitor Spectrometer construction on/off"); 
   ConstructLBNFNano->SetParameterName("ConstructLBNFNano",true); 
   ConstructLBNFNano->AvailableForStates(G4State_PreInit,G4State_Idle);
   }
   */
}

LBNEDetectorMessenger::~LBNEDetectorMessenger() 
{

   delete detDir;
   delete LBNEDir;
   delete ConstructSimpAbsorber;
   delete ConstructSculptedAbsorber;
   delete DwStrAbsConcreteWidth;
   delete DwStrAbsSteelWidth;
   delete ConstructTarget;
   delete SetBeamlineAngle;
   delete UpdateCmd;
   delete ConstructCmd;
   delete SetHornCurrent;
   delete SetSkinDepthInnerRad;
   for (size_t kH=0; kH != 3; kH++) {
     delete SetHornsDeltaEccentricityIO[kH];
     delete SetHornsDeltaEllipticityI[kH];
     delete SetHornsCurrentEqualizerLongAbsLength[kH];
     delete SetHornsCurrentEqualizerQuadAmpl[kH];
     delete SetFileNameFieldMapForCE[kH];
     delete SetHornsCurrentEqualizerOctAmpl[kH];
     delete SetHornsCurrentMultiplier[kH];
   }
   delete ZCoordForPerfectFocusing;
   delete DisableSpoiler;
   delete DisableSculptedLayers;
   delete DisableMask;
   delete ExpandAlLayers;
   delete RemoveLayers;
}


void LBNEDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newValue)
{

   //Turn on and off simple geometry
   //Note: ConstructSimpAbsorber defined in LBNEDetectorMessemger.hh
   //      SetConstructSimpAbsorber and fConstructSimpAbsorber defined in LBNEDetectorConstruction.hh
   if ( command == ConstructSimpAbsorber ) 
     {
       LBNEDetector->SetConstructSimpAbsorber(ConstructSimpAbsorber->GetNewBoolValue(newValue));
     }
   
   if ( command == ConstructSculptedAbsorber ) 
     {
       LBNEDetector->SetConstructSculptedAbsorber(ConstructSculptedAbsorber->GetNewBoolValue(newValue));
     }

   if ( command == DisableSpoiler ) 
     {
       LBNEDetector->SetDisableSpoiler(DisableSpoiler->GetNewBoolValue(newValue));
     }
  
   if ( command == DisableSculptedLayers ) 
     {
       LBNEDetector->SetDisableSculptedLayers(DisableSculptedLayers->GetNewBoolValue(newValue));
     }

   if ( command == DisableMask )
     {
       LBNEDetector->SetDisableMask(DisableMask->GetNewBoolValue(newValue));  
     }   

   if ( command == ExpandAlLayers )
     {
       LBNEDetector->SetExpandAlLayers(ExpandAlLayers->GetNewBoolValue(newValue));
     }

   if ( command == RemoveLayers )
     {
       LBNEDetector->SetRemoveLayers(RemoveLayers->GetNewIntValue(newValue));
     }

   if ( command == DwStrAbsSteelWidth )
     {
       G4UIcmdWithADoubleAndUnit* cmdSteelDwStr = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
       LBNEDetector->SetDwStrAbsSteelWidth(cmdSteelDwStr->GetNewDoubleValue(newValue));  
       //G4double SteelWidth = LBNEDetector->GetDwStrAbsSteelWidth();
     }
 

    if ( command == DwStrAbsConcreteWidth)
     {
       G4UIcmdWithADoubleAndUnit* cmdConcreteDwStr = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
       LBNEDetector->SetDwStrAbsConcreteWidth(cmdConcreteDwStr->GetNewDoubleValue(newValue));
       // G4double ConcWidth = LBNEDetector->GetDwStrAbsConcreteWidth();
     }

  if ( command == UpdateCmd ) 
   {
      LBNEDetector->Construct();
      return;
   }
  if ( command == ConstructCmd ) 
   {
      LBNEDetector->Construct();
      return;
   }
   if (command == SetHornCurrent) {
     G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
     LBNEDetector->SetHornCurrent(cmdWD->GetNewDoubleValue(newValue));
   }
   if (command == SetSkinDepthInnerRad) {
     G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
     LBNEDetector->SetSkinDepthInnerRad(cmdWD->GetNewDoubleValue(newValue));
   }
   for (size_t kH=0; kH != 3; kH++) {
     if (command == SetHornsDeltaEccentricityIO[kH]) {
       G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
       LBNEDetector->SetDeltaEccentricityIO(kH, cmdWD->GetNewDoubleValue(newValue));
     }
     if (command == SetHornsDeltaEllipticityI[kH]) {
       G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
       LBNEDetector->SetDeltaEllipticityI(kH, cmdWD->GetNewDoubleValue(newValue));
     }
     if (command == SetHornsCurrentEqualizerLongAbsLength[kH]) {
       G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
       LBNEDetector->SetCurrentEqualizerLongAbsLength(kH, cmdWD->GetNewDoubleValue(newValue));
     }
     if (command == SetHornsCurrentEqualizerQuadAmpl[kH]) {
       G4UIcmdWithADouble* cmdWD = dynamic_cast<G4UIcmdWithADouble*> (command);
       LBNEDetector->SetCurrentEqualizerQuadAmpl(kH, cmdWD->GetNewDoubleValue(newValue));
     }
     if (command == SetFileNameFieldMapForCE[kH]) {
       LBNEDetector->SetFileNameFieldMapForCE(kH, newValue);
     }
     if (command == SetHornsCurrentEqualizerOctAmpl[kH]) {
       G4UIcmdWithADouble* cmdWD = dynamic_cast<G4UIcmdWithADouble*> (command);
       LBNEDetector->SetCurrentEqualizerOctAmpl(kH, cmdWD->GetNewDoubleValue(newValue));
     }
     if (command == SetHornsCurrentMultiplier[kH]) {
       G4UIcmdWithADouble* cmdWD = dynamic_cast<G4UIcmdWithADouble*> (command);
       LBNEDetector->SetCurrentMultipiler(kH, cmdWD->GetNewDoubleValue(newValue));
     }
   }
   if (command == ZCoordForPerfectFocusing) {
     G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (command);
     double val = cmdWD->GetNewDoubleValue(newValue);
     LBNEDetector->SetZCoordForPerfectFocusing(val);
     if (val < 0.) {
       std::cout << " You have chosen a near perfect focus exiting the target Helium containment vessel " << std::endl;
       std::cout << " This only make sens if the real horn current is vanishlingly small.. " << std::endl;
       LBNEDetector->SetHornCurrent(2.0e-11);
       std::cout << " ...Done.. Also, Horn's conductors are made of air " << std::endl;
       LBNEVolumePlacements* volP=LBNEVolumePlacements::Instance();
       G4String matAir("Air");
       volP->SetHorn1InnerConductorMaterial(matAir);
       volP->SetHorn1AllConductorMaterial(matAir);
       volP->SetHorn2InnerConductorMaterial(matAir);
       volP->SetHorn2AllConductorMaterial(matAir);
     } else if (val > 22499.*CLHEP::mm) {
       std::cout << " You have chosen a near perfect focus entering the decay channel at Z =  " << val << std::endl;
     }
   }
   /*
   ** not yet deployed.. 
   **
   if (command == ConstructLBNFNano) {
      G4UIcmdWithABool* cmdWD = dynamic_cast<G4UIcmdWithABool*> (command);
       LBNEDetector->SetConstructLBNFNano(cmdWD->GetNewBoolValue(newValue));
   }
   */
}