Scripts/Python/misaligned_simulation.py

0001 #!/usr/bin/env python3
0002 """Simulate a detector that is placed differently from what reconstruction assumes.
0003
0004 The `AlignmentDecorator` writes its alignment into the simulation geometry context
0005 only, so Fatras and the digitization see a shifted layer while seeding, the CKF and
0006 the fit stay on the nominal geometry. The shift shows up as a bias in the track
0007 state residuals of that layer.
0008 """
0009
0010 import argparse
0011 import os
0012 from pathlib import Path
0013
0014 import acts
0015 import acts.examples
0016 from acts import UnitConstants as u
0017 from acts.examples import (
0018     GaussianVertexGenerator,
0019     RandomNumbers,
0020     Sequencer,
0021     StructureSelector,
0022     TelescopeDetector,
0023 )
0024 from acts.examples.alignment import (
0025     AlignmentDecorator,
0026     AlignmentGeneratorGlobalShift,
0027     GeoIdAlignmentStore,
0028 )
0029 from acts.examples.reconstruction import (
0030     CkfConfig,
0031     SeedFinderConfigArg,
0032     SeedFinderOptionsArg,
0033     SeedingAlgorithm,
0034     TrackSelectorConfig,
0035     addCKFTracks,
0036     addSeeding,
0037 )
0038 from acts.examples.simulation import (
0039     EtaConfig,
0040     MomentumConfig,
0041     ParticleConfig,
0042     ParticleSelectorConfig,
0043     PhiConfig,
0044     addDigiParticleSelection,
0045     addDigitization,
0046     addFatras,
0047     addParticleGun,
0048 )
0049
0050 # The telescope layout of millepede_alignment.py: square sensors in the x-z plane,
0051 # geo id layer running 2, 4, ..., 18 for the nine layers, all in volume 1.
0052 LAYER_POSITIONS = [30, 60, 90, 120, 150, 180, 210, 240, 270]
0053
0054
0055 def addMisalignment(
0056     s: Sequencer,
0057     trackingGeometry: acts.TrackingGeometry,
0058     layer: int,
0059     shift: acts.Vector3,
0060     target: AlignmentDecorator.Target,
0061     logLevel: acts.logging.Level = acts.logging.WARNING,
0062 ) -> AlignmentDecorator:
0063     """Shift one layer, and decorate the selected geometry context(s) with it."""
0064     geoId = acts.GeometryIdentifier(volume=1, layer=layer, sensitive=1)
0065
0066     generator = AlignmentGeneratorGlobalShift()
0067     generator.shift = shift
0068
0069     cfg = AlignmentDecorator.Config()
0070     cfg.target = target
0071     cfg.nominalStore = GeoIdAlignmentStore(
0072         StructureSelector(trackingGeometry).selectedTransforms(
0073             acts.GeometryContext.dangerouslyDefaultConstruct(), geoId
0074         )
0075     )
0076     # A single IOV covering the whole run, no time dependence
0077     cfg.iovGenerators = [((0, 10000000), generator)]
0078
0079     decorator = AlignmentDecorator(cfg, logLevel)
0080     s.addContextDecorator(decorator)
0081     return decorator
0082
0083
0084 def addTelescopeChain(
0085     s: Sequencer,
0086     trackingGeometry: acts.TrackingGeometry,
0087     field: acts.MagneticFieldProvider,
0088     outputDir: Path,
0089     rnd: RandomNumbers,
0090 ) -> None:
0091     srcdir = Path(__file__).resolve().parent.parent.parent.parent
0092
0093     addParticleGun(
0094         s,
0095         MomentumConfig(10 * u.GeV, 100 * u.GeV, transverse=True),
0096         EtaConfig(-0.3, 0.3),
0097         # the telescope points along +y
0098         PhiConfig(60 * u.degree, 120 * u.degree),
0099         ParticleConfig(1, acts.PdgParticle.eMuon, randomizeCharge=True),
0100         vtxGen=GaussianVertexGenerator(
0101             mean=acts.Vector4(0, 0, 0, 0),
0102             stddev=acts.Vector4(5.0 * u.mm, 0.0 * u.mm, 5.0 * u.mm, 0.0 * u.ns),
0103         ),
0104         multiplicity=1,
0105         rnd=rnd,
0106     )
0107     addFatras(
0108         s,
0109         trackingGeometry,
0110         field,
0111         enableInteractions=True,
0112         rnd=rnd,
0113     )
0114     addDigitization(
0115         s,
0116         trackingGeometry,
0117         field,
0118         digiConfigFile=srcdir / "Examples/Configs/telescope-digi-smearing-config.json",
0119         rnd=rnd,
0120     )
0121     addDigiParticleSelection(
0122         s,
0123         ParticleSelectorConfig(measurements=(3, None), removeNeutral=True),
0124     )
0125     addSeeding(
0126         s,
0127         trackingGeometry,
0128         field,
0129         seedFinderConfigArg=SeedFinderConfigArg(
0130             r=(20 * u.mm, 200 * u.mm),
0131             deltaR=(1 * u.mm, 300 * u.mm),
0132             collisionRegion=(-250 * u.mm, 250 * u.mm),
0133             z=(-100 * u.mm, 100 * u.mm),
0134             maxSeedsPerSpM=1,
0135             sigmaScattering=5,
0136             radLengthPerSeed=0.1,
0137             minPt=0.5 * u.GeV,
0138             impactMax=3 * u.mm,
0139         ),
0140         seedFinderOptionsArg=SeedFinderOptionsArg(bFieldInZ=2 * u.T),
0141         seedingAlgorithm=SeedingAlgorithm.GridTriplet,
0142         initialSigmas=[
0143             3 * u.mm,
0144             3 * u.mm,
0145             1 * u.degree,
0146             1 * u.degree,
0147             0 * u.e / u.GeV,
0148             1 * u.ns,
0149         ],
0150         initialSigmaQoverPt=0.1 * u.e / u.GeV,
0151         initialSigmaPtRel=0.1,
0152         initialVarInflation=[1.0] * 6,
0153         geoSelectionConfigFile=srcdir
0154         / "Examples/Configs/telescope-seeding-config.json",
0155     )
0156     addCKFTracks(
0157         s,
0158         trackingGeometry,
0159         field,
0160         TrackSelectorConfig(),
0161         CkfConfig(
0162             chi2CutOffMeasurement=150.0,
0163             chi2CutOffOutlier=250.0,
0164             numMeasurementsCutOff=50,
0165             seedDeduplication=True,
0166             stayOnSeed=True,
0167         ),
0168         outputDirRoot=outputDir,
0169         writeTrackStates=True,
0170         writePerformance=False,
0171     )
0172
0173
0174 def main():
0175     parser = argparse.ArgumentParser(description=__doc__)
0176     parser.add_argument(
0177         "--output", "-o", type=Path, default=Path.cwd() / "misaligned_output"
0178     )
0179     parser.add_argument("--events", "-n", type=int, default=1000)
0180     parser.add_argument(
0181         "--layer",
0182         type=int,
0183         default=4,
0184         help="geo id layer to misalign (2, 4, ..., 18)",
0185     )
0186     parser.add_argument(
0187         "--shift",
0188         type=float,
0189         default=0.2,
0190         help="global z shift of that layer in mm",
0191     )
0192     parser.add_argument(
0193         "--target",
0194         choices=["sim", "reco", "both"],
0195         default="sim",
0196         help=(
0197             "which geometry context sees the shift. 'sim' and 'reco' put the "
0198             "detector and the reconstruction hypothesis out of step, 'both' keeps "
0199             "the shift perfectly known and residuals centred."
0200         ),
0201     )
0202     args = parser.parse_args()
0203
0204     os.makedirs(args.output, exist_ok=True)
0205
0206     detector = TelescopeDetector(
0207         bounds=[200, 200],
0208         positions=LAYER_POSITIONS,
0209         stereos=[0] * len(LAYER_POSITIONS),
0210         binValue=1,
0211     )
0212     trackingGeometry = detector.trackingGeometry()
0213     field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
0214
0215     s = Sequencer(events=args.events, numThreads=1, outputDir=str(args.output))
0216
0217     addMisalignment(
0218         s,
0219         trackingGeometry,
0220         layer=args.layer,
0221         shift=acts.Vector3(0, 0, args.shift * u.mm),
0222         target={
0223             "sim": AlignmentDecorator.Target.eSim,
0224             "reco": AlignmentDecorator.Target.eReco,
0225             "both": AlignmentDecorator.Target.eBoth,
0226         }[args.target],
0227     )
0228
0229     addTelescopeChain(s, trackingGeometry, field, args.output, RandomNumbers(seed=42))
0230
0231     s.run()
0232
0233     print(
0234         f"\nWrote {args.output / 'trackstates_ckf.root'}. The residual "
0235         f"res_eLOC0/res_eLOC1 for volume 1 layer {args.layer} carries the "
0236         f"{args.shift} mm shift unless --target both."
0237     )
0238
0239
0240 if __name__ == "__main__":
0241     main()