diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index 27b30e25df16dd8baf1c4ce4b82b6a9e6a0f7b6f..8ea26ebb33fef5490110206cab8c940eecfa716b 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -44,8 +44,17 @@ parser.add_argument("--isMC", action='store_true',
parser.add_argument("--partial", action='store_true',
help="Allow partial input files")
+parser.add_argument("--trigFilt", action='store_true',
+ help="apply trigger event filter")
+parser.add_argument("--scintFilt", action='store_true',
+ help="apply scintillator event filter")
+parser.add_argument("--NoTrackFilt", action='store_true',
+ help="Don't apply track event filter (default: do)")
+
parser.add_argument("--unblind", action='store_true',
help="Don't apply signal blinding (default: do)")
+parser.add_argument("--onlyblind", action='store_true',
+ help="Only store events that were blinded (will override unblind arg)")
parser.add_argument("--fluka", action='store_true',
help="Add FLUKA weights to ntuple")
@@ -138,8 +147,12 @@ if filepath.is_dir():
print(f"Run = {runstr}")
print(f"First = {firstseg}")
print(f"Last = {lastseg}")
- print(f"Args = {args.tag}")
+ print(f"Tag = {args.tag}")
+ print(f"Trigger Filter = {args.trigFilt}")
+ print(f"Scintillator Filter = {args.scintFilt}")
+ print(f"Track Filter = {not args.NoTrackFilt}")
print(f"Blind = {not args.unblind}")
+ print(f"OnlyBlinded = {args.onlyblind}")
# Find any tags
tagstr = firststem.replace(f"{firstfaser}-{firstshort}-{runstr}-{firstseg}", "")
@@ -243,7 +256,7 @@ if args.isMC:
acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile))
else:
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind)))
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = (not args.NoTrackFilt), DoTrigFilter = args.trigFilt))
if not args.verbose:
from AthenaConfiguration.ComponentFactory import CompFactory
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 3f476edb40480e5871a2df39a3cd2695ea894e4c..c06051d2a5d9ee23b527fff81ec78e1ae47a3f20 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -18,9 +18,6 @@
#include <TH1F.h>
#include <numeric>
-constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
-
-
NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
ISvcLocator *pSvcLocator)
: AthReentrantAlgorithm(name, pSvcLocator),
@@ -167,10 +164,13 @@ StatusCode NtupleDumperAlg::initialize()
addWaveBranches("Preshower",2,12);
addWaveBranches("Calo",4,0);
+ m_tree->Branch("ScintHit", &m_scintHit);
+
addCalibratedBranches("Calo",4,0);
addBranch("Calo_total_nMIP", &m_calo_total_nMIP);
addBranch("Calo_total_E_dep", &m_calo_total_E_dep);
addBranch("Calo_total_E_EM", &m_calo_total_E_EM);
+ addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged);
addCalibratedBranches("Preshower",2,12);
addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP);
@@ -261,6 +261,11 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("Track_ThetaX_atCalo", &m_thetaxCalo);
m_tree->Branch("Track_ThetaY_atCalo", &m_thetayCalo);
+ m_tree->Branch("Track_x_atMaxRadius", &m_x_atMaxRadius);
+ m_tree->Branch("Track_y_atMaxRadius", &m_y_atMaxRadius);
+ m_tree->Branch("Track_z_atMaxRadius", &m_z_atMaxRadius);
+ m_tree->Branch("Track_r_atMaxRadius", &m_r_atMaxRadius);
+
//TRUTH
m_tree->Branch("t_pdg", &m_t_pdg);
m_tree->Branch("t_barcode", &m_t_barcode);
@@ -385,10 +390,10 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14]));
- m_MIP_sim_Edep_calo = 0.0585; // MIP deposits 0.0585 GeV of energy in calo
- m_MIP_sim_Edep_preshower = 0.004894; // MIP deposits 0.004894 GeV of energy in a preshower layer
-
- if (m_doBlinding) {
+ if (m_onlyBlinded){
+ ATH_MSG_INFO("Only events that would be blinded are saved in ntuple");
+ m_doBlinding = false;
+ } else if (m_doBlinding) {
ATH_MSG_INFO("Blinding will be enforced for real data.");
} else {
ATH_MSG_INFO("Blinding will NOT be enforced for real data.");
@@ -410,85 +415,97 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
isMC = true;
}
- // if real data, store charge in histograms from random events and only fill ntuple from coincidence events
+ // if real data, correct for diigitzer timing jitter with clock phase
if (!isMC) {
- SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx);
- m_tap=triggerData->tap();
- // unpack trigger word: 1=calo, 2=veotnu|veto1|preshower, 4=TimingLayer, 8=(VetoNu|Veto2)&Preshower, 16=random, 32=LED
- bool trig_random = false;
- if ( m_tap == 16 ) trig_random = true;
+ // correct waveform time with clock phase
+ // needs to be done before processing waveform values
+ SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_ClockWaveformContainer, ctx);
+ ATH_CHECK(clockHandle.isValid());
- bool trig_coincidence_preshower_and_vetoes = false;
- if ( (m_tap&8) != 0 ) trig_coincidence_preshower_and_vetoes = true;
+ if (clockHandle->phase() < -2.0) { // wrap around clock pahse so -pi goes to pi
+ m_clock_phase = ((clockHandle->phase() + 2*3.14159) / 3.14159) * 12.5;
+ } else {
+ m_clock_phase = (clockHandle->phase() / 3.14159) * 12.5;
+ }
+ }
- bool trig_coincidence_timing_and_vetoesORpreshower = false;
- if ( ((m_tap&4)!=0) && ((m_tap&2)!=0) ) trig_coincidence_timing_and_vetoesORpreshower = true;
+ // process all waveform data for all scintillator and calorimeter channels
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
+ ATH_CHECK(vetoNuContainer.isValid());
- bool trig_coincidence_timing_and_calo = false;
- if ( ((m_tap&4)!=0) && ((m_tap&1)!=0) ) trig_coincidence_timing_and_calo = true;
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
+ ATH_CHECK(vetoContainer.isValid());
- bool trig_coincidence_vetoesORpreshower_and_calo = false;
- if ( ((m_tap&2)!=0) && ((m_tap&1)!=0) ) trig_coincidence_vetoesORpreshower_and_calo = true;
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ ATH_CHECK(triggerContainer.isValid());
- bool trig_calo = (m_tap & 1);
+ SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
+ ATH_CHECK(preshowerContainer.isValid());
- // for random trigger, store charge of scintillators in histograms
- if (trig_random) {
- // Read in Waveform containers
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
- ATH_CHECK(vetoNuContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
+ ATH_CHECK(ecalContainer.isValid());
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
- ATH_CHECK(vetoContainer.isValid());
+ FillWaveBranches(*vetoNuContainer);
+ FillWaveBranches(*vetoContainer);
+ FillWaveBranches(*triggerContainer);
+ FillWaveBranches(*preshowerContainer);
+ FillWaveBranches(*ecalContainer);
- SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
- ATH_CHECK(triggerContainer.isValid());
+ // if real data, store charge in histograms from random events
+ if (!isMC) {
+ SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx);
+ m_tap=triggerData->tap();
+ // for random trigger, store charge of scintillators in histograms
+ if ((m_tap&16) != 0) {
+ // Fill histograms
+ for (int chan = 0; chan<15; chan++) {
+ m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]);
+ }
+ return StatusCode::SUCCESS; // finished with this randomly triggered event
+ }
- SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
- ATH_CHECK(preshowerContainer.isValid());
+ if (m_doTrigFilter) {
- SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
- ATH_CHECK(ecalContainer.isValid());
+ bool trig_coincidence_preshower_and_vetoes = ( (m_tap&8) != 0 );
+ bool trig_coincidence_timing_and_vetoesORpreshower = ( ((m_tap&4)!=0) && ((m_tap&2)!=0) );
+ bool trig_calo = (m_tap & 1);
- // Fill histograms
- if (vetoNuContainer.isValid()) {
- for (auto hit : *vetoNuContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (vetoContainer.isValid()) {
- for (auto hit : *vetoContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (triggerContainer.isValid()) {
- for (auto hit : *triggerContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (preshowerContainer.isValid()) {
- for (auto hit : *preshowerContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (ecalContainer.isValid()) {
- for (auto hit : *ecalContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
+ if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_calo) ) {
+ // don't process events that fail to activate coincidence triggers
+ ATH_MSG_DEBUG("event did not pass trigger filter");
+ return StatusCode::SUCCESS;
}
+ }
- return StatusCode::SUCCESS; // finished with this randomly triiggered event
-
- // 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;
+ if (m_doScintFilter) {
+ // filter events, but use waveform peak cuts instead of triggers, as triggers could miss signals slightly out of time
+ // digitizer channels described here: https://twiki.cern.ch/twiki/bin/viewauth/FASER/CaloScintillator
+ bool calo_trig = (m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0);
+ bool vetoNu_trig = (m_wave_raw_peak[4] > 25.0 && m_wave_raw_peak[5] > 25.0);
+ bool vetoSt2_trig = (m_wave_raw_peak[6] > 25.0 && m_wave_raw_peak[7] > 25.0);
+ bool timingScint_trig = ((m_wave_raw_peak[8] > 25.0 && m_wave_raw_peak[9] > 25.0) || (m_wave_raw_peak[10] > 25.0 && m_wave_raw_peak[11] > 25.0));
+ bool preshower_trig = (m_wave_raw_peak[12] > 3.0 && m_wave_raw_peak[13] > 3.0);
+ bool vetoSt1_trig = (m_wave_raw_peak[14] > 25.0);
+
+ bool veto_OR_trig = (vetoNu_trig || vetoSt1_trig || vetoSt2_trig);
+
+ bool passes_ScintFilter = false;
+ if (calo_trig) {
+ passes_ScintFilter = true;
+ } else if (veto_OR_trig && timingScint_trig) {
+ passes_ScintFilter = true;
+ } else if (veto_OR_trig && preshower_trig) {
+ passes_ScintFilter = true;
+ } else if (timingScint_trig && preshower_trig) {
+ passes_ScintFilter = true;
+ }
+
+ if (!passes_ScintFilter) {
+ ATH_MSG_DEBUG("event did not pass scint filter");
+ return StatusCode::SUCCESS; // only store events that pass filter
+ }
}
+
// store trigger data in ntuple variables
m_tbp=triggerData->tbp();
m_tap=triggerData->tap();
@@ -530,18 +547,34 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("LHC data distanceToTrainStart = " << lhcData->distanceToTrainStart() );
ATH_MSG_DEBUG("LHC data distanceToPreviousColliding = " << lhcData->distanceToPreviousColliding() );
- // correct waveform time with clock phase
- SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_ClockWaveformContainer, ctx);
- ATH_CHECK(clockHandle.isValid());
-
- if (clockHandle->phase() < -2.0) { // wrap around clock pahse so -pi goes to pi
- m_clock_phase = ((clockHandle->phase() + 2*3.14159) / 3.14159) * 12.5;
- } else {
- m_clock_phase = (clockHandle->phase() / 3.14159) * 12.5;
- }
-
} // done with processing only on real data
+ // fill scintHit word with bits that reflect if a scintillator was hit (1 = vetoNu0, 2 = vetoNu1, 4 = vetoSt1_1, 8 vetoSt2_0, 16 = vetoSt2_1, 32 = Timing scint, 64 = preshower0, 128 = preshower1)
+ if (m_wave_raw_charge[4] > 40.0) {
+ m_scintHit = m_scintHit | 1;
+ }
+ if (m_wave_raw_charge[5] > 40.0) {
+ m_scintHit = m_scintHit | 2;
+ }
+ if (m_wave_raw_charge[14] > 40.0) {
+ m_scintHit = m_scintHit | 4;
+ }
+ if (m_wave_raw_charge[6] > 40.0) {
+ m_scintHit = m_scintHit | 8;
+ }
+ if (m_wave_raw_charge[7] > 40.0) {
+ m_scintHit = m_scintHit | 16;
+ }
+ if (m_wave_raw_charge[8]+m_wave_raw_charge[9]+m_wave_raw_charge[10]+m_wave_raw_charge[11] > 25.0) {
+ m_scintHit = m_scintHit | 32;
+ }
+ if (m_wave_raw_charge[12] > 2.5) {
+ m_scintHit = m_scintHit | 64;
+ }
+ if (m_wave_raw_charge[13] > 2.5) {
+ m_scintHit = m_scintHit | 128;
+ }
+
m_run_number = ctx.eventID().run_number();
m_event_number = ctx.eventID().event_number();
m_event_time = ctx.eventID().time_stamp();
@@ -574,7 +607,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
for (auto particle : *truthParticleContainer)
{
- // loop over first 10 truth particles (for non A' samples)
+ // loop over first 10 truth particles (for non A' samples)
if (ipart++ < 10) {
@@ -590,9 +623,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
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);
+ m_truth_prod_x.push_back(999999);
+ m_truth_prod_y.push_back(999999);
+ m_truth_prod_z.push_back(999999);
}
if ( particle->hasDecayVtx()) {
@@ -600,9 +633,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
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);
+ m_truth_dec_x.push_back(999999);
+ m_truth_dec_y.push_back(999999);
+ m_truth_dec_z.push_back(999999);
}
}
@@ -621,9 +654,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
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);
+ m_truthM_x.push_back(999999);
+ m_truthM_y.push_back(999999);
+ m_truthM_z.push_back(999999);
}
}
@@ -639,9 +672,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_truthd0_y.push_back(particle->prodVtx()->y());
m_truthd0_z.push_back(particle->prodVtx()->z());
} else {
- m_truthd0_x.push_back(NaN);
- m_truthd0_y.push_back(NaN);
- m_truthd0_z.push_back(NaN);
+ m_truthd0_x.push_back(999999);
+ m_truthd0_y.push_back(999999);
+ m_truthd0_z.push_back(999999);
}
}
if ( particle->pdgId() == -11) // daughter particle (electron)
@@ -656,9 +689,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_truthd1_y.push_back(particle->prodVtx()->y());
m_truthd1_z.push_back(particle->prodVtx()->z());
} else {
- m_truthd1_x.push_back(NaN);
- m_truthd1_y.push_back(NaN);
- m_truthd1_z.push_back(NaN);
+ m_truthd1_x.push_back(999999);
+ m_truthd1_y.push_back(999999);
+ m_truthd1_z.push_back(999999);
}
}
}
@@ -667,7 +700,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
// load in calibrated calo container
- if (m_doCalib) {
SG::ReadHandle<xAOD::CalorimeterHitContainer> ecalCalibratedContainer { m_ecalCalibratedContainer, ctx };
ATH_CHECK(ecalCalibratedContainer.isValid());
for (auto hit : *ecalCalibratedContainer) {
@@ -694,6 +726,12 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
//}
}
+ // add a fudged energy variable that corrects the MC to agree with the testbeam results
+ m_calo_total_E_EM_fudged = m_calo_total_E_EM;
+ if (isMC) {
+ m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
+ }
+
// load in calibrated preshower container
SG::ReadHandle<xAOD::CalorimeterHitContainer> preshowerCalibratedContainer { m_preshowerCalibratedContainer, ctx };
ATH_CHECK(preshowerCalibratedContainer.isValid());
@@ -707,40 +745,24 @@ 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 };
- ATH_CHECK(vetoNuContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
- ATH_CHECK(vetoContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
- ATH_CHECK(triggerContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
- ATH_CHECK(preshowerContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
- ATH_CHECK(ecalContainer.isValid());
-
- FillWaveBranches(*vetoNuContainer);
- FillWaveBranches(*vetoContainer);
- FillWaveBranches(*triggerContainer);
- FillWaveBranches(*preshowerContainer);
- 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
// Only blind colliding BCIDs (+/- 1)
- if ( (!isMC) && m_doBlinding && abs(m_distanceToCollidingBCID) <= 1) {
+ bool blinded_event = false;
+ if ((!isMC) && 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;
+ if (m_wave_raw_charge[4] < 40.0 and m_wave_raw_charge[5] < 40.0 and m_wave_raw_charge[6] < 40.0 and m_wave_raw_charge[7] < 40.0 and m_wave_raw_charge[14] < 40.0) { // channels 4 and 5 are vetoNu, channels 6,7, and 14 are veto
+ blinded_event = true;
}
}
}
+ if (blinded_event && m_doBlinding) {
+ return StatusCode::SUCCESS;
+ } else if (!blinded_event && m_onlyBlinded) {
+ 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;
@@ -753,73 +775,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
auto gctx = faserGeometryContext.context();
- // loop over clusters and store how many clusters are in each tracking station
- SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
- ATH_CHECK(clusterContainer.isValid());
-
- for (auto collection : *clusterContainer)
- {
- Identifier id = collection->identify();
- int station = m_sctHelper->station(id);
- int clusters = (int) collection->size();
- switch (station)
- {
- case 0:
- m_station0Clusters += clusters;
- // following lines commented out depict how to access cluster position
- //for (auto cluster : *collection) {
- // if (cluster == nullptr) continue;
- // auto pos = cluster->globalPosition();
- // m_station0ClusterX.push_back(pos.x());
- //}
- break;
- case 1:
- m_station1Clusters += clusters;
- break;
- case 2:
- m_station2Clusters += clusters;
- break;
- case 3:
- m_station3Clusters += clusters;
- break;
- default:
- ATH_MSG_FATAL("Unknown tracker station number " << station);
- break;
- }
- }
-
- // loop over spacepoints and store each space point position
- SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
- ATH_CHECK(spacePointContainer.isValid());
- for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
- m_nspacepoints += spacePointCollection->size();
- for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
- auto pos = spacePoint->globalPosition();
- m_spacepointX.push_back(pos.x());
- m_spacepointY.push_back(pos.y());
- m_spacepointZ.push_back(pos.z());
- }
- }
-
- // loop over track segments and store position, momentum, chi2, and nDOF for each segment
- SG::ReadHandle<TrackCollection> trackSegmentCollection {m_trackSegmentCollection, ctx};
- ATH_CHECK(trackSegmentCollection.isValid());
- for (const Trk::Track* trackSeg : *trackSegmentCollection) {
- if (trackSeg == nullptr) continue;
- m_ntracksegs += 1;
- m_trackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
- m_trackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
- auto SegParameters = trackSeg->trackParameters()->front();
- const Amg::Vector3D SegPosition = SegParameters->position();
- const Amg::Vector3D SegMomentum = SegParameters->momentum();
- m_trackseg_x.push_back(SegPosition.x());
- m_trackseg_y.push_back(SegPosition.y());
- m_trackseg_z.push_back(SegPosition.z());
- m_trackseg_px.push_back(SegMomentum.x());
- m_trackseg_py.push_back(SegMomentum.y());
- m_trackseg_pz.push_back(SegMomentum.z());
- }
-
// Write out all truth particle barcodes that have a momentum larger than MinMomentum (default is 50 GeV)
std::map<int, size_t> truthParticleCount {};
if (isMC) {
@@ -836,7 +791,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// extrapolate track to all scintillator positions and store extrapolated position and angle
SG::ReadHandle<TrackCollection> trackCollection;
if (m_useIFT) {
- SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that excludes IFT
+ SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that includes IFT
trackCollection = tc;
} else {
SG::ReadHandle<TrackCollection> tc {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
@@ -896,6 +851,20 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_pzdown.push_back(downstreamParameters->momentum().z());
m_pdown.push_back(sqrt( pow(downstreamParameters->momentum().x(),2) + pow(downstreamParameters->momentum().y(),2) + pow(downstreamParameters->momentum().z(),2) ));
+ // find and store the position of the measurement on track with the largest radius
+ double maxRadius = 0;
+ Amg::Vector3D positionAtMaxRadius {};
+ for (const Trk::TrackParameters* trackParameters : *track->trackParameters()) {
+ if (radius(trackParameters->position()) > maxRadius) {
+ maxRadius = radius(trackParameters->position());
+ positionAtMaxRadius = trackParameters->position();
+ }
+ }
+ m_r_atMaxRadius.push_back(maxRadius);
+ m_x_atMaxRadius.push_back(positionAtMaxRadius.x());
+ m_y_atMaxRadius.push_back(positionAtMaxRadius.y());
+ m_z_atMaxRadius.push_back(positionAtMaxRadius.z());
+
if (isMC) { // if simulation, store track truth info as well
auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
if (truthParticle != nullptr) {
@@ -917,18 +886,18 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
m_t_prodVtx_z.push_back(truthParticle->prodVtx()->z());
} else {
- m_t_prodVtx_x.push_back(NaN);
- m_t_prodVtx_y.push_back(NaN);
- m_t_prodVtx_z.push_back(NaN);
+ m_t_prodVtx_x.push_back(999999);
+ m_t_prodVtx_y.push_back(999999);
+ m_t_prodVtx_z.push_back(999999);
}
if (truthParticle->hasDecayVtx()) {
m_t_decayVtx_x.push_back(truthParticle->decayVtx()->x());
m_t_decayVtx_y.push_back(truthParticle->decayVtx()->y());
m_t_decayVtx_z.push_back(truthParticle->decayVtx()->z());
} else {
- m_t_decayVtx_x.push_back(NaN);
- m_t_decayVtx_y.push_back(NaN);
- m_t_decayVtx_z.push_back(NaN);
+ m_t_decayVtx_x.push_back(999999);
+ m_t_decayVtx_y.push_back(999999);
+ m_t_decayVtx_z.push_back(999999);
}
m_t_px.push_back(truthParticle->px());
m_t_py.push_back(truthParticle->py());
@@ -940,167 +909,175 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_t_eta.push_back(truthParticle->p4().Eta());
} else {
ATH_MSG_WARNING("Can not find truthParticle.");
- setNaN();
+ clearTrackTruth();
}
} else {
- setNaN();
+ clearTrackTruth();
}
- // fill extrapolation vectors with NaN, will get set to real number if the track extrapolation succeeds
- m_xVetoNu.push_back(NaN);
- m_yVetoNu.push_back(NaN);
- m_thetaxVetoNu.push_back(NaN);
- m_thetayVetoNu.push_back(NaN);
- m_xVetoStation1.push_back(NaN);
- m_yVetoStation1.push_back(NaN);
- m_thetaxVetoStation1.push_back(NaN);
- m_thetayVetoStation1.push_back(NaN);
- m_xVetoStation2.push_back(NaN);
- m_yVetoStation2.push_back(NaN);
- m_thetaxVetoStation2.push_back(NaN);
- m_thetayVetoStation2.push_back(NaN);
- m_xTrig.push_back(NaN);
- m_yTrig.push_back(NaN);
- m_thetaxTrig.push_back(NaN);
- m_thetayTrig.push_back(NaN);
- m_xPreshower1.push_back(NaN);
- m_yPreshower1.push_back(NaN);
- m_thetaxPreshower1.push_back(NaN);
- m_thetayPreshower1.push_back(NaN);
- m_xPreshower2.push_back(NaN);
- m_yPreshower2.push_back(NaN);
- m_thetaxPreshower2.push_back(NaN);
- m_thetayPreshower2.push_back(NaN);
- m_xCalo.push_back(NaN);
- m_yCalo.push_back(NaN);
- m_thetaxCalo.push_back(NaN);
- m_thetayCalo.push_back(NaN);
-
- // extrapolate track from first station
- if (stationMap.count(1) > 0) { // extrapolation crashes if the track parameters are is too far away to extrapolate
- Amg::Vector3D position = upstreamParameters->position();
- Amg::Vector3D momentum = upstreamParameters->momentum();
- Acts::BoundVector params = Acts::BoundVector::Zero();
- params[Acts::eBoundLoc0] = -position.y();
- params[Acts::eBoundLoc1] = position.x();
- params[Acts::eBoundPhi] = momentum.phi();
- params[Acts::eBoundTheta] = momentum.theta();
- params[Acts::eBoundQOverP] = upstreamParameters->charge() / momentum.mag();
- params[Acts::eBoundTime] = 0;
- auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParameters(std::move(startSurface), params, upstreamParameters->charge());
-
- auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_VetoNu, Acts::backward);
- if (targetParameters_VetoNu != nullptr) {
- auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
- auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
- m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
- m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
- m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
- m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
- } else {
- ATH_MSG_INFO("vetoNu null targetParameters");
- }
-
- auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto1, Acts::backward);
- if (targetParameters_Veto1 != nullptr) {
- auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
- auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
- m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
- m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
- m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
- m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
- } else {
- ATH_MSG_INFO("veto1 null targetParameters");
- }
+ // fill extrapolation vectors with dummy values, will get set to real number if the track extrapolation succeeds
+ m_xVetoNu.push_back(999999);
+ m_yVetoNu.push_back(999999);
+ m_thetaxVetoNu.push_back(999999);
+ m_thetayVetoNu.push_back(999999);
+ m_xVetoStation1.push_back(999999);
+ m_yVetoStation1.push_back(999999);
+ m_thetaxVetoStation1.push_back(999999);
+ m_thetayVetoStation1.push_back(999999);
+ m_xVetoStation2.push_back(999999);
+ m_yVetoStation2.push_back(999999);
+ m_thetaxVetoStation2.push_back(999999);
+ m_thetayVetoStation2.push_back(999999);
+ m_xTrig.push_back(999999);
+ m_yTrig.push_back(999999);
+ m_thetaxTrig.push_back(999999);
+ m_thetayTrig.push_back(999999);
+ m_xPreshower1.push_back(999999);
+ m_yPreshower1.push_back(999999);
+ m_thetaxPreshower1.push_back(999999);
+ m_thetayPreshower1.push_back(999999);
+ m_xPreshower2.push_back(999999);
+ m_yPreshower2.push_back(999999);
+ m_thetaxPreshower2.push_back(999999);
+ m_thetayPreshower2.push_back(999999);
+ m_xCalo.push_back(999999);
+ m_yCalo.push_back(999999);
+ m_thetaxCalo.push_back(999999);
+ m_thetayCalo.push_back(999999);
+
+ // Define paramters for track extrapolation from most upstream measurement
+ Amg::Vector3D position_up = upstreamParameters->position();
+ Amg::Vector3D momentum_up = upstreamParameters->momentum();
+ Acts::BoundVector params_up = Acts::BoundVector::Zero();
+ params_up[Acts::eBoundLoc0] = -position_up.y();
+ params_up[Acts::eBoundLoc1] = position_up.x();
+ params_up[Acts::eBoundPhi] = momentum_up.phi();
+ params_up[Acts::eBoundTheta] = momentum_up.theta();
+ params_up[Acts::eBoundQOverP] = upstreamParameters->charge() / momentum_up.mag();
+ params_up[Acts::eBoundTime] = 0;
+ auto startSurface_up = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_up.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_up(std::move(startSurface_up), params_up, upstreamParameters->charge());
+
+ // Define paramters for track extrapolation from most downstream measurement
+ Amg::Vector3D position_down = downstreamParameters->position();
+ Amg::Vector3D momentum_down = downstreamParameters->momentum();
+ Acts::BoundVector params_down = Acts::BoundVector::Zero();
+ params_down[Acts::eBoundLoc0] = -position_down.y();
+ params_down[Acts::eBoundLoc1] = position_down.x();
+ params_down[Acts::eBoundPhi] = momentum_down.phi();
+ params_down[Acts::eBoundTheta] = momentum_down.theta();
+ params_down[Acts::eBoundQOverP] = downstreamParameters->charge() / momentum_down.mag();
+ params_down[Acts::eBoundTime] = 0;
+ auto startSurface_down = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_down.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_down(std::move(startSurface_down), params_down, downstreamParameters->charge());
+
+ // Extrapolate track to scintillators
+ auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_VetoNu, Acts::backward);
+ if (targetParameters_VetoNu != nullptr) {
+ auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
+ auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
+ m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
+ m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
+ m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
+ m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
+ } else {
+ ATH_MSG_INFO("vetoNu null targetParameters");
+ }
- auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto2, Acts::backward);
- if (targetParameters_Veto2 != nullptr) {
- auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
- auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
- m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
- m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
- m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
- m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
- } else {
- ATH_MSG_INFO("veto2 null targetParameters");
- }
+ auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto1, Acts::backward);
+ if (targetParameters_Veto1 != nullptr) {
+ auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
+ auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
+ m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
+ m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
+ m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
+ m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
+ } else {
+ ATH_MSG_INFO("veto1 null targetParameters");
+ }
- auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Trig, Acts::backward); // must extrapolate backsards to trig plane if track starts in station 1
- if (targetParameters_Trig != nullptr) {
- auto targetPosition_Trig = targetParameters_Trig->position(gctx);
- auto targetMomentum_Trig = targetParameters_Trig->momentum();
- m_xTrig[m_longTracks] = targetPosition_Trig.x();
- m_yTrig[m_longTracks] = targetPosition_Trig.y();
- m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
- m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
- } else {
- ATH_MSG_INFO("Trig null targetParameters");
- }
+ auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto2, Acts::backward);
+ if (targetParameters_Veto2 != nullptr) {
+ auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
+ auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
+ m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
+ m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
+ m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
+ m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
+ } else {
+ ATH_MSG_INFO("veto2 null targetParameters");
+ }
+ auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Trig, Acts::backward); // must extrapolate backsards to trig plane if track starts in station 1
+ if (targetParameters_Trig != nullptr) {
+ auto targetPosition_Trig = targetParameters_Trig->position(gctx);
+ auto targetMomentum_Trig = targetParameters_Trig->momentum();
+ m_xTrig[m_longTracks] = targetPosition_Trig.x();
+ m_yTrig[m_longTracks] = targetPosition_Trig.y();
+ m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
+ m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
+ } else {
+ ATH_MSG_INFO("Trig null targetParameters");
}
- // extrapolate track from tracking station 3
- if (stationMap.count(3) > 0) { // extrapolation crashes if the track does not end in the Station 3, as it is too far away to extrapolate
- Amg::Vector3D positionDown = downstreamParameters->position();
- Amg::Vector3D momentumDown = downstreamParameters->momentum();
- Acts::BoundVector paramsDown = Acts::BoundVector::Zero();
- paramsDown[Acts::eBoundLoc0] = -positionDown.y();
- paramsDown[Acts::eBoundLoc1] = positionDown.x();
- paramsDown[Acts::eBoundPhi] = momentumDown.phi();
- paramsDown[Acts::eBoundTheta] = momentumDown.theta();
- paramsDown[Acts::eBoundQOverP] = downstreamParameters->charge() / momentumDown.mag();
- paramsDown[Acts::eBoundTime] = 0;
- auto startSurfaceDown = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, positionDown.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParametersDown(std::move(startSurfaceDown), paramsDown, downstreamParameters->charge());
-
- auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower1, Acts::forward);
- if (targetParameters_Preshower1 != nullptr) {
- auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
- auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
- m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
- m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
- m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
- m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
- } else {
- ATH_MSG_INFO("Preshower1 null targetParameters");
- }
+ auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower1, Acts::forward);
+ if (targetParameters_Preshower1 != nullptr) {
+ auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
+ auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
+ m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
+ m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
+ m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
+ m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
+ } else {
+ ATH_MSG_INFO("Preshower1 null targetParameters");
+ }
- auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower2, Acts::forward);
- if (targetParameters_Preshower2 != nullptr) {
- auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
- auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
- m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
- m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
- m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
- m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
- } else {
- ATH_MSG_INFO("Preshower2 null targetParameters");
- }
+ auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower2, Acts::forward);
+ if (targetParameters_Preshower2 != nullptr) {
+ auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
+ auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
+ m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
+ m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
+ m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
+ m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
+ } else {
+ ATH_MSG_INFO("Preshower2 null targetParameters");
+ }
- auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Calo, Acts::forward);
- if (targetParameters_Calo != nullptr) {
- auto targetPosition_Calo = targetParameters_Calo->position(gctx);
- auto targetMomentum_Calo = targetParameters_Calo->momentum();
- m_xCalo[m_longTracks] = targetPosition_Calo.x();
- m_yCalo[m_longTracks] = targetPosition_Calo.y();
- m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
- m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
- } else {
- ATH_MSG_INFO("Calo null targetParameters");
- }
+ auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Calo, Acts::forward);
+ if (targetParameters_Calo != nullptr) {
+ auto targetPosition_Calo = targetParameters_Calo->position(gctx);
+ auto targetMomentum_Calo = targetParameters_Calo->momentum();
+ m_xCalo[m_longTracks] = targetPosition_Calo.x();
+ m_yCalo[m_longTracks] = targetPosition_Calo.y();
+ m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
+ m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
+ } else {
+ ATH_MSG_INFO("Calo null targetParameters");
}
m_longTracks++;
}
+ if (!isMC) {
+ if (m_doTrackFilter) { // filter events: colliding bunches have at least one long track, non-colliding bunches have at least one long track or one calo module with raw_peak > 3 mV
+ if (m_longTracks == 0 && abs(m_distanceToCollidingBCID) <= 1) {
+ return StatusCode::SUCCESS;
+ } else if (abs(m_distanceToCollidingBCID) > 1) {
+ if (!(m_longTracks > 0 || m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0 )) {
+ return StatusCode::SUCCESS;
+ }
+ }
+ }
+ }
+
if (isMC) {
for (auto &tp : truthParticleCount) {
m_truthParticleBarcode.push_back(tp.first);
@@ -1109,40 +1086,97 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
}
+ // loop over clusters and store how many clusters are in each tracking station
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+ ATH_CHECK(clusterContainer.isValid());
+
+ for (auto collection : *clusterContainer)
+ {
+ Identifier id = collection->identify();
+ int station = m_sctHelper->station(id);
+ int clusters = (int) collection->size();
+ switch (station)
+ {
+ case 0:
+ m_station0Clusters += clusters;
+ // following lines commented out depict how to access cluster position
+ //for (auto cluster : *collection) {
+ // if (cluster == nullptr) continue;
+ // auto pos = cluster->globalPosition();
+ // m_station0ClusterX.push_back(pos.x());
+ //}
+ break;
+ case 1:
+ m_station1Clusters += clusters;
+ break;
+ case 2:
+ m_station2Clusters += clusters;
+ break;
+ case 3:
+ m_station3Clusters += clusters;
+ break;
+ default:
+ ATH_MSG_FATAL("Unknown tracker station number " << station);
+ break;
+ }
+ }
+
+ // loop over spacepoints and store each space point position
+ SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
+ ATH_CHECK(spacePointContainer.isValid());
+ for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+ m_nspacepoints += spacePointCollection->size();
+ for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
+ auto pos = spacePoint->globalPosition();
+ m_spacepointX.push_back(pos.x());
+ m_spacepointY.push_back(pos.y());
+ m_spacepointZ.push_back(pos.z());
+ }
+ }
+
+ // loop over track segments and store position, momentum, chi2, and nDOF for each segment
+ SG::ReadHandle<TrackCollection> trackSegmentCollection {m_trackSegmentCollection, ctx};
+ ATH_CHECK(trackSegmentCollection.isValid());
+ for (const Trk::Track* trackSeg : *trackSegmentCollection) {
+ if (trackSeg == nullptr) continue;
+ m_ntracksegs += 1;
+ m_trackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
+ m_trackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
+ auto SegParameters = trackSeg->trackParameters()->front();
+ const Amg::Vector3D SegPosition = SegParameters->position();
+ const Amg::Vector3D SegMomentum = SegParameters->momentum();
+ m_trackseg_x.push_back(SegPosition.x());
+ m_trackseg_y.push_back(SegPosition.y());
+ m_trackseg_z.push_back(SegPosition.z());
+ m_trackseg_px.push_back(SegMomentum.x());
+ m_trackseg_py.push_back(SegMomentum.y());
+ m_trackseg_pz.push_back(SegMomentum.z());
+ }
+
// finished processing event, now fill ntuple tree
m_tree->Fill();
m_eventsPassed += 1;
return StatusCode::SUCCESS;
}
-
StatusCode NtupleDumperAlg::finalize()
{
ATH_MSG_INFO("Number of events passed Ntuple selectioon = " << m_eventsPassed);
return StatusCode::SUCCESS;
}
-bool NtupleDumperAlg::waveformHitOK(const xAOD::WaveformHit* hit) const
-{
- if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED) || hit->status_bit(xAOD::WaveformStatus::SECONDARY)) return false;
- return true;
-}
-
void
NtupleDumperAlg::clearTree() const
{
- // set all float variables to NaN
- // set all int variables to 999999 (can't set int to NaN, because NaN is a double)
- // set all counter variables to zero
- // set all trigger words to zero
- m_run_number = 999999;
- m_event_number = 999999;
- m_event_time = 999999;
- m_bcid = 999999;
-
- m_fillNumber = 999999;
- m_betaStar = NaN;
- m_crossingAngle = NaN;
+ // don't use NaN, they are annoying when doing cuts
+ m_run_number = 0;
+ m_event_number = 0;
+ m_event_time = 0;
+ m_bcid = 0;
+
+ m_fillNumber = 0;
+ m_betaStar = 0;
+ m_crossingAngle = 0;
m_distanceToCollidingBCID = 999999;
m_distanceToUnpairedB1 = 999999;
m_distanceToUnpairedB2 = 999999;
@@ -1156,25 +1190,28 @@ NtupleDumperAlg::clearTree() const
m_inputBitsNext=0;
for(int ii=0;ii<15;ii++) {
- m_wave_localtime[ii]=NaN;
- m_wave_peak[ii]=NaN;
- m_wave_width[ii]=NaN;
- m_wave_charge[ii]=NaN;
-
- m_wave_raw_peak[ii]=NaN;
- m_wave_raw_charge[ii]=NaN;
- m_wave_baseline_mean[ii]=NaN;
- m_wave_baseline_rms[ii]=NaN;
+ m_wave_localtime[ii]=0;
+ m_wave_peak[ii]=0;
+ m_wave_width[ii]=0;
+ m_wave_charge[ii]=0;
+
+ m_wave_raw_peak[ii]=0;
+ m_wave_raw_charge[ii]=0;
+ m_wave_baseline_mean[ii]=0;
+ m_wave_baseline_rms[ii]=0;
m_wave_status[ii]=1; // default = 1 means below threshold
- m_calibrated_nMIP[ii]=NaN;
- m_calibrated_E_dep[ii]=NaN;
- m_calibrated_E_EM[ii]=NaN;
+ m_calibrated_nMIP[ii]=0;
+ m_calibrated_E_dep[ii]=0;
+ m_calibrated_E_EM[ii]=0;
}
+ m_scintHit = 0;
+
m_calo_total_nMIP=0;
m_calo_total_E_dep=0;
m_calo_total_E_EM=0;
+ m_calo_total_E_EM_fudged=0;
m_preshower_total_nMIP=0;
m_preshower_total_E_dep=0;
@@ -1264,6 +1301,11 @@ NtupleDumperAlg::clearTree() const
m_thetaxCalo.clear();
m_thetayCalo.clear();
+ m_x_atMaxRadius.clear();
+ m_y_atMaxRadius.clear();
+ m_z_atMaxRadius.clear();
+ m_r_atMaxRadius.clear();
+
m_t_pdg.clear();
m_t_barcode.clear();
m_t_truthHitRatio.clear();
@@ -1292,7 +1334,7 @@ NtupleDumperAlg::clearTree() const
m_truthParticleIsFiducial.clear();
m_truthLeptonMomentum = 0;
- m_truthBarcode = 0;
+ m_truthBarcode = -1;
m_truthPdg = 0;
m_truth_P.clear();
@@ -1334,28 +1376,32 @@ NtupleDumperAlg::clearTree() const
}
-void NtupleDumperAlg::setNaN() const {
+void NtupleDumperAlg::clearTrackTruth() const {
m_t_pdg.push_back(0);
m_t_barcode.push_back(-1);
- m_t_truthHitRatio.push_back(NaN);
- m_t_prodVtx_x.push_back(NaN);
- m_t_prodVtx_y.push_back(NaN);
- m_t_prodVtx_z.push_back(NaN);
- m_t_decayVtx_x.push_back(NaN);
- m_t_decayVtx_y.push_back(NaN);
- m_t_decayVtx_z.push_back(NaN);
- m_t_px.push_back(NaN);
- m_t_py.push_back(NaN);
- m_t_pz.push_back(NaN);
- m_t_theta.push_back(NaN);
- m_t_phi.push_back(NaN);
- m_t_p.push_back(NaN);
- m_t_pT.push_back(NaN);
- m_t_eta.push_back(NaN);
+ m_t_truthHitRatio.push_back(0);
+ m_t_prodVtx_x.push_back(999999);
+ m_t_prodVtx_y.push_back(999999);
+ m_t_prodVtx_z.push_back(999999);
+ m_t_decayVtx_x.push_back(999999);
+ m_t_decayVtx_y.push_back(999999);
+ m_t_decayVtx_z.push_back(999999);
+ m_t_px.push_back(0);
+ m_t_py.push_back(0);
+ m_t_pz.push_back(0);
+ m_t_theta.push_back(999999);
+ m_t_phi.push_back(999999);
+ m_t_p.push_back(0);
+ m_t_pT.push_back(0);
+ m_t_eta.push_back(999999);
for (int station = 0; station < 4; ++station) {
- m_t_st_x[station].push_back(NaN);
- m_t_st_y[station].push_back(NaN);
- m_t_st_z[station].push_back(NaN);
+ m_t_st_x[station].push_back(999999);
+ m_t_st_y[station].push_back(999999);
+ m_t_st_z[station].push_back(999999);
}
m_isFiducial.push_back(false);
}
+
+double NtupleDumperAlg::radius(const Acts::Vector3 &position) const {
+ return std::sqrt(position.x() * position.x() + position.y() * position.y());
+}
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 6d4ccc381a4b65831678a657c035ecc4e79c2b65..9c4a90240c8cbc3fc21681e6b254aecc2064b35d 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -52,12 +52,13 @@ private:
bool waveformHitOK(const xAOD::WaveformHit* hit) const;
void clearTree() const;
- void setNaN() const;
+ void clearTrackTruth() const;
void addBranch(const std::string &name,float* var);
void addBranch(const std::string &name,unsigned int* var);
void addWaveBranches(const std::string &name, int nchannels, int first);
void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
void addCalibratedBranches(const std::string &name, int nchannels, int first);
+ double radius(const Acts::Vector3 &position) const;
ServiceHandle <ITHistSvc> m_histSvc;
@@ -101,17 +102,23 @@ 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" };
- IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
- DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
- 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"};
+ 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 no veto signal adn Calo signal > 20 GeV e-" };
+ BooleanProperty m_onlyBlinded { this, "OnlyBlinded", false, "Only events that would be blinded are saved" };
+ BooleanProperty m_doTrigFilter { this, "DoTrigFilter", false, "Only events that pass trigger cuts are passed" };
+ 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" };
+
+ 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" };
+ IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
+ DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
+ 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)"};
+ DoubleProperty m_caloMC_FudgeFactor {this, "CaloFudgeFactorMC", 1.088, "Correct the MC energy calibration by this fudge factor"};
double m_baseEventCrossSection {1.0};
const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -150,7 +157,9 @@ private:
mutable float m_wave_baseline_mean[15];
mutable float m_wave_baseline_rms[15];
mutable unsigned int m_wave_status[15];
-
+
+ mutable unsigned int m_scintHit;
+
mutable float m_calibrated_nMIP[15];
mutable float m_calibrated_E_dep[15];
mutable float m_calibrated_E_EM[15];
@@ -159,6 +168,8 @@ private:
mutable float m_calo_total_E_dep;
mutable float m_calo_total_E_EM;
+ mutable float m_calo_total_E_EM_fudged;
+
mutable float m_preshower_total_nMIP;
mutable float m_preshower_total_E_dep;
@@ -170,9 +181,6 @@ private:
mutable float m_Preshower12_Edep;
mutable float m_Preshower13_Edep;
- mutable float m_MIP_sim_Edep_calo;
- mutable float m_MIP_sim_Edep_preshower;
-
mutable float m_clock_phase;
mutable unsigned int m_station0Clusters;
@@ -247,6 +255,10 @@ private:
mutable std::vector<double> m_yCalo;
mutable std::vector<double> m_thetaxCalo;
mutable std::vector<double> m_thetayCalo;
+ mutable std::vector<double> m_x_atMaxRadius;
+ mutable std::vector<double> m_y_atMaxRadius;
+ mutable std::vector<double> m_z_atMaxRadius;
+ mutable std::vector<double> m_r_atMaxRadius;
mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
@@ -324,9 +336,6 @@ private:
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;