Feature/aran geometry

python/duneggd/larfd/CalculateVolume.C

@@ -0,0 +1,115 @@
+void CalculateVolume(std::string fname="larfd_rn30cm_noOpDet_v3.gdml",
+                    std::string topvol="volWorld",
+                    // std::string vol="volColdCryoLayer1",
+                    std::string vol="volWorld",
+                    // std::string vol="volWarmSkin",
+		     int npoint=pow(10, 9)) {
+
+  gSystem->Load("libGeom");
+  gSystem->Load("libGdml");
+
+  TGeoManager::Import(fname.c_str());
+  TGeoVolume* v = gGeoManager->GetVolume(vol.c_str());
+  TGeoVolume* tv = gGeoManager->GetVolume(topvol.c_str());
+  if (v==nullptr) {
+    std::cout << "Didn't find the volume " + vol << "\n";
+    return;
+  }
+  if (tv==nullptr) {
+    std::cout << "Didn't find the top volume " + topvol << "\n";
+    return;
+  }
+  tv->SetAsTopVolume();
+  TGeoShape* s = v->GetShape();
+  TGeoBBox *box = (TGeoBBox *)s;
+  double dx = box->GetDX();
+  double dy = box->GetDY();
+  double dz = box->GetDZ();
+
+  double worldVolume = (2*dx) * (2*dy) * (2*dz);
+  int    exponent    = (int)log10(fabs(worldVolume));
+  double mantissa    = worldVolume / pow(10, exponent);
+  printf("Volume of the world box: ");
+  printf("%f x 10^%d ", mantissa, exponent-6);
+  printf("[m^3]\n");
+
+  double ox = (box->GetOrigin())[0];
+  double oy = (box->GetOrigin())[1];
+  double oz = (box->GetOrigin())[2];
+
+  double *xyz = new double[3];
+  TGeoNode *node = 0;
+  int i=0;
+  std::map<TGeoMaterial*, int> materials;
+  std::map<TGeoVolume*,   int> volumes  ;
+  std::map<TGeoMedium*,   int> mediums  ;
+
+  while (i<npoint) {
+    if (i%10000==0) {
+      std::cout << "\r" << 100.*i/npoint << "\t%" << std::flush;
+    }
+    xyz[0] = ox-dx+2*dx*gRandom->Uniform();
+    xyz[1] = oy-dy+2*dy*gRandom->Uniform();
+    xyz[2] = oz-dz+2*dz*gRandom->Uniform();
+    gGeoManager->SetCurrentPoint(xyz);
+    node = gGeoManager->FindNode();
+    if (!node) continue;
+    if (node->IsOverlapping()) continue;
+    i++;
+
+    TGeoMedium  * med = node->GetMedium();
+    TGeoMaterial* mat = med ->GetMaterial();
+    TGeoVolume  * vol = node->GetVolume();
+    materials[mat]++;
+    mediums  [med]++;
+    volumes  [vol]++;
+  }
+
+  std::cout << "\n\n\n";
+  std::cout << npoint << " points thrown in " << vol << ".\n";
+  std::cout << "\n\n\n";
+  printf(" %-20s", "Medium");
+  printf(" %-20s", "Volume [%]");
+  printf(" %-30s", "Volume [cc]");
+  printf(" %-20s", "Error [%]");
+  printf("\n");
+  for (auto const& it: mediums) {
+    double number  = 100.*it.second / npoint;
+    double physVol = number * worldVolume;
+    double error   = number * sqrt(1. / it.second + 1. / npoint);
+    printf(" %-20s", it.first->GetName());
+    printf(" %-20f", number);
+    printf(" %-30f", physVol);
+    printf(" %-20f\n", error);
+  }
+
+  std::cout << "\n\n\n";
+  printf(" %-20s", "Material");
+  printf(" %-20s", "Volume [%]");
+  printf(" %-30s", "Volume [cc]");  
+  printf(" %-20s", "Error [%]");
+  printf("\n");
+  for (auto const& it: materials) {
+    double number  = 100.*it.second / npoint;
+    double physVol = number * worldVolume;
+    double error   = number * sqrt(1. / it.second + 1. / npoint);
+    printf(" %-20s", it.first->GetName());
+    printf(" %-20f", number); 
+    printf(" %-30f", physVol);   
+    printf(" %-20f\n", error);
+  }
+
+  std::cout << "\n\n\n";
+  // printf(" %-20s", "Volume");
+  // printf(" %-20s", "Volume [%]");
+  // printf(" %-20s", "Error [%]");
+  // printf("\n");
+  // for (auto const& it: volumes) {
+  //   double number = 100.*it.second / npoint;
+  //   double error  = number * sqrt(1. / it.second + 1. / npoint);
+  //   printf(" %-20s", it.first->GetName());
+  //   printf(" %-20.6f", number);
+  //   printf(" %-20.6f\n", error);
+  // }
+
+}

