Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:07

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