diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 9f62fd57b12fa450f6e2218b841501ca2351cf2a..c3f5313c0e52ddc68009ac2d6d3e41103b103b8a 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -6,6 +6,7 @@
#include <map>
#include <utility>
+#include <cmath>
CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name,
ISvcLocator* pSvcLocator)
@@ -99,19 +100,31 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
// Create structure to store pulse for each channel
std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
+ // Sum energy for each channel (i.e. identifier)
+ std::map<unsigned int, float> esum;
+ for (const auto& hit : *caloHitHandle) {
+ esum[hit.identify()] += hit.energyLoss();
+ }
+
+ // Loop over time samples
for (const auto& tk : m_timekernel) {
std::map<unsigned int, float> counts;
// Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
- for (const auto& hit : *caloHitHandle) {
- counts[hit.identify()] += tk * hit.energyLoss();
+ //for (const auto& hit : *caloHitHandle) {
+ // counts[hit.identify()] += tk * hit.energyLoss();
+ //}
+
+ // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+ for (const auto& e : esum) {
+ counts[e.first] = tk * e.second;
}
// Subtract count from basleine and add result to correct waveform vector
for (const auto& c : counts) {
- unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = baseline - c.second;
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = std::round(baseline - c.second);
if (value < 0) {
ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
diff --git a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
index 83839b9cba1a4c6a50eef52852ed8ee2127dfaa7..c49770ad2b7e9468eca196af582d2156b9ab977c 100644
--- a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
+++ b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
@@ -19,12 +19,11 @@ def CalorimeterReconstructionCfg(flags, **kwargs):
kwargs.setdefault("CaloHitContainerKey", "CaloHits")
kwargs.setdefault("PreshowerHitContainerKey", "PreshowerHits")
- recoAlg = CompFactory.CaloRecAlg("CaloRecAlg", **kwargs)
- acc.addEventAlgo(recoAlg)
-
- dbInstance = kwargs.get("dbInstance", "TRIGGER_OFL") # what should this be set to???
+ acc.merge(CaloRecToolCfg(flags, **kwargs))
+ kwargs.pop("MC_calibTag") # Remove this if it is specified so it does not get pased to CaloRecAlg
- acc.merge(CaloRecToolCfg(flags))
+ recoAlg = CompFactory.CaloRecAlg("CaloRecAlg", isMC = flags.Input.isMC, **kwargs)
+ acc.addEventAlgo(recoAlg)
return acc
diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
index f4f78d488d6c542637003e38ea0d8959f89b6b1c..b6df5fed6053f7499a6f45a95e4b51d2e6633ddd 100644
--- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
+++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
@@ -91,12 +91,11 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("calo_hit filled has charge of " << charge << " pC");
float gainRatio = 1.0;
- if (m_isMC) {
- gainRatio = 1.0; // put dummy value for now, this will end up being ratio of digi scale factors
- } else {
+ if (!m_isMC) { // MC already has correct MIP charge stored in MIPcharge_ref, so only need to to HV extrapolation with reral data
gainRatio = extrapolateHVgain(hit->channel());
- ATH_MSG_DEBUG("HV gain ratio = " << gainRatio );
}
+ ATH_MSG_DEBUG("HV gain ratio = " << gainRatio );
+
float Nmip = (charge * gainRatio) / MIPcharge_ref;
ATH_MSG_DEBUG("Nmip = " << Nmip );
diff --git a/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py b/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py
index 65be875da12b6deee3cbf9f35217cee6f40784ff..473a70f7327a39b4029b0af40045c84b1054cc35 100644
--- a/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py
+++ b/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py
@@ -19,8 +19,15 @@ def CaloRecToolCfg(flags, name="CaloRecTool", **kwargs):
# Probably need to figure this out!
dbInstance = kwargs.get("dbInstance", "TRIGGER_OFL")
+ MC_calibTag = kwargs.get("MC_calibTag", "")
+
+ # MC calibration db folder MIP_ref has version tags
+ if len(MC_calibTag) > 0:
+ acc.merge(addFolders(flags, "/WAVE/Calibration/MIP_ref", dbInstance, className="CondAttrListCollection", tag=MC_calibTag))
+ else:
+ acc.merge(addFolders(flags, "/WAVE/Calibration/MIP_ref", dbInstance, className="CondAttrListCollection"))
+
acc.merge(addFolders(flags, "/WAVE/Calibration/HV", dbInstance, className="CondAttrListCollection"))
- acc.merge(addFolders(flags, "/WAVE/Calibration/MIP_ref", dbInstance, className="CondAttrListCollection"))
return acc
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
index 03368e676c1e73caca7034aedc5fdcff20ed3e63..b945ffbbd19c024d549d0d0f84fdfca768759f68 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
@@ -29,6 +29,8 @@ parser.add_argument("-s", "--slice", type=int, default=0,
help="Specify file slice to produce")
parser.add_argument("-f", "--files", type=int, default=5,
help="Specify number of input files to run in one batch")
+parser.add_argument("--complete", action="store_true",
+ help="This is the complete run, remove segment numbers")
parser.add_argument("-t", "--tag", default="",
help="Specify digi tag (to append to output filename)")
parser.add_argument("--subtractTime", type=float, default=-999.,
@@ -69,24 +71,35 @@ if len(dirlist) == 0:
print("HITS files available:")
[print(file) for file in dirlist]
+# Get run number for pattern matching
+stem = dirlist[0].stem
+spl = stem.split('-')
+short = spl[1]
+run = spl[2]
+
for seg in seglist:
- # Assume these are in numerical order from 0
- if seg >= len(dirlist):
- print(f"Requested file segment {seg} but only {len(dirlist)} files found")
- if args.partial:
- break
- else:
- sys.exit(1) # Abort this job
+ # No longer assume all HITS files are complete
+ #if seg >= len(dirlist):
+ # print(f"Requested file segment {seg} but only {len(dirlist)} files found")
+ # if args.partial:
+ # break
+ # else:
+ # sys.exit(1) # Abort this job
# Check if segment number exists in hits file (this is not perfect)
- segstr = f"{seg:05d}"
- if segstr not in dirlist[seg].name:
- print(f"Segment {segstr} not in file {dirlist[seg].name}!")
+ segstr = f"FaserMC-{short}-{run}-{seg:05d}*HITS.root"
+ matchlist = list(dirpath.glob(segstr))
+ if len(matchlist) == 0:
+ print(f"Segment {segstr} not found!")
if not args.partial: sys.exit(1) # abort
+ continue
+ elif len(matchlist) > 1:
+ print(f"Segment {segstr} matched {len(matchlist)} files!")
+ sys.exit(1) # Should never happen!
else:
- print(f"Segment {segstr} found in file {dirlist[seg]}")
+ print(f"Segment {segstr} found!")
- filelist.append(dirlist[seg])
+ filelist.append(matchlist[0])
if len(filelist) == 0:
# Asked for range that doesn't exist
@@ -116,8 +129,11 @@ seghi = int(spl[3])
if len(args.short) > 0:
short += args.short
+# Shouldn't normally remove segments with MC as we may always add
+# more files later. --complete overrides this if wanted
+
# Build output filename
-if seglo == 0 and (seghi+1) == len(dirlist): # Full run
+if args.complete and (seglo == 0) and ((seghi+1) == len(dirlist)): # Full run
outfile = f"FaserMC-{short}-{run}"
elif seglo == seghi: # Single segment
outfile = f"FaserMC-{short}-{run}-{seglo:05}"
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
index 34e0275042d59ec162c44e8198925e9f6d098853..fcf5e2045b192ef913a9fc1ea0b75c8fae4074bd 100644
--- a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
@@ -2,10 +2,10 @@
# Used with a condor file to submit to vanilla universe
#
# Usage:
-# submit_faser_digi.sh [--highGain] filepath [release_directory] [working_directory]
+# submit_faser_digi.sh filepath [release_directory] [working_directory]
#
# Options:
-# --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
+# --digiTag <tag> - override digitization tag for calo gain
# --geom - geometry setting
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
@@ -186,6 +186,9 @@ cd "$file_stem"
# Run job
#
faser_digi.py $geomstr $gainstr $timestr $tagstr "$file_path"
+digi_code=$?
+echo "Return code: $digi_code"
+
#
# Print out ending time
date
@@ -196,25 +199,52 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-RDO.root to $outdest"
- mkdir -p $outdest
- eos cp *-RDO.root ${outdest}/ || true
+ thefile=`ls *-RDO.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-RDO.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
- mkdir -p $logdest
+ eos mkdir -p $logdest
eos cp $logfile $logdest/$logfile
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
- mkdir -p $outdest
+ eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $digi_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
index cb67ef8838ee26a01e15436bd44668e672b20379..50278c33e2435236fff3bc2b6a94c13aa5794817 100755
--- a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
@@ -2,10 +2,11 @@
# Used with a condor file to submit to vanilla universe
#
# Usage:
-# submit_faser_digi_merge.sh [--highGain] dirpath slice nfiles [release_directory] [working_directory]
+# submit_faser_digi_merge.sh dirpath slice nfiles [release_directory] [working_directory]
#
# Options:
-# --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
+# --digiTag <tag> - override digitization tag for calo gain
+# --partial - allow missing files in merge
# --geom - geometry setting
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
@@ -221,6 +222,8 @@ cd "$file_stem"
# Run job
#
faser_digi_merge.py $partialstr $geomstr $gainstr $timestr $tagstr --slice $slice --files $nfiles $dir_path
+digi_code=$?
+echo "Return code: $digi_code"
#
# Print out ending time
date
@@ -231,25 +234,52 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-RDO.root to $outdest"
- mkdir -p $outdest
- eos cp *-RDO.root ${outdest}/ || true
+ thefile=`ls *-RDO.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-RDO.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
- mkdir -p $logdest
+ eos mkdir -p $logdest
eos cp $logfile $logdest/$logfile
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
- mkdir -p $outdest
+ eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $digi_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 43f98387e7089cfdf258531ec57aeea072b70f53..055dd7200201ce69823d798609cb991d85f4e794 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -4,6 +4,13 @@
#
# Parser function for particle gun samples
#
+
+def float_or_none(arg):
+ if arg is None or arg == "None":
+ return None
+
+ return float(arg)
+
def faser_pgparser():
import sys
@@ -33,15 +40,24 @@ def faser_pgparser():
help="Specify PDG ID of particle (note plus/minus different) or list (e.g.: --pid -13 13)")
parser.add_argument("--mass", default=105.66, type=float,
help="Specify particle mass (in MeV)")
- parser.add_argument("--radius", default=100., type=float,
+ parser.add_argument("--radius", default=100., type=float, nargs="*",
help="Specify radius (in mm)")
- parser.add_argument("--angle", default=0.005, type=float,
+ parser.add_argument("--angle", default=0.005, type=float_or_none,
help="Specify angular width (in Rad)")
parser.add_argument("--zpos", default=None, type=float,
help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
+ parser.add_argument("--pidd1", default=None, type=int,
+ help="Specify PDG ID of daugther 1 for DIF generator")
+ parser.add_argument("--pidd2", default=None, type=int,
+ help="Specify PDG ID of daugther 2 for DIF generator")
+
+
parser.add_argument("--sampler", default="log",
- help="Specify energy sampling (log, lin, const)")
+ help="Specify energy sampling (log, lin, const, hist, hist2D)")
+ parser.add_argument("--hist_name", default="",
+ help="Specify energy sampling name for hist sampler file.root:hist")
+
parser.add_argument("--minE", default=10., type=float,
help="Minimum energy in GeV (for log or lin sampler)")
parser.add_argument("--maxE", default=1000., type=float,
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index cc1628efdb9245e32b7abbbc9ca60a60b2151fad..169badecf7cf2a392a95b098a858634c5a2523c6 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -101,24 +101,50 @@ if __name__ == '__main__':
from AthenaCommon.PhysicalConstants import pi
if isinstance(args.pid, list):
- # Note args.pid is a list, must make this a set for ParticleGun
- pidarg = set(args.pid)
+ # Note args.pid is a list, must make this a set for ParticleGun if have > 1
+ if len(args.pid) > 1:
+ pidarg = set(args.pid)
+ else:
+ pidarg = args.pid[0]
else:
# Just pass a single value
pidarg = args.pid
print(f"Using pid: {args.pid} => {pidarg}")
+ if isinstance(args.radius, list) and len(args.radius) == 1:
+ args.radius = args.radius[0]
+
+
# Create the simgun dictionary
# Negative radius gives uniform sampling
# Positive radius gives Gaussian sampling
- sg_dict = {
- "Generator" : "SingleParticle",
- "pid" : pidarg, "mass" : args.mass,
- "theta" : PG.GaussianSampler(0, args.angle, oneside = True),
- "phi" : [0, 2*pi], "radius" : args.radius,
- "randomSeed" : args.outfile
- }
+
+ if args.pidd1:
+ if args.pidd2 is None:
+ args.pidd2 = -args.pidd1
+
+ sg_dict = {
+ "Generator" : "DecayInFlight",
+ "mother_pid" : pidarg,
+ "daughter1_pid" : args.pidd1,
+ "daughter2_pid" : args.pidd2,
+ "mass" : args.mass,
+ "theta" : PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
+ "phi" : [0, 2*pi], "radius" : args.radius,
+ "randomSeed" : args.outfile
+ }
+
+ else:
+
+ sg_dict = {
+ "Generator" : "SingleParticle",
+ "pid" : pidarg, "mass" : args.mass,
+ "theta" : PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
+ "phi" : [0, 2*pi], "radius" : args.radius,
+ "randomSeed" : args.outfile
+ }
+
# -1000 is safely upstream of detector (to be checked)
# Note zpos is in mm!
@@ -132,6 +158,12 @@ if __name__ == '__main__':
sg_dict["energy"] = PG.LogSampler(args.minE*GeV, args.maxE*GeV)
elif args.sampler == "const":
sg_dict["energy"] = PG.ConstSampler(args.maxE*GeV)
+ elif args.sampler == "hist":
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+ elif args.sampler == "hist2D":
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH2Sampler(fname, hname)
else:
print(f"Sampler {args.sampler} not known!")
sys.exit(1)
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index b470cb3fd20ae482d4cd1aececd535191a35c00e..21e1d51a6af7d92ebac09d051e80c919d1c39687 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -33,6 +33,8 @@ parser.add_argument("-v", "--verbose", action='store_true',
help="Turn on DEBUG output")
parser.add_argument("--isMC", action='store_true',
help="Running on digitised MC rather than data")
+parser.add_argument("--MC_calibTag", default="",
+ help="Specify tag used to reconstruct MC calo energy: (WAVE-Calibration-01-LG-nofilt, WAVE-Calibration-01-LG, WAVE-Calibration-01-HG-nofilt, or WAVE-Calibration-01-HG) ")
parser.add_argument("--testBeam", action='store_true',
help="Set geometry for 2021 test beam")
@@ -129,8 +131,6 @@ elif runtype == "TI12Data02":
# Final 2022 TI12 geometry
elif runtype == "TI12Data03":
ConfigFlags.GeoModel.FaserVersion = "FASERNU-03"
- # ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"
- # Use the updated field map
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"
useCal = True
if not args.isMC:
@@ -141,7 +141,6 @@ else:
print("Specify correct type or update list")
sys.exit(1)
-
# Must use original input string here, as pathlib mangles double // in path names
ConfigFlags.Input.Files = [ args.file_path ]
@@ -194,14 +193,10 @@ if useLHC:
from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
acc.merge(WaveformReconstructionCfg(ConfigFlags))
-# Calorimeter reconstruction
-if args.isMC:
- # Not ready for MC quite yet
- pass
-
-elif useCal:
+# Calorimeter Energy reconstruction
+if useCal:
from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
- acc.merge(CalorimeterReconstructionCfg(ConfigFlags))
+ acc.merge(CalorimeterReconstructionCfg(ConfigFlags, MC_calibTag=args.MC_calibTag))
# Tracker clusters
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
@@ -295,6 +290,12 @@ if not args.isMC:
if args.verbose:
acc.foreach_component("*").OutputLevel = VERBOSE
ConfigFlags.dump()
+ # Print out POOL conditions
+ import os
+ print(f"ATLAS_POOLCOND_PATH: {os.environ['ATLAS_POOLCOND_PATH']}")
+ print(f"PoolSvc.ReadCatalog: {acc.getService('PoolSvc').ReadCatalog}")
+ print(f"PoolSvc.WriteCatalog: {acc.getService('PoolSvc').WriteCatalog}")
+
else:
acc.foreach_component("*").OutputLevel = INFO
@@ -314,5 +315,11 @@ b = time.time()
from AthenaCommon.Logging import log
log.info(f"Finish execution in {b-a} seconds")
-sys.exit(int(sc.isFailure()))
+# Signal errors
+if sc.isSuccess():
+ log.info("Execution succeeded")
+ sys.exit(0)
+else:
+ log.info("Execution failed, return 1")
+ sys.exit(1)
diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
index 4fe3e2d626796ab0ba52fddee0aee5e311cf59ed..f75b144a1bf2d0a725692dbaa82be34ad788cd5c 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
+++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
@@ -8,7 +8,10 @@
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
# --geom - specify geometry
+#
+# Monte Carlo options:
# --isMC - needed for MC reco
+# --caloTag <tag> - override MC reco tag for calo gain (to match digi tag)
#
# file_path - full file name (with path)
# release_directory - optional path to release install directory (default pwd)
@@ -25,6 +28,9 @@
# Keep track of time
SECONDS=0
#
+# Job option strings
+gainstr=""
+#
# Parse command-line options
while [ -n "$1" ]
do
@@ -45,9 +51,16 @@ do
shift;;
--isMC)
+ echo "Set isMC true"
ismc=1
shift;;
+ --caloTag)
+ echo "Override calo digi tag with $2"
+ gainstr="--MC_calibTag $2"
+ shift;
+ shift;;
+
--) # End of options
shift; # Eat this
break;; # And stop parsing
@@ -119,6 +132,7 @@ echo `date` - $HOSTNAME
echo "File: $file_name"
echo "Filepath: $file_path"
echo "Geom: $geom"
+echo "Gain: $gainstr"
echo "Release: $release_directory"
echo "Output: $output_directory"
echo "Starting: $starting_directory"
@@ -127,7 +141,6 @@ echo "Starting: $starting_directory"
export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase
source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh
#
-# Try automatic
# Always go back to the starting directory in case paths are relative
cd "$starting_directory"
cd "$release_directory"
@@ -140,6 +153,11 @@ source run/setup.sh
#
echo "ATLAS_POOLCOND_PATH = $ATLAS_POOLCOND_PATH"
#
+# Check if there are data overrides
+if [ -d "run/data/sqlite200" ]; then
+ cond_directory=`pwd -P`/run/data # Get absolute path
+fi
+#
# Try to find a release tag
cd calypso
recotag=`git describe --tags`
@@ -170,6 +188,29 @@ else
fi
cd "$file_stem"
#
+# Check if there were data overrides in the release directory
+if [[ -z "$cond_directory" ]]; then
+ echo "No local conditions directory found!"
+else
+ echo "Local conditions directory found! Copying to run directory..."
+ echo Copying $cond_directory
+ cp -r $cond_directory .
+ ls -R data
+fi
+#
+# Further check if there is a pool conditions override
+#if [[ -d "data/poolcond" ]]; then
+# echo "Local POOL directory found!"
+# echo "Change ATLAS_POOLCOND_PATH"
+# echo " from $ATLAS_POOLCOND_PATH"
+# export ATLAS_POOLCOND_PATH=`pwd -P`/data
+# echo " to $ATLAS_POOLCOND_PATH"
+#else
+# echo "No local pool files found, use default:"
+# echo " $ATLAS_POOLCOND_PATH"
+#fi
+echo "ATLAS_POOLCOND_PATH: $ATLAS_POOLCOND_PATH"
+#
# Run job
if [[ -z "$tag" ]]; then
tagstr=""
@@ -189,7 +230,9 @@ else
mcstr="--isMC"
fi
#
-faser_reco.py "--nevents=$nevents" $geomstr $tagstr $mcstr "$file_path"
+faser_reco.py "--nevents=$nevents" $geomstr $tagstr $mcstr $gainstr "$file_path"
+reco_code=$?
+echo "Return code: $reco_code"
#
# Print out ending time
date
@@ -203,11 +246,27 @@ if ! [ -z "$outdest" ]
then
echo "Output directory:"
ls -l
- echo "copy *-xAOD.root to $outdest"
- eos mkdir -p $outdest
- # Keep this line from stopping script, so we might get a log file
- # || true ensures script continues even if copy fails
- eos cp *-xAOD.root ${outdest}/ || true
+ thefile=`ls *-xAOD.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ # Keep this line from stopping script, so we might get a log file
+ # || true ensures script continues even if copy fails
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that it worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-xAOD.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Copy log file second
@@ -228,3 +287,11 @@ then
eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $reco_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh b/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
index 83a2470f1e93bd9703e0e44c5abca3a9aeebe7c1..6e4d4df2eac92d4836c2b77104093ab2dd81c01c 100755
--- a/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
+++ b/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
@@ -35,6 +35,11 @@ do
xangle=1
shift;; # This 'eats' the argument
+ --fshift)
+ echo "Applying FLUKA crossing-angle shift"
+ xangle=2
+ shift;;
+
-l | --log)
logdest="$2";
shift;
@@ -183,9 +188,16 @@ cd "${file_stem}"
#fi
if [[ -z "$xangle" ]]; then
faser_simulate.py --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
+elif [[ $xangle == 1 ]]; then
+ faser_simulate.py --yangle -0.000160 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
else
- faser_simulate.py --yangle -0.000150 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ echo "Applying double crossing-angle shift for FLUKA!"
+ faser_simulate.py --yangle -0.000310 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
fi
+echo "Return code: $sim_code"
#
# Print out ending time
date
@@ -196,16 +208,34 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-HITS.root to $outdest"
- eos mkdir -p $outdest
- eos cp *-HITS.root ${outdest}/ || true
+ thefile=`ls *-HITS.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-HITS.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
eos mkdir -p $logdest
@@ -213,9 +243,17 @@ then
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+# Make sure to return an error is calypso failed
+if [ $sim_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser02.py b/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser02.py
index 2e944e3fc88b76d6794a7df7b9ce4469fc763d51..56b856fa824564a1368675eb45bb1d6166a7949d 100644
--- a/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser02.py
+++ b/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser02.py
@@ -53,11 +53,11 @@ if __name__ == "__main__":
# Flags for this job
ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
- ConfigFlags.GeoModel.FaserVersion = "FASER-02" # Default FASER geometry
- ConfigFlags.IOVDb.GlobalTag = "OFLCOND-"+ ConfigFlags.GeoModel.FaserVersion # Always needed; must match FaserVersion
- ConfigFlags.IOVDb.DBConnection = "sqlite://;schema=" + ConfigFlags.GeoModel.FaserVersion + "_ALLP200.db;dbname=OFLP200"
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Default FASER geometry
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Old field map
+ ConfigFlags.IOVDb.DBConnection = "sqlite://;schema=FASER-02_ALLP200.db;dbname=OFLP200"
ConfigFlags.GeoModel.Align.Disable = True # Hack to avoid loading alignment when we want to create it from scratch
- ConfigFlags.addFlag("WriteAlignment.PoolFileName", ConfigFlags.GeoModel.FaserVersion + "_Align.pool.root")
+ ConfigFlags.addFlag("WriteAlignment.PoolFileName", "FASER-02_Align.pool.root")
# Parse flags from command line and lock
ConfigFlags.addFlag("AlignDbTool.AlignmentConstants", {})
diff --git a/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser03.py b/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser03.py
new file mode 100755
index 0000000000000000000000000000000000000000..2bbddc7cde737cfd5506c0c4bc6610ca820fb55b
--- /dev/null
+++ b/Control/CalypsoExample/WriteAlignment/python/WriteAlignmentConfig_Faser03.py
@@ -0,0 +1,82 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+#!/usr/bin/env python
+import sys
+from AthenaCommon.Constants import VERBOSE, INFO
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+def WriteAlignmentCfg(flags, name="WriteAlignmentAlg", alignmentConstants={}, **kwargs):
+
+ # Initialize GeoModel
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ a = FaserGeometryCfg(flags)
+
+ # This section is to allow alignment to be written to a conditions DB file
+ from IOVDbSvc.IOVDbSvcConfig import IOVDbSvcCfg
+ result = IOVDbSvcCfg(flags)
+ iovDbSvc = result.getPrimary()
+ iovDbSvc.dbConnection=flags.IOVDb.DBConnection
+ a.merge(result)
+
+ AthenaPoolCnvSvc=CompFactory.AthenaPoolCnvSvc
+ apcs=AthenaPoolCnvSvc()
+ a.addService(apcs)
+ EvtPersistencySvc=CompFactory.EvtPersistencySvc
+ a.addService(EvtPersistencySvc("EventPersistencySvc",CnvServices=[apcs.getFullJobOptName(),]))
+
+ a.addService(CompFactory.IOVRegistrationSvc(PayloadTable=False,OutputLevel=VERBOSE))
+
+ # Configure the algorithm itself
+ WriteAlignmentAlg = CompFactory.WriteAlignmentAlg
+ outputTool = CompFactory.AthenaOutputStreamTool("DbStreamTool", OutputFile = flags.WriteAlignment.PoolFileName,
+ PoolContainerPrefix="ConditionsContainer",
+ TopLevelContainerName = "<type>",
+ SubLevelBranchName= "<key>" )
+
+ trackerAlignDBTool = CompFactory.TrackerAlignDBTool("AlignDbTool", OutputTool = outputTool,
+ OutputLevel=VERBOSE,
+ AlignmentConstants = {})
+ kwargs.setdefault("AlignDbTool", trackerAlignDBTool)
+ trackerAlignDBTool.AlignmentConstants = alignmentConstants
+ a.addEventAlgo(WriteAlignmentAlg(name, **kwargs))
+
+ return a
+
+
+if __name__ == "__main__":
+ # from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Configurable import Configurable
+ # from AthenaConfiguration.ComponentFactory import CompFactory
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+ Configurable.configurableRun3Behavior = True
+
+# Flags for this job
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Default FASER geometry
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"
+ ConfigFlags.IOVDb.DBConnection = "sqlite://;schema=FASER-03_ALLP200.db;dbname=OFLP200"
+ ConfigFlags.GeoModel.Align.Disable = True # Hack to avoid loading alignment when we want to create it from scratch
+ ConfigFlags.addFlag("WriteAlignment.PoolFileName", "FASER-03_Align.pool.root")
+
+# Parse flags from command line and lock
+ ConfigFlags.addFlag("AlignDbTool.AlignmentConstants", {})
+ ConfigFlags.fillFromArgs(sys.argv[1:])
+ ConfigFlags.lock()
+
+# Configure components
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ acc = MainServicesCfg(ConfigFlags)
+
+# Set things up to create a conditions DB with neutral Tracker alignment transforms
+ acc.merge(WriteAlignmentCfg(ConfigFlags, alignmentConstants=ConfigFlags.AlignDbTool.AlignmentConstants, ValidRunStart=1, ValidEvtStart=0, ValidRunEnd=9999999, ValidEvtEnd=9999999, CondTag=ConfigFlags.GeoModel.FaserVersion.replace("FASER", "TRACKER-ALIGN"), ))
+
+# Configure verbosity
+ # ConfigFlags.dump()
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ acc.foreach_component("*").OutputLevel = VERBOSE
+ acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # log.setLevel(VERBOSE)
+
+# Execute and finish
+ sys.exit(int(acc.run(maxEvents=1).isFailure()))
diff --git a/Generators/DIFGenerator/python/DIFSampler.py b/Generators/DIFGenerator/python/DIFSampler.py
index 8ec22459aba93a430093c4bd7b2a3f13e2d0f206..007830b63c82ec06f19af6fbd3169cfc10ce8e2e 100644
--- a/Generators/DIFGenerator/python/DIFSampler.py
+++ b/Generators/DIFGenerator/python/DIFSampler.py
@@ -93,9 +93,8 @@ class DIFSampler(PG.ParticleSampler):
# def __init__(self, daughter1_pid = 11, daughter2_pid = -11, mother_pid = None, mother_mom = PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (10*MeV)**2),0,10*MeV,0),mother_pos = CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0)):
def __init__(self, daughter1_pid=13, daughter2_pid=-13, mother_pid=None,
mother_mom=PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (500*MeV)**2), [0, 0.0002], 500*MeV, [0, 2*pi]),
- my_z_position=-1500):
- # mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
- self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=CylinderSampler([0, 100**2], [0, 2*pi], my_z_position, 0))
+ mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
+ self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=mother_pos)
self.daughter1 = self.particle(daughter1_pid)
self.daughter2 = self.particle(daughter2_pid)
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index ac30494a4815f5bb06a5d67408b4c0551297a63b..d2e82a88d3fd6229b2e3ad33932ead7d5d650720 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -9,7 +9,7 @@ import ParticleGun as PG
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
# from AthenaCommon.Constants import VERBOSE, INFO
-from AthenaCommon.SystemOfUnits import TeV
+from AthenaCommon.SystemOfUnits import TeV, GeV, MeV
from AthenaCommon.PhysicalConstants import pi
### add radial pos sampler ### with gaussian beam implemented
@@ -144,7 +144,31 @@ def FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs) :
pg = cfg.getPrimary()
from DIFGenerator import DIFSampler
- pg.sampler = DIFSampler(**kwargs)
+
+ if "radius" in kwargs:
+ kwargs["mother_pos"] = RadialPosSampler(x = kwargs.setdefault("x", 0.0),
+ y = kwargs.setdefault("y", 0.0),
+ z = kwargs.setdefault("z", -3750.0),
+ r = kwargs.setdefault("radius", 1.0),
+ t = kwargs.setdefault("t", 0.0) )
+ else:
+ from DIFGenerator.DIFSampler import CylinderSampler
+ kwargs["mother_pos"] = CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)
+
+ if not "mother_mom" in kwargs:
+ kwargs["mother_mom"] = PG.EThetaMPhiSampler(kwargs.setdefault("energy", ((1*TeV)**2 + (500*MeV)**2)**0.5),
+ kwargs.setdefault("theta", [0, 0.0002]),
+ kwargs.setdefault("mass", 500*MeV),
+ kwargs.setdefault("phi", [0, 2*pi]) )
+
+
+
+ pg.sampler = DIFSampler(kwargs.setdefault("daughter1_pid", 13),
+ kwargs.setdefault("daughter2_pid", -13),
+ kwargs.setdefault("mother_pid", None),
+ kwargs["mother_mom"],
+ kwargs["mother_pos"] )
+
return cfg
@@ -190,9 +214,12 @@ def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
return cfg
def FaserParticleGunCfg(ConfigFlags) :
- # generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
- generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
+ generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
+ # generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
+
kwargs = ConfigFlags.Sim.Gun
+ del kwargs["Generator"]
+
if generator == "SingleEcalParticle" :
return FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs)
elif generator == "Cosmics" :
diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index 6fcd6829b29ea474d5c7380526d3c46bab9c6871..19b1db116991c082d0771699813df2fe94c0a1d9 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -35,12 +35,17 @@ class RadialPosSampler(Sampler):
def r(self):
"r position sampler"
- fwhm = 2*self.radius
- sig = fwhm/(2 * sqrt(2 * log(2)))
-
- if self.radius < 0:
- return sqrt(random.uniform(0, abs(self.radius**2)))
+ #fwhm = 2*self.radius
+ #sig = fwhm/(2 * sqrt(2 * log(2)))
+
+ if isinstance(self.radius, list) and len(self.radius) == 2:
+ # If a list of length 2, generate uniformally over an annulus from r[0] to r[1]
+ return sqrt(random.uniform(abs(self.radius[0]**2), abs(self.radius[1]**2)))
+ elif self.radius < 0:
+ # Generate uniformally up to |r| if r is < 0
+ return sqrt(random.uniform(0, abs(self.radius**2)))
else:
+ # Else generate as a Gaussian with widht r
x = random.gauss(0, self.radius)
y = random.gauss(0, self.radius)
return sqrt(x**2 + y**2)
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 3d84bb15354e0bd1ab8bf47eb1a380b27d08939a..6ef5d1d0377a147246b97b9e203fcc224355bae2 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -47,13 +47,13 @@ class ForeseeGenerator(object):
else:
self.foresee = Foresee(path = self.path)
- # Generate 6 cm high to account for translation from ATLAS to FASER coord. system
+ # Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
# TODO: relax this a bit as daughters may enter even if mother doesn't
# self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1",
# channels=[self.mode], distance=480, length=1.5 ,
# luminosity=1/1000.) # 1 pb-1
- self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.06)**2)< 0.1",
+ self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.1",
channels=[self.mode], distance=480, length=1.5 ,
luminosity=1/1000.) # 1 pb-1
@@ -459,7 +459,7 @@ def production(args):
f.write()
if args.outdir:
- p = proc.run(["cp", filename, outdir + "/" + filename.split("/")[-1]])
+ p = proc.run(["cp", filename, outdir + "/" + filename.split("/")[-1].replace(".hepmc", f"{nrun}.hepmc")])
if not p.returncode:
proc.run(["rm", filename])
@@ -492,7 +492,7 @@ if __name__ == "__main__":
parser.add_argument("--couplings", "-c", required = False, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
parser.add_argument("--pid1", default = None, required = False, type = int, help = "PID of daughter 1")
parser.add_argument("--pid2", default = None, type = int, help = "PID of daughter 2 (if not set then will be -PID1)")
- parser.add_argument("--Ecom", default = "14", help = "Center of mass energy [TeV]")
+ parser.add_argument("--Ecom", default = "13.6", help = "Center of mass energy [TeV]")
parser.add_argument("--outdir", "-o", default = None, help = "Output path")
parser.add_argument("--path", default = ".", help = "Path to foresee installation")
parser.add_argument("--hepmc", action = "store_true", help = "Write HepMC events")
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index d4ff9c8b2cb9274a7aade2e8b0a0f6ceb6bf21e9..b3cfc2b55dd58ddd921c4833f7191fef16e618c3 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -65,6 +65,7 @@ def plotn(f, args, configs, x, y, outname):
c._objs.append(plot(f, *cfg))
c.Print(f"{args.output}/{outname}.eps")
+ c.Print(f"{args.output}/{outname}.png")
return
diff --git a/Generators/ForeseeGenerator/share/validate_grid.py b/Generators/ForeseeGenerator/share/validate_grid.py
index e3d5302dff75a44f6cc2a77497f625e6659706db..88edb1194ee6ecebe37464906b26a03afde7b1da 100644
--- a/Generators/ForeseeGenerator/share/validate_grid.py
+++ b/Generators/ForeseeGenerator/share/validate_grid.py
@@ -1,21 +1,27 @@
import json
import subprocess as proc
-import os
+from glob import glob
+import os, sys
grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
with open(grid) as f:
data = json.load(f)
-name = "DarkPhotonGrid2"
+name = "DarkPhotonGrid13p6_65mm"
path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
-ecom = 14
+ecom = 13.6
pid1 = 11
pid2 = -11
for coup, masses in data["samples"].items():
for m in masses:
- infile = f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):e}to_{pid1}_{pid2}_{name}.hepmc"
+ files = glob(f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}*.hepmc")
+
+ if len(files) != 1:
+ continue
+
+ infile = files[0]
outfile = infile.replace(".hepmc", ".EVNT.pool.root")
valfile = infile.replace(".hepmc", ".VAL.pool.root")
valdir = infile.replace(".hepmc", "")
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
index c64112cc84e23963b63ea381817d1457f1c2a122..9f541da7c972ac01170b8acfafc8001b6eea5a14 100644
--- a/Generators/ParticleGun/python/histsampling.py
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -27,6 +27,9 @@ def load_hist(*args):
h = args[0].Get(args[1]).Clone()
if h is None:
raise Exception("Error in histogram loading from " + args)
+ else:
+ h.SetDirectory(0)
+
return h
@@ -40,7 +43,10 @@ def get_sampling_vars(h):
globalbin_to_axisbin = {} # for reverse axis bin lookup to get edges
globalbins = [] # because they aren't easily predicted, nor contiguous
cheights = [0] # cumulative "histogram" from which to uniformly sample
- if issubclass(type(h), ROOT.TH1):
+
+
+
+ if issubclass(type(h), ROOT.TH1) and not issubclass(type(h), ROOT.TH2):
for ix in range(1, h.GetNbinsX()+1):
iglobal = h.GetBin(ix)
globalbins.append(iglobal)
@@ -109,10 +115,6 @@ class TH1(object):
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th1)
return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin)
- def __getattr__(self, attr):
- "Forward all attributes to the contained TH1"
- return getattr(self.th1, attr)
-
class TH2(object):
"Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
@@ -127,6 +129,4 @@ class TH2(object):
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th2)
return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin)
- def __getattr__(self, attr):
- "Forward other attributes to the contained TH2"
- return getattr(self.th2, attr)
+
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index 90e9676a5ce20dc1b52e48c79f0b36a7da49f639..d6e540ea6cf7a747d9152202cd488c135fa27099 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -1,7 +1,7 @@
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
import ROOT, math, random
-from ParticleGun.histsampling import TH1
+from ParticleGun.histsampling import TH1,TH2
## For convenience
PI = math.pi
@@ -191,13 +191,26 @@ class TH1Sampler(ContinuousSampler):
def __init__(self, *args):
self.hist = TH1(*args)
- if self.hist.GetEntries() < 1:
- raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.GetName())
+ if self.hist.th1.GetEntries() < 1:
+ raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.th1.GetName())
def shoot(self):
return self.hist.GetRandom()
+class TH2Sampler(ContinuousSampler):
+ "Randomly sample from a 2D ROOT histogram."
+
+ def __init__(self, *args):
+ self.hist = TH2(*args)
+ if self.hist.th2.GetEntries() < 1:
+ raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.th2.GetName())
+
+ def shoot(self):
+ return self.hist.GetRandom()
+
+
+
########################################
@@ -572,9 +585,12 @@ class EThetaMPhiSampler(MomSampler):
# TODO: ensure that E >= m!
- def __init__(self, energy, theta, mass=0.0, phi=[0, TWOPI]):
+ def __init__(self, energy, theta=None, mass=0.0, phi=[0, TWOPI]):
self.energy = energy
- self.theta = theta
+ if theta is None:
+ self._theta = None
+ else:
+ self.theta = theta
self.mass = mass
self.phi = phi
@@ -616,10 +632,15 @@ class EThetaMPhiSampler(MomSampler):
pz = p cos(theta)
pt = p sin(theta)
"""
- e = self.energy()
+
+ if self._theta is None:
+ e,theta = self.energy()
+ else:
+ e = self.energy()
+ theta = self.theta()
+
m = self.mass()
p = math.sqrt( e**2 - m**2 )
- theta = self.theta()
pz = p * math.cos(theta)
pt = p * math.sin(theta)
phi = self.phi()
diff --git a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx
index e3d33d683e0b47d064216515cf83fb280cb9edd4..c353833fc5d100423c0a97c778cd5243db289184 100644
--- a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx
+++ b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx
@@ -750,4 +750,4 @@ NeutrinoSearchAlg::clearTree() const
m_syVetoNu = 0;
m_thetaxVetoNu = 0;
m_thetayVetoNu = 0;
-}
\ No newline at end of file
+}
diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
index dc2dae3126aaa00283775626da1b29bbb0dd3ac3..fab78641d7ae2a3842904122968a80170bca9a28 100644
--- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -1,4 +1,4 @@
-atlas_subdir(NtupleDumper)
+atlas_subdir( NtupleDumper )
atlas_add_component(
NtupleDumper
@@ -8,5 +8,5 @@ atlas_add_component(
LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib
)
-atlas_install_python_modules(python/*.py)
-atlas_install_scripts(scripts/*.py)
+atlas_install_python_modules( python/*.py )
+atlas_install_scripts( scripts/*.py scripts/*.sh )
diff --git a/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py
index 5bc52c41df03c934465f97658a39aeb342350631..da44697fc23f3318f599d4451f03fcc196a3f463 100644
--- a/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py
+++ b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py
@@ -2,103 +2,33 @@
Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
"""
-from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg
-def NtupleDumperAlgCfg(flags, **kwargs):
+def NtupleDumperAlgCfg(flags, OutName, **kwargs):
# Initialize GeoModel
from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
acc = FaserGeometryCfg(flags)
acc.merge(MagneticFieldSvcCfg(flags))
- # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
- # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+ #acc.merge(ActsTrackingGeometrySvcCfg(flags))
+ #acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags)
+ acc.merge(result)
actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
- actsExtrapolationTool.MaxSteps = 1000
- actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+ actsExtrapolationTool.MaxSteps = 10000
+ actsExtrapolationTool.TrackingGeometryTool = actsTrackingGeometryTool
NtupleDumperAlg = CompFactory.NtupleDumperAlg("NtupleDumperAlg",**kwargs)
NtupleDumperAlg.ExtrapolationTool = actsExtrapolationTool
acc.addEventAlgo(NtupleDumperAlg)
thistSvc = CompFactory.THistSvc()
- thistSvc.Output += ["HIST2 DATAFILE='Data-tuple.root' OPT='RECREATE'"]
+ thistSvc.Output += [f"HIST2 DATAFILE='{OutName}' OPT='RECREATE'"]
acc.addService(thistSvc)
return acc
-if __name__ == "__main__":
-
- import sys
- from AthenaCommon.Logging import log, logging
- from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
- from AthenaCommon.Configurable import Configurable
- from CalypsoConfiguration.AllConfigFlags import ConfigFlags
- from AthenaConfiguration.TestDefaults import defaultTestFiles
- from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
- from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
- # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
-
- # Set up logging and new style config
- log.setLevel(INFO)
- Configurable.configurableRun3Behavior = True
-
- # Configure
- ConfigFlags.Input.Files = [
- '/eos/experiment/faser/rec/2022/p0008//008119/Faser-Physics-008119-00168-p0008-xAOD.root',
-
-
- ]
- ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
- ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
- ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
- ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
- ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
- ConfigFlags.Common.isOnline = False
- ConfigFlags.GeoModel.Align.Dynamic = False
- ConfigFlags.Beam.NumberOfCollisions = 0.
-
- ConfigFlags.Detector.GeometryFaserSCT = True
-
- ConfigFlags.lock()
-
- # Core components
- acc = MainServicesCfg(ConfigFlags)
- acc.merge(PoolReadCfg(ConfigFlags))
-
- # algorithm
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, UseFlukaWeights=True))
-
- # silencio
- AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
- AthenaEventLoopMgr.EventPrintoutInterval=500
- acc.addService(AthenaEventLoopMgr)
-
- # # Hack to avoid problem with our use of MC databases when isMC = False
- replicaSvc = acc.getService("DBReplicaSvc")
- replicaSvc.COOLSQLiteVetoPattern = ""
- replicaSvc.UseCOOLSQLite = True
- replicaSvc.UseCOOLFrontier = False
- replicaSvc.UseGeomSQLite = True
-
- # Timing
- #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
-
- # Dump config
- # logging.getLogger('forcomps').setLevel(VERBOSE)
- # acc.foreach_component("*").OutputLevel = VERBOSE
- # acc.foreach_component("*ClassID*").OutputLevel = INFO
- # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
- # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
- # acc.getService("StoreGateSvc").Dump = True
- # acc.getService("ConditionStore").Dump = True
- # acc.printConfig(withDetails=True)
- # ConfigFlags.dump()
-
- # Execute and finish
- sc = acc.run(maxEvents=-1)
-
- # Success should be 0
- sys.exit(not sc.isSuccess())
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
index 8a47857614ebb70589ed7cb5ab6ca20aef91c878..8cb79460633e7cf16ba5818a43206c6dbe7c20ab 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
@@ -6,31 +6,9 @@
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
-from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
-
-def NtupleDumperAlgCfg(flags, OutName, **kwargs):
- # Initialize GeoModel
- from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
- acc = FaserGeometryCfg(flags)
-
- acc.merge(MagneticFieldSvcCfg(flags))
- # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
- # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
-
- actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
- actsExtrapolationTool.MaxSteps = 10000
- actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
-
- NtupleDumperAlg = CompFactory.NtupleDumperAlg("NtupleDumperAlg",**kwargs)
- NtupleDumperAlg.ExtrapolationTool = actsExtrapolationTool
- acc.addEventAlgo(NtupleDumperAlg)
-
- thistSvc = CompFactory.THistSvc()
- thistSvc.Output += [f"HIST2 DATAFILE='{OutName}' OPT='RECREATE'"]
- acc.addService(thistSvc)
-
- return acc
+# Move NtupleDumperAlgCfg to python/NtupleDumperConfig
+from NtupleDumper.NtupleDumperConfig import NtupleDumperAlgCfg
if __name__ == "__main__":
@@ -56,30 +34,41 @@ if __name__ == "__main__":
Configurable.configurableRun3Behavior = True
dataDir=f"/eos/experiment/faser/rec/2022/{ptag}/{runno:06d}"
- if inputIsMC: dataDir=f"/eos/experiment/faser/sim/mdc/foresee/{ui.runnum}/rec/r0009/"
+
files=sorted(glob.glob(f"{dataDir}/Faser*"))
- if num==-1: num=len(files)
- fileListInitial=files[0:num]#[num*filesPerJob:(num+1)*filesPerJob]
+ if filesPerJob <= 0: filesPerJob = len(files)
+
+ start = num*filesPerJob
+ end = (num+1)*filesPerJob
+ if end > len(files):
+ end = len(files)
+ fileListInitial=files[start:end]
fileList=[]
for fName in fileListInitial:
+ #if fName[:4] == '/eos':
+ # fName = f'root:/{fName}'
+
try:
fh=ROOT.TFile(fName)
fileList.append(fName)
except OSError:
print("Warning bad file: ",fName)
- log.info(f"Analyzing Run {runno} files {0} to {num} (num={num})")
+ log.info(f"Analyzing Run {runno} files {start} to {end} (num={num})")
log.info(f"Got {len(fileList)} files out of {len(fileListInitial)}")
- outName=f"Data-tuple-{runno:06d}-{num:05d}-{filesPerJob}.root"
+ if start == 0 and end == len(files):
+ outName=f"Faser-Physics-{runno:06d}-{ptag}-PHYS.root"
+ else:
+ outName=f"Faser-Physics-{runno:06d}-{start:05d}-{(end-1):05d}-{ptag}-PHYS.root"
# Configure
ConfigFlags.Input.Files = fileList
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
- ConfigFlags.Input.isMC = inputIsMC # Needed to bypass autoconfig
- ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
ConfigFlags.Common.isOnline = False
ConfigFlags.GeoModel.Align.Dynamic = False
ConfigFlags.Beam.NumberOfCollisions = 0.
@@ -93,36 +82,36 @@ if __name__ == "__main__":
acc.merge(PoolReadCfg(ConfigFlags))
# algorithm
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outName, UseFlukaWeights=True))
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outName))
AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
AthenaEventLoopMgr.EventPrintoutInterval=1000
acc.addService(AthenaEventLoopMgr)
# # Hack to avoid problem with our use of MC databases when isMC = False
- if not inputIsMC:
- replicaSvc = acc.getService("DBReplicaSvc")
- replicaSvc.COOLSQLiteVetoPattern = ""
- replicaSvc.UseCOOLSQLite = True
- replicaSvc.UseCOOLFrontier = False
- replicaSvc.UseGeomSQLite = True
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
# Timing
#acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
# Dump config
- # logging.getLogger('forcomps').setLevel(VERBOSE)
- # acc.foreach_component("*").OutputLevel = VERBOSE
- # acc.foreach_component("*ClassID*").OutputLevel = INFO
- # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
- # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
- # acc.getService("StoreGateSvc").Dump = True
- # acc.getService("ConditionStore").Dump = True
- # acc.printConfig(withDetails=True)
- # ConfigFlags.dump()
+ if False:
+ logging.getLogger('forcomps').setLevel(VERBOSE)
+ acc.foreach_component("*").OutputLevel = VERBOSE
+ acc.foreach_component("*ClassID*").OutputLevel = INFO
+ acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ acc.getService("StoreGateSvc").Dump = True
+ acc.getService("ConditionStore").Dump = True
+ acc.printConfig(withDetails=True)
+ ConfigFlags.dump()
# Execute and finish
- sc = acc.run(maxEvents=ui.events)
+ sc = acc.run(maxEvents=-1)
# Success should be 0
sys.exit(not sc.isSuccess())
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
new file mode 100755
index 0000000000000000000000000000000000000000..27b30e25df16dd8baf1c4ce4b82b6a9e6a0f7b6f
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -0,0 +1,284 @@
+#!/usr/bin/env python
+#
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+#
+# Run with
+# faser_ntuple_maker.py [options]
+#
+# Options:
+# --isMC - needed to reconstruct MC data
+#
+import sys
+import time
+import argparse
+
+a = time.time()
+
+parser = argparse.ArgumentParser(description="Run PHYS ntuple production")
+
+parser.add_argument("path",
+ help="Fully qualified input path (directory or file)")
+
+parser.add_argument("--slice", type=int, default=0,
+ help="Specify ordinal output file to produce")
+parser.add_argument("--files", type=int, default=1,
+ help="Specify reco files per slice")
+parser.add_argument("--merge", type=int, default=1,
+ help="Specify merged files per reco file (MC only)")
+parser.add_argument("--last", type=int, default=0,
+ help="Specify last file in slice (normally --files)")
+
+parser.add_argument("--outfile", default="",
+ help="Override output file name")
+
+parser.add_argument("-v", "--verbose", action='store_true',
+ help="Turn on DEBUG output")
+parser.add_argument("-n", "--nevents", type=int, default=-1,
+ help="Specify number of events to process (default: all)")
+
+parser.add_argument("-t", "--tag", default="",
+ help="Specify tag (to append to output filename)")
+
+parser.add_argument("--isMC", action='store_true',
+ help="Running on digitised MC rather than data")
+parser.add_argument("--partial", action='store_true',
+ help="Allow partial input files")
+
+parser.add_argument("--unblind", action='store_true',
+ help="Don't apply signal blinding (default: do)")
+
+parser.add_argument("--fluka", action='store_true',
+ help="Add FLUKA weights to ntuple")
+parser.add_argument("--genie", action='store_true',
+ help="Add Genie weights to ntuple")
+
+args = parser.parse_args()
+
+from pathlib import Path
+
+# Must figure out what we are doing here
+filepath = Path(args.path)
+
+filelist = []
+# If this is a directory, need to create file list
+if filepath.is_dir():
+
+ # Use expected data pattern to find files (data only)
+ runstr = filepath.stem
+
+ # Make list of segments to search for
+ seglist = []
+ if args.merge > 1:
+ start = args.slice * args.files * args.merge
+
+ # Number of files to combine
+ if args.last > 0:
+ num = args.last
+ else:
+ num = args.files
+
+ # Make segment list
+ for i in range(start, start+num*args.merge, args.merge):
+ seg = f"{i:05d}-{(i+args.merge-1):05d}"
+ seglist.append(seg)
+
+ else:
+ start = args.slice * args.files
+ if args.last > 0:
+ end = start + args.last
+ else:
+ end = start + args.files
+
+ seglist = [f'{seg:05d}' for seg in range(start, end)]
+
+
+ for seg in seglist:
+
+ if args.isMC:
+ searchstr = f"FaserMC-*-{seg}-*xAOD.root"
+ else:
+ searchstr = f"Faser-Physics-{runstr}-{seg}-*xAOD.root"
+
+ flist = list(filepath.glob(searchstr))
+ if len(flist) == 0:
+ print(f"Didn't find file {searchstr}!")
+ if args.partial: continue
+ sys.exit(1)
+
+ elif len(flist) > 1:
+ print(f"Found multiple matches for {searchstr}!")
+ print(flist)
+ sys.exit(1)
+
+ filestr = str(flist[0].resolve())
+ # Use proper EOS file path here?
+ if filestr[:4] == '/eos':
+ filestr = f"root://eospublic.cern.ch/{filestr}"
+ filelist.append(filestr)
+ # End of loop over segments
+
+ # Parse name to create outfile
+ firstfile = Path(filelist[0])
+ firststem = str(firstfile.stem)
+ firstfaser = firststem.split('-')[0]
+ firstshort = firststem.split('-')[1]
+ runstr = firststem.split('-')[2]
+ firstseg = firststem.split('-')[3]
+ if args.merge > 1:
+ firstseg2 = firststem.split('-')[4]
+
+ lastfile = Path(filelist[-1])
+ laststem = str(lastfile.stem)
+ lastseg = laststem.split('-')[3]
+ if args.merge > 1:
+ lastseg = laststem.split('-')[4]
+
+ print(f"Faser = {firstfaser}")
+ print(f"Short = {firstshort}")
+ print(f"Run = {runstr}")
+ print(f"First = {firstseg}")
+ print(f"Last = {lastseg}")
+ print(f"Args = {args.tag}")
+ print(f"Blind = {not args.unblind}")
+
+ # Find any tags
+ tagstr = firststem.replace(f"{firstfaser}-{firstshort}-{runstr}-{firstseg}", "")
+ if args.merge > 1:
+ tagstr = tagstr.replace(f"-{firstseg2}", "")
+
+ tagstr = tagstr.replace("-xAOD", "")
+
+ print(f"Tag = {tagstr}")
+
+ # Build output name
+ outfile = f"{firstfaser}-{firstshort}-{runstr}-{firstseg}-{lastseg}"
+
+ # This will include the leading -
+ if len(tagstr) > 0:
+ outfile += f"{tagstr}"
+
+ if len(args.tag) > 0 and args.tag not in tagstr:
+ outfile += f"-{args.tag}"
+
+ outfile += "-PHYS.root"
+
+
+# If this is a single file, just process that
+# Could be a url, so don't check if this is a file
+else:
+ filelist.append(args.path)
+
+ # Build output name
+ filestem = str(filepath.stem)
+
+ # Remove -xAOD from end
+ if filestem[-5:] == "-xAOD":
+ filestem = filestem[:-5]
+
+ # Do we want to add a tag?
+ if len(args.tag) > 0 and args.tag not in filestem:
+ filestem += f"-{args.tag}"
+
+ outfile = f"{filestem}-PHYS.root"
+
+# Print out what we found
+if len(filelist) == 0:
+ printf("Found no files for {args.path}!")
+ sys.exit(1)
+elif len(filelist) == 1:
+ print("Processing file:")
+else:
+ print("Processing files:")
+for f in filelist:
+ print(f)
+
+# Override output file?
+if len(args.outfile) > 0:
+ outfile = args.outfile
+
+print("Output file:")
+print(outfile)
+
+# OK, lets run the job here
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+
+Configurable.configurableRun3Behavior = True
+
+# Configure
+ConfigFlags.Input.Files = filelist
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ConfigFlags.Input.isMC = args.isMC
+if args.isMC:
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions
+else:
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions
+
+ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+
+ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ConfigFlags.Common.isOnline = False
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.Detector.GeometryFaserSCT = True
+ConfigFlags.lock()
+
+# Core components
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+
+# algorithm
+from NtupleDumper.NtupleDumperConfig import NtupleDumperAlgCfg
+if args.isMC:
+ if args.genie:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, UseGenieWeights=True))
+ elif args.fluka:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, UseFlukaWeights=True))
+ else:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile))
+
+else:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind)))
+
+if not args.verbose:
+ from AthenaConfiguration.ComponentFactory import CompFactory
+ AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+ AthenaEventLoopMgr.EventPrintoutInterval=1000
+ acc.addService(AthenaEventLoopMgr)
+
+# Hack to avoid problem with our use of MC databases when isMC = False
+if not args.isMC:
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+if args.verbose:
+ log.setLevel(VERBOSE)
+ acc.printConfig(withDetails=True)
+ ConfigFlags.dump()
+else:
+ log.setLevel(INFO)
+
+acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M"
+
+# Execute and finish
+sc = acc.run(maxEvents=args.nevents)
+
+b = time.time()
+log.info(f"Finish execution in {b-a} seconds")
+
+# Success should be 0
+if sc.isSuccess():
+ log.info("Execution succeeded")
+ sys.exit(0)
+else:
+ log.info("Execution failed, return 1")
+ sys.exit(1)
+
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh b/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9e3273dbb970f5b5d80958f67656490c045a7612
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh
@@ -0,0 +1,282 @@
+#!/bin/bash
+#
+# Used with a condor file to submit a physics ntuple job
+#
+# Usage:
+# submit_ntuple_maker.sh dirpath slice nfiles [release_directory] [working_directory]
+#
+# Options:
+# --out - specify output location (in EOS) to copy output HITS file
+# --log - specify output location (in EOS) for log file
+# --isMC - needed for MC reco
+# --partial - allow missing files
+# --merge - merge factor of reco files (for MC)
+# --fluka - create fluka weights
+# --genie - create genie weights
+#
+# dirpath - full directory path to HITS files
+# slice - ordinal output file number
+# nfiles - number of HITS files to process per slice
+# release_directory - optional path to release install directory (default pwd)
+# working_directory - optional path to output directory location (default pwd)
+# lastfile - number of files in this slice (for last file, default 0)
+#
+# The release directory must already be set up
+# (so an unqualified asetup can set up the release properly)
+#
+#----------------------------------------
+# Keep track of time
+SECONDS=0
+#
+# Defaults
+ismc=""
+partialstr=""
+mergestr=""
+flukastr=""
+geniestr=""
+unblindstr=""
+#
+# Parse command-line options
+while [ -n "$1" ]
+do
+ case "$1" in
+ -l | --log)
+ logdest="$2";
+ shift;
+ shift;; # Must eat 2 options here
+
+ -o | --out)
+ outdest="$2";
+ shift;
+ shift;;
+
+ --isMC)
+ ismc="--isMC"
+ shift;;
+
+ --partial)
+ echo "Allowing partial merge"
+ partialstr="--partial"
+ shift;;
+
+ --merge)
+ mergestr="--merge $2";
+ shift;
+ shift;;
+
+ --fluka)
+ flukastr="--fluka";
+ shift;;
+
+ --genie)
+ geniestr="--genie";
+ shift;;
+
+ --unblind)
+ unblindstr="--unblind";
+ shift;;
+
+ --) # End of options
+ shift; # Eat this
+ break;; # And stop parsing
+
+ -*)
+ echo "Unknown option $1"
+ shift;;
+
+ *) break;; # Not an option, don't shift
+ esac
+done
+#
+# Parse command-line options
+dir_path=${1}
+slice=${2}
+nfiles=${3}
+release_directory=${4}
+working_directory=${5}
+last_file=${6}
+#
+# Set defaults if arguments aren't provided
+if [ -z "$dir_path" ]
+then
+ echo "No directory specified!"
+ echo "Usage: submit_ntuple_maker.sh directory slice nfiles [release dir] [output dir]"
+ exit 1
+fi
+#
+if [ -z "$slice" ]
+then
+ echo "Slice number not specified!"
+ echo "Usage: submit_ntuple_maker.sh directory slice nfiles [release dir] [output dir]"
+ exit 1
+fi
+#
+if [ -z "$nfiles" ]
+then
+ echo "Files per slice not specified!"
+ echo "Usage: submit_ntuple_maker.sh directory slice nfiles [release dir] [output dir]"
+ exit 1
+fi
+#
+if [ -z "$release_directory" ]
+then
+ release_directory=`pwd`
+fi
+#
+if [ -z "$working_directory" ]
+then
+ working_directory=`pwd`
+fi
+#
+starting_directory=`pwd`
+#
+# Now extract the run number and file stem
+#
+# First, get an example filename
+file_name=`ls -1 $dir_path | head -1`
+#
+# Now split based on '.' to get stem
+defaultIFS=$IFS
+IFS='.'
+read file_stem ext <<< "$file_name"
+#
+# Finally extract the run number
+IFS='-'
+# Read the split words into an array based on delimiter
+read faser short run_number segment <<< "$file_stem"
+#
+# Set the IFS delimeter back or else echo doesn't work...
+IFS=$defaultIFS
+#
+# Make output directory if needed
+output_directory="$working_directory/$run_number"
+mkdir -p "$output_directory"
+#
+# Need to make up an output name
+printf -v slicestr "%04d" $slice
+#
+file_stem="$faser-$short-$run_number-$slicestr-PHYS"
+#
+# This magic redirects everything in this script to our log file
+logfile="${file_stem}.ntp.log"
+exec >& "$output_directory/$logfile"
+echo `date` - $HOSTNAME
+echo "Directory: $dir_path"
+echo "Slice: $slice"
+echo "NFiles: $nfiles"
+echo "Merge: $mergestr"
+echo "Release: $release_directory"
+echo "Output: $output_directory"
+echo "Starting: $starting_directory"
+echo "job: $file_stem"
+echo "Last: $last_file"
+#
+# Set up the release (do this automatically)?
+export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase
+source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh
+#
+# Always go back to the starting directory in case paths are relative
+cd "$starting_directory"
+cd "$release_directory"
+#
+# Do this by hand
+asetup --input=calypso/asetup.faser Athena,22.0.49
+source run/setup.sh
+#
+echo "ATLAS_POOLCOND_PATH = $ATLAS_POOLCOND_PATH"
+#
+# Check if there are data overrides
+if [ -d "run/data/sqlite200" ]; then
+ cond_directory=`pwd -P`/run/data # Get absolute path
+fi
+#
+# Try to find a release tag
+cd calypso
+recotag=`git describe --tags`
+if [[ "$recotag" == "reco/r"???? ]]; then
+ tag=`echo "$recotag" | cut -c 6-11`
+ echo "Found reco tag: $tag"
+fi
+if [[ "$recotag" == "digi/d"???? ]]; then
+ tag=`echo "$recotag" | cut -c 6-11`
+ echo "Found digi tag: $tag"
+fi
+if [[ "$recotag" == "sim/s"???? ]]; then
+ tag=`echo "$recotag" | cut -c 5-10`
+ echo "Found sim tag: $tag"
+fi
+#
+if [[ -z "$tag" ]]; then
+ tagstr=""
+else
+ tagstr="--tag $tag"
+fi
+#
+if [[ -z "$last_file" ]]; then
+ last_file_str=""
+else
+ last_file_str="--last $last_file"
+fi
+#
+# Move to the run directory
+cd "$starting_directory"
+cd "$output_directory"
+#
+# Make run directory
+if [[ -e "$file_stem" ]]; then
+ echo "Directory $file_stem already exists"
+else
+ mkdir "$file_stem"
+fi
+cd "$file_stem"
+#
+# Check if there are data overrides in the relese directory
+if [[ -z "$cond_directory" ]]; then
+ echo "No local conditions directory found!"
+else
+ echo "Local conditions directory found! Copying to run directory..."
+ echo Copying $cond_directory
+ cp -r $cond_directory .
+ ls -R data
+fi
+#
+export EOS_MGM_URL=root://eospublic.cern.ch
+#
+# Run job
+#
+faser_ntuple_maker.py $last_file_str $partialstr $tagstr $ismc --slice $slice --files $nfiles $mergestr $flukastr $geniestr $unblindstr $dir_path
+ntup_code=$?
+echo "Return code: $ntup_code"
+#
+# Print out ending time
+date
+echo "Job finished after $SECONDS seconds"
+#
+# Copy output to EOS if desired
+if ! [ -z "$outdest" ]
+then
+ ls -l
+ echo "copy *-PHYS.root to $outdest"
+ mkdir -p $outdest
+ eos cp *-PHYS.root ${outdest}/ || true
+fi
+#
+# Also copy log file
+if ! [ -z "$logdest" ]
+then
+ cd ..
+ ls -l
+ echo "copy $logfile to $logdest"
+ mkdir -p $logdest
+ eos cp $logfile $logdest/$logfile
+elif ! [ -z "$outdest" ]
+then
+ cd ..
+ ls -l
+ echo "copy $logfile to $outdest"
+ mkdir -p $outdest
+ eos cp $logfile $outdest/$logfile
+fi
+
+# Make sure to return an error if ntuple failed
+exit $ntup_code
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 1369684d16f1071c2b5eab6e5853b87459d0e06a..3f476edb40480e5871a2df39a3cd2695ea894e4c 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -106,6 +106,7 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_preshowerCalibratedContainer.initialize());
ATH_CHECK(m_ecalCalibratedContainer.initialize());
+ ATH_CHECK(m_eventInfoKey.initialize());
ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
@@ -160,8 +161,8 @@ StatusCode NtupleDumperAlg::initialize()
//WAVEFORMS
addWaveBranches("VetoNu",2,4);
- addWaveBranches("VetoSt1",2,6);
- addWaveBranches("VetoSt2",1,14);
+ addWaveBranches("VetoSt1",1,14);
+ addWaveBranches("VetoSt2",2,6);
addWaveBranches("Timing",4,8);
addWaveBranches("Preshower",2,12);
addWaveBranches("Calo",4,0);
@@ -197,6 +198,7 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("TrackSegment_pz", &m_trackseg_pz);
//TrackCollection
+ m_tree->Branch("Track_PropagationError", &m_propagationError, "Track_PropagationError/I");
m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
m_tree->Branch("Track_Chi2", &m_Chi2);
m_tree->Branch("Track_nDoF", &m_DoF);
@@ -304,6 +306,24 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I");
m_tree->Branch("CrossSection", &m_crossSection, "crossSection/D");
+
+ // first 10 truth particles
+
+ m_tree->Branch("truth_P", &m_truth_P);
+ m_tree->Branch("truth_px", &m_truth_px);
+ m_tree->Branch("truth_py", &m_truth_py);
+ m_tree->Branch("truth_pz", &m_truth_pz);
+ m_tree->Branch("truth_m", &m_truth_m);
+ m_tree->Branch("truth_pdg", &m_truth_pdg);
+
+ m_tree->Branch("truth_prod_x", &m_truth_prod_x);
+ m_tree->Branch("truth_prod_y", &m_truth_prod_y);
+ m_tree->Branch("truth_prod_z", &m_truth_prod_z);
+
+ m_tree->Branch("truth_dec_x", &m_truth_dec_x);
+ m_tree->Branch("truth_dec_y", &m_truth_dec_y);
+ m_tree->Branch("truth_dec_z", &m_truth_dec_z);
+
// for mother + daughter particle truth infomation
m_tree->Branch("truthM_P", &m_truthM_P);
@@ -410,6 +430,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
bool trig_coincidence_vetoesORpreshower_and_calo = false;
if ( ((m_tap&2)!=0) && ((m_tap&1)!=0) ) trig_coincidence_vetoesORpreshower_and_calo = true;
+ bool trig_calo = (m_tap & 1);
+
// for random trigger, store charge of scintillators in histograms
if (trig_random) {
// Read in Waveform containers
@@ -462,7 +484,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
return StatusCode::SUCCESS; // finished with this randomly triiggered event
- } else if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_coincidence_timing_and_calo || trig_coincidence_vetoesORpreshower_and_calo) ) {
+ // Try adding calo-only trigger
+ } else if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_coincidence_timing_and_calo || trig_coincidence_vetoesORpreshower_and_calo || trig_calo) ) {
// don't process events that fail to activate coincidence triggers
return StatusCode::SUCCESS;
}
@@ -546,30 +569,64 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
{
- for (auto particle : *truthParticleContainer)
- {
- if ( particle->barcode() <= 3)
- {
- if ( particle->barcode() == 1) // mother particle (A')
- {
-
- m_truthM_P.push_back(particle->p4().P());
- m_truthM_px.push_back(particle->p4().X());
- m_truthM_py.push_back(particle->p4().Y());
- m_truthM_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthM_x.push_back(particle->prodVtx()->x());
- m_truthM_y.push_back(particle->prodVtx()->y());
- m_truthM_z.push_back(particle->prodVtx()->z());
- } else {
+ int ipart(0);
+ for (auto particle : *truthParticleContainer)
+ {
+
+ // loop over first 10 truth particles (for non A' samples)
+
+ if (ipart++ < 10) {
+
+ m_truth_P.push_back(particle->p4().P());
+ m_truth_px.push_back(particle->p4().X());
+ m_truth_py.push_back(particle->p4().Y());
+ m_truth_pz.push_back(particle->p4().Z());
+ m_truth_m.push_back(particle->m());
+ m_truth_pdg.push_back(particle->pdgId());
+
+ if ( particle->hasProdVtx()) {
+ m_truth_prod_x.push_back(particle->prodVtx()->x());
+ m_truth_prod_y.push_back(particle->prodVtx()->y());
+ m_truth_prod_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truth_prod_x.push_back(NaN);
+ m_truth_prod_y.push_back(NaN);
+ m_truth_prod_z.push_back(NaN);
+ }
+
+ if ( particle->hasDecayVtx()) {
+ m_truth_dec_x.push_back(particle->decayVtx()->x());
+ m_truth_dec_y.push_back(particle->decayVtx()->y());
+ m_truth_dec_z.push_back(particle->decayVtx()->z());
+ } else {
+ m_truth_dec_x.push_back(NaN);
+ m_truth_dec_y.push_back(NaN);
+ m_truth_dec_z.push_back(NaN);
+ }
+ }
+
+
+ if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
+
+ if ( particle->pdgId() == 32) { // mother particle (A')
+
+ m_truthM_P.push_back(particle->p4().P());
+ m_truthM_px.push_back(particle->p4().X());
+ m_truthM_py.push_back(particle->p4().Y());
+ m_truthM_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasDecayVtx()) { // decay vertex for A' particle
+ m_truthM_x.push_back(particle->decayVtx()->x());
+ m_truthM_y.push_back(particle->decayVtx()->y());
+ m_truthM_z.push_back(particle->decayVtx()->z());
+ } else {
m_truthM_x.push_back(NaN);
m_truthM_y.push_back(NaN);
m_truthM_z.push_back(NaN);
}
+ }
- }
if ( particle->pdgId() == 11) // daughter particle (positron)
{
m_truthd0_P.push_back(particle->p4().P());
@@ -604,12 +661,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_truthd1_z.push_back(NaN);
}
}
- }
+ }
}
- }
+ }
}
// load in calibrated calo container
+ if (m_doCalib) {
SG::ReadHandle<xAOD::CalorimeterHitContainer> ecalCalibratedContainer { m_ecalCalibratedContainer, ctx };
ATH_CHECK(ecalCalibratedContainer.isValid());
for (auto hit : *ecalCalibratedContainer) {
@@ -649,6 +707,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("Calibrated preshower: ch is " << ch << ", edep is " << hit->E_dep() << ", E_EM is " << hit->E_EM() << ", Nmip is " << hit->Nmip() << ", fit_to_raw_ratio is " << hit->fit_to_raw_ratio());
}
+ }
// process all waveeform data fro all scintillator and calorimeter channels
SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
@@ -673,14 +732,23 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
FillWaveBranches(*ecalContainer);
// enforce blinding such that events that miss the veto layers and have a large calo signal are skipped and not in the output root file
- if ( (!isMC) && m_doBlinding) {
- if ( m_calo_total_E_EM/1000.0 > 10.0 ) { // energy is in MeV so divide by 1000 to compare to 10 GeV
+ // Only blind colliding BCIDs (+/- 1)
+ if ( (!isMC) && m_doBlinding && abs(m_distanceToCollidingBCID) <= 1) {
+ if (m_calo_total_E_EM > m_blindingCaloThreshold ) { // energy is in MeV
if (m_wave_status[4] == 1 and m_wave_status[5] == 1 and m_wave_status[6] == 1 and m_wave_status[7] == 1 and m_wave_status[14] == 1) { // hit status == 1 means it is below threshold. channles 4 and 5 are vetoNu, channels 6,7, and 14 are veto
return StatusCode::SUCCESS;
}
}
}
+ SG::ReadDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx);
+ if (eventInfo->errorState(xAOD::EventInfo_v1::SCT) == xAOD::EventInfo::Error) {
+ m_propagationError = 1;
+ ATH_MSG_DEBUG("NtupleDumper: xAOD::EventInfo::SCT::Error");
+ } else {
+ m_propagationError = 0;
+ }
+
// get geometry context
FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
auto gctx = faserGeometryContext.context();
@@ -766,8 +834,16 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// loop over all reconstructed tracks and use only the tracks that have hits in all three tracking stations (excludes IFT)
// store track parameters at most upstream measurement and at most downstream measurement
// extrapolate track to all scintillator positions and store extrapolated position and angle
- SG::ReadHandle<TrackCollection> trackCollection {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
+ SG::ReadHandle<TrackCollection> trackCollection;
+ if (m_useIFT) {
+ SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that excludes IFT
+ trackCollection = tc;
+ } else {
+ SG::ReadHandle<TrackCollection> tc {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
+ trackCollection = tc;
+ }
ATH_CHECK(trackCollection.isValid());
+
for (const Trk::Track* track : *trackCollection)
{
if (track == nullptr) continue;
@@ -1218,6 +1294,44 @@ NtupleDumperAlg::clearTree() const
m_truthLeptonMomentum = 0;
m_truthBarcode = 0;
m_truthPdg = 0;
+
+ m_truth_P.clear();
+ m_truth_px.clear();
+ m_truth_py.clear();
+ m_truth_pz.clear();
+ m_truth_m.clear();
+ m_truth_pdg.clear();
+ m_truth_prod_x.clear();
+ m_truth_prod_y.clear();
+ m_truth_prod_z.clear();
+ m_truth_dec_x.clear();
+ m_truth_dec_y.clear();
+ m_truth_dec_z.clear();
+
+ m_truthM_P.clear();
+ m_truthM_px.clear();
+ m_truthM_py.clear();
+ m_truthM_pz.clear();
+ m_truthM_x.clear();
+ m_truthM_y.clear();
+ m_truthM_z.clear();
+
+ m_truthd0_P.clear();
+ m_truthd0_px.clear();
+ m_truthd0_py.clear();
+ m_truthd0_pz.clear();
+ m_truthd0_x.clear();
+ m_truthd0_y.clear();
+ m_truthd0_z.clear();
+
+ m_truthd1_P.clear();
+ m_truthd1_px.clear();
+ m_truthd1_py.clear();
+ m_truthd1_pz.clear();
+ m_truthd1_x.clear();
+ m_truthd1_y.clear();
+ m_truthd1_z.clear();
+
}
void NtupleDumperAlg::setNaN() const {
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 5cc342021993ca07354e6a4924f7e5cfb07ae552..6d4ccc381a4b65831678a657c035ecc4e79c2b65 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -21,6 +21,8 @@
#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "StoreGate/ReadDecorHandle.h"
#include <vector>
@@ -84,7 +86,8 @@ private:
SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
SG::ReadHandleKey<xAOD::WaveformClock> m_ClockWaveformContainer { this, "WaveformClockKey", "WaveformClock", "ReadHandleKey for ClockWaveforms Container"};
- ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+ SG::ReadDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
@@ -98,6 +101,8 @@ private:
const PreshowerID* m_preshowerHelper;
const EcalID* m_ecalHelper;
+ BooleanProperty m_useIFT { this, "UseIFT", false, "Use IFT tracks" };
+ BooleanProperty m_doCalib { this, "DoCalib", true, "Fill calibrated calorimeter quantities" };
BooleanProperty m_doBlinding { this, "DoBlinding", true, "Blinding will not output events with Calo signal > 10 GeV e-" };
BooleanProperty m_useFlukaWeights { this, "UseFlukaWeights", false, "Flag to weight events according to value stored in HepMC::GenEvent" };
BooleanProperty m_useGenieWeights { this, "UseGenieWeights", false, "Flag to weight events according to Genie luminosity" };
@@ -106,6 +111,8 @@ private:
DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
+ DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"};
+
double m_baseEventCrossSection {1.0};
const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -189,6 +196,7 @@ private:
mutable std::vector<double> m_trackseg_pz;
mutable int m_longTracks;
+ mutable int m_propagationError;
mutable std::vector<double> m_Chi2;
mutable std::vector<double> m_DoF;
mutable std::vector<double> m_xup;
@@ -274,7 +282,7 @@ private:
mutable std::vector<double> m_truthM_py;
mutable std::vector<double> m_truthM_pz;
- mutable std::vector<double> m_truthM_x;
+ mutable std::vector<double> m_truthM_x; // decay vertex of A'
mutable std::vector<double> m_truthM_y;
mutable std::vector<double> m_truthM_z;
@@ -284,7 +292,7 @@ private:
mutable std::vector<double> m_truthd0_py;
mutable std::vector<double> m_truthd0_pz;
- mutable std::vector<double> m_truthd0_x;
+ mutable std::vector<double> m_truthd0_x; // production vertex for daughter particles
mutable std::vector<double> m_truthd0_y;
mutable std::vector<double> m_truthd0_z;
@@ -298,6 +306,27 @@ private:
mutable std::vector<double> m_truthd1_y;
mutable std::vector<double> m_truthd1_z;
+ // first 10 truth particles
+
+ mutable std::vector<double> m_truth_P;
+ mutable std::vector<double> m_truth_px;
+ mutable std::vector<double> m_truth_py;
+ mutable std::vector<double> m_truth_pz;
+ mutable std::vector<double> m_truth_m;
+
+ mutable std::vector<double> m_truth_dec_x; // components of decay vertex (mm)
+ mutable std::vector<double> m_truth_dec_y;
+ mutable std::vector<double> m_truth_dec_z;
+
+ mutable std::vector<double> m_truth_prod_x; // components of production vertex (mm)
+ mutable std::vector<double> m_truth_prod_y;
+ mutable std::vector<double> m_truth_prod_z;
+
+ mutable std::vector<int> m_truth_pdg; // pdg of first 10 truth particles
+
+
+
+
mutable double m_truthLeptonMomentum;
mutable int m_truthBarcode;
mutable int m_truthPdg;
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index ea03a77f16d12cd4cb45c99c96a20c1fbc8323e7..1869860e84a401fffbfe6a588d6b8e78ba702e9a 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -3,6 +3,7 @@
#include "ScintSimEvent/ScintHitIdHelper.h"
#include <map>
+#include <cmath>
ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name,
@@ -110,42 +111,68 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
waveforms = m_digiTool->create_waveform_map(m_preshowerID);
}
+ // Sum energy for each channel (i.e. identifier)
+ std::map<Identifier, float> esum;
+
+ Identifier id;
+ for (const auto& hit : *scintHitHandle) {
+ if (first.isTrigger()) {
+ Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+ // Need to do something better here, but for now just fill both
+ id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+ esum[id] += hit.energyLoss();
+ id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+ esum[id] += hit.energyLoss();
+ } else {
+ // All others have only 1 PMT
+ // Use detector id (not hit id) throughout
+ id = hit.getIdentifier();
+ esum[id] += hit.energyLoss();
+ }
+ }
+
// Loop over time samples
- // Should really do sums first, as repeating these in the loop is a waste
for (const auto& tk : m_timekernel) {
std::map<Identifier, float> counts;
- // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
- Identifier id;
- for (const auto& hit : *scintHitHandle) {
-
- // Special handling for trigger scintillator
- if (first.isTrigger()) {
- // Hits are created per plate
- // Need to split between 2 PMTs
-
- // Identifier from hit is plate ID
- Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
-
- // Need to do something better here, but for now just fill both
- id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
- counts[id] += tk * hit.energyLoss();
- id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
- counts[id] += tk * hit.energyLoss();
-
- } else {
- // All others have only 1 PMT
- // Use detector id (not hit id) throughout
- id = hit.getIdentifier();
- counts[id] += tk * hit.energyLoss();
- }
+// // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
+// Identifier id;
+// for (const auto& hit : *scintHitHandle) {
+//
+// // Special handling for trigger scintillator
+// if (first.isTrigger()) {
+// // Hits are created per plate
+// // Need to split between 2 PMTs
+//
+// // Identifier from hit is plate ID
+// Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+//
+// // Need to do something better here, but for now just fill both
+// id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+// counts[id] += tk * hit.energyLoss();
+// id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+// counts[id] += tk * hit.energyLoss();
+//
+// } else {
+// // All others have only 1 PMT
+// // Use detector id (not hit id) throughout
+// id = hit.getIdentifier();
+// counts[id] += tk * hit.energyLoss();
+// }
+// }
+
+ // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+ for (const auto& e : esum) {
+ // Convert hit id to Identifier and store
+ id = e.first;
+ counts[id] = tk * e.second;
}
// Subtract count from basleine and add result to correct waveform vector
for (const auto& c : counts) {
- unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = baseline - c.second;
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = std::round(baseline - c.second);
if (value < 0) {
ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
diff --git a/Tracker/TrackerAlignTools/TrackerAlignGenTools/src/TrackerAlignDBTool.cxx b/Tracker/TrackerAlignTools/TrackerAlignGenTools/src/TrackerAlignDBTool.cxx
index 4483eb5b4cd52b31aa70cb915ba6500641297685..227a1d7c2157aa09c54f2803b2654bec924610e0 100644
--- a/Tracker/TrackerAlignTools/TrackerAlignGenTools/src/TrackerAlignDBTool.cxx
+++ b/Tracker/TrackerAlignTools/TrackerAlignGenTools/src/TrackerAlignDBTool.cxx
@@ -209,8 +209,8 @@ StatusCode TrackerAlignDBTool::createDB() const
{
const auto iStation = m_sctid->station(ident);
const auto iLayer = m_sctid->layer(ident);
- const auto iModuleEta = m_sctid->eta_module(ident);
- const auto iModulePhi = m_sctid->phi_module(ident);
+ const auto iModuleEta = m_sctid->eta_module(ident);
+ const auto iModulePhi = m_sctid->phi_module(ident);
int iModule = iModulePhi;
if (iModuleEta < 0) iModule +=4;
@@ -258,26 +258,31 @@ StatusCode TrackerAlignDBTool::createDB() const
auto iStation = m_sctid->station(ident);
auto iLayer = m_sctid->layer(ident);
- const auto buildKey = [](auto iStation, auto iLayer) {
- std::stringstream ss;
- ss << iStation << iLayer;
- return ss.str();
- };
-
- const auto key = buildKey(iStation, iLayer);
- if (not (m_alignment.find(key) == m_alignment.end())) {
- const std::vector<double> c = m_alignment.value().find(key)->second;
- ATH_MSG_VERBOSE("Applying correction for " << key);
- ATH_MSG_VERBOSE(c[0] << " " << c[1] << " " << c[2] << " " << c[3] << " " << c[4] << " " << c[5]);
- Amg::Translation3D newtranslation(c[0], c[1], c[2]);
- Amg::Transform3D alignment = newtranslation * Amg::RotationMatrix3D::Identity();
- alignment *= Amg::AngleAxis3D(c[5], Amg::Vector3D(0.,0.,1.));
- alignment *= Amg::AngleAxis3D(c[4], Amg::Vector3D(0.,1.,0.));
- alignment *= Amg::AngleAxis3D(c[3], Amg::Vector3D(1.,0.,0.));
-
- pat->add(ident2, Amg::EigenTransformToCLHEP(alignment));
+ const auto buildKey = [](auto iStation, auto iLayer) {
+ std::stringstream ss;
+ ss << iStation << iLayer;
+ return ss.str();
+ };
+
+ const auto key = buildKey(iStation, iLayer);
+ if (not (m_alignment.find(key) == m_alignment.end())) {
+ const std::vector<double> c = m_alignment.value().find(key)->second;
+ ATH_MSG_VERBOSE("Applying correction for " << key);
+ ATH_MSG_VERBOSE(c[0] << " " << c[1] << " " << c[2] << " " << c[3] << " " << c[4] << " " << c[5]);
+ Amg::Translation3D newtranslation(c[0], c[1], c[2]);
+ Amg::Transform3D alignment = newtranslation * Amg::RotationMatrix3D::Identity();
+ alignment *= Amg::AngleAxis3D(c[5], Amg::Vector3D(0.,0.,1.));
+ alignment *= Amg::AngleAxis3D(c[4], Amg::Vector3D(0.,1.,0.));
+ alignment *= Amg::AngleAxis3D(c[3], Amg::Vector3D(1.,0.,0.));
+
+ Amg::Translation3D newtransl(0,0,element->transform().translation()[2]);
+ Amg::Transform3D newcorr= newtransl * Amg::RotationMatrix3D::Identity();
+
+ pat->add(ident2, Amg::EigenTransformToCLHEP(newcorr*alignment*newcorr.inverse()));
+
+// pat->add(ident2, Amg::EigenTransformToCLHEP(alignment));
} else {
- ATH_MSG_VERBOSE("No correction given for " << key);
+ ATH_MSG_VERBOSE("No correction given for " << key);
}
}
else
@@ -294,12 +299,6 @@ StatusCode TrackerAlignDBTool::createDB() const
std::string key1 = dirkey(ident1, 1);
if ((pat = getTransPtr(key1)))
{
- // Amg::Translation3D translation(0.1,0.2,0.3);
- //Amg::Transform3D globshift=translation*Amg::RotationMatrix3D::Identity();
- //std::cout<<"rotation"<<std::endl;
- //std::cout<<globshift.rotation()(0,0)<<" , "<<globshift.rotation()(1,1)<<" , "<<globshift.rotation()(2,2)<<std::endl;
- //std::cout<<"translation"<<std::endl;
- //std::cout<<globshift.translation()(0,0)<<" , "<<globshift.translation()(1,1)<<" , "<<globshift.translation()(2,2)<<std::endl;
Amg::Transform3D globshift;
globshift.setIdentity();
for (int station : m_stations)
diff --git a/Tracker/TrackerDetDescr/FaserSCT_GeoModel/python/FaserSCT_GeoModelConfig.py b/Tracker/TrackerDetDescr/FaserSCT_GeoModel/python/FaserSCT_GeoModelConfig.py
index 038c8c38bb3c6d0daed962726f310867873a4c02..b9a3215b3696127b083a3cd6a71fa25e10f227be 100644
--- a/Tracker/TrackerDetDescr/FaserSCT_GeoModel/python/FaserSCT_GeoModelConfig.py
+++ b/Tracker/TrackerDetDescr/FaserSCT_GeoModel/python/FaserSCT_GeoModelConfig.py
@@ -2,7 +2,7 @@
from AthenaConfiguration.ComponentFactory import CompFactory
from AthenaConfiguration.Enums import ProductionStep
-from IOVDbSvc.IOVDbSvcConfig import addFolders, addFoldersSplitOnline
+from IOVDbSvc.IOVDbSvcConfig import addFolders #, addFoldersSplitOnline
def FaserSCT_GeometryCfg( flags ):
from FaserGeoModel.GeoModelConfig import GeoModelCfg
@@ -27,17 +27,20 @@ def FaserSCT_GeometryCfg( flags ):
if flags.GeoModel.Align.Disable:
return acc
+ # Don't specify db= in addFolders calls below
+ # dbname set globally to OFLP200 for MC and CONDBR3 for data
+ # Now have proper data and MC DBs for the alignment
if flags.GeoModel.Align.Dynamic:
# acc.merge(addFoldersSplitOnline(flags,"INDET","/Indet/Onl/AlignL1/ID","/Indet/AlignL1/ID",className="CondAttrListCollection"))
# acc.merge(addFoldersSplitOnline(flags,"INDET","/Indet/Onl/AlignL2/SCT","/Indet/AlignL2/SCT",className="CondAttrListCollection"))
# acc.merge(addFoldersSplitOnline(flags,"INDET","/Indet/Onl/AlignL3","/Indet/AlignL3",className="AlignableTransformContainer"))
print("FaserSCT dynamic align flag is not supported!")
else:
- print("Override Alignment dbname to OFLP200, fix this when alignment available in CONDBR3")
if flags.Common.Project != "AthSimulation" and (flags.Common.ProductionStep != ProductionStep.Simulation or flags.Overlay.DataOverlay):
- acc.merge(addFolders(flags,"/Tracker/Align", "SCT_OFL", className="AlignableTransformContainer", db="OFLP200"))
+ acc.merge(addFolders(flags,"/Tracker/Align", "SCT_OFL", className="AlignableTransformContainer"))
else:
- acc.merge(addFolders(flags, "/Tracker/Align", "SCT_OFL", db="OFLP200"))
+ acc.merge(addFolders(flags, "/Tracker/Align", "SCT_OFL"))
+
if flags.Common.Project != "AthSimulation": # Protection for AthSimulation builds
if flags.Common.ProductionStep != ProductionStep.Simulation or flags.Overlay.DataOverlay:
FaserSCT_AlignCondAlg = CompFactory.FaserSCT_AlignCondAlg
diff --git a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SCT_DetectorManager.cxx b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SCT_DetectorManager.cxx
index 33aa094fc5abee616cd7713d054c510eed394960..96acb3c5f306e472c35af834f597f0c2a74c0f86 100755
--- a/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SCT_DetectorManager.cxx
+++ b/Tracker/TrackerDetDescr/TrackerReadoutGeometry/src/SCT_DetectorManager.cxx
@@ -180,6 +180,7 @@ namespace TrackerDD {
GeoVAlignmentStore* alignStore) const
{
if (level == 0) { // 0 - At the element level
+ std::cout<<"like SCT_DetectorManager::setAlignableTransformDelta 183 level "<<level<<std::endl;
// We retrieve it via a hashId.
IdentifierHash idHash = m_idHelper->wafer_hash(id);
@@ -194,6 +195,7 @@ namespace TrackerDD {
SiDetectorElement * element = m_elementCollection[idHash];
if (!element) return false;
+ std::cout<<"like SCT_DetectorManager::setAlignableTransformDelta 198"<<std::endl;
// Its a local transform
//See header file for definition of m_isLogical
@@ -213,6 +215,7 @@ namespace TrackerDD {
} else if (level == 1) { // module level
+ std::cout<<"like SCT_DetectorManager::setAlignableTransformDelta 218"<<std::endl;
// We retrieve it via a hashId.
IdentifierHash idHash = m_idHelper->wafer_hash(id);
if (!idHash.is_valid()) return false;
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
index aecebc372e415c6c58f1464706c6f25ea6001dc2..a194e065b23c83dc4b345a73a50f1ec499a95916 100644
--- a/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
@@ -23,6 +23,8 @@ StatusCode MyExtrapolationExample::initialize() {
}
StatusCode MyExtrapolationExample::execute(const EventContext &ctx) const {
+ // The following uses the nominal geometry (without alignment)
+ // To use the aligned version, use getGeometryContext() instead
const Acts::GeometryContext gctx =
m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext().context();
diff --git a/Tracker/TrackerRecAlgs/TrackerSpacePointFormation/python/TrackerSpacePointFormationConfig.py b/Tracker/TrackerRecAlgs/TrackerSpacePointFormation/python/TrackerSpacePointFormationConfig.py
index 4bd7402d130957407bdf9d3ca3a660b18530060a..d98f18b66a4f4d5a6eb8fa461382b021b10de7b9 100644
--- a/Tracker/TrackerRecAlgs/TrackerSpacePointFormation/python/TrackerSpacePointFormationConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerSpacePointFormation/python/TrackerSpacePointFormationConfig.py
@@ -67,7 +67,8 @@ def TrackerSpacePointFinderCfg(flags, **kwargs):
"""Return ComponentAccumulator for SCT SpacePoints and Output"""
acc=TrackerDDSiElementPropertiesTableCondAlgCfg(flags)
acc.merge(TrackerSpacePointFinderBasicCfg(flags, **kwargs))
- acc.merge(TrackerSpacePointFinder_OutputCfg(flags))
+ if flags.Output.doWriteESD:
+ acc.merge(TrackerSpacePointFinder_OutputCfg(flags))
return acc
def StatisticsCfg(flags, **kwargs):
diff --git a/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py b/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py
index d81e6fd90ce1bd9128b715857840fd93c9fcf126..9b13d866e9ef5808d3953f5a7529683341566aa4 100644
--- a/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py
+++ b/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py
@@ -37,6 +37,7 @@ def ActsTrackingGeometryToolCfg(configFlags, name = "ActsTrackingGeometryTool" )
result = ComponentAccumulator()
acc = ActsTrackingGeometrySvcCfg(configFlags)
result.merge(acc)
+ result.merge(ActsAlignmentCondAlgCfg(configFlags))
Acts_ActsTrackingGeometryTool = FaserActsTrackingGeometryTool
actsTrackingGeometryTool = Acts_ActsTrackingGeometryTool("TrackingGeometryTool")
result.addPublicTool(actsTrackingGeometryTool)
@@ -61,7 +62,7 @@ def ActsAlignmentCondAlgCfg(configFlags, name = "ActsAlignmentCondAlg", **kwargs
acc = ActsTrackingGeometrySvcCfg(configFlags)
result.merge(acc)
- Acts_ActsAlignmentCondAlg = CompFactory.ActsAlignmentCondAlg
+ Acts_ActsAlignmentCondAlg = CompFactory.FaserActsAlignmentCondAlg
actsAlignmentCondAlg = Acts_ActsAlignmentCondAlg(name, **kwargs)
result.addCondAlgo(actsAlignmentCondAlg)
@@ -98,17 +99,6 @@ def ActsMaterialTrackWriterSvcCfg(name="FaserActsMaterialTrackWriterSvc",
result.addService(ActsMaterialTrackWriterSvc, primary=True)
return result
-def ActsMaterialStepConverterToolCfg(name = "ActsMaterialStepConverterTool" ) :
- result=ComponentAccumulator()
-
- Acts_ActsMaterialStepConverterTool = CompFactory.ActsMaterialStepConverterTool
- ActsMaterialStepConverterTool = Acts_ActsMaterialStepConverterTool(name)
-
- from AthenaCommon.Constants import INFO
- ActsMaterialStepConverterTool.OutputLevel = INFO
-
- result.addPublicTool(ActsMaterialStepConverterTool, primary=True)
- return result
def ActsSurfaceMappingToolCfg(configFlags, name = "FaserActsSurfaceMappingTool" ) :
result=ComponentAccumulator()
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsAlignmentCondAlg.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsAlignmentCondAlg.cxx
index 5f0735b99c7192d7d761191046157b42340e7cf9..3e38e6fd02c81a67a6f90450b96d074d680006c7 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsAlignmentCondAlg.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsAlignmentCondAlg.cxx
@@ -108,10 +108,12 @@ StatusCode FaserActsAlignmentCondAlg::execute() {
size_t nElems = 0;
trkGeom->visitSurfaces([&actsAlignStore, &nElems](const Acts::Surface *srf) {
const Acts::DetectorElementBase *detElem =
- srf->associatedDetectorElement();
+ srf->associatedDetectorElement();
+ if(detElem){
const auto *gmde = static_cast<const FaserActsDetectorElement *>(detElem);
gmde->storeTransform(actsAlignStore.get());
nElems++;
+ }
});
ATH_MSG_DEBUG("FaserActsAlignmentStore populated for " << nElems
<< " detector elements");
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx
index 5545334849da0bb6b068c768773ae4878e52b5f5..59357dfc70de1a9b9d025d2c2bbf2b5a7d0a421b 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx
@@ -118,9 +118,10 @@ FaserActsDetectorElement::transform(const Acts::GeometryContext& anygctx) const
void
FaserActsDetectorElement::storeTransform(FaserActsAlignmentStore* gas) const
{
- Amg::Transform3D g2l
- = m_detElement->getMaterialGeom()->getDefAbsoluteTransform()
- * Amg::CLHEPTransformToEigen(m_detElement->recoToHitTransform());
+// Amg::Transform3D g2l
+// = m_detElement->getMaterialGeom()->getDefAbsoluteTransform()
+// * Amg::CLHEPTransformToEigen(m_detElement->recoToHitTransform());
+ Amg::Transform3D g2l = m_detElement->transform() ;
// need to make sure translation has correct units
g2l.translation() *= length_unit;
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
index 0e92508b08f01b792c15913a45f37cda8c0eb5c0..d31fe8d3c11e804e115aef2cc6dfc08a2354b9e7 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
@@ -126,8 +126,7 @@ FaserActsExtrapolationTool::propagationSteps(const EventContext& ctx,
Acts::MagneticFieldContext mctx = getMagneticFieldContext(ctx);
const FaserActsGeometryContext& gctx
- //= m_trackingGeometryTool->getGeometryContext(ctx);
- = m_trackingGeometryTool->getNominalGeometryContext();
+ = m_trackingGeometryTool->getGeometryContext(ctx);
auto anygctx = gctx.context();
@@ -200,8 +199,7 @@ FaserActsExtrapolationTool::propagate(const EventContext& ctx,
Acts::MagneticFieldContext mctx = getMagneticFieldContext(ctx);
const FaserActsGeometryContext& gctx
- //= m_trackingGeometryTool->getGeometryContext(ctx);
- = m_trackingGeometryTool->getNominalGeometryContext();
+ = m_trackingGeometryTool->getGeometryContext(ctx);
auto anygctx = gctx.context();
@@ -319,8 +317,7 @@ FaserActsExtrapolationTool::propagate(const EventContext& ctx,
ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " begin");
Acts::MagneticFieldContext mctx = getMagneticFieldContext(ctx);
- const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext();
-// const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext(ctx);
+ const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext(ctx);
auto anygctx = gctx.context();
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
index a45c7c06d37214248a19b189ca40886339c9f1b4..aef43cf4022b059edf059d614d22929b84644871 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx
@@ -63,7 +63,7 @@ FaserActsTrackingGeometrySvc::initialize()
//get all the subdetectors
const GeoModelExperiment* theExpt = nullptr;
- ATH_CHECK( m_detStore->retrieve( theExpt ,"FASER"));
+ ATH_CHECK( m_detStore->retrieve( theExpt ));
const GeoVDetectorManager *vetoManager = theExpt->getManager("Veto");
const GeoVDetectorManager *triggerManager = theExpt->getManager("Trigger");
const GeoVDetectorManager *dipoleManager = theExpt->getManager("Dipole");
diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
index f93c094287a4966f84689abaa2615da2b2ef8882..874d604249142448e76634c8dac9a8d0ef1db8f5 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
@@ -29,6 +29,7 @@ atlas_add_component(FaserActsKalmanFilter
CircleFit.h
CircleFitTrackSeedTool.h
CKF2.h
+ CKF2Alignment.h
CombinatorialKalmanFilterAlg.h
EffPlotTool.h
FASERSourceLink.h
@@ -73,6 +74,7 @@ atlas_add_component(FaserActsKalmanFilter
src/CircleFit.cxx
src/CircleFitTrackSeedTool.cxx
src/CKF2.cxx
+ src/CKF2Alignment.cxx
src/CreateTrkTrackTool.h
src/CreateTrkTrackTool.cxx
# src/ClusterTrackSeedTool.cxx
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Alignment.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Alignment.py
new file mode 100644
index 0000000000000000000000000000000000000000..eb5f84da0eb61bf7031090cf49911b9728ea82d0
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Alignment.py
@@ -0,0 +1,145 @@
+#!/usr/bin/env python
+
+import sys
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from AthenaConfiguration.TestDefaults import defaultTestFiles
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
+from FaserActsKalmanFilter.GhostBustersConfig import GhostBustersCfg
+from TruthSeededTrackFinder.TruthSeededTrackFinderConfig import TruthSeededTrackFinderCfg
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsExtrapolationToolCfg
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+THistSvc=CompFactory.getComps("THistSvc")
+
+CKF2Alignment, FaserActsExtrapolationTool, FaserActsTrackingGeometryTool=CompFactory.getComps("CKF2Alignment", "FaserActsExtrapolationTool","FaserActsTrackingGeometryTool")
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+
+
+def CKF2AlignmentCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags)
+# acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
+ acc.merge(result)
+ track_seed_tool = CompFactory.CircleFitTrackSeedTool()
+ #remove the IFT in the track finding or not
+ track_seed_tool.removeIFT = False
+ sigma_loc0 = 1.9e-2
+ sigma_loc1 = 9e-1
+ sigma_phi = 3.3e-2
+ sigma_theta = 2.9e-4
+ p = 1000
+ sigma_p = 0.1 * p
+ sigma_qop = sigma_p / (p * p)
+ initial_variance_inflation = [100, 100, 100, 100, 1000]
+ track_seed_tool.covLoc0 = initial_variance_inflation[0] * sigma_loc1 * sigma_loc1
+ track_seed_tool.covLoc1 = initial_variance_inflation[1] * sigma_loc0 * sigma_loc0
+ track_seed_tool.covPhi = initial_variance_inflation[2] * sigma_phi * sigma_phi
+ track_seed_tool.covTheta = initial_variance_inflation[3] * sigma_theta * sigma_theta
+ track_seed_tool.covQOverP = initial_variance_inflation[4] * sigma_qop * sigma_qop
+ track_seed_tool.std_cluster = 0.0231
+ track_seed_tool.TrackCollection = "Segments"
+ # useless in the alignment, will clean them up
+ trajectory_states_writer_tool = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool.noDiagnostics = kwargs.pop("noDiagnostics", True)
+ trajectory_states_writer_tool1 = CompFactory.RootTrajectoryStatesWriterTool()
+ trajectory_states_writer_tool1.noDiagnostics = kwargs.pop("noDiagnostics", True)
+ trajectory_states_writer_tool1.FilePath = "/pool/track_states_ckf1_10root"
+
+ trajectory_summary_writer_tool = CompFactory.RootTrajectorySummaryWriterTool()
+ trajectory_summary_writer_tool.noDiagnostics = kwargs.pop("noDiagnostics", True)
+ trajectory_summary_writer_tool1 = CompFactory.RootTrajectorySummaryWriterTool()
+ trajectory_summary_writer_tool1.FilePath = "/pool/track_summary_ckf1_0.root"
+ trajectory_summary_writer_tool1.noDiagnostics = kwargs.pop("noDiagnostics", True)
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = actsTrackingGeometryTool
+ acc.merge(result)
+
+ trajectory_performance_writer_tool = CompFactory.PerformanceWriterTool("PerformanceWriterTool")
+ trajectory_performance_writer_tool.ExtrapolationTool = actsExtrapolationTool
+ trajectory_performance_writer_tool.noDiagnostics = kwargs.pop("noDiagnostics", True)
+ ckf2fit = CompFactory.CKF2Alignment(**kwargs)
+ #biased residual or unbiased
+ ckf2fit.BiasedResidual = False
+ ckf2fit.ExtrapolationTool = actsExtrapolationTool
+ ckf2fit.TrackingGeometryTool=actsTrackingGeometryTool
+ kalman_fitter1 = CompFactory.KalmanFitterTool(name="fitterTool1")
+ kalman_fitter1.noDiagnostics = True
+ kalman_fitter1.ActsLogging = "INFO"
+ kalman_fitter1.SummaryWriter = False
+ kalman_fitter1.StatesWriter = False
+ kalman_fitter1.SeedCovarianceScale = 10
+ kalman_fitter1.isMC = False
+ kalman_fitter1.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool1
+ kalman_fitter1.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool1
+ ckf2fit.KalmanFitterTool1 = kalman_fitter1
+ ckf2fit.TrackSeed = track_seed_tool
+ ckf2fit.ActsLogging = "INFO"
+ ckf2fit.isMC = False
+ ckf2fit.nMax = 10
+ ckf2fit.chi2Max = 100000
+ ckf2fit.maxSteps = 5000
+ histSvc= CompFactory.THistSvc()
+ # output ntuple
+ histSvc.Output += [ "CKF2Alignment DATAFILE='kfalignment_mc.root' OPT='RECREATE'"]
+ acc.addService(histSvc)
+ acc.addEventAlgo(ckf2fit)
+ return acc
+
+if __name__ == "__main__":
+
+# log.setLevel(INFO)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [ '' ] #input RDO
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER cosmic ray geometry (station 2 only)
+ ConfigFlags.TrackingGeometry.MaterialSource = "material-maps.json" # material map
+ ConfigFlags.Input.isMC = True
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+ ConfigFlags.addFlag("Input.InitialTimeStamp", 0)
+ ConfigFlags.Exec.MaxEvents = -1
+ ConfigFlags.Output.doWriteESD = False
+ #ConfigFlags.Concurrency.NumThreads = 1
+# ConfigFlags.addFlag("Alignment.Output", "ckf_alignment.root")
+ ConfigFlags.fillFromArgs(sys.argv[1:])
+ ConfigFlags.dump()
+ ConfigFlags.lock()
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(FaserGeometryCfg(ConfigFlags))
+ acc.merge(PoolReadCfg(ConfigFlags))
+ acc.merge(PoolWriteCfg(ConfigFlags))
+ from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+ # Inner Detector
+ acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+ acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+ acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.61, MinClustersPerFit=5, TanThetaXZCut=0.083))
+ acc.merge(GhostBustersCfg(ConfigFlags))
+ acc.merge(CKF2AlignmentCfg(ConfigFlags))
+ logging.getLogger('forcomps').setLevel(INFO)
+ acc.foreach_component("*").OutputLevel = INFO
+ acc.foreach_component("*ClassID*").OutputLevel = INFO
+
+ # Execute and finish
+ sc = acc.run()
+ sys.exit(not sc.isSuccess())
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
index 1035484b832fec4f32ca36400309a36cd250ed91..5d6ed5c4f3170cfc57b058d6bc9432329c1da561 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
@@ -8,6 +8,7 @@ from AthenaConfiguration.ComponentFactory import CompFactory
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometrySvcCfg
+from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg
def FaserActsAlignmentCondAlgCfg(flags, **kwargs):
@@ -72,7 +73,11 @@ def CKF2Cfg(flags, **kwargs):
actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
actsExtrapolationTool.MaxSteps = 1000
- actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ #use the tracking geometry tool which has enabled the alignment condition algo
+ result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags)
+ actsExtrapolationTool.TrackingGeometryTool = actsTrackingGeometryTool
+ acc.merge(result)
trajectory_performance_writer_tool = CompFactory.PerformanceWriterTool("PerformanceWriterTool")
trajectory_performance_writer_tool.ExtrapolationTool = actsExtrapolationTool
@@ -95,15 +100,14 @@ def CKF2Cfg(flags, **kwargs):
ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool
ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool
ckf.PerformanceWriterTool = trajectory_performance_writer_tool
+ ckf.TrackingGeometryTool=actsTrackingGeometryTool
ckf.isMC = flags.Input.isMC
ckf.SummaryWriter = True
ckf.StatesWriter = False
ckf.PerformanceWriter = False
ckf.nMax = 10
- # Use larger chi2 cut until alignment improves
- # ckf.chi2Max = 25
- ckf.chi2Max = 100000
+ ckf.chi2Max = 25
# Protect against running out of memory on busy events
ckf.maxSteps = 5000
acc.addEventAlgo(ckf)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
index ebe56d179393fcf3f51cb8df510271fbf5a85ebe..57d7c035d2fe56f3ec62a4a9c0013421c37b57ff 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/GhostBustersConfig.py
@@ -26,7 +26,8 @@ def GhostBustersCfg(flags, **kwargs):
acts_tracking_geometry_svc = ActsTrackingGeometrySvcCfg(flags)
acc.merge(acts_tracking_geometry_svc )
acc.addEventAlgo(CompFactory.GhostBusters(**kwargs))
- acc.merge(GhostBusters_OutputCfg(flags))
+ if flags.Output.doWriteESD:
+ acc.merge(GhostBusters_OutputCfg(flags))
# thistSvc = CompFactory.THistSvc()
# thistSvc.Output += ["HIST2 DATAFILE='GhostBusters.root' OPT='RECREATE'"]
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/ActsTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/ActsTrackSeedTool.cxx
index 4161c19eecbc6b30b0436732f76498f71f3f4fa6..817fcfea171e76d70d75f49a540e0604f2f52984 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/ActsTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/ActsTrackSeedTool.cxx
@@ -36,7 +36,7 @@ StatusCode ActsTrackSeedTool::run(std::vector<int> /*maskedLayers*/) {
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
- const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext();
+ const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext();
Acts::GeometryContext geoctx = gctx.context();
const int kSize = 1;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 858b97d6b2655859d77f027b31a57ce5782b95ae..f990be9458e0ee2ac0a5bdef8b670703f1424c26 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -2,6 +2,7 @@
#include "StoreGate/ReadHandle.h"
#include "StoreGate/ReadCondHandleKey.h"
+#include "StoreGate/WriteDecorHandle.h"
#include "TrackerSpacePoint/FaserSCT_SpacePointCollection.h"
#include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
#include "TrackerIdentifier/FaserSCT_ID.h"
@@ -28,9 +29,9 @@
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
-
-
#include "Acts/EventData/Measurement.hpp"
+#include "Acts/Propagator/PropagatorError.hpp"
+#include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp"
size_t CKF2::TrajectoryInfo::nClusters {0};
@@ -49,6 +50,7 @@ StatusCode CKF2::initialize() {
ATH_CHECK(m_kalmanFitterTool1.retrieve());
ATH_CHECK(m_createTrkTrackTool.retrieve());
ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_eventInfoKey.initialize());
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool.retrieve());
}
@@ -78,13 +80,15 @@ StatusCode CKF2::execute() {
const EventContext& ctx = Gaudi::Hive::currentContext();
m_numberOfEvents++;
+ SG::WriteDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx);
+
SG::WriteHandle trackContainer{m_trackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
- const FaserActsGeometryContext& faserActsGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
+ const FaserActsGeometryContext& faserActsGeometryContext = m_trackingGeometryTool->getGeometryContext();
auto gctx = faserActsGeometryContext.context();
Acts::MagneticFieldContext magFieldContext = getMagneticFieldContext(ctx);
Acts::CalibrationContext calibContext;
@@ -150,6 +154,14 @@ StatusCode CKF2::execute() {
std::list<TrajectoryInfo> allTrajectories;
for (auto &result : results) {
if (not result.ok()) {
+ // TODO use status bits for different errors
+ // result.error() == Acts::CombinatorialKalmanFilterError::NoTrackFound
+ if (result.error() == Acts::PropagatorError::StepCountLimitReached ||
+ result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) {
+ if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) {
+ ATH_MSG_WARNING (" cannot set error state for SCT");
+ }
+ }
continue;
}
CKFResult ckfResult = result.value();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 2111afcf6b52d46b45e1e2539d105290b937689e..170f0353904bef82095146c50bf98d5721ef2138 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -21,6 +21,7 @@
#include "KalmanFitterTool.h"
#include <boost/dynamic_bitset.hpp>
#include "CreateTrkTrackTool.h"
+#include "xAODEventInfo/EventInfo.h"
using ConstTrackStateProxy = Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>;
using ClusterSet = boost::dynamic_bitset<>;
@@ -148,6 +149,7 @@ private:
ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
Gaudi::Property<bool> m_isMC {this, "isMC", false};
SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "CKFTrackCollection", "Output track collection name" };
+ SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
};
#endif // FASERACTSKALMANFILTER_CKF2_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..afd377fc5effaa67f7c9e949bf6ad31c2bddd665
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.cxx
@@ -0,0 +1,1405 @@
+#include "CKF2Alignment.h"
+
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/ReadCondHandleKey.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePointCollection.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrkPrepRawData/PrepRawData.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrkRIO_OnTrack/RIO_OnTrack.h"
+#include "TrkSurfaces/Surface.h"
+#include "Identifier/Identifier.h"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "FaserActsKalmanFilter/IndexSourceLink.h"
+#include "FaserActsKalmanFilter/Measurement.h"
+#include "FaserActsRecMultiTrajectory.h"
+#include "TrackSelection.h"
+#include <algorithm>
+
+#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
+
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
+#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+
+
+#include "Acts/EventData/Measurement.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "TFile.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TMatrixDSym.h"
+#include "TMatrixD.h"
+#include <fstream>
+#include <iostream>
+
+size_t CKF2Alignment::TrajectoryInfo::nClusters {0};
+
+using TrajectoriesContainer = std::vector<FaserActsRecMultiTrajectory>;
+
+using namespace Acts::UnitLiterals;
+
+CKF2Alignment::CKF2Alignment(
+ const std::string& name, ISvcLocator* pSvcLocator)
+ : AthAlgorithm(name, pSvcLocator), m_thistSvc("THistSvc", name) {}
+
+
+ StatusCode CKF2Alignment::initialize() {
+ ATH_CHECK(m_fieldCondObjInputKey.initialize());
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ ATH_CHECK(m_trackSeedTool.retrieve());
+ ATH_CHECK(m_kalmanFitterTool1.retrieve());
+ ATH_CHECK(m_createTrkTrackTool.retrieve());
+ ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
+ m_fit = makeTrackFinderFunction(m_trackingGeometryTool->trackingGeometry(),
+ m_resolvePassive, m_resolveMaterial, m_resolveSensitive);
+ m_kf = makeTrackFitterFunction(m_trackingGeometryTool->trackingGeometry());
+ // FIXME fix Acts logging level
+ if (m_actsLogging == "VERBOSE") {
+ m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::VERBOSE);
+ } else if (m_actsLogging == "DEBUG") {
+ m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::DEBUG);
+ } else {
+ m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::INFO);
+ }
+ CHECK(m_thistSvc.retrieve());
+ m_tree= new TTree("trackParam","tree");
+ initializeTree();
+ CHECK(m_thistSvc->regTree("/CKF2Alignment/trackParam",m_tree));
+ m_tree->AutoSave();
+ return StatusCode::SUCCESS;
+ }
+
+
+StatusCode CKF2Alignment::execute() {
+ clearVariables();
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+ m_numberOfEvents++;
+
+ ATH_MSG_DEBUG("get the tracking geometry");
+ bool save=false;
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
+ = m_trackingGeometryTool->trackingGeometry();
+
+ const FaserActsGeometryContext& faserActsGeometryContext = m_trackingGeometryTool->getGeometryContext();
+ auto gctx = faserActsGeometryContext.context();
+ Acts::MagneticFieldContext magFieldContext = getMagneticFieldContext(ctx);
+ Acts::CalibrationContext calibContext;
+ ATH_MSG_DEBUG("get the seed");
+
+ CHECK(m_trackSeedTool->run());
+ std::shared_ptr<const Acts::Surface> initialSurface =
+ m_trackSeedTool->initialSurface();
+ std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>> initialParameters =
+ m_trackSeedTool->initialTrackParameters();
+ std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks =
+ m_trackSeedTool->sourceLinks();
+ double origin = m_trackSeedTool->targetZPosition();
+
+ ATH_MSG_DEBUG("get the output of seeding");
+ std::shared_ptr<std::vector<Measurement>> measurements = m_trackSeedTool->measurements();
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters = m_trackSeedTool->clusters();
+ std::vector<const Tracker::FaserSCT_Cluster*> clusters_sta0;
+ clusters_sta0.clear();
+ std::shared_ptr<std::vector<const Tracker::FaserSCT_SpacePoint*>> spacePoints = m_trackSeedTool->spacePoints();
+ std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters = m_trackSeedTool->seedClusters();
+
+ TrajectoryInfo::nClusters = sourceLinks->size();
+ TrajectoriesContainer trajectories;
+ trajectories.reserve(initialParameters->size());
+ if(sourceLinks->size()<15||sourceLinks->size()>40)return StatusCode::SUCCESS;
+ ATH_MSG_DEBUG("size of clusters/sourcelinks "<<clusters->size()<<"/"<<sourceLinks->size());
+
+ Acts::PropagatorPlainOptions pOptions;
+ pOptions.maxSteps = m_maxSteps;
+
+ Acts::MeasurementSelector::Config measurementSelectorCfg = {
+ {Acts::GeometryIdentifier(), {m_chi2Max, m_nMax}},
+ };
+
+
+ // Set the CombinatorialKalmanFilter options
+ CKF2Alignment::TrackFinderOptions options(
+ gctx, magFieldContext, calibContext,
+ IndexSourceLinkAccessor(), MeasurementCalibrator(*measurements),
+ Acts::MeasurementSelector(measurementSelectorCfg),
+ Acts::LoggerWrapper{*m_logger}, pOptions, &(*initialSurface));
+
+ // Perform the track finding for all initial parameters
+ m_numberOfTrackSeeds += initialParameters->size();
+ ATH_MSG_DEBUG("Invoke track finding with " << initialParameters->size() << " seeds.");
+ IndexSourceLinkContainer tmp;
+ for (const auto& sl : *sourceLinks) {
+ tmp.emplace_hint(tmp.end(), sl);
+ }
+
+ for (const auto& init : *initialParameters) {
+ ATH_MSG_DEBUG(" position: " << init.position(gctx).transpose());
+ ATH_MSG_DEBUG(" momentum: " << init.momentum().transpose());
+ ATH_MSG_DEBUG(" charge: " << init.charge());
+ }
+
+ ATH_MSG_DEBUG("sourcelink size "<<tmp.size());
+ auto results = (*m_fit)(tmp, *initialParameters, options);
+
+ // results contain a MultiTrajectory for each track seed with a trajectory of each branch of the CKF.
+ // To simplify the ambiguity solving a list of MultiTrajectories is created, each containing only a single track.
+ std::list<TrajectoryInfo> allTrajectories;
+ for (auto &result : results) {
+ if (not result.ok()) {
+ continue;
+ }
+ CKFResult ckfResult = result.value();
+ for (size_t trackTip : ckfResult.lastMeasurementIndices) {
+ allTrajectories.emplace_back(TrajectoryInfo(FaserActsRecMultiTrajectory(
+ ckfResult.fittedStates, {trackTip}, {{trackTip, ckfResult.fittedParameters.at(trackTip)}})));
+ }
+ }
+ m_numberOfFittedTracks += allTrajectories.size();
+
+ // the list of MultiTrajectories is sorted by the number of measurements using the chi2 value as a tie-breaker
+ allTrajectories.sort([](const TrajectoryInfo &left, const TrajectoryInfo &right) {
+ if (left.nMeasurements > right.nMeasurements) return true;
+ if (left.nMeasurements < right.nMeasurements) return false;
+ if (left.chi2 < right.chi2) return true;
+ else return false;
+ });
+
+ // select all tracks with at least 13 heats and with 6 or less shared hits, starting from the best track
+ // TODO use Gaudi parameters for the number of hits and shared hits
+ // TODO allow shared hits only in the first station?
+ std::vector<FaserActsRecMultiTrajectory> selectedTrajectories {};
+ while (not allTrajectories.empty()) {
+ TrajectoryInfo selected = allTrajectories.front();
+ selectedTrajectories.push_back(selected.trajectory);
+ allTrajectories.remove_if([&](const TrajectoryInfo &p) {
+ return (p.nMeasurements <= 12) || ((p.clusterSet & selected.clusterSet).count() > 6);
+ });
+ }
+
+
+ //only use the first track
+ int trackid=0;
+ for (const FaserActsRecMultiTrajectory &traj : selectedTrajectories) {
+ if(trackid>0)continue;
+ trackid++;
+ const auto& mj = traj.multiTrajectory();
+ const auto& trackTips = traj.tips();
+ if(trackTips.empty())continue;
+ auto& trackTip = trackTips.front();
+ auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
+ if(trajState.nMeasurements<15)continue;
+
+ const auto params = traj.trackParameters(traj.tips().front());
+ ATH_MSG_DEBUG("Fitted parameters");
+ ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ ATH_MSG_DEBUG(" charge: " << params.charge());
+ ATH_MSG_VERBOSE("checkte track");
+ std::vector<double> t_align_stationId_sp;
+ std::vector<double> t_align_centery_sp;
+ std::vector<double> t_align_layerId_sp;
+ std::vector<double> t_align_moduleId_sp;
+ std::vector<double> t_align_clustersize_sp;
+ std::vector<double> t_align_stereoId;
+ std::vector<double> t_align_global_residual_y_sp;
+ std::vector<double> t_align_global_measured_y_sp;
+ std::vector<double> t_align_global_measured_x_sp;
+ std::vector<double> t_align_global_measured_z_sp;
+ std::vector<double> t_align_global_measured_ye_sp;
+ std::vector<double> t_align_global_fitted_ye_sp;
+ std::vector<double> t_align_global_fitted_y_sp;
+ std::vector<double> t_align_global_fitted_x_sp;
+ std::vector<double> t_align_local_residual_x_sp;
+ std::vector<double> t_align_unbiased;
+ std::vector<double> t_align_local_pull_x_sp;
+ std::vector<double> t_align_local_measured_x_sp;
+ std::vector<double> t_align_local_measured_xe_sp;
+ std::vector<double> t_align_local_fitted_xe_sp;
+ std::vector<double> t_align_local_fitted_x_sp;
+ std::vector<double> t_align_derivation_x_x;
+ std::vector<double> t_align_derivation_x_y;
+ std::vector<double> t_align_derivation_x_z;
+ std::vector<double> t_align_derivation_x_rx;
+ std::vector<double> t_align_derivation_x_ry;
+ std::vector<double> t_align_derivation_x_rz;
+ std::vector<double> t_align_derivation_global_y_x;
+ std::vector<double> t_align_derivation_global_y_y;
+ std::vector<double> t_align_derivation_global_y_z;
+ std::vector<double> t_align_derivation_global_y_rx;
+ std::vector<double> t_align_derivation_global_y_ry;
+ std::vector<double> t_align_derivation_global_y_rz;
+ bool foundift=false;
+ std::unique_ptr<const Acts::BoundTrackParameters> ini_Param=nullptr;
+ //get residual from ACTS
+ if(m_biased){
+ mj.visitBackwards(trackTip, [&](const auto &state) {
+ /// Only fill the track states with non-outlier measurement
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag) and not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ return true;
+ }
+
+ const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ Acts::Vector3 mom(1, 1, 1);
+ Acts::Vector3 global = surface.localToGlobal(gctx, local, dir);
+ if (state.hasSmoothed())
+ {
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance());
+ auto covariance = state.smoothedCovariance();
+
+ /// Local hit residual info
+ auto H = state.effectiveProjector();
+ auto resCov = state.effectiveCalibratedCovariance() +
+ H * covariance * H.transpose();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+
+ if (state.hasUncalibrated()) {
+ const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
+ auto trans2 = cluster->detectorElement()->surface().transform();
+ auto trans1 = trans2;
+ Identifier id = cluster->identify();
+ int istation = m_idHelper->station(id);
+ int ilayer = m_idHelper->layer(id);
+ int stereoid = m_idHelper->side(id);
+ if(stereoid==1){
+ trans1=cluster->detectorElement()->otherSide()->surface().transform();
+ }
+ ATH_MSG_DEBUG(" save the residual "<<residual(Acts::eBoundLoc0));
+ const auto iModuleEta = m_idHelper->eta_module(id);
+ const auto iModulePhi = m_idHelper->phi_module(id);
+ int iModule = iModulePhi;
+ if (iModuleEta < 0) iModule +=4;
+ ATH_MSG_DEBUG("ID "<<istation<<"/"<<ilayer<<"/"<<iModule<<"/"<<stereoid<<"/"<<ilayer);
+ t_align_stationId_sp.push_back(istation);
+ t_align_centery_sp.push_back(cluster->detectorElement()->center().y());
+ t_align_layerId_sp.push_back(ilayer);
+ t_align_moduleId_sp.push_back(iModule);
+ t_align_clustersize_sp.push_back(cluster->rdoList().size());
+ t_align_stereoId.push_back(stereoid);
+ t_align_global_measured_x_sp.push_back(cluster->globalPosition().x());
+ t_align_global_measured_y_sp.push_back(cluster->globalPosition().y());
+ t_align_global_measured_z_sp.push_back(cluster->globalPosition().z());
+ t_align_global_fitted_x_sp.push_back(global.x());
+ t_align_global_fitted_y_sp.push_back(global.y());
+ t_align_local_residual_x_sp.push_back(residual(Acts::eBoundLoc0));
+ t_align_unbiased.push_back(0);
+ t_align_local_pull_x_sp.push_back(residual(Acts::eBoundLoc0)/sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ t_align_local_measured_x_sp.push_back(cluster->localPosition().x());
+ t_align_local_measured_xe_sp.push_back(sqrt(cluster->localCovariance()(0,0)));
+ t_align_local_fitted_xe_sp.push_back(sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+ t_align_local_fitted_x_sp.push_back(meas[Acts::eBoundLoc0]);
+ save=true;
+ Amg::Vector2D local_mea=cluster->localPosition() ;
+ t_align_derivation_x_x.push_back(getDerivation(ctx, parameter, trans1,local_mea, 0, stereoid, trans2, 0));
+ t_align_derivation_x_y.push_back(getDerivation(ctx, parameter, trans1,local_mea, 1, stereoid, trans2, 0));
+ t_align_derivation_x_z.push_back(getDerivation(ctx, parameter, trans1,local_mea, 2, stereoid, trans2, 0));
+ t_align_derivation_x_rx.push_back(getDerivation(ctx, parameter, trans1,local_mea, 3, stereoid, trans2, 0));
+ t_align_derivation_x_ry.push_back(getDerivation(ctx, parameter, trans1,local_mea, 4, stereoid, trans2, 0));
+ t_align_derivation_x_rz.push_back(getDerivation(ctx, parameter, trans1,local_mea, 5, stereoid, trans2, 0));
+ t_align_derivation_global_y_x.push_back(getDerivation(ctx, parameter, trans1,local_mea, 0, stereoid, trans2, 1));
+ t_align_derivation_global_y_y.push_back(getDerivation(ctx, parameter, trans1,local_mea, 1, stereoid, trans2, 1));
+ t_align_derivation_global_y_z.push_back(getDerivation(ctx, parameter, trans1,local_mea, 2, stereoid, trans2, 1));
+ t_align_derivation_global_y_rx.push_back(getDerivation(ctx, parameter, trans1,local_mea, 3, stereoid, trans2, 1));
+ t_align_derivation_global_y_ry.push_back(getDerivation(ctx, parameter, trans1,local_mea, 4, stereoid, trans2, 1));
+ t_align_derivation_global_y_rz.push_back(getDerivation(ctx, parameter, trans1,local_mea, 5, stereoid, trans2, 1));
+ ATH_MSG_VERBOSE("local derivation "<<t_align_derivation_x_x.back()<<" "<<t_align_derivation_x_y.back()<<" "<<t_align_derivation_x_z.back()<<" "<<t_align_derivation_x_rx.back()<<" "<<t_align_derivation_x_ry.back()<<" "<<t_align_derivation_x_rz.back());
+ ATH_MSG_VERBOSE("global derivation "<<t_align_derivation_global_y_x.back()<<" "<<t_align_derivation_global_y_y.back()<<" "<<t_align_derivation_global_y_z.back()<<""<<t_align_derivation_global_y_rx.back()<<" "<<t_align_derivation_global_y_ry.back()<<" "<<t_align_derivation_global_y_rz.back());
+ //get the residual for IFT
+ if(istation==1&&ilayer==0&&!foundift){
+ Amg::Transform3D ini_trans = Amg::Translation3D(0,0, -1910.) * Amg::RotationMatrix3D::Identity();
+ auto ini_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(ini_trans, std::make_shared<const Acts::RectangleBounds>(1000.,1000.));
+ ini_Param = m_extrapolationTool->propagate( ctx, parameter, *ini_surface, Acts::backward);
+ //std::unique_ptr<const Acts::BoundTrackParameters> ini_Param = m_extrapolationTool->propagate( ctx, parameter, *ini_surface, Acts::backward);
+ //find the closet cluster
+ if(ini_Param!=nullptr){
+ for(int i=0;i<3;i++){
+ for(int st=0;st<2;st++){
+ double chisq=999.;
+ size_t iclus=999;
+ double fitx=-999.;
+ double fity=-999.;
+ double gfitx=-999.;
+ double gfity=-999.;
+ double resi=-999.;
+ for(size_t ic =0 ;ic <clusters->size();ic++){
+ auto cluster = clusters->at(ic);
+ Identifier id = cluster->identify();
+ int jstation = m_idHelper->station(id);
+ int jlayer = m_idHelper->layer(id);
+ if(jstation!=0||jlayer!=i)continue;
+ int jstereoid = m_idHelper->side(id);
+ if(st==0){if(jstereoid==1)continue;}
+ if(st==1){if(jstereoid!=1)continue;}
+ auto trans2 = cluster->detectorElement()->surface().transform();
+ Amg::Transform3D shift_trans = Amg::Translation3D(0,0, cluster->globalPosition().z()) * Amg::RotationMatrix3D::Identity();
+ auto shift_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(shift_trans, std::make_shared<const Acts::RectangleBounds>(1000.,1000.));
+ auto shift_Param = m_extrapolationTool->propagate( ctx, *ini_Param, *shift_surface);
+ if(shift_Param!=nullptr){
+ auto shift_parameter = shift_Param->parameters();
+ Amg::Vector2D local_m(shift_parameter[Acts::eBoundLoc0],shift_parameter[Acts::eBoundLoc1]);
+ Amg::Vector3D glopos(local_m.x(),local_m.y(),cluster->globalPosition().z());
+ auto localpos=trans2.inverse()*glopos;
+ ATH_MSG_DEBUG("looping cluster "<<local_m.x()<<" "<<local_m.y()<<" : "<<cluster->globalPosition());
+ ATH_MSG_DEBUG("local cluster "<<localpos.x()<<" "<<localpos.y()<<" : "<<cluster->localPosition());
+ double dist2 =(glopos.x()-cluster->globalPosition().x())*(glopos.x()-cluster->globalPosition().x())/40./40.+(glopos.y()-cluster->globalPosition().y())*(glopos.y()-cluster->globalPosition().y());
+ if(dist2<chisq){
+ chisq=dist2;
+ iclus=ic;
+ fitx=localpos.x();
+ fity=localpos.y();
+ gfitx=glopos.x();
+ gfity=glopos.y();
+ resi=cluster->localPosition().x() - localpos.x();
+ }
+ }
+ }
+ if(iclus>200)continue;
+ clusters_sta0.push_back(clusters->at(iclus));
+ foundift=true;
+ ATH_MSG_DEBUG(" local position "<<fitx<<" "<<fity<<" : "<<clusters->at(iclus)->localPosition());
+ ATH_MSG_DEBUG(" global position "<<gfitx<<" "<<gfity<<" : "<<clusters->at(iclus)->globalPosition());
+ Amg::Vector3D localresi(clusters->at(iclus)->localPosition().x()-fitx,clusters->at(iclus)->localPosition().y()-fity,0.);
+ auto gloresi=clusters->at(iclus)->detectorElement()->surface().transform()*localresi;
+ ATH_MSG_DEBUG(" local/global residual "<<localresi.x()<<" "<<localresi.y()<<" : "<<gloresi.x()<<" "<<gloresi.y()<<" "<<gloresi.z());
+
+ Identifier id = clusters->at(iclus)->identify();
+ int istation = m_idHelper->station(id);
+ int ilayer = m_idHelper->layer(id);
+ if(istation!=0||ilayer!=i)continue;
+ int stereoid = m_idHelper->side(id);
+ auto trans2 = clusters->at(iclus)->detectorElement()->surface().transform();
+ auto trans1 = trans2;
+ if(stereoid==1){
+ trans1=clusters->at(iclus)->detectorElement()->otherSide()->surface().transform();
+ }
+ ATH_MSG_DEBUG(" save the residual "<<resi);
+ const auto iModuleEta = m_idHelper->eta_module(id);
+ const auto iModulePhi = m_idHelper->phi_module(id);
+ int iModule = iModulePhi;
+ if (iModuleEta < 0) iModule +=4;
+ ATH_MSG_DEBUG("ID "<<istation<<"/"<<ilayer<<"/"<<iModule<<"/"<<stereoid<<"/"<<ilayer);
+ t_align_stationId_sp.push_back(istation);
+ t_align_centery_sp.push_back(clusters->at(iclus)->detectorElement()->center().y());
+ t_align_layerId_sp.push_back(ilayer);
+ t_align_moduleId_sp.push_back(iModule);
+ t_align_clustersize_sp.push_back(clusters->at(iclus)->rdoList().size());
+ t_align_stereoId.push_back(stereoid);
+ t_align_global_measured_x_sp.push_back(clusters->at(iclus)->globalPosition().x());
+ t_align_global_measured_y_sp.push_back(clusters->at(iclus)->globalPosition().y());
+ t_align_global_measured_z_sp.push_back(clusters->at(iclus)->globalPosition().z());
+ t_align_global_fitted_x_sp.push_back(gfitx);
+ t_align_global_fitted_y_sp.push_back(gfity);
+ t_align_local_residual_x_sp.push_back(resi);
+ t_align_unbiased.push_back(2);
+ t_align_local_pull_x_sp.push_back(resi/sqrt(clusters->at(iclus)->localCovariance()(0,0)));
+ t_align_local_measured_x_sp.push_back(clusters->at(iclus)->localPosition().x());
+ t_align_local_measured_xe_sp.push_back(sqrt(cluster->localCovariance()(0,0)));
+ t_align_local_fitted_xe_sp.push_back(sqrt(cluster->localCovariance()(0,0)));
+ t_align_local_fitted_x_sp.push_back(fitx);
+ Amg::Vector2D local_mea=clusters->at(iclus)->localPosition() ;
+ t_align_derivation_x_x.push_back(getDerivation(ctx, parameter, trans1,local_mea, 0, stereoid, trans2, 0));
+ t_align_derivation_x_y.push_back(getDerivation(ctx, parameter, trans1,local_mea, 1, stereoid, trans2, 0));
+ t_align_derivation_x_z.push_back(getDerivation(ctx, parameter, trans1,local_mea, 2, stereoid, trans2, 0));
+ t_align_derivation_x_rx.push_back(getDerivation(ctx, parameter, trans1,local_mea, 3, stereoid, trans2, 0));
+ t_align_derivation_x_ry.push_back(getDerivation(ctx, parameter, trans1,local_mea, 4, stereoid, trans2, 0));
+ t_align_derivation_x_rz.push_back(getDerivation(ctx, parameter, trans1,local_mea, 5, stereoid, trans2, 0));
+ t_align_derivation_global_y_x.push_back(getDerivation(ctx, parameter, trans1,local_mea, 0, stereoid, trans2, 1));
+ t_align_derivation_global_y_y.push_back(getDerivation(ctx, parameter, trans1,local_mea, 1, stereoid, trans2, 1));
+ t_align_derivation_global_y_z.push_back(getDerivation(ctx, parameter, trans1,local_mea, 2, stereoid, trans2, 1));
+ t_align_derivation_global_y_rx.push_back(getDerivation(ctx, parameter, trans1,local_mea, 3, stereoid, trans2, 1));
+ t_align_derivation_global_y_ry.push_back(getDerivation(ctx, parameter, trans1,local_mea, 4, stereoid, trans2, 1));
+ t_align_derivation_global_y_rz.push_back(getDerivation(ctx, parameter, trans1,local_mea, 5, stereoid, trans2, 1));
+ ATH_MSG_VERBOSE("local derivation "<<t_align_derivation_x_x.back()<<" "<<t_align_derivation_x_y.back()<<" "<<t_align_derivation_x_z.back()<<" "<<t_align_derivation_x_rx.back()<<" "<<t_align_derivation_x_ry.back()<<" "<<t_align_derivation_x_rz.back());
+ ATH_MSG_VERBOSE("global derivation "<<t_align_derivation_global_y_x.back()<<" "<<t_align_derivation_global_y_y.back()<<" "<<t_align_derivation_global_y_z.back()<<""<<t_align_derivation_global_y_rx.back()<<" "<<t_align_derivation_global_y_ry.back()<<" "<<t_align_derivation_global_y_rz.back());
+ }
+ }
+ }
+ }
+ }
+ }
+ return true;
+ });
+ ATH_MSG_VERBOSE("check track");
+
+ std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj);
+ if (track==nullptr||ini_Param==nullptr) continue;
+ origin = -1910.;
+ std::vector<TSOS4Residual> unbiase_resi= m_kalmanFitterTool1->getUnbiasedResidual(ctx, gctx, track.get(), Acts::BoundVector::Zero(), m_isMC, origin,clusters_sta0, *ini_Param);
+ ATH_MSG_DEBUG(" found tsos size "<<unbiase_resi.size());
+ if(unbiase_resi.size()>0){
+
+ ATH_MSG_DEBUG(" check the tsos");
+ for (auto& tsos : unbiase_resi){
+ ATH_MSG_DEBUG(" check the param");
+ auto cluster = tsos.cluster;
+ ATH_MSG_DEBUG("Positions "<<tsos.fit_global_x<<" "<<tsos.fit_global_y<<" "<<tsos.fit_global_z<<" , "<<cluster->globalPosition().x()<<" "<<cluster->globalPosition().y()<<" "<<cluster->globalPosition().z());
+ ATH_MSG_DEBUG("Found the ctsos ");
+ double localx = tsos.fit_local_x;
+ double localy = tsos.fit_local_y;
+ double globalx= tsos.fit_global_x;
+ double globaly= tsos.fit_global_y;
+ double globalz= tsos.fit_global_z;
+ //transfer to cluster surface
+ Amg::Vector3D glo_pos(globalx,globaly,globalz);
+ std::optional<Amg::Vector2D> loc_pos = cluster->detectorElement()->surface().globalToLocal(glo_pos);
+
+ ATH_MSG_DEBUG("local pos "<<localx<<" "<<localy<<" , "<<loc_pos->x()<<" "<<loc_pos->y()<<" :vs "<<cluster->localPosition().x()<<" "<<cluster->localPosition().y());
+ localx = loc_pos->x();
+ localy = loc_pos->y();
+ ATH_MSG_DEBUG("global pos "<<globalx<<" "<<globaly<<" vs "<<cluster->globalPosition().x()<<" "<<cluster->globalPosition().y());
+ // outputTracks->push_back(std::move(newtrack));
+ ATH_MSG_DEBUG("Global pos "<<globalx<<" , "<<globaly);
+ auto trans2 = cluster->detectorElement()->surface().transform();
+ auto trans1 = trans2;
+ Identifier id = cluster->identify();
+ int istation = m_idHelper->station(id);
+ int ilayer = m_idHelper->layer(id);
+ if(istation!=0)continue;
+ //int layerid = istation*3 + ilayer;
+ // if(layerid == 0 || layerid == 11)continue;
+ int stereoid = m_idHelper->side(id);
+ if(stereoid==1){
+ trans1=cluster->detectorElement()->otherSide()->surface().transform();
+ }
+ Acts::BoundTrackParameters fitParameter=tsos.parameter;
+ const auto iModuleEta = m_idHelper->eta_module(id);
+ const auto iModulePhi = m_idHelper->phi_module(id);
+ int iModule = iModulePhi;
+ if (iModuleEta < 0) iModule +=4;
+ ATH_MSG_DEBUG("ID "<<istation<<"/"<<ilayer<<"/"<<iModule<<"/"<<stereoid<<"/"<<ilayer);
+ t_align_stationId_sp.push_back(istation);
+ t_align_centery_sp.push_back(cluster->detectorElement()->center().y());
+ t_align_layerId_sp.push_back(ilayer);
+ t_align_moduleId_sp.push_back(iModule);
+ t_align_clustersize_sp.push_back(cluster->rdoList().size());
+ t_align_stereoId.push_back(stereoid);
+ Amg::Vector2D local_mea=cluster->localPosition() ;
+ Amg::Vector2D local_meae=cluster->localPosition() + Amg::Vector2D(sqrt(cluster->localCovariance()(0,0)),0);
+ Amg::Vector2D local_fit(localx, localy);
+ Amg::Vector2D local_fite(localx+0.0001, localy+0.0001);
+ ATH_MSG_DEBUG(" save the residual "<<local_mea.x()-local_fit.x()<<" vs "<<tsos.residual);
+ t_align_local_residual_x_sp.push_back(tsos.residual);
+ t_align_unbiased.push_back(1);
+ t_align_local_pull_x_sp.push_back((local_mea.x()-local_fit.x())/(sqrt(((cluster->localCovariance()(0,0)+(local_fite.x()-local_fit.x())*(local_fite.x()-local_fit.x()))))));
+ t_align_local_measured_x_sp.push_back(local_mea.x());
+ t_align_local_measured_xe_sp.push_back(fabs(local_meae.x()-local_mea.x()));
+ t_align_local_fitted_xe_sp.push_back(fabs(local_fit.x()-local_fite.x()));
+ t_align_local_fitted_x_sp.push_back(local_fit.x());
+
+ Amg::Vector3D global_fit= cluster->detectorElement()->surface().localToGlobal(local_fit);
+ Amg::Vector3D global_fite= cluster->detectorElement()->surface().localToGlobal( local_fite);
+ Amg::Vector3D global_mea= cluster->detectorElement()->surface().localToGlobal( local_mea);
+ Amg::Vector3D global_meae= cluster->detectorElement()->surface().localToGlobal( local_meae);
+ t_align_global_residual_y_sp.push_back(global_mea.y()-global_fit.y());
+ t_align_global_measured_x_sp.push_back(global_mea.x());
+ t_align_global_measured_y_sp.push_back(global_mea.y());
+ t_align_global_measured_z_sp.push_back(global_mea.z());
+ t_align_global_measured_ye_sp.push_back(fabs(global_mea.y()-global_meae.y()));
+ t_align_global_fitted_ye_sp.push_back(fabs(global_fite.y()-global_fite.y()));
+ t_align_global_fitted_x_sp.push_back(global_fit.x());
+ t_align_global_fitted_y_sp.push_back(global_fit.y());
+ t_align_derivation_x_x.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 0, stereoid, trans2, 0));
+ t_align_derivation_x_y.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 1, stereoid, trans2, 0));
+ t_align_derivation_x_z.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 2, stereoid, trans2, 0));
+ t_align_derivation_x_rx.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 3, stereoid, trans2, 0));
+ t_align_derivation_x_ry.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 4, stereoid, trans2, 0));
+ t_align_derivation_x_rz.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 5, stereoid, trans2, 0));
+ t_align_derivation_global_y_x.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 0, stereoid, trans2, 1));
+ t_align_derivation_global_y_y.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 1, stereoid, trans2, 1));
+ t_align_derivation_global_y_z.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 2, stereoid, trans2, 1));
+ t_align_derivation_global_y_rx.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 3, stereoid, trans2, 1));
+ t_align_derivation_global_y_ry.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 4, stereoid, trans2, 1));
+ t_align_derivation_global_y_rz.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 5, stereoid, trans2, 1));
+ ATH_MSG_VERBOSE("local derivation "<<t_align_derivation_x_x.back()<<" "<<t_align_derivation_x_y.back()<<" "<<t_align_derivation_x_z.back()<<" "<<t_align_derivation_x_rx.back()<<" "<<t_align_derivation_x_ry.back()<<" "<<t_align_derivation_x_rz.back());
+ ATH_MSG_VERBOSE("global derivation "<<t_align_derivation_global_y_x.back()<<" "<<t_align_derivation_global_y_y.back()<<" "<<t_align_derivation_global_y_z.back()<<""<<t_align_derivation_global_y_rx.back()<<" "<<t_align_derivation_global_y_ry.back()<<" "<<t_align_derivation_global_y_rz.back());
+ }
+ }
+ }
+ else{
+ std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj);
+ if (track!=nullptr) {
+ for (const Trk::MeasurementBase *meas : *track->measurementsOnTrack()) {
+ const auto* clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(meas);
+ if (clusterOnTrack) {
+ const Tracker::FaserSCT_Cluster* cluster = clusterOnTrack->prepRawData();
+ ATH_MSG_DEBUG(" Found the cluster to be remove "<<cluster->globalPosition());
+ //remove the cluster and re-fit to get unbiased residual
+ Identifier id = clusterOnTrack->identify();
+ Identifier waferId = m_idHelper->wafer_id(id);
+ ATH_MSG_DEBUG("like cluster ID "<<m_idHelper->layer(id)<<" "<<m_idHelper->station(id));
+ double localx(-999);
+ double localxe(999);
+ double localy(-999);
+ double localye(999);
+ double globalx(-999);
+ double globaly(999);
+ double globalz(999);
+// int charge =0;
+// Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Identity();
+ ATH_MSG_DEBUG(" make unbiased residual");
+ m_numberOfSelectedTracks++;
+ if((*track->measurementsOnTrack()).size()<15)continue;
+ //remove the cluster and re-fit to get unbiased residual
+ std::vector<TSOS4Residual> unbiase_resi= m_kalmanFitterTool1->getUnbiasedResidual(ctx, gctx, track.get(), cluster->globalPosition().z(), Acts::BoundVector::Zero(), m_isMC, origin);
+ if(unbiase_resi.size()>0){
+
+ ATH_MSG_DEBUG(" check the tsos");
+ for (auto& tsos : unbiase_resi){
+ ATH_MSG_DEBUG(" check the param");
+ if(tsos.cluster->identify()!=cluster->identify())continue;
+ ATH_MSG_DEBUG("Positions "<<tsos.fit_global_x<<" "<<tsos.fit_global_y<<" "<<tsos.fit_global_z<<" , "<<cluster->globalPosition().x()<<" "<<cluster->globalPosition().y()<<" "<<cluster->globalPosition().z());
+ ATH_MSG_DEBUG("Found the ctsos ");
+ localx = tsos.fit_local_x;
+ localy = tsos.fit_local_y;
+ globalx= tsos.fit_global_x;
+ globaly= tsos.fit_global_y;
+ globalz= tsos.fit_global_z;
+ //transfer to cluster surface
+ Amg::Vector3D glo_pos(globalx,globaly,globalz);
+ std::optional<Amg::Vector2D> loc_pos = cluster->detectorElement()->surface().globalToLocal(glo_pos);
+
+ ATH_MSG_DEBUG("local pos "<<localx<<" "<<localy<<" , "<<loc_pos->x()<<" "<<loc_pos->y()<<" :vs "<<cluster->localPosition().x()<<" "<<cluster->localPosition().y());
+ localx = loc_pos->x();
+ localy = loc_pos->y();
+ ATH_MSG_DEBUG("global pos "<<globalx<<" "<<globaly<<" vs "<<cluster->globalPosition().x()<<" "<<cluster->globalPosition().y());
+ save=true;
+ ATH_MSG_DEBUG("Global pos "<<globalx<<" , "<<globaly);
+ auto trans2 = cluster->detectorElement()->surface().transform();
+ auto trans1 = trans2;
+ Identifier id = cluster->identify();
+ int istation = m_idHelper->station(id);
+ int ilayer = m_idHelper->layer(id);
+ int stereoid = m_idHelper->side(id);
+ if(stereoid==1){
+ trans1=cluster->detectorElement()->otherSide()->surface().transform();
+ }
+ Acts::BoundTrackParameters fitParameter=tsos.parameter;
+ const auto iModuleEta = m_idHelper->eta_module(id);
+ const auto iModulePhi = m_idHelper->phi_module(id);
+ int iModule = iModulePhi;
+ if (iModuleEta < 0) iModule +=4;
+ ATH_MSG_DEBUG("ID "<<istation<<"/"<<ilayer<<"/"<<iModule<<"/"<<stereoid<<"/"<<ilayer);
+ t_align_centery_sp.push_back(cluster->detectorElement()->center().y());
+ t_align_stationId_sp.push_back(istation);
+ t_align_layerId_sp.push_back(ilayer);
+ t_align_moduleId_sp.push_back(iModule);
+ t_align_clustersize_sp.push_back(cluster->rdoList().size());
+ t_align_stereoId.push_back(stereoid);
+ Amg::Vector2D local_mea=cluster->localPosition() ;
+ Amg::Vector2D local_meae=cluster->localPosition() + Amg::Vector2D(sqrt(cluster->localCovariance()(0,0)),0);
+ Amg::Vector2D local_fit(localx, localy);
+ Amg::Vector2D local_fite(localx+localxe, localy+localye);
+ ATH_MSG_DEBUG(" save the residual "<<local_mea.x()-local_fit.x()<<" vs "<<tsos.residual);
+ t_align_local_residual_x_sp.push_back(tsos.residual);
+ t_align_unbiased.push_back(3);
+ t_align_local_pull_x_sp.push_back((local_mea.x()-local_fit.x())/(sqrt(((cluster->localCovariance()(0,0)+(local_fite.x()-local_fit.x())*(local_fite.x()-local_fit.x()))))));
+ t_align_local_measured_x_sp.push_back(local_mea.x());
+ t_align_local_measured_xe_sp.push_back(fabs(local_meae.x()-local_mea.x()));
+ t_align_local_fitted_xe_sp.push_back(fabs(local_fit.x()-local_fite.x()));
+ t_align_local_fitted_x_sp.push_back(local_fit.x());
+
+ Amg::Vector3D global_fit= cluster->detectorElement()->surface().localToGlobal(local_fit);
+ Amg::Vector3D global_fite= cluster->detectorElement()->surface().localToGlobal( local_fite);
+ Amg::Vector3D global_mea= cluster->detectorElement()->surface().localToGlobal( local_mea);
+ Amg::Vector3D global_meae= cluster->detectorElement()->surface().localToGlobal( local_meae);
+ t_align_global_residual_y_sp.push_back(global_mea.y()-global_fit.y());
+ t_align_global_measured_x_sp.push_back(global_mea.x());
+ t_align_global_measured_y_sp.push_back(global_mea.y());
+ t_align_global_measured_z_sp.push_back(global_mea.z());
+ t_align_global_measured_ye_sp.push_back(fabs(global_mea.y()-global_meae.y()));
+ t_align_global_fitted_ye_sp.push_back(fabs(global_fite.y()-global_fite.y()));
+ t_align_global_fitted_x_sp.push_back(global_fit.x());
+ t_align_global_fitted_y_sp.push_back(global_fit.y());
+ t_align_derivation_x_x.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 0, stereoid, trans2, 0));
+ t_align_derivation_x_y.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 1, stereoid, trans2, 0));
+ t_align_derivation_x_z.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 2, stereoid, trans2, 0));
+ t_align_derivation_x_rx.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 3, stereoid, trans2, 0));
+ t_align_derivation_x_ry.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 4, stereoid, trans2, 0));
+ t_align_derivation_x_rz.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 5, stereoid, trans2, 0));
+ t_align_derivation_global_y_x.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 0, stereoid, trans2, 1));
+ t_align_derivation_global_y_y.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 1, stereoid, trans2, 1));
+ t_align_derivation_global_y_z.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 2, stereoid, trans2, 1));
+ t_align_derivation_global_y_rx.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 3, stereoid, trans2, 1));
+ t_align_derivation_global_y_ry.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 4, stereoid, trans2, 1));
+ t_align_derivation_global_y_rz.push_back(getDerivation(ctx, fitParameter, trans1,local_mea, 5, stereoid, trans2, 1));
+ ATH_MSG_VERBOSE("local derivation "<<t_align_derivation_x_x.back()<<" "<<t_align_derivation_x_y.back()<<" "<<t_align_derivation_x_z.back()<<" "<<t_align_derivation_x_rx.back()<<" "<<t_align_derivation_x_ry.back()<<" "<<t_align_derivation_x_rz.back());
+ ATH_MSG_VERBOSE("global derivation "<<t_align_derivation_global_y_x.back()<<" "<<t_align_derivation_global_y_y.back()<<" "<<t_align_derivation_global_y_z.back()<<""<<t_align_derivation_global_y_rx.back()<<" "<<t_align_derivation_global_y_ry.back()<<" "<<t_align_derivation_global_y_rz.back());
+ }
+ }
+ }
+ }
+ }
+ }
+ ATH_MSG_VERBOSE("finish getting residuals");
+ if(save){
+ ATH_MSG_VERBOSE("save the results");
+ m_fitParam_nMeasurements.push_back(trajState.nMeasurements);
+ m_fitParam_nStates.push_back(trajState.nStates);
+ m_fitParam_nOutliers.push_back(trajState.nOutliers);
+ m_fitParam_nHoles.push_back(trajState.nHoles);
+ m_fitParam_chi2.push_back(trajState.chi2Sum);
+ m_fitParam_ndf.push_back(trajState.NDF);
+ m_fitParam_x.push_back(params.position(gctx).transpose().x());
+ m_fitParam_y.push_back(params.position(gctx).transpose().y());
+ m_fitParam_z.push_back(params.position(gctx).transpose().z());
+ m_fitParam_charge.push_back(params.charge());
+ m_fitParam_px.push_back(params.momentum().transpose().x());
+ m_fitParam_py.push_back(params.momentum().transpose().y());
+ m_fitParam_pz.push_back(params.momentum().transpose().z());
+ m_align_stationId_sp.push_back(t_align_stationId_sp);
+ m_align_centery_sp.push_back(t_align_centery_sp);
+ m_align_layerId_sp.push_back(t_align_layerId_sp);
+ m_align_moduleId_sp.push_back(t_align_moduleId_sp);
+ m_align_clustersize_sp.push_back(t_align_clustersize_sp);
+ m_align_stereoId.push_back(t_align_stereoId);
+ m_align_global_residual_y_sp.push_back(t_align_global_residual_y_sp);
+ m_align_global_measured_x_sp.push_back(t_align_global_measured_x_sp);
+ m_align_global_measured_y_sp.push_back(t_align_global_measured_y_sp);
+ m_align_global_measured_z_sp.push_back(t_align_global_measured_z_sp);
+ m_align_global_measured_ye_sp.push_back(t_align_global_measured_ye_sp);
+ m_align_global_fitted_ye_sp.push_back(t_align_global_fitted_ye_sp);
+ m_align_global_fitted_x_sp.push_back(t_align_global_fitted_x_sp);
+ m_align_global_fitted_y_sp.push_back(t_align_global_fitted_y_sp);
+ m_align_local_residual_x_sp.push_back(t_align_local_residual_x_sp);
+ m_align_unbiased_sp.push_back(t_align_unbiased);
+ m_align_local_pull_x_sp.push_back(t_align_local_pull_x_sp);
+ m_align_local_measured_x_sp.push_back(t_align_local_measured_x_sp);
+ m_align_local_measured_xe_sp.push_back(t_align_local_measured_xe_sp);
+ m_align_local_fitted_xe_sp.push_back(t_align_local_fitted_xe_sp);
+ m_align_local_fitted_x_sp.push_back(t_align_local_fitted_x_sp);
+ m_align_derivation_x_x.push_back(t_align_derivation_x_x);
+ m_align_derivation_x_y.push_back(t_align_derivation_x_y);
+ m_align_derivation_x_z.push_back(t_align_derivation_x_z);
+ m_align_derivation_x_rx.push_back(t_align_derivation_x_rx);
+ m_align_derivation_x_ry.push_back(t_align_derivation_x_ry);
+ m_align_derivation_x_rz.push_back(t_align_derivation_x_rz);
+ m_align_derivation_global_y_x.push_back(t_align_derivation_global_y_x);
+ m_align_derivation_global_y_y.push_back(t_align_derivation_global_y_y);
+ m_align_derivation_global_y_z.push_back(t_align_derivation_global_y_z);
+ m_align_derivation_global_y_rx.push_back(t_align_derivation_global_y_rx);
+ m_align_derivation_global_y_ry.push_back(t_align_derivation_global_y_ry);
+ m_align_derivation_global_y_rz.push_back(t_align_derivation_global_y_rz);
+ }
+ }
+
+ if(save)
+ m_tree->Fill();
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode CKF2Alignment::finalize() {
+ ATH_MSG_INFO("CombinatorialKalmanFilterAlg::finalize()");
+ ATH_MSG_INFO(m_numberOfEvents << " events processed.");
+ ATH_MSG_INFO(m_numberOfTrackSeeds << " seeds.");
+ ATH_MSG_INFO(m_numberOfFittedTracks << " fitted tracks.");
+ ATH_MSG_INFO(m_numberOfSelectedTracks << " selected and re-fitted tracks.");
+ //nominal geometry
+ double nominaly[4][3][8]={0.};
+ double sigma=0.064;
+ nominaly[1][0][4]= -95.61;
+ nominaly[1][0][0]=-95.61;
+ nominaly[1][0][5]= -31.87;
+ nominaly[1][0][1]=-31.87;
+ nominaly[1][0][6]= 31.873;
+ nominaly[1][0][2]=31.877;
+ nominaly[1][0][7]= 95.61;
+ nominaly[1][0][3]=95.61;
+ nominaly[1][1][4]= -100.61;
+ nominaly[1][1][0]=-100.61;
+ nominaly[1][1][5]= -36.87;
+ nominaly[1][1][1]=-36.87;
+ nominaly[1][1][6]= 26.87;
+ nominaly[1][1][2]=26.87;
+ nominaly[1][1][7]= 90.61;
+ nominaly[1][1][3]=90.61;
+ nominaly[1][2][4]= -90.61;
+ nominaly[1][2][0]=-90.61;
+ nominaly[1][2][5]= -26.87;
+ nominaly[1][2][1]=-26.87;
+ nominaly[1][2][6]= 36.87;
+ nominaly[1][2][2]=36.87;
+ nominaly[1][2][7]= 100.61;
+ nominaly[1][2][3]=100.61;
+ nominaly[2][0][4]= -95.61;
+ nominaly[2][0][0]=-95.61;
+ nominaly[2][0][5]= -31.87;
+ nominaly[2][0][1]=-31.87;
+ nominaly[2][0][6]= 31.87;
+ nominaly[2][0][2]=31.87;
+ nominaly[2][0][7]= 95.61;
+ nominaly[2][0][3]=95.61;
+ nominaly[2][1][4]= -100.61;
+ nominaly[2][1][0]=-100.61;
+ nominaly[2][1][5]= -36.87;
+ nominaly[2][1][1]=-36.87;
+ nominaly[2][1][6]= 26.87;
+ nominaly[2][1][2]=26.87;
+ nominaly[2][1][7]= 90.61;
+ nominaly[2][1][3]=90.61;
+ nominaly[2][2][4]= -90.61;
+ nominaly[2][2][0]=-90.61;
+ nominaly[2][2][5]= -26.87;
+ nominaly[2][2][1]=-26.87;
+ nominaly[2][2][6]= 36.87;
+ nominaly[2][2][2]=36.87;
+ nominaly[2][2][7]= 100.61;
+ nominaly[2][2][3]=100.61;
+ nominaly[3][0][4]= -95.61;
+ nominaly[3][0][0]=-95.61;
+ nominaly[3][0][5]= -31.87;
+ nominaly[3][0][1]=-31.87;
+ nominaly[3][0][6]= 31.87;
+ nominaly[3][0][2]=31.87;
+ nominaly[3][0][7]= 95.61;
+ nominaly[3][0][3]=95.61;
+ nominaly[3][1][4]= -100.61;
+ nominaly[3][1][0]=-100.61;
+ nominaly[3][1][5]= -36.87;
+ nominaly[3][1][1]=-36.87;
+ nominaly[3][1][6]= 26.87;
+ nominaly[3][1][2]=26.87;
+ nominaly[3][1][7]= 90.61;
+ nominaly[3][1][3]=90.61;
+ nominaly[3][2][4]= -90.61;
+ nominaly[3][2][0]=-90.61;
+ nominaly[3][2][5]= -26.87;
+ nominaly[3][2][1]=-26.87;
+ nominaly[3][2][6]= 36.87;
+ nominaly[3][2][2]=36.87;
+ nominaly[3][2][7]=100.61;
+ nominaly[3][2][3]=100.61;
+ //dump the matrix to a txt file for alignment
+ std::vector<double>* t_fitParam_chi2=0;
+ std::vector<double>* t_fitParam_pz=0;
+ std::vector<double>* t_fitParam_x=0;
+ std::vector<double>* t_fitParam_y=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_residual_x_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_measured_x_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_measured_xe_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_stationId_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_centery_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_layerId_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_moduleId_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_unbiased_sp=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_x=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_y=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_z=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_rx=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_ry=0;
+ std::vector<std::vector<double>>* t_fitParam_align_local_derivation_x_rz=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_x=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_y=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_z=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_rx=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_ry=0;
+ std::vector<std::vector<double>>* t_fitParam_align_global_derivation_y_rz=0;
+ std::vector<std::vector<double>>* t_fitParam_align_stereoId=0;
+ std::vector<int>* t_fitParam_nMeasurements=0;
+ m_tree->SetBranchAddress("fitParam_nMeasurements", &t_fitParam_nMeasurements);
+ m_tree->SetBranchAddress("fitParam_align_stereoId",&t_fitParam_align_stereoId);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_x",&t_fitParam_align_local_derivation_x_x);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_y",&t_fitParam_align_local_derivation_x_y);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_z",&t_fitParam_align_local_derivation_x_z);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_rx",&t_fitParam_align_local_derivation_x_rx);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_ry",&t_fitParam_align_local_derivation_x_ry);
+ m_tree->SetBranchAddress("fitParam_align_local_derivation_x_rz",&t_fitParam_align_local_derivation_x_rz);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_x",&t_fitParam_align_global_derivation_y_x);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_y",&t_fitParam_align_global_derivation_y_y);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_z",&t_fitParam_align_global_derivation_y_z);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_rx",&t_fitParam_align_global_derivation_y_rx);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_ry",&t_fitParam_align_global_derivation_y_ry);
+ m_tree->SetBranchAddress("fitParam_align_global_derivation_y_rz",&t_fitParam_align_global_derivation_y_rz);
+ m_tree->SetBranchAddress("fitParam_align_local_residual_x_sp",&t_fitParam_align_local_residual_x_sp);
+ m_tree->SetBranchAddress("fitParam_chi2",&t_fitParam_chi2);
+ m_tree->SetBranchAddress("fitParam_pz",&t_fitParam_pz);
+ m_tree->SetBranchAddress("fitParam_x",&t_fitParam_x);
+ m_tree->SetBranchAddress("fitParam_y",&t_fitParam_y);
+ m_tree->SetBranchAddress("fitParam_align_local_measured_x_sp",&t_fitParam_align_local_measured_x_sp);
+ m_tree->SetBranchAddress("fitParam_align_local_measured_xe_sp",&t_fitParam_align_local_measured_xe_sp);
+ m_tree->SetBranchAddress("fitParam_align_layerId_sp",&t_fitParam_align_layerId_sp);
+ m_tree->SetBranchAddress("fitParam_align_moduleId_sp",&t_fitParam_align_moduleId_sp);
+ m_tree->SetBranchAddress("fitParam_align_unbiased_sp",&t_fitParam_align_unbiased_sp);
+ m_tree->SetBranchAddress("fitParam_align_stationId_sp",&t_fitParam_align_stationId_sp);
+ m_tree->SetBranchAddress("fitParam_align_centery_sp",&t_fitParam_align_centery_sp);
+ int ntrk_mod[12][8]={0};
+ int ntrk_sta[4]={0};
+ int ntrk_lay[12]={0};
+ //define matrix
+ int Ndim=6;
+ std::vector<TMatrixD> denom_lay;
+ std::vector<TMatrixD> num_lay;
+ std::vector<std::vector<TMatrixD>> denom_mod;
+ std::vector<std::vector<TMatrixD>> num_mod;
+ std::vector<std::vector<TMatrixD>> denom_mod1;
+ std::vector<std::vector<TMatrixD>> num_mod1;
+ //initialize to 0
+ auto num_matrix = [](int n) {
+ TMatrixD num_t(n,1);
+ for(int i=0;i<n;i++){
+ num_t[i][0]=0.;
+ }
+ return num_t;
+ };
+ auto denom_matrix = [](int n) {
+ TMatrixD denom_t(n,n);
+ for(int i=0;i<n;i++){
+ for(int j=0;j<n;j++){
+ denom_t[i][j]=0.;
+ }
+ }
+ return denom_t;
+ };
+ std::vector<TMatrixD> denom_sta;
+ std::vector<TMatrixD> num_sta;
+ for(int i=0;i<4;i++){
+ denom_sta.push_back(denom_matrix(Ndim));
+ num_sta.push_back(num_matrix(Ndim));
+ }
+ for(int i=0;i<12;i++){
+ denom_lay.push_back(denom_matrix(Ndim));
+ num_lay.push_back(num_matrix(Ndim));
+ std::vector<TMatrixD> denom_l;
+ std::vector<TMatrixD> num_l;
+ std::vector<TMatrixD> denom_l1;
+ std::vector<TMatrixD> num_l1;
+ for(int j=0;j<8;j++){
+ denom_l.push_back(denom_matrix(Ndim));
+ num_l.push_back(num_matrix(Ndim));
+ denom_l1.push_back(denom_matrix(Ndim));
+ num_l1.push_back(num_matrix(Ndim));
+ }
+ denom_mod.push_back(denom_l1);
+ denom_mod1.push_back(denom_l1);
+ num_mod.push_back(num_l);
+ num_mod1.push_back(num_l1);
+ }
+ double chi_mod_y[12][8]={0};
+ std::cout<<"found "<<m_tree->GetEntries()<<" events "<<std::endl;
+ //loop over all the entries
+ for(int ievt=0;ievt<m_tree->GetEntries();ievt++){
+ m_tree->GetEntry(ievt);
+ if(ievt%10000==0)std::cout<<"processing "<<ievt<<" event"<<std::endl;
+ if(t_fitParam_chi2->size()<1)continue;
+ for(int i=0;i<1;i+=1){
+ if(t_fitParam_chi2->at(i)>100||t_fitParam_pz->at(i)<300||sqrt(t_fitParam_x->at(i)*t_fitParam_x->at(i)+t_fitParam_y->at(i)*t_fitParam_y->at(i))>95)continue;
+ if(t_fitParam_align_local_residual_x_sp->at(i).size()<15)continue;//only use good track
+ for(size_t j=0;j<t_fitParam_align_local_residual_x_sp->at(i).size();j++){
+ double resx1=999.;
+ double x1e=999.;
+ if(m_biased&&t_fitParam_align_stationId_sp->at(i).at(j)==0){
+ if(t_fitParam_align_unbiased_sp->at(i).at(j)==1){
+ resx1=(t_fitParam_align_local_residual_x_sp->at(i).at(j));
+ x1e=sqrt(t_fitParam_align_local_measured_xe_sp->at(i).at(j));
+ if(fabs(resx1)>1.0)continue;
+ }
+ else continue;
+ }
+ else{
+ resx1=t_fitParam_align_local_residual_x_sp->at(i).at(j);
+ x1e=sqrt(t_fitParam_align_local_measured_xe_sp->at(i).at(j));
+ if(fabs(resx1)>0.2)continue;
+ }
+ TMatrixD m1(Ndim,1);
+ m1[0][0]=t_fitParam_align_local_derivation_x_x->at(i).at(j)/x1e;
+ m1[1][0]=t_fitParam_align_local_derivation_x_y->at(i).at(j)/x1e;
+ m1[2][0]=t_fitParam_align_local_derivation_x_z->at(i).at(j)/x1e;
+ m1[3][0]=t_fitParam_align_local_derivation_x_rx->at(i).at(j)/x1e;
+ m1[4][0]=t_fitParam_align_local_derivation_x_ry->at(i).at(j)/x1e;
+ m1[5][0]=t_fitParam_align_local_derivation_x_rz->at(i).at(j)/x1e;
+ TMatrixD m2(1,Ndim);
+ TMatrixD m3(Ndim,Ndim);
+ m2.Transpose(m1);
+ m3=m1*m2;
+ TMatrixD m4(Ndim,1);
+ m4[0][0]=t_fitParam_align_local_derivation_x_x->at(i).at(j)/x1e*resx1/x1e;
+ m4[1][0]=t_fitParam_align_local_derivation_x_y->at(i).at(j)/x1e*resx1/x1e;
+ m4[2][0]=t_fitParam_align_local_derivation_x_z->at(i).at(j)/x1e*resx1/x1e;
+ m4[3][0]=t_fitParam_align_local_derivation_x_rx->at(i).at(j)/x1e*resx1/x1e;
+ m4[4][0]=t_fitParam_align_local_derivation_x_ry->at(i).at(j)/x1e*resx1/x1e;
+ m4[5][0]=t_fitParam_align_local_derivation_x_rz->at(i).at(j)/x1e*resx1/x1e;
+ TMatrixD m6(Ndim,1);
+ m6=m4;
+ TMatrixD mg1(Ndim,1);
+ mg1[0][0]=t_fitParam_align_global_derivation_y_x->at(i).at(j)/x1e;
+ mg1[1][0]=t_fitParam_align_global_derivation_y_y->at(i).at(j)/x1e;
+ mg1[2][0]=t_fitParam_align_global_derivation_y_z->at(i).at(j)/x1e;
+ mg1[3][0]=t_fitParam_align_global_derivation_y_rx->at(i).at(j)/x1e;
+ mg1[4][0]=t_fitParam_align_global_derivation_y_ry->at(i).at(j)/x1e;
+ mg1[5][0]=t_fitParam_align_global_derivation_y_rz->at(i).at(j)/x1e;
+ TMatrixD mg2(1,Ndim);
+ TMatrixD mg3(Ndim,Ndim);
+ mg2.Transpose(mg1);
+ mg3=mg1*mg2;
+ TMatrixD mg4(Ndim,1);
+ mg4[0][0]=t_fitParam_align_global_derivation_y_x->at(i).at(j)/x1e*resx1/x1e;
+ mg4[1][0]=t_fitParam_align_global_derivation_y_y->at(i).at(j)/x1e*resx1/x1e;
+ mg4[2][0]=t_fitParam_align_global_derivation_y_z->at(i).at(j)/x1e*resx1/x1e;
+ mg4[3][0]=t_fitParam_align_global_derivation_y_rx->at(i).at(j)/x1e*resx1/x1e;
+ mg4[4][0]=t_fitParam_align_global_derivation_y_ry->at(i).at(j)/x1e*resx1/x1e;
+ mg4[5][0]=t_fitParam_align_global_derivation_y_rz->at(i).at(j)/x1e*resx1/x1e;
+ TMatrixD mg6(Ndim,1);
+ mg6=mg4;
+ int istation=t_fitParam_align_stationId_sp->at(i).at(j);
+ int ilayer=t_fitParam_align_layerId_sp->at(i).at(j);
+ int imodule=t_fitParam_align_moduleId_sp->at(i).at(j);
+ //station: 1,2,3
+ if(istation>=0&&istation<4){
+ denom_sta[istation]+=mg3;
+ num_sta[istation]+=mg6;
+ ntrk_sta[istation]+=1;
+ }
+ int ilayer_tot = ilayer + (istation )*3;
+ double deltay=t_fitParam_align_centery_sp->at(i).at(j)-nominaly[istation][ilayer][imodule];
+ ATH_MSG_DEBUG("nominal/new geometry Y position = "<<nominaly[istation][ilayer][imodule]<<"/"<<t_fitParam_align_centery_sp->at(i).at(j));
+ chi_mod_y[ilayer_tot][imodule]+=deltay/sigma/sigma;
+ if(ilayer_tot>=0&&ilayer_tot<12){
+ denom_lay[ilayer_tot]+=mg3;
+ num_lay[ilayer_tot]+=mg6;
+ ntrk_lay[ilayer_tot]+=1;
+ denom_mod[ilayer_tot][imodule]+=m3;
+ num_mod[ilayer_tot][imodule]+=m6;
+ ntrk_mod[ilayer_tot][imodule]+=1;
+ }
+ }
+ }
+ }
+
+ std::cout<<"Get the nominal transform"<<std::endl;
+ //get the alignment constants
+ std::ofstream align_output("alignment_matrix.txt",std::ios::out);
+ for(int i=0;i<4;i++){
+ for(int ii=0;ii<Ndim;ii++){
+ align_output<<i<<" "<<ii<<" "<<denom_sta[i][ii][0]<<" "<<denom_sta[i][ii][1]<<" "<<denom_sta[i][ii][2]<<" "<<denom_sta[i][ii][3]<<" "<<denom_sta[i][ii][4]<<" "<<denom_sta[i][ii][5]<<" "<<0<<std::endl;
+ if(ii==0)
+ std::cout<<i<<" "<<ii<<" "<<denom_sta[i][ii][0]<<" "<<denom_sta[i][ii][1]<<" "<<denom_sta[i][ii][2]<<" "<<denom_sta[i][ii][3]<<" "<<denom_sta[i][ii][4]<<" "<<denom_sta[i][ii][5]<<" "<<0<<std::endl;
+ }
+ align_output<<i<<" "<<6<<" "<<num_sta[i][0][0]<<" "<<num_sta[i][1][0]<<" "<<num_sta[i][2][0]<<" "<<num_sta[i][3][0]<<" "<<num_sta[i][4][0]<<" "<<num_sta[i][5][0]<<" "<<ntrk_sta[i]<<std::endl;
+ }
+ std::cout<<"Get the deltas for stations"<<std::endl;
+ //layers
+ for(size_t i=0;i<denom_lay.size();i++){
+ for(int ii=0;ii<Ndim;ii++){
+ align_output<<i/3<<i%3<<" "<<ii<<" "<<denom_lay[i][ii][0]<<" "<<denom_lay[i][ii][1]<<" "<<denom_lay[i][ii][2]<<" "<<denom_lay[i][ii][3]<<" "<<denom_lay[i][ii][4]<<" "<<denom_lay[i][ii][5]<<" "<<0<<std::endl;
+ }
+ align_output<<i/3<<i%3<<" "<<6<<" "<<num_lay[i][0][0]<<" "<<num_lay[i][1][0]<<" "<<num_lay[i][2][0]<<" "<<num_lay[i][3][0]<<" "<<num_lay[i][4][0]<<" "<<num_lay[i][5][0]<<" "<<ntrk_lay[i]<<std::endl;
+ }
+
+ std::cout<<"Get the deltas for layers"<<std::endl;
+ for(size_t i=0;i<denom_lay.size();i++){
+ for( size_t j=0;j<num_mod[i].size();j++){
+ for(int ii=0;ii<Ndim;ii++){
+ align_output<<i/3<<i%3<<j<<" "<<ii<<" "<<denom_mod[i][j][ii][0]<<" "<<denom_mod[i][j][ii][1]<<" "<<denom_mod[i][j][ii][2]<<" "<<denom_mod[i][j][ii][3]<<" "<<denom_mod[i][j][ii][4]<<" "<<denom_mod[i][j][ii][5]<<" "<<chi_mod_y[i][j]<<std::endl;
+ }
+ align_output<<i/3<<i%3<<j<<" "<<6<<" "<<num_mod[i][j][0][0]<<" "<<num_mod[i][j][1][0]<<" "<<num_mod[i][j][2][0]<<" "<<num_mod[i][j][3][0]<<" "<<num_mod[i][j][4][0]<<" "<<num_mod[i][j][5][0]<<" "<<ntrk_mod[i][j]<<std::endl;
+ }
+ }
+ std::cout<<"closing the output"<<std::endl;
+ align_output.close();
+ return StatusCode::SUCCESS;
+}
+
+
+Acts::MagneticFieldContext CKF2Alignment::getMagneticFieldContext(const EventContext& ctx) const {
+ SG::ReadCondHandle<FaserFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
+ if (!readHandle.isValid()) {
+ std::stringstream msg;
+ msg << "Failed to retrieve magnetic field condition data " << m_fieldCondObjInputKey.key() << ".";
+ throw std::runtime_error(msg.str());
+ }
+ const FaserFieldCacheCondObj* fieldCondObj{*readHandle};
+ return Acts::MagneticFieldContext(fieldCondObj);
+}
+
+
+
+
+void CKF2Alignment::computeSharedHits(std::vector<IndexSourceLink>* sourceLinks, TrackFinderResult& results) const {
+ // Compute shared hits from all the reconstructed tracks
+ // Compute nSharedhits and Update ckf results
+ // hit index -> list of multi traj indexes [traj, meas]
+ static_assert(Acts::SourceLinkConcept<IndexSourceLink>,
+ "Source link does not fulfill SourceLinkConcept");
+
+ std::vector<std::size_t> firstTrackOnTheHit(
+ sourceLinks->size(), std::numeric_limits<std::size_t>::max());
+ std::vector<std::size_t> firstStateOnTheHit(
+ sourceLinks->size(), std::numeric_limits<std::size_t>::max());
+
+ for (unsigned int iresult = 0; iresult < results.size(); iresult++) {
+ if (not results.at(iresult).ok()) {
+ continue;
+ }
+
+ auto& ckfResult = results.at(iresult).value();
+ auto& measIndexes = ckfResult.lastMeasurementIndices;
+
+ for (auto measIndex : measIndexes) {
+ ckfResult.fittedStates.visitBackwards(measIndex, [&](const auto& state) {
+ if (not state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag))
+ return;
+
+ std::size_t hitIndex = state.uncalibrated().index();
+
+ // Check if hit not already used
+ if (firstTrackOnTheHit.at(hitIndex) ==
+ std::numeric_limits<std::size_t>::max()) {
+ firstTrackOnTheHit.at(hitIndex) = iresult;
+ firstStateOnTheHit.at(hitIndex) = state.index();
+ return;
+ }
+
+ // if already used, control if first track state has been marked
+ // as shared
+ int indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
+ int indexFirstState = firstStateOnTheHit.at(hitIndex);
+ if (not results.at(indexFirstTrack).value().fittedStates.getTrackState(indexFirstState).typeFlags().test(Acts::TrackStateFlag::SharedHitFlag))
+ results.at(indexFirstTrack).value().fittedStates.getTrackState(indexFirstState).typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
+
+ // Decorate this track
+ results.at(iresult).value().fittedStates.getTrackState(state.index()).typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
+ });
+ }
+ }
+}
+
+
+namespace {
+
+ using Updater = Acts::GainMatrixUpdater;
+ using Smoother = Acts::GainMatrixSmoother;
+
+ using Stepper = Acts::EigenStepper<>;
+ using Navigator = Acts::Navigator;
+ using Propagator = Acts::Propagator<Stepper, Navigator>;
+ using CKF = Acts::CombinatorialKalmanFilter<Propagator, Updater, Smoother>;
+
+ using TrackParametersContainer = std::vector<CKF2Alignment::TrackParameters>;
+
+ struct TrackFinderFunctionImpl
+ : public CKF2Alignment::TrackFinderFunction {
+ CKF trackFinder;
+
+ TrackFinderFunctionImpl(CKF&& f) : trackFinder(std::move(f)) {}
+
+ CKF2Alignment::TrackFinderResult operator()(
+ const IndexSourceLinkContainer& sourcelinks,
+ const TrackParametersContainer& initialParameters,
+ const CKF2Alignment::TrackFinderOptions& options)
+ const override {
+ return trackFinder.findTracks(sourcelinks, initialParameters, options);
+ };
+ };
+
+} // namespace
+
+std::shared_ptr<CKF2Alignment::TrackFinderFunction>
+CKF2Alignment::makeTrackFinderFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ bool resolvePassive, bool resolveMaterial, bool resolveSensitive) {
+ auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
+ Stepper stepper(std::move(magneticField));
+ Navigator::Config cfg{trackingGeometry};
+ cfg.resolvePassive = resolvePassive;
+ cfg.resolveMaterial = resolveMaterial;
+ cfg.resolveSensitive = resolveSensitive;
+ Navigator navigator(cfg);
+ Propagator propagator(std::move(stepper), std::move(navigator));
+ CKF trackFinder(std::move(propagator));
+
+ // build the track finder functions. owns the track finder object.
+ return std::make_shared<TrackFinderFunctionImpl>(std::move(trackFinder));
+}
+
+
+namespace {
+
+ using Updater = Acts::GainMatrixUpdater;
+ using Smoother = Acts::GainMatrixSmoother;
+ using Stepper = Acts::EigenStepper<>;
+ using Propagator = Acts::Propagator<Stepper, Acts::Navigator>;
+ using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>;
+
+ struct TrackFitterFunctionImpl
+ : public CKF2Alignment::TrackFitterFunction {
+ Fitter trackFitter;
+
+ TrackFitterFunctionImpl(Fitter &&f) : trackFitter(std::move(f)) {}
+
+ CKF2Alignment::KFResult operator()(
+ const std::vector<IndexSourceLink> &sourceLinks,
+ const Acts::BoundTrackParameters &initialParameters,
+ const CKF2Alignment::TrackFitterOptions &options)
+ const override {
+ return trackFitter.fit(sourceLinks, initialParameters, options);
+ };
+ };
+
+} // namespace
+
+
+std::shared_ptr<CKF2Alignment::TrackFitterFunction>
+CKF2Alignment::makeTrackFitterFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry) {
+ auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
+ auto stepper = Stepper(std::move(magneticField));
+ Acts::Navigator::Config cfg{trackingGeometry};
+ cfg.resolvePassive = false;
+ cfg.resolveMaterial = true;
+ cfg.resolveSensitive = true;
+ Acts::Navigator navigator(cfg);
+ Propagator propagator(std::move(stepper), std::move(navigator));
+ Fitter trackFitter(std::move(propagator));
+ return std::make_shared<TrackFitterFunctionImpl>(std::move(trackFitter));
+}
+
+
+void CKF2Alignment::initializeTree(){
+ m_tree->Branch("evtId",&m_numberOfEvents);
+ m_tree->Branch("fitParam_x", &m_fitParam_x);
+ m_tree->Branch("fitParam_charge", &m_fitParam_charge);
+ m_tree->Branch("fitParam_y", &m_fitParam_y);
+ m_tree->Branch("fitParam_z", &m_fitParam_z);
+ m_tree->Branch("fitParam_px", &m_fitParam_px);
+ m_tree->Branch("fitParam_py", &m_fitParam_py);
+ m_tree->Branch("fitParam_pz", &m_fitParam_pz);
+ m_tree->Branch("fitParam_chi2", &m_fitParam_chi2);
+ m_tree->Branch("fitParam_ndf", &m_fitParam_ndf);
+ m_tree->Branch("fitParam_nHoles", &m_fitParam_nHoles);
+ m_tree->Branch("fitParam_nOutliers", &m_fitParam_nOutliers);
+ m_tree->Branch("fitParam_nStates", &m_fitParam_nStates);
+ m_tree->Branch("fitParam_nMeasurements", &m_fitParam_nMeasurements);
+ m_tree->Branch("fitParam_align_stereoId", &m_align_stereoId);
+ m_tree->Branch("fitParam_align_stationId_sp", &m_align_stationId_sp);
+ m_tree->Branch("fitParam_align_centery_sp", &m_align_centery_sp);
+ m_tree->Branch("fitParam_align_layerId_sp", &m_align_layerId_sp);
+ m_tree->Branch("fitParam_align_moduleId_sp", &m_align_moduleId_sp);
+ m_tree->Branch("fitParam_align_clustersize_sp", &m_align_clustersize_sp);
+ m_tree->Branch("fitParam_align_local_derivation_x_x", &m_align_derivation_x_x);
+ m_tree->Branch("fitParam_align_local_derivation_x_y", &m_align_derivation_x_y);
+ m_tree->Branch("fitParam_align_local_derivation_x_z", &m_align_derivation_x_z);
+ m_tree->Branch("fitParam_align_local_derivation_x_rx", &m_align_derivation_x_rx);
+ m_tree->Branch("fitParam_align_local_derivation_x_ry", &m_align_derivation_x_ry);
+ m_tree->Branch("fitParam_align_local_derivation_x_rz", &m_align_derivation_x_rz);
+ m_tree->Branch("fitParam_align_global_derivation_y_x", &m_align_derivation_global_y_x);
+ m_tree->Branch("fitParam_align_global_derivation_y_y", &m_align_derivation_global_y_y);
+ m_tree->Branch("fitParam_align_global_derivation_y_z", &m_align_derivation_global_y_z);
+ m_tree->Branch("fitParam_align_global_derivation_y_rx", &m_align_derivation_global_y_rx);
+ m_tree->Branch("fitParam_align_global_derivation_y_ry", &m_align_derivation_global_y_ry);
+ m_tree->Branch("fitParam_align_global_derivation_y_rz", &m_align_derivation_global_y_rz);
+ m_tree->Branch("fitParam_align_local_residual_x_sp", &m_align_local_residual_x_sp);
+ m_tree->Branch("fitParam_align_unbiased_sp", &m_align_unbiased_sp);
+ m_tree->Branch("fitParam_align_local_pull_x_sp", &m_align_local_pull_x_sp);
+ m_tree->Branch("fitParam_align_local_measured_x_sp", &m_align_local_measured_x_sp);
+ m_tree->Branch("fitParam_align_local_measured_xe_sp", &m_align_local_measured_xe_sp);
+ m_tree->Branch("fitParam_align_local_fitted_xe_sp", &m_align_local_fitted_xe_sp);
+ m_tree->Branch("fitParam_align_local_fitted_x_sp", &m_align_local_fitted_x_sp);
+ m_tree->Branch("fitParam_align_global_residual_y_sp", &m_align_global_residual_y_sp);
+ m_tree->Branch("fitParam_align_global_measured_x_sp", &m_align_global_measured_x_sp);
+ m_tree->Branch("fitParam_align_global_measured_y_sp", &m_align_global_measured_y_sp);
+ m_tree->Branch("fitParam_align_global_measured_ye_sp", &m_align_global_measured_ye_sp);
+ m_tree->Branch("fitParam_align_global_measured_z_sp", &m_align_global_measured_z_sp);
+ m_tree->Branch("fitParam_align_global_fitted_ye_sp", &m_align_global_fitted_ye_sp);
+ m_tree->Branch("fitParam_align_global_fitted_x_sp", &m_align_global_fitted_x_sp);
+ m_tree->Branch("fitParam_align_global_fitted_y_sp", &m_align_global_fitted_y_sp);
+}
+
+void CKF2Alignment::clearVariables(){
+ m_fitParam_x.clear();
+ m_fitParam_charge.clear();
+ m_fitParam_y.clear();
+ m_fitParam_z.clear();
+ m_fitParam_px.clear();
+ m_fitParam_py.clear();
+ m_fitParam_pz.clear();
+ m_fitParam_chi2.clear();
+ m_fitParam_ndf.clear();
+ m_fitParam_nHoles.clear();
+ m_fitParam_nOutliers.clear();
+ m_fitParam_nStates.clear();
+ m_fitParam_nMeasurements.clear();
+ m_align_derivation_x_x.clear();
+ m_align_derivation_x_y.clear();
+ m_align_derivation_x_z.clear();
+ m_align_derivation_x_rx.clear();
+ m_align_derivation_x_ry.clear();
+ m_align_derivation_x_rz.clear();
+ m_align_derivation_global_y_x.clear();
+ m_align_derivation_global_y_y.clear();
+ m_align_derivation_global_y_z.clear();
+ m_align_derivation_global_y_rx.clear();
+ m_align_derivation_global_y_ry.clear();
+ m_align_derivation_global_y_rz.clear();
+ m_align_local_residual_x_sp.clear();
+ m_align_unbiased_sp.clear();
+ m_align_local_pull_x_sp.clear();
+ m_align_local_measured_x_sp.clear();
+ m_align_local_measured_xe_sp.clear();
+ m_align_local_fitted_xe_sp.clear();
+ m_align_local_fitted_x_sp.clear();
+ m_align_global_residual_y_sp.clear();
+ m_align_global_measured_x_sp.clear();
+ m_align_global_measured_y_sp.clear();
+ m_align_global_measured_z_sp.clear();
+ m_align_global_measured_ye_sp.clear();
+ m_align_global_fitted_ye_sp.clear();
+ m_align_global_fitted_x_sp.clear();
+ m_align_global_fitted_y_sp.clear();
+ m_align_stationId_sp.clear();
+ m_align_centery_sp.clear();
+ m_align_layerId_sp.clear();
+ m_align_moduleId_sp.clear();
+ m_align_clustersize_sp.clear();
+ m_align_stereoId.clear();
+}
+
+
+double CKF2Alignment::getLocalDerivation(Amg::Transform3D trans1, Amg::Vector2D loc_pos, int ia, int side , Amg::Transform3D trans2)const {
+ double delta[6]={0.001,0.001,0.001,0.0001,0.0001,0.0001};
+ auto translation_m = Amg::Translation3D(0,0,0);
+ Amg::Transform3D transform_m= translation_m* Amg::RotationMatrix3D::Identity();
+ auto translation_p = Amg::Translation3D(0,0,0);
+ Amg::Transform3D transform_p= translation_p* Amg::RotationMatrix3D::Identity();
+ switch (ia)
+ {
+ case 0:
+ transform_m *= Amg::Translation3D(0-delta[ia],0,0);
+ transform_p *= Amg::Translation3D(0+delta[ia],0,0);
+ break;
+ case 1:
+ transform_m *= Amg::Translation3D(0,0-delta[ia],0);
+ transform_p *= Amg::Translation3D(0,0+delta[ia],0);
+ break;
+ case 2:
+ transform_m *= Amg::Translation3D(0,0,0-delta[ia]);
+ transform_p *= Amg::Translation3D(0,0,0+delta[ia]);
+ break;
+ case 3:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(1,0,0));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(1,0,0));
+ break;
+ case 4:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(0,1,0));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(0,1,0));
+ break;
+ case 5:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(0,0,1));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(0,0,1));
+ break;
+ default:
+ ATH_MSG_FATAL("Unknown alignment parameter" );
+ break;
+ }
+
+ auto local3 = Amg::Vector3D(loc_pos.x(),loc_pos.y(),0.);
+ Amg::Vector3D glo_pos = trans1 * local3;
+ Amg::Vector3D local_m = (trans1 * transform_m).inverse()*glo_pos;
+ Amg::Vector3D local_p = (trans1 * transform_p).inverse()*glo_pos;
+ if(side!=1)
+ return (local_p.x()-local_m.x())/2./delta[ia];
+ else{
+ Amg::Vector3D local_m2=trans2.inverse()*trans1*local_m;
+ Amg::Vector3D local_p2=trans2.inverse()*trans1*local_p;
+ return (local_p2.x()-local_m2.x())/2./delta[ia];
+ }
+}
+
+
+double CKF2Alignment::getDerivation(const EventContext& ctx, Acts::BoundTrackParameters& fitParameters, Amg::Transform3D trans1, Amg::Vector2D loc_pos, int ia, int side , Amg::Transform3D trans2, int global)const {
+ double delta[6]={0.001,0.001,0.001,0.0001,0.0001,0.0001};
+ auto translation_m = Amg::Translation3D(0,0,0);
+ Amg::Transform3D transform_m= translation_m* Amg::RotationMatrix3D::Identity();
+ auto translation_p = Amg::Translation3D(0,0,0);
+ Amg::Transform3D transform_p= translation_p* Amg::RotationMatrix3D::Identity();
+ switch (ia)
+ {
+ case 0:
+ transform_m *= Amg::Translation3D(0-delta[ia],0,0);
+ transform_p *= Amg::Translation3D(0+delta[ia],0,0);
+ break;
+ case 1:
+ transform_m *= Amg::Translation3D(0,0-delta[ia],0);
+ transform_p *= Amg::Translation3D(0,0+delta[ia],0);
+ break;
+ case 2:
+ transform_m *= Amg::Translation3D(0,0,0-delta[ia]);
+ transform_p *= Amg::Translation3D(0,0,0+delta[ia]);
+ break;
+ case 3:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(1,0,0));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(1,0,0));
+ break;
+ case 4:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(0,1,0));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(0,1,0));
+ break;
+ case 5:
+ transform_m *= Amg::AngleAxis3D(0-delta[ia], Amg::Vector3D(0,0,1));
+ transform_p *= Amg::AngleAxis3D(0+delta[ia], Amg::Vector3D(0,0,1));
+ break;
+ default:
+ ATH_MSG_FATAL("Unknown alignment parameter" );
+ break;
+ }
+
+ auto local3 = Amg::Vector3D(loc_pos.x(),loc_pos.y(),0.);
+ Amg::Vector3D glo_pos = trans2 * local3;
+ auto trans_m = trans1 * transform_m;
+ auto trans_p = trans1 * transform_p;
+ if(side==1 && global!=1){
+ trans_m = trans1 * transform_m * trans1.inverse() * trans2;
+ trans_p = trans1 * transform_p * trans1.inverse() * trans2;
+ }
+ if (global==1){
+ Amg::Transform3D trans_g = Amg::Translation3D(0,0, glo_pos.z()) * Amg::RotationMatrix3D::Identity();
+ auto trans_l2g = trans_g.inverse()*trans1;
+ trans_m = trans1 * trans_l2g.inverse() * transform_m * trans_l2g;
+ trans_p = trans1 * trans_l2g.inverse() * transform_p * trans_l2g;
+ if(side==1){
+ auto tmpm = trans_m;
+ auto tmpp = trans_p;
+ trans_m = tmpm * trans1.inverse() * trans2;
+ trans_p = tmpp * trans1.inverse() * trans2;
+ }
+ }
+
+ Amg::Transform3D ini_trans = Amg::Translation3D(0,0, glo_pos.z()-10) * Amg::RotationMatrix3D::Identity();
+ auto ini_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(ini_trans, std::make_shared<const Acts::RectangleBounds>(1000.,1000.));
+ std::unique_ptr<const Acts::BoundTrackParameters> ini_Param = m_extrapolationTool->propagate( ctx, fitParameters, *ini_surface, Acts::backward);
+ if(ini_Param==nullptr)return -9999.;
+
+
+ auto shift_m_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(trans_m, std::make_shared<const Acts::RectangleBounds>(1000.,1000.));
+ auto shift_p_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(trans_p, std::make_shared<const Acts::RectangleBounds>(1000.,1000.));
+ auto shift_m_Param = m_extrapolationTool->propagate( ctx, *ini_Param, *shift_m_surface, Acts::forward);
+ auto shift_p_Param = m_extrapolationTool->propagate( ctx, *ini_Param, *shift_p_surface, Acts::forward);
+ if(shift_m_Param!=nullptr||shift_p_Param!=nullptr){
+ auto shift_m_parameter = shift_m_Param->parameters();
+ auto shift_p_parameter = shift_p_Param->parameters();
+ Amg::Vector2D local_m(shift_m_parameter[Acts::eBoundLoc0],shift_m_parameter[Acts::eBoundLoc1]);
+ Amg::Vector2D local_p(shift_p_parameter[Acts::eBoundLoc0],shift_p_parameter[Acts::eBoundLoc1]);
+ return (local_p.x()-local_m.x())/2./delta[ia];
+ }
+ return -99999.;
+}
+
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.h
new file mode 100755
index 0000000000000000000000000000000000000000..1cc5f368942797b2a193ed69fdebf4bf2db2df4c
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2Alignment.h
@@ -0,0 +1,199 @@
+#ifndef FASERACTSKALMANFILTER_CKF2ALIGNMENT_H
+#define FASERACTSKALMANFILTER_CKF2ALIGNMENT_H
+
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "Acts/TrackFitting/KalmanFitter.hpp"
+#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
+#include "Acts/TrackFinding/MeasurementSelector.hpp"
+#include "FaserActsKalmanFilter/Measurement.h"
+#include "MagFieldConditions/FaserFieldCacheCondObj.h"
+#include "TrkTrack/TrackCollection.h"
+#include "FaserActsKalmanFilter/ITrackSeedTool.h"
+#include "KalmanFitterTool.h"
+#include <boost/dynamic_bitset.hpp>
+#include "GaudiKernel/ITHistSvc.h"
+#include "CreateTrkTrackTool.h"
+
+using ConstTrackStateProxy = Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>;
+using ClusterSet = boost::dynamic_bitset<>;
+
+class TTree;
+class FaserSCT_ID;
+
+
+namespace TrackerDD {
+ class SCT_DetectorManager;
+}
+
+class CKF2Alignment : public AthAlgorithm {
+ public:
+ CKF2Alignment(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~CKF2Alignment() = default;
+
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+ void initializeTree();
+ void clearVariables();
+
+ using TrackFinderOptions =
+ Acts::CombinatorialKalmanFilterOptions<IndexSourceLinkAccessor,
+ MeasurementCalibrator,
+ Acts::MeasurementSelector>;
+ using CKFResult = Acts::CombinatorialKalmanFilterResult<IndexSourceLink>;
+ using TrackFitterResult = Acts::Result<CKFResult>;
+ using TrackFinderResult = std::vector<TrackFitterResult>;
+
+ using KFResult =
+ Acts::Result<Acts::KalmanFitterResult<IndexSourceLink>>;
+
+ using TrackFitterOptions =
+ Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder,
+ Acts::VoidReverseFilteringLogic>;
+
+ using TrackParameters = Acts::CurvilinearTrackParameters;
+ using TrackParametersContainer = std::vector<TrackParameters>;
+
+ // Track Finding
+ class TrackFinderFunction {
+ public:
+ virtual ~TrackFinderFunction() = default;
+ virtual TrackFinderResult operator()(const IndexSourceLinkContainer&,
+ const TrackParametersContainer&,
+ const TrackFinderOptions&) const = 0;
+ };
+
+ static std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ bool resolvePassive, bool resolveMaterial, bool resolveSensitive);
+
+ // Track Fitting
+ class TrackFitterFunction {
+ public:
+ virtual ~TrackFitterFunction() = default;
+ virtual KFResult operator()(const std::vector<IndexSourceLink>&,
+ const Acts::BoundTrackParameters&,
+ const TrackFitterOptions&) const = 0;
+ };
+
+ struct TrajectoryInfo {
+ TrajectoryInfo(const FaserActsRecMultiTrajectory &traj) :
+ trajectory{traj}, clusterSet{nClusters} {
+ auto state = Acts::MultiTrajectoryHelpers::trajectoryState(traj.multiTrajectory(), traj.tips().front());
+ traj.multiTrajectory().visitBackwards(traj.tips().front(), [&](const ConstTrackStateProxy& state) {
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
+ return true;
+ }
+ clusterSet.set(state.uncalibrated().index());
+ return true;
+ });
+ nMeasurements = state.nMeasurements;
+ chi2 = state.chi2Sum;
+ }
+
+ static size_t nClusters;
+ FaserActsRecMultiTrajectory trajectory;
+ ClusterSet clusterSet;
+ size_t nMeasurements;
+ double chi2;
+ };
+
+ static std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry);
+
+ virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const;
+
+ private:
+ int m_numberOfEvents {0};
+ int m_numberOfTrackSeeds {0};
+ int m_numberOfFittedTracks {0};
+ int m_numberOfSelectedTracks {0};
+
+ void computeSharedHits(std::vector<IndexSourceLink>* sourceLinks, TrackFinderResult& results) const;
+ std::shared_ptr<TrackFinderFunction> m_fit;
+ std::shared_ptr<TrackFitterFunction> m_kf;
+ std::unique_ptr<const Acts::Logger> m_logger;
+ const FaserSCT_ID* m_idHelper {nullptr};
+ const TrackerDD::SCT_DetectorManager* m_detManager {nullptr};
+
+ Gaudi::Property<bool> m_biased{this, "BiasedResidual", true};
+ Gaudi::Property<std::string> m_actsLogging {this, "ActsLogging", "VERBOSE"};
+ Gaudi::Property<int> m_minNumberMeasurements {this, "MinNumberMeasurements", 12};
+ Gaudi::Property<bool> m_resolvePassive {this, "resolvePassive", false};
+ Gaudi::Property<bool> m_resolveMaterial {this, "resolveMaterial", true};
+ Gaudi::Property<bool> m_resolveSensitive {this, "resolveSensitive", true};
+ Gaudi::Property<double> m_maxSteps {this, "maxSteps", 10000};
+ Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15};
+ Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10};
+ SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
+ ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"};
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ ToolHandle<KalmanFitterTool> m_kalmanFitterTool1 {this, "KalmanFitterTool1", "KalmanFitterTool"};
+ ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
+ Gaudi::Property<bool> m_isMC {this, "isMC", false};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool{this, "ExtrapolationTool", "FaserActsExtrapolationTool"};
+
+ double getDerivation(const EventContext& ctx, Acts::BoundTrackParameters& fitParameters, Amg::Transform3D trans1, Amg::Vector2D loc_pos, int ia, int side , Amg::Transform3D trans2, int global)const ;
+ double getLocalDerivation(Amg::Transform3D trans1, Amg::Vector2D loc_pos, int ia, int side , Amg::Transform3D trans2)const;
+ //output ntuple
+ ServiceHandle<ITHistSvc> m_thistSvc;
+ TTree* m_tree;
+ mutable int m_eventId=0;
+ mutable int m_trackId=0;
+ mutable std::vector<double> m_fitParam_x;
+ mutable std::vector<double> m_fitParam_charge;
+ mutable std::vector<double> m_fitParam_y;
+ mutable std::vector<double> m_fitParam_z;
+ mutable std::vector<double> m_fitParam_px;
+ mutable std::vector<double> m_fitParam_py;
+ mutable std::vector<double> m_fitParam_pz;
+ mutable std::vector<double> m_fitParam_chi2;
+ mutable std::vector<int> m_fitParam_ndf;
+ mutable std::vector<int> m_fitParam_nHoles;
+ mutable std::vector<int> m_fitParam_nOutliers;
+ mutable std::vector<int> m_fitParam_nStates;
+ mutable std::vector<int> m_fitParam_nMeasurements;
+ mutable std::vector<std::vector<double>> m_align_stationId_sp;
+ mutable std::vector<std::vector<double>> m_align_centery_sp;
+ mutable std::vector<std::vector<double>> m_align_layerId_sp;
+ mutable std::vector<std::vector<double>> m_align_moduleId_sp;
+ mutable std::vector<std::vector<double>> m_align_clustersize_sp;
+ mutable std::vector<std::vector<double>> m_align_global_residual_y_sp;
+ mutable std::vector<std::vector<double>> m_align_global_measured_x_sp;
+ mutable std::vector<std::vector<double>> m_align_global_measured_y_sp;
+ mutable std::vector<std::vector<double>> m_align_global_measured_z_sp;
+ mutable std::vector<std::vector<double>> m_align_global_measured_ye_sp;
+ mutable std::vector<std::vector<double>> m_align_global_fitted_ye_sp;
+ mutable std::vector<std::vector<double>> m_align_global_fitted_y_sp;
+ mutable std::vector<std::vector<double>> m_align_global_fitted_x_sp;
+ mutable std::vector<std::vector<double>> m_align_local_residual_x_sp;
+ mutable std::vector<std::vector<double>> m_align_unbiased_sp;
+ mutable std::vector<std::vector<double>> m_align_local_pull_x_sp;
+ mutable std::vector<std::vector<double>> m_align_local_measured_x_sp;
+ mutable std::vector<std::vector<double>> m_align_local_measured_xe_sp;
+ mutable std::vector<std::vector<double>> m_align_local_fitted_xe_sp;
+ mutable std::vector<std::vector<double>> m_align_local_fitted_x_sp;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_x;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_y;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_z;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_rx;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_ry;
+ mutable std::vector<std::vector<double>> m_align_derivation_x_rz;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_x;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_y;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_z;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_rx;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_ry;
+ mutable std::vector<std::vector<double>> m_align_derivation_global_y_rz;
+ mutable std::vector<std::vector<double>> m_align_stereoId;
+};
+
+#endif // FASERACTSKALMANFILTER_CKF2Alignment_H
+
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index 41611bad5882de7e15ce3de15f6c4833fbf10581..952f75287f594e5df31feb2b6a493bfe235a4086 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -7,7 +7,6 @@
#include "Identifier/Identifier.h"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "CircleFit.h"
-#include "Identifier/Identifier.h"
#include "LinearFit.h"
#include "TrackClassification.h"
#include <array>
@@ -80,12 +79,14 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
if (std::find(maskedLayers.begin(), maskedLayers.end(), 3 * m_idHelper->station(id) + m_idHelper->layer(id)) != maskedLayers.end()) {
continue;
}
+// if(m_idHelper->station(id)==0) continue;
CircleFitTrackSeedTool::s_indexMap[cluster->identify()] = measurements.size();
if (identifierMap->count(id) != 0) {
Acts::GeometryIdentifier geoId = identifierMap->at(id);
IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
// create measurement
const auto& par = cluster->localPosition();
+ auto cova = cluster->localCovariance();
Eigen::Matrix<double, 1, 1> pos {par.x(),};
Eigen::Matrix<double, 1, 1> cov {0.0231 * 0.0231,};
ThisMeasurement meas(sourceLink, Indices, pos, cov);
@@ -101,7 +102,13 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
std::array<std::vector<Segment>, 4> segments {};
for (const Trk::Track* track : *trackCollection) {
auto s = Segment(track, m_idHelper, maskedLayers);
+ //remove IFT in seeding
+ if(m_removeIFT){
+ if (s.station != -1 &&s.station != 0) segments[s.station].push_back(s);
+ }
+ else{
if (s.station != -1) segments[s.station].push_back(s);
+ }
}
std::vector<Segment> combination {};
@@ -151,7 +158,7 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
selectedSeeds.begin(), selectedSeeds.end(), [](const Seed &lhs, const Seed &rhs) {
return lhs.minZ < rhs.minZ;
});
- double origin = !selectedSeeds.empty() ? minSeed->minZ - 1 : 0;
+ double origin = !selectedSeeds.empty() ? minSeed->minZ - 10 : 0;
m_targetZPosition = origin;
std::vector<Acts::CurvilinearTrackParameters> initParams {};
for (const Seed &seed : selectedSeeds) {
@@ -212,6 +219,7 @@ CircleFitTrackSeedTool::Segment::Segment(const Trk::Track* track, const FaserSCT
}
}
station = idHelper->station(id);
+// if(station==0)continue;
clusterSet.set(CircleFitTrackSeedTool::s_indexMap.at(id));
auto fitParameters = track->trackParameters()->front();
position = fitParameters->position();
@@ -265,12 +273,11 @@ CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
size = clusters.size();
}
+
void CircleFitTrackSeedTool::Seed::fit() {
CircleFit::CircleData circleData(positions);
CircleFit::Circle circle = CircleFit::circleFit(circleData);
momentum = circle.r > 0 ? circle.r * 0.001 * 0.3 * 0.55 : 9999999.;
- // the magnetic field bends a positively charged particle downwards, thus the
- // y-component of the center of a fitted circle is negative.
charge = circle.cy < 0 ? 1 : -1;
}
@@ -281,8 +288,6 @@ void CircleFitTrackSeedTool::Seed::fakeFit(double B) {
cy = circle.cy;
r = circle.r;
momentum = r * 0.3 * B * m_MeV2GeV;
- // the magnetic field bends a positively charged particle downwards, thus the
- // y-component of the center of a fitted circle is negative.
charge = circle.cy < 0 ? 1 : -1;
}
@@ -292,6 +297,7 @@ void CircleFitTrackSeedTool::Seed::linearFit(const std::vector<Acts::Vector2> &p
c0 = origin[1] - origin[0] * c1;
}
+
double CircleFitTrackSeedTool::Seed::getY(double z) {
double sqt = std::sqrt(-cx*cx + 2*cx*z + r*r - z*z);
return abs(cy - sqt) < abs(cy + sqt) ? cy - sqt : cy + sqt;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
index 7866b89fc32fe2ccd8a0391f421b1d203ff1e8e0..0467c92620eebd955f0a78f11cfdea3fd0cb0b1d 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
@@ -116,6 +116,7 @@ private:
Gaudi::Property<double> m_covTheta {this, "covTheta", 1};
Gaudi::Property<double> m_covQOverP {this, "covQOverP", 1};
Gaudi::Property<double> m_covTime {this, "covTime", 1};
+ Gaudi::Property<bool> m_removeIFT{this, "removeIFT", false};
};
inline const std::shared_ptr<std::vector<Acts::CurvilinearTrackParameters>>
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
index ee2a864b51d0cd632322128256666679c4beb490..dcb167496790d2b5af143b55b86cc4698f902976 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx
@@ -73,7 +73,7 @@ StatusCode CombinatorialKalmanFilterAlg::execute() {
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
- const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext();
+ const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext();
auto geoctx = gctx.context();
Acts::MagneticFieldContext magFieldContext = getMagneticFieldContext(ctx);
Acts::CalibrationContext calibContext;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
index f88b4e15f59fc7e7d3a4b434c8328942982a67ba..5dd831f9592b5323280d5cc09e10798fa342966f 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
@@ -109,7 +109,7 @@ StatusCode FaserActsKalmanFilterAlg::execute() {
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
- const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext();
+ const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext();
auto geoctx = gctx.context();
Acts::MagneticFieldContext magctx = getMagneticFieldContext(ctx);
Acts::CalibrationContext calctx;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index 051687e21c6d0b1c1874a26ac4ddeb5d55c55124..c0d02183f5616df5b624152fee61d4c1bf6bba37 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -39,6 +39,329 @@ StatusCode KalmanFitterTool::finalize() {
return StatusCode::SUCCESS;
}
+//get the unbiased residual for the cluster at Z = clusz
+//output is a vector of the outliers
+//
+std::vector<TSOS4Residual>
+KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track* inputTrack, double clusz, const Acts::BoundVector& inputVector,
+ bool isMC, double origin) const {
+ std::vector<TSOS4Residual> resi;
+ resi.clear();
+ std::vector<FaserActsRecMultiTrajectory> myTrajectories;
+ ATH_MSG_DEBUG("start kalmanfilter "<<inputTrack->measurementsOnTrack()->size()<<" "<<origin<<" "<<clusz);
+
+ if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
+ ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
+ return resi;
+ }
+
+ if (!inputTrack->trackParameters() || inputTrack->trackParameters()->empty()) {
+ ATH_MSG_DEBUG("Input track has no track parameters and cannot be fitted");
+ return resi;
+ }
+
+
+ auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
+
+ Acts::MagneticFieldContext mfContext = getMagneticFieldContext(ctx);
+ Acts::CalibrationContext calibContext = Acts::CalibrationContext();
+
+ auto [sourceLinks, measurements] = getMeasurementsFromTrack(inputTrack);
+ auto trackParameters = getParametersFromTrack(inputTrack->trackParameters()->front(), inputVector, origin);
+ ATH_MSG_DEBUG("trackParameters: " << trackParameters.parameters().transpose());
+ ATH_MSG_DEBUG("position: " << trackParameters.position(gctx).transpose());
+ ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
+ ATH_MSG_DEBUG("charge: " << trackParameters.charge());
+ ATH_MSG_DEBUG("size of measurements : " << sourceLinks.size());
+
+ FaserActsOutlierFinder faserActsOutlierFinder{0};
+ faserActsOutlierFinder.cluster_z=clusz;
+ faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
+ Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
+ kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
+ faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ //FaserActsOutlierFinder(), Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ //Acts::VoidOutlierFinder(), Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ Acts::PropagatorPlainOptions(), &(*pSurface));
+ kfOptions.multipleScattering = false;
+ kfOptions.energyLoss = false;
+
+ ATH_MSG_DEBUG("Invoke fitter");
+ auto result = (*m_fit)(sourceLinks, trackParameters, kfOptions);
+
+ if (result.ok()) {
+ const auto& fitOutput = result.value();
+ if (fitOutput.fittedParameters) {
+ const auto& params = fitOutput.fittedParameters.value();
+ ATH_MSG_DEBUG("ReFitted parameters");
+ ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ ATH_MSG_DEBUG(" charge: " << params.charge());
+ using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
+ IndexedParams indexedParams;
+ indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
+ std::vector<size_t> trackTips;
+ trackTips.reserve(1);
+ trackTips.emplace_back(fitOutput.lastMeasurementIndex);
+ myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
+ std::move(trackTips),
+ std::move(indexedParams));
+ for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
+ const auto& mj = traj.multiTrajectory();
+ const auto& trackTips = traj.tips();
+ if(trackTips.empty())continue;
+ auto& trackTip = trackTips.front();
+ mj.visitBackwards(trackTip, [&](const auto &state) {
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ return true;
+ }
+
+ if (state.hasUncalibrated()&&state.hasSmoothed()) {
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance());
+ auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
+ const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ ATH_MSG_DEBUG(" residual/global position: " << residual(Acts::eBoundLoc0)<<" "<<parameter.position(gctx).x()<<" "<<parameter.position(gctx).y()<<" "<<parameter.position(gctx).z());
+ }
+ return true;
+ });
+
+ }
+
+ } else {
+ ATH_MSG_DEBUG("No fitted parameters for track");
+ }
+ }
+ return resi;
+
+}
+
+//get the unbiased residual for IFT, i.e. remove whole IFT in the track fitting
+//output is a vector of the outliers
+//
+std::vector<TSOS4Residual>
+KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track* inputTrack, const Acts::BoundVector& inputVector,
+ bool isMC, double origin, std::vector<const Tracker::FaserSCT_Cluster*>& clusters, const Acts::BoundTrackParameters ini_Param) const {
+ std::vector<TSOS4Residual> resi;
+ resi.clear();
+ std::vector<FaserActsRecMultiTrajectory> myTrajectories;
+
+ if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
+ ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
+ return resi;
+ }
+
+ if (!inputTrack->trackParameters() || inputTrack->trackParameters()->empty()) {
+ ATH_MSG_DEBUG("Input track has no track parameters and cannot be fitted");
+ return resi;
+ }
+
+
+ auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
+
+ Acts::MagneticFieldContext mfContext = getMagneticFieldContext(ctx);
+ Acts::CalibrationContext calibContext = Acts::CalibrationContext();
+
+ auto [sourceLinks, measurements] = getMeasurementsFromTrack(inputTrack/*, wafer_id*/, clusters);
+ ATH_MSG_DEBUG("trackParameters: " << ini_Param.parameters().transpose());
+ ATH_MSG_DEBUG("position: " << ini_Param.position(gctx).transpose());
+ ATH_MSG_DEBUG("momentum: " << ini_Param.momentum().transpose());
+ ATH_MSG_DEBUG("charge: " << ini_Param.charge());
+
+ FaserActsOutlierFinder faserActsOutlierFinder{0};
+ faserActsOutlierFinder.cluster_z=-1000000.;
+ faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
+ Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
+ kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
+ faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ Acts::PropagatorPlainOptions(), &(*pSurface));
+ kfOptions.multipleScattering = false;
+ kfOptions.energyLoss = false;
+
+ ATH_MSG_DEBUG("Invoke fitter");
+ auto result = (*m_fit)(sourceLinks, ini_Param, kfOptions);
+
+ if (result.ok()) {
+ const auto& fitOutput = result.value();
+ if (fitOutput.fittedParameters) {
+ const auto& params = fitOutput.fittedParameters.value();
+ ATH_MSG_DEBUG("ReFitted parameters");
+ ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ ATH_MSG_DEBUG(" charge: " << params.charge());
+ using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
+ IndexedParams indexedParams;
+ indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
+ std::vector<size_t> trackTips;
+ trackTips.reserve(1);
+ trackTips.emplace_back(fitOutput.lastMeasurementIndex);
+ myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
+ std::move(trackTips),
+ std::move(indexedParams));
+ for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
+ const auto& mj = traj.multiTrajectory();
+ const auto& trackTips = traj.tips();
+ if(trackTips.empty())continue;
+ auto& trackTip = trackTips.front();
+ mj.visitBackwards(trackTip, [&](const auto &state) {
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ return true;
+ }
+
+ if (state.hasUncalibrated()&&state.hasSmoothed()) {
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance());
+ auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
+ const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ }
+ return true;
+ });
+
+ }
+
+ } else {
+ ATH_MSG_DEBUG("No fitted parameters for track");
+ }
+ }
+ return resi;
+
+}
+
+//get the unbiased residual for IFT, i.e. remove whole IFT in the track fitting
+//output is a vector of the outliers
+//
+std::vector<TSOS4Residual>
+KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track* inputTrack, const Acts::BoundVector& inputVector,
+ bool isMC, double origin) const {
+ std::vector<TSOS4Residual> resi;
+ resi.clear();
+ std::vector<FaserActsRecMultiTrajectory> myTrajectories;
+
+ if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
+ ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
+ return resi;
+ }
+
+ if (!inputTrack->trackParameters() || inputTrack->trackParameters()->empty()) {
+ ATH_MSG_DEBUG("Input track has no track parameters and cannot be fitted");
+ return resi;
+ }
+
+
+ auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
+
+ Acts::MagneticFieldContext mfContext = getMagneticFieldContext(ctx);
+ Acts::CalibrationContext calibContext = Acts::CalibrationContext();
+
+ auto [sourceLinks, measurements] = getMeasurementsFromTrack(inputTrack/*, wafer_id*/);
+ auto trackParameters = getParametersFromTrack(inputTrack->trackParameters()->front(), inputVector, origin);
+ ATH_MSG_DEBUG("trackParameters: " << trackParameters.parameters().transpose());
+ ATH_MSG_DEBUG("position: " << trackParameters.position(gctx).transpose());
+ ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
+ ATH_MSG_DEBUG("charge: " << trackParameters.charge());
+
+
+ FaserActsOutlierFinder faserActsOutlierFinder{0};
+ faserActsOutlierFinder.cluster_z=-1000000.;
+ faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
+ Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
+ kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
+ faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ Acts::PropagatorPlainOptions(), &(*pSurface));
+ kfOptions.multipleScattering = false;
+ kfOptions.energyLoss = false;
+
+ ATH_MSG_DEBUG("Invoke fitter");
+ auto result = (*m_fit)(sourceLinks, trackParameters, kfOptions);
+
+ if (result.ok()) {
+ const auto& fitOutput = result.value();
+ if (fitOutput.fittedParameters) {
+ const auto& params = fitOutput.fittedParameters.value();
+ ATH_MSG_DEBUG("ReFitted parameters");
+ ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ ATH_MSG_DEBUG(" charge: " << params.charge());
+ using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
+ IndexedParams indexedParams;
+ indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
+ std::vector<size_t> trackTips;
+ trackTips.reserve(1);
+ trackTips.emplace_back(fitOutput.lastMeasurementIndex);
+ myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
+ std::move(trackTips),
+ std::move(indexedParams));
+ for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
+ const auto& mj = traj.multiTrajectory();
+ const auto& trackTips = traj.tips();
+ if(trackTips.empty())continue;
+ auto& trackTip = trackTips.front();
+ mj.visitBackwards(trackTip, [&](const auto &state) {
+ auto typeFlags = state.typeFlags();
+ if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ return true;
+ }
+
+ if (state.hasUncalibrated()&&state.hasSmoothed()) {
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance());
+ auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
+ const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ }
+ return true;
+ });
+
+ }
+
+ } else {
+ ATH_MSG_DEBUG("No fitted parameters for track");
+ }
+ }
+ return resi;
+
+}
+
std::unique_ptr<Trk::Track>
KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx,
Trk::Track* inputTrack, const Acts::BoundVector& inputVector,
@@ -70,9 +393,12 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
ATH_MSG_DEBUG("charge: " << trackParameters.charge());
- Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic>
+ FaserActsOutlierFinder faserActsOutlierFinder{0};
+ faserActsOutlierFinder.cluster_z=-1000000.;
+ faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
+ Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
- Acts::VoidOutlierFinder(), Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
+ faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
Acts::PropagatorPlainOptions(), &(*pSurface));
kfOptions.multipleScattering = false;
kfOptions.energyLoss = false;
@@ -95,9 +421,6 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
std::vector<size_t> trackTips;
trackTips.reserve(1);
trackTips.emplace_back(fitOutput.lastMeasurementIndex);
- // trajectories.emplace_back(std::move(fitOutput.fittedStates),
- // std::move(trackTips),
- // std::move(indexedParams));
myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
std::move(trackTips),
std::move(indexedParams));
@@ -172,6 +495,84 @@ Acts::MagneticFieldContext KalmanFitterTool::getMagneticFieldContext(const Event
}
+std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
+KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, std::vector<const Tracker::FaserSCT_Cluster*>& clusters) const {
+ const int kSize = 1;
+ ATH_MSG_DEBUG("clusters in refit "<<clusters.size());
+ std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0};
+ using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>;
+ using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
+
+ std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
+ std::vector<IndexSourceLink> sourceLinks;
+ std::vector<Measurement> measurements;
+ for(auto cluster : clusters){
+ Identifier id = cluster->identify();
+ int ista = m_idHelper->station(id);
+ int ilay = m_idHelper->layer(id);
+ ATH_MSG_DEBUG("station/layer "<<ista<<" "<<ilay);
+ Identifier waferId = m_idHelper->wafer_id(id);
+ if (identifierMap->count(waferId) != 0) {
+ ATH_MSG_DEBUG("found the identifier "<<ista<<" "<<ilay);
+ Acts::GeometryIdentifier geoId = identifierMap->at(waferId);
+ IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
+ Eigen::Matrix<double, 1, 1> pos {cluster->localPosition().x()};
+ Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
+ ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ sourceLinks.push_back(sourceLink);
+ measurements.emplace_back(std::move(actsMeas));
+ }
+ }
+ for (const Trk::MeasurementBase *meas : *track->measurementsOnTrack()) {
+ const auto* clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(meas);
+ const Tracker::FaserSCT_Cluster* cluster = clusterOnTrack->prepRawData();
+ if (clusterOnTrack) {
+ Identifier id = clusterOnTrack->identify();
+ Identifier waferId = m_idHelper->wafer_id(id);
+ if (identifierMap->count(waferId) != 0) {
+ Acts::GeometryIdentifier geoId = identifierMap->at(waferId);
+ IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
+ Eigen::Matrix<double, 1, 1> pos {meas->localParameters()[Trk::locX],};
+ Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
+ ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ sourceLinks.push_back(sourceLink);
+ measurements.emplace_back(std::move(actsMeas));
+ }
+ }
+ }
+ return std::make_tuple(sourceLinks, measurements);
+}
+
+std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
+KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, Identifier& wafer_id) const {
+ const int kSize = 1;
+ std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0};
+ using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>;
+ using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
+
+ std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
+ std::vector<IndexSourceLink> sourceLinks;
+ std::vector<Measurement> measurements;
+ for (const Trk::MeasurementBase *meas : *track->measurementsOnTrack()) {
+ const auto* clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(meas);
+ const Tracker::FaserSCT_Cluster* cluster = clusterOnTrack->prepRawData();
+ if (clusterOnTrack) {
+ Identifier id = clusterOnTrack->identify();
+ Identifier waferId = m_idHelper->wafer_id(id);
+
+ if (identifierMap->count(waferId) != 0) {
+ Acts::GeometryIdentifier geoId = identifierMap->at(waferId);
+ IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
+ Eigen::Matrix<double, 1, 1> pos {meas->localParameters()[Trk::locX],};
+ Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
+ ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ sourceLinks.push_back(sourceLink);
+ measurements.emplace_back(std::move(actsMeas));
+ }
+ }
+ }
+ return std::make_tuple(sourceLinks, measurements);
+}
std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track) const {
const int kSize = 1;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
index 0894e3ab2139decc1f0950a89f6200dc9ea2930f..584cda2c98bf32d7050150a690a0926439e11f91 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
@@ -15,9 +15,54 @@
#include "TrkTrack/Track.h"
#include "TrkTrack/TrackCollection.h"
#include "CreateTrkTrackTool.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "Identifier/Identifier.h"
+struct TSOS4Residual{
+ double fit_local_x;
+ double fit_local_y;
+ double fit_global_x;
+ double fit_global_y;
+ double fit_global_z;
+ const Tracker::FaserSCT_Cluster* cluster;
+ double residual;
+ Acts::BoundTrackParameters parameter;
+
+};
class FaserSCT_ID;
+//set the cluster to be removed as outlier in order to get the unbiased residual
+/// Outlier finder using a Chi2 cut.
+struct FaserActsOutlierFinder {
+ double StateChiSquaredPerNumberDoFCut = 10000.;
+ double cluster_z = -10000.;
+ template <typename track_state_t>
+ bool operator()(const track_state_t& state) const {
+ //remove the whole IFT
+ if(cluster_z<-10000){
+ if(state.uncalibrated().hit()->globalPosition().z()<-100)return true;
+ }
+
+ if (not state.hasCalibrated() or not state.hasPredicted()) {
+ return false;
+ }
+ return Acts::visit_measurement(
+ state.calibrated(), state.calibratedCovariance(),
+ state.calibratedSize(),
+ [&](const auto calibrated, const auto calibratedCovariance) {
+ //remove the cluster
+ if(fabs(state.uncalibrated().hit()->globalPosition().z()-cluster_z)<3)return true;
+ constexpr size_t kMeasurementSize = decltype(calibrated)::RowsAtCompileTime;
+ const auto H =
+ state.projector()
+ .template topLeftCorner<kMeasurementSize, Acts::BoundIndices::eBoundSize>()
+ .eval();
+ const auto residual = calibrated - H * state.predicted();
+ double chi2 = (residual.transpose() * ((calibratedCovariance + H * state.predictedCovariance() * H.transpose())).inverse() * residual).value();
+ return bool(chi2 > StateChiSquaredPerNumberDoFCut * kMeasurementSize);
+ });
+ }
+};
class KalmanFitterTool : virtual public AthAlgTool {
public:
@@ -29,7 +74,8 @@ public:
using TrackParameters = Acts::BoundTrackParameters;
using IndexedParams = std::unordered_map<size_t, TrackParameters>;
using TrackFitterOptions =
- Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic>;
+ Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>;
+ //Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic>;
using TrackFitterResult = Acts::Result<Acts::KalmanFitterResult<IndexSourceLink>>;
class TrackFitterFunction {
public:
@@ -45,12 +91,27 @@ public:
Trk::Track *inputTrack,
const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(),
bool isMC=false, double origin=0) const;
+ std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track *inputTrack,
+ const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(),
+ bool isMC=false, double origin=0) const;
+ std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track *inputTrack, double clusz,
+ const Acts::BoundVector& inputVector ,
+ bool isMC, double origin) const;
+ std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
+ Trk::Track *inputTrack,
+ const Acts::BoundVector& inputVector ,
+ bool isMC, double origin, std::vector<const Tracker::FaserSCT_Cluster*>& clus,const Acts::BoundTrackParameters ini_Param) const;
private:
const FaserSCT_ID* m_idHelper {nullptr};
std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
getMeasurementsFromTrack(Trk::Track *track) const;
- // Acts::BoundTrackParameters getParametersFromTrack(const Acts::BoundVector& params, const Trk::TrackParameters *inputParameters) const;
+ std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
+ getMeasurementsFromTrack(Trk::Track *track, Identifier& wafer_id) const;
+ std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
+ getMeasurementsFromTrack(Trk::Track *track, std::vector<const Tracker::FaserSCT_Cluster*>& clusters) const;
Acts::BoundTrackParameters getParametersFromTrack(const Trk::TrackParameters *inputParameters, const Acts::BoundVector& inputVector, double origin) const;
std::shared_ptr<TrackFitterFunction> m_fit;
std::unique_ptr<const Acts::Logger> m_logger;
@@ -67,6 +128,12 @@ private:
ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
+
+// Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> getExtensions();
+
+// FaserActsOutlierFinder m_outlierFinder{0};
+// ReverseFilteringLogic m_reverseFilteringLogic{0};
+// Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> m_kfExtensions;
};
#endif //FASERACTSKALMANFILTER_KALMANFITTERTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TruthSeededTrackFinderTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TruthSeededTrackFinderTool.cxx
index f0d9796f1a16ce52aee2bf52dbc11ec28c96a7a5..7f9f427960ba221d8fd3f489cbb0e38d627ffc9f 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TruthSeededTrackFinderTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TruthSeededTrackFinderTool.cxx
@@ -46,7 +46,7 @@ StatusCode TruthSeededTrackFinderTool::run() {
SG::ReadHandle<SpacePointForSeedCollection> spacePointCollection {m_spacePointCollectionKey};
ATH_CHECK(spacePointCollection.isValid());
- Acts::GeometryContext geoctx = m_trackingGeometryTool->getNominalGeometryContext().context();
+ Acts::GeometryContext geoctx = m_trackingGeometryTool->getGeometryContext().context();
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TruthTrackFinderTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TruthTrackFinderTool.cxx
index bfd5e9b740bbeb53529fc2c6e1048656187e1d5a..672be6fa934ead4725a43d81bc2f5e156612449e 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TruthTrackFinderTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TruthTrackFinderTool.cxx
@@ -39,7 +39,7 @@ StatusCode TruthTrackFinderTool::run() {
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
- Acts::GeometryContext gctx = m_trackingGeometryTool->getNominalGeometryContext().context();
+ Acts::GeometryContext gctx = m_trackingGeometryTool->getGeometryContext().context();
const int kSize = 2;
using ParametersVector = Acts::ActsVector<kSize>;
using CovarianceMatrix = Acts::ActsSymMatrix<kSize>;
@@ -188,4 +188,4 @@ StatusCode TruthTrackFinderTool::run() {
StatusCode TruthTrackFinderTool::finalize() {
return StatusCode::SUCCESS;
-}
\ No newline at end of file
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
index d9fefbddfb3b5e651f5ac697ac231bb9e41918d3..be32f4959409c57903c674bf7fd5d28209ad295d 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
@@ -23,6 +23,7 @@
#include "../TrackSeedWriterTool.h"
#include "../ActsTrackSeedTool.h"
#include "../CKF2.h"
+#include "../CKF2Alignment.h"
#include "../KalmanFitterTool.h"
#include "../MyTrackSeedTool.h"
#include "../SeedingAlg.h"
@@ -53,6 +54,7 @@ DECLARE_COMPONENT(PerformanceWriterTool)
DECLARE_COMPONENT(TrackSeedWriterTool)
DECLARE_COMPONENT(ActsTrackSeedTool)
DECLARE_COMPONENT(CKF2)
+DECLARE_COMPONENT(CKF2Alignment)
DECLARE_COMPONENT(KalmanFitterTool)
DECLARE_COMPONENT(MyTrackSeedTool)
DECLARE_COMPONENT(SeedingAlg)
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index c30351902a1b79e9a58ad33d5f58a5d04ed13faa..37cb7af183c2d89b9ada3579f76820bfe6487c5c 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -51,7 +51,7 @@ public:
virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0;
/// Generate random baseline
- virtual unsigned int generate_baseline(int mean, int rms) const = 0;
+ virtual float generate_baseline(float mean, float rms) const = 0;
/// Create structure to store pulse for each channel
template <class T>
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index c62d7f8a753490907bd1e8a513cad5341a2a46be..3a48027d57e3710a2ec592aed6a9ed749b9f6848 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -22,7 +22,7 @@ WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
m_nsamples = 600;
- m_random = new TRandom3();
+ m_random = new TRandom3(0);
return StatusCode::SUCCESS;
}
@@ -41,7 +41,7 @@ WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
return timekernel;
}
-unsigned int
-WaveformDigitisationTool::generate_baseline(int mean, int rms) const {
+float
+WaveformDigitisationTool::generate_baseline(float mean, float rms) const {
return m_random->Gaus(mean, rms);
}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index e2dd5169152845824927baeeae7ce8fc36ab46f8..b6554c337830fbd2183752cc1b05123ea827b4f1 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -32,7 +32,7 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
std::vector<float> evaluate_timekernel(TF1* kernel) const;
/// Generate random baseline
- unsigned int generate_baseline(int mean, int rms) const;
+ float generate_baseline(float mean, float rms) const;
private:
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
index 8882027bfc13a7b28a7043f2b9feb70ed8e25877..f7ef9339ae32919ef6600cf099da2d42c5eb4a69 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
@@ -759,6 +759,12 @@ WaveformReconstructionTool::fitGaussian(const WaveformFitResult& raw, const std:
// Find time here
gfit.time = gfunc.GetX( gfit.peak * m_timingPeakFraction, time.front(), gfit.mean );
+ if (isnan(gfit.time)) {
+ ATH_MSG_WARNING(" Gaussian fit returned time as NaN, assume fit failed");
+ gfit = raw;
+ gfit.fit_status = false;
+ }
+
ATH_MSG_DEBUG("G Fit Mean: " << gfit.mean << " RMS: " << gfit.sigma
<< " Peak: " << gfit.peak << " Int: " << gfit.integral << " Time: " << gfit.time);
}
@@ -839,6 +845,12 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit,
// Find time here
cbfit.time = cbfunc.GetX( cbfit.peak * m_timingPeakFraction, time.front(), cbfit.mean );
+ if (isnan(cbfit.time)) {
+ ATH_MSG_WARNING(" Crystal Ball fit returned time as NaN, assume fit failes");
+ cbfit.time = gfit.time;
+ cbfit.fit_status = false;
+ }
+
ATH_MSG_DEBUG("CB Fit Mean: " << cbfit.mean << " RMS: " << cbfit.sigma
<< " Peak: " << cbfit.peak << " Int: " << cbfit.integral
<< " Time: " << cbfit.time
diff --git a/Waveform/WaveformConditions/WaveCondUtils/scripts/makeDigiDB.py b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeDigiDB.py
index 202c5e266107d9e753c26b79d4de470d2f35653e..40e4247398e516c884152f5e30de34fd26977a8c 100755
--- a/Waveform/WaveformConditions/WaveCondUtils/scripts/makeDigiDB.py
+++ b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeDigiDB.py
@@ -108,7 +108,13 @@ preshower_record['cb_sigma'] = 4.0
preshower_record['cb_alpha'] = -0.32
preshower_record['cb_n'] = 1000.
+# Dec. 10, 2022
+# Values: 0.56, 5.6, 16.8, 168
+# Dec. 19, 2022
+# Lowered by 10% to 0.51, 5.1, 15.3, 153
+
# Low gain, filters in default
+calo_record['cb_norm'] = 0.51
tag = "WAVE-Digitization-01-LG"
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, calo_record, 0, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, vetonu_record, 1, tag)
@@ -126,7 +132,7 @@ wave_folder.createTagRelation("OFLCOND-FASER-03", "WAVE-01")
# Low gain, filters out
tag = "WAVE-Digitization-01-LG-nofilt"
-calo_record['cb_norm'] = 5.6
+calo_record['cb_norm'] = 5.1
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, calo_record, 0, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, vetonu_record, 1, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, veto_record, 2, tag)
@@ -135,7 +141,7 @@ folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, preshower_record,
# High gain
tag = "WAVE-Digitization-01-HG"
-calo_record['cb_norm'] = 16.8
+calo_record['cb_norm'] = 15.3
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, calo_record, 0, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, vetonu_record, 1, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, veto_record, 2, tag)
@@ -144,7 +150,7 @@ folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, preshower_record,
# High gain, filters out
tag = "WAVE-Digitization-01-HG-nofilt"
-calo_record['cb_norm'] = 168
+calo_record['cb_norm'] = 153
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, calo_record, 0, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, vetonu_record, 1, tag)
folder.storeObject( cool.ValidityKeyMin, cool.ValidityKeyMax, veto_record, 2, tag)