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