diff --git a/Calorimeter/CaloDetDescr/CaloIdDictFiles/data/IdDictCalorimeter.xml b/Calorimeter/CaloDetDescr/CaloIdDictFiles/data/IdDictCalorimeter.xml
index 006fabfa2116e1390b3c2f612e3b378e45f5562d..37bed2379f814331055b953646980bbd9e3b444e 100644
--- a/Calorimeter/CaloDetDescr/CaloIdDictFiles/data/IdDictCalorimeter.xml
+++ b/Calorimeter/CaloDetDescr/CaloIdDictFiles/data/IdDictCalorimeter.xml
@@ -21,6 +21,6 @@
<range field="part" value="Ecal" />
<range field="row" values="Bottom Top" wraparound="FALSE" />
<range field="module" values="Starboard Port" wraparound="FALSE" />
- <range field="pmt" minvalue="0" maxvalue="0" />
+ <range field="pmt" minvalue="0" maxvalue="1" />
</region>
-</IdDictionary>
\ No newline at end of file
+</IdDictionary>
diff --git a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
index c49770ad2b7e9468eca196af582d2156b9ab977c..bfa25550254ef72e1499e420d69a0ccde3a9f547 100644
--- a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
+++ b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py
@@ -15,8 +15,11 @@ def CalorimeterReconstructionCfg(flags, **kwargs):
acc = ComponentAccumulator()
kwargs.setdefault("CaloWaveHitContainerKey", "CaloWaveformHits")
+ kwargs.setdefault("Calo2WaveHitContainerKey", "Calo2WaveformHits")
kwargs.setdefault("PreshowerWaveHitContainerKey", "PreshowerWaveformHits")
+
kwargs.setdefault("CaloHitContainerKey", "CaloHits")
+ kwargs.setdefault("Calo2HitContainerKey", "Calo2Hits")
kwargs.setdefault("PreshowerHitContainerKey", "PreshowerHits")
acc.merge(CaloRecToolCfg(flags, **kwargs))
diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
index b6df5fed6053f7499a6f45a95e4b51d2e6633ddd..c8844946d2f8b1242c22ec0540a07d23ff2c7ada 100644
--- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
+++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx
@@ -12,23 +12,19 @@ StatusCode CaloRecAlg::initialize() {
// Set key to read calo hits from
ATH_CHECK( m_caloWaveHitContainerKey.initialize() );
+ ATH_CHECK( m_calo2WaveHitContainerKey.initialize() );
// Set key to read preshower hits from
ATH_CHECK( m_preshowerWaveHitContainerKey.initialize() );
// Set key to write container
ATH_CHECK( m_caloHitContainerKey.initialize() );
+ ATH_CHECK( m_calo2HitContainerKey.initialize() );
ATH_CHECK( m_preshowerHitContainerKey.initialize() );
// Initalize tools
ATH_CHECK( m_recoCalibTool.retrieve() );
- // Store calibrattiion factos in a vector for ease of access
- m_EM_mu_Map[0] = m_calo_ch0_EM_mu;
- m_EM_mu_Map[1] = m_calo_ch1_EM_mu;
- m_EM_mu_Map[2] = m_calo_ch2_EM_mu;
- m_EM_mu_Map[3] = m_calo_ch3_EM_mu;
-
return StatusCode::SUCCESS;
}
//----------------------------------------------------------------------
@@ -54,6 +50,15 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Calorimeter Waveform Hit container found with zero length!");
}
+ 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() );
@@ -63,116 +68,55 @@ StatusCode CaloRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Preshower Waveform Hit container found with zero length!");
}
- // Find the output waveform container
+ // Only correct gain in real data
+ bool correct_gain = !m_isMC;
+
+ // Reconstruct calorimeter hits
SG::WriteHandle<xAOD::CalorimeterHitContainer> caloHitContainerHandle(m_caloHitContainerKey, ctx);
ATH_CHECK( caloHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(),
std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) );
-
- SG::WriteHandle<xAOD::CalorimeterHitContainer> preshowerHitContainerHandle(m_preshowerHitContainerKey, ctx);
- ATH_CHECK( preshowerHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(),
- std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) );
-
ATH_MSG_DEBUG("WaveformsHitContainer '" << caloHitContainerHandle.name() << "' initialized");
- ATH_MSG_DEBUG("WaveformsHitContainer '" << preshowerHitContainerHandle.name() << "' initialized");
// Loop over calo hits and calibrate each primary hit
for( const auto& hit : *caloWaveHitHandle ) {
if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue;
-
- // Create a new calo hit
xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit();
caloHitContainerHandle->push_back(calo_hit);
+ calo_hit->addHit(caloWaveHitHandle.get(), hit);
+ m_recoCalibTool->reconstruct(ctx, calo_hit, correct_gain);
+ }
+ ATH_MSG_DEBUG("CaloHitContainer '" << caloHitContainerHandle.name() << "' filled with "<< caloHitContainerHandle->size() <<" items");
- ATH_MSG_DEBUG("calo_hit in channel " << hit->channel() );
-
- float MIPcharge_ref = m_recoCalibTool->getMIPcharge_ref(hit->channel()); // get reference MIP charge from database
-
- float charge = hit->integral()/50.0; // divide by 50 ohms to get charge
- ATH_MSG_DEBUG("calo_hit filled has charge of " << charge << " pC");
-
- float gainRatio = 1.0;
- if (!m_isMC) { // MC already has correct MIP charge stored in MIPcharge_ref, so only need to to HV extrapolation with reral data
- gainRatio = extrapolateHVgain(hit->channel());
- }
- ATH_MSG_DEBUG("HV gain ratio = " << gainRatio );
-
-
- float Nmip = (charge * gainRatio) / MIPcharge_ref;
- ATH_MSG_DEBUG("Nmip = " << Nmip );
- calo_hit->set_Nmip(Nmip); // set Nmip value
-
- float E_dep = Nmip * m_MIP_sim_Edep_calo;
- ATH_MSG_DEBUG("E_dep in MeV = " << E_dep );
- calo_hit->set_E_dep(E_dep); // set calibrated E_dep value
-
- float E_EM = Nmip * m_EM_mu_Map[hit->channel()];
- ATH_MSG_DEBUG("Em E in MeV = " << E_EM );
- calo_hit->set_E_EM(E_EM); // set calibrated E_EM value
-
- float fit_to_raw_ratio = 1.0;
- if (hit->integral() != 0.0) { // avoid possibility of division by zero error
- fit_to_raw_ratio = hit->raw_integral() / hit->integral();
- }
- calo_hit->set_fit_to_raw_ratio(fit_to_raw_ratio); // set fit-to-raw-ratio that can be used to take any of the calibrated values to what they would be if we used the raw integral instead of the fit integral
-
- calo_hit->set_channel(hit->channel()); // set channel number
+ // 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");
- calo_hit->clearWaveformLinks();
- calo_hit->addHit(caloWaveHitHandle.get(), hit); // create link to calo waveform hit
+ 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");
- ATH_MSG_DEBUG("CaloHitContainer '" << caloHitContainerHandle.name() << "' filled with "<< caloHitContainerHandle->size() <<" items");
+ // Reconstruct preshower hits
+ SG::WriteHandle<xAOD::CalorimeterHitContainer> preshowerHitContainerHandle(m_preshowerHitContainerKey, ctx);
+ ATH_CHECK( preshowerHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(),
+ std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) );
+ ATH_MSG_DEBUG("WaveformsHitContainer '" << preshowerHitContainerHandle.name() << "' initialized");
for( const auto& hit : *preshowerWaveHitHandle ) {
if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue;
-
- // Create a new preshower hit
- xAOD::CalorimeterHit* preshower_hit = new xAOD::CalorimeterHit();
- preshowerHitContainerHandle->push_back(preshower_hit);
-
- ATH_MSG_DEBUG("preshower_hit in channel " << hit->channel() );
-
- float MIPcharge_ref = m_recoCalibTool->getMIPcharge_ref(hit->channel()); // get reference MIP charge from database
-
- float charge = hit->integral()/50.0; // divide by 50 ohms to get charge
- ATH_MSG_DEBUG("preshower_hit filled has charge of " << charge << " pC");
-
- float Nmip = charge / MIPcharge_ref;
- ATH_MSG_DEBUG("Nmip = " << Nmip );
- preshower_hit->set_Nmip(Nmip); // set Nmip value
-
- float E_dep = Nmip * m_MIP_sim_Edep_preshower;
- ATH_MSG_DEBUG("E_dep in GeV = " << E_dep );
- preshower_hit->set_E_dep(E_dep); // set calibrated E_dep value
-
- float fit_to_raw_ratio = 1.0;
- if (hit->integral() != 0.0) { // avoid possibility of division by zero error
- fit_to_raw_ratio = hit->raw_integral() / hit->integral();
- }
- preshower_hit->set_fit_to_raw_ratio(fit_to_raw_ratio); // set fit-to-raw-ratio that can be used to take any of the calibrated values to what they would be if we used the raw integral instead of the fit integral
-
- preshower_hit->set_channel(hit->channel()); // set channel number
-
- preshower_hit->clearWaveformLinks();
- preshower_hit->addHit(preshowerWaveHitHandle.get(), hit); // create link to preshower waveform hit
- }
-
+ xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit();
+ preshowerHitContainerHandle->push_back(calo_hit);
+ calo_hit->addHit(preshowerWaveHitHandle.get(), hit);
+ m_recoCalibTool->reconstruct(ctx, calo_hit, correct_gain=false);
+ }
ATH_MSG_DEBUG("PreshowerHitContainer '" << preshowerHitContainerHandle.name() << "' filled with "<< preshowerHitContainerHandle->size() <<" items");
-
-
+
return StatusCode::SUCCESS;
}
-//----------------------------------------------------------------------
-float CaloRecAlg::extrapolateHVgain(int channel) const {
- float PMT_hv = m_recoCalibTool->getHV(channel);
- float PMT_hv_ref = m_recoCalibTool->getHV_ref(channel);
- TF1 gaincurve = m_recoCalibTool->get_PMT_HV_curve(channel);
-
- float gaincurve_atHV = gaincurve.Eval(PMT_hv);
- float gaincurve_atHVref = gaincurve.Eval(PMT_hv_ref);
-
- return ( gaincurve_atHVref / gaincurve_atHV ) * pow( PMT_hv_ref / PMT_hv , 6.6);
-}
-//----------------------------------------------------------------------
-
diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h
index 721f521621ae519778a8ddadd8b5711f1beeac9b..b11fc05ae5cf91b8e5b4aed67074d07674a511a7 100644
--- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h
+++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h
@@ -71,6 +71,10 @@ class CaloRecAlg : public AthReentrantAlgorithm {
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_caloWaveHitContainerKey {this, "CaloWaveHitContainerKey", "CaloWaveformHits"};
//@}
+ //@{
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_calo2WaveHitContainerKey {this, "Calo2WaveHitContainerKey", "Calo2WaveformHits"};
+ //@}
+
//@{
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerWaveHitContainerKey {this, "PreshowerWaveHitContainerKey", "PreshowerWaveformHits"};
//@}
@@ -80,21 +84,10 @@ class CaloRecAlg : public AthReentrantAlgorithm {
*/
//@{
SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_caloHitContainerKey {this, "CaloHitContainerKey", "CaloHits"};
+ SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_calo2HitContainerKey {this, "Calo2HitContainerKey", "Calo2Hits"};
SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_preshowerHitContainerKey {this, "PreshowerHitContainerKey", "PreshowerHits"};
//@}
- float extrapolateHVgain(int channel) const;
-
- 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_preshower {this, "MIP_sim_Edep_preshower", 4.894}; // MIP deposits 4.894 MeV of energy in a preshower layer
-
- FloatProperty m_calo_ch0_EM_mu {this, "m_calo_ch0_EM_mu", 330.0}; // factor used to do rough calibration of calo ch0 to EM energy: 0.33 GeV or 330 MeV
- FloatProperty m_calo_ch1_EM_mu {this, "m_calo_ch1_EM_mu", 330.0}; // factor used to do rough calibration of calo ch1 to EM energy: 0.33 GeV or 330 MeV
- FloatProperty m_calo_ch2_EM_mu {this, "m_calo_ch2_EM_mu", 330.0}; // factor used to do rough calibration of calo ch2 to EM energy: 0.33 GeV or 330 MeV
- FloatProperty m_calo_ch3_EM_mu {this, "m_calo_ch3_EM_mu", 330.0}; // factor used to do rough calibration of calo ch3 to EM energy: 0.33 GeV or 330 MeV
-
- float m_EM_mu_Map[4]; // vector that holds EM_mu calibration factors for calo channels
-
Gaudi::Property<bool> m_isMC {this, "isMC", false};
};
diff --git a/Calorimeter/CaloRecTools/CMakeLists.txt b/Calorimeter/CaloRecTools/CMakeLists.txt
index f0e9f07b53f5426b360c2639d4f7953a44d57792..abfdf9927c6467cbd3f8e847c09a1b799b0b327b 100644
--- a/Calorimeter/CaloRecTools/CMakeLists.txt
+++ b/Calorimeter/CaloRecTools/CMakeLists.txt
@@ -13,7 +13,8 @@ atlas_add_library( CaloRecToolsLib
CaloRecTools/*.h src/*.cxx src/*.h
PUBLIC_HEADERS CaloRecTools
PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES AthenaBaseComps AthenaPoolUtilities AthenaKernel GaudiKernel
+ LINK_LIBRARIES xAODFaserWaveform xAODFaserCalorimeter
+ AthenaBaseComps AthenaPoolUtilities AthenaKernel GaudiKernel
PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
)
diff --git a/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h b/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h
index feb141d915d6214a681f437d94b8118e14c33820..e52efd93dd560bc38a83b9a23c10c16e6df533d9 100644
--- a/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h
+++ b/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h
@@ -17,6 +17,10 @@
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/EventContext.h"
+// Data Classes
+// #include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserCalorimeter/CalorimeterHit.h"
+
#include "TF1.h"
///Interface for Calo reco algorithms
@@ -40,6 +44,11 @@ class ICaloRecTool : virtual public IAlgTool
virtual float getMIPcharge_ref(int channel) const = 0;
virtual TF1 get_PMT_HV_curve(int channel) const = 0;
+
+ // Reconstruct one waveform hit and put the resulting Calorimeter hit in the container
+ virtual void reconstruct(const EventContext& ctx,
+ xAOD::CalorimeterHit* hit,
+ bool correct_gain) const = 0;
};
#endif // CALORECTOOL_ICALORECTOOL_H
diff --git a/Calorimeter/CaloRecTools/src/CaloRecTool.cxx b/Calorimeter/CaloRecTools/src/CaloRecTool.cxx
index 916011019805ab9087c57427859086dec456eb80..01a81e287ca312ecde8616cf22a1c23a29186767 100644
--- a/Calorimeter/CaloRecTools/src/CaloRecTool.cxx
+++ b/Calorimeter/CaloRecTools/src/CaloRecTool.cxx
@@ -40,6 +40,29 @@ CaloRecTool::initialize() {
HVgaincurves_rootFile->Close(); // close the root file
+ // These should be in DB, but just hardcode this for now
+ m_MIP_sim_Edep[0] = m_MIP_sim_Edep_calo.value();
+ m_MIP_sim_Edep[1] = m_MIP_sim_Edep_calo.value();
+ m_MIP_sim_Edep[2] = m_MIP_sim_Edep_calo.value();
+ m_MIP_sim_Edep[3] = m_MIP_sim_Edep_calo.value();
+
+ m_MIP_sim_Edep[12] = m_MIP_sim_Edep_preshower.value();
+ m_MIP_sim_Edep[13] = m_MIP_sim_Edep_preshower.value();
+
+ m_MIP_sim_Edep[16] = m_MIP_sim_Edep_calo2.value();
+ m_MIP_sim_Edep[17] = m_MIP_sim_Edep_calo2.value();
+ m_MIP_sim_Edep[18] = m_MIP_sim_Edep_calo2.value();
+ m_MIP_sim_Edep[19] = m_MIP_sim_Edep_calo2.value();
+
+ m_EM_mu_Map[0] = m_calo_EM_mu.value();
+ m_EM_mu_Map[1] = m_calo_EM_mu.value();
+ m_EM_mu_Map[2] = m_calo_EM_mu.value();
+ m_EM_mu_Map[3] = m_calo_EM_mu.value();
+
+ m_EM_mu_Map[16] = m_calo_EM_mu.value();
+ m_EM_mu_Map[17] = m_calo_EM_mu.value();
+ m_EM_mu_Map[18] = m_calo_EM_mu.value();
+ m_EM_mu_Map[19] = m_calo_EM_mu.value();
return StatusCode::SUCCESS;
}
@@ -150,7 +173,7 @@ float CaloRecTool::getMIPcharge_ref(const EventContext& ctx, int channel) const
ATH_MSG_DEBUG("in getMIPcharge_ref("<<channel<<")");
- float MIP_charge_ref =0.;
+ float MIP_charge_ref =1.; // Default for no calibration parameters
// Read Cond Handle
SG::ReadCondHandle<CondAttrListCollection> readHandle{m_MIP_ref_ReadKey, ctx};
@@ -187,8 +210,63 @@ float CaloRecTool::getMIPcharge_ref(int channel) const {
return CaloRecTool::getMIPcharge_ref(ctx, channel);
}
+//----------------------------------------------------------------------
+float CaloRecTool::extrapolateHVgain(const EventContext& ctx, int channel) const {
+ float PMT_hv = getHV(ctx, channel);
+ float PMT_hv_ref = getHV_ref(ctx, channel);
+ TF1 gaincurve = get_PMT_HV_curve(channel);
+
+ float gaincurve_atHV = gaincurve.Eval(PMT_hv);
+ float gaincurve_atHVref = gaincurve.Eval(PMT_hv_ref);
+
+ return ( gaincurve_atHVref / gaincurve_atHV ) * pow( PMT_hv_ref / PMT_hv , 6.6);
+}
+
//--------------------------------------------------------------
+//----------------------------------------------------------------------
+// Reconstruct one waveform hit
+//xAOD::CalorimeterHit*
+void
+CaloRecTool::reconstruct(const EventContext& ctx,
+ xAOD::CalorimeterHit* calo_hit,
+ bool correct_gain=false) const {
+
+ // Get the waveform hit attached to this calo hit
+ const xAOD::WaveformHit* wave_hit = calo_hit->Hit(0);
+ ATH_MSG_DEBUG("calo_hit in channel " << wave_hit->channel() );
+
+ float MIPcharge_ref = getMIPcharge_ref(ctx, wave_hit->channel()); // get reference MIP charge from database
+
+ float charge = wave_hit->integral()/50.0; // divide by 50 ohms to get charge
+ ATH_MSG_DEBUG("calo_hit filled has charge of " << charge << " pC");
+
+ float gainRatio = 1.0;
+ if (correct_gain) { // MC already has correct MIP charge stored in MIPcharge_ref, so only need to to HV extrapolation with reral data
+ gainRatio = extrapolateHVgain(ctx, wave_hit->channel());
+ }
+ ATH_MSG_DEBUG("HV gain ratio = " << gainRatio );
+
+ float Nmip = (charge * gainRatio) / MIPcharge_ref;
+ ATH_MSG_DEBUG("Nmip = " << Nmip );
+ calo_hit->set_Nmip(Nmip); // set Nmip value
+
+ float E_dep = Nmip * m_MIP_sim_Edep[wave_hit->channel()];
+ ATH_MSG_DEBUG("E_dep in MeV = " << E_dep );
+ calo_hit->set_E_dep(E_dep); // set calibrated E_dep value
+
+ float E_EM = Nmip * m_EM_mu_Map[wave_hit->channel()];
+ ATH_MSG_DEBUG("Em E in MeV = " << E_EM );
+ calo_hit->set_E_EM(E_EM); // set calibrated E_EM value
+
+ float fit_to_raw_ratio = 1.0;
+ if (wave_hit->integral() != 0.0) { // avoid possibility of division by zero error
+ fit_to_raw_ratio = wave_hit->raw_integral() / wave_hit->integral();
+ }
+ calo_hit->set_fit_to_raw_ratio(fit_to_raw_ratio); // set fit-to-raw-ratio that can be used to take any of the calibrated values to what they would be if we used the raw integral instead of the fit integral
+
+}
+
diff --git a/Calorimeter/CaloRecTools/src/CaloRecTool.h b/Calorimeter/CaloRecTools/src/CaloRecTool.h
index 6ab891ad42784e3388d933f42de9267545fab3d9..de736af74393bf424b2fd471795aea7bc2808f52 100644
--- a/Calorimeter/CaloRecTools/src/CaloRecTool.h
+++ b/Calorimeter/CaloRecTools/src/CaloRecTool.h
@@ -25,6 +25,10 @@
#include "Gaudi/Property.h"
#include "GaudiKernel/EventContext.h"
+// Data Classes
+#include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserCalorimeter/CalorimeterHit.h"
+
// Include ROOT classes
#include "TF1.h"
#include "TFile.h"
@@ -57,6 +61,11 @@ class CaloRecTool: public extends<AthAlgTool, ICaloRecTool> {
// method for returning PMT HV calibration curves from root file
virtual TF1 get_PMT_HV_curve(int channel) const override;
+ // Reconstruct one waveform hit and put the resulting Calorimeter hit in the container
+ virtual void reconstruct(const EventContext& ctx,
+ xAOD::CalorimeterHit* hit,
+ bool correct_gain) const override;
+
TFile* HVgaincurves_rootFile;
TF1* chan0_HVgaincurve_pntr;
@@ -80,6 +89,18 @@ class CaloRecTool: public extends<AthAlgTool, ICaloRecTool> {
SG::ReadCondHandleKey<CondAttrListCollection> m_PMT_HV_ReadKey{this, "PMT_HV_ReadKey", "/WAVE/Calibration/HV", "Key of folder for PMT HV reading"};
SG::ReadCondHandleKey<CondAttrListCollection> m_MIP_ref_ReadKey{this, "MIP_ref_ReadKey", "/WAVE/Calibration/MIP_ref", "Key of folder for MIP charge calibration measurment, also stores PMT HV used to measure the reference MIP charge"};
+ // 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_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
+
+ float m_MIP_sim_Edep[32]; // vector that holds Edep factors for calo and preshower
+ float m_EM_mu_Map[32]; // vector that holds EM_mu calibration factors for calo channels
+
+ float extrapolateHVgain(const EventContext& ctx, int channel) const;
+
};
#endif // CALORECTOOLS_CALORECTOOL_H
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index ac3c57ce5fe602f9f817f8092ef0112ae1217a6e..de2d757a56a56b0c1f782a9c8add40fd8a28e062 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -27,7 +27,8 @@ parser.add_argument("--merge", type=int, default=1,
help="Specify merged files per reco file (MC only)")
parser.add_argument("--last", type=int, default=0,
help="Specify last file in slice (normally --files)")
-
+parser.add_argument("-g", "--geom", default="",
+ help="Specify geometry (if it can't be parsed from run number\n Values: Ti12Data03 (2022 TiT12)")
parser.add_argument("--outfile", default="",
help="Override output file name")
@@ -244,9 +245,9 @@ if len(args.geom) > 0:
runtype = args.geom
else:
if args.isMC:
- runtype = TI12MC04
+ runtype = "TI12MC04"
else:
- runtype = TI12Data04
+ runtype = "TI12Data04"
if runtype in ["TI12Data04", "TI12MC04"]:
ConfigFlags.GeoModel.FaserVersion = "FASERNU-04" # FASER geometry
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index def572ed1dbd7ad5e6af132eaeb17fbfc5a57f11..1692580349dde4d73684761e4a9dd33f81530274 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -179,13 +179,15 @@ StatusCode NtupleDumperAlg::initialize()
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",1,14);
+ addWaveBranches("VetoSt1",2,14);
addWaveBranches("VetoSt2",2,6);
addWaveBranches("Timing",4,8);
addWaveBranches("Preshower",2,12);
addWaveBranches("Calo",4,0);
-
+ addWaveBranches("CaloHi", 4,16);
+
m_tree->Branch("ScintHit", &m_scintHit);
addCalibratedBranches("Calo",4,0);
@@ -417,6 +419,12 @@ StatusCode NtupleDumperAlg::initialize()
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]));
@@ -433,6 +441,11 @@ StatusCode NtupleDumperAlg::initialize()
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]));
if (m_onlyBlinded){
ATH_MSG_INFO("Only events that would be blinded are saved in ntuple");
@@ -577,7 +590,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// for random (only) triggers, store charge of scintillators in histograms
if (m_tap == 16) {
// Fill histograms
- for (int chan = 0; chan<15; chan++) {
+ for (unsigned int chan = 0; chan<nchan; chan++) {
m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]);
}
}
@@ -1381,7 +1394,7 @@ NtupleDumperAlg::clearTree() const
m_inputBits=0;
m_inputBitsNext=0;
- for(int ii=0;ii<15;ii++) {
+ for(unsigned int ii=0;ii<nchan;ii++) {
m_wave_localtime[ii]=0;
m_wave_peak[ii]=0;
m_wave_width[ii]=0;
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 5f7dedda2c7801f58ecdd5365dc152dbdab763e7..d636d8ef70dcdd7fbcb90fd44eddccd8c2ac36ac 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -84,8 +84,10 @@ private:
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
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::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_preshowerCalibratedContainer { this, "preshowerCalibratedContainer", "PreshowerHits", "Preshower Calibrated hit container name" };
SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer", "Tracker cluster container name" };
@@ -146,7 +148,9 @@ private:
mutable TTree* m_tree;
- mutable TH1* m_HistRandomCharge[15];
+ const static unsigned int nchan=20;
+
+ mutable TH1* m_HistRandomCharge[nchan];
mutable unsigned int m_run_number;
mutable unsigned int m_event_number;
@@ -170,22 +174,22 @@ private:
mutable unsigned int m_inputBits;
mutable unsigned int m_inputBitsNext;
- mutable float m_wave_localtime[15];
- mutable float m_wave_peak[15];
- mutable float m_wave_width[15];
- mutable float m_wave_charge[15];
+ 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[15];
- mutable float m_wave_raw_charge[15];
- mutable float m_wave_baseline_mean[15];
- mutable float m_wave_baseline_rms[15];
- mutable unsigned int m_wave_status[15];
+ 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 unsigned int m_scintHit;
- mutable float m_calibrated_nMIP[15];
- mutable float m_calibrated_E_dep[15];
- mutable float m_calibrated_E_EM[15];
+ mutable float m_calibrated_nMIP[nchan];
+ mutable float m_calibrated_E_dep[nchan];
+ mutable float m_calibrated_E_EM[nchan];
mutable float m_calo_total_nMIP;
mutable float m_calo_total_E_dep;
@@ -204,6 +208,7 @@ private:
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;
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
index f7ef9339ae32919ef6600cf099da2d42c5eb4a69..1caf49d4a3be081da8c8b8339324125ee6df788b 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
@@ -96,6 +96,14 @@ WaveformReconstructionTool::reconstructPrimary(
int lo_edge = int((trigger_time+offset)/2.) + m_windowStart;
int hi_edge = int((trigger_time+offset)/2.) + m_windowStart + m_windowWidth;
+ if (hi_edge >= int(wave.size())) {
+ // This likely means we have the wrong digitizer range
+ ATH_MSG_WARNING("Found channel " << wave.channel() << " with low edge: " << lo_edge << " hi edge: " << hi_edge << " > wave.size() " << wave.size());
+ ATH_MSG_WARNING(" trigger_time + offset: " << (trigger_time+offset) << " => " << int((trigger_time+offset)/2.));
+ newhit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_INVALID);
+ return StatusCode::SUCCESS;
+ }
+
// Fill raw hit values
fillRawHitValues(wave, lo_edge, hi_edge, newhit);
@@ -367,7 +375,11 @@ WaveformReconstructionTool::fillRawHitValues(const RawWaveform& wave,
// First, make sure we don't overflow the waveform range
if (lo_edge < 0) lo_edge = 0;
if (hi_edge >= int(wave.size())) hi_edge = wave.size() - 1;
-
+ if (hi_edge <= lo_edge) {
+ ATH_MSG_WARNING("Channel " << wave.channel() << " waveform with lo: " << lo_edge << " hi: " << hi_edge << "!");
+ return;
+ }
+
ATH_MSG_DEBUG("Fill channel " << wave.channel()
<< " waveform from sample " << lo_edge << " to " << hi_edge);
diff --git a/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py
index 57c3390c087b7f41850716d6ec24aa7d7e7044f6..0074c2b69906eda40667defdebd2596ff96d8652 100755
--- a/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py
+++ b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py
@@ -15,7 +15,7 @@ nominal_data = {
6525: 820.
}
-offset_channels = 16
+# offset_channels = 20
# Run
# 0 - initial data
@@ -53,7 +53,10 @@ offset_data = {
5042: ehn1_offsets,
5050: ti12_offsets,
# IFT and VetoNu installed
- 6525: [ -10., -10., -10., -10., -25., -25., 0., 0., 0., 0., 0., 0., 18., 18., 0., 0. ]
+ 6525: [ -10., -10., -10., -10., -25., -25., 0., 0., 0., 0., 0., 0., 18., 18., 0., 0. ],
+# 2024, add 2nd digitizer
+ 13847: [ -10., -10., -10., -10., -25., -25., 0., 0., 0., 0., 0., 0., 18., 18., 0., 0., -10., -10., -10., -10. ]
+
}
attr_list_desc = '<timeStamp>run-lumi</timeStamp><addrHeader><address_header service_type="71" clid="40774348" /></addrHeader><typeName>AthenaAttributeList</typeName>'
@@ -82,13 +85,12 @@ for run, data in offset_data.items():
assert run > lastRun, 'Run numbers out of order'
assert run <= maxInt32, 'Run number out of range'
lastRun = run
- assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
- for i in range(offset_channels):
+ #assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
+ for i in range(len(data)):
assert isinstance(data[i], float), 'Offset time is not float'
# Data looks OK
-
from PyCool import cool
from CoolConvUtilities.AtlCoolLib import indirectOpen
@@ -135,7 +137,7 @@ offsetFolder = db.createFolder('/WAVE/DAQ/TimingOffset', offsetFolderSpec, cond_
lastValid = cool.ValidityKeyMax
for firstValidRun, offset_list in reversed(offset_data.items()):
firstValid = (firstValidRun << 32)
- for channel in range(offset_channels):
+ for channel in range(len(offset_list)):
offsetRecord = cool.Record(offsetSpec)
offsetRecord[ 'TriggerOffset' ] = float(offset_list[channel])
offsetFolder.storeObject( firstValid, lastValid, offsetRecord, cool.ChannelId(channel) )
@@ -154,8 +156,9 @@ nominal_data = {
0: 820.
}
# No offsets by default
+# Expand this to 20 channels for 2024 data
offset_data = {
- 0: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
+ 0: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
}
# Validate again just in case
@@ -176,8 +179,8 @@ for run, data in offset_data.items():
assert run > lastRun, 'Run numbers out of order'
assert run <= maxInt32, 'Run number out of range'
lastRun = run
- assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
- for i in range(offset_channels):
+ #assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
+ for i in range(len(data)):
assert isinstance(data[i], float), 'Offset time is not float'
# Data looks OK
@@ -224,7 +227,7 @@ offsetFolder = db.createFolder('/WAVE/DAQ/TimingOffset', offsetFolderSpec, cond_
lastValid = cool.ValidityKeyMax
for firstValidRun, offset_list in reversed(offset_data.items()):
firstValid = (firstValidRun << 32)
- for channel in range(offset_channels):
+ for channel in range(len(offset_list)):
offsetRecord = cool.Record(offsetSpec)
offsetRecord[ 'TriggerOffset' ] = float(offset_list[channel])
offsetFolder.storeObject( firstValid, lastValid, offsetRecord, cool.ChannelId(channel) )
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
index a7de79076eddf9fa132849deb24089dacbf996dd..2828113ed35bf6d24c19ccc31ca3470fd0a39d98 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
@@ -94,10 +94,17 @@ WaveformCableMappingTool::getCableMapping(const EventContext& ctx) const {
// Ugh, cant use switch statement with strings
// Must do this using an if ladder
if (det_type == "calo") {
- identifier = m_ecalID->pmt_id(rowVal, moduleVal, pmtVal);
+ // Do checks of PMT identifier
+ identifier = m_ecalID->pmt_id(rowVal, moduleVal, pmtVal, true);
+ // Test a few things
+ ATH_MSG_DEBUG("Calo ID:" << identifier);
+ ATH_MSG_DEBUG("PMT:" << m_ecalID->pmt(identifier) << " PMT Max:" << m_ecalID->pmt_max(identifier));
}
else if (det_type == "calo2") {
- identifier = m_ecalID->pmt_id(rowVal, moduleVal, pmtVal);
+ // Do checks of PMT identifier
+ identifier = m_ecalID->pmt_id(rowVal, moduleVal, pmtVal, true);
+ ATH_MSG_DEBUG("Calo2 ID:" << identifier);
+ ATH_MSG_DEBUG("PMT:" << m_ecalID->pmt(identifier) << " PMT Max:" << m_ecalID->pmt_max(identifier));
}
else if (det_type == "veto") {
identifier = m_vetoID->pmt_id(stationVal, plateVal, pmtVal);