python/duneggd/larfd/CheckOverlaps.C

@@ -0,0 +1,41 @@
+
+
+#ifndef __CINT__
+// PUT HEADERS HERE
+#endif
+
+void CheckOverlaps() {
+
+  std::string GeoWithWires    = "larfd_rn200cm_noOpDet_Wood.gdml";
+  std::string GeoWithoutWires = "larfd_rn200cm_noOpDet_Wood_nowires.gdml";
+
+
+  gSystem->Load("libGeom");
+  gSystem->Load("libGdml");
+
+  TGeoManager::Import(GeoWithWires.c_str());
+
+  std::string topVol ="volWorld";
+
+  // gGeoManager->GetTopNode();
+  gGeoManager->CheckOverlaps(1e-5,"d");
+  gGeoManager->PrintOverlaps();
+
+  TFile *tf1 = new TFile("overlaps.root", "RECREATE");
+  gGeoManager->Write();
+  tf1->Close();
+
+  TGeoManager::Import(GeoWithoutWires.c_str());
+
+  // gGeoManager->GetTopNode();
+  gGeoManager->CheckOverlaps(1e-5,"d");
+  gGeoManager->PrintOverlaps();
+
+  TFile *tf2 = new TFile("overlaps_nowires.root", "RECREATE");
+  gGeoManager->Write();
+  tf2->Close();
+
+
+
+}
+

python/duneggd/larfd/Configurations/MakeConfigs.sh

@@ -0,0 +1,5 @@
+cp larfd_nowires_AllAir_20cm.cfg larfd_nowires_AllAir_40cm.cfg
+cp larfd_nowires_AllAir_20cm.cfg larfd_nowires_AllAir_60cm.cfg
+cp larfd_nowires_AllAir_20cm.cfg larfd_nowires_AllAir_80cm.cfg
+cp larfd_nowires_AllAir_20cm.cfg larfd_nowires_AllAir_100cm.cfg
+cp larfd_nowires_AllAir_20cm.cfg larfd_nowires_AllAir_200cm.cfg

python/duneggd/larfd/Configurations/larfd_rn200cm_noOpDet.cfg

