Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:10:58

0001 #!/usr/bin/env python3
0002 
0003 import os
0004 import argparse
0005 import time
0006 from datetime import datetime
0007 
0008 from orion.client import build_experiment
0009 from orion.storage.base import get_storage
0010 
0011 import acts
0012 from acts import (
0013     SurfaceMaterialMapper,
0014     VolumeMaterialMapper,
0015     Navigator,
0016     Propagator,
0017     StraightLineStepper,
0018     MaterialMapJsonConverter,
0019 )
0020 from acts.examples import (
0021     Sequencer,
0022     WhiteBoard,
0023     AlgorithmContext,
0024     ProcessCode,
0025     RootMaterialTrackReader,
0026     MaterialMapping,
0027     JsonMaterialWriter,
0028     JsonFormat,
0029 )
0030 from acts.examples.odd import getOpenDataDetector
0031 
0032 
0033 def runMaterialMappingNoTrack(
0034     trackingGeometry,
0035     decorators,
0036     outputDir,
0037     inputDir,
0038     mapName="material-map",
0039     mapSurface=True,
0040     mapVolume=True,
0041     format=JsonFormat.Json,
0042     readCachedSurfaceInformation=False,
0043     s=None,
0044 ):
0045     """
0046     Implementation of the material mapping that doesn't write the material tracks.
0047     Used to create the material map that will then be used to compute the material variance.
0048 
0049     trackingGeometry : The tracking geometry
0050     decorators : The decorators for the tracking geometry
0051     outputDir : Output directory for the material map
0052     inputDir : Input directory containing the Material track
0053     mapName : Name of the material map
0054     mapSurface : Is material being mapped onto surfaces ?
0055     mapVolume : Is material being mapped onto volumes ?
0056     format : Json format used to write the material map (json, cbor, ...)
0057     readCachedSurfaceInformation : If set to true it will be assumed that the surface has already been associated with each material interaction in the input file.
0058     """
0059 
0060     s = s or Sequencer(numThreads=1)
0061     for decorator in decorators:
0062         s.addContextDecorator(decorator)
0063     wb = WhiteBoard(acts.logging.INFO)
0064     context = AlgorithmContext(0, 0, wb)
0065     for decorator in decorators:
0066         assert decorator.decorate(context) == ProcessCode.SUCCESS
0067 
0068     # Read material step information from a ROOT TTRee
0069     s.addReader(
0070         RootMaterialTrackReader(
0071             level=acts.logging.INFO,
0072             outputMaterialTracks="material-tracks",
0073             fileList=[
0074                 os.path.join(
0075                     inputDir,
0076                     "optimised-material-map_tracks.root"
0077                     if readCachedSurfaceInformation
0078                     else "geant4_material_tracks.root",
0079                 )
0080             ],
0081             readCachedSurfaceInformation=readCachedSurfaceInformation,
0082         )
0083     )
0084 
0085     stepper = StraightLineStepper()
0086     mmAlgCfg = MaterialMapping.Config(context.geoContext, context.magFieldContext)
0087     mmAlgCfg.trackingGeometry = trackingGeometry
0088     mmAlgCfg.inputMaterialTracks = "material-tracks"
0089 
0090     if mapSurface:
0091         navigator = Navigator(
0092             trackingGeometry=trackingGeometry,
0093             resolveSensitive=True,
0094             resolveMaterial=True,
0095             resolvePassive=True,
0096         )
0097         propagator = Propagator(stepper, navigator)
0098         mapper = SurfaceMaterialMapper(level=acts.logging.INFO, propagator=propagator)
0099         mmAlgCfg.materialSurfaceMapper = mapper
0100 
0101     if mapVolume:
0102         navigator = Navigator(
0103             trackingGeometry=trackingGeometry,
0104         )
0105         propagator = Propagator(stepper, navigator)
0106         mapper = VolumeMaterialMapper(
0107             level=acts.logging.INFO, propagator=propagator, mappingStep=999
0108         )
0109         mmAlgCfg.materialVolumeMapper = mapper
0110 
0111     jmConverterCfg = MaterialMapJsonConverter.Config(
0112         processSensitives=True,
0113         processApproaches=True,
0114         processRepresenting=True,
0115         processBoundaries=True,
0116         processVolumes=True,
0117         context=context.geoContext,
0118     )
0119 
0120     jmw = JsonMaterialWriter(
0121         level=acts.logging.VERBOSE,
0122         converterCfg=jmConverterCfg,
0123         fileName=os.path.join(outputDir, mapName),
0124         writeFormat=format,
0125     )
0126 
0127     mmAlgCfg.materialWriters = [jmw]
0128 
0129     s.addAlgorithm(MaterialMapping(level=acts.logging.INFO, config=mmAlgCfg))
0130 
0131     return s
0132 
0133 
0134 def runMaterialMappingVariance(
0135     binMap,
0136     events,
0137     job,
0138     inputPath,
0139     pathExp,
0140     pipeResult,
0141     readCachedSurfaceInformation=False,
0142 ):
0143     """
0144     Run the material mapping and compute the variance for each bin of each surface
0145     Return a dict with the GeometryId value of the surface as a key that stores
0146     a list of pairs corresponding to the variance and number of tracks associated with each bin of the surface
0147 
0148     binMap : Map containing the binning for each surface
0149     events : Number of event to use in the mapping
0150     job : ID of the job
0151     inputPath : Directory containing the input geantino track and the json geometry
0152     pathExp : Material mapping optimisation path
0153     pipeResult : Pipe to send back the score to the main python instance
0154     readCachedSurfaceInformation : Are surface information stored in the material track. Switch to true if the mapping was already performed to improve the speed.
0155     """
0156     print(
0157         datetime.now().strftime("%H:%M:%S") + "    Start mapping for job " + str(job),
0158         flush=True,
0159     )
0160     mapName = "material-map-" + str(job)
0161     mapSurface = True
0162     mapVolume = False
0163 
0164     # Create a MappingMaterialDecorator based on the tracking geometry
0165     matDeco = acts.IMaterialDecorator.fromFile(
0166         str(os.path.join(inputPath, "geometry-map.json"))
0167     )
0168     detectorTemp, trackingGeometryTemp, decoratorsTemp = getOpenDataDetector(matDeco)
0169     matMapDeco = acts.MappingMaterialDecorator(
0170         tGeometry=trackingGeometryTemp, level=acts.logging.ERROR
0171     )
0172     # Update the binning using the bin map corresponding to this trial
0173     matMapDeco.setBinningMap(binMap)
0174 
0175     del detectorTemp
0176     del trackingGeometryTemp
0177     del decoratorsTemp
0178 
0179     # Decorate the detector with the MappingMaterialDecorator
0180     detector, trackingGeometry, decorators = getOpenDataDetector(matMapDeco)
0181 
0182     # Sequence for the mapping, only use one thread when mapping material
0183     sMap = acts.examples.Sequencer(
0184         events=events, numThreads=1, logLevel=acts.logging.INFO
0185     )
0186 
0187     # Run the material mapping without writing the material track (as they take too much space)
0188     runMaterialMappingNoTrack(
0189         trackingGeometry,
0190         decorators,
0191         outputDir=pathExp,
0192         inputDir=inputPath,
0193         mapName=mapName,
0194         format=JsonFormat.Cbor,
0195         mapVolume=mapVolume,
0196         readCachedSurfaceInformation=readCachedSurfaceInformation,
0197         s=sMap,
0198     )
0199     sMap.run()
0200     del sMap  # Need to be deleted to write the material map to cbor
0201     del detector
0202     del trackingGeometry
0203     del decorators
0204 
0205     # Compute the variance by rerunning the mapping
0206     print(
0207         datetime.now().strftime("%H:%M:%S")
0208         + "    Job "
0209         + str(job)
0210         + ": second pass to compute the variance",
0211         flush=True,
0212     )
0213     # Use the material map from the previous mapping as an input
0214     cborMap = os.path.join(pathExp, (mapName + ".cbor"))
0215     matDecoVar = acts.IMaterialDecorator.fromFile(cborMap)
0216     detectorVar, trackingGeometryVar, decoratorsVar = getOpenDataDetector(matDecoVar)
0217     s = acts.examples.Sequencer(events=events, numThreads=1, logLevel=acts.logging.INFO)
0218     for decorator in decoratorsVar:
0219         s.addContextDecorator(decorator)
0220     wb = WhiteBoard(acts.logging.ERROR)
0221     context = AlgorithmContext(0, 0, wb)
0222     for decorator in decoratorsVar:
0223         assert decorator.decorate(context) == ProcessCode.SUCCESS
0224 
0225     # Read material step information from a ROOT TTRee
0226     reader = RootMaterialTrackReader(
0227         level=acts.logging.ERROR,
0228         outputMaterialTracks="material-tracks",
0229         fileList=[
0230             os.path.join(
0231                 inputPath,
0232                 "optimised-material-map_tracks.root"
0233                 if readCachedSurfaceInformation
0234                 else "geant4_material_tracks.root",
0235             )
0236         ],
0237         readCachedSurfaceInformation=readCachedSurfaceInformation,
0238     )
0239     s.addReader(reader)
0240     stepper = StraightLineStepper()
0241     mmAlgCfg = MaterialMapping.Config(context.geoContext, context.magFieldContext)
0242     mmAlgCfg.trackingGeometry = trackingGeometryVar
0243     mmAlgCfg.inputMaterialTracks = "material-tracks"
0244 
0245     if mapSurface:
0246         navigator = Navigator(
0247             trackingGeometry=trackingGeometryVar,
0248             resolveSensitive=True,
0249             resolveMaterial=True,
0250             resolvePassive=True,
0251         )
0252         propagator = Propagator(stepper, navigator)
0253         surfaceCfg = SurfaceMaterialMapper.Config(
0254             computeVariance=True
0255         )  # Don't forget to turn the `computeVariance` to true
0256         mapper = SurfaceMaterialMapper(
0257             config=surfaceCfg, level=acts.logging.ERROR, propagator=propagator
0258         )
0259         mmAlgCfg.materialSurfaceMapper = mapper
0260 
0261     if mapVolume:
0262         navigator = Navigator(
0263             trackingGeometry=trackingGeometryVar,
0264         )
0265         propagator = Propagator(stepper, navigator)
0266         mapper = VolumeMaterialMapper(
0267             level=acts.logging.ERROR, propagator=propagator, mappingStep=999
0268         )
0269         mmAlgCfg.materialVolumeMapper = mapper
0270 
0271     mapping = MaterialMapping(level=acts.logging.ERROR, config=mmAlgCfg)
0272     s.addAlgorithm(mapping)
0273     s.run()
0274 
0275     # Compute the scoring parameters
0276     score = dict()
0277     for key in binMap:
0278         objective = 0
0279         nonZero = 0
0280         binParameters = mapping.scoringParameters(key)
0281         # Objective : Sum of variance in all bin divided by the number of bin
0282         # The variance is scaled by (1.0 + 1.0/(number of hits in the bin)) to encourage larger bin at equal score
0283         for parameters in binParameters:
0284             if parameters[1] != 0:
0285                 objective += parameters[0] * (1.0 + 1.0 / parameters[1])
0286                 nonZero += 1
0287         if nonZero != 0:
0288             objective = objective / nonZero
0289         score[key] = [dict(name="surface_score", type="objective", value=objective)]
0290     print(
0291         datetime.now().strftime("%H:%M:%S")
0292         + "    Mapping over for job "
0293         + str(job)
0294         + " : now sending score",
0295         flush=True,
0296     )
0297     pipeResult.send(score)
0298 
0299     del mapping
0300     del s
0301     del detectorVar
0302     del trackingGeometryVar
0303     del decoratorsVar
0304     os.remove(cborMap)
0305 
0306 
0307 def surfaceExperiment(key, nbJobs, pathDB, pathResult, pipeBin, pipeResult, doPloting):
0308     """
0309     This function create an experiment for a given single surface
0310     Due to how Orion is implemented only one DB can exist per job, this thus need to be call using pythons multiprocessing to circumvent the issue.
0311 
0312     key : Id of the surface corresponding to this experiment
0313     nbJobs : Total number of jobs to be executed simultaneously
0314     pathDB : Path to the databases
0315     pathResult : Path to write the result of the optimisation
0316     pipeBin : Pipe use to send the experiment binning to the main python instance
0317     pipeResult : Pipe to receive the result of the optimisation
0318     doPloting : true if we want to plot the result of the optimisation and obtain the optimal material map
0319     """
0320     # Create the database
0321     storage = {
0322         "database": {
0323             "name": "database_" + str(key),
0324             "type": "pickleddb",
0325             "host": os.path.join(pathDB, "database_" + str(key) + ".pkl"),
0326             "timeout": 43200,
0327         },
0328     }
0329     # Create the search space, the range of the binning can be chosen here
0330     # x represent X or phi depending on the type of surface
0331     # y represent Y, R or Z depending on the type of surface
0332     space = {
0333         "x": "uniform(1, 240, discrete=True)",
0334         "y": "uniform(1, 240, discrete=True)",
0335     }
0336     # Build the experiment
0337     experiments = build_experiment(
0338         "s_" + str(key),
0339         version="1",
0340         space=space,
0341         algorithms="random",
0342         storage=storage,
0343         max_idle_time=43200,
0344     )
0345     # Clean trial that haven't been completed
0346     store = get_storage()
0347     store.delete_trials(uid=experiments.id, where={"status": {"$ne": "completed"}})
0348     store.release_algorithm_lock(uid=experiments.id)
0349 
0350     # Suggest one binning per job and then send them via the pipe
0351     trials = dict()
0352     binMap = dict()
0353     for job in range(nbJobs):
0354         trials[job] = experiments.suggest()
0355         binMap[job] = (trials[job].params["x"], trials[job].params["y"])
0356         print(
0357             datetime.now().strftime("%H:%M:%S")
0358             + "    Binning for job "
0359             + str(job)
0360             + " and surface "
0361             + str(key)
0362             + " has been selected",
0363             flush=True,
0364         )
0365         pipeBin.send(binMap[job])
0366     print(
0367         datetime.now().strftime("%H:%M:%S")
0368         + "    All binning for surface "
0369         + str(key)
0370         + " has been sent",
0371         flush=True,
0372     )
0373     # Store the score resulting for the jobs in the database
0374     for job in range(nbJobs):
0375         score = pipeResult.recv()
0376         print(
0377             datetime.now().strftime("%H:%M:%S")
0378             + "    Received score for job "
0379             + str(job)
0380             + " and surface "
0381             + str(key),
0382             flush=True,
0383         )
0384         experiments.observe(trials[job], score)
0385         print(
0386             datetime.now().strftime("%H:%M:%S")
0387             + "    Score for job "
0388             + str(job)
0389             + " and surface "
0390             + str(key)
0391             + " has been written",
0392             flush=True,
0393         )
0394 
0395     # Create some performances plots for each surface
0396     if doPloting:
0397         print(
0398             datetime.now().strftime("%H:%M:%S")
0399             + "    All the jobs are over. Now creating the optimisation plots",
0400             flush=True,
0401         )
0402 
0403         pathExpSurface = os.path.join(pathResult, "b_" + str(key))
0404         if not os.path.isdir(pathExpSurface):
0405             os.makedirs(pathExpSurface)
0406 
0407         regret = experiments.plot.regret()
0408         regret.write_html(pathExpSurface + "/regret.html")
0409 
0410         parallel_coordinates = experiments.plot.parallel_coordinates()
0411         parallel_coordinates.write_html(pathExpSurface + "/parallel_coordinates.html")
0412 
0413         lpi = experiments.plot.lpi()
0414         lpi.write_html(pathExpSurface + "/lpi.html")
0415 
0416         partial_dependencies = experiments.plot.partial_dependencies()
0417         partial_dependencies.write_html(pathExpSurface + "/partial_dependencies.html")
0418 
0419         # Select the optimal binning and send it via the pipe
0420         df = experiments.to_pandas()
0421         best = df.iloc[df.objective.idxmin()]
0422         print(
0423             datetime.now().strftime("%H:%M:%S")
0424             + "    Best score for surface "
0425             + str(key)
0426             + " : "
0427             + str(best),
0428             flush=True,
0429         )
0430         resultBinMap = (best.x, best.y)
0431         pipeBin.send(resultBinMap)
0432 
0433 
0434 if "__main__" == __name__:
0435     print(datetime.now().strftime("%H:%M:%S") + "    Starting")
0436     # Optimiser arguments
0437     parser = argparse.ArgumentParser()
0438     parser.add_argument(
0439         "--numberOfJobs", nargs="?", default=2, type=int
0440     )  # number of parallele jobs
0441     parser.add_argument(
0442         "--topNumberOfEvents", nargs="?", default=100, type=int
0443     )  # number of events per trials
0444     parser.add_argument(
0445         "--inputPath", nargs="?", default=os.getcwd(), type=str
0446     )  # path to the input
0447     parser.add_argument(
0448         "--outputPath", nargs="?", default="", type=str
0449     )  # path to the output
0450     parser.add_argument(
0451         "--doPloting", action="store_true"
0452     )  # Return the optimisation plot and create the optimal material map
0453     parser.add_argument(
0454         "--readCachedSurfaceInformation", action="store_true"
0455     )  # Use surface information from the material track
0456     parser.set_defaults(doPloting=False)
0457     parser.set_defaults(readCachedSurfaceInformation=False)
0458     args = parser.parse_args()
0459 
0460     # Define the useful path and create them if they do not exist
0461     pathExp = os.path.join(args.outputPath, "Mapping")
0462     pathDB = os.path.join(pathExp, "Database")
0463     pathResult = os.path.join(pathExp, "Result")
0464     if not os.path.isdir(pathExp):
0465         os.makedirs(pathExp)
0466     if not os.path.isdir(pathDB):
0467         os.makedirs(pathDB)
0468     if not os.path.isdir(pathResult):
0469         os.makedirs(pathResult)
0470 
0471     # Create the tracking geometry, uses the json file to configure the proto-surfaces
0472     matDeco = acts.IMaterialDecorator.fromFile(
0473         str(os.path.join(args.inputPath, "geometry-map.json"))
0474     )
0475     detector, trackingGeometry, decorators = getOpenDataDetector(matDeco)
0476 
0477     # Use the MappingMaterialDecorator to create a binning map that can be optimised
0478     matMapDeco = acts.MappingMaterialDecorator(
0479         tGeometry=trackingGeometry, level=acts.logging.WARNING
0480     )
0481     binDict = matMapDeco.binningMap()
0482     del detector, decorators
0483 
0484     # Create the pipes that will be used to transfer data to/from the jobs
0485     from multiprocessing import Process, Pipe
0486 
0487     binPipes_child = dict()
0488     resultPipes_child = dict()
0489     scorePipes_child = dict()
0490 
0491     binPipes_parent = dict()
0492     resultPipes_parent = dict()
0493     scorePipes_parent = dict()
0494 
0495     expJob = dict()
0496     OptiJob = dict()
0497 
0498     # Build one experiment per surface
0499     # The binning of the surfaces are independent so we split
0500     # an optimisation problem with a large number of variable into a lot of optimisation with 2
0501     for key in binDict:
0502         binPipes_parent[key], binPipes_child[key] = Pipe()
0503         scorePipes_parent[key], scorePipes_child[key] = Pipe()
0504         expJob[key] = Process(
0505             target=surfaceExperiment,
0506             args=(
0507                 key,
0508                 args.numberOfJobs,
0509                 pathDB,
0510                 pathResult,
0511                 binPipes_child[key],
0512                 scorePipes_child[key],
0513                 args.doPloting,
0514             ),
0515         )
0516         expJob[key].start()
0517 
0518     # Prepare `args.numberOfJobs` material mapping jobs
0519     for job in range(args.numberOfJobs):
0520         resultPipes_parent[job], resultPipes_child[job] = Pipe()
0521         binMap = dict()
0522         # Collect the binning for all the surfaces and create a bin map
0523         for key in binDict:
0524             binMap[key] = binPipes_parent[key].recv()
0525         print(
0526             datetime.now().strftime("%H:%M:%S")
0527             + "    Binning for job "
0528             + str(job)
0529             + " have been selected, now running the mapping",
0530             flush=True,
0531         )
0532         # Launch the material mapping with the bin map
0533         OptiJob[job] = Process(
0534             target=runMaterialMappingVariance,
0535             args=(
0536                 binMap,
0537                 args.topNumberOfEvents,
0538                 job,
0539                 args.inputPath,
0540                 pathExp,
0541                 resultPipes_child[job],
0542                 args.readCachedSurfaceInformation,
0543             ),
0544         )
0545         OptiJob[job].start()
0546 
0547     # Collect the score from the material mapping, this pauses the script until all the jobs have been completed
0548     for job in range(args.numberOfJobs):
0549         scores = resultPipes_parent[job].recv()
0550         print(
0551             datetime.now().strftime("%H:%M:%S")
0552             + "    Retried score for job "
0553             + str(job),
0554             flush=True,
0555         )
0556         for key in binDict:
0557             score = scores[key]
0558             scorePipes_parent[key].send(score)
0559         print(
0560             datetime.now().strftime("%H:%M:%S")
0561             + "    Job number "
0562             + str(job)
0563             + " is over",
0564             flush=True,
0565         )
0566 
0567     if args.doPloting:
0568         # The optimal binning has been found.
0569         # Run the material mapping one last to obtain a usable material map
0570         print(
0571             datetime.now().strftime("%H:%M:%S")
0572             + "    Running the material mapping to obtain the optimised material map",
0573             flush=True,
0574         )
0575         resultBinMap = dict()
0576         for key in binDict:
0577             resultBinMap[key] = binPipes_parent[key].recv()
0578         matMapDeco.setBinningMap(resultBinMap)
0579 
0580         # Decorate the detector with the MappingMaterialDecorator
0581         resultDetector, resultTrackingGeometry, resultDecorators = getOpenDataDetector(
0582             matMapDeco
0583         )
0584 
0585         # Sequence for the mapping, only use one thread when mapping material
0586         rMap = acts.examples.Sequencer(
0587             events=args.topNumberOfEvents, numThreads=1, logLevel=acts.logging.INFO
0588         )
0589 
0590         # Run the material mapping
0591         from material_mapping import runMaterialMapping
0592 
0593         runMaterialMapping(
0594             resultTrackingGeometry,
0595             resultDecorators,
0596             outputDir=args.outputPath,
0597             inputDir=args.inputPath,
0598             mapName="optimised-material-map",
0599             format=JsonFormat.Json,
0600             mapVolume=False,
0601             s=rMap,
0602         )
0603         rMap.run()
0604         del rMap  # Need to be deleted to write the material map to cbor
0605         del resultDetector, resultTrackingGeometry, resultDecorators
0606     print(
0607         datetime.now().strftime("%H:%M:%S")
0608         + "    Waiting for all the score to have been stored",
0609         flush=True,
0610     )
0611 
0612     # Create a timer that will be used to terminate the Process
0613     deadline = time.time() + 600
0614 
0615     for key in binDict:
0616         timeout = max(deadline - time.time(), 0)
0617         expJob[key].join(timeout=timeout)
0618         expJob[key].terminate()
0619         print(
0620             datetime.now().strftime("%H:%M:%S")
0621             + "    Experiment for surface "
0622             + str(key)
0623             + " is over",
0624             flush=True,
0625         )