From 31d8595ebcc32204dfd224e93cc87b3a8da4110e Mon Sep 17 00:00:00 2001
From: Peter Alan Steinberg <peter.steinberg@bnl.gov>
Date: Sat, 29 Oct 2016 14:52:40 +0200
Subject: [PATCH] 'CalibEnergy and CalibTime are in EDM even when calibrations
off, to maintain format, but set to -1000 each to avoid anyone using them'
(ZdcAnalysis-00-00-16)
---
.../ZDC/ZdcAnalysis/CMakeLists.txt | 28 +
.../ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx | 315 ++++++
.../ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx | 150 +++
.../ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx | 506 ++++++++++
.../ZdcAnalysis/Root/ZDCTriggerEfficiency.cxx | 63 ++
.../ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx | 947 ++++++++++++++++++
.../ZdcAnalysis/Root/ZdcAnalysis_entries.cxx | 16 +
.../ZDC/ZdcAnalysis/Root/ZdcSincInterp.cxx | 27 +
.../ZdcAnalysis/IZdcAnalysisTool.h | 29 +
.../ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h | 153 +++
.../ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h | 255 +++++
.../ZdcAnalysis/ZDCPulseAnalyzer.h | 290 ++++++
.../ZdcAnalysis/ZDCTriggerEfficiency.h | 101 ++
.../ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h | 139 +++
.../ZdcAnalysis/ZdcAnalysis/ZdcSincInterp.h | 12 +
.../ZDC/ZdcAnalysis/cmt/Makefile.RootCore | 60 ++
.../ZDC/ZdcAnalysis/cmt/requirements | 24 +
.../ZdcAnalysis/data/ZdcAnalysisConfig.conf | 3 +
.../ZdcAnalysis/tools/RunZDCTreeAnalysis.C | 216 ++++
.../ZdcAnalysis/tools/ZDCNLCalibration.cxx | 603 +++++++++++
.../ZDC/ZdcAnalysis/tools/ZDCNLCalibration.h | 286 ++++++
.../ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.C | 412 ++++++++
.../ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.h | 512 ++++++++++
23 files changed, 5147 insertions(+)
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/CMakeLists.txt
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCTriggerEfficiency.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysis_entries.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcSincInterp.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/IZdcAnalysisTool.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCTriggerEfficiency.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcSincInterp.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/cmt/requirements
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/data/ZdcAnalysisConfig.conf
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/tools/RunZDCTreeAnalysis.C
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.cxx
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.h
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.C
create mode 100644 ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.h
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/CMakeLists.txt b/ForwardDetectors/ZDC/ZdcAnalysis/CMakeLists.txt
new file mode 100644
index 000000000000..f9fdbeb3d980
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/CMakeLists.txt
@@ -0,0 +1,28 @@
+atlas_subdir( ZDCAnalysis )
+
+
+atlas_depends_on_subdirs( PUBLIC
+ Control/AthToolSupport/AsgTools
+ External/GaudiInterface
+ Event/xAOD/xAODForward
+ Event/xAOD/xAODTrigL1Calo
+ PRIVATE
+ Event/xAOD/xAODEventInfo
+ Tools/PathResolver)
+
+find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread MathMore Minuit Minuit2 Matrix Physics HistPainter Rint Graf Graf3d Gpad Html Postscript Gui GX11TTF GX11 )
+
+atlas_add_component( ZDCAnalysis
+ Root/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODForward xAODTrigL1Calo xAODEventInfo
+ )
+
+atlas_add_library (ZDCAnalysisLib
+ Root/*.cxx
+ PUBLIC_HEADERS ZDCAnalysis
+ PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODForward xAODTrigL1Calo xAODEventInfo
+ )
+
+atlas_install_joboptions( share/*.py)
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
new file mode 100644
index 000000000000..8f989dcfb888
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
@@ -0,0 +1,315 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <ZdcAnalysis/ZDCDataAnalyzer.h>
+#include <ZdcAnalysis/ZDCPulseAnalyzer.h>
+
+#include <sstream>
+
+ZDCDataAnalyzer::ZDCDataAnalyzer(int nSample, float deltaTSample, size_t preSampleIdx, std::string fitFunction,
+ const ZDCModuleFloatArray& peak2ndDerivMinSamples,
+ const ZDCModuleFloatArray& peak2ndDerivMinThresholdsHG,
+ const ZDCModuleFloatArray& peak2ndDerivMinThresholdsLG,
+ bool forceLG) :
+ _nSample(nSample), _deltaTSample(deltaTSample), _preSampleIdx(preSampleIdx),
+ _fitFunction(fitFunction),
+ _forceLG(forceLG),
+ _debugLevel(-1),
+ _eventCount(0),
+ _moduleMask(0),
+ _moduleSum({0, 0}),
+ _moduleSumErrSq({0, 0}),
+ _moduleSumPreSample({0,0}),
+ _calibModuleSum({0, 0}),
+ _calibModuleSumErrSq({0,0}),
+ _averageTime({0, 0}),
+ _fail({false, false}),
+ _haveT0Calib(false),
+ _haveECalib(false),
+ _currentLB(-1)
+{
+ _moduleAnalyzers[0] = {0, 0, 0, 0};
+ _moduleAnalyzers[1] = {0, 0, 0, 0};
+
+ _calibAmplitude[0] = {0, 0, 0, 0};
+ _calibAmplitude[1] = {0, 0, 0, 0};
+
+ _calibTime[0] = {0, 0, 0, 0};
+ _calibTime[1] = {0, 0, 0, 0};
+
+ // For now we are using hard-coded gain factors and pedestals
+ //
+ _HGGains[0] = {9.51122, 9.51980, 9.51122, 9.51122};
+ _HGGains[1] = {9.50842, 9.49662, 9.50853, 9.50842};
+
+ _pedestals[0] = {100, 100, 100, 100};
+ _pedestals[1] = {100, 100, 100, 100};
+
+ // Default "calibrations"
+ //
+ _currentECalibCoeff = {1, 1, 1, 1, 1, 1, 1, 1};
+
+ _currentT0OffsetsHG = {0, 0, 0, 0, 0, 0, 0, 0};
+ _currentT0OffsetsLG = {0, 0, 0, 0, 0, 0, 0, 0};
+
+ // Construct the per-module pulse analyzers
+ //
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ std::ostringstream moduleTag;
+ moduleTag << "_s" << side << "_m" << module;
+
+ _moduleAnalyzers[side][module] = new ZDCPulseAnalyzer(moduleTag.str().c_str(), _nSample, _deltaTSample, _preSampleIdx,
+ _pedestals[side][module], _HGGains[side][module], _fitFunction,
+ peak2ndDerivMinSamples[side][module],
+ peak2ndDerivMinThresholdsHG[side][module],
+ peak2ndDerivMinThresholdsLG[side][module]);
+ if (_forceLG) _moduleAnalyzers[side][module]->SetForceLG(true);
+ }
+ }
+
+ // By default we perform quiet pulse fits
+ //
+ ZDCPulseAnalyzer::SetQuietFits(true);
+}
+
+ZDCDataAnalyzer::~ZDCDataAnalyzer()
+{
+ // Delete the per-module pulse analyzers
+ //
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ if (!_moduleAnalyzers[side][module]) delete _moduleAnalyzers[side][module];
+ }
+ }
+}
+
+void ZDCDataAnalyzer::SetTauT0Values(const ZDCModuleBoolArray& fixTau1, const ZDCModuleBoolArray& fixTau2,
+ const ZDCModuleFloatArray& tau1, const ZDCModuleFloatArray& tau2,
+ const ZDCModuleFloatArray& t0HG, const ZDCModuleFloatArray& t0LG)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetTauT0Values(fixTau1[side][module], fixTau2[side][module],
+ tau1[side][module], tau2[side][module], t0HG[side][module], t0LG[side][module]);
+ }
+ }
+}
+
+void ZDCDataAnalyzer::SetADCOverUnderflowValues(const ZDCModuleFloatArray& HGOverflowADC, const ZDCModuleFloatArray& HGUnderflowADC,
+ const ZDCModuleFloatArray& LGOverflowADC)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetADCOverUnderflowValues(HGOverflowADC[side][module], HGUnderflowADC[side][module], LGOverflowADC[side][module]);
+ }
+ }
+}
+
+void ZDCDataAnalyzer::SetCutValues(const ZDCModuleFloatArray& chisqDivAmpCutHG, const ZDCModuleFloatArray& chisqDivAmpCutLG,
+ const ZDCModuleFloatArray& deltaT0MinHG, const ZDCModuleFloatArray& deltaT0MaxHG,
+ const ZDCModuleFloatArray& deltaT0MinLG, const ZDCModuleFloatArray& deltaT0MaxLG)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetCutValues(chisqDivAmpCutHG[side][module], chisqDivAmpCutLG[side][module],
+ deltaT0MinHG[side][module], deltaT0MaxHG[side][module],
+ deltaT0MinLG[side][module], deltaT0MaxLG[side][module]);
+ }
+ }
+}
+
+void ZDCDataAnalyzer::SetTimingCorrParams(const std::array<std::array<std::vector<float>, 4>, 2>& HGParamArr,
+ const std::array<std::array<std::vector<float>, 4>, 2>& LGParamArr)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetTimingCorrParams(HGParamArr.at(side).at(module), LGParamArr.at(side).at(module));
+ }
+ }
+
+}
+
+void ZDCDataAnalyzer::SetNonlinCorrParams(const std::array<std::array<std::vector<float>, 4>, 2>& HGNonlinCorrParams)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetNonlinCorrParams(HGNonlinCorrParams.at(side).at(module));
+ }
+ }
+}
+
+void ZDCDataAnalyzer::StartEvent(int lumiBlock)
+{
+ if (_debugLevel > 0) {
+ std::cout << "Starting new event, event index = " << _eventCount << std::endl;
+ }
+
+ // See if we have to load up new calibrations
+ //
+ if (lumiBlock != _currentLB) {
+ if (_debugLevel > 0) {
+ std::cout << "Starting new luminosity block " << lumiBlock << std::endl;
+ }
+
+ if (_haveECalib) {
+ if (_debugLevel > 1) {
+ std::cout << "Loading energy calibrations for event " << _eventCount << ", lumi block " << lumiBlock << std::endl;
+ }
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ float splineLBMin = _LBDepEcalibSplines[side][module]->GetXmin();
+ float splineLBMax = _LBDepEcalibSplines[side][module]->GetXmax();
+
+ if (lumiBlock >= splineLBMin && lumiBlock <= splineLBMax) {
+ _currentECalibCoeff[side][module] = _LBDepEcalibSplines[side][module]->Eval(lumiBlock);
+ }
+ else if (lumiBlock < splineLBMin) {
+ _currentECalibCoeff[side][module] = _LBDepEcalibSplines[side][module]->Eval(splineLBMin);
+ }
+ else {
+ _currentECalibCoeff[side][module] = _LBDepEcalibSplines[side][module]->Eval(splineLBMax);
+ }
+ }
+ }
+ } // end of if (_haveEcalib) {
+
+ if (_haveT0Calib) {
+ if (_debugLevel > 1) {
+ std::cout << "Loading timing calibrations for event " << _eventCount << ", lumi block " << lumiBlock << std::endl;
+ }
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ float splineLBMin = _T0HGOffsetSplines[side][module]->GetXmin();
+ float splineLBMax = _T0HGOffsetSplines[side][module]->GetXmax();
+
+ if (lumiBlock >= splineLBMin && lumiBlock <= splineLBMax) {
+ _currentT0OffsetsHG[side][module] = _T0HGOffsetSplines[side][module]->Eval(lumiBlock);
+ _currentT0OffsetsLG[side][module] = _T0LGOffsetSplines[side][module]->Eval(lumiBlock);
+ }
+ else if (lumiBlock < splineLBMin) {
+ _currentT0OffsetsHG[side][module] = _T0HGOffsetSplines[side][module]->Eval(splineLBMin);
+ _currentT0OffsetsLG[side][module] = _T0LGOffsetSplines[side][module]->Eval(splineLBMin);
+ }
+ else {
+ _currentT0OffsetsHG[side][module] = _T0HGOffsetSplines[side][module]->Eval(splineLBMax);
+ _currentT0OffsetsLG[side][module] = _T0LGOffsetSplines[side][module]->Eval(splineLBMax);
+ }
+ }
+ }
+ } // end of if (_haveT0Calib)
+ }
+
+ // Initialize transient results
+ //
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _dataLoaded[side][module] = false;
+ _moduleStatus[side][module] = 0;
+ _calibAmplitude[side][module] = 0;
+ _calibTime[side][module] = 0;
+ // _moduleFail[side][module] = false;
+ }
+
+ _moduleSum[side] = 0;
+ _moduleSumErrSq[side] = 0;
+ _moduleSumPreSample[side] = 0;
+
+ _calibModuleSum[side] = 0;
+ _calibModuleSumErrSq[side] = 0;
+
+ _averageTime[side] = 0;
+ _fail[side] = false;
+ }
+
+ _moduleMask = 0;
+ _currentLB = lumiBlock;
+}
+
+void ZDCDataAnalyzer::LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples)
+{
+ if (_debugLevel > 1) {
+ std::cout << "Loading data for event index " << _eventCount << ", side, module = " << side << ", " << module << std::endl;
+
+ if (_debugLevel > 2) {
+ std::cout << " Number of HG and LG samples = " << HGSamples.size() << ", " << LGSamples.size() << std::endl;
+ if (_debugLevel > 3) {
+ for (size_t sample = 0; sample < HGSamples.size(); sample++) {
+ std::cout << "HGSample[" << sample << "] = " << HGSamples[sample] << std::endl;
+ }
+ }
+ }
+ }
+
+ ZDCPulseAnalyzer* pulseAna_p = _moduleAnalyzers[side][module];
+
+ bool result = pulseAna_p->LoadAndAnalyzeData(HGSamples, LGSamples);
+ _dataLoaded[side][module] = true;
+
+ if (result) {
+ if (!pulseAna_p->BadT0() && !pulseAna_p->BadChisq()) {
+ int moduleMaskBit = 4*side + module;
+ _moduleMask |= 1<< moduleMaskBit;
+
+ float amplitude = pulseAna_p->GetAmplitude();
+ float ampError = pulseAna_p->GetAmpError();
+
+ _calibAmplitude[side][module] = amplitude*_currentECalibCoeff[side][module];
+
+ float calibAmpError = ampError * _currentECalibCoeff[side][module];
+
+ float timeCalib = pulseAna_p->GetT0Corr();
+ if (pulseAna_p->UseLowGain()) timeCalib -= _currentT0OffsetsLG[side][module];
+ else timeCalib -= _currentT0OffsetsHG[side][module];
+
+ _calibTime[side][module] = timeCalib;
+
+ _moduleSum[side] += amplitude;
+ _moduleSumErrSq[side] += ampError*ampError;
+
+ _moduleSumPreSample[side] += pulseAna_p->GetPreSampleAmp();
+
+ _calibModuleSum[side] += _calibAmplitude[side][module];
+ _calibModuleSumErrSq[side] += calibAmpError*calibAmpError;
+
+ _averageTime[side] += _calibTime[side][module]*_calibAmplitude[side][module];
+ }
+ }
+ else {
+ if (pulseAna_p->Failed()) {
+ if (_debugLevel >= 0) {
+ std::cout << "ZDCPulseAnalyzer::LoadData() returned fail for event " << _eventCount << ", side, module = " << side << ", " << module << std::endl;
+ }
+
+ _fail[side] = true;
+ }
+ }
+
+ _moduleStatus[side][module] = pulseAna_p->GetStatusMask();
+}
+
+bool ZDCDataAnalyzer::FinishEvent()
+{
+ // First make sure that all data is loaded
+ //
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ if (!_dataLoaded[side][module]) return false;
+ }
+
+ // Divide the average times by the calibrated module sums
+ //
+ if (_calibModuleSum[side] > 0) {
+ _averageTime[side] /= _calibModuleSum[side];
+ }
+ else {
+ _averageTime[side] = 0;
+ }
+ }
+
+ _eventCount++;
+ return true;
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
new file mode 100644
index 000000000000..0d146a734436
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
@@ -0,0 +1,150 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ZdcAnalysis/ZDCFitWrapper.h"
+
+ZDCFitExpFermiVariableTaus::ZDCFitExpFermiVariableTaus(std::string tag, float tmin, float tmax, bool fixTau1, bool fixTau2, float tau1, float tau2) :
+ ZDCFitWrapper(new TF1(("ExpFermiVariableTaus" + tag).c_str(), ZDCFermiExpFit, tmin, tmax, 5)),
+ _fixTau1(fixTau1), _fixTau2(fixTau2), _tau1(tau1), _tau2(tau2)
+{
+ TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+
+ theTF1->SetParName(0, "Amp");
+ theTF1->SetParName(1, "T0");
+ theTF1->SetParName(2, "#tau_{1}");
+ theTF1->SetParName(3, "#tau_{2}");
+ theTF1->SetParName(4, "s_{b}");
+
+ theTF1->SetParLimits(0, 0, 1024);
+ theTF1->SetParLimits(1, tmin, tmax);
+ theTF1->SetParLimits(2, 2.5, 6);
+ theTF1->SetParLimits(3, 10, 60);
+ theTF1->SetParLimits(4, -0.005, 10);
+
+ if (_fixTau1) theTF1->FixParameter(2, _tau1);
+ if (_fixTau2) theTF1->FixParameter(3, _tau2);
+}
+
+void ZDCFitExpFermiVariableTaus::Initialize(float initialAmp, float initialT0)
+{
+ TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+
+ theTF1->SetParameter(0, initialAmp);
+ theTF1->SetParameter(1, initialT0);
+
+ if (!_fixTau1) theTF1->SetParameter(2, _tau1);
+ if (!_fixTau2) theTF1->SetParameter(3, _tau2);
+
+ theTF1->SetParameter(4, 0);
+}
+
+ZDCFitExpFermiFixedTaus::ZDCFitExpFermiFixedTaus(std::string tag, float tmin, float tmax, float tau1, float tau2) :
+ ZDCFitWrapper(new TF1(("ExpFermiFixedTaus" + tag).c_str(), this, tmin, tmax, 3)),
+ _tau1(tau1), _tau2(tau2)
+{
+ TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+
+ theTF1->SetParLimits(0, 0, 1024);
+ theTF1->SetParLimits(1, tmin, tmax);
+ theTF1->SetParLimits(2, -0.005, 10);
+
+ theTF1->SetParName(0, "Amp");
+ theTF1->SetParName(1, "T0");
+ theTF1->SetParName(2, "s_{b}");
+
+ // Now create the reference function that we use to evaluate ExpFermiFit more efficiently
+ //
+ std::string funcNameRefFunc = "ExpFermiFixedTausRefFunc" + tag;
+
+ _expFermiFunc = new TF1(funcNameRefFunc.c_str(), ZDCFermiExpFit, -50, 100, 5);
+
+ _expFermiFunc->SetParameter(0, 1);
+ _expFermiFunc->SetParameter(1, 0);
+ _expFermiFunc->SetParameter(2, _tau1);
+ _expFermiFunc->SetParameter(3, _tau2);
+ _expFermiFunc->SetParameter(4, 0);
+
+ _norm = 1./_expFermiFunc->GetMaximum();
+ _timeCorr = 0.5*_tau1*std::log(_tau2/_tau1 - 1.0);
+}
+
+void ZDCFitExpFermiFixedTaus::Initialize(float initialAmp, float initialT0)
+{
+ GetWrapperTF1()->SetParameter(0, initialAmp);
+ GetWrapperTF1()->SetParameter(1, initialT0);
+ GetWrapperTF1()->SetParameter(2, 0);
+}
+
+ZDCFitExpFermiPrePulse::ZDCFitExpFermiPrePulse(std::string tag, float tmin, float tmax, float tau1, float tau2) :
+ ZDCPrePulseFitWrapper(new TF1(("ExpFermiPrePulse" + tag).c_str(), this, tmin, tmax, 5)),
+ _tau1(tau1), _tau2(tau2)
+{
+ // Create the reference function that we use to evaluate ExpFermiFit more efficiently
+ //
+ std::string funcNameRefFunc = "ExpFermiPerPulseRefFunc" + tag;
+
+ _expFermiFunc = new TF1(funcNameRefFunc.c_str(), ZDCFermiExpFit, -50, 100, 5);
+
+ _expFermiFunc->SetParameter(0, 1);
+ _expFermiFunc->SetParameter(1, 0);
+ _expFermiFunc->SetParameter(2, _tau1);
+ _expFermiFunc->SetParameter(3, _tau2);
+ _expFermiFunc->SetParameter(4, 0);
+
+ _norm = 1./_expFermiFunc->GetMaximum();
+ _timeCorr = 0.5*_tau1*std::log(_tau2/_tau1 - 1.0);
+
+ // Now set up the actual TF1
+ //
+ TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+
+ theTF1->SetParLimits(0, 0, 1024);
+ theTF1->SetParLimits(1, tmin, tmax);
+ theTF1->SetParLimits(2, 0, 1024);
+ theTF1->SetParLimits(3, -20, 10);
+ theTF1->SetParLimits(4, -0.005, 10);
+
+ theTF1->SetParName(0, "Amp");
+ theTF1->SetParName(1, "T0");
+ theTF1->SetParName(2, "Amp_{pre}");
+ theTF1->SetParName(3, "T0_{pre}");
+ theTF1->SetParName(4, "s_{b}");
+
+}
+
+void ZDCFitExpFermiPrePulse::Initialize(float initialAmp, float initialT0)
+{
+ GetWrapperTF1()->SetParameter(0, initialAmp);
+ GetWrapperTF1()->SetParameter(1, initialT0);
+ GetWrapperTF1()->SetParameter(2, 5);
+ GetWrapperTF1()->SetParameter(3, 0);
+ GetWrapperTF1()->SetParameter(4, 0);
+}
+
+double ZDCFermiExpFit(double* xvec, double* pvec)
+{
+ static const float offsetScale = 0;
+
+ double t = xvec[0];
+
+ double amp = pvec[0];
+ double t0 = pvec[1];
+ double tau1 = pvec[2];
+ double tau2 = pvec[3];
+ double linSlope = pvec[4];
+
+ double tauRatio = tau2/tau1;
+ double tauRatioMinunsOne = tauRatio - 1;
+
+ double norm = ( std::exp(-offsetScale/tauRatio)*pow(1./tauRatioMinunsOne, 1./(1 + tauRatio))/
+ ( 1 + pow(1./tauRatioMinunsOne, 1./(1 + 1/tauRatio))) );
+
+ double deltaT = t - (t0 - offsetScale*tau1);
+ if (deltaT < 0) deltaT = 0;
+
+ double expTerm = std::exp(-deltaT/tau2);
+ double fermiTerm = 1./(1. + std::exp(-(t - t0)/tau1));
+
+ return amp*expTerm*fermiTerm/norm + linSlope*t;
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
new file mode 100644
index 000000000000..141c2859f98a
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
@@ -0,0 +1,506 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ZdcAnalysis/ZDCPulseAnalyzer.h"
+
+#include <numeric>
+#include <algorithm>
+
+#include "TFitResult.h"
+#include "TFitResultPtr.h"
+
+bool ZDCPulseAnalyzer::_quietFits = true;
+std::string ZDCPulseAnalyzer::_fitOptions = "";
+
+ZDCPulseAnalyzer::ZDCPulseAnalyzer(std::string tag, int Nsample, float deltaTSample, size_t preSampleIdx, int pedestal,
+ float gainHG, std::string fitFunction, int peak2ndDerivMinSample,
+ float peak2ndDerivMinThreshHG, float peak2ndDerivMinThreshLG) :
+ _tag(tag), _Nsample(Nsample), _deltaTSample(deltaTSample), _preSampleIdx(preSampleIdx),
+ _pedestal(pedestal), _gainHG(gainHG), _forceLG(false), _fitFunction(fitFunction),
+ _peak2ndDerivMinSample(peak2ndDerivMinSample),
+ _peak2ndDerivMinThreshHG(peak2ndDerivMinThreshHG),
+ _peak2ndDerivMinThreshLG(peak2ndDerivMinThreshLG),
+ _haveTimingCorrections(false), _haveNonlinCorr(false), _initializedFits(false),
+ _ADCSamplesHGSub(Nsample, 0), _ADCSamplesLGSub(Nsample, 0),
+ _ADCSSampSigHG(Nsample, 1), _ADCSSampSigLG(Nsample, 1), // By default the ADC uncertainties are set to one
+ _samplesSub(Nsample, 0),
+ _defaultFitWrapper(0), _prePulseFitWrapper(0)
+{
+ // Create the histogram used for fitting
+ //
+ _tmin = -deltaTSample/2;
+ _tmax = _tmin + ((float) Nsample)*deltaTSample;
+
+ std::string histName = "ZDCFitHist" + tag;
+
+ _fitHist = new TH1F(histName.c_str(), "", _Nsample, _tmin, _tmax);
+ _fitHist->SetDirectory(0);
+
+ SetDefaults();
+ Reset();
+}
+
+ZDCPulseAnalyzer::~ZDCPulseAnalyzer()
+{
+ if (!_defaultFitWrapper) delete _defaultFitWrapper;
+ if (!_prePulseFitWrapper) delete _prePulseFitWrapper;
+}
+
+void ZDCPulseAnalyzer::SetDefaults()
+{
+ _nominalTau1 = 4;
+ _nominalTau2 = 21;
+
+ _fixTau1 = false;
+ _fixTau2 = false;
+
+ _HGOverflowADC = 900;
+ _HGUnderflowADC = 20;
+ _LGOverflowADC = 1000;
+
+ _chisqDivAmpCutLG = 100;
+ _chisqDivAmpCutHG = 100;
+
+ _T0CutLowLG = _tmin;
+ _T0CutLowHG = _tmin;
+
+ _T0CutHighLG = _tmax;
+ _T0CutHighHG = _tmax;
+
+ _LGT0CorrParams.assign(4, 0);
+ _HGT0CorrParams.assign(4, 0);
+
+ _fitOptions = "s";
+}
+
+void ZDCPulseAnalyzer::Reset()
+{
+ _haveData = false;
+ _havePulse = false;
+ _fail = false;
+
+ _HGOverflow = false;
+ _HGUnderflow = false;
+
+ _LGOverflow = false;
+ _LGUnderflow = false;
+
+ _PSHGOverUnderflow = false;
+
+ _prePulse = false;
+ _fitFailed = false;
+
+ _useLowGain = false;
+
+ _badT0 = false;
+ _badChisq = false;
+
+ _fitAmplitude = 0;
+ _fitT0 = -100;
+ _fitChisq = 0;
+
+ _fitT0Corr = -100;
+ _amplitude = 0;
+ _ampError = 0;
+ _preSampleAmp = 0;
+ _bkgdMaxFraction = 0;
+
+ _samplesSub.clear();
+ _samplesDeriv.clear();
+ _samplesDeriv2nd.clear();
+}
+
+void ZDCPulseAnalyzer::SetTauT0Values(bool fixTau1, bool fixTau2, float tau1, float tau2, float t0HG, float t0LG)
+{
+ _fixTau1 = fixTau1;
+ _fixTau2 = fixTau2;
+ _nominalTau1 = tau1;
+ _nominalTau2 = tau2;
+
+ _nominalT0HG = t0HG;
+ _nominalT0LG = t0LG;
+
+ _initializedFits = false;
+}
+
+void ZDCPulseAnalyzer::SetADCOverUnderflowValues(int HGOverflowADC, int HGUnderflowADC, int LGOverflowADC)
+{
+ _HGOverflowADC = HGOverflowADC;
+ _LGOverflowADC = LGOverflowADC;
+ _HGUnderflowADC = HGUnderflowADC;
+}
+
+void ZDCPulseAnalyzer::SetCutValues(float chisqDivAmpCutHG, float chisqDivAmpCutLG,
+ float deltaT0MinHG, float deltaT0MaxHG,
+ float deltaT0MinLG, float deltaT0MaxLG)
+{
+ _chisqDivAmpCutHG = chisqDivAmpCutHG;
+ _chisqDivAmpCutLG = chisqDivAmpCutLG;
+
+ _T0CutLowHG = deltaT0MinHG;
+ _T0CutLowLG = deltaT0MinLG;
+
+ _T0CutHighHG = deltaT0MaxHG;
+ _T0CutHighLG = deltaT0MaxLG;
+}
+
+void ZDCPulseAnalyzer::SetupFitFunctions()
+{
+ if (!_defaultFitWrapper) delete _defaultFitWrapper;
+ if (!_prePulseFitWrapper) delete _prePulseFitWrapper;
+
+ if (_fitFunction == "FermiExp") {
+ if (!_fixTau1 || !_fixTau2) {
+ //
+ // Use the variable tau version of the expFermiFit
+ //
+ _defaultFitWrapper = new ZDCFitExpFermiVariableTaus(_tag, _tmin, _tmax, _fixTau1, _fixTau2, _nominalTau1, _nominalTau2);
+ }
+ else {
+ _defaultFitWrapper = new ZDCFitExpFermiFixedTaus(_tag, _tmin, _tmax, _nominalTau1, _nominalTau2);
+ }
+
+ _prePulseFitWrapper = new ZDCFitExpFermiPrePulse(_tag, _tmin, _tmax, _nominalTau1, _nominalTau2);
+ }
+ else {
+ throw;
+ }
+
+ _initializedFits = true;
+}
+
+bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG)
+{
+ if (!_initializedFits) SetupFitFunctions();
+
+ // Clear any transient data
+ //
+ Reset();
+
+ // Make sure we have the right number of samples. Should never fail. necessry?
+ //
+ if (ADCSamplesHG.size() != _Nsample || ADCSamplesLG.size() != _Nsample) {
+ _fail = true;
+ return false;
+ }
+
+ // Now do pedestal subtraction and check for overflows
+ //
+ for (size_t isample = 0; isample < _Nsample; isample++) {
+ float ADCHG = ADCSamplesHG[isample];
+ float ADCLG = ADCSamplesLG[isample];
+
+ if (ADCHG > _HGOverflowADC) {
+ _HGOverflow = true;
+
+ if (isample == _preSampleIdx) _PSHGOverUnderflow = true;
+ }
+ else if (ADCHG < _HGUnderflowADC) {
+ _HGUnderflow = true;
+
+ if (isample == _preSampleIdx) _PSHGOverUnderflow = true;
+ }
+
+ if (ADCLG > _LGOverflowADC) {
+ _LGOverflow = true;
+ _fail = true;
+ _amplitude = 1024 * _gainHG; // Give a vale here even though we know it's wrong because
+ // the user may not check the return value and we know that
+ // amplitude is bigger than this
+ }
+
+ if (ADCLG == 0) {
+ _LGUnderflow = true;
+ _fail = true;
+ }
+
+ _ADCSamplesHGSub[isample] = ADCHG - _pedestal;
+ _ADCSamplesLGSub[isample] = ADCLG - _pedestal;
+ }
+
+ // if (_fail) return false;
+
+ // Now we actually do the analysis using the high or low gain samples as determined
+ // on the basis of the above analysis
+ //
+ _useLowGain = _HGUnderflow || _HGOverflow || _forceLG;
+
+ if (_useLowGain) {
+ bool result = AnalyzeData(_ADCSamplesLGSub, _ADCSSampSigLG, _LGT0CorrParams, _chisqDivAmpCutLG,
+ _T0CutLowLG, _T0CutHighLG, _peak2ndDerivMinThreshLG);
+ if (result) {
+ //
+ // Multiply amplitude by gain factor
+ //
+ _amplitude = _fitAmplitude * _gainHG;
+ _ampError = _fitAmpError * _gainHG;
+ _preSampleAmp = _preSample * _gainHG;
+ }
+
+ return result;
+ }
+ else {
+ bool result = AnalyzeData(_ADCSamplesHGSub, _ADCSSampSigHG, _HGT0CorrParams, _chisqDivAmpCutHG,
+ _T0CutLowHG, _T0CutHighHG, _peak2ndDerivMinThreshHG);
+ if (result) {
+ _preSampleAmp = _preSample;
+ _amplitude = _fitAmplitude;
+ _ampError = _fitAmpError;
+
+ // If we have a non-linear correction, apply it here
+ //
+ if (_haveNonlinCorr) {
+ float ampCorrFact = _amplitude/1000. - 0.5;
+ float corrFactor = 1./(1. + _nonLinCorrParams[0]*ampCorrFact + _nonLinCorrParams[1]*ampCorrFact*ampCorrFact);
+
+ _amplitude *= corrFactor;
+ _ampError *= corrFactor;
+ }
+ }
+
+ return result;
+ }
+}
+
+bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, // The samples used for this event
+ const std::vector<float>& samplesSig, // The "resolution" on the ADC value
+ const std::vector<float>& t0CorrParams, // The parameters used to correct the t0
+ float maxChisqDivAmp, // The maximum chisq / amplitude ratio
+ float minT0Corr, float maxT0Corr, // The minimum and maximum corrected T0 values
+ float peak2ndDerivMinThresh
+ )
+{
+ // Do the presample subtraction using a lambda function
+ //
+ float preSample = samples[_preSampleIdx];
+ _preSample = preSample;
+ _samplesSub = samples;
+
+ std::for_each(_samplesSub.begin(), _samplesSub.end(), [=] (float& adcUnsub) {return adcUnsub -= preSample;} );
+
+ // Find maximum and minimum values
+ //
+ std::pair<SampleCIter, SampleCIter> minMaxIters = std::minmax_element(_samplesSub.begin(), _samplesSub.end());
+ SampleCIter minIter = minMaxIters.first;
+ SampleCIter maxIter = minMaxIters.second;
+
+ _maxADCValue = *maxIter;
+ _minADCValue = *minIter;
+
+ float maxMinADCDiff = _maxADCValue - _minADCValue;
+
+ _maxSampl = std::distance(_samplesSub.cbegin(), maxIter);
+ _minSampl = std::distance(_samplesSub.cbegin(), minIter);
+
+ // Calculate first and 2nd "derivatives"
+ //
+ std::vector<float> deriv(_Nsample, 0);
+ std::adjacent_difference(samples.begin(), samples.end(), deriv.begin());
+
+ // Save the derivatives, dropping the useless first element
+ //
+ _samplesDeriv.insert(_samplesDeriv.begin(), deriv.begin() + 1, deriv.end());
+
+ std::vector<float> deriv2nd(_Nsample, 0);
+ std::adjacent_difference(_samplesDeriv.begin(), _samplesDeriv.end(), deriv2nd.begin());
+
+ // Save the second derivatives, dropping the useless first element
+ //
+ _samplesDeriv2nd.insert(_samplesDeriv2nd.begin(), deriv2nd.begin() + 1, deriv2nd.end());
+
+ // Find the sampl which has the lowest 2nd derivative
+ //
+ SampleCIter minDeriv2ndIter = std::min_element(_samplesDeriv2nd.begin(), _samplesDeriv2nd.end());
+
+ _minDeriv2nd = *minDeriv2ndIter;
+ _minDeriv2ndIndex = std::distance(_samplesDeriv2nd.cbegin(), minDeriv2ndIter);
+
+ // The ket test for presence of a pulse: is the minimum 2nd derivative in the right place and is it
+ // large enough (sufficiently negative)
+ //
+ if (std::abs(_minDeriv2ndIndex - _peak2ndDerivMinSample) <= 1 && _minDeriv2nd <= peak2ndDerivMinThresh) {
+ _havePulse = true;
+ }
+ else {
+ _havePulse = false;
+ }
+
+ // Now decide whether we have a preceeding pulse or not
+ //
+
+ // Have to be careful with the derivative indices, but test up to minDeriv2ndIndex - 1
+ //
+ bool haveNegDeriv = false;
+
+ for (int isampl = 0; isampl < _minDeriv2ndIndex; isampl++) {
+
+ // If any of the derivatives prior to the peak (as determined by the 2nd derivative) is negative
+ // Then we are on the tailing edge of a previous pulse.
+ //
+ if (_samplesDeriv[isampl] < 0) {
+
+ //
+ // Put a threshold on the derivative as a fraction of the difference between max and min so that
+ // we don't waste lots of time for insignificant prepulses.
+ //
+ if (std::abs(_samplesDeriv[isampl]) / (_maxADCValue - _minADCValue + 1) > 0.05) _prePulse = true;
+
+ // Also, if we have two negative derivatives in a row, we have a prepulse
+ //
+ if (haveNegDeriv) _prePulse = true;
+ else haveNegDeriv = true;
+ }
+
+ // If any of the 2nd derivatives before the minium are negative, we have the possibility of positive
+ // preceeding pulse
+ //
+ if (_samplesDeriv2nd[isampl] < 0 && _samplesDeriv2nd[isampl + 1] > 0) _prePulse = true;
+ }
+
+ // Additional specific tests for pre pulses etc.
+ //
+ if (_deltaTSample < 20) {
+ //
+ // Don't do this check for the 40 MHz data as the pulse often drops below baseline
+ // due to overshoot, not pre-pulses
+ //
+ if (_samplesSub[_Nsample - 1] <= 0 ||
+ _samplesSub[_Nsample - 1] <= _samplesSub[_preSampleIdx+1] ) _prePulse = true;
+
+ }
+
+
+ FillHistogram(_samplesSub, samplesSig);
+
+ // Stop now if we have no pulse
+ //
+ if (_fail || !_havePulse) return false;
+
+ DoFit();
+
+ if (!FitFailed()) {
+ _fitT0Corr = _fitT0;
+
+ if (_haveTimingCorrections) {
+
+ //
+ // We correct relative to the middle of the amplitude range, divided by 100 to make corrections smaller
+ //
+ float t0CorrFact = (_fitAmplitude - 500)/100.;
+
+ float correction = (t0CorrParams[0] + t0CorrParams[1]*t0CorrFact + t0CorrParams[2]*t0CorrFact*t0CorrFact +
+ t0CorrParams[3]*t0CorrFact*t0CorrFact*t0CorrFact);
+
+ _fitT0Corr -= correction;
+ }
+
+
+ // Now check for valid chisq and valid time
+ //
+ if (_fitChisq / _fitAmplitude > maxChisqDivAmp) _badChisq = true;
+ if (_fitT0Corr < minT0Corr || _fitT0Corr > maxT0Corr) _badT0 = true;
+
+ // std::cout << "Testing t0 cut: t0corr = " << _fitT0Corr << ", min max cuts = " << minT0Corr << ", " << maxT0Corr << std::endl;
+ }
+
+ return !_fitFailed;
+}
+
+void ZDCPulseAnalyzer::DoFit()
+{
+ // Set the initial values
+ //
+ float ampInitial = _maxADCValue - _minADCValue;
+ float t0Initial = (_useLowGain ? _nominalT0LG : _nominalT0HG);
+
+ ZDCFitWrapper* fitWrapper = _defaultFitWrapper;
+ if (PrePulse()) fitWrapper = _prePulseFitWrapper;
+
+ fitWrapper->Initialize(ampInitial, t0Initial);
+ if (PrePulse() && _samplesDeriv[0] < 0) {
+ (static_cast<ZDCPrePulseFitWrapper*>(_prePulseFitWrapper))->SetInitialPrePulseT0(-10);
+ }
+
+ // Now perform the fit
+ //
+ std::string options = _fitOptions;
+ if (QuietFits()) options += "Q0";
+ options += "s";
+
+ TFitResultPtr result_ptr = _fitHist->Fit(fitWrapper->GetWrapperTF1(), options.c_str());
+
+ int fitStatus = result_ptr;
+ if (fitStatus != 0) {
+ TFitResultPtr try2Result_ptr = _fitHist->Fit(fitWrapper->GetWrapperTF1(), options.c_str());
+ if ((int) try2Result_ptr != 0) _fitFailed = true;
+ }
+ else _fitFailed = false;
+
+ // if (!_fitFailed) {
+ _bkgdMaxFraction = fitWrapper->GetBkgdMaxFraction();
+ if (std::abs(_bkgdMaxFraction) > 0.25) {
+ options += "e";
+ _fitHist->Fit(fitWrapper->GetWrapperTF1(), options.c_str());
+ _bkgdMaxFraction = fitWrapper->GetBkgdMaxFraction();
+ }
+
+ _fitAmplitude = fitWrapper->GetAmplitude();
+ _fitT0 = fitWrapper->GetTime() - t0Initial;
+ _fitChisq = result_ptr->Chi2();
+
+ _fitTau1 = fitWrapper->GetTau1();
+ _fitTau2 = fitWrapper->GetTau2();
+
+ _fitAmpError = fitWrapper->GetAmpError();
+ // }
+}
+
+void ZDCPulseAnalyzer::Dump() const
+{
+ std::cout << "ZDCPulseAnalyzer dump for tag = " << _tag << std::endl;
+
+ std::cout << "Presample index, value = " << _preSampleIdx << ", " << _preSample << std::endl;
+
+ std::cout << "samplesSub ";
+ for (size_t sample = 0; sample < _samplesSub.size(); sample++) {
+ std::cout << ", [" << sample << "] = " << _samplesSub[sample];
+ }
+ std::cout << std::endl;
+
+ std::cout << "samplesDeriv ";
+ for (size_t sample = 0; sample < _samplesDeriv.size(); sample++) {
+ std::cout << ", [" << sample << "] = " << _samplesDeriv[sample];
+ }
+ std::cout << std::endl;
+
+ std::cout << "samplesDeriv2nd ";
+ for (size_t sample = 0; sample < _samplesDeriv2nd.size(); sample++) {
+ std::cout << ", [" << sample << "] = " << _samplesDeriv2nd[sample];
+ }
+ std::cout << std::endl;
+
+ std::cout << "minimum 2nd deriv sample, value = " << _minDeriv2ndIndex << ", " << _minDeriv2nd << std::endl;
+}
+
+unsigned int ZDCPulseAnalyzer::GetStatusMask() const
+{
+ unsigned int statusMask = 0;
+
+ if (HavePulse()) statusMask |= 1 << PulseBit;
+ if (UseLowGain()) statusMask |= 1 << LowGainBit;
+ if (Failed()) statusMask |= 1 << FailBit;
+ if (HGOverflow()) statusMask |= 1 << HGOverflowBit;
+
+ if (HGUnderflow()) statusMask |= 1 << HGUnderflowBit;
+ if (PSHGOverUnderflow()) statusMask |= 1 << PSHGOverUnderflowBit;
+ if (LGOverflow()) statusMask |= 1 << LGOverflowBit;
+ if (LGUnderflow()) statusMask |= 1 << LGUnderflowBit;
+
+ if (FitFailed()) statusMask |= 1 << FitFailedBit;
+ if (PrePulse()) statusMask |= 1 << PrePulseBit;
+ if (BadChisq()) statusMask |= 1 << BadChisqBit;
+ if (BadT0()) statusMask |= 1 << BadT0Bit;
+
+ return statusMask;
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCTriggerEfficiency.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCTriggerEfficiency.cxx
new file mode 100644
index 000000000000..d1f702893b66
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCTriggerEfficiency.cxx
@@ -0,0 +1,63 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ZdcAnalysis/ZDCTriggerEfficiency.h"
+#include <vector>
+
+using namespace std;
+float ZDCTriggerEfficiency::GetEfficiency(int side, float ADCSum){
+ double alpha = _currentParams[side][0];
+ double beta = _currentParams[side][1];
+ double theta = _currentParams[side][2];
+ if(alpha<10e-5||beta<10e-5||theta<10e-5)
+ {
+ return -1;
+ }
+ float m = exp(-ADCSum/theta);
+ float p = exp(-(ADCSum - alpha) / beta);
+ float efficiency = (1 - m) / (1 + p);
+ return efficiency;
+}
+
+std::pair<float, float> ZDCTriggerEfficiency::GetEfficiencyAndError(int side, float ADCSum){
+ double alpha = _currentParams[side][0];
+ double beta = _currentParams[side][1];
+ double theta = _currentParams[side][2];
+ if(alpha<10e-5||beta<10e-5||theta<10e-5)
+ {
+ return std::make_pair(-1,-1);
+ }
+
+ double alphaErr = _currentParamErrors[side][0];
+ double betaErr = _currentParamErrors[side][1];
+ double thetaErr = _currentParamErrors[side][2];
+
+ double cov_alpha_beta = _currentCorrCoefff[side][0];
+ double cov_alpha_theta = _currentCorrCoefff[side][1];
+ double cov_beta_theta = _currentCorrCoefff[side][2];
+
+ float m = exp(-ADCSum/theta);
+ float n = exp(alpha/beta);
+ float p = exp(-(ADCSum - alpha) / beta);
+ float efficiency = (1 - m) / (1 + p);
+ float dda = (1 - m) * -n / beta / (1 + p) / (1 + p);
+ float ddb = (1 - m) * -p / pow(1 + p, 2.0) * ADCSum / beta / beta;
+ float ddt = -ADCSum * m / theta / theta / (1 + p);
+
+ float efficiencyErr = sqrt(
+ alphaErr * alphaErr * dda * dda +
+ betaErr * betaErr * ddb * ddb +
+ thetaErr * thetaErr * ddt * ddt +
+ 2 * cov_alpha_beta * dda * ddb +
+ 2 * cov_alpha_theta * dda * ddt +
+ 2 * cov_beta_theta * ddb * ddt);
+ if(alphaErr * alphaErr * dda * dda +
+ betaErr * betaErr * ddb * ddb +
+ thetaErr * thetaErr * ddt * ddt +
+ 2 * cov_alpha_beta * dda * ddb +
+ 2 * cov_alpha_theta * dda * ddt +
+ 2 * cov_beta_theta * ddb * ddt<0) efficiencyErr = 0;
+
+ return std::make_pair(efficiency, efficiencyErr);
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
new file mode 100644
index 000000000000..535aa78f78b8
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
@@ -0,0 +1,947 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ZdcAnalysis/ZdcAnalysisTool.h"
+#include "TGraph.h"
+#include "TEnv.h"
+#include "TSystem.h"
+#include "ZdcAnalysis/ZdcSincInterp.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "TFile.h"
+#include <sstream>
+#include <memory>
+#include <xAODForward/ZdcModuleAuxContainer.h>
+#include "PathResolver/PathResolver.h"
+
+namespace ZDC
+{
+
+ ZdcAnalysisTool::ZdcAnalysisTool(const std::string& name) : asg::AsgTool(name), m_name(name), m_init(false),
+ m_writeAux(false), m_eventReady(false),
+ m_runNumber(0), m_lumiBlock(0),
+ m_splines{0,0,0,0,0,0,0,0}, m_zdcTriggerEfficiency(0)
+{
+
+#ifndef XAOD_STANDALONE
+ declareInterface<IZdcAnalysisTool>(this);
+#endif
+
+ declareProperty("ZdcModuleContainerName",m_zdcModuleContainerName="ZdcModules","Location of ZDC processed data");
+ declareProperty("Configuration", m_configuration = "PbPb2015");
+ declareProperty("FlipEMDelay",m_flipEMDelay=false);
+ declareProperty("LowGainOnly",m_lowGainOnly=false);
+ declareProperty("WriteAux",m_writeAux=true);
+ declareProperty("AuxSuffix",m_auxSuffix="");
+
+ // The following job properties enable/disable and affect the calibration of the ZDC energies
+ //
+ declareProperty("DoCalib",m_doCalib=true);
+ declareProperty("ZdcAnalysisConfigPath",m_zdcAnalysisConfigPath="$ROOTCOREBIN/data/ZdcAnalysis","ZDC Analysis config file path");
+ //declareProperty("ForceCalibRun",m_forceCalibRun=287931); // last run of Pb+Pb 2015
+ declareProperty("ForceCalibRun",m_forceCalibRun=-1); // last run of Pb+Pb 2015
+ declareProperty("ForceCalibLB",m_forceCalibLB=814); // last LB of Pb+Pb 2015
+
+ // The following parameters are primarily used for the "default" configuration, but also may be
+ // use to modify/tailor other configurations
+ //
+ declareProperty("NumSampl", m_numSample = 7);
+ declareProperty("DeltaTSample", m_deltaTSample = 25);
+ declareProperty("Presample", m_presample = 0);
+ declareProperty("PeakSample", m_peakSample = 2);
+ declareProperty("Peak2ndDerivThresh", m_Peak2ndDerivThresh = 10);
+
+ declareProperty("T0", m_t0 = 50);
+ declareProperty("Tau1", m_tau1 = 5);
+ declareProperty("Tau2", m_tau2 = 25);
+ declareProperty("FixTau1", m_fixTau1 = false);
+ declareProperty("FixTau2", m_fixTau2 = false);
+
+ declareProperty("DeltaTCut", m_deltaTCut = 25);
+ declareProperty("ChisqRatioCut", m_ChisqRatioCut = 10);
+
+ //ATH_MSG_INFO("Creating ZdcAnalysisoTool named " << m_name);
+ //ATH_MSG_INFO("ZDC config file path " << m_zdcAnalysisConfigPath);
+
+}
+
+ZdcAnalysisTool::~ZdcAnalysisTool()
+{
+ ATH_MSG_DEBUG("Deleting ZdcAnalysisTool named " << m_name);
+ //SafeDelete( m_zdcDataAnalyzer );
+}
+
+void ZdcAnalysisTool::initializeTriggerEffs(unsigned int runNumber)
+{
+ if (!m_zdcTriggerEfficiency)
+ {
+ m_zdcTriggerEfficiency = new ZDCTriggerEfficiency();
+ }
+
+ std::string filename = PathResolverFindCalibFile( ("ZdcAnalysis/"+m_zdcTriggerEffParamsFileName) );
+ //std::cout << "Found trigger config file " << filename << std::endl;
+
+ TFile* file = TFile::Open(filename.c_str());
+ if (!file->IsOpen())
+ {
+ ATH_MSG_ERROR("No trigger efficiencies at " << filename);
+ return;
+ }
+
+ //file->Print();
+
+ stringstream Aalpha_name;
+ Aalpha_name<<"A_alpha_"<<runNumber;
+ TSpline3* par_A_alpha = (TSpline3*)file->GetObjectChecked(Aalpha_name.str().c_str(),"TSpline3");
+
+ if (!par_A_alpha)
+ {
+ ATH_MSG_ERROR("No trigger efficiencies for run number " << runNumber);
+ m_doCalib = false;
+ }
+
+ stringstream Abeta_name;
+ Abeta_name<<"A_beta_"<<runNumber;
+ TSpline3* par_A_beta = (TSpline3*)file->GetObjectChecked(Abeta_name.str().c_str(),"TSpline3");
+ stringstream Atheta_name;
+ Atheta_name<<"A_theta_"<<runNumber;
+ TSpline3* par_A_theta = (TSpline3*)file->GetObjectChecked(Atheta_name.str().c_str(),"TSpline3");
+
+ stringstream Calpha_name;
+ Calpha_name<<"C_alpha_"<<runNumber;
+ TSpline3* par_C_alpha = (TSpline3*)file->GetObjectChecked(Calpha_name.str().c_str(),"TSpline3");
+ stringstream Cbeta_name;
+ Cbeta_name<<"C_beta_"<<runNumber;
+ TSpline3* par_C_beta = (TSpline3*)file->GetObjectChecked(Cbeta_name.str().c_str(),"TSpline3");
+ stringstream Ctheta_name;
+ Ctheta_name<<"C_theta_"<<runNumber;
+ TSpline3* par_C_theta = (TSpline3*)file->GetObjectChecked(Ctheta_name.str().c_str(),"TSpline3");
+
+ stringstream Err_Aalpha_name;
+ Err_Aalpha_name<<"A_alpha_error_"<<runNumber;
+ TSpline3* parErr_A_alpha = (TSpline3*)file->GetObjectChecked(Err_Aalpha_name.str().c_str(),"TSpline3");
+ stringstream Err_Abeta_name;
+ Err_Abeta_name<<"A_beta_error_"<<runNumber;
+ TSpline3* parErr_A_beta = (TSpline3*)file->GetObjectChecked(Err_Abeta_name.str().c_str(),"TSpline3");
+ stringstream Err_Atheta_name;
+ Err_Atheta_name<<"A_theta_error_"<<runNumber;
+ TSpline3* parErr_A_theta = (TSpline3*)file->GetObjectChecked(Err_Atheta_name.str().c_str(),"TSpline3");
+
+ stringstream Err_Calpha_name;
+ Err_Calpha_name<<"C_alpha_error_"<<runNumber;
+ TSpline3* parErr_C_alpha = (TSpline3*)file->GetObjectChecked(Err_Calpha_name.str().c_str(),"TSpline3");
+ stringstream Err_Cbeta_name;
+ Err_Cbeta_name<<"C_beta_error_"<<runNumber;
+ TSpline3* parErr_C_beta = (TSpline3*)file->GetObjectChecked(Err_Cbeta_name.str().c_str(),"TSpline3");
+ stringstream Err_Ctheta_name;
+ Err_Ctheta_name<<"C_theta_error_"<<runNumber;
+ TSpline3* parErr_C_theta = (TSpline3*)file->GetObjectChecked(Err_Ctheta_name.str().c_str(),"TSpline3");
+
+
+ stringstream Cov_A_alpha_beta_name;
+ Cov_A_alpha_beta_name<<"cov_A_alpha_beta_"<<runNumber;
+ TSpline3* cov_A_alpha_beta = (TSpline3*)file->GetObjectChecked(Cov_A_alpha_beta_name.str().c_str(),"TSpline3");
+ stringstream Cov_A_alpha_theta_name;
+ Cov_A_alpha_theta_name<<"cov_A_alpha_theta_"<<runNumber;
+ TSpline3* cov_A_alpha_theta = (TSpline3*)file->GetObjectChecked(Cov_A_alpha_theta_name.str().c_str(),"TSpline3");
+ stringstream Cov_A_beta_theta_name;
+ Cov_A_beta_theta_name<<"cov_A_beta_theta_"<<runNumber;
+ TSpline3* cov_A_beta_theta = (TSpline3*)file->GetObjectChecked(Cov_A_beta_theta_name.str().c_str(),"TSpline3");
+
+ stringstream Cov_C_alpha_beta_name;
+ Cov_C_alpha_beta_name<<"cov_C_alpha_beta_"<<runNumber;
+ TSpline3* cov_C_alpha_beta = (TSpline3*)file->GetObjectChecked(Cov_C_alpha_beta_name.str().c_str(),"TSpline3");
+ stringstream Cov_C_alpha_theta_name;
+ Cov_C_alpha_theta_name<<"cov_C_alpha_theta_"<<runNumber;
+ TSpline3* cov_C_alpha_theta = (TSpline3*)file->GetObjectChecked(Cov_C_alpha_theta_name.str().c_str(),"TSpline3");
+ stringstream Cov_C_beta_theta_name;
+ Cov_C_beta_theta_name<<"cov_C_beta_theta_"<<runNumber;
+ TSpline3* cov_C_beta_theta = (TSpline3*)file->GetObjectChecked(Cov_C_beta_theta_name.str().c_str(),"TSpline3");
+
+ std::array<std::vector<TSpline3*>, 2> effparams;
+ std::array<std::vector<TSpline3*>, 2> effparamErrors;
+ std::array<std::vector<TSpline3*>, 2> effparamsCorrCoeffs;
+ //side0: C; side1: A
+ effparams[0] = {par_C_alpha, par_C_beta, par_C_theta};
+ effparams[1] = {par_A_alpha, par_A_beta, par_A_theta};
+ effparamErrors[0] = {parErr_C_alpha, parErr_C_beta, parErr_C_theta};
+ effparamErrors[1] = {parErr_A_alpha, parErr_A_beta, parErr_A_theta};
+ effparamsCorrCoeffs[0] = {cov_C_alpha_beta, cov_C_alpha_theta, cov_C_beta_theta};
+ effparamsCorrCoeffs[1] = {cov_A_alpha_beta, cov_A_alpha_theta, cov_A_beta_theta};
+
+ m_zdcTriggerEfficiency->SetEffParamsAndErrors(effparams, effparamErrors);
+ m_zdcTriggerEfficiency->SetEffParamCorrCoeffs(effparamsCorrCoeffs);
+
+ return;
+
+}
+
+ZDCDataAnalyzer* ZdcAnalysisTool::initializeDefault()
+{
+ // We rely completely on the default parameters specified in the job properties to control:
+ // # samples
+ // frequency (more precisely, time/sample)
+ // which sample to use as the pre-sample
+ // where to expact the maxim of the peak (min 2nd derivative)
+ // thresholds on the 2nd derivative for valid pulses
+ // whether to fix the tau values in the pulse fitting
+ // the default tau values
+ // the nominal T0
+ // delta T and chisq/amp cuts
+ //
+ // For now, we continue to use hard-coded values for the maximum and minimum ADC values
+ // For now we also use the FermiExp pulse model.
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau1, tau2, peak2ndDerivMinSamples, t0;
+ ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG;
+ ZDCDataAnalyzer::ZDCModuleFloatArray deltaT0CutLow, deltaT0CutHigh, chisqDivAmpCut;
+ ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ fixTau1Arr[side][module] = m_fixTau1;
+ fixTau2Arr[side][module] = m_fixTau2;
+ tau1[side][module] = m_tau1;
+ tau2[side][module] = m_tau2;
+
+ peak2ndDerivMinSamples[side][module] = m_peakSample;
+ peak2ndDerivMinThresholdsHG[side][module] = -m_Peak2ndDerivThresh;
+ peak2ndDerivMinThresholdsLG[side][module] = -m_Peak2ndDerivThresh/2;
+
+ t0[side][module] = m_t0;
+ deltaT0CutLow[side][module] = -m_deltaTCut;
+ deltaT0CutHigh[side][module] = m_deltaTCut;
+ chisqDivAmpCut[side][module] = m_ChisqRatioCut;
+ }
+ }
+
+ ATH_MSG_INFO( "Default: delta t cut, value low = " << deltaT0CutLow[0][0] << ", high = " << deltaT0CutHigh[0][0] );
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {800, 800, 800, 800, 800, 800, 800, 800};
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {1020, 1020, 1020, 1020, 1020, 1020, 1020, 1020};
+
+ // Construct the data analyzer
+ //
+ ZDCDataAnalyzer* zdcDataAnalyzer = new ZDCDataAnalyzer(m_numSample, m_deltaTSample, m_presample, "FermiExp", peak2ndDerivMinSamples,
+ peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG, m_lowGainOnly);
+
+ zdcDataAnalyzer->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+ zdcDataAnalyzer->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0, t0);
+ zdcDataAnalyzer->SetCutValues(chisqDivAmpCut, chisqDivAmpCut, deltaT0CutLow, deltaT0CutHigh, deltaT0CutLow, deltaT0CutHigh);
+
+ return zdcDataAnalyzer;
+}
+
+
+void ZdcAnalysisTool::initialize40MHz()
+{
+ // We have a complete configuration and so we override all of the default parameters
+ //
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau1 = {4.2, 3.8, 5.2, 5.0,
+ 5.0, 3.7, 3.5, 3.5};
+
+ // identical to 80 MHz -- is this right
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau2 = {20.0, 20.4, 18.9, 20.8,
+ 19.1, 21.9, 22.6, 23.4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinSamples = {1, 1, 2, 1,
+ 1, 1, 1, 1};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsHG = {-8, -8, -8, -8,
+ -8, -8, -8, -8};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsLG = {-4, -4, -4, -4,
+ -4, -4, -4, -4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {800, 800, 800, 800, 800, 800, 800, 800};
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {1020, 1020, 1020, 1020, 1020, 1020, 1020, 1020};
+
+ // Set Tau and nominal timing offsets
+ ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr;
+
+ bool fixTau1 = true;
+ bool fixTau2 = true;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ fixTau1Arr[side][module] = fixTau1;
+ fixTau2Arr[side][module] = fixTau2;
+ }
+ }
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0HG = {53.942, 49.887, 59.633, 46.497,
+ 46.314, 42.267, 50.327, 41.605};
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0LG = {51.771, 47.936, 57.438, 44.191,
+ 44.295, 41.755, 48.081, 40.175};
+
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutHG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutLG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowHG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighHG = {8, 8, 8, 11, 8, 10, 8, 12};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowLG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighLG = {8, 8, 8, 11, 8, 10, 8, 12};
+
+ std::array<std::array<std::vector<float>, 4>, 2> slewingParamsHG, slewingParamsLG;
+
+ slewingParamsHG[0][0] = {0, -7.904e-02, 4.686e-02, 1.530e-03 };
+ slewingParamsHG[0][1] = {0, 2.250e-02, 4.732e-02, 6.050e-03 };
+ slewingParamsHG[0][2] = {0, 4.388e-02, 6.707e-02, -5.526e-05 };
+ slewingParamsHG[0][3] = {0, 1.205e-01, 2.726e-02, 2.610e-03 };
+
+ slewingParamsHG[1][0] = {0, 6.861e-02, 5.175e-03, -9.018e-04 };
+ slewingParamsHG[1][1] = {0, 3.855e-01, -4.442e-02, -2.022e-02 };
+ slewingParamsHG[1][2] = {0, -4.337e-03, 3.841e-02, 4.661e-03 };
+ slewingParamsHG[1][3] = {0, 3.623e-01, -3.882e-02, -1.805e-02 };
+
+ slewingParamsLG[0][0] = {0, 1.708e-02, 7.929e-02, 5.079e-03 };
+ slewingParamsLG[0][1] = {0, 1.406e-01, 1.209e-01, -1.922e-04 };
+ slewingParamsLG[0][2] = {0, 1.762e-01, 1.118e-01, 1.679e-04 };
+ slewingParamsLG[0][3] = {0, 1.361e-02, -2.685e-02, -4.168e-02 };
+
+ slewingParamsLG[1][0] = {0, 1.962e-01, -5.025e-03, -2.001e-02 };
+ slewingParamsLG[1][1] = {0, 3.258e-01, 1.229e-02, -2.925e-02 };
+ slewingParamsLG[1][2] = {0, 1.393e-01, 8.113e-02, -2.594e-03 };
+ slewingParamsLG[1][3] = {0, 1.939e-01, 2.188e-02, -5.579e-02 };
+
+ std::array<std::array<std::vector<float>, 4>, 2> moduleHGNonLinCorr;
+ moduleHGNonLinCorr[0][0] = {-3.76800e-02, 4.63597e-02};
+ moduleHGNonLinCorr[0][1] = {-1.02185e-01, -1.17548e-01};
+ moduleHGNonLinCorr[0][2] = {-8.78451e-02, -1.52174e-01};
+ moduleHGNonLinCorr[0][3] = {-1.04835e-01, -1.96514e-01};
+ moduleHGNonLinCorr[1][0] = {-6.83115e-02, 3.57802e-02};
+ moduleHGNonLinCorr[1][1] = {-1.08162e-01, -1.91413e-01};
+ moduleHGNonLinCorr[1][2] = {-7.82514e-02, -1.21218e-01};
+ moduleHGNonLinCorr[1][3] = {-2.34354e-02, -2.52033e-01};
+
+ m_zdcDataAnalyzer_40MHz = new ZDCDataAnalyzer(7,25,0,"FermiExp",peak2ndDerivMinSamples,
+ peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG,m_lowGainOnly);
+
+ m_zdcDataAnalyzer_40MHz->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+ m_zdcDataAnalyzer_40MHz->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0HG, t0LG);
+ m_zdcDataAnalyzer_40MHz->SetCutValues(chisqDivAmpCutHG, chisqDivAmpCutLG, DeltaT0CutLowHG, DeltaT0CutHighHG, DeltaT0CutLowLG, DeltaT0CutHighLG);
+ m_zdcDataAnalyzer_40MHz->SetTimingCorrParams(slewingParamsHG, slewingParamsLG);
+ m_zdcDataAnalyzer_40MHz->SetNonlinCorrParams(moduleHGNonLinCorr);
+
+}
+
+void ZdcAnalysisTool::initialize80MHz()
+{
+ // We have a complete configuration and so we override all of the default parameters
+ //
+
+ _peak2ndDerivMinSamples = {3, 2, 3, 2,
+ 2, 2, 2, 2};
+
+ _peak2ndDerivMinThresholdsHG = {-8, -8, -8, -8,
+ -8, -8, -8, -8};
+
+ _peak2ndDerivMinThresholdsLG = {-4, -4, -4, -4,
+ -4, -4, -4, -4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {800, 800, 800, 800, 800, 800, 800, 800};
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {950, 950, 950, 950, 950, 950, 950, 950};
+
+ // Set Tau and nominal timing offsets
+ ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr;
+
+ bool fixTau1 = true;
+ bool fixTau2 = true;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ fixTau1Arr[side][module] = fixTau1;
+ fixTau2Arr[side][module] = fixTau2;
+ }
+ }
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau1 = {3.9, 3.4, 4.1, 4.2,
+ 4.2, 3.6, 3.3, 3.4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau2 = {20.0, 20.4, 18.9, 20.8,
+ 19.1, 21.9, 22.6, 23.4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0HG = {44.24, 40.35, 49.3, 36.0,
+ 36.0, 31.1, 40.75, 30.5};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0LG = {42.65, 38.5, 47.4, 34,
+ 33.7, 29.9, 39.0, 29.3};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutHG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutLG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowHG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighHG = {8, 8, 8, 11, 8, 10, 8, 12};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowLG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighLG = {8, 8, 8, 11, 8, 10, 8, 12};
+
+ std::array<std::array<std::vector<float>, 4>, 2> slewingParamsHG, slewingParamsLG;
+
+ slewingParamsHG[0][0] = {0, -6.5e-2, 2.85e-2, -2.83e-3};
+ slewingParamsHG[0][1] = {0, -5.5e-2, 5.13e-2, 5.6e-3};
+ slewingParamsHG[0][2] = {0, -1.45e-3, 9.3e-2, 3.9e-3};
+ slewingParamsHG[0][3] = {0, -2.36e-2, 8.3e-2, 1.1e-3};
+
+ slewingParamsHG[1][0] = {0, -6.5e-2, 4.84e-2, -3.7e-3};
+ slewingParamsHG[1][1] = {0, 1.34e-2, 6.57e-2, 5.37e-3};
+ slewingParamsHG[1][2] = {0, -5.37e-2, 3.49e-2, 3.8e-3};
+ slewingParamsHG[1][3] = {0, -3.3e-2, 3.9e-2, 2.2e-3};
+
+ slewingParamsLG[0][0] = {0, -9.6e-2, 4.39e-2, 2.93e-3 };
+ slewingParamsLG[0][1] = {0, -5.0e-2, 14.9e-2, 20.6e-3 };
+ slewingParamsLG[0][2] = {0, -4.4e-2, 5.3e-2, 0, };
+ slewingParamsLG[0][3] = {0, -9.90e-2, 4.08e-2, 0, };
+
+ slewingParamsLG[1][0] = {0, -8.7e-2, 4.2e-2, -3.2e-3 };
+ slewingParamsLG[1][1] = {0, -3.26e-2, 3.84e-2, -2.32e-3};
+ slewingParamsLG[1][2] = {0, -26.8e-2, -2.64e-2, -5.3e-3 };
+ slewingParamsLG[1][3] = {0, -13.2e-2, 0.45e-2, -2.4e-3 };
+
+ std::array<std::array<std::vector<float>, 4>, 2> moduleHGNonLinCorr;
+ moduleHGNonLinCorr[0][0] = {-3.76800e-02, 4.63597e-02};
+ moduleHGNonLinCorr[0][1] = {-1.02185e-01, -1.17548e-01};
+ moduleHGNonLinCorr[0][2] = {-8.78451e-02, -1.52174e-01};
+ moduleHGNonLinCorr[0][3] = {-1.04835e-01, -1.96514e-01};
+ moduleHGNonLinCorr[1][0] = {-6.83115e-02, 3.57802e-02};
+ moduleHGNonLinCorr[1][1] = {-1.08162e-01, -1.91413e-01};
+ moduleHGNonLinCorr[1][2] = {-7.82514e-02, -1.21218e-01};
+ moduleHGNonLinCorr[1][3] = {-2.34354e-02, -2.52033e-01};
+
+ m_zdcDataAnalyzer_80MHz = new ZDCDataAnalyzer(7 , 12.5, 0, "FermiExp",_peak2ndDerivMinSamples,
+ _peak2ndDerivMinThresholdsHG, _peak2ndDerivMinThresholdsLG,m_lowGainOnly);
+
+ m_zdcDataAnalyzer_80MHz->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+ m_zdcDataAnalyzer_80MHz->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0HG, t0LG);
+ m_zdcDataAnalyzer_80MHz->SetCutValues(chisqDivAmpCutHG, chisqDivAmpCutLG, DeltaT0CutLowHG, DeltaT0CutHighHG, DeltaT0CutLowLG, DeltaT0CutHighLG);
+ m_zdcDataAnalyzer_80MHz->SetTimingCorrParams(slewingParamsHG, slewingParamsLG);
+ m_zdcDataAnalyzer_80MHz->SetNonlinCorrParams(moduleHGNonLinCorr);
+}
+
+StatusCode ZdcAnalysisTool::initializeTool()
+{
+ tf1SincInterp = new TF1("SincInterp",ZDC::SincInterp,-5.,160.,8);
+ tf1SincInterp->SetNpx(300);
+
+ // Set up calibrations
+ //
+ std::string filename = PathResolverFindCalibFile( "ZdcAnalysis/ZdcAnalysisConfig.conf" );
+ TEnv env(filename.c_str());
+
+ m_zdcEnergyCalibFileName = std::string(env.GetValue("ZdcEnergyCalibFileName","ZdcCalibrations_v1.root"));
+ ATH_MSG_INFO("ZDC energy calibration filename " << m_zdcEnergyCalibFileName);
+ m_zdcTimeCalibFileName = std::string(env.GetValue("ZdcTimeCalibFileName","ZdcTimeCalibrations_v1.root"));
+ ATH_MSG_INFO("ZDC time calibration filename " << m_zdcTimeCalibFileName);
+ m_zdcTriggerEffParamsFileName = std::string(env.GetValue("ZdcTriggerEffFileName", "ZdcTriggerEffParameters_v4.root"));
+ ATH_MSG_INFO("ZDC trigger efficiencies filename " << m_zdcTriggerEffParamsFileName);
+
+ ATH_MSG_INFO("Configuration: "<<m_configuration);
+ ATH_MSG_INFO("FlipEMDelay: "<<m_flipEMDelay);
+ ATH_MSG_INFO("LowGainOnly: "<<m_lowGainOnly);
+ ATH_MSG_INFO("WriteAux: "<<m_writeAux);
+ ATH_MSG_INFO("AuxSuffix: "<<m_auxSuffix);
+ ATH_MSG_INFO("DoCalib: "<<m_doCalib);
+ ATH_MSG_INFO("ForceCalibRun: "<<m_forceCalibRun);
+ ATH_MSG_INFO("ForceCalibLB: "<<m_forceCalibLB);
+ ATH_MSG_INFO("NumSampl: "<<m_numSample);
+ ATH_MSG_INFO("DeltaTSample: "<<m_deltaTSample);
+ ATH_MSG_INFO("Presample: "<<m_presample);
+ ATH_MSG_INFO("PeakSample: "<<m_peakSample);
+ ATH_MSG_INFO("Peak2ndDerivThresh: "<<m_Peak2ndDerivThresh);
+ ATH_MSG_INFO("T0: "<<m_t0);
+ ATH_MSG_INFO("Tau1: "<<m_tau1);
+ ATH_MSG_INFO("Tau2: "<<m_tau2);
+ ATH_MSG_INFO("FixTau1: "<<m_fixTau1);
+ ATH_MSG_INFO("FixTau2: "<<m_fixTau2);
+ ATH_MSG_INFO("DeltaTCut: "<<m_deltaTCut);
+ ATH_MSG_INFO("ChisqRatioCut: "<<m_ChisqRatioCut);
+
+ if (m_forceCalibRun>-1) {
+ ATH_MSG_INFO("CAREFUL: forcing calibration run/LB =" << m_forceCalibRun << "/" << m_forceCalibLB);
+
+ if (m_forceCalibLB < 0) {
+ ATH_MSG_ERROR("Invalid settings: Forced run > 0 but lumi block < 0");
+ return StatusCode::FAILURE;
+ }
+ }
+
+ // Use configuration to direct initialization
+ //
+ if (m_configuration == "default") {
+ m_zdcDataAnalyzer = initializeDefault();
+ }
+ else if (m_configuration == "PbPb2015") {
+ initialize80MHz();
+ initialize40MHz();
+
+ m_zdcDataAnalyzer = m_zdcDataAnalyzer_80MHz; // default
+ }
+ else {
+ ATH_MSG_ERROR("Unknown configuration: " << m_configuration);
+ return StatusCode::FAILURE;
+ }
+
+ // If an aux suffix is provided, prepend it with "_" so we don't have to do so at each use
+ //
+ if (m_writeAux && m_auxSuffix != "") {
+ m_auxSuffix = "_" + m_auxSuffix;
+ }
+
+
+ m_init = true;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode ZdcAnalysisTool::configureNewRun(unsigned int runNumber)
+{
+ ATH_MSG_INFO("Setting up new run " << runNumber);
+
+ // We do nothing for the default configuration
+ //
+ if (m_configuration != "default") {
+ if (m_configuration == "PbPb2015") {
+ //
+ // Two periods, 40 MHz and 80 MHz readout
+ //
+ if (runNumber < 287222) m_zdcDataAnalyzer = m_zdcDataAnalyzer_40MHz;
+ else m_zdcDataAnalyzer = m_zdcDataAnalyzer_80MHz;
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode ZdcAnalysisTool::recoZdcModule(const xAOD::ZdcModule& module)
+{
+ ATH_MSG_DEBUG("Processing ZDC module S/T/M/C = "
+ << module.side() << " "
+ << module.type() << " "
+ << module.zdcModule() << " "
+ << module.channel()
+ );
+ const std::vector<unsigned short>* adc0;
+ const std::vector<unsigned short>* adc1;
+
+
+ if (module.type()==0 && module.zdcModule()==0 && m_flipEMDelay) // flip delay/non-delay for EM big tube
+ {
+ adc0 = &(*(module.TTg0d1Link()))->adc();
+ adc1 = &(*(module.TTg1d1Link()))->adc();
+ }
+ else
+ {
+ adc0 = &(*(module.TTg0d0Link()))->adc();
+ adc1 = &(*(module.TTg1d0Link()))->adc();
+ }
+
+ float amp;
+ float time;
+ float qual;
+
+ float deltaT = m_deltaTSample;
+
+ sigprocMaxFinder(*adc0,deltaT,amp,time,qual);
+ if (m_writeAux) {
+ module.auxdecor<float>("amplitudeG0_mf" + m_auxSuffix)=amp;
+ module.auxdecor<float>("timeG0_mf" + m_auxSuffix)=time;
+ }
+
+ sigprocMaxFinder(*adc1,deltaT,amp,time,qual);
+ module.auxdecor<float>("amplitudeG1_mf" + m_auxSuffix)=amp;
+ module.auxdecor<float>("timeG1_mf" + m_auxSuffix)=time;
+
+ if (module.type()==0)
+ {
+ sigprocSincInterp(*adc0,deltaT,amp,time,qual);
+ module.auxdecor<float>("amplitudeG0_si" + m_auxSuffix)=amp;
+ module.auxdecor<float>("timeG0_si" + m_auxSuffix)=time;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& moduleContainer)
+{
+ if (!m_eventReady)
+ {
+ ATH_MSG_INFO("Event not ready for ZDC reco!");
+ return StatusCode::FAILURE;
+ }
+
+ const xAOD::EventInfo* eventInfo = 0;
+ ATH_CHECK(evtStore()->retrieve(eventInfo,"EventInfo"));
+
+ // check for new run number, if new, possibly update configuration and/or calibrations
+ //
+ unsigned int thisRunNumber = eventInfo->runNumber();
+ if (thisRunNumber != m_runNumber) {
+ ATH_MSG_DEBUG("ZDC analysis tool will be configured for run " << thisRunNumber);
+
+ ATH_CHECK(configureNewRun(thisRunNumber)); // ALWAYS check methods that return StatusCode
+
+ ATH_MSG_DEBUG("Setting up calibrations");
+
+ if (m_doCalib) {
+ //
+ // Check for calibration override
+ //
+ unsigned int calibRunNumber = thisRunNumber;
+ if (m_forceCalibRun > -1) calibRunNumber = m_forceCalibRun;
+
+ setEnergyCalibrations(calibRunNumber);
+ initializeTriggerEffs(calibRunNumber); // if energy calibrations fail to load, then so will trigger efficiencies
+ //setTimeCalibrations();
+ }
+
+ m_runNumber = thisRunNumber;
+ }
+
+ m_lumiBlock = eventInfo->lumiBlock();
+
+ unsigned int calibLumiBlock = m_lumiBlock;
+ if (m_doCalib) {
+ if (m_forceCalibRun > 0) calibLumiBlock = m_forceCalibLB;
+ }
+
+ ATH_MSG_DEBUG("Starting event processing");
+ m_zdcDataAnalyzer->StartEvent(calibLumiBlock);
+
+ const std::vector<unsigned short>* adc0;
+ const std::vector<unsigned short>* adc1;
+
+ ATH_MSG_DEBUG("Processing modules");
+ for (const auto zdcModule : moduleContainer)
+ {
+
+ if (zdcModule->type()!=0) continue;
+
+ if (zdcModule->zdcModule()==0 && m_flipEMDelay) // flip delay/non-delay for EM big tube
+ {
+ adc0 = &(*(zdcModule->TTg0d1Link()))->adc();
+ adc1 = &(*(zdcModule->TTg1d1Link()))->adc();
+ }
+ else
+ {
+ adc0 = &(*(zdcModule->TTg0d0Link()))->adc();
+ adc1 = &(*(zdcModule->TTg1d0Link()))->adc();
+ }
+
+ static std::vector<float> HGADCSamples(m_numSample);
+ static std::vector<float> LGADCSamples(m_numSample);
+
+ std::copy(adc0->begin(),adc0->end(),LGADCSamples.begin());
+ std::copy(adc1->begin(),adc1->end(),HGADCSamples.begin());
+
+ int side = (zdcModule->side()==-1) ? 0 : 1 ;
+
+ //std::cout << "LG: ";for (int i = 0;i<7;i++){std::cout<<LGADCSamples.at(i)<< " ";};std::cout<<std::endl;
+ //std::cout << "HG: ";for (int i = 0;i<7;i++){std::cout<<HGADCSamples.at(i)<< " ";};std::cout<<std::endl;
+
+ m_zdcDataAnalyzer->LoadAndAnalyzeData(side,zdcModule->zdcModule(),HGADCSamples, LGADCSamples);
+ }
+
+ ATH_MSG_DEBUG("Finishing event processing");
+
+ m_zdcDataAnalyzer->FinishEvent();
+
+ ATH_MSG_DEBUG("Adding variables");
+
+ for (const auto zdcModule : moduleContainer)
+ {
+
+ if (zdcModule->type()!=0) continue;
+
+ int side = (zdcModule->side()==-1) ? 0 : 1 ;
+ int mod = zdcModule->zdcModule();
+
+ if (m_writeAux) {
+ if (m_doCalib) {
+ float calibEnergy = m_zdcDataAnalyzer->GetModuleCalibAmplitude(side, mod);
+ zdcModule->auxdecor<float>("CalibEnergy" + m_auxSuffix) = calibEnergy;
+ zdcModule->auxdecor<float>("CalibTime" + m_auxSuffix) = m_zdcDataAnalyzer->GetModuleCalibTime(side, mod);
+ }
+ else
+ {
+ zdcModule->auxdecor<float>("CalibEnergy" + m_auxSuffix) = -1000;
+ zdcModule->auxdecor<float>("CalibTime" + m_auxSuffix) = -1000;
+ }
+
+ zdcModule->auxdecor<unsigned int>("Status" + m_auxSuffix) = m_zdcDataAnalyzer->GetModuleStatus(side, mod);
+ zdcModule->auxdecor<float>("Amplitude" + m_auxSuffix) = m_zdcDataAnalyzer->GetModuleAmplitude(side, mod);
+ zdcModule->auxdecor<float>("Time" + m_auxSuffix) = m_zdcDataAnalyzer->GetModuleTime(side, mod);
+
+ const ZDCPulseAnalyzer* pulseAna_p = m_zdcDataAnalyzer->GetPulseAnalyzer(side, mod);
+ zdcModule->auxdecor<float>("Chisq" + m_auxSuffix) = pulseAna_p->GetChisq();
+ zdcModule->auxdecor<float>("FitAmp" + m_auxSuffix) = pulseAna_p->GetFitAmplitude();
+ zdcModule->auxdecor<float>("FitAmpError" + m_auxSuffix) = pulseAna_p->GetAmpError();
+ zdcModule->auxdecor<float>("FitT0" + m_auxSuffix) = pulseAna_p->GetFitT0();
+ zdcModule->auxdecor<float>("BkgdMaxFraction" + m_auxSuffix) = pulseAna_p->GetBkgdMaxFraction();
+ }
+ /*
+ std::cout << "side = " << side << " module=" << zdcModule->zdcModule() << " CalibEnergy=" << zdcModule->auxdecor<float>("CalibEnergy")
+ << " should be " << m_zdcDataAnalyzer->GetModuleCalibAmplitude(side,mod) << std::endl;
+ */
+ }
+
+ // Record sum objects
+
+ //xAOD::ZdcModuleContainer* newModuleContainer = new xAOD::ZdcModuleContainer();
+ //xAOD::ZdcModuleAuxContainer* newModuleAuxContainer = new xAOD::ZdcModuleAuxContainer() ;
+
+ std::unique_ptr<xAOD::ZdcModuleContainer> newModuleContainer( new xAOD::ZdcModuleContainer() );
+ std::unique_ptr<xAOD::ZdcModuleAuxContainer> newModuleAuxContainer( new xAOD::ZdcModuleAuxContainer() );
+
+ newModuleContainer->setStore( newModuleAuxContainer.get() );
+
+ for (int iside=0;iside<2;iside++)
+ {
+ xAOD::ZdcModule* zdc_sum = new xAOD::ZdcModule;
+ newModuleContainer.get()->push_back(zdc_sum);
+ zdc_sum->setSide(iside);
+
+ float calibEnergy = getCalibModuleSum(iside);
+ zdc_sum->auxdecor<float>("CalibEnergy") = calibEnergy;
+
+ float finalEnergy = calibEnergy;
+ if (iside==0) finalEnergy = finalEnergy*(1 + 7e-7 * finalEnergy); // nonlinear correction for side C
+
+ zdc_sum->auxdecor<float>("FinalEnergy") = finalEnergy;
+ zdc_sum->auxdecor<float>("AverageTime") = getAverageTime(iside);
+ zdc_sum->auxdecor<unsigned int>("Status") = !sideFailed(iside);
+ zdc_sum->auxdecor<unsigned int>("ModuleMask") = (getModuleMask()>>(4*iside)) & 0xF;
+ }
+
+ ATH_CHECK( evtStore()->record( newModuleContainer.release() , "ZdcSums" + m_auxSuffix) ) ;
+ ATH_CHECK( evtStore()->record( newModuleAuxContainer.release() , "ZdcSums" + m_auxSuffix +"Aux.") );
+
+ return StatusCode::SUCCESS;
+}
+
+void ZdcAnalysisTool::setEnergyCalibrations(unsigned int runNumber)
+{
+
+ char name[128];
+ /*
+ sprintf(name,"%s/%s",m_zdcAnalysisConfigPath.c_str(),m_zdcEnergyCalibFileName.c_str());
+ ATH_MSG_INFO("Opening energy calibration file " << name);
+ TFile* fCalib = TFile::Open(name);
+ */
+
+ std::string filename = PathResolverFindCalibFile( ("ZdcAnalysis/"+m_zdcEnergyCalibFileName).c_str() );
+ ATH_MSG_INFO("Opening energy calibration file " << filename);
+ TFile* fCalib = TFile::Open(filename.c_str());
+
+ //std::array<std::array<TSpline*, 4>, 2> splines;
+
+ for (int iside=0;iside<2;iside++)
+ {
+ for(int imod=0;imod<4;imod++)
+ {
+ // need to delete old splines
+ //std::cout << "m_splines[iside][imod] = " << m_splines[iside][imod] << std::endl;
+ SafeDelete( m_splines[iside][imod] );
+
+ sprintf(name,"ZDC_Gcalib_run%d_s%d_m%d",runNumber,iside,imod);
+ ATH_MSG_INFO("Searching for graph " << name);
+ TGraph* g = (TGraph*) fCalib->GetObjectChecked(name,"TGraph");
+ if (!g && m_doCalib)
+ {
+ ATH_MSG_ERROR("No calibrations for run " << runNumber);
+ m_doCalib = false;
+ }
+
+ if (g)
+ {
+ m_splines[iside][imod] = new TSpline3(name,g);
+ }
+ else
+ {
+ ATH_MSG_ERROR("No graph " << name);
+ }
+ }
+ }
+ fCalib->Close();
+ if (m_doCalib) m_zdcDataAnalyzer->LoadEnergyCalibrations(m_splines);
+
+ return;
+}
+
+void ZdcAnalysisTool::setTimeCalibrations(unsigned int runNumber)
+{
+ char name[128];
+ sprintf(name,"%s/%s",m_zdcAnalysisConfigPath.c_str(),m_zdcTimeCalibFileName.c_str());
+ ATH_MSG_INFO("Opening time calibration file " << name);
+ TFile* fCalib = TFile::Open(name);
+ std::array<std::array<TSpline*, 4>, 2> T0HGOffsetSplines;
+ std::array<std::array<TSpline*, 4>, 2> T0LGOffsetSplines;
+ TSpline3* spline;
+ for (int iside=0;iside<2;iside++)
+ {
+ for(int imod=0;imod<4;imod++)
+ {
+ sprintf(name,"ZDC_T0calib_run%d_s%d_m%d_HG",runNumber,iside,imod);
+ spline = (TSpline3*) fCalib->GetObjectChecked(name,"TSpline3");
+ T0HGOffsetSplines[iside][imod] = spline;
+ sprintf(name,"ZDC_T0calib_run%d_s%d_m%d_LG",runNumber,iside,imod);
+ spline = (TSpline3*) fCalib->GetObjectChecked(name,"TSpline3");
+ T0LGOffsetSplines[iside][imod] = spline;
+ }
+ }
+ m_zdcDataAnalyzer->LoadT0Calibrations(T0HGOffsetSplines,T0LGOffsetSplines);
+ fCalib->Close();
+}
+
+StatusCode ZdcAnalysisTool::reprocessZdc()
+{
+ if (!m_init)
+ {
+ ATH_MSG_INFO("Tool not initialized!");
+ return StatusCode::FAILURE;
+ }
+ m_eventReady = false;
+ //std::cout << "Trying to retrieve " << m_zdcModuleContainerName << std::endl;
+ m_zdcModules = 0;
+ ATH_CHECK(evtStore()->retrieve(m_zdcModules,m_zdcModuleContainerName));
+ m_eventReady = true;
+
+ ATH_CHECK(recoZdcModules(*m_zdcModules));
+
+ return StatusCode::SUCCESS;
+}
+
+ bool ZdcAnalysisTool::sigprocMaxFinder(const std::vector<unsigned short>& adc, float deltaT, float& amp, float& time, float& qual)
+ {
+ size_t nsamp = adc.size();
+ float presamp = adc.at(0);
+ unsigned short max_adc = 0;
+ int max_index = -1;
+ for (size_t i = 0;i<nsamp;i++)
+ {
+ if (adc[i]>max_adc)
+ {
+ max_adc = adc[i];
+ max_index = i;
+ }
+ }
+ amp = max_adc - presamp;
+ time = max_index*deltaT;
+ qual = 1.;
+
+ if(max_index==-1)
+ {
+ qual=0.;
+ return false;
+ }
+
+ return true;
+ }
+
+ bool ZdcAnalysisTool::sigprocSincInterp(const std::vector<unsigned short>& adc, float deltaT, float& amp, float& time, float& qual)
+ {
+ size_t nsamp = adc.size();
+ float presamp = adc.at(0);
+ tf1SincInterp->SetParameter(0,deltaT);
+ for (size_t i = 0;i<nsamp;i++)
+ {
+ tf1SincInterp->SetParameter(i+1,adc.at(i)-presamp);
+ }
+ amp = tf1SincInterp->GetMaximum();
+ time = tf1SincInterp->GetMaximumX();
+ qual = 1.;
+ return true;
+ }
+
+ float ZdcAnalysisTool::getModuleSum(int side)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleSum(side);
+ }
+
+/*
+ float ZdcAnalysisTool::getModuleAmplitude(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetPulseAnalyzer(side,mod)->GetAmplitude();
+ }
+
+ float ZdcAnalysisTool::getModuleFitAmplitude(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetPulseAnalyzer(side,mod)->GetFitAmplitude();
+ }
+
+ float ZdcAnalysisTool::getModuleFitT0(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetPulseAnalyzer(side,mod)->GetFitT0();
+ }
+
+ float ZdcAnalysisTool::getModuleTime(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleTime(side,mod);
+ }
+
+ float ZdcAnalysisTool::getModuleCalibAmplitude(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleCalibAmplitude(side,mod);
+ }
+
+ float ZdcAnalysisTool::getModuleCalibTime(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleCalibTime(side,mod);
+ }
+
+ float ZdcAnalysisTool::getModuleStatus(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleStatus(side,mod);
+ }
+
+ float ZdcAnalysisTool::getModuleChisq(int side, int mod)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleChisq(side,mod);
+ }
+*/
+
+ float ZdcAnalysisTool::getCalibModuleSum(int side)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetCalibModuleSum(side);
+ }
+
+ float ZdcAnalysisTool::getAverageTime(int side)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetAverageTime(side);
+ }
+
+ bool ZdcAnalysisTool::sideFailed(int side)
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->SideFailed(side);
+ }
+
+ unsigned int ZdcAnalysisTool::getModuleMask()
+ {
+ if (!m_zdcDataAnalyzer) return 0;
+ return m_zdcDataAnalyzer->GetModuleMask();
+ }
+
+ float ZdcAnalysisTool::getTriggerEfficiency(int side)
+ {
+ if (!m_doCalib) return -1;
+
+ m_zdcTriggerEfficiency->UpdatelumiBlock(m_lumiBlock);
+ float adcSum = getModuleSum(side);
+ float eff = m_zdcTriggerEfficiency->GetEfficiency(side, adcSum);
+ //if(eff<0) std::cout << "eff= " << eff << " --> Invalid lumiBlock input: m_lumiBlock = " << m_lumiBlock << std::endl; // looking for -1's
+ return eff;
+ }
+
+} // namespace ZDC
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysis_entries.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysis_entries.cxx
new file mode 100644
index 000000000000..ff4c36e6867d
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysis_entries.cxx
@@ -0,0 +1,16 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef XAOD_STANDALONE
+
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "ZdcAnalysis/ZdcAnalysisTool.h"
+
+DECLARE_TOOL_FACTORY (ZDC::ZdcAnalysisTool)
+
+DECLARE_FACTORY_ENTRIES (ZdcAnalysis) {
+ DECLARE_TOOL (ZDC::ZdcAnalysisTool)
+}
+
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcSincInterp.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcSincInterp.cxx
new file mode 100644
index 000000000000..7593c7e6a67d
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcSincInterp.cxx
@@ -0,0 +1,27 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ZdcAnalysis/ZdcSincInterp.h"
+#include "TMath.h"
+#include <cmath>
+
+namespace ZDC
+{
+ double SincInterp(double* xvec, double* pvec)
+ {
+ // pvec are the sample values
+ double ret = 0;
+ double T = pvec[0]; // deltaT
+ double t = xvec[0];
+ for (int isamp = 0;isamp<7;isamp++)
+ {
+ double arg = (t - isamp*T)/T;
+ if (arg!=0.0)
+ {
+ ret += pvec[isamp+1] * std::sin(TMath::Pi()*arg)/(TMath::Pi()*arg);
+ }
+ }
+ return ret;
+ }
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/IZdcAnalysisTool.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/IZdcAnalysisTool.h
new file mode 100644
index 000000000000..c4a9e9bdeb9f
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/IZdcAnalysisTool.h
@@ -0,0 +1,29 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef __IZDCRECTOOL_H__
+#define __IZDCRECTOOL_H__
+
+#include "AsgTools/IAsgTool.h"
+#include "xAODForward/ZdcModuleContainer.h"
+
+namespace ZDC
+{
+
+class IZdcAnalysisTool : virtual public asg::IAsgTool
+{
+ ASG_TOOL_INTERFACE( ZDC::IZdcAnalysisTool )
+
+ public:
+
+ /// Initialize the tool.
+ virtual StatusCode initializeTool() = 0;
+ virtual StatusCode recoZdcModule(const xAOD::ZdcModule& module) = 0;
+ virtual StatusCode recoZdcModules(const xAOD::ZdcModuleContainer& moduleContainer) = 0;
+ virtual StatusCode reprocessZdc() = 0;
+};
+
+}
+
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
new file mode 100644
index 000000000000..9bcafdd3cd29
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
@@ -0,0 +1,153 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef _ZDCDataAnalyzer_h
+#define _ZDCDataAnalyzer_h
+
+#include <ZdcAnalysis/ZDCPulseAnalyzer.h>
+#include <TSpline.h>
+
+#include <array>
+#include <string>
+
+class ZDCDataAnalyzer
+{
+public:
+ typedef std::array<std::array<float, 4>, 2> ZDCModuleFloatArray;
+ typedef std::array<std::array<bool, 4>, 2> ZDCModuleBoolArray;
+
+private:
+ size_t _nSample;
+ float _deltaTSample;
+ size_t _preSampleIdx;
+ std::string _fitFunction;
+ bool _forceLG;
+
+ std::array<std::array<ZDCPulseAnalyzer*, 4>, 2> _moduleAnalyzers;
+
+ int _debugLevel;
+ int _eventCount;
+
+ ZDCModuleFloatArray _HGGains;
+ ZDCModuleFloatArray _pedestals;
+
+ bool _haveECalib;
+ std::array<std::array<TSpline*, 4>, 2> _LBDepEcalibSplines;
+
+ bool _haveT0Calib;
+ std::array<std::array<TSpline*, 4>, 2> _T0HGOffsetSplines;
+ std::array<std::array<TSpline*, 4>, 2> _T0LGOffsetSplines;
+
+ // Transient data that is updated each LB or each event
+ //
+ int _currentLB;
+ ZDCModuleFloatArray _currentECalibCoeff;
+ ZDCModuleFloatArray _currentT0OffsetsHG;
+ ZDCModuleFloatArray _currentT0OffsetsLG;
+
+ std::array<std::array<bool, 4>, 2> _dataLoaded;
+ // std::array<std::array<bool, 4>, 2> _moduleFail;
+
+ unsigned int _moduleMask;
+
+ std::array<std::array<unsigned int, 4>, 2> _moduleStatus;
+ std::array<std::array<float, 4>, 2> _calibAmplitude;
+ std::array<std::array<float, 4>, 2> _calibTime;
+
+ std::array<float, 2> _moduleSum;
+ std::array<float, 2> _moduleSumErrSq;
+ std::array<float, 2> _moduleSumPreSample;
+
+ std::array<float, 2> _calibModuleSum;
+ std::array<float, 2> _calibModuleSumErrSq;
+
+ std::array<float, 2> _averageTime;
+ std::array<bool, 2> _fail;
+
+public:
+
+ ZDCDataAnalyzer(int nSample, float deltaTSample, size_t preSampleIdx, std::string fitFunction,
+ const ZDCModuleFloatArray& peak2ndDerivMinSamples,
+ const ZDCModuleFloatArray& peak2ndDerivMinThresholdsHG,
+ const ZDCModuleFloatArray& peak2ndDerivMinThresholdsLG,
+ bool forceLG = false);
+
+ ~ZDCDataAnalyzer();
+
+ void SetDebugLevel(int level = 0) {
+ _debugLevel = level;
+ if (level < 2) ZDCPulseAnalyzer::SetQuietFits(true);
+ else ZDCPulseAnalyzer::SetQuietFits(false);
+ }
+
+ unsigned int GetModuleMask() const {return _moduleMask;}
+
+ float GetModuleSum(size_t side) const {return _moduleSum.at(side);}
+ float GetModuleSumErr(size_t side) const {return std::sqrt(_moduleSumErrSq.at(side));}
+
+ float GetCalibModuleSum(size_t side) const {return _calibModuleSum.at(side);}
+ float GetCalibModuleSumErr(size_t side) const {return std::sqrt(_calibModuleSumErrSq.at(side));}
+
+ float GetModuleSumPreSample(size_t side) const {return _moduleSumPreSample.at(side);}
+
+ float GetAverageTime(size_t side) const {return _averageTime.at(side);}
+ bool SideFailed(size_t side) const {return _fail.at(side);}
+
+ float GetModuleAmplitude(size_t side, size_t module) const {return _moduleAnalyzers.at(side).at(module)->GetAmplitude();}
+ float GetModuleTime(size_t side, size_t module) const {return _moduleAnalyzers.at(side).at(module)->GetT0Corr();}
+ float GetModuleChisq(size_t side, size_t module) const {return _moduleAnalyzers.at(side).at(module)->GetChisq();}
+
+ float GetModuleCalibAmplitude(size_t side, size_t module) const {return _calibAmplitude.at(side).at(module);}
+ float GetModuleCalibTime(size_t side, size_t module) const {return _calibTime.at(side).at(module);}
+ float GetModuleStatus(size_t side, size_t module) const {return _moduleStatus.at(side).at(module);}
+
+ const ZDCPulseAnalyzer* GetPulseAnalyzer(size_t side, size_t module) const {return _moduleAnalyzers.at(side).at(module);}
+
+ void SetADCOverUnderflowValues(const ZDCModuleFloatArray& HGOverflowADC, const ZDCModuleFloatArray& HGUnderflowADC,
+ const ZDCModuleFloatArray& LGOverflowADC);
+
+ void SetTauT0Values(const ZDCModuleBoolArray& fxiTau1, const ZDCModuleBoolArray& fxiTau2,
+ const ZDCModuleFloatArray& tau1, const ZDCModuleFloatArray& tau2,
+ const ZDCModuleFloatArray& t0HG, const ZDCModuleFloatArray& t0LG);
+
+ void SetCutValues(const ZDCModuleFloatArray& chisqDivAmpCutHG, const ZDCModuleFloatArray& chisqDivAmpCutLG,
+ const ZDCModuleFloatArray& deltaT0MinHG, const ZDCModuleFloatArray& deltaT0MaxHG,
+ const ZDCModuleFloatArray& deltaT0MinLG, const ZDCModuleFloatArray& deltaT0MaxLG);
+
+
+ void SetTimingCorrParams(const std::array<std::array<std::vector<float>, 4>, 2>& HGParamArr,
+ const std::array<std::array<std::vector<float>, 4>, 2>& LGParamArr);
+
+ void SetNonlinCorrParams(const std::array<std::array<std::vector<float>, 4>, 2>& HGNonlinCorrParams);
+
+ void LoadEnergyCalibrations(const std::array<std::array<TSpline*, 4>, 2>& calibSplines)
+ {
+ if (_debugLevel > 0) {
+ std::cout << "Loading energy calibrations" << std::endl;
+ }
+
+ _LBDepEcalibSplines = calibSplines;
+ _haveECalib = true;
+ }
+
+ void LoadT0Calibrations(const std::array<std::array<TSpline*, 4>, 2>& T0HGOffsetSplines,
+ const std::array<std::array<TSpline*, 4>, 2>& T0LGOffsetSplines)
+ {
+ if (_debugLevel > 0) {
+ std::cout << "Loading timing calibrations" << std::endl;
+ }
+ _T0HGOffsetSplines = T0HGOffsetSplines;
+ _T0LGOffsetSplines = T0LGOffsetSplines;
+
+ _haveT0Calib = true;
+ }
+
+ void StartEvent(int lumiBlock);
+
+ void LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples);
+
+ bool FinishEvent();
+
+};
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
new file mode 100644
index 000000000000..e24e80bb0ee3
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
@@ -0,0 +1,255 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef _ZDCFitWrapper_h
+#define _ZDCFitWrapper_h
+
+// Base class that defines the interface
+//
+#include <TF1.h>
+
+double ZDCFermiExpFit(double* xvec, double* pvec);
+
+class ZDCFitWrapper
+{
+ TF1* _wrapperTF1;
+
+public:
+ ZDCFitWrapper(TF1* wrapperTF1) : _wrapperTF1(wrapperTF1) {}
+
+ virtual ~ZDCFitWrapper() {delete _wrapperTF1;}
+
+ virtual void Initialize(float initialAmp, float initialT0) = 0;
+
+ virtual float GetAmplitude() const = 0;
+ virtual float GetAmpError() const = 0;
+ virtual float GetTime() const = 0;
+ virtual float GetTau1() const = 0;
+ virtual float GetTau2() const = 0;
+
+ virtual float GetBkgdMaxFraction() const = 0;
+
+ virtual float GetShapeParameter(size_t index) const = 0;
+
+ virtual double operator()(double *x, double *p) = 0;
+
+ virtual TF1* GetWrapperTF1() {return _wrapperTF1;}
+ virtual const TF1* GetWrapperTF1() const {return _wrapperTF1;}
+};
+
+class ZDCPrePulseFitWrapper : public ZDCFitWrapper
+{
+public:
+ ZDCPrePulseFitWrapper(TF1* wrapperTF1) : ZDCFitWrapper(wrapperTF1) {}
+
+ virtual void SetInitialPrePulseT0(float t0) = 0;
+};
+
+class ZDCFitExpFermiVariableTaus : public ZDCFitWrapper
+{
+ bool _fixTau1;
+ bool _fixTau2;
+
+ float _tau1;
+ float _tau2;
+
+public:
+
+ ZDCFitExpFermiVariableTaus(std::string tag, float tmin, float tmax, bool fixTau1, bool fixTau2, float tau1, float tau2);
+
+ virtual void Initialize(float initialAmp, float initialT0);
+
+ virtual float GetAmplitude() const {return GetWrapperTF1()->GetParameter(0); }
+ virtual float GetAmpError() const {return GetWrapperTF1()->GetParError(0); }
+
+ virtual float GetTau1() const {return GetWrapperTF1()->GetParameter(2);}
+ virtual float GetTau2() const {return GetWrapperTF1()->GetParameter(3);}
+
+ virtual float GetTime() const {
+ const TF1* theTF1 = GetWrapperTF1();
+
+ float fitT0 = theTF1->GetParameter(1);
+
+ float tau1 = theTF1->GetParameter(2);
+ float tau2 = theTF1->GetParameter(3);
+
+ // Correct the time to the maximum (the factor of 1/2 is still not fully understood)
+ //
+ if (tau2 > tau1) fitT0 += 0.5*tau1*std::log(tau2/tau1 - 1.0);
+ return fitT0;
+ }
+
+ virtual float GetShapeParameter(size_t index) const
+ {
+ if (index == 0) return GetWrapperTF1()->GetParameter(2);
+ else if (index == 1) return GetWrapperTF1()->GetParameter(3);
+ else throw;
+ }
+
+ virtual float GetBkgdMaxFraction() const
+ {
+ const TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+ double amp = theTF1->GetParameter(0);
+ double slope = theTF1->GetParameter(4);
+
+ double background = slope*GetTime();
+ return background/amp;
+ }
+
+ // virtual float GetNDOF() const {return _fitFunc->GetNDF(); }
+
+ virtual double operator()(double *x, double *p) {
+ return ZDCFermiExpFit(x, p);
+ }
+
+};
+
+class ZDCFitExpFermiFixedTaus : public ZDCFitWrapper
+{
+ float _tau1;
+ float _tau2;
+
+ float _norm;
+ float _timeCorr;
+
+ TF1* _expFermiFunc;
+
+public:
+
+ ZDCFitExpFermiFixedTaus(std::string tag, float tmin, float tmax, float tau1, float tau2);
+
+ ~ZDCFitExpFermiFixedTaus() {
+ delete _expFermiFunc;
+ }
+
+ virtual void Initialize(float initialAmp, float initialT0);
+
+ virtual float GetAmplitude() const {return GetWrapperTF1()->GetParameter(0); }
+ virtual float GetAmpError() const {return GetWrapperTF1()->GetParError(0); }
+
+ virtual float GetTau1() const {return _tau1;}
+ virtual float GetTau2() const {return _tau2;}
+
+ virtual float GetTime() const {
+ float fitT0 = GetWrapperTF1()->GetParameter(1);
+
+ // Correct the time to the maximum (the factor of 1/2 is still not fully understood)
+ //
+ fitT0 += _timeCorr;
+ return fitT0;
+ }
+
+ virtual float GetShapeParameter(size_t index) const
+ {
+ if (index == 0) return _tau1;
+ else if (index == 1) return _tau2;
+ else throw;
+ }
+
+ virtual float GetBkgdMaxFraction() const
+ {
+ const TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+ double amp = theTF1->GetParameter(0);
+ double slope = theTF1->GetParameter(2);
+
+ double background = slope*GetTime();
+ return background/amp;
+ }
+
+ virtual double operator() (double *x, double *p)
+ {
+ double amp = p[0];
+ double t0 = p[1];
+ double deltaT = x[0] - t0;
+
+ double bckgd = p[2]*x[0];
+
+ double expFermi = amp*_norm*_expFermiFunc->operator()(deltaT);
+
+ return expFermi + bckgd;
+ }
+};
+
+class ZDCFitExpFermiPrePulse : public ZDCPrePulseFitWrapper
+{
+ float _tau1;
+ float _tau2;
+ float _norm;
+ float _timeCorr;
+ TF1* _expFermiFunc;
+
+public:
+ ZDCFitExpFermiPrePulse(std::string tag, float tmin, float tmax, float tau1, float tau2);
+
+ virtual void Initialize(float initialAmp, float initialT0);
+ // void InitializePrePulseT0(float initPreT0) { GetWrapperTF1()->SetParameter(3, initPreT0);}
+
+ virtual void SetInitialPrePulseT0(float t0) {GetWrapperTF1()->SetParameter(3, t0);}
+
+ virtual float GetAmplitude() const {return GetWrapperTF1()->GetParameter(0); }
+ virtual float GetAmpError() const {return GetWrapperTF1()->GetParError(0); }
+
+ virtual float GetTau1() const {return GetWrapperTF1()->GetParameter(2);}
+ virtual float GetTau2() const {return GetWrapperTF1()->GetParameter(3);}
+
+ virtual float GetTime() const {
+ float fitT0 = GetWrapperTF1()->GetParameter(1);
+
+ // Correct the time to the maximum (the factor of 1/2 is still not fully understood)
+ //
+ fitT0 += _timeCorr;
+ return fitT0;
+ }
+
+ virtual float GetShapeParameter(size_t index) const
+ {
+ if (index == 0) return _tau1;
+ else if (index == 1) return _tau2;
+ else if (index < 5) return GetWrapperTF1()->GetParameter(index);
+ else throw;
+ }
+
+ virtual float GetBkgdMaxFraction() const
+ {
+ const TF1* theTF1 = ZDCFitWrapper::GetWrapperTF1();
+
+ double maxTime = GetTime();
+
+ double amp = theTF1->GetParameter(0);
+ double preAmp = theTF1->GetParameter(2);
+ double preT0 = theTF1->GetParameter(3);
+ double slope = theTF1->GetParameter(4);
+
+ double deltaTPre = maxTime - preT0;
+
+ double background = slope*maxTime + preAmp*_norm*(_expFermiFunc->operator()(deltaTPre) -
+ _expFermiFunc->operator()(-preT0));
+
+ return background/amp;
+ }
+
+ virtual double operator() (double *x, double *p)
+ {
+ double t = x[0];
+
+ double amp = p[0];
+ double t0 = p[1];
+ double preAmp = p[2];
+ double preT0 = p[3];
+ double linSlope = p[4];
+
+ double deltaT = t - t0;
+ double deltaTPre = t - preT0;
+
+ double bckgd = linSlope*t;
+
+ double pulse1 = amp*_norm*_expFermiFunc->operator()(deltaT);
+ double pulse2 = preAmp*_norm*(_expFermiFunc->operator()(deltaTPre) -
+ _expFermiFunc->operator()(-preT0));
+
+ return pulse1 + pulse2 + bckgd;
+ }
+};
+
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
new file mode 100644
index 000000000000..8939f8170715
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
@@ -0,0 +1,290 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef _ZDCPulseAnalyzer_h
+#define _ZDCPulseAnalyzer_h
+
+//#include <valarray>
+#include <vector>
+#include <string>
+
+#include <TF1.h>
+#include <TH1.h>
+
+#include "ZdcAnalysis/ZDCFitWrapper.h"
+
+class ZDCPulseAnalyzer
+{
+public:
+ enum {PulseBit = 0, LowGainBit, FailBit, HGOverflowBit,
+ HGUnderflowBit, PSHGOverUnderflowBit, LGOverflowBit, LGUnderflowBit,
+ PrePulseBit, PostPulseBit, FitFailedBit, BadChisqBit, BadT0Bit};
+
+private:
+ typedef std::vector<float>::const_iterator SampleCIter;
+
+ // Static data
+ //
+ static std::string _fitOptions;
+ static bool _quietFits;
+
+ // Quantities provided/set in the constructor
+ //
+ std::string _tag;
+ size_t _Nsample;
+ size_t _preSampleIdx;
+ float _deltaTSample;
+ int _pedestal;
+ float _gainHG;
+ bool _forceLG;
+ float _tmin;
+ float _tmax;
+ std::string _fitFunction;
+ float _peak2ndDerivMinSample;
+ float _peak2ndDerivMinThreshLG;
+ float _peak2ndDerivMinThreshHG;
+
+ // Default fit values and cuts that can be set via modifier methods
+ //
+ int _HGOverflowADC;
+ int _HGUnderflowADC;
+ int _LGOverflowADC;
+
+ float _nominalT0HG;
+ float _nominalT0LG;
+
+ float _nominalTau1;
+ float _nominalTau2;
+
+ bool _fixTau1;
+ bool _fixTau2;
+
+ float _chisqDivAmpCutLG; // maximum good LG chisq / amplitude
+ float _chisqDivAmpCutHG; // maximum good HG chisq / amplitude
+
+ float _T0CutLowLG; // minimum good corrected time for LG fits
+ float _T0CutHighLG; // maximum good corrected time for LG fits
+
+ float _T0CutLowHG; // minimum good corrected time for HG fits
+ float _T0CutHighHG; // maximum good corrected time for HG fits
+
+ bool _haveTimingCorrections;
+ std::vector<float> _LGT0CorrParams; // Parameters used to correct the fit LG times
+ std::vector<float> _HGT0CorrParams; // Parameters used to correct the fit HG times
+
+ bool _haveNonlinCorr;
+ std::vector<float> _nonLinCorrParams;
+
+ // Histogram used to perform the fits and function wrappers
+ //
+ mutable TH1* _fitHist;
+
+ bool _initializedFits;
+ ZDCFitWrapper* _defaultFitWrapper;
+ ZDCFitWrapper* _prePulseFitWrapper;
+
+ // TH1* _LGFitHist;
+ // TH1* _fitHist;
+ // TF1* _FitFuncNoPrePulse;
+ // TF1* _FitFuncPrePulse;
+
+ // Dynamic data loaded for each pulse (event)
+ // ==========================================
+
+ // Statuses
+ //
+ bool _haveData;
+ bool _havePulse;
+ bool _fail;
+ bool _useLowGain;
+ bool _HGOverflow;
+ bool _HGUnderflow;
+ bool _LGOverflow;
+ bool _LGUnderflow;
+ bool _PSHGOverUnderflow;
+ bool _prePulse;
+ bool _postPulse;
+ bool _fitFailed;
+ bool _badT0;
+ bool _badChisq;
+
+ // Pulse analysis
+ //
+ float _preSample;
+
+ float _maxADCValue;
+ float _minADCValue;
+ float _maxDelta;
+ float _minDelta;
+
+ int _maxSampl;
+ int _minSampl;
+ int _maxDeltaSampl;
+ int _minDeltaSampl;
+
+ float _minDeriv2nd;
+ int _minDeriv2ndIndex;
+
+ float _fitAmplitude;
+ float _fitAmpError;
+ float _fitT0;
+ float _fitT0Corr;
+ float _fitTau1;
+ float _fitTau2;
+ float _fitChisq;
+ float _amplitude;
+ float _ampError;
+ float _preSampleAmp;
+ float _bkgdMaxFraction;
+
+ std::vector<float> _ADCSamplesHGSub;
+ std::vector<float> _ADCSamplesLGSub;
+
+ std::vector<float> _ADCSSampSigHG;
+ std::vector<float> _ADCSSampSigLG;
+
+ std::vector<float> _samplesSub;
+ std::vector<float> _samplesDeriv;
+ std::vector<float> _samplesDeriv2nd;
+
+ // Private methods
+ //
+ void Reset();
+ void SetDefaults();
+ void SetupFitFunctions();
+
+ bool AnalyzeData(const std::vector<float>& samples, // The samples used for this event
+ const std::vector<float>& samplesSig, // The "resolution" on the ADC value
+ const std::vector<float>& toCorrParams, // The parameters used to correct the t0
+ float maxChisqDivAmp, // The maximum chisq / amplitude ratio
+ float minT0Corr, float maxT0Corr, // The minimum and maximum corrected T0 values
+ float peak2ndDerivMinThresh
+ );
+
+
+ void FillHistogram(const std::vector<float>& samples, const std::vector<float>& samplesSig) const
+ {
+ // Set the data and errors in the histogram object
+ //
+ for (size_t isample = 0; isample < _Nsample; isample++) {
+ _fitHist->SetBinContent(isample + 1, samples[isample]);
+ _fitHist->SetBinError(isample + 1, samplesSig[isample]);
+ }
+ }
+
+public:
+
+ ZDCPulseAnalyzer(std::string tag, int Nsample, float deltaTSample, size_t preSampleIdx, int pedestal, float gainHG,
+ std::string fitFunction, int peak2ndDerivMinSample, float peak2DerivMinThreshHG, float peak2DerivMinThreshLG);
+
+ ~ZDCPulseAnalyzer();
+
+ static void SetFitOPtions(std::string fitOptions) { _fitOptions = fitOptions;}
+ static void SetQuietFits(bool quiet) {_quietFits = quiet;}
+ static bool QuietFits() {return _quietFits;}
+
+ void SetForceLG(bool forceLG) {_forceLG = forceLG;}
+ bool ForceLG() const {return _forceLG;}
+
+ void SetCutValues(float chisqDivAmpCutHG, float chisqDivAmpCutLG,
+ float deltaT0MinHG, float deltaT0MaxHG,
+ float deltaT0MinLG, float deltaT0MaxLG) ;
+
+ void SetTauT0Values(bool fixTau1, bool fixTau2, float tau1, float tau2, float t0HG, float t0LG);
+
+ void SetADCOverUnderflowValues(int HGOverflowADC, int HGUnderflowADC, int LGOverflowADC);
+
+ void SetTimingCorrParams(std::vector<float> HGT0CorrParams, std::vector<float> LGT0CorrParams)
+ {
+ if (HGT0CorrParams.size() == 4 && LGT0CorrParams.size() == 4) {
+ _HGT0CorrParams = HGT0CorrParams;
+ _LGT0CorrParams = LGT0CorrParams;
+ _haveTimingCorrections = true;
+ }
+ else throw;
+ }
+
+ void SetNonlinCorrParams(const std::vector<float>& params)
+ {
+ // Check for valid length
+ //
+ if (params.size() != 2) throw;
+
+ std::cout << "Setting non-linear parameters for module: " << _tag << ", vlues = "
+ << params[0] << ", " << params[1] << std::endl;
+
+ _nonLinCorrParams = params;
+ }
+
+ bool LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG);
+
+ bool HaveData() const {return _haveData;}
+ bool HavePulse() const {return _havePulse;}
+ bool Failed() const {return _fail;}
+ bool BadChisq() const {return _badChisq;}
+ bool BadT0() const {return _badT0;}
+ bool FitFailed() const {return _fitFailed;}
+ bool PrePulse() const {return _prePulse;}
+ bool UseLowGain() const {return _useLowGain;}
+
+ bool HGOverflow() const {return _HGOverflow;}
+ bool HGUnderflow() const {return _HGOverflow;}
+ bool LGOverflow() const {return _LGOverflow;}
+ bool LGUnderflow() const {return _LGUnderflow;}
+
+ bool PSHGOverUnderflow() const {return _PSHGOverUnderflow;}
+
+ float GetFitAmplitude() const {return _fitAmplitude;}
+ float GetFitT0() const {return _fitT0;}
+ float GetT0Corr() const {return _fitT0Corr;}
+ float GetChisq() const {return _fitChisq;}
+
+ float GetFitTau1() const {return _fitTau1;}
+ float GetFitTau2() const {return _fitTau2;}
+
+ float GetAmplitude() const {return _amplitude;}
+ float GetAmpError() const {return _ampError;}
+
+ float GetPresample() const {return _preSample;}
+
+ float GetMaxADC() const {return _maxADCValue;}
+ float GetMinADC() const {return _minADCValue;}
+
+ int GetMaxADCSample() const {return _maxSampl;}
+ int GetMinADCSample() const {return _minSampl;}
+
+ float GetMaxDelta() const {return _maxDelta;}
+ float GetMinDelta() const {return _minDelta;}
+
+ float GetMinDeriv2nd() const {return _minDeriv2nd;}
+ float GetMinDeriv2ndIndex() const {return _minDeriv2ndIndex;}
+
+ unsigned int GetStatusMask() const;
+
+ float GetPreSampleAmp() const {return _preSampleAmp;}
+ float GetBkgdMaxFraction() const {return _bkgdMaxFraction;}
+
+ const TH1* GetHistogramPtr() const {
+ //
+ // We defer filling the histogram if we don't have a pulse until the histogram is requested
+ //
+ if (!_havePulse) {
+ if (UseLowGain()) FillHistogram(_samplesSub, _ADCSSampSigLG);
+ else FillHistogram(_samplesSub, _ADCSSampSigHG);
+ }
+
+ return _fitHist;
+ }
+
+ void DoFit();
+
+ void Dump() const;
+
+ const std::vector<float>& GetSamplesSub() const {return _samplesSub;}
+ const std::vector<float>& GetSamplesDeriv() const {return _samplesDeriv;}
+ const std::vector<float>& GetSamplesDeriv2nd() const {return _samplesDeriv2nd;}
+};
+
+
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCTriggerEfficiency.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCTriggerEfficiency.h
new file mode 100644
index 000000000000..0bc68d012e02
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCTriggerEfficiency.h
@@ -0,0 +1,101 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef _ZDCTriggerEfficiency_
+#define _ZDCTriggerEfficiency_
+
+#include <TSpline.h>
+#include <iostream>
+#include <vector>
+#include <array>
+#include <algorithm>
+using namespace std;
+class ZDCTriggerEfficiency{
+ public:
+ bool _haveParams;
+ bool _haveCorrCoeffs;
+
+ unsigned int _minLB;
+ unsigned int _maxLB;
+
+ unsigned int _currentLB;
+
+ std::array<std::vector<TSpline3*>, 2> _effParams;
+ std::array<std::vector<TSpline3*>, 2> _effParamErrors;
+ std::array<std::vector<TSpline3*>, 2> _effParamCorrCoeffs;
+
+ std::array<std::vector<double>, 2> _currentParams;
+ std::array<std::vector<double>, 2> _currentParamErrors;
+ std::array<std::vector<double>, 2> _currentCorrCoefff;
+
+ void UpdatelumiBlock(unsigned int lumiBlock)
+ {
+ if (!_haveParams) throw;
+
+ if (lumiBlock != _currentLB) {
+ int newLB = std::max(std::min(lumiBlock, _maxLB), _minLB);
+ for (int side : {0, 1}) {
+ _currentParams[side].clear();
+ _currentParamErrors[side].clear();
+ _currentCorrCoefff[side].clear();
+
+ for (size_t ipar = 0; ipar < _effParams[side].size(); ipar++) {
+ _currentParams[side].push_back(_effParams[side][ipar]->Eval(newLB));
+ _currentParamErrors[side].push_back(_effParamErrors[side][ipar]->Eval(newLB));
+
+ if (_haveCorrCoeffs) {
+ _currentCorrCoefff[side].push_back(_effParamCorrCoeffs[side][ipar]->Eval(newLB));
+ }
+ }
+ _currentLB = lumiBlock;
+ }
+ }
+ }
+
+public:
+ ZDCTriggerEfficiency() : // in order, alpha-beta, alpha-theta, beta-theta
+ _haveParams(false), _haveCorrCoeffs(false),
+ _minLB(0), _maxLB(0), _currentLB(0)
+ {
+
+ }
+
+ void SetEffParamsAndErrors(std::array<std::vector<TSpline3*>, 2> effParams,
+ std::array<std::vector<TSpline3*>, 2> effParamErrors)
+ {
+ for (int side : {0, 1}) {
+ _effParams[side] = effParams[side];
+ _effParamErrors[side] = effParamErrors[side];
+ for (int iarr=0;iarr<3;iarr++)
+ {
+ _effParams[side].at(iarr)->Print();
+ _effParamErrors[side].at(iarr)->Print();
+ }
+ }
+
+ _minLB = _effParams[0][0]->GetXmin();
+ _maxLB = _effParams[0][0]->GetXmax();
+ _haveParams = true;
+ }
+
+ void SetEffParamCorrCoeffs(std::array<std::vector<TSpline3*>, 2> effParamsCorrCoeffs)
+ {
+ for (int side : {0, 1}) {
+ _effParamCorrCoeffs[side] = effParamsCorrCoeffs[side];
+ for (int iarr=0;iarr<3;iarr++)
+ {
+ _effParamCorrCoeffs[side].at(iarr)->Print();
+ }
+ }
+
+ _haveCorrCoeffs = true;
+ }
+
+ float GetEfficiency(int side, float ADCSum);
+
+ std::pair<float, float> GetEfficiencyAndError(int side, float ADCSum);
+
+};
+
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
new file mode 100644
index 000000000000..278e115d3683
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
@@ -0,0 +1,139 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef __ZDCANALYSISTOOL_H__
+#define __ZDCANALYSISTOOL_H__
+
+#include "AsgTools/AsgTool.h"
+#include "xAODForward/ZdcModuleContainer.h"
+#include "xAODTrigL1Calo/TriggerTowerContainer.h"
+
+#include "ZdcAnalysis/ZDCDataAnalyzer.h"
+#include "ZdcAnalysis/ZDCTriggerEfficiency.h"
+#include "ZdcAnalysis/IZdcAnalysisTool.h"
+
+#include "TF1.h"
+#include "TMath.h"
+
+namespace ZDC
+{
+
+class ZdcAnalysisTool : public virtual IZdcAnalysisTool, public asg::AsgTool
+{
+
+ ASG_TOOL_CLASS(ZdcAnalysisTool, ZDC::IZdcAnalysisTool)
+
+ public:
+ ZdcAnalysisTool(const std::string& name);
+ virtual ~ZdcAnalysisTool();
+
+ //interface from AsgTool
+ StatusCode initializeTool();
+ StatusCode initialize() {return initializeTool();}
+ void initialize80MHz();
+ void initialize40MHz();
+ void initializeTriggerEffs(unsigned int runNumber);
+
+ StatusCode recoZdcModule(const xAOD::ZdcModule& module);
+ StatusCode recoZdcModules(const xAOD::ZdcModuleContainer& moduleContainer);
+ StatusCode reprocessZdc();
+
+ // methods for processing, used for decoration
+ bool sigprocMaxFinder(const std::vector<unsigned short>& adc, float deltaT, float& amp, float& time, float& qual);
+ bool sigprocSincInterp(const std::vector<unsigned short>& adc, float deltaT, float& amp, float& time, float& qual);
+
+ void setEnergyCalibrations(unsigned int runNumber);
+ void setTimeCalibrations(unsigned int runNumber);
+
+ float getModuleSum(int side);
+
+ //float getModuleFitAmplitude(int side, int imod);
+ //float getModuleFitT0(int side, int imod);
+ //float getModuleAmplitude(int side, int imod);
+ //float getModuleTime(int side, int imod);
+ //float getModuleChisq(int side, int imod);
+ //float getModuleStatus(int side, int imod);
+
+ //float getModuleCalibAmplitude(int side, int imod);
+ //float getModuleCalibTime(int side, int imod);
+
+ float getCalibModuleSum(int side);
+ float getAverageTime(int side);
+ bool sideFailed(int side);
+ unsigned int getModuleMask();
+
+ float getTriggerEfficiency(int side);
+
+ const ZDCDataAnalyzer* getDataAnalyzer() {return m_zdcDataAnalyzer;}
+
+ private:
+ // Private methods
+ //
+ ZDCDataAnalyzer* initializeDefault();
+
+ StatusCode configureNewRun(unsigned int runNumber);
+
+ // Data members
+ //
+ std::string m_name;
+ bool m_init;
+ std::string m_configuration;
+ std::string m_zdcAnalysisConfigPath;
+ std::string m_zdcEnergyCalibFileName;
+ std::string m_zdcTimeCalibFileName;
+ std::string m_zdcTriggerEffParamsFileName;
+
+ bool m_writeAux;
+ std::string m_auxSuffix;
+
+ bool m_eventReady;
+ unsigned int m_runNumber;
+ unsigned int m_lumiBlock;
+
+ // internal functions
+ TF1* tf1SincInterp;
+
+ std::string m_zdcModuleContainerName;
+ const xAOD::ZdcModuleContainer* m_zdcModules;
+ bool m_flipEMDelay;
+ bool m_lowGainOnly;
+ bool m_doCalib;
+ int m_forceCalibRun;
+ int m_forceCalibLB;
+
+ // Parameters that control the pulse fitting analysis
+ //
+ unsigned int m_numSample;
+ float m_deltaTSample;
+ unsigned int m_presample;
+ unsigned int m_peakSample;
+ float m_Peak2ndDerivThresh;
+ float m_t0;
+ float m_tau1;
+ float m_tau2;
+ bool m_fixTau1;
+ bool m_fixTau2;
+ float m_deltaTCut;
+ float m_ChisqRatioCut;
+
+ std::array<std::array<TSpline*, 4>, 2> m_splines;
+
+ ZDCDataAnalyzer* m_zdcDataAnalyzer;
+ ZDCDataAnalyzer* m_zdcDataAnalyzer_40MHz;
+ ZDCDataAnalyzer* m_zdcDataAnalyzer_80MHz;
+ ZDCDataAnalyzer::ZDCModuleFloatArray _peak2ndDerivMinSamples;
+ ZDCDataAnalyzer::ZDCModuleFloatArray _peak2ndDerivMinThresholdsHG;
+ ZDCDataAnalyzer::ZDCModuleFloatArray _peak2ndDerivMinThresholdsLG;
+
+
+ ZDCTriggerEfficiency* m_zdcTriggerEfficiency;
+
+};
+
+} // namespace ZDC
+
+#endif
+
+
+
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcSincInterp.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcSincInterp.h
new file mode 100644
index 000000000000..ec4f60a5d69a
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcSincInterp.h
@@ -0,0 +1,12 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef _ZDC_SINC_INTERP_
+#define _ZDC_SINC_INTERP_
+
+namespace ZDC
+{
+double SincInterp(double* xvec, double* pvec);
+}
+#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
new file mode 100644
index 000000000000..4b26464a70b3
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
@@ -0,0 +1,60 @@
+# this makefile also gets parsed by shell scripts
+# therefore it does not support full make syntax and features
+# edit with care
+
+# for full documentation check:
+# https://twiki.cern.ch/twiki/bin/viewauth/Atlas/RootCore#Package_Makefile
+
+
+# the name of the package:
+PACKAGE = ZdcAnalysis
+
+# the libraries to link with this one:
+PACKAGE_PRELOAD =
+
+# additional compilation flags to pass (not propagated to dependent packages):
+PACKAGE_CXXFLAGS =
+
+# additional compilation flags to pass (propagated to dependent packages):
+PACKAGE_OBJFLAGS =
+
+# additional linker flags to pass (for compiling the library):
+PACKAGE_LDFLAGS =
+
+# additional linker flags to pass (for compiling binaries):
+PACKAGE_BINFLAGS =
+
+# additional linker flags to pass (propagated to client libraries):
+PACKAGE_LIBFLAGS =
+
+# the list of packages we depend on:
+PACKAGE_DEP = xAODRootAccess xAODBase xAODForward xAODEventInfo AsgTools xAODTrigL1Calo PathResolver
+
+# the list of packages we use if present, but that we can work without :
+PACKAGE_TRYDEP =
+
+# list pattern of scripts to link directly into binary path:
+PACKAGE_SCRIPTS =
+
+# whether to use pedantic compilation:
+PACKAGE_PEDANTIC = 1
+
+# whether to turn *off* optimization (set to dict to do it only for
+# dictionaries):
+PACKAGE_NOOPT = 0
+
+# whether to build no library (needs to be set if no source files are
+# present):
+PACKAGE_NOCC = 0
+
+# whether we build a reflex dictionary:
+PACKAGE_REFLEX = 0
+
+# the list of all unit tests that should be called in recursive testing,
+# i.e. in unit tests that call other unit tests
+# for that unit tests need to pass on all machines, and run very fast
+PACKAGE_RECURSIVE_UT =
+
+
+
+include $(ROOTCOREDIR)/Makefile-common
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/cmt/requirements b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/requirements
new file mode 100644
index 000000000000..937d3d8020ad
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/requirements
@@ -0,0 +1,24 @@
+package ZdcAnalysis
+
+author Peter Steinberg <peter.steinberg@bnl.gov>
+
+use AtlasPolicy AtlasPolicy-*
+
+use AsgTools AsgTools-* Control/AthToolSupport
+use AtlasROOT AtlasROOT-* External
+use GaudiInterface GaudiInterface-* External
+use xAODForward xAODForward-* Event/xAOD
+use xAODTrigL1Calo xAODTrigL1Calo-* Event/xAOD
+
+private
+use xAODEventInfo xAODEventInfo-* Event/xAOD
+use PathResolver PathResolver-* Tools
+end_private
+
+library ZdcAnalysis ../Root/*.cxx
+apply_pattern component_library
+
+apply_pattern declare_python_modules files="*.py"
+apply_pattern declare_joboptions files="*.py"
+
+
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/data/ZdcAnalysisConfig.conf b/ForwardDetectors/ZDC/ZdcAnalysis/data/ZdcAnalysisConfig.conf
new file mode 100644
index 000000000000..2378e74959fc
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/data/ZdcAnalysisConfig.conf
@@ -0,0 +1,3 @@
+ZdcEnergyCalibFileName ZdcCalibration_comb_v4.root
+ZdcTriggerEffFileName ZdcTriggerEffParameters_v3.root
+
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/tools/RunZDCTreeAnalysis.C b/ForwardDetectors/ZDC/ZdcAnalysis/tools/RunZDCTreeAnalysis.C
new file mode 100644
index 000000000000..fc9c53e45090
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/tools/RunZDCTreeAnalysis.C
@@ -0,0 +1,216 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <string>
+#include <iostream>
+#include <sstream>
+
+void SetT0Values(double LGt0Array[2][4], double HGt0Array[2][4])
+{
+
+ HGt0Array[0][0] = 44.0 + 0.24;
+ HGt0Array[0][1] = 40.7 - 0.34;
+ HGt0Array[0][2] = 48.2 + 0.13;
+ HGt0Array[0][3] = 35.5 + 0.53;
+
+ HGt0Array[1][0] = 35.2 + 0.04;
+ HGt0Array[1][1] = 31.3 - 0.22;
+ HGt0Array[1][2] = 40.9 - 0.15;
+ HGt0Array[1][3] = 30.7 - 0.17;
+
+ LGt0Array[0][0] = 42.5 + 0.15;
+ LGt0Array[0][1] = 39.0 - 0.52 + 0.027;
+ LGt0Array[0][2] = 46.3 - 0.06 - 0.06;
+ LGt0Array[0][3] = 33.7 - 0.37 - 0.10;
+
+ LGt0Array[1][0] = 33.4 + 0.33 - 0.66;
+ LGt0Array[1][1] = 30.1 - 0.21 - 0.04;
+ LGt0Array[1][2] = 39.5 - 0.37 - 0.05;
+ LGt0Array[1][3] = 29.6 - 0.28 - 0.058;
+
+}
+
+void LoadCalibrations(ZDCTreeAnalysis* ana, std::string filename, int runNumber, bool eCalib = true, bool t0Calib = true)
+{
+ TFile* file = TFile::Open(filename.c_str());
+ if (!file->IsOpen()) {
+ std::cout << "Error opening calibration file, quitting" << std::endl;
+ }
+
+ std::ostringstream runString;
+ runString << "run" << runNumber;
+
+ if (eCalib) {
+
+ std::string calibNamePrefix = "ZDC_Ecalib_" + runString.str();
+
+ std::array<std::array<TSpline*, 4>, 2> ECalibSplinePtrs;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ std::ostringstream sideModStr;
+ sideModStr << "_s" << side << "_m" << module;
+
+ TSpline3* spline_p = (TSpline3*) file->GetObjectChecked((calibNamePrefix + sideModStr.str()).c_str(), "TSpline3");
+
+ ECalibSplinePtrs[side][module] = spline_p;
+ }
+ }
+
+ ana->LoadEnergyCalibrations(ECalibSplinePtrs);
+ }
+
+ if (t0Calib) {
+ std::string calibNamePrefix = "ZDC_T0calib_" + runString.str();
+
+ std::array<std::array<TSpline*, 4>, 2> t0CalibSplinePtrsHG, t0CalibSplinePtrsLG;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ std::ostringstream sideModStr;
+ sideModStr << "_s" << side << "_m" << module;
+
+ std::string HGName = calibNamePrefix + "_HG" + sideModStr.str();
+ std::string LGName = calibNamePrefix + "_LG" + sideModStr.str();
+
+ TSpline3* HGSpline_p = (TSpline3*) file->GetObjectChecked(HGName.c_str(), "TSpline3");
+ TSpline3* LGSpline_p = (TSpline3*) file->GetObjectChecked(LGName.c_str(), "TSpline3");
+
+ if (!HGSpline_p || !LGSpline_p) {
+ std::cout << "Error reading t0 calibrations for side, mod = " << side << ", " << module << " -- quitting" << std::endl;
+ return;
+ }
+
+ t0CalibSplinePtrsHG[side][module] = HGSpline_p;
+ t0CalibSplinePtrsLG[side][module] = LGSpline_p;
+ }
+ }
+
+ ana->LoadT0Calibrations(t0CalibSplinePtrsHG, t0CalibSplinePtrsLG);
+ }
+
+ file->Close();
+}
+
+ZDCTreeAnalysis* InitZDCAnalysis(std::string inputFile, std::string calibFile, bool fixTau1, bool fixTau2,
+ bool eCalib = true, bool t0Calib = true, bool doSlew = false)
+{
+
+ ZDCTreeAnalysis* ana = new ZDCTreeAnalysis(inputFile, 7, 12.5, 0);
+
+ ana->SetDebugLevel(0);
+
+ // ADC over/underflow thresholds
+ //
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGOverFlowADC = {800, 800, 800, 800, 800, 800, 800, 800};
+ ZDCDataAnalyzer::ZDCModuleFloatArray HGUnderFlowADC = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray LGOverFlowADC = {950, 950, 950, 950, 950, 950, 950, 950};
+
+ ana->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+
+ // Set Tau and nominal timing offsets
+ //
+ ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr;
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ fixTau1Arr[side][module] = fixTau1;
+ fixTau2Arr[side][module] = fixTau2;
+ }
+ }
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau1 = {3.9, 3.4, 4.1, 4.2,
+ 4.2, 3.6, 3.3, 3.4};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray tau2 = {20.0, 20.4, 18.9, 20.8,
+ 19.1, 21.9, 22.6, 23.4};
+
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0HG = {44.0 + 0.24, 40.7 - 0.34, 48.2 + 0.13, 35.5 + 0.53,
+ 35.2 + 0.04, 31.3 - 0.22, 40.9 - 0.15, 30.7 - 0.17};
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray t0LG = {42.5 + 0.15, 39.0 - 0.52 + 0.027, 46.3 - 0.06 - 0.06, 33.7 - 0.37 - 0.10,
+ 33.4 + 0.33 - 0.66, 30.1 - 0.21 - 0.04, 39.5 - 0.37 - 0.05, 29.6 - 0.28 - 0.058};
+
+
+ ana->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0HG, t0LG);
+
+ int runNumber = ana->GetRunNumber();
+
+ if (doSlew) {
+ // float HGSlewCoeff[2][4][3] = {-6.5e-2, 2.85e-2, -2.83e-3,
+ // -5.5e-2, 5.13e-2, 5.6e-3,
+ // -1.45e-3, 9.3e-2, 3.9e-3,
+ // -2.36e-2, 8.3e-2, 1.1e-3,
+ // -6.5e-2, 4.84e-2, -3.7e-3,
+ // 1.34e-2, 6.57e-2, 5.37e-3,
+ // -5.37e-2, 3.49e-2, 3.8e-3,
+ // -3.3e-2,3.9e-2,2.2e-3};
+
+ // float LGSlewCoeff[2][4][3] = {-9.6e-2, 4.39e-2, 2.93e-3, //[0][0]
+ // -5.0e-2, 14.9e-2, 20.6e-3, //[0][1]
+ // -4.4e-2, 5.3e-2, 0, //[0][2]
+ // -9.90e-2, 4.08e-2, 0, //[0][3]
+ // -8.7e-2, 4.2e-2, -3.2e-3, //[1][0]
+ // -3.26e-2, 3.84e-2, -2.32e-3, //[1][1]
+ // -26.8e-2, -2.64e-2, -5.3e-3,//[1][2]
+ // -13.2e-2, 0.45e-2, -2.4e-3}; //[1][3]
+ //
+ //
+ std::array<std::array<std::vector<float>, 4>, 2> slewingParamsHG, slewingParamsLG;
+
+ slewingParamsHG[0][0] = {0, -6.5e-2, 2.85e-2, -2.83e-3};
+ slewingParamsHG[0][1] = {0, -5.5e-2, 5.13e-2, 5.6e-3};
+ slewingParamsHG[0][2] = {0, -1.45e-3, 9.3e-2, 3.9e-3};
+ slewingParamsHG[0][3] = {0, -2.36e-2, 8.3e-2, 1.1e-3};
+
+ slewingParamsHG[1][0] = {0, -6.5e-2, 4.84e-2, -3.7e-3};
+ slewingParamsHG[1][1] = {0, 1.34e-2, 6.57e-2, 5.37e-3};
+ slewingParamsHG[1][2] = {0, -5.37e-2, 3.49e-2, 3.8e-3};
+ slewingParamsHG[1][3] = {0, -3.3e-2, 3.9e-2, 2.2e-3};
+
+ slewingParamsLG[0][0] = {0, -9.6e-2, 4.39e-2, 2.93e-3 };
+ slewingParamsLG[0][1] = {0, -5.0e-2, 14.9e-2, 20.6e-3 };
+ slewingParamsLG[0][2] = {0, -4.4e-2, 5.3e-2, 0, };
+ slewingParamsLG[0][3] = {0, -9.90e-2, 4.08e-2, 0, };
+
+ slewingParamsLG[1][0] = {0, -8.7e-2, 4.2e-2, -3.2e-3 };
+ slewingParamsLG[1][1] = {0, -3.26e-2, 3.84e-2, -2.32e-3};
+ slewingParamsLG[1][2] = {0, -26.8e-2, -2.64e-2, -5.3e-3 };
+ slewingParamsLG[1][3] = {0, -13.2e-2, 0.45e-2, -2.4e-3 };
+
+ ana->SetSlewingCoeff(slewingParamsHG, slewingParamsLG);
+ }
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutHG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray chisqDivAmpCutLG = {10, 10, 10, 10, 10, 10, 10, 10};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowHG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighHG = {8, 8, 8, 11, 8, 10, 8, 12};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutLowLG = {-6, -5, -5, -5, -5, -5, -5, -5};
+ ZDCDataAnalyzer::ZDCModuleFloatArray DeltaT0CutHighLG = {8, 8, 8, 11, 8, 10, 8, 12};
+
+ ana->SetCutValues(chisqDivAmpCutHG, chisqDivAmpCutLG, DeltaT0CutLowHG, DeltaT0CutHighHG, DeltaT0CutLowLG, DeltaT0CutHighLG);
+
+
+ // float moduleECalib[2][4] = {1., 1, 1, 1, 1., 1, 1, 1};
+
+ // if (eCalib) ana->EnableECalib(moduleECalib);
+
+ if (eCalib || t0Calib) LoadCalibrations(ana, calibFile, runNumber, eCalib, t0Calib);
+
+ return ana;
+}
+
+void RunZDCAnalysis(std::string inputFile, std::string outputFile, std::string calibFile, int nevent = -1,
+ bool fixTau1 = true, bool fixTau2 = true, bool eCalib = false, bool t0Calib = false, bool doSlew = false)
+{
+
+ ZDCTreeAnalysis* ana = InitZDCAnalysis(inputFile, calibFile, fixTau1, fixTau2, eCalib, t0Calib, doSlew);
+
+ ana->OpenOutputTree(outputFile);
+ ana->Loop(nevent);
+ ana->CloseOutputTree();
+}
+
+
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.cxx
new file mode 100644
index 000000000000..b251b010fca1
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.cxx
@@ -0,0 +1,603 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#define ZDCNLCalibration_cxx
+#include "ZDCNLCalibration.h"
+
+#include <TF1.h>
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+#include <iomanip>
+#include <numeric>
+
+void ZDCNLCalibration::FillLumiBlockEvtMap()
+{
+ Long64_t numEntries = m_tree->GetEntriesFast();
+
+ unsigned int countPassed = 0;
+
+ int lastLumiBlock = -1;
+ int LBStartEntry = -1;
+
+ for (Long64_t entry = 0; entry < numEntries; entry++) {
+ Long64_t numBytes = m_tree->GetEntry(entry);
+
+ if (m_useGRL && passBits !=0) {
+ //std::cout << "lumiBlock = " << lumiBlock << " passBits = " << passBits << std::endl;
+ continue; // reject based on GRL or errors
+ }
+
+ if (lumiBlock != lastLumiBlock) {
+ if (lastLumiBlock > -1) {
+ m_LumiBlockEvtMap.insert(LBEvtMap::value_type(lumiBlock, std::pair<unsigned int, unsigned int>(LBStartEntry, entry - 1)));
+ }
+
+ LBStartEntry = entry;
+ lastLumiBlock = lumiBlock;
+ }
+
+ countPassed++;
+ }
+
+ if (m_debugLevel >= 0) {
+ std::cout << "From a total of " << numEntries << " entries, " << countPassed << " were selected for analysis" << std::endl;
+ }
+}
+
+void ZDCNLCalibration::AddCalibration(size_t side, std::string tag, const CalibData& calib)
+{
+ std::map<std::string, CalibData>::iterator iter = m_calibrations.at(side).find(tag);
+ if (iter != m_calibrations[side].end()) (*iter).second = calib;//*iter = std::pair<std::string, CalibData>(tag, calib);
+ else {
+ m_calibrations[side].insert(std::pair<std::string, CalibData>(tag, calib));
+ }
+}
+
+CalibData ZDCNLCalibration::GetCalibration(size_t side, std::string tag)
+{
+ CalibData null;
+
+ auto iter = m_calibrations.at(side).find(tag);
+ if (iter != m_calibrations.at(side).end()) {
+ return iter->second;
+ }
+ else return null;
+}
+
+std::pair<float, float> ZDCNLCalibration::FindSNRange(size_t LBLow, size_t LBHigh, size_t side)
+{
+ static TH1D* snHist = new TH1D("ZDCNLCalibSNHist", "", 200, 0, 2000);
+ snHist->Reset();
+
+ LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
+ LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
+
+ // Select which trigger we're going to look at based on the side
+ //
+ const bool& trigger = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
+
+ for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; iter++) {
+ unsigned int entryStart = iter->second.first;
+ unsigned int entryEnd = iter->second.second;
+
+ for (unsigned int entry = entryStart; entry < entryEnd; entry++) {
+ Long64_t nb = m_tree->GetEntry(entry);
+ if (nb < 0) {
+ std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
+ throw;
+ }
+
+ if (!trigger) continue;
+
+ int testMask = zdc_ZdcModuleMask >> side*4;
+ if ((testMask & 0xf) != 0xf) continue;
+
+ float sumAmp = zdc_ZdcAmp[side];
+ snHist->Fill(sumAmp);
+ }
+ }
+
+ int ibin = snHist->GetMaximumBin();
+ double center = snHist->GetBinCenter(ibin);
+ double maxCounts = snHist->GetBinContent(ibin);
+
+ static TF1* gaussFunc = new TF1("gaussFunc", "gaus", 100, 2000);
+
+ gaussFunc->SetParameters(maxCounts, center, 100);
+ snHist->Fit("gaussFunc", "", "", center*0.7, center*1.3);
+
+ double gaussMean = gaussFunc->GetParameter(1);
+ double gaussWidth = gaussFunc->GetParameter(2);
+
+ double cutLow = gaussMean -2.*gaussWidth;
+ double cutHigh = gaussMean + 2*gaussWidth;
+
+ if (m_debugLevel >= 0) {
+ std::cout << "SN cut range = " << cutLow << " to " << cutHigh << std::endl;
+ }
+
+ return std::pair<float, float>(cutLow, cutHigh);
+}
+
+std::pair<std::pair<float, float>,std::pair<float, float> >
+ZDCNLCalibration::FindSNTwoNRanges(size_t LBLow, size_t LBHigh, size_t side)
+{
+ static TH1D* sntnHist = new TH1D("ZDCNLCalibSNTNHist", "", 200, 0, 2000);
+ sntnHist->Reset();
+
+ LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
+ LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
+
+ // Select which trigger we're going to look at based on the side
+ //
+ const bool& trigger = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
+
+ for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; iter++) {
+ unsigned int entryStart = iter->second.first;
+ unsigned int entryEnd = iter->second.second;
+
+ for (unsigned int entry = entryStart; entry < entryEnd; entry++) {
+ Long64_t nb = m_tree->GetEntry(entry);
+ if (nb < 0) {
+ std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
+ throw;
+ }
+
+ if (!trigger) continue;
+
+ int testMask = zdc_ZdcModuleMask >> side*4;
+ if ((testMask & 0xf) != 0xf) continue;
+
+ float sumAmp = zdc_ZdcAmp[side];
+ sntnHist->Fill(sumAmp);
+ }
+ }
+
+ int ibin = sntnHist->GetMaximumBin();
+ double center = sntnHist->GetBinCenter(ibin);
+ double maxCounts = sntnHist->GetBinContent(ibin);
+
+ static TF1* twogaussFunc = new TF1("twogaussFunc", "[0]*exp(-pow((x-[1])/[2],2)/2.) + [3]*exp(-pow((x-[4])/[5],2)/2.)", 100, 2000);
+
+ twogaussFunc->SetParameters(maxCounts, center, 100, maxCounts/2, 2*center, 200);
+ sntnHist->Fit("twogaussFunc", "", "", center*0.7, 2*center*1.3);
+
+ double gaussMean1 = twogaussFunc->GetParameter(1);
+ double gaussMean2 = twogaussFunc->GetParameter(4);
+
+ double gaussWidth1 = std::abs(twogaussFunc->GetParameter(2));
+ double gaussWidth2 = std::abs(twogaussFunc->GetParameter(5));
+
+ double cutLow1 = gaussMean1 -2.*gaussWidth1;
+ double cutHigh1 = gaussMean1 + 2*gaussWidth1
+;
+ double cutLow2 = gaussMean2 -2.*gaussWidth2;
+ double cutHigh2 = gaussMean2 + 2*gaussWidth2;
+
+ if (cutHigh1 > cutLow2) {
+ double middle = 0.5*(cutHigh1 + cutLow2);
+ cutHigh1 = middle;
+ cutLow2 = middle;
+ }
+
+ if (m_debugLevel >= 0) {
+ std::cout << "SN cut range = " << cutLow1 << " to " << cutHigh1 << std::endl;
+ std::cout << "TN cut range = " << cutLow2 << " to " << cutHigh2 << std::endl;
+ }
+
+ return std::pair<std::pair<float, float>, std::pair<float, float> >(std::pair<float, float>(cutLow1, cutHigh1),
+ std::pair<float, float>(cutLow2, cutHigh2));
+}
+
+void ZDCNLCalibration::Calibrate(size_t side, std::string calibInput, std::string calibOutput,
+ size_t LBLow, size_t LBHigh, std::array<int, 4> maxPowerModule,
+ std::vector<std::pair<double, double> > nNeutERange,
+ bool excludeHE, float heSumThresh, float HEDeweight)
+{
+ std::for_each(maxPowerModule.begin(), maxPowerModule.end(),
+ [](int power) {
+ if (power == 0) {
+ std::cout << "Found max power = 0, quitting" << std::endl;
+ return;
+ }
+ } );
+
+ // Set up container for sums that we will use to calcualting the weights
+ //
+ size_t numWeights = std::accumulate(maxPowerModule.begin(), maxPowerModule.end(), 0);
+ size_t arrayDim = 4*m_maxNLPower;
+
+ size_t arrayDim2D = arrayDim*arrayDim;
+
+ size_t numNeutMax = nNeutERange.size();
+
+ std::vector<double> sumsHE(arrayDim, 0);
+ std::vector<double> sums2DHE(arrayDim2D, 0);
+
+ std::vector<std::vector<double> > sumsNeutVec(numNeutMax, std::vector<double>(arrayDim, 0));
+ std::vector<std::vector<double> > sums2DNeutVec(numNeutMax, std::vector<double>(arrayDim2D, 0));
+
+ LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
+ LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
+
+ // Select which trigger we're going to look at based on the side
+ //
+ const bool& triggerOpp = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
+ const bool& triggerSame = (side == 1 ? L1_ZDC_A : L1_ZDC_C);
+
+ std::vector<int> numNeutEvent(numNeutMax, 0);
+ int numHEEvent = 0;
+
+ //
+ //
+ std::string calibName = (calibInput != "" ? calibInput : "default");
+ CalibData calib = GetCalibration(side, calibName);
+
+ for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; iter++) {
+ unsigned int entryStart = iter->second.first;
+ unsigned int entryEnd = iter->second.second;
+
+ for (unsigned int entry = entryStart; entry < entryEnd; entry++) {
+ Long64_t nb = m_tree->GetEntry(entry);
+ if (nb < 0) {
+ std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
+ throw;
+ }
+
+ int testMask = zdc_ZdcModuleMask >> side*4;
+ if ((testMask & 0xf) != 0xf) continue;
+
+ double sumAmp = CalculateEnergy(zdc_ZdcModuleAmp[side], calib);
+
+ if (triggerSame && sumAmp > heSumThresh && !excludeHE) {
+ AddToSums(sumsHE, sums2DHE, zdc_ZdcModuleAmp[side]);
+ numHEEvent++;
+ }
+
+ // For 1N we require the opposite side trigger
+ //
+ if (triggerOpp && sumAmp > nNeutERange[0].first && sumAmp < nNeutERange[0].second) {
+ AddToSums(sumsNeutVec[0], sums2DNeutVec[0], zdc_ZdcModuleAmp[side]);
+ numNeutEvent[0]++;
+ }
+
+ for (size_t iNNeut = 1; iNNeut < numNeutMax; iNNeut++) {
+ if (sumAmp > nNeutERange[iNNeut].first && sumAmp < nNeutERange[iNNeut].second) {
+ AddToSums(sumsNeutVec[iNNeut], sums2DNeutVec[iNNeut], zdc_ZdcModuleAmp[side]);
+ numNeutEvent[iNNeut]++;
+ }
+ }
+ }
+ }
+
+ if (m_debugLevel >= 0) {
+ std::cout << "Statistics for side " << side << ", lb range = " << LBLow << ", " << LBHigh
+ << " : # high energy events = " << numHEEvent;
+
+ for (size_t iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
+ std::cout << ", # of " + std::to_string(iNNeut) + " N events = " << numNeutEvent[iNNeut];
+ }
+
+ std::cout << std::endl;
+ }
+
+
+ // Normalize the averages
+ //
+ for (size_t iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
+ if (numNeutEvent[iNNeut] == 0) {
+ std::cout << "Found no " + std::to_string(iNNeut) + "N events, quitting" << std::endl;
+ return;
+ }
+
+ std::for_each(sumsNeutVec[iNNeut].begin(), sumsNeutVec[iNNeut].end(), [&](double& sum) {sum /= numNeutEvent[iNNeut];});
+ std::for_each(sums2DNeutVec[iNNeut].begin(), sums2DNeutVec[iNNeut].end(), [&](double& sum) {sum /= numNeutEvent[iNNeut];});
+ }
+
+ if (!excludeHE) {
+ if (numHEEvent == 0) {
+ std::cout << "Found no high-energy events, quitting" << std::endl;
+ return;
+ }
+
+ std::for_each(sumsHE.begin(), sumsHE.end(), [&](double& sum) {sum /= numHEEvent;});
+ std::for_each(sums2DHE.begin(), sums2DHE.end(), [&](double& sum) {sum /= numHEEvent;});
+ }
+
+
+ // Now we have the averages, set up to do the minimization
+ //
+ size_t numFitParams = numWeights;
+ TMatrixD minimMatrix(numFitParams, numFitParams);
+ TVectorD minimVector(numFitParams);
+
+ FillMinimizationData(minimMatrix, minimVector, maxPowerModule, HEDeweight,
+ sumsNeutVec, sumsHE, sums2DNeutVec, sums2DHE);
+
+ if (m_debugLevel >= 2) {
+ std::cout << "Dumping matrix equation: " << std::endl;
+
+ for (size_t index1 = 0; index1 < numFitParams; index1++) {
+ std::cout << "| ";
+
+ for (size_t index2 = 0; index2 < numFitParams; index2++) {
+ std::cout << std::scientific << std::setw(9) << std::setprecision(2) << minimMatrix(index1, index2) << " ";
+ }
+
+ std::cout << "| | w_" << index1 << " = | " << minimVector(index1) << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+
+ TDecompLU lu(minimMatrix);
+ lu.Decompose();
+
+ bool status = lu.Solve(minimVector);
+
+ if (!status) {
+ std::cout << "Minimization failed for side " << side << ", lb range "
+ << LBLow << ", " << LBHigh << endl;
+ }
+ else {
+ CalibData calibData;
+ calibData.LBStart = LBLow;
+ calibData.LBEnd = LBHigh;
+
+ int weightIndex = 0;
+ for (size_t module : {0, 1, 2, 3}) {
+ calibData.weights[module].assign(maxPowerModule[module], 0);
+
+ for (size_t power = 0; power < maxPowerModule[module]; power++) {
+ calibData.weights.at(module).at(power) = minimVector(weightIndex++);
+
+ if (m_debugLevel >= 0) {
+ std::cout << "Weight for module,power " << module << ", " << power << " = " << calibData.weights.at(module).at(power)
+ << std::endl;
+ }
+ }
+ }
+
+ // Now store this calibration
+ //
+ AddCalibration(side, calibOutput, calibData);
+ }
+}
+
+void ZDCNLCalibration::TestCalibration(int side, std::string name)
+{
+ std::string calibName = (name != "" ? name : "default");
+ CalibData calib = GetCalibration(side, calibName);
+
+ if (m_debugLevel >= 0) {
+ std::cout << "Testing calibration for side " << side << "name: " << calibName << std::endl;
+ }
+
+ m_testCalibSNHist = new TH1D("testCalibSNHist" , "", 4000, 0, 200000);
+ for (int module : {0, 1, 2, 3}) {
+ std::string fracHistName = "testCalibHEFracHist_" + std::to_string(module);
+ std::string energyHistName = "testCalibEnergyHist_" + std::to_string(module);
+
+ m_testCalibHEFracHist[module] = new TH1D(fracHistName.c_str(), "", 200, 0., 1.);
+ m_testCalibEnergyHist[module] = new TH1D(energyHistName.c_str(), "", 1000, 0., 100000);
+ }
+
+ static std::array<float, 4> treeRawAmp;
+ static std::array<float, 4> treeCalibAmp;
+ static float treeRawSum, treeCalibSum, treeRawSumOpp;
+ static unsigned int treeEventNumber;
+ static bool passedSN, passedTN, passedHE;
+
+ if (m_testTree != 0) delete m_testTree;
+
+ m_testTree = new TTree("ZDCNLCalibTest", "ZDCNLCalibTest");
+ m_testTree->SetDirectory(0);
+
+ m_testTree->Branch("EventNumber", &treeEventNumber, "eventNumber/I");
+ m_testTree->Branch("lumiBlock", &lumiBlock, "lumiBlock/I");
+ m_testTree->Branch("L1_ZDC_A", &L1_ZDC_A, "L1_ZDC_A/b");
+ m_testTree->Branch("L1_ZDC_C", &L1_ZDC_C, "L1_ZDC_C/b");
+ m_testTree->Branch("L1_ZDC_A_C", &L1_ZDC_A_C, "L1_ZDC_A_C/b");
+
+ m_testTree->Branch("RawSum", &treeRawSum, "RawSum/f");
+ m_testTree->Branch("CalibSum", &treeCalibSum, "CalibSum/f");
+ m_testTree->Branch("RawAmp", &treeRawAmp[0], "RawAmp[4]/f");
+ m_testTree->Branch("CalibAmp", &treeCalibAmp[0], "CalibAmp[4]/f");
+ m_testTree->Branch("RawSumOpp", &treeRawSumOpp, "RawSumOpp/f");
+
+ LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(calib.LBStart);
+ LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(calib.LBEnd);
+
+ int count = 0;
+ for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; iter++) {
+ unsigned int entryStart = iter->second.first;
+ unsigned int entryEnd = iter->second.second;
+
+ for (unsigned int entry = entryStart; entry < entryEnd; entry++) {
+ Long64_t nb = m_tree->GetEntry(entry);
+ if (nb < 0) {
+ std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
+ throw;
+ }
+
+ int testMask = zdc_ZdcModuleMask >> side*4;
+ if ((testMask & 0xf) != 0xf) continue;
+
+ double sumAmp = zdc_ZdcAmp[side];
+
+ // Calculate the calibrated energy
+ //
+ float sumCalibAmp = 0;
+ float weightIndex = 0;
+
+ std::array<double, 4> calibAmpModule = {0, 0, 0, 0};
+
+ if (m_debugLevel > 2 && count < 20) {
+ std::cout << "TestCalibration:: dumping calibration results for event " << count << std::endl;
+ }
+
+ for (size_t module : {0, 1, 2, 3}) {
+ float amp = zdc_ZdcModuleAmp[side][module];
+ float ampPow = amp;
+
+ for (size_t power = 0; power < calib.weights[module].size(); power++) {
+ calibAmpModule[module] += calib.weights[module][power]*ampPow;
+ ampPow *= amp;;
+ }
+
+ double energy = calibAmpModule[module];
+ sumCalibAmp += energy;
+ m_testCalibEnergyHist[module]->Fill(energy);
+
+ treeRawAmp[module] = amp;
+ treeCalibAmp[module] = energy;
+
+ if (m_debugLevel > 2 && count < 20) {
+ std::cout << "Module " << module << " : raw, calibrated energy = " << std::fixed << amp << ", " << energy << std::endl;
+ }
+ }
+
+ if (m_debugLevel > 2 && count < 20) {
+ std::cout << "Total: raw, calibrated energy = " << std::fixed << sumAmp << ", " << sumCalibAmp << std::endl;
+ }
+
+ treeRawSum = sumAmp;
+ treeRawSumOpp = (side == 0 ? zdc_ZdcAmp[1] : zdc_ZdcAmp[0]);
+ treeCalibSum = sumCalibAmp;
+
+ m_testTree->Fill();
+
+ m_testCalibSNHist->Fill(sumCalibAmp);
+
+ count++;
+ }
+ }
+
+ m_haveTest = true;
+}
+
+void ZDCNLCalibration::FillMinimizationData(TMatrixD& minimMatrix, TVectorD& minimVector,
+ std::array<int, 4> maxPowerModule, float HEDeweight,
+ const std::vector<std::vector<double> >& sums1DVec, const std::vector<double>& sumsHE,
+ const std::vector<std::vector<double> >& sums2DVec, const std::vector<double>& sumsHE2D)
+{
+ size_t numWeights = minimVector.GetNrows();
+ size_t numColumns = minimMatrix.GetNcols();
+ size_t numRows = minimMatrix.GetNrows();
+ size_t numNeutMax = sums1DVec.size();
+
+ if (numRows != numWeights || numColumns != numWeights) throw;
+
+ // if (m_debugLevel > 0) std::cout << "Filling single neutorn entries in matrix" << std::endl;
+
+ double sumFracSqr = 0;
+ for (size_t module : {0, 1, 2, 3}) {
+ sumFracSqr += m_HEFraction.at(module)*m_HEFraction.at(module);
+ }
+
+ // Start with the single neutron contributions
+ //
+ {
+ size_t index1 = 0;
+
+ for (size_t module1 : {0, 1, 2, 3}) {
+ for (size_t power1 = 0; power1 < maxPowerModule[module1]; power1++) {
+
+ size_t index2 = 0;
+ size_t sumIndexOffset = (module1*m_maxNLPower + power1)*4*m_maxNLPower;
+
+ for (size_t module2 : {0, 1, 2, 3}) {
+ for (size_t power2 = 0; power2 < maxPowerModule[module2]; power2++) {
+
+ // The N neutron contributions
+ //
+ unsigned int sum2DIndex = sumIndexOffset + power2;
+ double element = 0;
+ std::for_each(sums2DVec.begin(), sums2DVec.end(), [&] (const std::vector<double>& sums) {element += 2*sums.at(sum2DIndex);});
+
+ // The TN contribution
+ //
+ // element += 2*sumsTN2D.at(sumIndexOffset + power2);
+
+ // The HE fraction contribution
+ //
+ double sum = sumsHE2D.at(sumIndexOffset + power2);
+ sum *= HEDeweight;
+
+ if (module2 == module1) element += 2*sum;
+ element -= 2*m_HEFraction.at(module2)*sum;
+ element -= 2*m_HEFraction.at(module1)*sum;
+ element += 2*sumFracSqr*sum;
+
+ minimMatrix(index1, index2) += element;
+
+ if (m_debugLevel > 2) {
+ std::cout << "Filling module 1, power 1 = " << module1 << ", " << power1 << ", module 2, power 2 = "
+ << module2 << ", " << power2 << " with sum at index " << sumIndexOffset + power2
+ << ", value = " << minimMatrix(index1, index2) << std::endl;
+ }
+
+ index2++;
+ }
+
+ sumIndexOffset += m_maxNLPower;
+ }
+
+ size_t sum1DIndex = module1*m_maxNLPower + power1;
+
+ minimVector[index1] = 0;
+ for (unsigned int iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
+ minimVector[index1] += 2*((double) iNNeut + 1)*m_SNEnergy*sums1DVec[iNNeut].at(sum1DIndex);
+ }
+
+ if (m_debugLevel > 2) {
+ std::cout << "Filling vector module , power = " << module1 << ", " << power1 << " with sum at index " << sum1DIndex
+ << ", value = " << minimVector[index1] << std::endl;
+ }
+
+ // minimVector[index1] = 2*m_SNEnergy*sumsSN.at(sum1DIndex) + 2*2*m_SNEnergy*sumsTN.at(sum1DIndex);
+ index1++;
+ }
+ }
+ }
+
+ // if (m_debugLevel > 0) std::cout << "Filling high energy entries in matrix" << std::endl;
+
+ // {
+
+ // size_t sumIndexOffset = 0;
+ // size_t index1 = 0;
+
+ // for (size_t module1 : {0, 1, 2, 3}) {
+ // for (size_t power1 = 0; power1 < maxPowerModule[module1]; power1++) {
+ // size_t index2 = 0;
+
+ // for (size_t module2 : {0, 1, 2, 3}) {
+ // for (size_t power2 = 0; power2 < maxPowerModule[module2]; power2++) {
+
+ // double element = 0;
+ // double sum = sumsHE2D.at(sumIndexOffset + power2);
+
+ // // Now assemble to terms contributing to this matrix element
+ // //
+ // if (module2 == module1) element += 2*sum;
+ // element -= 2*m_HEFraction.at(module2)*sum;
+ // element -= 2*m_HEFraction.at(module1)*sum;
+ // element += 2*sumFracSqr*sum;
+
+ // minimMatrix(index1, index2) += m_HEDeweight*element;
+ // index2++;
+ // }
+
+ // sumIndexOffset += m_maxNLPower;
+ // }
+
+ // index1++;
+ // }
+ // }
+ // }
+
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.h b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.h
new file mode 100644
index 000000000000..6628fc7792db
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCNLCalibration.h
@@ -0,0 +1,286 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Aug 26 05:51:09 2016 by ROOT version 6.06/02
+// from TTree zdcTree/ZDC Tree
+// found on file: ../user.steinber.data15_hi.00287931.calibration_zdcCalib.merge.AOD.c932_m1533.42testNL3_output.root
+//////////////////////////////////////////////////////////
+
+#ifndef ZDCNLCalibration_h
+#define ZDCNLCalibration_h
+
+#include <TROOT.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TH1.h>
+
+#include <array>
+#include <map>
+#include <algorithm>
+#include <iostream>
+
+#include <TVectorD.h>
+#include <TMatrixD.h>
+#include <TDecompLU.h>
+
+
+struct CalibData
+{
+ std::array<std::vector<float>, 4> weights;
+ unsigned int LBStart;
+ unsigned int LBEnd;
+
+ CalibData()
+ {
+ std::vector<float> unity(1, 1);
+ weights = {unity, unity, unity, unity};
+ LBStart = 0;
+ LBEnd = 10000;
+ };
+
+CalibData(int inLBStart, int inLBEnd, const std::array<std::vector<float>, 4>& inWeights) :
+ LBStart(inLBStart), LBEnd(inLBEnd), weights(inWeights)
+ {};
+};
+
+// Header file for the classes stored in the TTree if any.
+
+class ZDCNLCalibration
+{
+ TFile* m_TFile;
+ TTree* m_tree;
+
+ size_t m_maxNLPower;
+ // float m_HEDeweight;
+ bool m_useGRL;
+ int m_debugLevel;
+
+ // size_t m_numWeights;
+ // size_t m_numFitParams;
+
+ const float m_SNEnergy;
+ const std::vector<float> m_HEFraction;
+
+ // Declaration of leaf types
+ UInt_t runNumber;
+ UInt_t eventNumber;
+ UInt_t lumiBlock;
+ UInt_t bcid;
+ UInt_t passBits;
+
+ Float_t zdc_ZdcAmp[2];
+ // Float_t zdc_ZdcEnergy[2];
+ // Float_t zdc_ZdcTime[2];
+ // Short_t zdc_ZdcStatus[2];
+ // Float_t zdc_ZdcTrigEff[2];
+ UInt_t zdc_ZdcModuleMask;
+ Float_t zdc_ZdcModuleAmp[2][4];
+ // Float_t zdc_ZdcModuleTime[2][4];
+ // Float_t zdc_ZdcModuleFitAmp[2][4];
+ // Float_t zdc_ZdcModuleFitT0[2][4];
+ // UInt_t zdc_ZdcModuleStatus[2][4];
+ // Float_t zdc_ZdcModuleChisq[2][4];
+ // Float_t zdc_ZdcModuleCalibAmp[2][4];
+ // Float_t zdc_ZdcModuleCalibTime[2][4];
+ // Float_t zdc_ZdcModuleBkgdMaxFraction[2][4];
+ // Float_t zdc_ZdcModuleAmpError[2][4];
+ Bool_t L1_ZDC_A;
+ Bool_t L1_ZDC_C;
+ Bool_t L1_ZDC_AND;
+ Bool_t L1_ZDC_A_C;
+
+ // List of branches
+ TBranch *b_runNumber; //!
+ TBranch *b_eventNumber; //!
+ TBranch *b_lumiBlock; //!
+ TBranch *b_bcid; //!
+
+ TBranch *b_passBits; //!
+
+ TBranch *b_zdc_ZdcAmp; //!
+ // TBranch *b_zdc_ZdcEnergy; //!
+ // TBranch *b_zdc_ZdcTime; //!
+ // TBranch *b_zdc_ZdcStatus; //!
+ // TBranch *b_zdc_ZdcTrigEff; //!
+ TBranch *b_zdc_ZdcModuleMask; //!
+ TBranch *b_zdc_ZdcModuleAmp; //!
+ // TBranch *b_zdc_ZdcModuleTime; //!
+ // TBranch *b_zdc_ZdcModuleFitAmp; //!
+ // TBranch *b_zdc_ZdcModuleFitT0; //!
+ // TBranch *b_zdc_ZdcModuleStatus; //!
+ // TBranch *b_zdc_ZdcModuleChisq; //!
+ // TBranch *b_zdc_ZdcModuleCalibAmp; //!
+ // TBranch *b_zdc_ZdcModuleCalibTime; //!
+ // TBranch *b_zdc_ZdcModuleBkgdMaxFraction; //!
+ // TBranch *b_zdc_ZdcModuleAmpError; //!
+ TBranch *b_L1_ZDC_A; //!
+ TBranch *b_L1_ZDC_C; //!
+ TBranch *b_L1_ZDC_AND; //!
+ TBranch *b_L1_ZDC_A_C; //!
+
+ // Map that keeps lumi block and event associations to optimize
+ // processing of lumi block ranges
+ //
+ typedef multimap<unsigned int, std::pair<unsigned int, unsigned int> > LBEvtMap;
+ LBEvtMap m_LumiBlockEvtMap;
+
+public:
+
+ std::array<std::map<std::string, CalibData>, 2> m_calibrations;
+
+ bool m_haveTest;
+ TH1D* m_testCalibSNHist;
+ std::array<TH1D*, 4> m_testCalibHEFracHist;
+ std::array<TH1D*, 4> m_testCalibEnergyHist;
+ TTree* m_testTree;
+
+
+ ZDCNLCalibration(std::string file, int maxNLPower = 3, bool useGRL = true, int debugLevel = 0);
+ virtual ~ZDCNLCalibration() {}
+
+ std::array<float, 4> FindSNPeaks(size_t LBLow, size_t LBHigh, size_t side);
+
+ std::pair<float, float> FindSNRange(size_t LBLow, size_t LBHigh, size_t side);
+
+ std::pair<std::pair<float, float>,std::pair<float, float> > FindSNTwoNRanges(size_t LBLow, size_t LBHigh, size_t side);
+
+ void SetDefaultCalibration(size_t side, const CalibData& calib) {
+ AddCalibration(side, "default", calib);
+ }
+
+ void Calibrate(size_t side, std::string calibInput, std::string calibOutput,
+ size_t LBLow, size_t LBHigh, std::array<int, 4> maxPowerModule,
+ std::vector<std::pair<double, double> > nNeutERange,
+ bool excludeHE, float heSumThresh, float HEDeweight);
+
+ void TestCalibration(int side, std::string calibName);
+
+ TH1* GetTestSNHist() {return (m_haveTest ? m_testCalibSNHist : 0);}
+ TH1* GetTestFracHist(size_t module) {return (m_haveTest ? m_testCalibHEFracHist.at(module) : 0);}
+
+ TTree* GetTestTree() {return m_testTree;}
+
+private:
+
+ void FillLumiBlockEvtMap();
+
+ void FillMinimizationData(TMatrixD& minimMatrix, TVectorD& minimVector,
+ std::array<int, 4> maxPowerModule, float HEDeweight,
+ const std::vector<std::vector<double> >& sums1DVec, const std::vector<double>& sumsHE,
+ const std::vector<std::vector<double> >& sums2DVec, const std::vector<double>& sumsHE2D);
+
+ void AddCalibration(size_t side, std::string tag, const CalibData& calib);
+ CalibData GetCalibration(size_t side, std::string tag);
+
+ // Add the four amplitudes from a given event to the set of sums used in the NL weight fitting
+ // we need to make all combinations of amplitude products and for all combinations of powers
+ //
+ void AddToSums(std::vector<double>& sums1D, std::vector<double>& sums2D, float* amps)
+ {
+ size_t index1D = 0;
+ size_t index2D = 0;
+
+ for (size_t module1 : {0, 1, 2, 3}) {
+ double amp1 = amps[module1];
+ double amp1Pow = amp1;
+
+ for (size_t power1 = 0; power1 < m_maxNLPower; power1++) {
+ for (size_t module2 : {0, 1, 2, 3}) {
+ double amp2 = amps[module2];
+ double amp2Pow = amp2;
+
+ for (size_t power2 = 0; power2 < m_maxNLPower; power2++) {
+ sums2D.at(index2D++) += amp1Pow*amp2Pow;
+ amp2Pow *= amp2;
+ }
+ }
+
+ sums1D.at(index1D++) += amp1Pow;
+ amp1Pow *= amp1;
+ }
+ }
+ }
+
+ static double CalculateEnergy(const float* moduleAmps, const CalibData& calib)
+ {
+ double energy = 0;
+
+ for (size_t module : {0, 1, 2, 3}) {
+ float amp = moduleAmps[module];
+ float ampPow = amp;
+
+ for (size_t power = 0; power < calib.weights[module].size(); power++) {
+ energy += calib.weights[module][power]*ampPow;
+ ampPow *= amp;;
+ }
+ }
+
+ return energy;
+ }
+
+};
+
+#endif
+
+#ifdef ZDCNLCalibration_cxx
+ZDCNLCalibration::ZDCNLCalibration(std::string file, int maxNLPower, bool useGRL, int debugLevel) :
+ m_TFile(0), m_tree(0), m_maxNLPower(maxNLPower), //m_HEDeweight(heDeweight),
+ m_useGRL(useGRL), m_debugLevel(debugLevel),
+ // m_numWeights(m_maxNLPower*4), m_numFitParams(m_maxNLPower *4),
+ m_SNEnergy(2510), m_HEFraction({0.31, 0.27, 0.21, 0.21}),
+ m_haveTest(false)
+{
+ std::cout << "Initializing ZDCNLCalibration with debug level " << m_debugLevel << std::endl;
+
+ TFile* temp_p = new TFile(file.c_str());
+ if (!temp_p->IsOpen()) {
+ std::cout << "Error opening input file " << file << std::endl;
+ return;
+ }
+
+ m_TFile = temp_p;
+ m_tree = static_cast<TTree*>(m_TFile->GetObjectChecked("zdcTree","TTree"));
+ if (!m_tree) {
+ std::cout << "Error reading tree from input file " << file << std::endl;
+ return;
+ }
+
+ m_tree->SetBranchAddress("runNumber", &runNumber, &b_runNumber);
+ m_tree->SetBranchAddress("eventNumber", &eventNumber, &b_eventNumber);
+ m_tree->SetBranchAddress("lumiBlock", &lumiBlock, &b_lumiBlock);
+ m_tree->SetBranchAddress("bcid", &bcid, &b_bcid);
+ m_tree->SetBranchAddress("passBits", &passBits, &b_passBits);
+
+ if (m_tree->FindBranch("zdc_ZdcModuleMask")) {
+ m_tree->SetBranchAddress("zdc_ZdcAmp", zdc_ZdcAmp, &b_zdc_ZdcAmp);
+ m_tree->SetBranchAddress("zdc_ZdcModuleMask", &zdc_ZdcModuleMask, &b_zdc_ZdcModuleMask);
+ m_tree->SetBranchAddress("zdc_ZdcModuleAmp", zdc_ZdcModuleAmp, &b_zdc_ZdcModuleAmp);
+ }
+ else if (m_tree->FindBranch("zdc_ModuleMask")) {
+ m_tree->SetBranchAddress("zdc_OptSumAmp", zdc_ZdcAmp, &b_zdc_ZdcAmp);
+ m_tree->SetBranchAddress("zdc_ModuleMask", (int*) &zdc_ZdcModuleMask, &b_zdc_ZdcModuleMask);
+ m_tree->SetBranchAddress("zdc_OptAmp", zdc_ZdcModuleAmp, &b_zdc_ZdcModuleAmp);
+ }
+ else throw;
+
+ m_tree->SetBranchAddress("L1_ZDC_A", &L1_ZDC_A, &b_L1_ZDC_A);
+ m_tree->SetBranchAddress("L1_ZDC_C", &L1_ZDC_C, &b_L1_ZDC_C);
+ m_tree->SetBranchAddress("L1_ZDC_AND", &L1_ZDC_AND, &b_L1_ZDC_AND);
+ m_tree->SetBranchAddress("L1_ZDC_A_C", &L1_ZDC_A_C, &b_L1_ZDC_A_C);
+
+ if (m_debugLevel > 0) {
+ std::cout << "Filling LumiBlock-event map" << std::endl;
+ }
+
+ FillLumiBlockEvtMap();
+
+ // Provide "null" default calibration
+ //
+ AddCalibration(0, "default", CalibData());
+ AddCalibration(1, "default", CalibData());
+}
+
+#endif // #ifdef ZDCNLCalibration_cxx
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.C b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.C
new file mode 100644
index 000000000000..223f0137777e
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.C
@@ -0,0 +1,412 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#define ZDCTreeAnalysis_cxx
+#include "ZDCTreeAnalysis.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TLatex.h>
+#include <TPaveStats.h>
+
+#include <cmath>
+
+template <typename T> T Sqr(const T& inT) {return inT*inT;}
+
+void ZDCTreeAnalysis::InitInternal()
+{
+ for (int iside = 0; iside < 2; iside++) {
+ for (int imod = 0; imod < 4; imod++) {
+ _fixTau1[iside][imod] = false;
+ _fixTau2[iside][imod] = false;
+ _moduleTau1[iside][imod] = 4;
+ _moduleTau2[iside][imod] = 20;
+
+ _chisqDivAmpCutHG[iside][imod] = 1e6;
+ _chisqDivAmpCutLG[iside][imod] = 1e6;
+ _DeltaT0CutLow[iside][imod] = -100;
+ _DeltaT0CutHigh[iside][imod] = 100;
+ _HGOverFlowADC[iside][imod] = 800;
+ _HGUnderFlowADC[iside][imod] = 25;
+
+ for (size_t ipar = 0; ipar < 3; ipar++) {
+ _T0SlewCoeffHG[iside][imod][ipar] = 0;
+ _T0SlewCoeffLG[iside][imod][ipar] = 0;
+ }
+
+ _modECalib[iside][imod] = 1;
+ _haveLBDepECalib[iside][imod] = false;
+ _haveLBDepT0[iside][imod] = false;
+ }
+ }
+
+ zdc_samplesSub = new std::vector<float>(7, 0);
+ zdc_samplesDeriv = new std::vector<float>(7, 0);
+ zdc_samplesDeriv2nd = new std::vector<float>(7, 0);
+}
+
+
+void ZDCTreeAnalysis::SetupBunchTrains()
+{
+ GetEntry(0);
+
+ _BCIDGap.assign(3300, 0);
+ _BCIDPosInTrain.assign(3300, 0);
+
+ std::cout << "Run = " << runNumber << std::endl;
+
+ if (runNumber >= 287706) {
+
+ int trainStartBCIDS[21] = {5, 155, 305, 446, 596, 746, 896, 1046, 1196, 1487, 1637, 1787, 1937, 2087, 2228, 2378, 2528, 2678, 2828, 2978, 3119};
+ int trainNumBCIDS[21] = {24, 24, 20, 24, 24, 24, 24, 24, 20, 24, 24, 24, 24, 22, 24, 24, 24, 24, 24, 22, 24};
+
+ int lastBCID = -1;
+ int numBCID = 0;
+
+ for (int train = 0; train < 21; train++) {
+ int bcid = trainStartBCIDS[train];
+
+ std::cout << "Train " << train << ", bcids = ";
+
+ for (int PIT = 0; PIT < trainNumBCIDS[train]; PIT++) {
+ int gapNext = PIT % 2 == 0 ? 4 : 6;
+
+ _BCIDGap[bcid] = bcid - lastBCID;
+ _BCIDPosInTrain[bcid] = PIT;
+
+ std::cout << bcid << ", ";
+
+ numBCID++;
+ lastBCID = bcid;
+ bcid += gapNext;
+ }
+
+ std::cout << std::endl;
+ }
+
+ std::cout << "Found " << numBCID << "bcids" << std::endl;
+ }
+}
+
+void ZDCTreeAnalysis::OpenOutputTree(std::string file)
+{
+ TFile* outputFile = new TFile(file.c_str(), "recreate");
+ if (!outputFile->IsOpen()) {
+ std::cout << "Error opening output file " << file << ", quitting" << std::endl;
+ return;
+ }
+
+ _outTFile = outputFile;
+
+ _outTree = fChain->CloneTree(0);
+ _outTree->Branch("bcidGap", &bcidGap, "bcidGap/I");
+ _outTree->Branch("bcidPosInTrain", &bcidPosInTrain, "bcidPosInTrain/I");
+
+ _outTree->Branch("zdc_SumAmp", zdc_SumAmp, "zdc_SumAmp[2]/F");
+ _outTree->Branch("zdc_SumCalibAmp", zdc_SumCalibAmp, "zdc_SumCalibAmp[2]/F");
+ _outTree->Branch("zdc_AvgTime", zdc_AvgTime, "zdc_AvgTime[2]/F");
+ _outTree->Branch("zdc_ModuleMask", &zdc_ModuleMask, "zdc_ModuleMask/I");
+
+ _outTree->Branch("zdc_Status", zdc_Status, "zdc_Status[2][4]/I");
+ _outTree->Branch("zdc_Amp", zdc_Amp, "zdc_Amp[2][4]/F");
+ _outTree->Branch("zdc_FitT0", zdc_FitT0, "zdc_FitT0[2][4]/F");
+ _outTree->Branch("zdc_T0Corr", zdc_T0Corr, "zdc_T0Corr[2][4]/F");
+ _outTree->Branch("zdc_FitChisq", zdc_FitChisq, "zdc_FitChisq[2][4]/F");
+
+ _outTree->Branch("zdc_CalibAmp", zdc_CalibAmp, "zdc_CalibAmp[2][4]/F");
+ _outTree->Branch("zdc_CalibTime", zdc_CalibTime, "zdc_CalibTime[2][4]/F");
+
+ _outTree->Branch("zdc_MaxADC", zdc_MaxADC, "zdc_MaxADC[2][4]/f");
+ _outTree->Branch("zdc_MinADC", zdc_MinADC, "zdc_MinADC[2][4]/f");
+ _outTree->Branch("zdc_MaxADCSample", zdc_MaxADCSample, "zdc_MaxADCSample[2][4]/I");
+ _outTree->Branch("zdc_MinADCSample", zdc_MaxADCSample, "zdc_MinADCSample[2][4]/I");
+ _outTree->Branch("zdc_Min2ndDeriv", zdc_Min2ndDeriv, "zdc_Min2ndDeriv[2][4]/F");
+ _outTree->Branch("zdc_Min2ndDerivSample", zdc_Min2ndDerivSample, "zdc_Min2ndDerivSample[2][4]/I");
+
+ _outTree->Branch("zdc_samplesSub", &zdc_samplesSub);
+ _outTree->Branch("zdc_samplesDeriv", &zdc_samplesDeriv);
+ _outTree->Branch("zdc_samplesDeriv2nd", &zdc_samplesDeriv2nd);
+
+ _doOutput = true;
+}
+
+void ZDCTreeAnalysis::CloseOutputTree()
+{
+ if (!_doOutput) return;
+
+ _outTree->AutoSave();
+ delete _outTFile;
+ _outTFile = 0;
+ _outTree = 0;
+ _doOutput = false;
+}
+
+void ZDCTreeAnalysis::Loop(int numEntries, int startEntry)
+{
+ if (fChain == 0) return;
+
+ int loopStart = (startEntry == -1 ? _currentEntry : startEntry);
+
+ Long64_t nentries = fChain->GetEntriesFast() - loopStart;
+
+ Long64_t numLoopEntry = (numEntries != -1 ? numEntries : nentries);
+
+ Long64_t nbytes = 0, nb = 0;
+
+ _inLoop = true;
+
+ for (Long64_t jentry = loopStart; jentry < loopStart + numLoopEntry; jentry++) {
+ Long64_t ientry = LoadTree(jentry);
+ if (ientry < 0) break;
+ nb = fChain->GetEntry(jentry);
+ _currentEntry = jentry;
+ nbytes += nb;
+
+ DoAnalysis();
+
+ if (jentry > 1) {
+ float log10Entry = std::floor(std::log ((float) jentry)/std::log(10));
+ int modEntry = jentry % (int) pow(10, log10Entry);
+
+ if (modEntry == 0) {
+ std::cout << "Processed event " << jentry << std::endl;
+ }
+ }
+ }
+
+ _inLoop = false;
+}
+
+void ZDCTreeAnalysis::DoAnalysis()
+{
+ // Starting a new event
+ //
+ _dataAnalyzer_p->StartEvent(lumiBlock);
+
+ bcidGap = _BCIDGap[bcid];
+ bcidPosInTrain = _BCIDPosInTrain[bcid];
+
+ if (_doOutput) {
+ // //
+ // // Initialize outout arrays
+ // //
+ }
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+
+ static std::vector<float> HGADCSamples(7);
+ static std::vector<float> LGADCSamples(7);
+
+ // Copy data into a proper vector of floats
+ //
+ if (module == 0) {
+ std::copy(&zdc_raw[side][module][0][1][0], &zdc_raw[side][module][0][1][7], LGADCSamples.begin());
+ std::copy(&zdc_raw[side][module][1][1][0], &zdc_raw[side][module][1][1][7], HGADCSamples.begin());
+ }
+ else {
+ std::copy(&zdc_raw[side][module][0][0][0], &zdc_raw[side][module][0][0][7], LGADCSamples.begin());
+ std::copy(&zdc_raw[side][module][1][0][0], &zdc_raw[side][module][1][0][7], HGADCSamples.begin());
+ }
+
+ // Load the data
+ //
+ _dataAnalyzer_p->LoadAndAnalyzeData(side, module, HGADCSamples, LGADCSamples);
+ }
+ }
+
+
+ // We're done
+ //
+ _dataAnalyzer_p->FinishEvent();
+
+ zdc_ModuleMask = _dataAnalyzer_p->GetModuleMask();
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+
+ const ZDCPulseAnalyzer* pulseAna_p = _dataAnalyzer_p->GetPulseAnalyzer(side, module);
+
+ *zdc_samplesSub = pulseAna_p->GetSamplesSub();
+ *zdc_samplesDeriv = pulseAna_p->GetSamplesDeriv();
+ *zdc_samplesDeriv2nd = pulseAna_p->GetSamplesDeriv2nd();
+
+ zdc_Min2ndDeriv[side][module] = pulseAna_p->GetMinDeriv2nd();
+ zdc_Min2ndDerivSample[side][module] = pulseAna_p->GetMinDeriv2ndIndex();
+
+ zdc_MaxADC[side][module] = pulseAna_p->GetMaxADC();
+ zdc_MinADC[side][module] = pulseAna_p->GetMinADC();
+
+ zdc_MaxADCSample[side][module] = pulseAna_p->GetMaxADCSample();
+ zdc_MinADCSample[side][module] = pulseAna_p->GetMinADCSample();
+
+ zdc_Status[side][module] = pulseAna_p->GetStatusMask();
+ zdc_Amp[side][module] = pulseAna_p->GetAmplitude();
+ zdc_FitT0[side][module] = pulseAna_p->GetFitT0();
+ zdc_T0Corr[side][module] = pulseAna_p->GetT0Corr();
+ zdc_FitChisq[side][module] = pulseAna_p->GetChisq();
+
+ zdc_CalibAmp[side][module] = _dataAnalyzer_p->GetModuleCalibAmplitude(side, module);
+ zdc_CalibTime[side][module] = _dataAnalyzer_p->GetModuleCalibTime(side, module);
+ }
+
+ zdc_SumAmp[side] = _dataAnalyzer_p->GetModuleSum(side);
+ zdc_SumCalibAmp[side] = _dataAnalyzer_p->GetCalibModuleSum(side);
+ zdc_AvgTime[side] = _dataAnalyzer_p->GetAverageTime(side);
+ zdc_sideFailed[side] = _dataAnalyzer_p->SideFailed(side);
+
+
+ }
+
+
+ if (_doOutput) {
+ _outTree->Fill();
+ }
+}
+
+
+void ZDCTreeAnalysis::PlotFits(std::string canvasSavePath)
+{
+ static bool first = true;
+ static std::array<TCanvas*, 2> plotCanvases;
+ static TLatex* statusLabel;
+ static TLatex* sideModuleLabel;
+ static TLatex* highLowGainLabel;
+ static TLatex* eventLabel;
+ static TLatex* warningLabel;
+
+ if (first) {
+ plotCanvases[0] = new TCanvas("ZDCFitsSide0","ZDCFitsSide0", 1000, 800);
+ plotCanvases[1] = new TCanvas("ZDCFitsSide1","ZDCFitsSide1", 1000, 800);
+
+ plotCanvases[0]->Divide(2, 2, 0.001, 0.0001);
+ plotCanvases[1]->Divide(2, 2, 0.001, 0.0001);
+
+ first = false;
+
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(111);
+
+ sideModuleLabel = new TLatex;
+ sideModuleLabel->SetTextFont(42);
+ sideModuleLabel->SetTextSize(0.05);
+ sideModuleLabel->SetTextAlign(12);
+ sideModuleLabel->SetTextColor(2);
+ sideModuleLabel->SetNDC(1);
+
+ statusLabel = new TLatex;
+ statusLabel->SetTextFont(42);
+ statusLabel->SetTextSize(0.045);
+ statusLabel->SetTextAlign(12);
+ statusLabel->SetNDC(1);
+
+ highLowGainLabel = new TLatex;
+ highLowGainLabel->SetTextFont(42);
+ highLowGainLabel->SetTextSize(0.045);
+ highLowGainLabel->SetTextAlign(12);
+ highLowGainLabel->SetNDC(1);
+
+ warningLabel = new TLatex;
+ warningLabel->SetTextFont(62);
+ warningLabel->SetTextColor(2);
+ warningLabel->SetTextSize(0.045);
+ warningLabel->SetTextAlign(12);
+ warningLabel->SetNDC(1);
+ }
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ TVirtualPad* pad_p = plotCanvases[side]->cd(module + 1);
+
+ pad_p->SetTopMargin(0.03);
+ pad_p->SetRightMargin(0.03);
+
+ const ZDCPulseAnalyzer* pulseAna_p = _dataAnalyzer_p->GetPulseAnalyzer(side, module);
+
+ const TH1* hist_p = pulseAna_p->GetHistogramPtr();
+ TH1* cloneHist_p = (TH1*) hist_p->Clone((std::string(hist_p->GetName()) + "_clone").c_str());
+
+ // Set the maximum to make sure we include the peak of the function
+ //
+ TF1* func_p = (TF1*) cloneHist_p->GetListOfFunctions()->First();
+ double funcMax = func_p->GetMaximum();
+
+ if (pulseAna_p->HavePulse()) {
+ if (cloneHist_p->GetMaximum() < funcMax*1.1) cloneHist_p->SetMaximum(funcMax*1.1);
+ if (cloneHist_p->GetMaximum() < 20) cloneHist_p->SetMaximum(20);
+ }
+ else {
+ float maxADC = pulseAna_p->GetMaxADC();
+ float minADC = pulseAna_p->GetMinADC();
+
+ cloneHist_p->SetMaximum(maxADC*1.1);
+ cloneHist_p->SetMinimum(std::min(0.0,minADC*1.1));
+ }
+
+ if (pulseAna_p->HavePulse()) cloneHist_p->Draw();
+ else cloneHist_p->Draw("hist");
+
+ gPad->Modified();
+ gPad->Update();
+
+ TPaveStats* stats_p = (TPaveStats *)cloneHist_p->GetListOfFunctions()->FindObject("stats");
+ if (stats_p) {
+ stats_p->SetX1NDC(0.135);
+ stats_p->SetX2NDC(0.4);
+ stats_p->SetFitFormat("5.3f");
+ }
+
+ unsigned int status = pulseAna_p->GetStatusMask();
+
+ std::ostringstream sideModule, statusHex;
+ sideModule << "Side " << side << ", module " << module;
+ statusHex << "0x" << std::hex << status;
+
+ sideModuleLabel->DrawLatex(0.135, 0.6, sideModule.str().c_str());
+
+ std::string statusString = "status = " + statusHex.str();
+ statusLabel->DrawLatex(0.135, 0.54, statusString.c_str());
+
+ if (pulseAna_p->UseLowGain()) {
+ highLowGainLabel->DrawLatex(0.135, 0.48, "low gain");
+ }
+ else {
+ highLowGainLabel->DrawLatex(0.135, 0.48, "high gain");
+ }
+
+ std::string warning = "";
+
+ if (pulseAna_p->Failed()) warning += "Failed";
+
+ if (pulseAna_p->BadT0()) {
+ if (warning != "") warning += ",";
+ warning += "BadT0";
+ }
+
+ if (pulseAna_p->BadChisq()) {
+ if (warning != "") warning += ",";
+ warning += "BadChisq";
+ }
+
+ if (warning != "") warningLabel->DrawLatex(0.135, 0.42, warning.c_str());
+
+ gPad->Modified();
+ gPad->Update();
+ }
+
+ plotCanvases[side]->Modified();
+ plotCanvases[side]->Update();
+ }
+}
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.h b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.h
new file mode 100644
index 000000000000..19be5440a2cb
--- /dev/null
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/tools/ZDCTreeAnalysis.h
@@ -0,0 +1,512 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Wed Dec 16 10:23:25 2015 by ROOT version 6.04/12
+// from TTree zdcTree/ZDC Tree
+// found on file: data15_hi.00287259.calibration_zdcCalib.recon.AOD.c931.root
+//////////////////////////////////////////////////////////
+
+#ifndef ZDCTreeAnalysis_h
+#define ZDCTreeAnalysis_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TSpline.h>
+
+// Header file for the classes stored in the TTree if any.
+#include <vector>
+#include <iostream>
+#include <sstream>
+
+#include "ZDCDataAnalyzer.h"
+
+class ZDCTreeAnalysis {
+public :
+ TTree *fChain; //!pointer to the analyzed TTree or TChain
+ Int_t fCurrent; //!current Tree number in a TChain
+
+ TFile* _outTFile;
+ TTree* _outTree;
+ bool _doOutput;
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+ // Declaration of leaf types
+ UInt_t runNumber;
+ UInt_t eventNumber;
+ UInt_t lumiBlock;
+ UInt_t bcid;
+ ULong64_t trigger;
+ UInt_t trigger_TBP;
+ UInt_t tbp[16];
+ UInt_t tav[16];
+
+ UInt_t passBits;
+ UShort_t zdc_raw[2][4][2][2][7];
+
+ Float_t fcalEt;
+ Int_t ntrk;
+ Int_t nvx;
+ Float_t vx[3];
+ Int_t vxntrk;
+ Float_t vxcov[6];
+ UShort_t mbts_countA;
+ UShort_t mbts_countC;
+ Float_t mbts_timeA;
+ Float_t mbts_timeC;
+ Float_t mbts_timeDiff;
+
+ Bool_t L1_ZDC_A;
+ Float_t ps_L1_ZDC_A;
+ Bool_t L1_ZDC_C;
+ Float_t ps_L1_ZDC_C;
+ Bool_t L1_ZDC_AND;
+ Float_t ps_L1_ZDC_AND;
+ Bool_t L1_ZDC_A_C;
+ Float_t ps_L1_ZDC_A_C;
+
+ // Data for branches that are to be added to the tree
+ //
+ int bcidGap;
+ int bcidPosInTrain;
+
+ // Data obtained from the analysis
+ //
+ // Summary data
+ //
+ int zdc_ModuleMask;
+ float zdc_SumAmp[2];
+ float zdc_SumCalibAmp[2];
+ float zdc_AvgTime[2];
+ bool zdc_sideFailed[2];
+
+ // Per-module data
+ //
+ std::vector<float>* zdc_samplesSub;
+ std::vector<float>* zdc_samplesDeriv;
+ std::vector<float>* zdc_samplesDeriv2nd;
+
+ float zdc_Presample[2][4];
+
+ float zdc_MaxADC[2][4];
+ int zdc_MaxADCSample[2][4];
+
+ float zdc_MinADC[2][4];
+ int zdc_MinADCSample[2][4];
+
+ float zdc_Min2ndDeriv[2][4];
+ int zdc_Min2ndDerivSample[2][4];
+
+ int zdc_Status[2][4];
+ float zdc_Amp[2][4];
+ float zdc_FitT0[2][4];
+ float zdc_FitChisq[2][4];
+ float zdc_T0Corr[2][4];
+
+ float zdc_CalibAmp[2][4];
+ float zdc_CalibTime[2][4];
+
+ // The object responsible for the actual analysis
+ //
+ ZDCDataAnalyzer* _dataAnalyzer_p;
+
+private:
+ // Critical construction parameters
+ //
+ int _nSample;
+ float _deltaTSample;
+ int _preSampleIdx;
+ std::string _fitFunction;
+
+ ZDCDataAnalyzer::ZDCModuleFloatArray _peak2ndDerivMinSamples;
+ ZDCDataAnalyzer::ZDCModuleFloatArray _peak2ndDerivMinThresholds;
+
+ // Cut quantities
+ //
+ float _HGOverFlowADC[2][4];
+ float _HGUnderFlowADC[2][4];
+ float _DeltaT0CutLow[2][4];
+ float _DeltaT0CutHigh[2][4];
+ float _chisqDivAmpCutHG[2][4];
+ float _chisqDivAmpCutLG[2][4];
+
+ // Per-module calibration factors
+ //
+ // Time-dependent
+ //
+ bool _haveLBDepT0[2][4];
+ TSpline* _moduleT0HGLB[2][4];
+ TSpline* _moduleT0LGLB[2][4];
+
+ bool _haveCalibrations;
+ float _modECalib[2][4];
+
+ bool _haveLBDepECalib[2][4];
+ TSpline* _modECalibLB[2][4];
+
+ // Allow per-module Tau1, tau2 settings
+ //
+ bool _haveModuleSettings[2][4];
+ bool _fixTau1[2][4];
+ bool _fixTau2[2][4];
+ float _moduleTau1[2][4];
+ float _moduleTau2[2][4];
+ float _moduleT0LG[2][4];
+ float _moduleT0HG[2][4];
+
+ std::array<std::array<float, 4>, 2> _moduleHGGains;
+
+ // Time slewing corrections
+ //
+ float _T0SlewCoeffHG[2][4][3];
+ float _T0SlewCoeffLG[2][4][3];
+
+ // Bunch information
+ //
+ std::vector<std::set<int> > _trains;
+ std::vector<int> _BCIDGap;
+ std::vector<int> _BCIDPosInTrain;
+
+ // Which entry we're reading
+ //
+ bool _inLoop;
+ int _currentEntry;
+ unsigned int _runNumber;
+
+private:
+ // List of branches
+ TBranch *b_runNumber; //!
+ TBranch *b_eventNumber; //!
+ TBranch *b_lumiBlock; //!
+ TBranch *b_bcid; //!
+ TBranch *b_trigger; //!
+ TBranch *b_trigger_TBP; //!
+ TBranch *b_tbp; //!
+ TBranch *b_tav; //!
+ TBranch *b_decisions; //!
+ TBranch *b_prescales; //!
+ TBranch *b_passBits; //!
+
+ TBranch *b_zdc_raw; //!
+ // TBranch *b_fcalEt; //!
+ // TBranch *b_ntrk; //!
+ // TBranch *b_nvx; //!
+ // TBranch *b_vx; //!
+ // TBranch *b_vxntrk; //!
+ // TBranch *b_vxcov; //!
+ // TBranch *b_mbts_countA; //!
+ // TBranch *b_mbts_countC; //!
+ // TBranch *b_mbts_timeA; //!
+ // TBranch *b_mbts_timeC; //!
+ // TBranch *b_mbts_timeDiff; //!
+
+ TBranch *b_L1_ZDC_A; //!
+ TBranch *b_ps_L1_ZDC_A; //!
+ TBranch *b_L1_ZDC_C; //!
+ TBranch *b_ps_L1_ZDC_C; //!
+ TBranch *b_L1_ZDC_AND; //!
+ TBranch *b_ps_L1_ZDC_AND; //!
+ TBranch *b_L1_ZDC_A_C; //!
+ TBranch *b_ps_L1_ZDC_A_C; //!
+
+ void InitInternal();
+ void SetupBunchTrains();
+
+public:
+
+ ZDCTreeAnalysis(std::string filename, int nSample = 7, double deltaT = 12.5, int preSamplIdx = 1, std::string fitFunction = "FermiExp");
+
+ void SetADCOverUnderflowValues(const ZDCDataAnalyzer::ZDCModuleFloatArray& HGOverFlowADC,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& HGUnderFlowADC,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& LGOverFlowADC)
+ {
+ _dataAnalyzer_p->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+ }
+
+ void SetCutValues(const ZDCDataAnalyzer::ZDCModuleFloatArray& chisqDivAmpCutHG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& chisqDivAmpCutLG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& DeltaT0CutLowHG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& DeltaT0CutHighHG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& DeltaT0CutLowLG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& DeltaT0CutHighLG)
+
+ {
+ _dataAnalyzer_p->SetCutValues(chisqDivAmpCutHG, chisqDivAmpCutLG,
+ DeltaT0CutLowHG, DeltaT0CutHighHG,
+ DeltaT0CutLowLG, DeltaT0CutHighLG);
+ }
+
+ void SetTauT0Values(const ZDCDataAnalyzer::ZDCModuleBoolArray& fixTau1,
+ const ZDCDataAnalyzer::ZDCModuleBoolArray& fixTau2,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& tau1,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& tau2,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& t0HG,
+ const ZDCDataAnalyzer::ZDCModuleFloatArray& t0LG)
+ {
+ _dataAnalyzer_p->SetTauT0Values(fixTau1, fixTau2, tau1, tau2, t0HG, t0LG);
+ }
+
+
+ void SetDebugLevel(int debugLevel = 0)
+ {
+ _dataAnalyzer_p->SetDebugLevel(debugLevel);
+ }
+
+ void LoadEnergyCalibrations(const std::array<std::array<TSpline*, 4>, 2>& calibSplines)
+ {
+ _dataAnalyzer_p->LoadEnergyCalibrations(calibSplines);
+ }
+
+ void LoadT0Calibrations(const std::array<std::array<TSpline*, 4>, 2>& calibSplinesHG,
+ const std::array<std::array<TSpline*, 4>, 2>& calibSplinesLG)
+ {
+ _dataAnalyzer_p->LoadtT0Calibrations(calibSplinesHG, calibSplinesLG);
+ }
+
+ void SetLBDepT0(int iside, int imod, TSpline* t0SplineLG, TSpline* t0SplineHG);
+
+ void SetSlewingCoeff(const std::array<std::array<std::vector<float>, 4>, 2>& HGParamArr,
+ const std::array<std::array<std::vector<float>, 4>, 2>& LGParamArr)
+ {
+ _dataAnalyzer_p->SetTimingCorrParams(HGParamArr, LGParamArr);
+ }
+
+ void OpenOutputTree(std::string file);
+ void CloseOutputTree();
+
+ void PlotFits(std::string canvasSavePath = "");
+
+ virtual ~ZDCTreeAnalysis();
+ virtual Int_t Cut(Long64_t entry);
+ virtual Int_t GetEntry(Long64_t entry);
+ virtual Long64_t LoadTree(Long64_t entry);
+ virtual void Init(TTree *tree);
+ virtual void Loop(int numEntries = -1, int startEntry = 0);
+ virtual Bool_t Notify();
+ virtual void Show(Long64_t entry = -1);
+
+ void LoadEntry(int entry)
+ {
+ GetEntry(entry);
+ DoAnalysis();
+ }
+
+ void LoadNextEntry()
+ {
+ GetEntry(_currentEntry + 1);
+ DoAnalysis();
+ }
+
+ unsigned int GetRunNumber() const {return _runNumber;}
+
+ void DoAnalysis();
+};
+
+
+#endif
+
+#ifdef ZDCTreeAnalysis_cxx
+
+ZDCTreeAnalysis::ZDCTreeAnalysis(std::string filename, int nSample, double deltaT, int preSamplIdx, std::string fitFunction) :
+ fChain(0), _outTFile(0), _outTree(0),
+ _nSample(nSample), _deltaTSample(deltaT), _preSampleIdx(preSamplIdx),
+ _doOutput(false), _currentEntry(-1), _inLoop(false),
+ _haveCalibrations(false)
+{
+
+ // Open the fiel, extract and initialize the tree
+ //
+ TTree* tree;
+ TFile* f = new TFile(filename.c_str());
+ f->GetObject("zdcTree",tree);
+ Init(tree);
+
+ // Capture the run number
+ //
+ GetEntry(1);
+ _runNumber = runNumber;
+
+ InitInternal();
+
+ // Figure out the bunch spacing
+ //
+ SetupBunchTrains();
+
+ // Set the HV gains
+ //
+ _moduleHGGains = {9.51122, 9.51980, 9.51122, 9.51122,
+ 9.51415, 9.5049, 9.51659, 9.51415};
+
+ _peak2ndDerivMinSamples = {3, 3, 3, 2,
+ 2, 2, 2, 2};
+
+ _peak2ndDerivMinThresholds = {-7, -7, -7, -7,
+ -7, -7, -7, -7};
+
+ _dataAnalyzer_p = new ZDCDataAnalyzer(_nSample, _deltaTSample, _preSampleIdx, fitFunction,
+ _peak2ndDerivMinSamples, _peak2ndDerivMinThresholds);
+
+/* 81Undelayed */
+/* 81H1: 9.5049 */
+/* 81H2: 9.51659 */
+/* 81H3: 9.50119 */
+/* 81EM: 9.51415 */
+/* 81Delayed */
+/* 81H1: 9.52632 */
+/* 81H2: 9.49662 */
+/* 81H3: 9.50853 */
+/* 81EM: 9.50842 */
+/* 12Undelayed */
+/* 12H1: 9.5198 */
+/* 12H2: 9.51122 */
+/* 12H3: 9.51122 */
+/* 12EM: 9.51122 */
+/* 12Delayed */
+/* 12H1: 9.51237 */
+/* 12H2: 9.50178 */
+/* 12H3: 9.50178 */
+/* 12EM: 9.50178 */
+}
+
+ZDCTreeAnalysis::~ZDCTreeAnalysis()
+{
+ if (!fChain) return;
+ delete fChain->GetCurrentFile();
+}
+
+Int_t ZDCTreeAnalysis::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+ if (!fChain) return 0;
+ int result = fChain->GetEntry(entry);
+ if (result > 0) _currentEntry = entry;
+ return result;
+}
+Long64_t ZDCTreeAnalysis::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+ if (!fChain) return -5;
+ Long64_t centry = fChain->LoadTree(entry);
+ if (centry < 0) return centry;
+ if (fChain->GetTreeNumber() != fCurrent) {
+ fCurrent = fChain->GetTreeNumber();
+ Notify();
+ }
+ return centry;
+}
+
+void ZDCTreeAnalysis::Init(TTree *tree)
+{
+ // The Init() function is called when the selector needs to initialize
+ // a new tree or chain. Typically here the branch addresses and branch
+ // pointers of the tree will be set.
+ // It is normally not necessary to make changes to the generated
+ // code, but the routine can be extended by the user if needed.
+ // Init() will be called many times when running on PROOF
+ // (once per file to be processed).
+
+ // Set object pointer
+ // decisions = 0;
+ // prescales = 0;
+ // Set branch addresses and branch pointers
+ if (!tree) return;
+ fChain = tree;
+ fCurrent = -1;
+ // fChain->SetMakeClass(1);
+
+ fChain->SetBranchAddress("runNumber", &runNumber, &b_runNumber);
+ fChain->SetBranchAddress("eventNumber", &eventNumber, &b_eventNumber);
+ fChain->SetBranchAddress("lumiBlock", &lumiBlock, &b_lumiBlock);
+ fChain->SetBranchAddress("bcid", &bcid, &b_bcid);
+ // fChain->SetBranchAddress("trigger", &trigger, &b_trigger);
+ // fChain->SetBranchAddress("trigger_TBP", &trigger_TBP, &b_trigger_TBP);
+ fChain->SetBranchAddress("tbp", tbp, &b_tbp);
+
+ // fChain->SetBranchAddress("tav", tav, &b_tav);
+ // // fChain->SetBranchAddress("decisions", &decisions, &b_decisions);
+ // // fChain->SetBranchAddress("prescales", &prescales, &b_prescales);
+ fChain->SetBranchAddress("passBits", &passBits, &b_passBits);
+
+ // fChain->SetBranchAddress("zdc_amp", zdc_amp, &b_zdc_amp);
+ // fChain->SetBranchAddress("zdc_amp_rp", zdc_amp_rp, &b_zdc_amp_rp);
+ // fChain->SetBranchAddress("zdc_time", zdc_time, &b_zdc_time);
+ fChain->SetBranchAddress("zdc_raw", zdc_raw, &b_zdc_raw);
+ // fChain->SetBranchAddress("zdc_ampHG", zdc_ampHG, &b_zdc_ampHG);
+ // fChain->SetBranchAddress("zdc_ampLG", zdc_ampLG, &b_zdc_ampLG);
+ // fChain->SetBranchAddress("zdc_sumHG", zdc_sumHG, &b_zdc_sumHG);
+ // fChain->SetBranchAddress("zdc_sumLG", zdc_sumLG, &b_zdc_sumLG);
+ // fChain->SetBranchAddress("zdc_ampHG_rp", zdc_ampHG_rp, &b_zdc_ampHG_rp);
+ // fChain->SetBranchAddress("zdc_ampLG_rp", zdc_ampLG_rp, &b_zdc_ampLG_rp);
+ // fChain->SetBranchAddress("zdc_sumHG_rp", zdc_sumHG_rp, &b_zdc_sumHG_rp);
+ // fChain->SetBranchAddress("zdc_sumLG_rp", zdc_sumLG_rp, &b_zdc_sumLG_rp);
+ // fChain->SetBranchAddress("fcalEt", &fcalEt, &b_fcalEt);
+ // fChain->SetBranchAddress("ntrk", &ntrk, &b_ntrk);
+ // fChain->SetBranchAddress("nvx", &nvx, &b_nvx);
+ // fChain->SetBranchAddress("vx", vx, &b_vx);
+ // fChain->SetBranchAddress("vxntrk", &vxntrk, &b_vxntrk);
+ // fChain->SetBranchAddress("vxcov", vxcov, &b_vxcov);
+ // fChain->SetBranchAddress("mbts_countA", &mbts_countA, &b_mbts_countA);
+ // fChain->SetBranchAddress("mbts_countC", &mbts_countC, &b_mbts_countC);
+ // fChain->SetBranchAddress("mbts_timeA", &mbts_timeA, &b_mbts_timeA);
+ // fChain->SetBranchAddress("mbts_timeC", &mbts_timeC, &b_mbts_timeC);
+ // fChain->SetBranchAddress("mbts_timeDiff", &mbts_timeDiff, &b_mbts_timeDiff);
+
+ fChain->SetBranchAddress("L1_ZDC_A", &L1_ZDC_A, &b_L1_ZDC_A);
+ // fChain->SetBranchAddress("ps_L1_ZDC_A", &ps_L1_ZDC_A, &b_ps_L1_ZDC_A);
+ fChain->SetBranchAddress("L1_ZDC_C", &L1_ZDC_C, &b_L1_ZDC_C);
+ // fChain->SetBranchAddress("ps_L1_ZDC_C", &ps_L1_ZDC_C, &b_ps_L1_ZDC_C);
+ fChain->SetBranchAddress("L1_ZDC_AND", &L1_ZDC_AND, &b_L1_ZDC_AND);
+ // fChain->SetBranchAddress("ps_L1_ZDC_AND", &ps_L1_ZDC_AND, &b_ps_L1_ZDC_AND);
+ fChain->SetBranchAddress("L1_ZDC_A_C", &L1_ZDC_A_C, &b_L1_ZDC_A_C);
+ // fChain->SetBranchAddress("ps_L1_ZDC_A_C", &ps_L1_ZDC_A_C, &b_ps_L1_ZDC_A_C);
+
+ fChain->SetBranchStatus("*", 0);
+ fChain->SetBranchStatus("runNumber", 1);
+ fChain->SetBranchStatus("eventNumber", 1);
+ fChain->SetBranchStatus("lumiBlock", 1);
+ fChain->SetBranchStatus("bcid", 1);
+ fChain->SetBranchStatus("zdc_raw", 1);
+ fChain->SetBranchStatus("tbp", 1);
+
+ fChain->SetBranchStatus("passBits", 1);
+
+ fChain->SetBranchStatus("L1_ZDC_A", 1);
+ fChain->SetBranchStatus("L1_ZDC_C", 1);
+ fChain->SetBranchStatus("L1_ZDC_AND", 1);
+ fChain->SetBranchStatus("L1_ZDC_A_C", 1);
+
+ Notify();
+}
+
+Bool_t ZDCTreeAnalysis::Notify()
+{
+ // The Notify() function is called when a new file is opened. This
+ // can be either for a new TTree in a TChain or when when a new TTree
+ // is started when using PROOF. It is normally not necessary to make changes
+ // to the generated code, but the routine can be extended by the
+ // user if needed. The return value is currently not used.
+
+ return kTRUE;
+}
+
+void ZDCTreeAnalysis::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+ if (!fChain) return;
+ fChain->Show(entry);
+}
+Int_t ZDCTreeAnalysis::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns 1 if entry is accepted.
+// returns -1 otherwise.
+ return 1;
+}
+#endif // #ifdef ZDCTreeAnalysis_cxx
--
GitLab