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