diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 4ac399e08ce06c47fc792e9f32d0c45806cabafc..7af8f3d219ce0095d580f92d98ddd9b1a52c3478 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -75,6 +75,8 @@ def faser_pgparser():
help="Write out full configuration")
parser.add_argument("--noexec", action='store_true',
help="Exit after parsing configuration (no execution)")
+ parser.add_argument('--filter',default="",type=str,
+ help='Specify "type" of filter, e.g. "muon_conversion" ,etc')
pg_args = parser.parse_args()
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index bf3cd11cb1e31e182b48372dc476c3766a2deee5..6a3c34308187d4407f3f01d63f88897a88006716 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -265,6 +265,16 @@ if __name__ == '__main__':
#
from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg
cfg.merge(G4FaserAlgCfg(ConfigFlags))
+
+
+#Event Filtering?
+ if args.filter == 'muon_conversion':
+ import McParticleEvent.Pythonizations
+ from GeneratorUtils.MuonFilters import ConversionFilter
+ filt = ConversionFilter("MuonFilter") # , muonEnergyLoss = 5000, minConversionEnergy = 2000) # , muonEnergyLoss = 5000, maxRadius = 100)
+ cfg.addEventAlgo(filt, sequenceName = "AthAlgSeq")
+ OutputStreamHITS = cfg.getEventAlgo("OutputStreamHITS")
+ OutputStreamHITS.RequireAlgs = ["MuonFilter"]
#
# Dump config
#
diff --git a/Generators/GeneratorUtils/python/MuonFilters.py b/Generators/GeneratorUtils/python/MuonFilters.py
new file mode 100644
index 0000000000000000000000000000000000000000..2540d3905ada00ccba8ffacf2f3f56e43d4b8248
--- /dev/null
+++ b/Generators/GeneratorUtils/python/MuonFilters.py
@@ -0,0 +1,132 @@
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+from AthenaPython import PyAthena
+from AthenaPython.PyAthena import StatusCode, McEventCollection, CLHEP
+from AthenaCommon.SystemOfUnits import cm, m, GeV
+import ROOT
+from math import sqrt
+
+try:
+ from AthenaPython.PyAthena import HepMC3 as HepMC
+except ImportError:
+ from AthenaPython.PyAthena import HepMC as HepMC
+
+class ConversionFilter(PyAthena.Alg):
+ def __init__(self, name="SimFilter", InputMCEventKey="TruthEvent", minConversionEnergy = None,
+ maxRadius = None, muonEnergyLoss = None):
+
+ super(ConversionFilter,self).__init__(name=name)
+
+ self.InputMCEventKey = InputMCEventKey
+ self.minConversionEnergy = minConversionEnergy
+ self.maxRadius = maxRadius
+ self.muonEnergyLoss = muonEnergyLoss
+
+ self.zEmulsionEnd = -1.6 * m
+
+ return
+
+ def findConversionFromMuon(self, evt):
+
+ for i, p in enumerate(evt.particles):
+
+ # Check incoming particle (muon) is within radius to
+ # avoid those conversions from muons coming in to magnet
+ if i == 1 and self.maxRadius is not None:
+ pvtx = p.production_vertex()
+ if pvtx:
+ ppos = pvtx.position()
+ if sqrt(ppos.x()**2 + ppos.y()**2) > self.maxRadius:
+ return False
+
+ # Find photon
+ if not abs(p.pdg_id()) == 22: continue
+
+ # Decaying to 2 children
+ ovtx = p.end_vertex()
+ if not ovtx: continue
+ if ovtx.particles_out_size() != 2: continue
+
+ # That are e+e- pair
+ children = list(ovtx.particles_out)
+ if not (abs(children[0].pdg_id()) == 11 and children[0].pdg_id() == -children[1].pdg_id()): continue
+ self.msg.debug("Found Photon conversion")
+
+ # Which comes from a single mother
+ ivtx = p.production_vertex()
+ if not ivtx: continue
+ if ivtx.particles_in_size() != 1: continue
+ mother = list(ivtx.particles_in)[0]
+
+ # That is a muon (allowing for another photon inbetween)
+ isFromMuon = False
+ if abs(mother.pdg_id()) == 13:
+ isFromMuon = True
+ elif abs(mother.pdg_id()) == 22:
+ ivtx = mother.production_vertex()
+ if ivtx and ivtx.particles_in_size() == 1:
+ gran = list(ivtx.particles_in)[0]
+ if abs(gran.pdg_id()) == 13:
+ isFromMuon = True
+
+ if isFromMuon:
+ # Optionally require a certain energy for at least one of the e+e- pair
+ if (self.minConversionEnergy is None
+ or children[0].momentum().e() > self.minConversionEnergy
+ or children[1].momentum().e() > self.minConversionEnergy):
+
+ self.msg.debug("... from muon")
+ self.setFilterPassed(True)
+ return True
+
+ return False
+
+ def findLargeELossMuon(self, evt):
+
+ lastMuonE = None
+ firstMuonE = None
+
+ for p in evt.particles:
+
+ # Find muons
+ if abs(p.pdg_id()) != 13: continue
+
+ # With a vertex inside the specified radius
+ pvtx = p.production_vertex()
+ if not pvtx: continue
+ ppos = pvtx.position()
+
+ if self.maxRadius is not None and sqrt(ppos.x()**2 + ppos.y()**2) > self.maxRadius:
+ continue
+
+ # Find the final muon produced in the emulsion
+ if ppos.z() < self.zEmulsionEnd:
+ firstMuonE = p.momentum().e()
+
+ # Find the final muon
+ lastMuonE = p.momentum().e()
+
+ # Look for E loss between end of emulsion and final
+ if lastMuonE is not None and firstMuonE is not None:
+ if firstMuonE - lastMuonE > self.muonEnergyLoss:
+ self.msg.debug(f"Found muon with E loss = {firstMuonE - lastMuonE}")
+ self.setFilterPassed(True)
+ return True
+
+ return False
+
+
+
+ def execute(self):
+ self.msg.debug(f"Executing {self.getName()}")
+
+ self.msg.debug(f"Reading {self.InputMCEventKey}")
+ evt = self.evtStore[self.InputMCEventKey][0]
+
+ self.setFilterPassed(False)
+
+ if self.muonEnergyLoss:
+ self.findLargeELossMuon(evt)
+ else:
+ self.findConversionFromMuon(evt)
+
+ return StatusCode.Success
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index 79200c089a8c439300df6937ff62d4dee6e144f3..a88fd42c68121ffea75f92c551ee466125359935 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -61,6 +61,12 @@ parser.add_argument("--NoTrackFilt", action='store_true',
parser.add_argument("--no_stable", action='store_true',
help="Don't apply stable beam requirement (default: do)")
+parser.add_argument("--randomTrigFilt", action='store_true',
+ help="Only store events that have a random trigger Events")
+
+parser.add_argument("--randomOnlyTrigFilt", action='store_true',
+ help="Only store events that only have a random trigger Events")
+
parser.add_argument("--unblind", action='store_true',
help="Don't apply signal blinding (default: do)")
parser.add_argument("--onlyblind", action='store_true',
@@ -204,7 +210,7 @@ else:
# Print out what we found
if len(filelist) == 0:
- printf("Found no files for {args.path}!")
+ print("Found no files for {args.path}!")
sys.exit(1)
elif len(filelist) == 1:
print("Processing file:")
@@ -228,6 +234,8 @@ print(f"OnlyBlinded = {args.onlyblind}")
print(f"Stable Beams = {not args.no_stable}")
print(f"Backward = {args.backward}")
print(f"GRL = {args.grl}")
+print(f"Random Filter = {args.randomTrigFilt}")
+print(f"Random Only Filter = {args.randomOnlyTrigFilt}")
# OK, lets run the job here
from AthenaCommon.Logging import log, logging
@@ -317,7 +325,7 @@ if args.isMC:
acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, **mc_kwargs))
else:
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = (not args.NoTrackFilt), DoTrigFilter = args.trigFilt, StableOnly = (not args.no_stable), **grl_kwargs) )
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = (not args.NoTrackFilt), DoTrigFilter = args.trigFilt, StableOnly = (not args.no_stable),DoRandomFilter = args.randomTrigFilt, DoRandomOnlyFilter=args.randomOnlyTrigFilt ,**grl_kwargs))
if not args.verbose:
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index c751d105649810a1b0db80f7b5fadc8d8bf624d6..2842156823b2d6e33d50ac5fd1faa58fb3082b18 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -200,6 +200,8 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_trackCollection.initialize());
ATH_CHECK(m_trackCollectionWithoutIFT.initialize());
ATH_CHECK(m_trackSegmentCollection.initialize());
+ ATH_CHECK(m_trackTrueSegmentCollection.initialize());
+
ATH_CHECK(m_vetoNuContainer.initialize());
ATH_CHECK(m_vetoContainer.initialize());
ATH_CHECK(m_triggerContainer.initialize());
@@ -323,6 +325,16 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("TrackSegment_py", &m_trackseg_py);
m_tree->Branch("TrackSegment_pz", &m_trackseg_pz);
+ addBranch("TrueTrackSegments",&m_nTruetracksegs);
+ m_tree->Branch("TrueTrackSegment_Chi2", &m_Truetrackseg_Chi2);
+ m_tree->Branch("TrueTrackSegment_nDoF", &m_Truetrackseg_DoF);
+ m_tree->Branch("TrueTrackSegment_x", &m_Truetrackseg_x);
+ m_tree->Branch("TrueTrackSegment_y", &m_Truetrackseg_y);
+ m_tree->Branch("TrueTrackSegment_z", &m_Truetrackseg_z);
+ m_tree->Branch("TrueTrackSegment_px", &m_Truetrackseg_px);
+ m_tree->Branch("TrueTrackSegment_py", &m_Truetrackseg_py);
+ m_tree->Branch("TrueTrackSegment_pz", &m_Truetrackseg_pz);
+
//TrackCollection
m_tree->Branch("Track_PropagationError", &m_propagationError, "Track_PropagationError/I");
m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
@@ -731,6 +743,25 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
return StatusCode::SUCCESS;
}
}
+ if (m_doRandomFilter) {
+ bool trig_random = m_tap & 16;
+ if ( !(trig_random)) {
+ // don't process events that fail to activate random trigger
+ ATH_MSG_DEBUG("event did not pass random trigger filter");
+ return StatusCode::SUCCESS;
+ }
+
+ }
+
+ if (m_doRandomOnlyFilter) {
+ bool trig_random = m_tap == 16;
+ if ( !(trig_random)) {
+ // don't process events that fail to activate random only trigger
+ ATH_MSG_DEBUG("event did not pass random trigger filter");
+ return StatusCode::SUCCESS;
+ }
+
+ }
if (m_doScintFilter) {
// filter events, but use waveform peak cuts instead of triggers, as triggers could miss signals slightly out of time
@@ -1505,6 +1536,26 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_trackseg_pz.push_back(SegMomentum.z());
}
+ SG::ReadHandle<TrackCollection> trackTrueSegmentCollection {m_trackTrueSegmentCollection, ctx};
+ ATH_CHECK(trackTrueSegmentCollection.isValid());
+ for (const Trk::Track* trackSeg : *trackTrueSegmentCollection) {
+ if (trackSeg == nullptr) continue;
+ m_nTruetracksegs += 1;
+ m_Truetrackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
+ m_Truetrackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
+ auto SegParameters = trackSeg->trackParameters()->front();
+ const Amg::Vector3D SegPosition = SegParameters->position();
+ const Amg::Vector3D SegMomentum = SegParameters->momentum();
+ m_Truetrackseg_x.push_back(SegPosition.x());
+ m_Truetrackseg_y.push_back(SegPosition.y());
+ m_Truetrackseg_z.push_back(SegPosition.z());
+ m_Truetrackseg_px.push_back(SegMomentum.x());
+ m_Truetrackseg_py.push_back(SegMomentum.y());
+ m_Truetrackseg_pz.push_back(SegMomentum.z());
+ }
+
+
+
// finished processing event, now fill ntuple tree
m_tree->Fill();
m_eventsPassed += 1;
@@ -1602,6 +1653,16 @@ NtupleDumperAlg::clearTree() const
m_trackseg_py.clear();
m_trackseg_pz.clear();
+ m_nTruetracksegs = 0;
+ m_Truetrackseg_Chi2.clear();
+ m_Truetrackseg_DoF.clear();
+ m_Truetrackseg_x.clear();
+ m_Truetrackseg_y.clear();
+ m_Truetrackseg_z.clear();
+ m_Truetrackseg_px.clear();
+ m_Truetrackseg_py.clear();
+ m_Truetrackseg_pz.clear();
+
m_xup.clear();
m_yup.clear();
m_zup.clear();
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 36729773eec08aa8674dfebb43bec17953731c92..c5051ae797e6eee967abf9de69394b88f86d6650 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -83,6 +83,8 @@ private:
SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollection", "Input track collection name" };
SG::ReadHandleKey<TrackCollection> m_trackCollectionWithoutIFT { this, "TrackCollectionWithoutIFT", "CKFTrackCollectionWithoutIFT", "Input track collection name (without IFT)" };
SG::ReadHandleKey<TrackCollection> m_trackSegmentCollection {this, "TrackSegmentCollection", "SegmentFit", "Input track segment collection name"};
+ SG::ReadHandleKey<TrackCollection> m_trackTrueSegmentCollection {this, "TrackTrueSegmentCollection", "Segments", "Input true track segment collection name"};
+
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoNuContainer { this, "VetoNuContainer", "VetoNuWaveformHits", "VetoNu hit container name" };
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoContainer { this, "VetoContainer", "VetoWaveformHits", "Veto hit container name" };
@@ -123,6 +125,9 @@ private:
BooleanProperty m_doScintFilter { this, "DoScintFilter", false, "Only events that pass scintillator coincidence cuts are passed" };
BooleanProperty m_doTrackFilter { this, "DoTrackFilter", true, "Only events that have >= 1 long track are passed, also non-colliding events with a track or calo signal are passed (default)" };
BooleanProperty m_stableOnly { this, "StableOnly", true, "Only events recorded during stable beams are saved (default)" };
+
+ BooleanProperty m_doRandomFilter { this, "DoRandomFilter", false, "Only events also recorded with random triggers are stored." };
+ BooleanProperty m_doRandomOnlyFilter { this, "DoRandomOnlyFilter", false, "Only events recorded with only random triggers are stored." };
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" };
@@ -240,6 +245,16 @@ private:
mutable std::vector<double> m_trackseg_py;
mutable std::vector<double> m_trackseg_pz;
+ mutable unsigned int m_nTruetracksegs;
+ mutable std::vector<double> m_Truetrackseg_Chi2;
+ mutable std::vector<double> m_Truetrackseg_DoF;
+ mutable std::vector<double> m_Truetrackseg_x;
+ mutable std::vector<double> m_Truetrackseg_y;
+ mutable std::vector<double> m_Truetrackseg_z;
+ mutable std::vector<double> m_Truetrackseg_px;
+ mutable std::vector<double> m_Truetrackseg_py;
+ mutable std::vector<double> m_Truetrackseg_pz;
+
mutable int m_longTracks;
mutable int m_propagationError;
mutable std::vector<int> m_module_eta0;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
index 6eab8b1b84f2a745e8d02395f8063306428dbd58..f191ba2c72d04c5668590abcdd320a0846ba41be 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/GhostBusters.cxx
@@ -34,6 +34,7 @@ StatusCode GhostBusters::initialize() {
StatusCode GhostBusters::execute(const EventContext &ctx) const {
m_event_number = ctx.eventID().event_number();
+ ATH_MSG_INFO("GHOSTBUSTING Event Number = " << m_event_number);
SG::WriteHandle outputTrackCollection {m_outputTrackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
@@ -62,6 +63,25 @@ StatusCode GhostBusters::execute(const EventContext &ctx) const {
double meanX = 0.5 * (minX + maxX);
double meanY = 0.5 * (minY + maxY);
double deltaX = deltaY / (2 * std::tan(0.02));
+
+ ATH_MSG_INFO("maxX = " << maxX);
+ ATH_MSG_INFO("maxY = " << maxY);
+ ATH_MSG_INFO("minX = " << minX);
+ ATH_MSG_INFO("minY = " << minY);
+ ATH_MSG_INFO("deltaY = " << deltaY);
+ ATH_MSG_INFO("meanX = " << meanX);
+ ATH_MSG_INFO("meanY = " << meanY);
+ ATH_MSG_INFO("deltaX = " << deltaX);
+
+ ATH_MSG_INFO("yLowerBound = " << meanY - m_yTolerance * deltaY);
+ ATH_MSG_INFO("yUpperBound = " << meanY + m_yTolerance * deltaY);
+ ATH_MSG_INFO("xLeftLowerBound = " << meanX - (1 + m_xTolerance) * deltaX);
+ ATH_MSG_INFO("xLeftUpperBound = " << meanX - (1 - m_xTolerance) * deltaX);
+ ATH_MSG_INFO("xRightLowerBound = " << meanX + (1 - m_xTolerance) * deltaX);
+ ATH_MSG_INFO("xRightUpperBound = " << meanX + (1 + m_xTolerance) * deltaX);
+
+
+
for (Segment &segment : stationSegments) {
bool isGhost = (segment.y() > meanY - m_yTolerance * deltaY) && (segment.y() < meanY + m_yTolerance * deltaY) && (
((segment.x() > meanX - (1 + m_xTolerance) * deltaX) &&
@@ -69,6 +89,10 @@ StatusCode GhostBusters::execute(const EventContext &ctx) const {
((segment.x() > meanX + (1 - m_xTolerance) * deltaX) && (segment.x() < meanX + (1 + m_xTolerance) * deltaX)));
if (isGhost) segment.setGhost();
if (not isGhost && segment.size() >= 5) nGoodSegments++;
+ ATH_MSG_INFO("Segment Y = " << segment.y());
+ ATH_MSG_INFO("Segment X = " << segment.x());
+ ATH_MSG_INFO("Segment isGhost = " << isGhost);
+ ATH_MSG_INFO("nGoodSegments = " << nGoodSegments);
}
for (const Segment &segment : stationSegments) {
m_x = segment.x();
@@ -77,17 +101,31 @@ StatusCode GhostBusters::execute(const EventContext &ctx) const {
m_chi2 = segment.chi2();
m_station = segment.station();
m_size = segment.size();
+ ATH_MSG_INFO("m_x = " << m_x);
+ ATH_MSG_INFO("m_y = " << m_y);
+ ATH_MSG_INFO("m_z = " << m_z);
// m_majorityHits = segment.majorityHits();
m_isGhost = segment.isGhost();
- if (nGoodSegments >= 2 && segment.size() == 4) m_isGhost = true;
+ ATH_MSG_INFO("m_isGhost = " << m_isGhost);
+ ATH_MSG_INFO("nGoodSegments = " << nGoodSegments);
+ ATH_MSG_INFO("segment.size() = " << segment.size());
+ // if (nGoodSegments >= 2 && segment.size() == 4) m_isGhost = true;
+ ATH_MSG_INFO("m_isGhost after if statement = " << m_isGhost);
m_tree->Fill();
if (segment.isGhost()) continue;
- if (nGoodSegments >= 2 && segment.size() == 4) continue;
+ // if (nGoodSegments >= 2 && segment.size() == 4) continue;
outputTracks->push_back(new Trk::Track(*segment.track()));
}
}
ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+
+
+
+
+
+
return StatusCode::SUCCESS;
}