diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx index 4d122f2e8ec6844496c4bff025779fca9045da1c..12f5ee5eed4a08a8ce0d3119f0bfe6a8430cb4d9 100644 --- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx +++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx @@ -98,11 +98,18 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const { // Loop over wavefrom vectors to make and store waveform unsigned int digitizer_samples = m_digiTool->digitizer_samples(); for (const auto& w : waveforms) { - RawWaveform* wfm = new RawWaveform(); + RawWaveform* wfm = new RawWaveform(); wfm->setWaveform(0, w.second); wfm->setIdentifier(w.first); - wfm->setChannel(m_mappingTool->getChannelMapping(w.first)); wfm->setSamples(digitizer_samples); + // Don't write out waveforms with no channel defined + // Causes problems in downstream reco (not known in cable map + int channel = m_mappingTool->getChannelMapping(w.first); + if (channel < 0) { + ATH_MSG_DEBUG("Skipping waveform with unknown channel!"); + continue; + } + wfm->setChannel(channel); waveformContainerHandle->push_back(wfm); } diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx index c8844946d2f8b1242c22ec0540a07d23ff2c7ada..5b132d2f9ce057254aedf9937f0359eea0303196 100644 --- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx +++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx @@ -51,14 +51,7 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const { } SG::ReadHandle<xAOD::WaveformHitContainer> calo2WaveHitHandle(m_calo2WaveHitContainerKey, ctx); - - ATH_CHECK( calo2WaveHitHandle.isValid() ); - ATH_MSG_DEBUG("Found ReadHandle for WaveformHitContainer " << m_calo2WaveHitContainerKey); - - if (calo2WaveHitHandle->size() == 0) { - ATH_MSG_DEBUG("Calorimeter 2 Waveform Hit container found with zero length!"); - } - + SG::ReadHandle<xAOD::WaveformHitContainer> preshowerWaveHitHandle(m_preshowerWaveHitContainerKey, ctx); ATH_CHECK( preshowerWaveHitHandle.isValid() ); @@ -87,21 +80,6 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const { } ATH_MSG_DEBUG("CaloHitContainer '" << caloHitContainerHandle.name() << "' filled with "<< caloHitContainerHandle->size() <<" items"); - // Reconstruct high-gain calorimeter hits - SG::WriteHandle<xAOD::CalorimeterHitContainer> calo2HitContainerHandle(m_calo2HitContainerKey, ctx); - ATH_CHECK( calo2HitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(), - std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) ); - ATH_MSG_DEBUG("WaveformsHitContainer '" << calo2HitContainerHandle.name() << "' initialized"); - - for( const auto& hit : *calo2WaveHitHandle ) { - if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue; - xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit(); - calo2HitContainerHandle->push_back(calo_hit); - calo_hit->addHit(calo2WaveHitHandle.get(), hit); - m_recoCalibTool->reconstruct(ctx, calo_hit, correct_gain=false); - } - ATH_MSG_DEBUG("Calo2HitContainer '" << calo2HitContainerHandle.name() << "' filled with "<< calo2HitContainerHandle->size() <<" items"); - // Reconstruct preshower hits SG::WriteHandle<xAOD::CalorimeterHitContainer> preshowerHitContainerHandle(m_preshowerHitContainerKey, ctx); ATH_CHECK( preshowerHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(), @@ -117,6 +95,35 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const { } ATH_MSG_DEBUG("PreshowerHitContainer '" << preshowerHitContainerHandle.name() << "' filled with "<< preshowerHitContainerHandle->size() <<" items"); + + // High-gain calo isn't guaranteed to exist + if ( calo2WaveHitHandle.isValid() ) { + ATH_MSG_DEBUG("Found ReadHandle for WaveformHitContainer " << m_calo2WaveHitContainerKey); + + if (calo2WaveHitHandle->size() == 0) { + ATH_MSG_DEBUG("Calorimeter 2 Waveform Hit container found with zero length!"); + } + + // Reconstruct high-gain calorimeter hits + SG::WriteHandle<xAOD::CalorimeterHitContainer> calo2HitContainerHandle(m_calo2HitContainerKey, ctx); + ATH_CHECK( calo2HitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(), + std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) ); + ATH_MSG_DEBUG("WaveformsHitContainer '" << calo2HitContainerHandle.name() << "' initialized"); + + for( const auto& hit : *calo2WaveHitHandle ) { + if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue; + xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit(); + calo2HitContainerHandle->push_back(calo_hit); + calo_hit->addHit(calo2WaveHitHandle.get(), hit); + m_recoCalibTool->reconstruct(ctx, calo_hit, correct_gain=false); + } + ATH_MSG_DEBUG("Calo2HitContainer '" << calo2HitContainerHandle.name() << "' filled with "<< calo2HitContainerHandle->size() <<" items"); + + } else { + ATH_MSG_DEBUG("No ReadHandle for WaveformHitContainer " << m_calo2WaveHitContainerKey); + } + + return StatusCode::SUCCESS; } diff --git a/Calorimeter/CaloRecTools/src/CaloRecTool.h b/Calorimeter/CaloRecTools/src/CaloRecTool.h index de736af74393bf424b2fd471795aea7bc2808f52..82ba2b60740dade0e64e91273e089d2ced98780a 100644 --- a/Calorimeter/CaloRecTools/src/CaloRecTool.h +++ b/Calorimeter/CaloRecTools/src/CaloRecTool.h @@ -91,7 +91,7 @@ class CaloRecTool: public extends<AthAlgTool, ICaloRecTool> { // Could also put these in DB, but just hardcode them for now FloatProperty m_MIP_sim_Edep_calo {this, "MIP_sim_Edep_calo", 58.5}; // MIP deposits 5.85 MeV of energy in calo - FloatProperty m_MIP_sim_Edep_calo2 {this, "MIP_sim_Edep_calo", 58.5}; // MIP deposits 5.85 MeV of energy in calo + FloatProperty m_MIP_sim_Edep_calo2 {this, "MIP_sim_Edep_calo2", 58.5}; // MIP deposits 5.85 MeV of energy in calo FloatProperty m_MIP_sim_Edep_preshower {this, "MIP_sim_Edep_preshower", 4.894}; // MIP deposits 4.894 MeV of energy in a preshower layer FloatProperty m_calo_EM_mu {this, "m_calo_EM_mu", 330.0}; // factor used to do rough calibration of calo to EM energy: 0.33 GeV or 330 MeV diff --git a/Neutrino/NeutrinoRecAlgs/src/NeutrinoRecAlgs.cxx b/Neutrino/NeutrinoRecAlgs/src/NeutrinoRecAlgs.cxx index 9aa8340b23742dca635703df70be9a9f6af45918..f6bee3da6c914ab6318a0ba9635292617bd4830e 100644 --- a/Neutrino/NeutrinoRecAlgs/src/NeutrinoRecAlgs.cxx +++ b/Neutrino/NeutrinoRecAlgs/src/NeutrinoRecAlgs.cxx @@ -4,6 +4,7 @@ #include "NeutrinoReadoutGeometry/EmulsionDetectorManager.h" #include "StoreGate/StoreGateSvc.h" #include "NeutrinoIdentifier/EmulsionID.h" +#include "NeutrinoSimEvent/NeutrinoHitIdHelper.h" #include "GeoPrimitives/CLHEPtoEigenConverter.h" NeutrinoRecAlgs::NeutrinoRecAlgs(const std::string& name, ISvcLocator* pSvcLocator) @@ -231,8 +232,17 @@ StatusCode NeutrinoRecAlgs::execute() m_trackid_begin_out_particle = 0; m_trackid_end_in_particle = 0; m_trackid_end_out_particle = 0; - m_num_in_particle = 0; - m_num_out_particle = 0; + + m_num_in_particle = -1; + m_num_out_particle = -1; + + m_vx_prod = -1; + m_vy_prod = -1; + m_vz_prod = -1; + + m_vx_decay = -1; + m_vy_decay = -1; + m_vz_decay = -1; m_pdg_in_particle.clear(); m_pdg_out_particle.clear(); @@ -321,14 +331,22 @@ StatusCode NeutrinoRecAlgs::execute() // hit.print(); Double_t thickness_base = 0.21;//mm - Double_t thickness_plate = 1.;//mm - Double_t thickness_film = 0.07;//mm + Double_t thickness_plate = 1.087;//mm + Double_t thickness_film = 0.065;//mm //module_num * base_num = 770 plates //Int_t module_num = 35; - Int_t base_num = 22; + Int_t base_num = 10; + if (m_sID == nullptr || m_sID->dictionaryVersion() == "FASERNU-03-770" || m_sID->dictionaryVersion() == "") + { + thickness_base = 0.21;//mm + thickness_plate = 1.;//mm + thickness_film = 0.07;//mm + base_num = 22; + } + //x, y : local coordinate //z : converted to the global coorfinate. m_x_start = hit.localStartPosition()[0]; diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt index 24b788d23b88e95ed37a20e86a109f198e836851..cb082e543f613fa613ef23105d5b5deebf97f481 100644 --- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt +++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt @@ -7,7 +7,7 @@ atlas_add_component( src/NtupleDumperAlg.h src/NtupleDumperAlg.cxx src/component/NtupleDumper_entries.cxx - LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib FaserActsmanVertexingLib AtlasHepMCLib + LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib FaserActsmanVertexingLib AtlasHepMCLib WaveformConditionsToolsLib PRIVATE_LINK_LIBRARIES nlohmann_json::nlohmann_json ) diff --git a/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py index da44697fc23f3318f599d4451f03fcc196a3f463..cc68970f9447e8b93433d93153b349247bcdad8a 100644 --- a/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py +++ b/PhysicsAnalysis/NtupleDumper/python/NtupleDumperConfig.py @@ -5,6 +5,7 @@ from AthenaConfiguration.ComponentFactory import CompFactory from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg +from WaveformConditionsTools.WaveformCableMappingConfig import WaveformCableMappingCfg def NtupleDumperAlgCfg(flags, OutName, **kwargs): # Initialize GeoModel @@ -14,6 +15,7 @@ def NtupleDumperAlgCfg(flags, OutName, **kwargs): acc.merge(MagneticFieldSvcCfg(flags)) #acc.merge(ActsTrackingGeometrySvcCfg(flags)) #acc.merge(FaserActsAlignmentCondAlgCfg(flags)) + acc.merge(WaveformCableMappingCfg(flags, **kwargs)) result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags) acc.merge(result) diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx index 468a6739bdd9606cc8f34c4ff0996cbda393273d..c751d105649810a1b0db80f7b5fadc8d8bf624d6 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx @@ -31,17 +31,76 @@ NtupleDumperAlg::NtupleDumperAlg(const std::string &name, void NtupleDumperAlg::addBranch(const std::string &name, - float* var) { + float* var) const { m_tree->Branch(name.c_str(),var,(name+"/F").c_str()); } void NtupleDumperAlg::addBranch(const std::string &name, - unsigned int* var) { + unsigned int* var) const { m_tree->Branch(name.c_str(),var,(name+"/I").c_str()); } +// Have to declare this cost to call during execute +void NtupleDumperAlg::defineWaveBranches() const { + + ATH_MSG_DEBUG("defineWaveBranches called"); + + // Use waveform map to find all defined waveform channels + auto mapping = m_mappingTool->getCableMapping(); + ATH_MSG_DEBUG("Cable mapping contains " << mapping.size() << " entries"); + + // mapping is std::map<int, std::pair<std::string, Identifier> > + // use this to fill my own map of channel lists keyed by type + + std::map<std::string, std::list<int>> wave_map; + + for (const auto& [key, value] : mapping) { + wave_map[value.first].push_back(key); + ATH_MSG_DEBUG("Found mapping " << value.first << " chan " << key); + } + + // Now go through found types and define ntuple entries + // Keys are defined in cable map and used by RawWaveformDecoder + for (const auto& [key, value] : wave_map) { + + if (key == std::string("calo")) { + addWaveBranches("CaloLo", value); + } + else if (key == std::string("calo2")) { + addWaveBranches("CaloHi", value); + } + else if (key == std::string("calonu")) { + addWaveBranches("CaloNu", value); + } + else if (key == std::string("veto")) { + addWaveBranches("Veto", value); + } + else if (key == std::string("vetonu")) { + addWaveBranches("VetoNu", value); + } + else if (key == std::string("trigger")) { + addWaveBranches("Timing", value); + } + else if (key == std::string("preshower")) { + addWaveBranches("Preshower", value); + } + else if (key == std::string("clock")) { + // ignore this + } + else if (key == std::string("none")) { + // ignore this + } + else { + ATH_MSG_WARNING("Found unknown mapping type " << key); + } + } + + ATH_MSG_DEBUG("defineWaveBranches done"); + +} + void NtupleDumperAlg::addWaveBranches(const std::string &name, int nchannels, - int first) { + int first) const { for(int ch=0;ch<nchannels;ch++) { std::string base=name+std::to_string(ch)+"_"; addBranch(base+"time",&m_wave_localtime[first]); @@ -57,22 +116,64 @@ void NtupleDumperAlg::addWaveBranches(const std::string &name, } } -void NtupleDumperAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave) const { +// Use channel list to add branches +void NtupleDumperAlg::addWaveBranches(const std::string &name, + std::list<int> channel_list) const { + + ATH_MSG_DEBUG("Adding " << name << " with channels " << channel_list); + + int nchannels = channel_list.size(); + for(int i=0;i<nchannels;i++) { + int ch = channel_list.front(); + channel_list.pop_front(); + + std::string base=name+std::to_string(i); + addBranch(base+"_triggertime",&m_wave_triggertime[ch]); + addBranch(base+"_localtime",&m_wave_localtime[ch]); + addBranch(base+"_bcidtime",&m_wave_bcidtime[ch]); + addBranch(base+"_peak",&m_wave_peak[ch]); + addBranch(base+"_width",&m_wave_width[ch]); + addBranch(base+"_charge",&m_wave_charge[ch]); + addBranch(base+"_raw_peak",&m_wave_raw_peak[ch]); + addBranch(base+"_raw_charge",&m_wave_raw_charge[ch]); + addBranch(base+"_baseline",&m_wave_baseline_mean[ch]); + addBranch(base+"_baseline_rms",&m_wave_baseline_rms[ch]); + addBranch(base+"_status",&m_wave_status[ch]); + ATH_MSG_DEBUG("Added " << base << " ch " << ch); + + // Also add the random histogram + std::string hname = "hRandomCharge"+std::to_string(ch); + base = name + std::to_string(i) + " Ch" + std::to_string(ch); + std::string title = base+" Charge from RandomEvents;charge (pC);Events/bin"; + m_HistRandomCharge[ch] = new TH1F(hname.c_str(), title.c_str(), 100, -1.0, 1.0); + StatusCode sc = histSvc()->regHist("/HIST2/RandomCharge"+std::to_string(ch), m_HistRandomCharge[ch]); + } +} + +void NtupleDumperAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave, bool isMC) const { + ATH_MSG_DEBUG("FillWaveBranches called"); for (auto hit : wave) { if ((hit->hit_status()&2)==0) { // dont store secoondary hits as they can overwrite the primary hit int ch=hit->channel(); - m_wave_localtime[ch]=hit->localtime()+m_clock_phase; + ATH_MSG_DEBUG("FillWaveBranches filling channel " << ch); + + m_wave_localtime[ch]=hit->localtime(); + m_wave_bcidtime[ch] = hit->bcid_time(); // Now corrected for offsets + // Not sure why this doesn't exist in MC... + if (!isMC) + m_wave_triggertime[ch] = hit->trigger_time(); // + m_clock_phase; m_wave_peak[ch]=hit->peak(); m_wave_width[ch]=hit->width(); - m_wave_charge[ch]=hit->integral()/50; + m_wave_charge[ch]=hit->integral()/50; // In pC m_wave_raw_peak[ch]=hit->raw_peak(); - m_wave_raw_charge[ch]=hit->raw_integral()/50; + m_wave_raw_charge[ch]=hit->raw_integral()/50; // In pC m_wave_baseline_mean[ch]=hit->baseline_mean(); m_wave_baseline_rms[ch]=hit->baseline_rms(); m_wave_status[ch]=hit->hit_status(); } } + ATH_MSG_DEBUG("FillWaveBranches done"); } void NtupleDumperAlg::addCalibratedBranches(const std::string &name, @@ -89,6 +190,9 @@ void NtupleDumperAlg::addCalibratedBranches(const std::string &name, StatusCode NtupleDumperAlg::initialize() { + + m_first = true; + ATH_CHECK(m_truthEventContainer.initialize()); ATH_CHECK(m_mcEventContainer.initialize()); ATH_CHECK(m_truthParticleContainer.initialize()); @@ -102,6 +206,7 @@ StatusCode NtupleDumperAlg::initialize() ATH_CHECK(m_preshowerContainer.initialize()); ATH_CHECK(m_ecalContainer.initialize()); ATH_CHECK(m_ecal2Container.initialize()); + ATH_CHECK(m_caloNuContainer.initialize()); ATH_CHECK(m_clusterContainer.initialize()); ATH_CHECK(m_simDataCollection.initialize()); ATH_CHECK(m_FaserTriggerData.initialize()); @@ -111,14 +216,10 @@ StatusCode NtupleDumperAlg::initialize() ATH_CHECK(m_preshowerCalibratedContainer.initialize()); ATH_CHECK(m_ecalCalibratedContainer.initialize()); ATH_CHECK(m_ecal2CalibratedContainer.initialize()); + ATH_CHECK(m_calonuCalibratedContainer.initialize()); ATH_CHECK(m_eventInfoKey.initialize()); ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID")); - ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID")); - ATH_CHECK(detStore()->retrieve(m_vetoHelper, "VetoID")); - ATH_CHECK(detStore()->retrieve(m_triggerHelper, "TriggerID")); - ATH_CHECK(detStore()->retrieve(m_preshowerHelper, "PreshowerID")); - ATH_CHECK(detStore()->retrieve(m_ecalHelper, "EcalID")); ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT")); ATH_CHECK(m_extrapolationTool.retrieve()); @@ -126,7 +227,8 @@ StatusCode NtupleDumperAlg::initialize() ATH_CHECK(m_trackTruthMatchingTool.retrieve()); ATH_CHECK(m_fiducialParticleTool.retrieve()); ATH_CHECK(m_vertexingTool.retrieve()); - + ATH_CHECK(m_mappingTool.retrieve()); + ATH_CHECK(m_spacePointContainerKey.initialize()); // Read GRL if requested @@ -164,7 +266,8 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("eventTimeNSOffset", &m_event_timeNSOffset, "eventTimeNSOffset/I"); m_tree->Branch("BCID", &m_bcid, "BCID/I"); m_tree->Branch("inGRL", &m_in_grl, "inGRL/I"); - + m_tree->Branch("clockPhase", &m_clock_phase, "clockPhase"); + m_tree->Branch("fillNumber", &m_fillNumber, "fillNumber/I"); m_tree->Branch("betaStar", &m_betaStar, "betaStar/F"); m_tree->Branch("crossingAngle", &m_crossingAngle, "crossingAngle/F"); @@ -180,15 +283,8 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("inputBits", &m_inputBits, "inputBits/I"); m_tree->Branch("inputBitsNext", &m_inputBitsNext, "inputBitsNext/I"); - //WAVEFORMS - // Need to put an option in here for pre/post 2024 data - addWaveBranches("VetoNu",2,4); - addWaveBranches("VetoSt1",2,14); - addWaveBranches("VetoSt2",2,6); - addWaveBranches("Timing",4,8); - addWaveBranches("Preshower",2,12); - addWaveBranches("Calo",4,0); - addWaveBranches("CaloHi", 4,16); + // WAVEFORMS + // Now defined on first event so we can use cable mapping tool m_tree->Branch("ScintHit", &m_scintHit); @@ -375,6 +471,33 @@ StatusCode NtupleDumperAlg::initialize() 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); + m_tree->Branch("truth_isFiducial", &m_truth_isFiducial); + + + m_tree->Branch("truth_st0_x", &m_truth_st_x[0]); + m_tree->Branch("truth_st0_y", &m_truth_st_y[0]); + m_tree->Branch("truth_st0_z", &m_truth_st_z[0]); + m_tree->Branch("truth_st1_x", &m_truth_st_x[1]); + m_tree->Branch("truth_st1_y", &m_truth_st_y[1]); + m_tree->Branch("truth_st1_z", &m_truth_st_z[1]); + m_tree->Branch("truth_st2_x", &m_truth_st_x[2]); + m_tree->Branch("truth_st2_y", &m_truth_st_y[2]); + m_tree->Branch("truth_st2_z", &m_truth_st_z[2]); + m_tree->Branch("truth_st3_x", &m_truth_st_x[3]); + m_tree->Branch("truth_st3_y", &m_truth_st_y[3]); + m_tree->Branch("truth_st3_z", &m_truth_st_z[3]); + m_tree->Branch("truth_st0_px", &m_truth_st_px[0]); + m_tree->Branch("truth_st0_py", &m_truth_st_py[0]); + m_tree->Branch("truth_st0_pz", &m_truth_st_pz[0]); + m_tree->Branch("truth_st1_px", &m_truth_st_px[1]); + m_tree->Branch("truth_st1_py", &m_truth_st_py[1]); + m_tree->Branch("truth_st1_pz", &m_truth_st_pz[1]); + m_tree->Branch("truth_st2_px", &m_truth_st_px[2]); + m_tree->Branch("truth_st2_py", &m_truth_st_py[2]); + m_tree->Branch("truth_st2_pz", &m_truth_st_pz[2]); + m_tree->Branch("truth_st3_px", &m_truth_st_px[3]); + m_tree->Branch("truth_st3_py", &m_truth_st_py[3]); + m_tree->Branch("truth_st3_pz", &m_truth_st_pz[3]); // for mother + daughter particle truth infomation @@ -393,6 +516,7 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("truthd0_x", &m_truthd0_x); m_tree->Branch("truthd0_y", &m_truthd0_y); m_tree->Branch("truthd0_z", &m_truthd0_z); + m_tree->Branch("truthd0_isFiducial", &m_truthd0_isFiducial); m_tree->Branch("truthd1_P", &m_truthd1_P); m_tree->Branch("truthd1_px", &m_truthd1_px); @@ -401,6 +525,7 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("truthd1_x", &m_truthd1_x); m_tree->Branch("truthd1_y", &m_truthd1_y); m_tree->Branch("truthd1_z", &m_truthd1_z); + m_tree->Branch("truthd1_isFiducial", &m_truthd1_isFiducial); m_tree->Branch("vertex_x", &m_vertex_x, "vertex_x/D"); m_tree->Branch("vertex_y", &m_vertex_y, "vertex_y/D"); @@ -418,48 +543,7 @@ StatusCode NtupleDumperAlg::initialize() ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree)); // Register histograms - m_HistRandomCharge[0] = new TH1F("hRandomCharge0", "Calo ch0 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[1] = new TH1F("hRandomCharge1", "Calo ch1 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[2] = new TH1F("hRandomCharge2", "Calo ch2 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[3] = new TH1F("hRandomCharge3", "Calo ch3 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[4] = new TH1F("hRandomCharge4", "VetoNu ch4 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[5] = new TH1F("hRandomCharge5", "VetoNu ch5 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[6] = new TH1F("hRandomCharge6", "Veto ch6 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[7] = new TH1F("hRandomCharge7", "Veto ch7 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[8] = new TH1F("hRandomCharge8", "Trig ch8 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[9] = new TH1F("hRandomCharge9", "Trig ch9 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[10] = new TH1F("hRandomCharge10", "Trig ch10 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[11] = new TH1F("hRandomCharge11", "Trig ch11 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[12] = new TH1F("hRandomCharge12", "Preshower ch12 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[13] = new TH1F("hRandomCharge13", "Preshower ch13 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[14] = new TH1F("hRandomCharge14", "Veto ch14 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[15] = new TH1F("hRandomCharge15", "Veto ch15 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - - m_HistRandomCharge[16] = new TH1F("hRandomCharge16", "CaloHi ch0 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[17] = new TH1F("hRandomCharge17", "CaloHi ch1 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[18] = new TH1F("hRandomCharge18", "CaloHi ch2 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - m_HistRandomCharge[19] = new TH1F("hRandomCharge19", "CaloHi ch3 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0); - - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge0", m_HistRandomCharge[0])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge1", m_HistRandomCharge[1])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge2", m_HistRandomCharge[2])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge3", m_HistRandomCharge[3])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge4", m_HistRandomCharge[4])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge5", m_HistRandomCharge[5])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge6", m_HistRandomCharge[6])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge7", m_HistRandomCharge[7])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge8", m_HistRandomCharge[8])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge9", m_HistRandomCharge[9])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge10", m_HistRandomCharge[10])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge11", m_HistRandomCharge[11])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge12", m_HistRandomCharge[12])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge15", m_HistRandomCharge[15])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge16", m_HistRandomCharge[16])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge17", m_HistRandomCharge[17])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge18", m_HistRandomCharge[18])); - ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge19", m_HistRandomCharge[19])); + // Also done on first event if (m_onlyBlinded){ ATH_MSG_INFO("Only events that would be blinded are saved in ntuple"); @@ -478,6 +562,12 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const { clearTree(); + // Define waveform channels on first event + if (m_first) { + defineWaveBranches(); + m_first = false; + } + // check if real data or simulation data bool isMC = false; SG::ReadHandle<xAOD::TruthEventContainer> truthEventContainer { m_truthEventContainer, ctx }; @@ -575,44 +665,60 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const } // process all waveform data for all scintillator and calorimeter channels + // Not all are guaranteed to exist, so fill the ones that are valid SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx }; - ATH_CHECK(vetoNuContainer.isValid()); + if (vetoNuContainer.isValid()) { + FillWaveBranches(*vetoNuContainer, isMC); + } SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx }; - ATH_CHECK(vetoContainer.isValid()); + if (vetoContainer.isValid()) { + FillWaveBranches(*vetoContainer, isMC); + } SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx }; - ATH_CHECK(triggerContainer.isValid()); + if (triggerContainer.isValid()) { + FillWaveBranches(*triggerContainer, isMC); + } SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx }; - ATH_CHECK(preshowerContainer.isValid()); + if (preshowerContainer.isValid()) { + FillWaveBranches(*preshowerContainer, isMC); + } SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx }; - ATH_CHECK(ecalContainer.isValid()); + if (ecalContainer.isValid()) { + FillWaveBranches(*ecalContainer, isMC); + } SG::ReadHandle<xAOD::WaveformHitContainer> ecal2Container { m_ecal2Container, ctx }; - ATH_CHECK(ecal2Container.isValid()); + if (ecal2Container.isValid()) { + FillWaveBranches(*ecal2Container, isMC); + } - FillWaveBranches(*vetoNuContainer); - FillWaveBranches(*vetoContainer); - FillWaveBranches(*triggerContainer); - FillWaveBranches(*preshowerContainer); - FillWaveBranches(*ecalContainer); - FillWaveBranches(*ecal2Container); + SG::ReadHandle<xAOD::WaveformHitContainer> caloNuContainer { m_caloNuContainer, ctx }; + if (caloNuContainer.isValid()) { + FillWaveBranches(*caloNuContainer, isMC); + } - // if real data, store charge in histograms from random events + // Some things for real data if (!isMC) { + // if real data, store charge in histograms from random events SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx); m_tap=triggerData->tap(); // for random (only) triggers, store charge of scintillators in histograms - if (m_tap == 16) { + if (m_tap == 16) { + ATH_MSG_DEBUG("Filling random charge histogram"); // Fill histograms - for (unsigned int chan = 0; chan<nchan; chan++) { - m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]); + for (unsigned int chan = 0; chan<max_chan; chan++) { + // Only fill histograms that have been defined + if (m_HistRandomCharge[chan]) + m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]); } } + // Apply ntuple selection for real data if (m_doTrigFilter) { bool trig_coincidence_preshower_and_vetoes = ( (m_tap&8) != 0 ); @@ -759,6 +865,21 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const m_truth_pz.push_back(particle->p4().Z()); m_truth_m.push_back(particle->m()); m_truth_pdg.push_back(particle->pdgId()); + m_truth_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode())); + + auto positions = m_fiducialParticleTool->getTruthPositions(particle->barcode()); // Add truth position information at each station + for (int station = 0; station < 4; ++station) { + m_truth_st_x[station].push_back(positions[station].x()); + m_truth_st_y[station].push_back(positions[station].y()); + m_truth_st_z[station].push_back(positions[station].z()); + } + + auto momenta = m_fiducialParticleTool->getTruthMomenta(particle->barcode()); // Add truth momentum information at each station + for (int station = 0; station < 4; ++station) { + m_truth_st_px[station].push_back(momenta[station].x()); + m_truth_st_py[station].push_back(momenta[station].y()); + m_truth_st_pz[station].push_back(momenta[station].z()); + } if ( particle->hasProdVtx()) { m_truth_prod_x.push_back(particle->prodVtx()->x()); @@ -788,6 +909,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const 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()); @@ -805,6 +927,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const m_truthd0_px.push_back(particle->p4().X()); m_truthd0_py.push_back(particle->p4().Y()); m_truthd0_pz.push_back(particle->p4().Z()); + m_truthd0_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode())); if ( particle->hasProdVtx()) { m_truthd0_x.push_back(particle->prodVtx()->x()); @@ -826,6 +949,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const m_truthd1_px.push_back(particle->p4().X()); m_truthd1_py.push_back(particle->p4().Y()); m_truthd1_pz.push_back(particle->p4().Z()); + m_truthd1_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode())); if ( particle->hasProdVtx()) { m_truthd1_x.push_back(particle->prodVtx()->x()); @@ -922,7 +1046,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const return StatusCode::SUCCESS; } - SG::ReadDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx); + SG::ReadHandle<xAOD::EventInfo> 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"); @@ -1404,7 +1528,8 @@ NtupleDumperAlg::clearTree() const m_event_timeNSOffset = 0; m_bcid = 0; m_in_grl = 0; - + m_clock_phase = 0; + m_fillNumber = 0; m_betaStar = 0; m_crossingAngle = 0; @@ -1420,8 +1545,10 @@ NtupleDumperAlg::clearTree() const m_inputBits=0; m_inputBitsNext=0; - for(unsigned int ii=0;ii<nchan;ii++) { + for(unsigned int ii=0;ii<max_chan;ii++) { m_wave_localtime[ii]=0; + m_wave_triggertime[ii]=0; + m_wave_bcidtime[ii]=0; m_wave_peak[ii]=0; m_wave_width[ii]=0; m_wave_charge[ii]=0; @@ -1595,6 +1722,18 @@ NtupleDumperAlg::clearTree() const m_truth_dec_x.clear(); m_truth_dec_y.clear(); m_truth_dec_z.clear(); + m_truth_isFiducial.clear(); + + for (int station = 0; station < 4; ++station) { + m_truth_st_x[station].clear(); + m_truth_st_y[station].clear(); + m_truth_st_z[station].clear(); + } + for (int station = 0; station < 4; ++station) { + m_truth_st_px[station].clear(); + m_truth_st_py[station].clear(); + m_truth_st_pz[station].clear(); + } m_truthM_P.clear(); m_truthM_px.clear(); @@ -1611,6 +1750,7 @@ NtupleDumperAlg::clearTree() const m_truthd0_x.clear(); m_truthd0_y.clear(); m_truthd0_z.clear(); + m_truthd0_isFiducial.clear(); m_truthd1_P.clear(); m_truthd1_px.clear(); @@ -1619,6 +1759,7 @@ NtupleDumperAlg::clearTree() const m_truthd1_x.clear(); m_truthd1_y.clear(); m_truthd1_z.clear(); + m_truthd1_isFiducial.clear(); m_vertex_x = NaN; m_vertex_y = NaN; diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h index d4fc22b9754617ea8c68aba3c365ef334abd0d5d..36729773eec08aa8674dfebb43bec17953731c92 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h @@ -22,11 +22,12 @@ #include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h" #include "TrackerSimEvent/FaserSiHitCollection.h" #include "xAODEventInfo/EventInfo.h" -#include "StoreGate/ReadDecorHandle.h" #include "FaserActsVertexing/IVertexingTool.h" #include "GeneratorObjects/McEventCollection.h" #include <boost/dynamic_bitset.hpp> +#include "WaveformConditionsTools/IWaveformCableMappingTool.h" + #include <vector> #include <nlohmann/json.hpp> @@ -56,13 +57,17 @@ public: private: + mutable bool m_first; + bool waveformHitOK(const xAOD::WaveformHit* hit) const; void clearTree() 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 addBranch(const std::string &name,float* var) const; + void addBranch(const std::string &name,unsigned int* var) const; + void defineWaveBranches() const; + void addWaveBranches(const std::string &name, int nchannels, int first) const; + void addWaveBranches(const std::string &name, std::list<int> channel_list) const; + void FillWaveBranches(const xAOD::WaveformHitContainer &wave, bool isMC) const; void addCalibratedBranches(const std::string &name, int nchannels, int first); double radius(const Acts::Vector3 &position) const; @@ -85,9 +90,11 @@ private: SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerContainer { this, "PreshowerContainer", "PreshowerWaveformHits", "Preshower hit container name" }; SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecalContainer { this, "EcalContainer", "CaloWaveformHits", "Ecal hit container name" }; SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecal2Container { this, "Ecal2Container", "Calo2WaveformHits", "Ecal hit container name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_caloNuContainer { this, "CaloNuContainer", "CaloNuWaveformHits", "CaloNu hit container name" }; SG::ReadHandleKey<xAOD::CalorimeterHitContainer> m_ecalCalibratedContainer { this, "EcalCalibratedContainer", "CaloHits", "Ecal Calibrated hit container name" }; SG::ReadHandleKey<xAOD::CalorimeterHitContainer> m_ecal2CalibratedContainer { this, "Ecal2CalibratedContainer", "Calo2Hits", "Ecal Calibrated hit container name" }; + SG::ReadHandleKey<xAOD::CalorimeterHitContainer> m_calonuCalibratedContainer { this, "CaloNuCalibratedContainer", "CaloNuHits", "CaloNu Calibrated hit container name" }; SG::ReadHandleKey<xAOD::CalorimeterHitContainer> m_preshowerCalibratedContainer { this, "preshowerCalibratedContainer", "PreshowerHits", "Preshower Calibrated hit container name" }; SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer", "Tracker cluster container name" }; @@ -96,21 +103,17 @@ 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"}; - SG::ReadDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"}; + SG::ReadHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo", "ReadHandleKey for 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"}; ToolHandle<FaserTracking::IVertexingTool> m_vertexingTool { this, "VertexingTool", "FaserTracking::PointOfClosestApproachSearchTool"}; + ToolHandle<IWaveformCableMappingTool> m_mappingTool { this, "WaveformCableMappingTool", "WaveformCableMappingTool" }; const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr}; const FaserSCT_ID* m_sctHelper; - const VetoNuID* m_vetoNuHelper; - const VetoID* m_vetoHelper; - const TriggerID* m_triggerHelper; - 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" }; @@ -148,9 +151,9 @@ private: mutable TTree* m_tree; - const static unsigned int nchan=20; - - mutable TH1* m_HistRandomCharge[nchan]; + //mutable unsigned int n_wave_chan; // Actual number of waveform channels + const static unsigned int max_chan=32; + mutable TH1* m_HistRandomCharge[max_chan]; mutable unsigned int m_run_number; mutable unsigned int m_event_number; @@ -174,22 +177,24 @@ private: mutable unsigned int m_inputBits; mutable unsigned int m_inputBitsNext; - mutable float m_wave_localtime[nchan]; - mutable float m_wave_peak[nchan]; - mutable float m_wave_width[nchan]; - mutable float m_wave_charge[nchan]; - - mutable float m_wave_raw_peak[nchan]; - mutable float m_wave_raw_charge[nchan]; - mutable float m_wave_baseline_mean[nchan]; - mutable float m_wave_baseline_rms[nchan]; - mutable unsigned int m_wave_status[nchan]; + mutable float m_wave_localtime[max_chan]; + mutable float m_wave_triggertime[max_chan]; + mutable float m_wave_bcidtime[max_chan]; + mutable float m_wave_peak[max_chan]; + mutable float m_wave_width[max_chan]; + mutable float m_wave_charge[max_chan]; + + mutable float m_wave_raw_peak[max_chan]; + mutable float m_wave_raw_charge[max_chan]; + mutable float m_wave_baseline_mean[max_chan]; + mutable float m_wave_baseline_rms[max_chan]; + mutable unsigned int m_wave_status[max_chan]; mutable unsigned int m_scintHit; - mutable float m_calibrated_nMIP[nchan]; - mutable float m_calibrated_E_dep[nchan]; - mutable float m_calibrated_E_EM[nchan]; + mutable float m_calibrated_nMIP[max_chan]; + mutable float m_calibrated_E_dep[max_chan]; + mutable float m_calibrated_E_EM[max_chan]; mutable float m_calo_total_nMIP; mutable float m_calo_total_E_dep; @@ -340,6 +345,7 @@ private: mutable std::vector<double> m_truthM_z; mutable std::vector<double> m_truthd0_P; + mutable std::vector<bool> m_truthd0_isFiducial; //Boolean for if electron is fiducial mutable std::vector<double> m_truthd0_px; mutable std::vector<double> m_truthd0_py; @@ -350,6 +356,7 @@ private: mutable std::vector<double> m_truthd0_z; mutable std::vector<double> m_truthd1_P; + mutable std::vector<bool> m_truthd1_isFiducial; // Boolean for if positron is fiducial mutable std::vector<double> m_truthd1_px; mutable std::vector<double> m_truthd1_py; @@ -376,6 +383,14 @@ private: mutable std::vector<double> m_truth_prod_z; mutable std::vector<int> m_truth_pdg; // pdg of first 10 truth particles + mutable std::vector<bool> m_truth_isFiducial; // Boolean for if first 10 truth particles are fiducial + + mutable std::array<std::vector<double>, 4> m_truth_st_x; // vector of the x components of the simulated hits of the truth particle for each station + mutable std::array<std::vector<double>, 4> m_truth_st_y; // vector of the y components of the simulated hits of the truth particle for each station + mutable std::array<std::vector<double>, 4> m_truth_st_z; // vector of the z components of the simulated hits of the truth particle for each station + mutable std::array<std::vector<double>, 4> m_truth_st_px; // x components of the true momentum at each station + mutable std::array<std::vector<double>, 4> m_truth_st_py; // y components of the true momentum at each station + mutable std::array<std::vector<double>, 4> m_truth_st_pz; // z components of the true momentum at each station mutable double m_truthLeptonMomentum; mutable int m_truthBarcode; diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py index 6af8ed0e1c78771a1121bce792b4d69105dc46e1..8600247ed919a34df3f81dfb6d059eb8fba28b87 100644 --- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py +++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py @@ -16,7 +16,9 @@ def WaveformReconstructionCfg(flags): """ Return all algorithms and tools for Waveform reconstruction """ acc = ComponentAccumulator() - if not flags.Input.isMC: + if flags.Input.isMC: + print("ClockRecAlg diabled for MC in WaveformReconstructionCfg!") + else: acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg")) if "TB" in flags.GeoModel.FaserVersion: @@ -28,7 +30,10 @@ def WaveformReconstructionCfg(flags): else: acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo")) - acc.merge(WaveformHitRecCfg(flags, "Calo2WaveformRecAlg", "Calo2")) + if flags.Input.isMC: + print("Calo2WaveformRecAlg diabled for MC in WaveformReconstructionCfg!") + else: + acc.merge(WaveformHitRecCfg(flags, "Calo2WaveformRecAlg", "Calo2")) # Make preshower window 200 ns wide (value in digitizer clock ticks) acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto", FitWindowWidth=100 ))