@@ -0,0 +1,218 @@
+[TPCPlaneZ]
+class                   = TPCPlane.TPCPlaneBuilder
+view                    = 'Z'
+wirePitch               = Q('0.479cm')
+wireAngle               = Q('0deg')
+nChannels               = 960
+nowires                 = False
+wireDiam                = Q('152um')
+# wireDiam                = Q('1520um')
+G10ThicknessFoot        = Q('1mm')
+G10ThicknessSide        = None
+# the position of the wire wrt the wire coming out of the boards
+# Z direction is 2.20mm(bottom left) - 44.8mm (see 8760104)
+HeadBoardScrewCentre    = [Q('0m'), Q('45mm'), Q('42.6mm')]
+HeadAPAFrameScrewCentre = [Q('0m'), Q('55.73mm'), Q('45mm')]
+wrapCover               =  Q('0.0625in')
+APAFrameDim             = {APAFrame:size}
+# planeDim                = [{TPCPlaneZ:wireDiam}, Q('606cm')-Q('7.61cm'), Q('229.4562cm')]
+planeDim                = [{TPCPlaneZ:wireDiam}, Q('606cm')-Q('7.61cm'), Q('229.593cm')]
+
+[TPCPlaneV]
+class                   = TPCPlane.TPCPlaneBuilder
+view                    = 'V'
+wirePitch               = Q('0.4669cm')
+wireAngle               = Q('-35.7deg')
+nChannels               = 800
+nowires                 = {TPCPlaneZ:nowires}
+wireDiam                = {TPCPlaneZ:wireDiam} 
+G10ThicknessFoot        = Q('1mm')
+G10ThicknessSide        = Q('3.18mm')
+# the position of the wire wrt the wire coming out of the boards
+# Z direction is ( 45mm - 30.96mm ) * tan (35.7deg) - 44.8mm (see 8760108)
+# HeadBoardScrewCentre    = [Q('0m'), Q('45mm'), Q('40.71mm')]
+HeadBoardScrewCentre    = [Q('0m'), Q('45mm'), Q('34.71mm')]
+HeadAPAFrameScrewCentre = {TPCPlaneZ:HeadAPAFrameScrewCentre}
+SideWrappingBoardOffset = Q('12.98mm') - Q('11.08mm')
+wrapCover               = {TPCPlaneZ:wrapCover}
+APAFrameDim             = {APAFrame:size}
+# planeDim                = [{TPCPlaneV:wireDiam}, Q('606cm')-Q('7.61cm')+Q('0.125in'), Q('230.0175cm')]
+planeDim                = [{TPCPlaneV:wireDiam}, Q('606cm')-Q('7.61cm')+Q('0.125in'), Q('230.0025cm')]
+
+[TPCPlaneU]
+class                   = TPCPlane.TPCPlaneBuilder
+view                    = 'U'
+wirePitch               = Q('0.4669cm')
+wireAngle               = Q('35.7deg')
+nChannels               = 800
+nowires                 = {TPCPlaneZ:nowires}
+wireDiam                = {TPCPlaneZ:wireDiam} 
+G10ThicknessFoot        = Q('1mm')
+G10ThicknessSide        = Q('3.18mm')
+HeadBoardScrewCentre    = {TPCPlaneZ:HeadBoardScrewCentre}
+HeadAPAFrameScrewCentre = {TPCPlaneZ:HeadAPAFrameScrewCentre}
+SideWrappingBoardOffset = Q('12.98mm') - Q('11.08mm')
+wrapCover               = {TPCPlaneZ:wrapCover}
+APAFrameDim             = {APAFrame:size}
+# planeDim                = [{TPCPlaneU:wireDiam}, Q('606cm')-Q('7.61cm')+Q('0.25in'), Q('230.6525cm')]
+planeDim                = [{TPCPlaneU:wireDiam}, Q('606cm')-Q('7.61cm')+Q('0.25in'), Q('230.6375cm')]
+
+[APAFrame]
+# https://edms.cern.ch/ui/file/2112698/1/8760050_ASSEMBLY_REV_WIP.PDF and ref therein.
+# Y,Z (there is no x extent, its the size of the foot beam)
+class                      = APAFrame.APAFrameBuilder
+size                       = [Q('4in'),Q('6060mm'), Q('2315.9mm')]
+footsize                   = [{APAFrame:size}[0], Q('4in'), {APAFrame:size}[2]]
+headsize                   = [{APAFrame:size}[0], Q('6in'), {APAFrame:size}[2]]
+footthickness              = Q('3.05mm')
+nribs                      = 4
+ribsize                    = [Q('5.08cm'), Q('4in')]
+ribthickness               = Q('3.05mm')
+APAGap_y                   = {TPC:APAGap_y}
+APAGap_z                   = {TPC:APAGap_z}
+
+APAFrameZSide_y            = Q('4in')
+APAFrameZSide_z            = Q('4in')
+
+ArapucaIn_x                = Q("1.10cm")
+ArapucaIn_y                = Q("4.65cm")
+ArapucaIn_z                = Q("23.4cm")
+
+ArapucaAcceptanceWindow_x  = Q("2.2cm")
+ArapucaAcceptanceWindow_y  = Q("9.3cm")
+ArapucaAcceptanceWindow_z  = Q("46.8cm")
+gapCenter_arapuca_z        = Q("16.0cm")
+gapBetween_arapuca_z       = Q("2.0cm")
+
+LightPaddle_x              = Q('3.3cm')
+LightPaddle_y              = Q('10.16cm')
+LightPaddle_z              = Q('219.8425cm')
+nLightPaddlePerAPA         = {Cryostat:nLightPaddlePerAPA}
+LightPaddleHeadOffset      = Q('4in')+Q('181.5mm')+Q('68mm')
+LightPaddleVerticalSpacing = Q('592mm')
+NoOpticalDetectors         = False
+
+
+[TPC]
+subbuilders      = ['TPCPlaneZ', 'TPCPlaneV', 'TPCPlaneU']
+class            = TPC.TPCBuilder
+driftDistance    = Q('3.6m')
+wirePlanePitch   = Q('0.476cm')
+APAGap_y         = Q('0.4cm')
+APAGap_z         = Q('0.8cm')
+cryoInnerDim     = {Cryostat:CryostatInnerDim}
+APAFrameDim      = {APAFrame:size}
+cathodeThickness = {Cryostat:cathodeThickness}
+nAPAs            = {Cryostat:nAPAs}
+
+
+[TPCOuter] # for when the apa is on the outside
+subbuilders      = {TPC:subbuilders}
+class            = {TPC:class}
+driftDistance    = Q('20cm')-Q('2.0282cm')-Q('0.4436cm')-Q('0.88cm') # "drift" is the dist to cryo wall
+# driftDistance    = Q('159.52mm')-Q('0.489754cm')-Q('0.203846cm') # "drift" is the dist to cryo wall
+wirePlanePitch   = {TPC:wirePlanePitch}
+APAGap_y         = {TPC:APAGap_y}
+APAGap_z         = {TPC:APAGap_z}
+cryoInnerDim     = {Cryostat:CryostatInnerDim}
+APAFrameDim      = {APAFrame:size}
+cathodeThickness = {Cryostat:cathodeThickness}
+nAPAs            = {Cryostat:nAPAs}
+
+
+[Cryostat]
+subbuilders		= ['TPC', 'TPCOuter', 'APAFrame']
+class			= Cryostat.CryostatBuilder
+CryostatInnerDim	= [Q('15100mm'),Q('14000mm'),Q('62000mm')]
+membraneThickness	= Q('800mm') #https://edms.cern.ch/ui/file/1834156/1/CENBNFCR0075.pdf
+# cathodeThickness	= Q('0.016cm')
+cathodeThickness	= Q('0.13in')
+# nAPAs			= [1, 2, 6]
+nAPAs			= [3, 2, 25]
+IgnoredAPAs		= [[]]
+outerAPAs		= True
+sideLAr			= Q('1cm')
+APAToFloor		= Q('49.2cm')
+LArLevel		= Q('13m')
+APAToUpstreamWall	= Q('301.2cm')
+APAToDownstreamWall	= Q('49.2cm')
+Layer1Thickness		= Q('0.12cm')
+Layer2Thickness		= Q('0.024m')
+Layer3Thickness		= Q('0.7748m')
+Layer1Material		= "SS304L"		
+Layer2Material		= "Wood"	
+Layer3Material		= "PolyurethaneFoam"
+WarmSkinMaterial        = "S460ML"
+DSSClearance		= [Q('402mm'),Q('400mm'),Q('1403mm')-Q('1600mm')/2] 
+nDSSBeam		= 0
+DSSBeamHeight		= Q('203.2mm')
+DSSBeamBase		= Q('101.6mm')
+DSSBeamThickW		= Q('6.9mm')
+DSSBeamThickF		= Q('10.4mm')
+#https://www.structural-drafting-net-expert.com/steel-sections-Europe-HL.html
+# or CENBNFCR0122.pdf
+IPEBeamHeight		= Q('1108mm')
+IPEBeamBase		= Q('402mm')
+IPEBeamThickW		= Q('22mm')
+IPEBeamThickF		= Q('40mm')
+#  CENBNFCR0129.pdf
+IPEBeamRoofHeight	= Q('600mm')
+IPEBeamRoofBase		= Q('300mm')
+IPEBeamRoofThickW	= Q('15.5mm')
+IPEBeamRoofThickF	= Q('30mm')
+# elevation of the transverse beams
+# CENBNFCR0080.pdf in
+BeamFloors		= [Q('2002mm')+Q('3608mm'), Q('6002mm')+Q('3608mm'), Q('10002mm')+Q('3608mm')]
+HoleDiam		= Q('800mm')
+TopBeam			= [0,2,3,5,6,8]
+#numbers of beams (same ref)
+nBeamX			= 39
+BeamSeparationX		= Q('1600mm')
+nBeamY			= 39
+BeamSeparationY		= Q('1600mm')
+nBeamZ			= 9
+BeamSeparationZ		= Q('1600mm')
+BeamMaterial            = "RadioS460ML"
+# CENBNFCR0214.pdf
+SteelThickness		= Q('10mm')
+nLightPaddlePerAPA	= 10
+
+FieldCageBarWidth	= Q("1.8125in")
+FieldCageBarHeight	= Q("0.4375in")
+FieldCageMaterial       = "ALUMINUM_AL"
+
+doWaterShielding	= False
+waterThickness		= Q("50cm")
+
+
+[DetEnclosureLAr]
+#https://docs.dunescience.org/cgi-bin/private/RetrieveFile?docid=14242&filename=ARUP-DWG-000006.00.IFI.00.01.pdf&version=1
+# Enjoy if you have to read that. Page 53.
+subbuilders     = ['Cryostat']
+class           = DetEncLAr.DetEncLArBuilder
+# detEncDim       = [Q('19.80m'),Q('23.15m'),Q('3.3m')+Q('60.45m')+Q('5.55m')+Q('12m')+Q('5.55m')+Q('60.45m')+Q('3.3m')-Q('84m')]
+detEncDim       = [Q('19.80m'),Q('23.15m'),Q('3.3m')+Q('60.45m')+Q('5.55m')+Q('12m')+Q('5.55m')+Q('60.45m')+Q('3.3m')]
+archRadius      = Q('12.84m')
+archHalfAngle   = Q('50deg')
+
+# https://docs.dunescience.org/cgi-bin/private/RetrieveFile?docid=464&filename=LBNF%20-%20Underground%20-%20North%20Cavern%20Clearance%20Envelopes%20-%202018-12-07.pdf&version=10
+# these numbers come from here. This document is not quite compatible with 100% CF but...
+# note the y number means the detector is on the floor
+# the x number is the number it *should be* if everything has the perfect size and the detector is in the centre
+# y and z are actually used.
+ConcreteBeamGap = [Q('427mm'),Q('0mm'),Q('427mm')]
+RadioRockThickness = Q('200cm')
+RadioRockMaterial  = 'RadioAverageDuneRock'
+ShotCreteThickness = Q('4in')
+ShotCreteMaterial  = "RadioShotcretePeteLein"
+ConcreteThickness  = Q('11in')
+ConcreteMaterial   = "RadioConcretePeteLein"
+GroutThickness     = Q('1in')
+GroutMaterial      = "RadioConcretePeteLein"
+
+[World]
+subbuilders         = ['DetEnclosureLAr']
+class               = World.WorldBuilder
+worldMat            = 'AverageDuneRock'
+worldDim            = [Q('400m'),Q('400m'),Q('400m')] # enough room for CRY, otherwise arbitrary
+encBackWallToHall_z = Q('46.25ft')