python/bin/g4GraphicalScan.py

0001 #!/usr/bin/env python3
0002 # ==========================================================================
0003 #  AIDA Detector description implementation
0004 # --------------------------------------------------------------------------
0005 # Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0006 # All rights reserved.
0007 #
0008 # For the licensing terms see $DD4hepINSTALL/LICENSE.
0009 # For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0010 #
0011 # ==========================================================================
0012 """
0013 Make a 2d slice through a dd4hep detector model.
0014 Present results as a series of TH2F
0015
0016 - density (1/X0), materials
0017
0018 in each bin, we consider the material along two lines,
0019 parallel to the histogram axes and passing through the bin centre.
0020
0021 The G4 geometry is scanned by shooting a geantino along various paths
0022
0023 D. Jeans, KEK. 2025/2/3
0024
0025 for usage instructions:
0026
0027   python3 g4GraphicalScan.py -h
0028
0029 e.g. for a scan at z=1000mm, in the range -10mm < x,y < 10mm, with 100x100 bins:
0030
0031   python3 g4GraphicalScan.py -c myModel.xml -s XY -x -10,10 -y -10,10 -z 1000 -n 100 -o scanOutput.root
0032 """
0033 import os
0034 import sys
0035 import optparse
0036 import subprocess
0037 import ROOT
0038
0039 # define the input parameters
0040
0041 parser = optparse.OptionParser()
0042 parser.formatter.width = 132
0043 parser.description = '2-dimensional material scan using Geant4.'
0044 parser.add_option('-c', '--compact', dest='compact', default=None,
0045                   help='compact xml input file',
0046                   metavar='<FILE>')
0047 parser.add_option('-s', '--sliceType',
0048                   dest='sliceType', default='ZX',
0049                   help='slice plane [XY, ZX, or ZY]',
0050                   metavar='<string>')
0051 parser.add_option('-x', '--xRange',
0052                   dest='xRange', default='-1000.,1000',
0053                   help='range to scan in x [in mm; give tuple "min,max" as string, or just "val" in case of ZY]',
0054                   metavar='<tuple>')
0055 parser.add_option('-y', '--yRange',
0056                   dest='yRange', default='-1000.,1000',
0057                   help='range to scan in y [in mm; give tuple "min,max" as string, or just "val" in case of ZX]',
0058                   metavar='<tuple>')
0059 parser.add_option('-z', '--zRange',
0060                   dest='zRange', default='-1000.,1000',
0061                   help='range to scan in z [in mm; give tuple "min,max" as string, or just "val" in case of XY]',
0062                   metavar='<tuple>')
0063 parser.add_option('-n', '--nBins',
0064                   dest='nBins', default='100',
0065                   help='number of bins in output histograms',
0066                   metavar='<int>')
0067 parser.add_option('-o', '--outputFile',
0068                   dest='outFile', default='output.root',
0069                   help='name of ouput root file',
0070                   metavar='<string>')
0071 parser.add_option("-P", "--noPilot",
0072                   action="store_true", dest="noPilot", default=False,
0073                   help="don't run the pilot job (e.g. if you're sure the geometry is good)")
0074
0075 (opts, args) = parser.parse_args()
0076 #
0077 # check that the requested inputs are valid
0078 #
0079 infileName = str(opts.compact)
0080 if not os.path.isfile(infileName):
0081     print('ERROR: cannot find requested input geometry file', infileName, file=sys.stderr)
0082     exit(1)
0083 print(infileName)
0084
0085 sliceType = str(opts.sliceType)
0086 if sliceType != 'XY' and sliceType != 'ZX' and sliceType != 'ZY':
0087     print('ERROR: unknown slice Type', sliceType, '. Choose XY, ZX or ZY.', file=sys.stderr)
0088     exit(1)
0089 print(sliceType)
0090
0091 planePos = -99999.
0092 planeAxis = ''
0093
0094 aa = str(opts.xRange).split(',')
0095 if len(aa) == 2 and sliceType != 'ZY':
0096     xRange = (float(aa[0]), float(aa[1]))
0097     if xRange[1] <= xRange[0]:
0098         print('ERROR, xmin is larger than xmax', file=sys.stderr)
0099         exit(1)
0100 elif len(aa) == 1 and sliceType == 'ZY':
0101     xRange = (float(aa[0]))
0102     planePos = xRange
0103     planeAxis = 'X'
0104 else:
0105     print('ERROR: could not determine xRange, or inconsistent with sliceType', file=sys.stderr)
0106     exit(1)
0107 print('xRange', xRange)
0108
0109 aa = str(opts.yRange).split(',')
0110 if len(aa) == 2 and sliceType != 'ZX':
0111     yRange = (float(aa[0]), float(aa[1]))
0112     if yRange[1] <= yRange[0]:
0113         print('ERROR, ymin is larger than ymax', file=sys.stderr)
0114         exit(1)
0115 elif len(aa) == 1 and sliceType == 'ZX':
0116     yRange = (float(aa[0]))
0117     planePos = yRange
0118     planeAxis = 'Y'
0119 else:
0120     print('ERROR: could not determine yRange, or inconsistent with sliceType', file=sys.stderr)
0121     exit(1)
0122 print('yRange', yRange)
0123
0124 aa = str(opts.zRange).split(',')
0125 if len(aa) == 2 and sliceType != 'XY':
0126     zRange = (float(aa[0]), float(aa[1]))
0127     if zRange[1] <= zRange[0]:
0128         print('ERROR, zmin is larger than zmax', file=sys.stderr)
0129         exit(1)
0130 elif len(aa) == 1 and sliceType == 'XY':
0131     zRange = (float(aa[0]))
0132     planePos = zRange
0133     planeAxis = 'Z'
0134 else:
0135     print('ERROR: could not determine zRange, or inconsistent with sliceType', file=sys.stderr)
0136     exit(1)
0137 print('zRange', zRange)
0138
0139 nBins = int(opts.nBins)
0140 print('nBins', nBins)
0141
0142 if nBins < 1:
0143     print('ERROR: crazy number of bins requested', nBins, file=sys.stderr)
0144     exit(1)
0145
0146 noPilot = bool(opts.noPilot)
0147 print('noPilot', noPilot)
0148
0149 outFileName = str(opts.outFile)
0150 #
0151 # define the "mother" histogram according to the requested slice type, ranges, number of bins
0152 #
0153 fout = ROOT.TFile(outFileName, 'recreate')
0154
0155 if sliceType == 'XY':
0156     h2 = ROOT.TH2F('hMat', 'h2', nBins, xRange[0], xRange[1], nBins, yRange[0], yRange[1])
0157     h2.GetXaxis().SetTitle('x [mm]')
0158     h2.GetYaxis().SetTitle('y [mm]')
0159 elif sliceType == 'ZX':
0160     h2 = ROOT.TH2F('hMat', 'h2', nBins, zRange[0], zRange[1], nBins, xRange[0], xRange[1])
0161     h2.GetXaxis().SetTitle('z [mm]')
0162     h2.GetYaxis().SetTitle('x [mm]')
0163 elif sliceType == 'ZY':
0164     h2 = ROOT.TH2F('hMat', 'h2', nBins, zRange[0], zRange[1], nBins, yRange[0], yRange[1])
0165     h2.GetXaxis().SetTitle('z [mm]')
0166     h2.GetYaxis().SetTitle('y [mm]')
0167 h2.Fill(0, 0, 0.0)  # to ensure there is at least one entry...otherwise doesn't get drawn...
0168
0169 #
0170 # this is where the materials will be stored
0171 #
0172 mats = {}
0173 #
0174 # we make scans along the X and Y axes of the mother histogram
0175 # prepare the input macro to ddsim
0176 #
0177 pilotName = '_pilot_' + outFileName + '.mac'
0178 pilotMac = open(pilotName, 'w')
0179 pilotMac.write('/gun/particle geantino' + '\n')
0180 pilotMac.write('/gun/energy 20 GeV' + '\n')
0181 pilotMac.write('/gun/number 1' + '\n')
0182
0183 steerName = '_' + outFileName + '.mac'
0184 steerMac = open(steerName, 'w')
0185 steerMac.write('/gun/particle geantino' + '\n')
0186 steerMac.write('/gun/energy 20 GeV' + '\n')
0187 steerMac.write('/gun/number 1' + '\n')
0188
0189 requestedStartPositions = {}
0190
0191 for iDir in range(0, 2):
0192     npilot = 0
0193     mats[iDir] = {}
0194     requestedStartPositions[iDir] = []
0195     if iDir == 0:
0196         axis = h2.GetXaxis()
0197     else:
0198         axis = h2.GetYaxis()
0199     #
0200     # the direction of the gun
0201     #
0202     if sliceType == 'XY':
0203         if iDir == 0:
0204             dirn = '0 1 0'
0205         else:
0206             dirn = '1 0 0'
0207     elif sliceType == 'ZX':
0208         if iDir == 0:
0209             dirn = '1 0 0'
0210         else:
0211             dirn = '0 0 1'
0212     elif sliceType == 'ZY':
0213         if iDir == 0:
0214             dirn = '0 1 0'
0215         else:
0216             dirn = '0 0 1'
0217     steerMac.write('/gun/direction ' + dirn + '\n')
0218     #
0219     # loop over the bins in this axis
0220     #
0221     for iX in range(1, nBins + 1):
0222
0223         mats[iDir][iX] = {}
0224         #
0225         # define the starting position of the gun
0226         #
0227         X = axis.GetBinCenter(iX)
0228         if iDir == 0:
0229             if sliceType == 'XY':
0230                 startPos = str(X) + ' '
0231                 startPos += str(yRange[0]) + ' '
0232                 startPos += str(zRange)
0233             elif sliceType == 'ZX':
0234                 startPos = str(xRange[0]) + ' '
0235                 startPos += str(yRange) + ' '
0236                 startPos += str(X)
0237             elif sliceType == 'ZY':
0238                 startPos = str(xRange) + ' '
0239                 startPos += str(yRange[0]) + ' '
0240                 startPos += str(X)
0241         else:
0242             if sliceType == 'XY':
0243                 startPos = str(xRange[0]) + ' '
0244                 startPos += str(X) + ' '
0245                 startPos += str(zRange)
0246             elif sliceType == 'ZX':
0247                 startPos = str(X) + ' '
0248                 startPos += str(yRange) + ' '
0249                 startPos += str(zRange[0])
0250             elif sliceType == 'ZY':
0251                 startPos = str(xRange) + ' '
0252                 startPos += str(X) + ' '
0253                 startPos += str(zRange[0])
0254
0255         steerMac.write('/gun/position ' + startPos + ' mm \n')
0256         steerMac.write('/run/beamOn' + '\n')
0257         if npilot < 1:
0258             pilotMac.write('/gun/position ' + startPos + ' mm \n')
0259             pilotMac.write('/run/beamOn' + '\n')
0260             npilot += 1
0261         requestedStartPositions[iDir].append(startPos)
0262
0263 steerMac.write('exit')
0264 steerMac.close()
0265
0266 pilotMac.write('exit')
0267 pilotMac.close()
0268 #
0269 # first try a pilot run to check model is OK
0270 #  pilot jobs has 2 events, one in each direction
0271 #
0272 if not noPilot:
0273     cmd = ['ddsim', '--compactFile', infileName, '--runType', 'run', '--enableG4Gun',
0274            '--action.step', 'Geant4MaterialScanner/MaterialScan', '-M', pilotName]
0275     print('running test pilot job...\n')
0276     for cc in cmd:
0277         print(cc, end=' ')
0278     print('\n')
0279     pilotresult = subprocess.run(cmd, capture_output=True, text=True, timeout=90)
0280     print('done, checking pilot result')
0281     has_Material_scan_between = 0
0282     has_Finished_run = 0
0283     for ll in pilotresult.stdout.splitlines():
0284         if 'Material scan between' in ll:
0285             has_Material_scan_between += 1
0286         if 'Finished run' in ll:
0287             has_Finished_run += 1
0288     if has_Material_scan_between != 2 or has_Finished_run != 2:
0289         print('ERROR, pilot job seems not to have finished successfully')
0290         for ll in pilotresult.stdout.splitlines():
0291             print(ll)
0292         print('ERROR, pilot job seems not to have finished successfully')
0293     print('pilot job seems OK')
0294 #
0295 # run ddsim with the full macro
0296 #
0297 cmd = ['ddsim', '--compactFile', infileName, '--runType', 'run', '--enableG4Gun',
0298        '--action.step', 'Geant4MaterialScanner/MaterialScan', '-M', steerName]
0299 print('now running main ddsim job..this may take some time')
0300 result = subprocess.run(cmd, capture_output=True, text=True, timeout=600)
0301 #
0302 # parse the results
0303 #
0304 iscan = 1
0305 inScan = False
0306 iDir = 0
0307 for line in result.stdout.splitlines():
0308     if 'Material scan between' in line:
0309         gg = line.split(')')[0].split('(')[1].split(',')
0310         startx = 10 * float(gg[0])  # convert cm -> mm
0311         starty = 10 * float(gg[1])
0312         startz = 10 * float(gg[2])
0313         inScan = True
0314         # check consistency with requested position.
0315         # if a gun position outside the world volume is requested,
0316         #  this can cause an inconsistency (it is started at 0,0,0)
0317         pp = requestedStartPositions[iDir][iscan - 1].split()
0318         rx = float(pp[0])
0319         ry = float(pp[1])
0320         rz = float(pp[2])
0321         if abs(rx - startx) > 1. or abs(ry - starty) > 1. or abs(rz - startz) > 1.:
0322             print('ERROR inconsistent starting gun position')
0323             print('  REQUESTED:', pp)
0324             print('  USED:', gg)
0325             print('The requested range probably lies partially outside the world volume')
0326             print('  use a more reasonable range and try again!')
0327             exit(1)
0328     elif 'Finished run' in line:
0329         iscan += 1
0330         if iscan == nBins + 1:   # now move to the second set of scans
0331             iDir = 1
0332             iscan = 1
0333         inScan = False
0334     elif r"+-----------------" in line:  # comment line
0335         continue
0336     elif r"|     \   Material" in line:  # comment line
0337         continue
0338     elif r"| Num. \  Name" in line:      # comment line
0339         continue
0340     elif r"| Layer \ " in line:          # comment line
0341         continue
0342     elif inScan and \
0343          '(' in line and \
0344          len(line.split('(')[0].split()) == 12 and \
0345          line.split()[0] == '|':  # this line contains material information
0346         index = int(line.split()[1])
0347         material = line.split()[2]
0348         radlen = 10 * float(line.split()[6])     # cm->mm
0349         thickness = 10 * float(line.split()[8])  # cm->mm
0350         endpos = line.split('(')[1].split(')')[0].split(',')
0351         endx = 10 * float(endpos[0])             # cm -> mm
0352         endy = 10 * float(endpos[1])
0353         endz = 10 * float(endpos[2])
0354         mats[iDir][iscan][index] = [material, radlen, thickness, endx, endy, endz]
0355 #
0356 # now all data is collected: fill the histograms
0357 #
0358 h2.SetTitle('materialScan at ' + planeAxis + '=' + str(planePos) + ' mm : 1/X_{0}')
0359
0360 hists = {}
0361 hists['x0'] = h2
0362
0363 for iDir in range(0, 2):   # the two directions
0364     if iDir == 1:
0365         mainaxis = h2.GetYaxis()  # perpendicular to the scan direction
0366         scanaxis = h2.GetXaxis()  # parallel to the scan direction
0367     elif iDir == 0:
0368         mainaxis = h2.GetXaxis()
0369         scanaxis = h2.GetYaxis()
0370
0371     for jj in range(1, mainaxis.GetNbins() + 1):  # loop over the scan lines
0372         mainaxis.GetBinCenter(jj)
0373         scandat = mats[iDir][jj]
0374
0375         hxn = scanaxis.GetNbins()
0376         hxl = scanaxis.GetBinLowEdge(1)
0377         hxh = scanaxis.GetBinUpEdge(hxn)
0378         hxbw = scanaxis.GetBinWidth(1)
0379
0380         curpos = hxl
0381
0382         for value in scandat.values():
0383             begpos = curpos
0384
0385             endx = value[3]
0386             endy = value[4]
0387             endz = value[5]
0388
0389             if sliceType == 'XY':
0390                 if iDir == 0:
0391                     endpos = endy
0392                 elif iDir == 1:
0393                     endpos = endx
0394             elif sliceType == 'ZX':
0395                 if iDir == 0:
0396                     endpos = endx
0397                 elif iDir == 1:
0398                     endpos = endz
0399             elif sliceType == 'ZY':
0400                 if iDir == 0:
0401                     endpos = endy
0402                 elif iDir == 1:
0403                     endpos = endz
0404
0405             radlen = value[1]
0406             thick = value[2]
0407             matStr = value[0]
0408
0409             if begpos < hxl and endpos < hxl:  # not in histo range (below)
0410                 pass
0411             elif begpos > hxh and endpos > hxh:  # not in histo range (above)
0412                 pass
0413             else:
0414                 if matStr not in hists.keys():  # not yet seen material: make a new histogram for it
0415                     hists[matStr] = h2.Clone(h2.GetName() + '_' + matStr)
0416                     hists[matStr].SetTitle(hists[matStr].GetTitle().replace('1/X_{0}', matStr))
0417                     hists[matStr].Reset()
0418                     hists[matStr].Fill(0, 0, 0.0)  # ensure at least one entry...otherwise doesn't get drawn..
0419
0420                 iy1 = scanaxis.FindBin(begpos)
0421                 iy2 = scanaxis.FindBin(endpos)
0422
0423                 if iy1 == iy2:   # this step entirely within one histo bin
0424                     if iDir == 1:
0425                         hists['x0']  .AddBinContent(iy1, jj, thick / radlen / hxbw)
0426                         hists[matStr].AddBinContent(iy1, jj, thick / hxbw)
0427                     else:
0428                         hists['x0']  .AddBinContent(jj, iy1, thick / radlen / hxbw)
0429                         hists[matStr].AddBinContent(jj, iy1, thick / hxbw)
0430                 else:   # this step extends over two or more bins
0431                     firstBinStubLength = scanaxis.GetBinUpEdge(iy1) - begpos
0432                     lastBinStubLength = endpos - scanaxis.GetBinLowEdge(iy2)
0433                     if iDir == 1:
0434                         hists['x0'].AddBinContent(iy1, jj, firstBinStubLength / radlen / hxbw)
0435                         hists['x0'].AddBinContent(iy2, jj, lastBinStubLength / radlen / hxbw)
0436                         hists[matStr].AddBinContent(iy1, jj, firstBinStubLength / hxbw)
0437                         hists[matStr].AddBinContent(iy2, jj, lastBinStubLength / hxbw)
0438                     else:
0439                         hists['x0'].AddBinContent(jj, iy1, firstBinStubLength / radlen / hxbw)
0440                         hists['x0'].AddBinContent(jj, iy2, lastBinStubLength / radlen / hxbw)
0441                         hists[matStr].AddBinContent(jj, iy1, firstBinStubLength / hxbw)
0442                         hists[matStr].AddBinContent(jj, iy2, lastBinStubLength / hxbw)
0443
0444                     if iy2 - iy1 > 1:  # fill the bins in between first and last
0445                         for i in range(iy1 + 1, iy2):
0446                             if iDir == 1:
0447                                 hists['x0']  .AddBinContent(i, jj, 1. / radlen)
0448                                 hists[matStr].AddBinContent(i, jj, 1.)
0449                             else:
0450                                 hists['x0']  .AddBinContent(jj, i, 1. / radlen)
0451                                 hists[matStr].AddBinContent(jj, i, 1.)
0452             curpos = endpos
0453
0454 for mm in hists.keys():
0455     hists[mm].Scale(1. / 2)  # average of the two direction scans
0456
0457 print('done filling histograms, now closing root file')
0458 fout.Write()
0459 fout.Close()
0460
0461 # clean up the macro files
0462 os.remove(steerName)
0463 os.remove(pilotName)