1 #ifndef __JCALIBRATE_JHVINTERPOLATOR__
2 #define __JCALIBRATE_JHVINTERPOLATOR__
12 #include "TGraphErrors.h"
13 #include "TMultiGraph.h"
14 #include "TAttMarker.h"
22 namespace JCALIBRATE {
58 const Int_t
N = g1->GetN();
60 g1->SetPoint (N, fabs(HV), gain);
61 g1->SetPointError(N, 0.0, gainError);
73 void SetPoint(Int_t i, Double_t HV, Double_t
gain, Double_t gainError)
77 g1->SetPoint (i, fabs(HV), gain);
78 g1->SetPointError(i, 0.0, gainError);
102 const bool checkHV(
const double HV1,
const double HV2)
const
104 return (fabs(HV1 - HV2) >
dHVmin);
133 if ((i >= 0 && i < g1->GetN()) &&
134 (j >= 0 && i < g1->GetN())) {
136 return ((g1->GetPointX(i) > g1->GetPointX(j) && g1->GetPointY(i) > g1->GetPointY(j)) ||
137 (g1->GetPointX(i) < g1->GetPointX(j) && g1->GetPointY(i) < g1->GetPointY(j)));
157 if ((i >= 0 && i < g1->GetN()) &&
158 (j >= 0 && i < g1->GetN())) {
160 const double HV_i = -fabs(g1->GetPointX(i));
161 const double HV_j = -fabs(g1->GetPointX(j));
162 const double gain_i = g1->GetPointY(i);
163 const double gain_j = g1->GetPointY(j);
185 if (g1->GetN() < 2 || !
checkGain(gainTarget)) {
194 for (Int_t
k = 1;
k < g1->GetN(); ++
k) {
196 const double dGain_i = fabs(g1->GetPointY(i) - gainTarget);
197 const double dGain_j = fabs(g1->GetPointY(j) - gainTarget);
198 const double dGain_k = fabs(g1->GetPointY(
k) - gainTarget);
200 if (dGain_k < dGain_i) {
205 }
else if ((dGain_k < dGain_j || !
areValid(i, j)) &&
215 const double logHV0 = log(fabs(g1->GetPointX(i)));
216 const double logHV1 = log(fabs(g1->GetPointX(j)));
218 const double logG0 = log(g1->GetPointY(i));
219 const double logG1 = log(g1->GetPointY(j));
220 const double elogG0 = g1->GetErrorY(i) / g1->GetPointY(i);
221 const double elogG1 = g1->GetErrorY(j) / g1->GetPointY(j);
223 const double dlogG0 = log(gainTarget) - logG0;
224 const double dlogG1 = log(gainTarget) - logG1;
226 const double slope = (logG1 - logG0) / (logHV1 - logHV0);
228 const double HVnom =
exp(dlogG0 / slope + logHV0);
229 const double eHVnom = HVnom * sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 +
230 dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0));
232 const double distance = fabs((log(HVnom) - logHV0) / (logHV0 - logHV1));
233 static const double maxDist = 2.0;
235 if (
checkHV(-HVnom) && distance < maxDist) {
239 g0->SetPoint (0, HVnom, gainTarget);
240 g0->SetPointError(0, eHVnom, 0.0);
259 static const double precision = 1e-3;
263 if ((g0->GetN() > 0 && fabs(g0->GetPointY(0) - gainTarget) < precision) ||
interpolate(gainTarget)) {
265 return g0->GetPointX(0);
269 THROW(
JValueOutOfRange,
"JHVInterpolator::getTargetError(): No valid inter-/extrapolation candidate found for given target gain of " << gainTarget <<
'.');
283 static const double precision = 1e-3;
287 if ((g0->GetN() > 0 && fabs(g0->GetPointY(0) - gainTarget) < precision) ||
interpolate(gainTarget)) {
289 return g0->GetErrorX(0);
293 THROW(
JValueOutOfRange,
"JHVInterpolator::getTargetError(): No valid inter-/extrapolation candidate found for given target gain of " << gainTarget <<
'.');
307 TList* list =
data->GetListOfGraphs();
315 g1 =
new TGraphErrors();
319 g1->SetMarkerStyle(kFullDotSmall);
339 TList* list =
data->GetListOfGraphs();
347 g0 =
new TGraphErrors();
351 g0->SetLineColor (kRed);
352 g0->SetMarkerColor(kRed);
353 g0->SetMarkerStyle(kFullDotSmall);
void SetPoint(Int_t i, Double_t HV, Double_t gain, Double_t gainError)
Set point with index i.
static const char * HVINTERPOLATOR_RESULT
TGraphErrors * getResult()
Get graph with the interpolation result.
double getTargetHVError(const double gainTarget)
Get error estimate on high-voltage corresponding to given target gain value.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
TGraphErrors * getData()
Get graph with the input data for the interpolation.
static const double FITTOT_GAIN_MAX
Default maximal gain.
static const double FITTOT_GAIN_MIN
Default minimal gain.
const bool areValid(const Int_t i, const Int_t j)
Checks whether two points are valid for inter-/extrapolation.
then fatal No sound hydrophone file $HYDROPHONE_TXT fi JGraph f $HYDROPHONE_TXT o $HYDROPHONE_ROOT sort gr k
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
static void setGainRange(const JRange< double > range)
Set valid gain range.
static const std::string HVINTERPOLATOR_SUFFIX
static JRange< double > hvRange
Allowed high-voltage range [V].
double getTargetHV(const double gainTarget)
Get high-voltage corresponding to given target gain value.
then set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))*exp(-0.5 *(y-[1])*(y-[1])/([2]*[2]))" JF2 -o $WORKDIR/f2.root -F "$FORMULA" -@ "p0
static void setHVRange(const JRange< double > range)
Set valid gain range.
static JRange< double > gainRange
Allowed gain range.
const bool checkHV(const double HV) const
Checks whether high-voltage is within range.
TMultiGraph * data
HV-versus-gain data.
static void setMinHVDistance(const double minDist)
Set minimal separation distance for high-voltage.
bool interpolate(const double gainTarget)
Inter-/Extrapolate the high-voltage value corresponding to the target gain value. ...
const bool checkGain(const double gain) const
Checks if gain is within range.
Auxiliary data structure to store high-voltage versus gain data and interpolate the nominal high-volt...
const bool areIncreasing(const Int_t i, const Int_t j)
Checks whether the gains of two points are strictly increasing as function of their absolute high-vol...
static const char * HVINTERPOLATOR_DATA
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
then $JPP_DIR examples JDetector JToT o $OUTPUT_FILE n N $NPE P gain
Auxiliary class to define a range between two values.
static double dHVmin
Minimal high-voltage difference between two points [V].
JHVInterpolator(TMultiGraph &object)
Constructor.
Exception for accessing a value in a collection that is outside of its range.
then usage $script[input file[working directory[option]]] nWhere option can be N
const bool checkHV(const double HV1, const double HV2) const
Checks whether two high-voltage values are different.
void AddPoint(Double_t HV, Double_t gain, Double_t gainError)
Add point to diagram.
Double_t g1(const Double_t x)
Function.