diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx index 9bb264022f15eac85bdc2576b48a53c5243ef96a..796491358924bd5e22978f1a923fb3f6d45d851d 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx @@ -79,25 +79,37 @@ void NtupleDumperAlg::defineWaveBranches() const { for (const auto& [key, value] : wave_map) { if (key == std::string("calo")) { - addWaveBranches("CaloLo", value); + m_calolo_channels = value; + // Figure out ntuple below } else if (key == std::string("calo2")) { - addWaveBranches("CaloHi", value); + m_calohi_channels = value; + // Figure out ntuple below } else if (key == std::string("calonu")) { + m_calonu_channels = value; addWaveBranches("CaloNu", value); + addCalibratedBranches("CaloNu", value); + addCaloSumBranches("CaloNu"); } else if (key == std::string("veto")) { + m_veto_channels = value; addWaveBranches("Veto", value); } else if (key == std::string("vetonu")) { + m_vetonu_channels = value; addWaveBranches("VetoNu", value); } else if (key == std::string("trigger")) { + m_trigger_channels = value; addWaveBranches("Timing", value); } else if (key == std::string("preshower")) { + m_preshower_channels = value; addWaveBranches("Preshower", value); + addCalibratedBranches("Preshower", value); + addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP); + addBranch("Preshower_total_E_dep", &m_preshower_total_E_dep); } else if (key == std::string("clock")) { // ignore this @@ -108,10 +120,79 @@ void NtupleDumperAlg::defineWaveBranches() const { else { ATH_MSG_WARNING("Found unknown mapping type " << key); } + } // End of loop over wave_map entries + + // Figure out calorimeter + if (m_calolo_channels.size() && !m_calohi_channels.size()) { + // Only "calo" type found, 2022/23 configuration + ATH_MSG_INFO("defineWaveBranches found " << m_calolo_channels.size() << " calo channels => 2022/23 configuration"); + m_dual_gain = false; + m_calo_channels = m_calolo_channels; + m_calolo_channels.clear(); + addWaveBranches("Calo", m_calo_channels); + addCalibratedBranches("Calo", m_calo_channels); + addCaloSumBranches("Calo"); } - + else if (m_calolo_channels.size() && m_calohi_channels.size()) { + // Both found, 2024/25 configuration (both defined) + ATH_MSG_INFO("defineWaveBranches found " << m_calolo_channels.size() << " calo and " << m_calohi_channels.size() << " calo2 channels => 2024/25 configuration"); + m_dual_gain = true; + addWaveBranches("CaloLo", m_calolo_channels); + addWaveBranches("CaloHo", m_calohi_channels); + addCalibratedBranches("CaloLo", m_calolo_channels); + addCalibratedBranches("CaloHo", m_calohi_channels); + addCaloSumBranches("CaloLo"); + addCaloSumBranches("CaloHi"); + } + else { + // Unknown configuration + ATH_MSG_WARNING("defineWaveBranches found " << m_calolo_channels.size() << " calo channels and " << m_calohi_channels.size() << " calo2 channels!"); + m_dual_gain = false; + return; + } + ATH_MSG_DEBUG("defineWaveBranches done"); +} +void +NtupleDumperAlg::addCaloSumBranches(const std::string &name) const { + + // Utility to define the calorimeter sum branches + if (name == std::string("Calo")) { + addBranch("Calo_total_nMIP", &m_calo_total_nMIP); + addBranch("Calo_total_E_dep", &m_calo_total_E_dep); + addBranch("Calo_total_fit_E_EM", &m_calo_total_fit_E_EM); + addBranch("Calo_total_raw_E_EM", &m_calo_total_raw_E_EM); + addBranch("Calo_total_E_EM", &m_calo_total_E_EM); + + addBranch("Calo_total_fit_E_EM_fudged", &m_calo_total_fit_E_EM_fudged); + addBranch("Calo_total_raw_E_EM_fudged", &m_calo_total_raw_E_EM_fudged); + addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged); + } + else if (name == std::string("CaloLo")) { + addBranch("CaloLo_total_nMIP", &m_calolo_total_nMIP); + addBranch("CaloLo_total_E_dep", &m_calolo_total_E_dep); + addBranch("CaloLo_total_fit_E_EM", &m_calolo_total_fit_E_EM); + addBranch("CaloLo_total_raw_E_EM", &m_calolo_total_raw_E_EM); + addBranch("CaloLo_total_E_EM", &m_calolo_total_E_EM); + } + else if (name == std::string("CaloHi")) { + addBranch("CaloHi_total_nMIP", &m_calohi_total_nMIP); + addBranch("CaloHi_total_E_dep", &m_calohi_total_E_dep); + addBranch("CaloHi_total_fit_E_EM", &m_calohi_total_fit_E_EM); + addBranch("CaloHi_total_raw_E_EM", &m_calohi_total_raw_E_EM); + addBranch("CaloHi_total_E_EM", &m_calohi_total_E_EM); + } + else if (name == std::string("CaloNu")) { + addBranch("CaloNu_total_nMIP", &m_calonu_total_nMIP); + addBranch("CaloNu_total_E_dep", &m_calonu_total_E_dep); + addBranch("CaloNu_total_fit_E_EM", &m_calonu_total_fit_E_EM); + addBranch("CaloNu_total_raw_E_EM", &m_calonu_total_raw_E_EM); + addBranch("CaloNu_total_E_EM", &m_calonu_total_E_EM); + } + else { + ATH_MSG_WARNING("addCaloSumBranches - Unknown type " << name << "!"); + } } // Use channel list to add branches @@ -175,17 +256,25 @@ void NtupleDumperAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave, b } void NtupleDumperAlg::addCalibratedBranches(const std::string &name, - int nchannels, - int first) { - for(int ch=0;ch<nchannels;ch++) { - std::string base=name+std::to_string(ch)+"_"; - addBranch(base+"nMIP",&m_calibrated_nMIP[first]); - addBranch(base+"E_dep",&m_calibrated_E_dep[first]); - addBranch(base+"E_EM",&m_calibrated_E_EM[first]); - first++; + 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+"_nMIP", &m_calibrated_nMIP[ch]); + addBranch(base+"_E_dep",&m_calibrated_E_dep[ch]); + addBranch(base+"_E_EM", &m_calibrated_E_EM[ch]); + ATH_MSG_DEBUG("Added " << base << " ch " << ch); + } } + StatusCode NtupleDumperAlg::initialize() { @@ -289,19 +378,6 @@ StatusCode NtupleDumperAlg::initialize() 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_fit_E_EM", &m_calo_total_fit_E_EM); - addBranch("Calo_total_raw_E_EM", &m_calo_total_raw_E_EM); - addBranch("Calo_total_E_EM", &m_calo_total_E_EM); - addBranch("Calo_total_fit_E_EM_fudged", &m_calo_total_fit_E_EM_fudged); - addBranch("Calo_total_raw_E_EM_fudged", &m_calo_total_raw_E_EM_fudged); - addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged); - - addCalibratedBranches("Preshower",2,12); - addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP); - addBranch("Preshower_total_E_dep", &m_preshower_total_E_dep); //TRACKER addBranch("nClusters0",&m_station0Clusters); @@ -324,7 +400,7 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("TrackSegment_py", &m_trackseg_py); m_tree->Branch("TrackSegment_pz", &m_trackseg_pz); - addBranch("TrueTrackSegments",&m_nTruetracksegs); + addBranch("TrueTrackSegments",&m_nTruetracksegs); m_tree->Branch("TrueTrackSegment_Chi2", &m_Truetrackseg_Chi2); m_tree->Branch("TrueTrackSegment_nDoF", &m_Truetrackseg_DoF); m_tree->Branch("TrueTrackSegment_x", &m_Truetrackseg_x); @@ -746,26 +822,24 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const return StatusCode::SUCCESS; } } + if (m_doRandomFilter) { - bool trig_random = m_tap & 16; - if ( !(trig_random)) { + bool trig_random = m_tap & 16; + if ( !(trig_random)) { // don't process events that fail to activate random trigger ATH_MSG_DEBUG("event did not pass random trigger filter"); return StatusCode::SUCCESS; } + } - } - - if (m_doRandomOnlyFilter) { - bool trig_random = m_tap == 16; - if ( !(trig_random)) { + if (m_doRandomOnlyFilter) { + bool trig_random = m_tap == 16; + if ( !(trig_random)) { // don't process events that fail to activate random only trigger ATH_MSG_DEBUG("event did not pass random trigger filter"); return StatusCode::SUCCESS; - } - - } - + } + } if (m_doScintFilter) { // Get channel mapping @@ -904,30 +978,33 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const } // 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; - } + + // This either needs to be updated to use cable map, or skipped + // Skip for now + // 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; + // } if (isMC) { // if simulation find MC cross section and primary lepton // Work out effective cross section for MC @@ -1076,18 +1153,48 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const m_calibrated_E_dep[ch] = hit->E_dep(); m_calibrated_E_EM[ch] = hit->E_EM(); - m_calo_total_nMIP += hit->Nmip(); - m_calo_total_E_dep += hit->E_dep(); - m_calo_total_fit_E_EM += hit->E_EM(); - m_calo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + ATH_MSG_DEBUG("Calibrated calo: 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()); + + if (m_dual_gain) { + m_calolo_total_nMIP += hit->Nmip(); + m_calolo_total_E_dep += hit->E_dep(); + m_calolo_total_fit_E_EM += hit->E_EM(); + m_calolo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + + if (m_wave_status[ch]&4) { // Overflows + m_calolo_total_E_EM += hit->E_EM(); + } else { + m_calolo_total_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + } + } else { + m_calo_total_nMIP += hit->Nmip(); + m_calo_total_E_dep += hit->E_dep(); + m_calo_total_fit_E_EM += hit->E_EM(); + m_calo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); - if (m_wave_status[ch]&4) { + if (m_wave_status[ch]&4) { // Overflows m_calo_total_E_EM += hit->E_EM(); - } else { + } else { m_calo_total_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); - } + } - ATH_MSG_DEBUG("Calibrated calo: 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()); + // add a fudged energy variable that corrects the MC to agree with the testbeam results + if (isMC) { + // Add fudge factor for MC + m_calo_total_fit_E_EM_fudged = m_calo_total_fit_E_EM * m_caloMC_FudgeFactor; + m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM * m_caloMC_FudgeFactor; + m_calo_total_E_EM_fudged = m_calo_total_E_EM * m_caloMC_FudgeFactor; + } else { + // Unchanged in data + m_calo_total_fit_E_EM_fudged = m_calo_total_fit_E_EM; + m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM; + m_calo_total_E_EM_fudged = m_calo_total_E_EM; + } + } //// the following is an example of how to access the linked waveform data from the calibrated data //auto measurements = &(hit->WaveformLinks())[0]; @@ -1101,16 +1208,60 @@ 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_fit_E_EM_fudged = m_calo_total_fit_E_EM; - m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM; - m_calo_total_E_EM_fudged = m_calo_total_E_EM; - if (isMC) { - m_calo_total_fit_E_EM_fudged = m_calo_total_fit_E_EM_fudged * m_caloMC_FudgeFactor; - m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM_fudged * m_caloMC_FudgeFactor; - m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor; + // load in calibrated calo2 container + SG::ReadHandle<xAOD::CalorimeterHitContainer> ecal2CalibratedContainer { m_ecal2CalibratedContainer, ctx }; + ATH_CHECK(ecal2CalibratedContainer.isValid()); + for (auto hit : *ecal2CalibratedContainer) { + int ch=hit->channel(); + m_calibrated_nMIP[ch] = hit->Nmip(); + m_calibrated_E_dep[ch] = hit->E_dep(); + m_calibrated_E_EM[ch] = hit->E_EM(); + + ATH_MSG_DEBUG("Calibrated calo: 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()); + + m_calohi_total_nMIP += hit->Nmip(); + m_calohi_total_E_dep += hit->E_dep(); + m_calohi_total_fit_E_EM += hit->E_EM(); + m_calohi_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + + if (m_wave_status[ch]&4) { // Overflows + m_calohi_total_E_EM += hit->E_EM(); + } else { + m_calohi_total_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + } } + // load in calibrated calo2 container + SG::ReadHandle<xAOD::CalorimeterHitContainer> calonuCalibratedContainer { m_calonuCalibratedContainer, ctx }; + ATH_CHECK(calonuCalibratedContainer.isValid()); + for (auto hit : *calonuCalibratedContainer) { + int ch=hit->channel(); + m_calibrated_nMIP[ch] = hit->Nmip(); + m_calibrated_E_dep[ch] = hit->E_dep(); + m_calibrated_E_EM[ch] = hit->E_EM(); + + ATH_MSG_DEBUG("Calibrated calo: 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()); + + m_calonu_total_nMIP += hit->Nmip(); + m_calonu_total_E_dep += hit->E_dep(); + m_calonu_total_fit_E_EM += hit->E_EM(); + m_calonu_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + + if (m_wave_status[ch]&4) { // Overflows + m_calonu_total_E_EM += hit->E_EM(); + } else { + m_calonu_total_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio()); + } + } + // load in calibrated preshower container SG::ReadHandle<xAOD::CalorimeterHitContainer> preshowerCalibratedContainer { m_preshowerCalibratedContainer, ctx }; ATH_CHECK(preshowerCalibratedContainer.isValid()); @@ -1127,6 +1278,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const // 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) + // This needs to be fixed!!! bool blinded_event = false; if ((!isMC) && abs(m_distanceToCollidingBCID) <= 1) { if (m_calo_total_E_EM > m_blindingCaloThreshold ) { // energy is in MeV @@ -1727,6 +1879,24 @@ NtupleDumperAlg::clearTree() const m_calo_total_raw_E_EM_fudged=0; m_calo_total_E_EM_fudged=0; + m_calolo_total_nMIP=0; + m_calolo_total_E_dep=0; + m_calolo_total_fit_E_EM=0; + m_calolo_total_raw_E_EM=0; + m_calolo_total_E_EM=0; + + m_calohi_total_nMIP=0; + m_calohi_total_E_dep=0; + m_calohi_total_fit_E_EM=0; + m_calohi_total_raw_E_EM=0; + m_calohi_total_E_EM=0; + + m_calonu_total_nMIP=0; + m_calonu_total_E_dep=0; + m_calonu_total_fit_E_EM=0; + m_calonu_total_raw_E_EM=0; + m_calonu_total_E_EM=0; + m_preshower_total_nMIP=0; m_preshower_total_E_dep=0; diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h index ef4c227c2661a1c1faa115d19b6bb9d2d59f09c8..36b9d17dfc7c5988ac7c4c98a35f5a25224c6261 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h @@ -65,10 +65,11 @@ private: 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); + void addCalibratedBranches(const std::string &name, std::list<int> channel_list) const; + void addCaloSumBranches(const std::string &name) const ; + double radius(const Acts::Vector3 &position) const; // Return map with string as key and list<int> as value @@ -157,6 +158,26 @@ private: double m_baseEventCrossSection {1.0}; const double kfemtoBarnsPerMilliBarn {1.0e12}; + // Flag for dual-gain calorimeter configuration + mutable bool m_dual_gain {false}; + + // Helper lists to hold channel lists + mutable std::list<int> m_calo_channels; + mutable std::list<int> m_calolo_channels; + mutable std::list<int> m_calohi_channels; + mutable std::list<int> m_calonu_channels; + mutable std::list<int> m_veto_channels; + mutable std::list<int> m_vetonu_channels; + mutable std::list<int> m_trigger_channels; + mutable std::list<int> m_preshower_channels; + + // Output statistics + mutable int m_eventsPassed = 0; + mutable int m_eventsFailedGRL = 0; + mutable int m_two_tracks = 0; + mutable int m_verticies = 0; + + // Output tree mutable TTree* m_tree; //mutable unsigned int n_wave_chan; // Actual number of waveform channels @@ -213,19 +234,28 @@ private: mutable float m_calo_total_fit_E_EM_fudged; mutable float m_calo_total_raw_E_EM_fudged; mutable float m_calo_total_E_EM_fudged; + + mutable float m_calolo_total_nMIP; + mutable float m_calolo_total_E_dep; + mutable float m_calolo_total_fit_E_EM; + mutable float m_calolo_total_raw_E_EM; + mutable float m_calolo_total_E_EM; + + mutable float m_calohi_total_nMIP; + mutable float m_calohi_total_E_dep; + mutable float m_calohi_total_fit_E_EM; + mutable float m_calohi_total_raw_E_EM; + mutable float m_calohi_total_E_EM; + + mutable float m_calonu_total_nMIP; + mutable float m_calonu_total_E_dep; + mutable float m_calonu_total_fit_E_EM; + mutable float m_calonu_total_raw_E_EM; + mutable float m_calonu_total_E_EM; mutable float m_preshower_total_nMIP; mutable float m_preshower_total_E_dep; - mutable float m_Calo0_Edep; - mutable float m_Calo1_Edep; - mutable float m_Calo2_Edep; - mutable float m_Calo3_Edep; - - mutable float m_Calo_Total_Edep; - mutable float m_Preshower12_Edep; - mutable float m_Preshower13_Edep; - mutable float m_clock_phase; mutable unsigned int m_station0Clusters; @@ -416,12 +446,6 @@ private: mutable double m_crossSection; mutable std::vector<double> m_genWeights; - - mutable int m_eventsPassed = 0; - mutable int m_eventsFailedGRL = 0; - mutable int m_two_tracks = 0; - mutable int m_verticies = 0; - mutable double m_vertex_x; // components of reconstructed vertex in mm mutable double m_vertex_y; mutable double m_vertex_z;