From f1c44d5e402ec8aea36248cf4381f20b6c034536 Mon Sep 17 00:00:00 2001
From: Peter Alan Steinberg <peter.steinberg@bnl.gov>
Date: Fri, 25 Nov 2016 22:31:58 +0100
Subject: [PATCH] 'added pPb2016 mode to analyze delayed samples. PreSampleAmp
added to output.' (ZdcAnalysis-00-00-17)
---
.../ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx | 131 +++-
.../ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx | 69 ++-
.../ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx | 575 +++++++++++++++++-
.../ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx | 188 +++++-
.../ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h | 10 +
.../ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h | 82 ++-
.../ZdcAnalysis/ZDCPulseAnalyzer.h | 110 +++-
.../ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h | 3 +
.../ZDC/ZdcAnalysis/cmt/Makefile.RootCore | 2 +-
9 files changed, 1083 insertions(+), 87 deletions(-)
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
index 8f989dcfb888..4bf56bd14c9b 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCDataAnalyzer.cxx
@@ -29,6 +29,9 @@ ZDCDataAnalyzer::ZDCDataAnalyzer(int nSample, float deltaTSample, size_t preSamp
_haveECalib(false),
_currentLB(-1)
{
+ _moduleDisabled[0] = {false, false, false, false};
+ _moduleDisabled[1] = {false, false, false, false};
+
_moduleAnalyzers[0] = {0, 0, 0, 0};
_moduleAnalyzers[1] = {0, 0, 0, 0};
@@ -85,6 +88,40 @@ ZDCDataAnalyzer::~ZDCDataAnalyzer()
}
}
+bool ZDCDataAnalyzer::DisableModule(size_t side, size_t module)
+{
+ if (side < 2 && module < 4) {
+ //
+ // Can't disable in the middle of analysis
+ //
+ if (_dataLoaded[side][module]) return false;
+ else {
+ _moduleDisabled[side][module] = true;
+ return true;
+ }
+ }
+ else {
+ return false;
+ }
+}
+
+void ZDCDataAnalyzer::EnableDelayed(float deltaT, const ZDCModuleFloatArray& undelayedDelayedPedestalDiff)
+{
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->EnableDelayed(deltaT, undelayedDelayedPedestalDiff[side][module]);
+ }
+ }
+}
+
+void ZDCDataAnalyzer::SetPeak2ndDerivMinTolerances(size_t tolerance) {
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ _moduleAnalyzers[side][module]->SetPeak2ndDerivMinTolerance(tolerance);
+ }
+ }
+}
+
void ZDCDataAnalyzer::SetTauT0Values(const ZDCModuleBoolArray& fixTau1, const ZDCModuleBoolArray& fixTau2,
const ZDCModuleFloatArray& tau1, const ZDCModuleFloatArray& tau2,
const ZDCModuleFloatArray& t0HG, const ZDCModuleFloatArray& t0LG)
@@ -231,8 +268,18 @@ void ZDCDataAnalyzer::StartEvent(int lumiBlock)
void ZDCDataAnalyzer::LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples)
{
+ // We immediately return if this module is disabled
+ //
+ if (_moduleDisabled[side][module]) {
+ if (_debugLevel > 2) {
+ std::cout << "Skipping analysis of disabled mofule for event index " << _eventCount << ", side, module = " << side << ", " << module << std::endl;
+ }
+
+ return;
+ }
+
if (_debugLevel > 1) {
- std::cout << "Loading data for event index " << _eventCount << ", side, module = " << side << ", " << module << std::endl;
+ std::cout << "/n 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;
@@ -291,18 +338,96 @@ void ZDCDataAnalyzer::LoadAndAnalyzeData(size_t side, size_t module, const std::
_moduleStatus[side][module] = pulseAna_p->GetStatusMask();
}
+void ZDCDataAnalyzer::LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples,
+ const std::vector<float> HGSamplesDelayed, const std::vector<float> LGSamplesDelayed)
+{
+ // We immediately return if this module is disabled
+ //
+ if (_moduleDisabled[side][module]) {
+ if (_debugLevel > 2) {
+ std::cout << "Skipping analysis of disabled mofule for event index " << _eventCount << ", side, module = " << side << ", " << module << std::endl;
+ }
+
+ return;
+ }
+
+ if (_debugLevel > 1) {
+ std::cout << "Loading undelayed and delayed 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;
+ }
+
+ for (size_t sample = 0; sample < HGSamples.size(); sample++) {
+ std::cout << "HGSampleDelayed[" << sample << "] = " << HGSamplesDelayed[sample] << std::endl;
+ }
+ }
+
+ }
+ }
+
+ ZDCPulseAnalyzer* pulseAna_p = _moduleAnalyzers[side][module];
+
+ bool result = pulseAna_p->LoadAndAnalyzeData(HGSamples, LGSamples, HGSamplesDelayed, LGSamplesDelayed);
+ _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;
+ if (!_dataLoaded[side][module] && !_moduleDisabled[side][module]) return false;
}
// Divide the average times by the calibrated module sums
//
- if (_calibModuleSum[side] > 0) {
+ if (_calibModuleSum[side] > 1e-6) {
_averageTime[side] /= _calibModuleSum[side];
}
else {
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
index 0d146a734436..bce7327df6c6 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCFitWrapper.cxx
@@ -20,7 +20,7 @@ ZDCFitExpFermiVariableTaus::ZDCFitExpFermiVariableTaus(std::string tag, float tm
theTF1->SetParLimits(1, tmin, tmax);
theTF1->SetParLimits(2, 2.5, 6);
theTF1->SetParLimits(3, 10, 60);
- theTF1->SetParLimits(4, -0.005, 10);
+ theTF1->SetParLimits(4, -0.05, 0.05);
if (_fixTau1) theTF1->FixParameter(2, _tau1);
if (_fixTau2) theTF1->FixParameter(3, _tau2);
@@ -47,7 +47,7 @@ ZDCFitExpFermiFixedTaus::ZDCFitExpFermiFixedTaus(std::string tag, float tmin, fl
theTF1->SetParLimits(0, 0, 1024);
theTF1->SetParLimits(1, tmin, tmax);
- theTF1->SetParLimits(2, -0.005, 10);
+ theTF1->SetParLimits(2, -0.05, 0.05);
theTF1->SetParName(0, "Amp");
theTF1->SetParName(1, "T0");
@@ -103,14 +103,13 @@ ZDCFitExpFermiPrePulse::ZDCFitExpFermiPrePulse(std::string tag, float tmin, floa
theTF1->SetParLimits(1, tmin, tmax);
theTF1->SetParLimits(2, 0, 1024);
theTF1->SetParLimits(3, -20, 10);
- theTF1->SetParLimits(4, -0.005, 10);
+ theTF1->SetParLimits(4, -0.05, 0.05);
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)
@@ -122,6 +121,68 @@ void ZDCFitExpFermiPrePulse::Initialize(float initialAmp, float initialT0)
GetWrapperTF1()->SetParameter(4, 0);
}
+ZDCFitExpFermiPulseSequence::ZDCFitExpFermiPulseSequence(std::string tag, float tmin, float tmax, float nominalT0, float deltaT, float tau1, float tau2) :
+ ZDCFitWrapper(0), _tau1(tau1), _tau2(tau2)
+{
+ // Create the reference function that we use to evaluate ExpFermiFit more efficiently
+ //
+ std::string funcNameRefFunc = "ExpFermiPulseSequenceRefFunc" + 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);
+
+ // Calculate how many pulses we have to analyze
+ //
+ size_t numPrePulse = std::ceil((nominalT0 - tmin)/deltaT) + 1;
+ for (size_t ipre = 0; ipre < numPrePulse; ipre++) {
+ _pulseDeltaT.push_back(((float) ipre) * -deltaT);
+ }
+
+ size_t numPostPulse = std::floor((tmax - nominalT0)/deltaT);
+ for (size_t ipost = 0; ipost < numPrePulse; ipost++) {
+ _pulseDeltaT.push_back(((float) ipost) * deltaT);
+ }
+
+ _numPulses = 1 + numPrePulse + numPostPulse;
+
+ // Now set up the actual TF1
+ //
+
+ _fitFunc = new TF1(("ExpFermiPulseSequence" + tag).c_str(), this, tmin, tmax, _numPulses + 1),
+
+ _fitFunc->SetParName(0, "Amp");
+ _fitFunc->SetParName(1, "T0");
+
+ _fitFunc->SetParLimits(0, 0, 1024);
+ _fitFunc->SetParLimits(1, nominalT0 - deltaT/2, nominalT0 + deltaT/2);
+
+ for (size_t ipulse = 0; ipulse < _numPulses - 1; ipulse++) {
+ std::string name = "Amp" + std::to_string(ipulse + 1);
+ _fitFunc->SetParName(ipulse + 2, name.c_str());
+ _fitFunc->SetParLimits(ipulse + 2, 0, 1024);
+ }
+}
+
+void ZDCFitExpFermiPulseSequence::Initialize(float initialAmp, float initialT0)
+{
+ _fitFunc->SetParameter(0, initialAmp);
+ _fitFunc->SetParameter(1, initialT0);
+
+ for (size_t ipulse = 0; ipulse < _numPulses - 1; ipulse++) {
+ _fitFunc->SetParameter(ipulse + 2, 0);
+ }
+}
+
+
+
double ZDCFermiExpFit(double* xvec, double* pvec)
{
static const float offsetScale = 0;
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
index 141c2859f98a..e731d3b73916 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZDCPulseAnalyzer.cxx
@@ -9,9 +9,71 @@
#include "TFitResult.h"
#include "TFitResultPtr.h"
+#include "TVirtualFitter.h"
+#include "TFitter.h"
+#include "TList.h"
+#include "TMinuit.h"
+
+
+extern int gErrorIgnoreLevel;
bool ZDCPulseAnalyzer::_quietFits = true;
std::string ZDCPulseAnalyzer::_fitOptions = "";
+TH1* ZDCPulseAnalyzer::_undelayedFitHist = 0;
+TH1* ZDCPulseAnalyzer::_delayedFitHist = 0;
+TF1* ZDCPulseAnalyzer::_combinedFitFunc = 0;
+
+void ZDCPulseAnalyzer::CombinedPulsesFCN(int& numParam, double*, double& f, double* par, int flag)
+{
+ // The first parameter is a correction factor to account for decrease in beam intensity between x
+ // and y scan. It is applied here and not passed to the actual fit function
+ //
+
+ if (flag >=1 && flag <= 3) {
+ f = 0;
+ return;
+ }
+
+ double chiSquare = 0;
+
+ float delayBaselineAdjust = par[0];
+
+ size_t nSamples = _undelayedFitHist->GetNbinsX();
+ if (_delayedFitHist->GetNbinsX() != nSamples) throw;
+
+ // undelayed
+ //
+ for (size_t isample = 0; isample < nSamples; isample++) {
+ double histValue = _undelayedFitHist->GetBinContent(isample + 1);
+ double histError = std::max(_undelayedFitHist->GetBinError(isample + 1), 1.0);
+ double t = _undelayedFitHist->GetBinCenter(isample + 1);
+
+ double funcVal = _combinedFitFunc->EvalPar(&t, &par[1]);
+
+ // std::cout << "Calculating chis^2, undelayed sample " << isample << "t = " << t << ", data = " << histValue << ", func = " << funcVal << std::endl;
+
+ double pull = (histValue - funcVal)/histError;
+ chiSquare += pull*pull;
+ }
+
+ // delayed
+ //
+ for (size_t isample = 0; isample < nSamples; isample++) {
+ double histValue = _delayedFitHist->GetBinContent(isample + 1);
+ double histError = std::max(_delayedFitHist->GetBinError(isample + 1), 1.0);
+ double t = _delayedFitHist->GetBinCenter(isample + 1);
+
+ double funcVal = _combinedFitFunc->EvalPar(&t, &par[1]) + delayBaselineAdjust;
+ double pull = (histValue - funcVal)/histError;
+
+ // std::cout << "Calculating chis^2, delayed sample " << isample << "t = " << t << ", data = " << histValue << ", func = " << funcVal << std::endl;
+
+ chiSquare += pull*pull;
+ }
+
+ f = chiSquare;
+}
+
ZDCPulseAnalyzer::ZDCPulseAnalyzer(std::string tag, int Nsample, float deltaTSample, size_t preSampleIdx, int pedestal,
float gainHG, std::string fitFunction, int peak2ndDerivMinSample,
@@ -19,9 +81,11 @@ ZDCPulseAnalyzer::ZDCPulseAnalyzer(std::string tag, int Nsample, float deltaTSam
_tag(tag), _Nsample(Nsample), _deltaTSample(deltaTSample), _preSampleIdx(preSampleIdx),
_pedestal(pedestal), _gainHG(gainHG), _forceLG(false), _fitFunction(fitFunction),
_peak2ndDerivMinSample(peak2ndDerivMinSample),
+ _peak2ndDerivMinTolerance(1),
_peak2ndDerivMinThreshHG(peak2ndDerivMinThreshHG),
_peak2ndDerivMinThreshLG(peak2ndDerivMinThreshLG),
_haveTimingCorrections(false), _haveNonlinCorr(false), _initializedFits(false),
+ _useDelayed(false), _delayedHist(0), _prePulseCombinedFitter(0), _defaultCombinedFitter(0),
_ADCSamplesHGSub(Nsample, 0), _ADCSamplesLGSub(Nsample, 0),
_ADCSSampSigHG(Nsample, 1), _ADCSSampSigLG(Nsample, 1), // By default the ADC uncertainties are set to one
_samplesSub(Nsample, 0),
@@ -47,6 +111,17 @@ ZDCPulseAnalyzer::~ZDCPulseAnalyzer()
if (!_prePulseFitWrapper) delete _prePulseFitWrapper;
}
+void ZDCPulseAnalyzer::EnableDelayed(float deltaT, float pedestalShift)
+{
+ _useDelayed = true;
+ _delayedDeltaT = deltaT;
+ _delayedPedestalDiff = pedestalShift;
+
+ _delayedHist = new TH1F((std::string(_fitHist->GetName()) + "delayed").c_str(), "", _Nsample, _tmin + _delayedDeltaT, _tmax + _delayedDeltaT);
+
+ _ADCSamplesHGSub.assign(2*_Nsample, 0);
+}
+
void ZDCPulseAnalyzer::SetDefaults()
{
_nominalTau1 = 4;
@@ -97,10 +172,12 @@ void ZDCPulseAnalyzer::Reset()
_badChisq = false;
_fitAmplitude = 0;
- _fitT0 = -100;
+ _fitTime = -100;
+ _fitTimeSub = -100;
+ _fitTimeCorr = -100;
+
_fitChisq = 0;
- _fitT0Corr = -100;
_amplitude = 0;
_ampError = 0;
_preSampleAmp = 0;
@@ -109,6 +186,12 @@ void ZDCPulseAnalyzer::Reset()
_samplesSub.clear();
_samplesDeriv.clear();
_samplesDeriv2nd.clear();
+
+ int sampleVecSize = _Nsample;
+ if (_useDelayed) sampleVecSize *= 2;
+
+ _ADCSamplesHGSub.assign(sampleVecSize, 0);
+ _ADCSamplesLGSub.assign(sampleVecSize, 0);
}
void ZDCPulseAnalyzer::SetTauT0Values(bool fixTau1, bool fixTau2, float tau1, float tau2, float t0HG, float t0LG)
@@ -172,6 +255,9 @@ void ZDCPulseAnalyzer::SetupFitFunctions()
bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG)
{
+ // std::cout << "Use delayed = " << _useDelayed << std::endl;
+
+ if (_useDelayed) throw;
if (!_initializedFits) SetupFitFunctions();
// Clear any transient data
@@ -185,6 +271,8 @@ bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::
return false;
}
+
+
// Now do pedestal subtraction and check for overflows
//
for (size_t isample = 0; isample < _Nsample; isample++) {
@@ -219,6 +307,162 @@ bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::
_ADCSamplesLGSub[isample] = ADCLG - _pedestal;
}
+ // std::cout << "Loaded data" << std::endl;
+
+ // 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(_Nsample, _preSampleIdx, _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(_Nsample, _preSampleIdx, _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::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG,
+ std::vector<float> ADCSamplesHGDelayed, std::vector<float> ADCSamplesLGDelayed)
+{
+ if (!_useDelayed) throw;
+ 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 ||
+ ADCSamplesHGDelayed.size() != _Nsample || ADCSamplesLGDelayed.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];
+
+ float ADCHGDelay = ADCSamplesHGDelayed[isample];
+ float ADCLGDelay = ADCSamplesLGDelayed[isample];
+
+ if (ADCHG > _HGOverflowADC) {
+ _HGOverflow = true;
+
+ if (isample == _preSampleIdx) _PSHGOverUnderflow = true;
+ }
+ else if (ADCHG < _HGUnderflowADC) {
+ _HGUnderflow = true;
+
+ if (isample == _preSampleIdx) _PSHGOverUnderflow = true;
+ }
+
+ if (ADCHGDelay > _HGOverflowADC) {
+ _HGOverflow = true;
+ }
+ else if (ADCHGDelay < _HGUnderflowADC) {
+ _HGUnderflow = true;
+ }
+
+ if (ADCLG > _LGOverflowADC) {
+ _LGOverflow = true;
+ _fail = true;
+ _amplitude = 1024 * _gainHG; // Give a value 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
+ }
+ else if (ADCLG == 0) {
+ _LGUnderflow = true;
+ _fail = true;
+ }
+
+ if (ADCLGDelay > _LGOverflowADC) {
+ _LGOverflow = true;
+ _fail = true;
+ _amplitude = 1024 * _gainHG; // Give a value 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
+ }
+ else if (ADCLGDelay == 0) {
+ _LGUnderflow = true;
+ _fail = true;
+ }
+
+ _ADCSamplesHGSub[isample*2] = ADCHG - _pedestal;
+ _ADCSamplesLGSub[isample*2] = ADCLG - _pedestal;
+
+ _ADCSamplesHGSub[isample*2 + 1] = ADCHGDelay - _pedestal - _delayedPedestalDiff;
+ _ADCSamplesLGSub[isample*2 + 1] = ADCLGDelay - _pedestal - _delayedPedestalDiff;
+ }
+
+ // if (!_PSHGOverUnderflow) {
+
+ // std::cout << "_ADCSamplesHGSub = ";
+ // for (size_t sample = 0; sample < _ADCSamplesHGSub.size(); sample++) {
+ // std::cout << ", [" << sample << "] = " << _ADCSamplesHGSub[sample];
+ // }
+
+ // // Attempt to address up front cases where we have significant offsets between the delayed and undelayed
+ // //
+ // float delta10 = _ADCSamplesHGSub[1] - _ADCSamplesHGSub[0];
+ // float delta21 = _ADCSamplesHGSub[2] - _ADCSamplesHGSub[1];
+ // float delta32 = _ADCSamplesHGSub[3] - _ADCSamplesHGSub[2];
+
+ // std::cout << "testing delayed shift: delta10 = " << delta10 << ", delta21 = " << delta21 << ", delta32 = " << delta32 << std::endl;
+
+ // if (std::abs(delta10 + delta21) <= 2 && std::abs(delta10 - delta32) <= 2) {
+ // // Very likely we have an offset between the delayed and undelayed
+ // //
+
+ // std::cout << "correcting HG delayed ssamples down by " << 0.5*(delta10 + delta32) << std::endl;
+
+ // for (int isample = 0; isample < _Nsample; isample++) {
+ // _ADCSamplesHGSub[isample*2 + 1] -= 0.5*(delta10 + delta32);
+ // }
+
+ // std::cout << " corrected _ADCSamplesHGSub = ";
+ // for (size_t sample = 0; sample < _ADCSamplesHGSub.size(); sample++) {
+ // std::cout << ", [" << sample << "] = " << _ADCSamplesHGSub[sample];
+ // }
+
+ // }
+ // }
+
// if (_fail) return false;
// Now we actually do the analysis using the high or low gain samples as determined
@@ -227,7 +471,8 @@ bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::
_useLowGain = _HGUnderflow || _HGOverflow || _forceLG;
if (_useLowGain) {
- bool result = AnalyzeData(_ADCSamplesLGSub, _ADCSSampSigLG, _LGT0CorrParams, _chisqDivAmpCutLG,
+
+ bool result = AnalyzeData(_Nsample*2, _preSampleIdx*2, _ADCSamplesLGSub, _ADCSSampSigLG, _LGT0CorrParams, _chisqDivAmpCutLG,
_T0CutLowLG, _T0CutHighLG, _peak2ndDerivMinThreshLG);
if (result) {
//
@@ -241,7 +486,7 @@ bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::
return result;
}
else {
- bool result = AnalyzeData(_ADCSamplesHGSub, _ADCSSampSigHG, _HGT0CorrParams, _chisqDivAmpCutHG,
+ bool result = AnalyzeData(_Nsample*2, _preSampleIdx*2, _ADCSamplesHGSub, _ADCSSampSigHG, _HGT0CorrParams, _chisqDivAmpCutHG,
_T0CutLowHG, _T0CutHighHG, _peak2ndDerivMinThreshHG);
if (result) {
_preSampleAmp = _preSample;
@@ -263,7 +508,8 @@ bool ZDCPulseAnalyzer::LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::
}
}
-bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, // The samples used for this event
+bool ZDCPulseAnalyzer::AnalyzeData(size_t nSamples, size_t preSampleIdx,
+ 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
@@ -273,10 +519,44 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
{
// Do the presample subtraction using a lambda function
//
- float preSample = samples[_preSampleIdx];
+ float preSample = samples[preSampleIdx];
_preSample = preSample;
_samplesSub = samples;
+
+ // std::cout << "_samplesSub = ";
+ // for (size_t sample = 0; sample < _samplesSub.size(); sample++) {
+ // std::cout << ", [" << sample << "] = " << _samplesSub[sample];
+ // }
+
+ if (_useDelayed) {
+ // Attempt to address up front cases where we have significant offsets between the delayed and undelayed
+ //
+ // float delta10 = _samplesSub[1] - _samplesSub[0];
+ // float delta21 = _samplesSub[2] - _samplesSub[1];
+ // float delta32 = _samplesSub[3] - _samplesSub[2];
+
+ // // std::cout << "testing delayed shift: delta10 = " << delta10 << ", delta21 = " << delta21 << ", delta32 = " << delta32 << std::endl;
+
+ // if (std::abs(delta10 + delta21) <= 0.1*std::abs(delta10 - delta21)/2 && std::abs(delta10 - delta32) <= 0.1*std::abs(delta10 + delta32)/2) {
+ // // Very likely we have an offset between the delayed and undelayed
+ // //
+
+ // std::cout << "correcting HG delayed ssamples down by " << 0.5*(delta10 + delta32) << std::endl;
+
+
+ double baselineCorr = 0.5*(_samplesSub[1] - _samplesSub[0] + _samplesSub[3] - _samplesSub[4]);
+
+ for (int isample = 0; isample < nSamples; isample++) {
+ if (isample%2) _samplesSub[isample] -= baselineCorr;
+ }
+
+ // std::cout << " corrected _samplesSub = ";
+ // for (size_t sample = 0; sample < _samplesSub.size(); sample++) {
+ // std::cout << ", [" << sample << "] = " << _samplesSub[sample];
+ // }
+ }
+
std::for_each(_samplesSub.begin(), _samplesSub.end(), [=] (float& adcUnsub) {return adcUnsub -= preSample;} );
// Find maximum and minimum values
@@ -295,14 +575,14 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
// Calculate first and 2nd "derivatives"
//
- std::vector<float> deriv(_Nsample, 0);
- std::adjacent_difference(samples.begin(), samples.end(), deriv.begin());
+ std::vector<float> deriv(nSamples, 0);
+ std::adjacent_difference(_samplesSub.begin(), _samplesSub.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::vector<float> deriv2nd(nSamples, 0);
std::adjacent_difference(_samplesDeriv.begin(), _samplesDeriv.end(), deriv2nd.begin());
// Save the second derivatives, dropping the useless first element
@@ -319,7 +599,7 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
// 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) {
+ if (std::abs(_minDeriv2ndIndex - (int) _peak2ndDerivMinSample) <= _peak2ndDerivMinTolerance && _minDeriv2nd <= peak2ndDerivMinThresh) {
_havePulse = true;
}
else {
@@ -328,34 +608,55 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
// Now decide whether we have a preceeding pulse or not
//
+ // std::cout << "_samplesSDeriv = ";
+ // for (size_t sample = 0; sample < _samplesDeriv.size(); sample++) {
+ // std::cout << ", [" << sample << "] = " << _samplesDeriv[sample];
+ // }
+
// Have to be careful with the derivative indices, but test up to minDeriv2ndIndex - 1
//
- bool haveNegDeriv = false;
+ bool lastNegDeriv = 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) {
+ if (_samplesDeriv[isampl] < -0.5) {
//
// 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;
+ if (std::abs(_samplesDeriv[isampl]) / (_maxADCValue - _minADCValue + 1) > 0.05) {
+
+ // std:: cout << "Setting prepulse due to large negative derivative: sample = " << isampl << ", value = " << _samplesDeriv[isampl]
+ // << ", threshold = " << (_maxADCValue - _minADCValue + 1)*0.05 << std::endl;
+ _prePulse = true;
+ }
// Also, if we have two negative derivatives in a row, we have a prepulse
//
- if (haveNegDeriv) _prePulse = true;
- else haveNegDeriv = true;
+ if (lastNegDeriv) {
+ _prePulse = true;
+ // std:: cout << "Setting prepulse due to successive negative derivatives" << std::endl;
+ }
+ else lastNegDeriv = true;
+ }
+ else {
+ lastNegDeriv = false;
}
// 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;
+ if (_samplesDeriv2nd[isampl] < 0 && _samplesDeriv2nd[isampl + 1] > 0) {
+ if (_samplesDeriv2nd[isampl] < 0.1*_minDeriv2nd) {
+ _prePulse = true;
+ // std:: cout << "Setting prepulse due to negative 2nd derivative: sample = " << isampl << ", value = " << _samplesDeriv2nd[isampl] << ", threshold = " << 0.1*_minDeriv2nd << std::endl;
+ }
+ }
}
// Additional specific tests for pre pulses etc.
@@ -365,8 +666,8 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
// 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;
+ //if //_samplesSub[_Nsample - 1] <= -_maxADCValue*0.05 ||
+ // if (_samplesSub[_Nsample - 1] - _samplesSub[_preSampleIdx+1] < -_maxADCValue*0.05) _prePulse = true;
}
@@ -377,10 +678,11 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
//
if (_fail || !_havePulse) return false;
- DoFit();
+ if (!_useDelayed) DoFit();
+ else DoFitCombined();
if (!FitFailed()) {
- _fitT0Corr = _fitT0;
+ _fitTimeCorr = _fitTimeSub;
if (_haveTimingCorrections) {
@@ -392,14 +694,14 @@ bool ZDCPulseAnalyzer::AnalyzeData(const std::vector<float>& samples, //
float correction = (t0CorrParams[0] + t0CorrParams[1]*t0CorrFact + t0CorrParams[2]*t0CorrFact*t0CorrFact +
t0CorrParams[3]*t0CorrFact*t0CorrFact*t0CorrFact);
- _fitT0Corr -= correction;
+ _fitTimeCorr -= correction;
}
// Now check for valid chisq and valid time
//
- if (_fitChisq / _fitAmplitude > maxChisqDivAmp) _badChisq = true;
- if (_fitT0Corr < minT0Corr || _fitT0Corr > maxT0Corr) _badT0 = true;
+ if (_fitChisq / (_fitAmplitude + 1.0e-6) > maxChisqDivAmp) _badChisq = true;
+ if (_fitTimeCorr < minT0Corr || _fitTimeCorr > maxT0Corr) _badT0 = true;
// std::cout << "Testing t0 cut: t0corr = " << _fitT0Corr << ", min max cuts = " << minT0Corr << ", " << maxT0Corr << std::endl;
}
@@ -440,13 +742,14 @@ void ZDCPulseAnalyzer::DoFit()
// if (!_fitFailed) {
_bkgdMaxFraction = fitWrapper->GetBkgdMaxFraction();
if (std::abs(_bkgdMaxFraction) > 0.25) {
- options += "e";
- _fitHist->Fit(fitWrapper->GetWrapperTF1(), options.c_str());
+ std::string tempOptions = options + "e";
+ _fitHist->Fit(fitWrapper->GetWrapperTF1(), tempOptions.c_str());
_bkgdMaxFraction = fitWrapper->GetBkgdMaxFraction();
}
_fitAmplitude = fitWrapper->GetAmplitude();
- _fitT0 = fitWrapper->GetTime() - t0Initial;
+ _fitTime = fitWrapper->GetTime();
+ _fitTimeSub = _fitTime - t0Initial;
_fitChisq = result_ptr->Chi2();
_fitTau1 = fitWrapper->GetTau1();
@@ -456,12 +759,230 @@ void ZDCPulseAnalyzer::DoFit()
// }
}
+void ZDCPulseAnalyzer::DoFitCombined()
+{
+ // 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()) {
+ if (_samplesDeriv[0] < 0) {
+ (static_cast<ZDCPrePulseFitWrapper*>(_prePulseFitWrapper))->SetInitialPrePulseT0(-10);
+ }
+ else {
+ (static_cast<ZDCPrePulseFitWrapper*>(_prePulseFitWrapper))->SetInitialPrePulseT0(0);
+ }
+ }
+
+ // Set up the virtual fitter
+ //
+ TFitter* theFitter = 0;
+
+ if (PrePulse()) {
+ if (!_prePulseCombinedFitter) _prePulseCombinedFitter = MakeCombinedFitter(fitWrapper->GetWrapperTF1());
+ theFitter = _prePulseCombinedFitter;
+
+ // std::cout << "Testing prepulse wrapper, npar = " << theFitter->GetNumberTotalParameters() << std::endl;
+ }
+ else {
+ if (!_defaultCombinedFitter) {
+ _defaultCombinedFitter = MakeCombinedFitter(fitWrapper->GetWrapperTF1());
+ }
+ theFitter = _defaultCombinedFitter;
+
+ // std::cout << "Testing default wrapper, npar = " << theFitter->GetNumberTotalParameters() << std::endl;
+ }
+
+ // Set the static pointers to histograms and function for use in FCN
+ //
+ _undelayedFitHist = _fitHist;
+ _delayedFitHist = _delayedHist;
+ _combinedFitFunc = fitWrapper->GetWrapperTF1();
+
+ size_t numFitPar = theFitter->GetNumberTotalParameters();
+
+ theFitter->GetMinuit()->fISW[4] = -1;
+
+ double arglist[100];
+ // arglist[0] = 1;
+ // arglist[1] = 0;
+
+ // theFitter->ExecuteCommand("SET PAR ", arglist, 2);
+
+ theFitter->SetParameter(0, "delayBaselineAdjust", 0, 0.01, 0, 0);
+
+ for (size_t ipar = 0; ipar < numFitPar - 1; ipar++) {
+ double parLimitLow, parLimitHigh;
+
+ _combinedFitFunc->GetParLimits(ipar, parLimitLow, parLimitHigh);
+
+ // std::cout << "Setting minuit parameter " << ipar + 1 << " to value " << _combinedFitFunc->GetParameter(ipar) << ", limits = "
+ // << parLimitLow << " to " << parLimitHigh << std::endl;
+
+ theFitter->SetParameter(ipar + 1, _combinedFitFunc->GetParName(ipar), _combinedFitFunc->GetParameter(ipar), 0.01, parLimitLow, parLimitHigh);
+
+ // arglist[0] = ipar + 2;;
+ // arglist[1] = _combinedFitFunc->GetParameter(ipar);
+ // theFitter->ExecuteCommand("SET PAR ", arglist, 2);
+
+ // theFitter->GetMinuit()->Command(("SET PAR " + std::to_string(ipar + 2) + " " + std::to_string(_combinedFitFunc->GetParameter(ipar))).c_str() );
+ }
+
+ // Now perform the fit
+ //
+ // double arglist[100];
+
+ if (_quietFits) {
+ theFitter->GetMinuit()->fISW[4] = -1;
+
+ int ierr= 0;
+ theFitter->GetMinuit()->mnexcm("SET NOWarnings",0,0,ierr);
+ }
+ else theFitter->GetMinuit()->fISW[4] = 0;
+
+ arglist[0] = 2000; // number of function calls
+ arglist[1] = 0.01; // tolerance
+ int status = theFitter->ExecuteCommand("MIGRAD",arglist,2);
+
+ double fitAmp = theFitter->GetParameter(1);
+ if (status || fitAmp < 1) {
+
+ // We first fit with no baseline adjust
+ //
+ theFitter->SetParameter(0, "delayBaselineAdjust", 0, 0.01, 0, 0);
+ theFitter->FixParameter(0);
+
+ for (size_t ipar = 0; ipar < numFitPar - 1; ipar++) {
+ double parLimitLow, parLimitHigh;
+
+ _combinedFitFunc->GetParLimits(ipar, parLimitLow, parLimitHigh);
+ theFitter->SetParameter(ipar + 1, _combinedFitFunc->GetParName(ipar), _combinedFitFunc->GetParameter(ipar), 0.01, parLimitLow, parLimitHigh);
+ }
+
+ arglist[0] = 2000; // number of function calls
+ arglist[1] = 0.01; // tolerance
+ status = theFitter->ExecuteCommand("MIGRAD",arglist,2);
+ if (!status) {
+ theFitter->ReleaseParameter(0);
+ _fitFailed = true;
+ }
+ else {
+ theFitter->ReleaseParameter(0);
+ status = theFitter->ExecuteCommand("MIGRAD",arglist,2);
+
+ if (status || fitAmp < 1) {_fitFailed = false;}
+ }
+ }
+ else _fitFailed = false;
+
+ if (!_quietFits) theFitter->GetMinuit()->fISW[4] = -1;
+
+ std::vector<double> funcParams(numFitPar - 1);
+ std::vector<double> funcParamErrs(numFitPar - 1);
+
+ // Save the baseline shift between delayed and undelayed samples
+ //
+ _delayedBaselineShift = theFitter->GetParameter(0);
+
+ // Capture and store the fit function parameteds and errors
+ //
+ for (size_t ipar = 1; ipar < numFitPar; ipar++) {
+ funcParams[ipar - 1] = theFitter->GetParameter(ipar);
+ funcParamErrs[ipar - 1] = theFitter->GetParError(ipar);
+ }
+
+ _combinedFitFunc->SetParameters(&funcParams[0]);
+ _combinedFitFunc->SetParErrors(&funcParamErrs[0]);
+
+ // Cpature the chi-square etc.
+ //
+ double chi2, edm, errdef;
+ int nvpar, nparx;
+
+ theFitter->GetStats(chi2, edm, errdef, nvpar, nparx);
+
+ int ndf = 2*_Nsample - nvpar;
+
+ _combinedFitFunc->SetChisquare(chi2);
+ _combinedFitFunc->SetNDF(ndf);
+
+ // add to list of functions
+ _undelayedFitHist->GetListOfFunctions()->Add(_combinedFitFunc);
+ _delayedFitHist->GetListOfFunctions()->Add(_combinedFitFunc);
+
+ // if (QuietFits()) options += "Q0";
+
+ // 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) {
+ // // std::string tempOptions = options + "e";
+ // _fitHist->Fit(fitWrapper->GetWrapperTF1(), tempOptions.c_str());
+ // _bkgdMaxFraction = fitWrapper->GetBkgdMaxFraction();
+ // }
+
+ _fitAmplitude = fitWrapper->GetAmplitude();
+ _fitTime = fitWrapper->GetTime();
+ _fitTimeSub = _fitTime - t0Initial;
+ _fitChisq = chi2;
+
+ _fitTau1 = fitWrapper->GetTau1();
+ _fitTau2 = fitWrapper->GetTau2();
+
+ _fitAmpError = fitWrapper->GetAmpError();
+ // }
+}
+
+
+TFitter* ZDCPulseAnalyzer::MakeCombinedFitter(TF1* func)
+{
+ // gErrorIgnoreLevel = 1001;
+ TVirtualFitter::SetDefaultFitter("Minuit");
+
+ size_t nFitParams = func->GetNpar() + 1;
+ TFitter* fitter = new TFitter(nFitParams);
+
+ fitter->GetMinuit()->fISW[4] = -1;
+ fitter->SetParameter(0, "delayBaselineAdjust", 0, 0.01, 0, 0);
+
+ for (size_t ipar = 0; ipar < nFitParams - 1; ipar++) {
+ double parLimitLow, parLimitHigh;
+
+ func->GetParLimits(ipar, parLimitLow, parLimitHigh);
+ fitter->SetParameter(ipar + 1, func->GetParName(ipar), func->GetParameter(ipar), 0.01, parLimitLow, parLimitHigh);
+ }
+
+ fitter->SetFCN(ZDCPulseAnalyzer::CombinedPulsesFCN);
+
+ double arglist[100];
+ arglist[0] = -1;
+
+ return fitter;
+}
+
void ZDCPulseAnalyzer::Dump() const
{
std::cout << "ZDCPulseAnalyzer dump for tag = " << _tag << std::endl;
std::cout << "Presample index, value = " << _preSampleIdx << ", " << _preSample << std::endl;
+ if (_useDelayed) {
+ std::cout << "using delayed samples with delta T = " << _delayedDeltaT << ", and pedestalDiff == " << _delayedPedestalDiff << std::endl;
+ }
+
std::cout << "samplesSub ";
for (size_t sample = 0; sample < _samplesSub.size(); sample++) {
std::cout << ", [" << sample << "] = " << _samplesSub[sample];
@@ -497,8 +1018,8 @@ unsigned int ZDCPulseAnalyzer::GetStatusMask() const
if (LGOverflow()) statusMask |= 1 << LGOverflowBit;
if (LGUnderflow()) statusMask |= 1 << LGUnderflowBit;
- if (FitFailed()) statusMask |= 1 << FitFailedBit;
if (PrePulse()) statusMask |= 1 << PrePulseBit;
+ if (FitFailed()) statusMask |= 1 << FitFailedBit;
if (BadChisq()) statusMask |= 1 << BadChisqBit;
if (BadT0()) statusMask |= 1 << BadT0Bit;
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
index 535aa78f78b8..cdd0a11d1514 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/Root/ZdcAnalysisTool.cxx
@@ -48,6 +48,9 @@ namespace ZDC
declareProperty("NumSampl", m_numSample = 7);
declareProperty("DeltaTSample", m_deltaTSample = 25);
declareProperty("Presample", m_presample = 0);
+ declareProperty("CombineDelay", m_combineDelay = false);
+ declareProperty("DelayDeltaT", m_delayDeltaT = -12.5);
+
declareProperty("PeakSample", m_peakSample = 2);
declareProperty("Peak2ndDerivThresh", m_Peak2ndDerivThresh = 10);
@@ -230,9 +233,89 @@ ZDCDataAnalyzer* ZdcAnalysisTool::initializeDefault()
zdcDataAnalyzer->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1, tau2, t0, t0);
zdcDataAnalyzer->SetCutValues(chisqDivAmpCut, chisqDivAmpCut, deltaT0CutLow, deltaT0CutHigh, deltaT0CutLow, deltaT0CutHigh);
+ if (m_combineDelay) {
+ ZDCDataAnalyzer::ZDCModuleFloatArray defaultPedestalShifts = {0, 0, 0, 0, 0, 0, 0, 0};
+
+ zdcDataAnalyzer->EnableDelayed(m_delayDeltaT, defaultPedestalShifts);
+ }
+
return zdcDataAnalyzer;
}
+ZDCDataAnalyzer* ZdcAnalysisTool::initializepPb2016()
+{
+ //
+ // 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 tau1Arr, tau2Arr, peak2ndDerivMinSamples, t0;
+ ZDCDataAnalyzer::ZDCModuleFloatArray peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG;
+ ZDCDataAnalyzer::ZDCModuleFloatArray deltaT0CutLow, deltaT0CutHigh, chisqDivAmpCut;
+ ZDCDataAnalyzer::ZDCModuleBoolArray fixTau1Arr, fixTau2Arr;
+
+ // For now we allow the tau values to be controlled by the job properties until they are better determined
+ //
+ const int peakSample = 5;
+ const float peak2ndDerivThreshHG = -12;
+ const float peak2ndDerivThreshLG = -10;
+ const float tau1 = 4.5;
+ const float tau2 = 22.;
+
+
+ for (size_t side : {0, 1}) {
+ for (size_t module : {0, 1, 2, 3}) {
+ fixTau1Arr[side][module] = m_fixTau1;
+ fixTau2Arr[side][module] = m_fixTau2;
+ tau1Arr[side][module] = tau1;
+ tau2Arr[side][module] = tau2;
+
+ peak2ndDerivMinSamples[side][module] = peakSample;
+ peak2ndDerivMinThresholdsHG[side][module] = peak2ndDerivThreshHG;
+ peak2ndDerivMinThresholdsLG[side][module] = peak2ndDerivThreshLG;
+
+ 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
+ //
+ // We adopt hard-coded values for the number of samples and the frequency which we kept fixed for all physics data
+ //
+ ZDCDataAnalyzer* zdcDataAnalyzer = new ZDCDataAnalyzer(7, 25, 1, "FermiExp", peak2ndDerivMinSamples,
+ peak2ndDerivMinThresholdsHG, peak2ndDerivMinThresholdsLG, m_lowGainOnly);
+
+ // Open up tolerances on the position of the peak for now
+ //
+ zdcDataAnalyzer->SetPeak2ndDerivMinTolerances(2);
+
+ // We alwyas disable the 12EM (sideC) module which was not present (LHCf)
+ //
+ zdcDataAnalyzer->DisableModule(0,0);
+
+ zdcDataAnalyzer->SetADCOverUnderflowValues(HGOverFlowADC, HGUnderFlowADC, LGOverFlowADC);
+ zdcDataAnalyzer->SetTauT0Values(fixTau1Arr, fixTau2Arr, tau1Arr, tau2Arr, t0, t0);
+ zdcDataAnalyzer->SetCutValues(chisqDivAmpCut, chisqDivAmpCut, deltaT0CutLow, deltaT0CutHigh, deltaT0CutLow, deltaT0CutHigh);
+
+ // We allow the combineDelay to be controlled by the properties
+ //
+ // if (m_combineDelay) {
+ m_combineDelay = true;
+ ZDCDataAnalyzer::ZDCModuleFloatArray defaultPedestalShifts = {0, 0, 0, 0, 0, 0, 0, 0};
+
+ zdcDataAnalyzer->EnableDelayed(-12.5, defaultPedestalShifts);
+ // }
+
+ return zdcDataAnalyzer;
+}
void ZdcAnalysisTool::initialize40MHz()
{
@@ -437,26 +520,6 @@ StatusCode ZdcAnalysisTool::initializeTool()
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);
@@ -478,6 +541,9 @@ StatusCode ZdcAnalysisTool::initializeTool()
m_zdcDataAnalyzer = m_zdcDataAnalyzer_80MHz; // default
}
+ else if (m_configuration == "pPb2016") {
+ m_zdcDataAnalyzer = initializepPb2016();
+ }
else {
ATH_MSG_ERROR("Unknown configuration: " << m_configuration);
return StatusCode::FAILURE;
@@ -489,6 +555,32 @@ StatusCode ZdcAnalysisTool::initializeTool()
m_auxSuffix = "_" + m_auxSuffix;
}
+ ATH_MSG_INFO("Configuration: "<<m_configuration);
+ ATH_MSG_INFO("FlipEMDelay: "<<m_flipEMDelay);
+ ATH_MSG_INFO("LowGainOnly: "<<m_lowGainOnly);
+
+ ATH_MSG_INFO("Using Combined delayed and undelayed samples: "<< m_combineDelay);
+
+ 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);
+
+ if (m_combineDelay) ATH_MSG_INFO("DelayDeltaT: "<< m_delayDeltaT);
+
+ 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);
m_init = true;
return StatusCode::SUCCESS;
@@ -608,8 +700,11 @@ StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& modul
ATH_MSG_DEBUG("Starting event processing");
m_zdcDataAnalyzer->StartEvent(calibLumiBlock);
- const std::vector<unsigned short>* adc0;
- const std::vector<unsigned short>* adc1;
+ const std::vector<unsigned short>* adcUndelayLG = 0;
+ const std::vector<unsigned short>* adcUndelayHG = 0;
+
+ const std::vector<unsigned short>* adcDelayLG = 0;
+ const std::vector<unsigned short>* adcDelayHG = 0;
ATH_MSG_DEBUG("Processing modules");
for (const auto zdcModule : moduleContainer)
@@ -619,27 +714,56 @@ StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& modul
if (zdcModule->zdcModule()==0 && m_flipEMDelay) // flip delay/non-delay for EM big tube
{
- adc0 = &(*(zdcModule->TTg0d1Link()))->adc();
- adc1 = &(*(zdcModule->TTg1d1Link()))->adc();
+ adcUndelayLG = &(*(zdcModule->TTg0d1Link()))->adc();
+ adcUndelayHG = &(*(zdcModule->TTg1d1Link()))->adc();
+
+ adcDelayLG = &(*(zdcModule->TTg0d0Link()))->adc();
+ adcDelayHG = &(*(zdcModule->TTg1d0Link()))->adc();
}
else
{
- adc0 = &(*(zdcModule->TTg0d0Link()))->adc();
- adc1 = &(*(zdcModule->TTg1d0Link()))->adc();
+ adcUndelayLG = &(*(zdcModule->TTg0d0Link()))->adc();
+ adcUndelayHG = &(*(zdcModule->TTg1d0Link()))->adc();
+
+ adcDelayLG = &(*(zdcModule->TTg0d1Link()))->adc();
+ adcDelayHG = &(*(zdcModule->TTg1d1Link()))->adc();
}
- static std::vector<float> HGADCSamples(m_numSample);
- static std::vector<float> LGADCSamples(m_numSample);
+ static std::vector<float> HGUndelADCSamples(m_numSample);
+ static std::vector<float> LGUndelADCSamples(m_numSample);
- std::copy(adc0->begin(),adc0->end(),LGADCSamples.begin());
- std::copy(adc1->begin(),adc1->end(),HGADCSamples.begin());
+ std::copy(adcUndelayLG->begin(), adcUndelayLG->end(), LGUndelADCSamples.begin());
+ std::copy(adcUndelayHG->begin(), adcUndelayHG->end(), HGUndelADCSamples.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);
+ if (!m_combineDelay) {
+ m_zdcDataAnalyzer->LoadAndAnalyzeData(side,zdcModule->zdcModule(), HGUndelADCSamples, LGUndelADCSamples);
+ }
+ else {
+ std::vector<float> HGDelayADCSamples(m_numSample);
+ std::vector<float> LGDelayADCSamples(m_numSample);
+
+ std::copy(adcDelayLG->begin(),adcDelayLG->end(), LGDelayADCSamples.begin());
+ std::copy(adcDelayHG->begin(),adcDelayHG->end(), HGDelayADCSamples.begin());
+
+ // If the delayed channels actually come earlier (as in the pPb in 2016), we invert the meaning of delayed and undelayed
+ // see the initialization sections for similar inversion on the sign of the pedestal difference
+ //
+ if (m_delayDeltaT > 0) {
+ m_zdcDataAnalyzer->LoadAndAnalyzeData(side,zdcModule->zdcModule(),
+ HGUndelADCSamples, LGUndelADCSamples,
+ HGDelayADCSamples, LGDelayADCSamples);
+ }
+ else {
+ m_zdcDataAnalyzer->LoadAndAnalyzeData(side,zdcModule->zdcModule(),
+ HGDelayADCSamples, LGDelayADCSamples,
+ HGUndelADCSamples, LGUndelADCSamples);
+ }
+ }
}
ATH_MSG_DEBUG("Finishing event processing");
@@ -678,6 +802,8 @@ StatusCode ZdcAnalysisTool::recoZdcModules(const xAOD::ZdcModuleContainer& modul
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();
+ zdcModule->auxdecor<float>("PreSampleAmp" + m_auxSuffix) = pulseAna_p->GetPreSampleAmp();
+ //zdcModule->auxdecor<float>("Presample" + m_auxSuffix) = pulseAna_p->GetPresample();
}
/*
std::cout << "side = " << side << " module=" << zdcModule->zdcModule() << " CalibEnergy=" << zdcModule->auxdecor<float>("CalibEnergy")
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
index 9bcafdd3cd29..e45dc44efcdc 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCDataAnalyzer.h
@@ -24,6 +24,7 @@ private:
std::string _fitFunction;
bool _forceLG;
+ ZDCModuleBoolArray _moduleDisabled;
std::array<std::array<ZDCPulseAnalyzer*, 4>, 2> _moduleAnalyzers;
int _debugLevel;
@@ -81,6 +82,8 @@ public:
else ZDCPulseAnalyzer::SetQuietFits(false);
}
+ void EnableDelayed(float deltaT, const ZDCModuleFloatArray& undelayedDelayedPedestalDiff);
+
unsigned int GetModuleMask() const {return _moduleMask;}
float GetModuleSum(size_t side) const {return _moduleSum.at(side);}
@@ -104,6 +107,10 @@ public:
const ZDCPulseAnalyzer* GetPulseAnalyzer(size_t side, size_t module) const {return _moduleAnalyzers.at(side).at(module);}
+ bool DisableModule(size_t side, size_t module);
+
+ void SetPeak2ndDerivMinTolerances(size_t tolerance);
+
void SetADCOverUnderflowValues(const ZDCModuleFloatArray& HGOverflowADC, const ZDCModuleFloatArray& HGUnderflowADC,
const ZDCModuleFloatArray& LGOverflowADC);
@@ -147,6 +154,9 @@ public:
void LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples);
+ void LoadAndAnalyzeData(size_t side, size_t module, const std::vector<float> HGSamples, const std::vector<float> LGSamples,
+ const std::vector<float> HGSamplesDelayed, const std::vector<float> LGSamplesDelayed);
+
bool FinishEvent();
};
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
index e24e80bb0ee3..e4760b15688a 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCFitWrapper.h
@@ -190,8 +190,8 @@ public:
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 GetTau1() const {return _tau1;}
+ virtual float GetTau2() const {return _tau2;}
virtual float GetTime() const {
float fitT0 = GetWrapperTF1()->GetParameter(1);
@@ -252,4 +252,82 @@ public:
}
};
+class ZDCFitExpFermiPulseSequence : public ZDCFitWrapper
+{
+ float _tau1;
+ float _tau2;
+ float _norm;
+ float _timeCorr;
+
+ size_t _numPulses;
+ std::vector<float> _pulseDeltaT;
+
+ TF1* _fitFunc;
+ TF1* _expFermiFunc;
+
+
+public:
+ ZDCFitExpFermiPulseSequence(std::string tag, float tmin, float tmax, float nominalT0, float deltaT, float tau1, float tau2);
+
+ ~ZDCFitExpFermiPulseSequence()
+ {
+ delete _fitFunc;
+ delete _expFermiFunc;
+ }
+
+ virtual TF1* GetWrapperTF1() {return _fitFunc;}
+ virtual const TF1* GetWrapperTF1() const {return _fitFunc;}
+
+ virtual void Initialize(float initialAmp, float initialT0);
+ // void InitializePrePulseT0(float initPreT0) { GetWrapperTF1()->SetParameter(3, initPreT0);}
+
+ virtual float GetAmplitude() const {return _fitFunc->GetParameter(0); }
+ virtual float GetAmpError() const {return _fitFunc->GetParError(0); }
+
+ virtual float GetTau1() const {return _tau1;}
+ virtual float GetTau2() const {return _tau2;}
+
+ virtual float GetTime() const {
+ float fitT0 = _fitFunc->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 < _numPulses + 2) return _fitFunc->GetParameter(index);
+ else throw;
+ }
+
+ virtual float GetBkgdMaxFraction() const
+ {
+ return 0;
+ }
+
+ virtual double operator() (double *x, double *p)
+ {
+ double t = x[0];
+
+ double mainAmp = p[0];
+ double t0 = p[1];
+
+ double mainDeltaT = t - t0;
+ double funcValue = mainAmp*_norm*_expFermiFunc->operator()(mainDeltaT);
+
+ for (size_t ipulse = 0; ipulse < _numPulses - 1; ipulse++) {
+ double deltaT = mainDeltaT - _pulseDeltaT[ipulse];
+ double amp = p[2 + ipulse];
+
+ double pulse = amp*_norm*_expFermiFunc->operator()(deltaT);
+ funcValue += pulse;
+ }
+
+ return funcValue;
+ }
+};
#endif
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
index 8939f8170715..7da2b6605167 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZDCPulseAnalyzer.h
@@ -11,15 +11,19 @@
#include <TF1.h>
#include <TH1.h>
+#include <TGraph.h>
+#include <TList.h>
#include "ZdcAnalysis/ZDCFitWrapper.h"
+class TFitter;
+
class ZDCPulseAnalyzer
{
public:
- enum {PulseBit = 0, LowGainBit, FailBit, HGOverflowBit,
- HGUnderflowBit, PSHGOverUnderflowBit, LGOverflowBit, LGUnderflowBit,
- PrePulseBit, PostPulseBit, FitFailedBit, BadChisqBit, BadT0Bit};
+ enum {PulseBit = 0, LowGainBit = 1, FailBit = 2, HGOverflowBit = 3,
+ HGUnderflowBit = 4, PSHGOverUnderflowBit = 5, LGOverflowBit = 6, LGUnderflowBit = 7,
+ PrePulseBit = 8, PostPulseBit = 9, FitFailedBit = 10, BadChisqBit = 11, BadT0Bit = 12};
private:
typedef std::vector<float>::const_iterator SampleCIter;
@@ -28,6 +32,9 @@ private:
//
static std::string _fitOptions;
static bool _quietFits;
+ static TH1* _undelayedFitHist;
+ static TH1* _delayedFitHist;
+ static TF1* _combinedFitFunc;
// Quantities provided/set in the constructor
//
@@ -40,8 +47,10 @@ private:
bool _forceLG;
float _tmin;
float _tmax;
+
std::string _fitFunction;
- float _peak2ndDerivMinSample;
+ size_t _peak2ndDerivMinSample;
+ size_t _peak2ndDerivMinTolerance;
float _peak2ndDerivMinThreshLG;
float _peak2ndDerivMinThreshHG;
@@ -84,10 +93,15 @@ private:
ZDCFitWrapper* _defaultFitWrapper;
ZDCFitWrapper* _prePulseFitWrapper;
- // TH1* _LGFitHist;
- // TH1* _fitHist;
- // TF1* _FitFuncNoPrePulse;
- // TF1* _FitFuncPrePulse;
+ // Delayed pulse emembers
+ //
+ bool _useDelayed;
+ float _delayedDeltaT;
+ float _delayedPedestalDiff;
+ mutable TH1* _delayedHist;
+
+ TFitter* _prePulseCombinedFitter;
+ TFitter* _defaultCombinedFitter;
// Dynamic data loaded for each pulse (event)
// ==========================================
@@ -128,8 +142,9 @@ private:
float _fitAmplitude;
float _fitAmpError;
- float _fitT0;
- float _fitT0Corr;
+ float _fitTime;
+ float _fitTimeSub;
+ float _fitTimeCorr;
float _fitTau1;
float _fitTau2;
float _fitChisq;
@@ -137,6 +152,7 @@ private:
float _ampError;
float _preSampleAmp;
float _bkgdMaxFraction;
+ float _delayedBaselineShift;
std::vector<float> _ADCSamplesHGSub;
std::vector<float> _ADCSamplesLGSub;
@@ -154,7 +170,8 @@ private:
void SetDefaults();
void SetupFitFunctions();
- bool AnalyzeData(const std::vector<float>& samples, // The samples used for this event
+ bool AnalyzeData(size_t nSamples, size_t preSample,
+ 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
@@ -165,14 +182,36 @@ private:
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]);
+ if (!_useDelayed) {
+ // 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]);
+ }
+ }
+ else {
+ // Set the data and errors in the histogram object
+ //
+ for (size_t isample = 0; isample < _Nsample; isample++) {
+ _fitHist->SetBinContent(isample + 1, samples[isample*2]);
+ _delayedHist->SetBinContent(isample + 1, samples[isample*2 + 1]);
+
+ _fitHist->SetBinError(isample + 1, samplesSig[isample]); // ***** horrible hack: fix ASAP
+ }
+
}
}
+ void DoFit();
+ void DoFitCombined();
+
+ static TFitter* MakeCombinedFitter(TF1* func);
+
+ // The minuit FCN used for fitting combined undelayed and delayed pulses
+ //
+ static void CombinedPulsesFCN(int& numParam, double*, double& f, double* par, int flag);
+
public:
ZDCPulseAnalyzer(std::string tag, int Nsample, float deltaTSample, size_t preSampleIdx, int pedestal, float gainHG,
@@ -184,6 +223,10 @@ public:
static void SetQuietFits(bool quiet) {_quietFits = quiet;}
static bool QuietFits() {return _quietFits;}
+ void EnableDelayed(float deltaT, float pedestalShift);
+
+ void SetPeak2ndDerivMinTolerance(size_t tolerance) {_peak2ndDerivMinTolerance = tolerance;}
+
void SetForceLG(bool forceLG) {_forceLG = forceLG;}
bool ForceLG() const {return _forceLG;}
@@ -219,6 +262,9 @@ public:
bool LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG);
+ bool LoadAndAnalyzeData(std::vector<float> ADCSamplesHG, std::vector<float> ADCSamplesLG,
+ std::vector<float> ADCSamplesHGDelayed, std::vector<float> ADCSamplesLGDelayed);
+
bool HaveData() const {return _haveData;}
bool HavePulse() const {return _havePulse;}
bool Failed() const {return _fail;}
@@ -236,8 +282,9 @@ public:
bool PSHGOverUnderflow() const {return _PSHGOverUnderflow;}
float GetFitAmplitude() const {return _fitAmplitude;}
- float GetFitT0() const {return _fitT0;}
- float GetT0Corr() const {return _fitT0Corr;}
+ float GetFitT0() const {return _fitTime;}
+ float GetT0Sub() const {return _fitTimeSub;}
+ float GetT0Corr() const {return _fitTimeCorr;}
float GetChisq() const {return _fitChisq;}
float GetFitTau1() const {return _fitTau1;}
@@ -277,7 +324,32 @@ public:
return _fitHist;
}
- void DoFit();
+ TGraph* GetCombinedGraph() const {
+ //
+ // We defer filling the histogram if we don't have a pulse until the histogram is requested
+ //
+ GetHistogramPtr();
+
+ TGraph* theGraph = new TGraph(2*_Nsample);
+ size_t npts = 0;
+
+ for (size_t ipt = 0; ipt < _fitHist->GetNbinsX(); ipt++) {
+ theGraph->SetPoint(npts++, _fitHist->GetBinCenter(ipt + 1), _fitHist->GetBinContent(ipt + 1));
+ }
+
+ for (size_t iDelayPt = 0; iDelayPt < _delayedHist->GetNbinsX(); iDelayPt++) {
+ theGraph->SetPoint(npts++, _delayedHist->GetBinCenter(iDelayPt + 1), _delayedHist->GetBinContent(iDelayPt + 1) - _delayedBaselineShift);
+ }
+
+ TF1* func_p = (TF1*) _fitHist->GetListOfFunctions()->Last();
+ theGraph->GetListOfFunctions()->Add(func_p);
+ theGraph->SetName(( std::string(_fitHist->GetName()) + "combinaed").c_str());
+
+ theGraph->SetMarkerStyle(20);
+ theGraph->SetMarkerColor(1);
+
+ return theGraph;
+ }
void Dump() const;
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
index 278e115d3683..ee752e0ff923 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/ZdcAnalysis/ZdcAnalysisTool.h
@@ -71,6 +71,7 @@ class ZdcAnalysisTool : public virtual IZdcAnalysisTool, public asg::AsgTool
// Private methods
//
ZDCDataAnalyzer* initializeDefault();
+ ZDCDataAnalyzer* initializepPb2016();
StatusCode configureNewRun(unsigned int runNumber);
@@ -98,6 +99,7 @@ class ZdcAnalysisTool : public virtual IZdcAnalysisTool, public asg::AsgTool
const xAOD::ZdcModuleContainer* m_zdcModules;
bool m_flipEMDelay;
bool m_lowGainOnly;
+ bool m_combineDelay;
bool m_doCalib;
int m_forceCalibRun;
int m_forceCalibLB;
@@ -110,6 +112,7 @@ class ZdcAnalysisTool : public virtual IZdcAnalysisTool, public asg::AsgTool
unsigned int m_peakSample;
float m_Peak2ndDerivThresh;
float m_t0;
+ float m_delayDeltaT;
float m_tau1;
float m_tau2;
bool m_fixTau1;
diff --git a/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
index 4b26464a70b3..9f3ee50d3b7b 100644
--- a/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
+++ b/ForwardDetectors/ZDC/ZdcAnalysis/cmt/Makefile.RootCore
@@ -10,7 +10,7 @@
PACKAGE = ZdcAnalysis
# the libraries to link with this one:
-PACKAGE_PRELOAD =
+PACKAGE_PRELOAD = Minuit
# additional compilation flags to pass (not propagated to dependent packages):
PACKAGE_CXXFLAGS =
--
GitLab