#ifndef __JCALIBRATE_JHVINTERPOLATOR__ #define __JCALIBRATE_JHVINTERPOLATOR__ #include "JLang/JException.hh" #include "JTools/JRange.hh" #include "JCalibrate/JFitToT.hh" #include "TKey.h" #include "TGraph.h" #include "TGraphErrors.h" #include "TMultiGraph.h" #include "TAttMarker.h" #include "TColor.h" /** * \author bjung */ namespace JCALIBRATE { using JTOOLS::JRange; using JLANG::JNoValue; using JLANG::JValueOutOfRange; using JLANG::JNullPointerException; static const std::string HVINTERPOLATOR_SUFFIX = ".HVxG"; static const char* HVINTERPOLATOR_DATA = "HVINTERPOLATOR_DATA"; static const char* HVINTERPOLATOR_RESULT = "HVINTERPOLATOR_RESULT"; /** * Auxiliary data structure to store high-voltage versus gain data and interpolate the nominal high-voltage. */ struct JHVInterpolator { /** * Constructor. * * \param object TMultiGraph object */ JHVInterpolator(TMultiGraph& object) : data(&object) { TGraphErrors* g0 = new TGraphErrors(); g0->SetName(HVINTERPOLATOR_RESULT); g0->SetLineColor (kRed); g0->SetMarkerColor(kRed); g0->SetMarkerStyle(kFullDotSmall); g0->Set(0); data->Add(g0); TGraphErrors* g1 = new TGraphErrors(); g1->SetName(HVINTERPOLATOR_DATA); g1->SetMarkerStyle(kFullDotSmall); g1->Set(0); data->Add(g1); } /** * Add point to diagram. * * \param HV high-voltage [V] * \param gain gain value * \param gainError error on gain value */ void AddPoint(Double_t HV, Double_t gain, Double_t gainError) { TGraphErrors* g1 = getData(); const Int_t N = g1->GetN(); g1->SetPoint (N, fabs(HV), gain); g1->SetPointError(N, 0.0, gainError); } /** * Set point with index i. * * \param i index of point * \param HV high-voltage [V] * \param gain gain value * \param gainError error on gain value */ void SetPoint(Int_t i, Double_t HV, Double_t gain, Double_t gainError) { TGraphErrors* g1 = getData(); g1->SetPoint (i, fabs(HV), gain); g1->SetPointError(i, 0.0, gainError); } /** * Checks whether high-voltage is within range * * \param HV high-voltage [V] * \return true if high-voltage is within range; else false */ const bool checkHV(const double HV) const { return (HV > hvRange.getLowerLimit() && HV < hvRange.getUpperLimit()); } /** * Checks whether two high-voltage values are different. * * \param HV1 first high-voltage [V] * \param HV2 second high-voltage [V] * \return true if absolute difference is greater than minimal required difference; else false */ const bool checkHV(const double HV1, const double HV2) const { return (fabs(HV1 - HV2) > dHVmin); } /** * Checks if gain is within range. * * \param gain gain value * \return true if gain within allowed range; else false */ const bool checkGain(const double gain) const { return (gain > gainRange.getLowerLimit() && gain < gainRange.getUpperLimit()); } /** * Checks whether the gains of two points are strictly increasing as function of their absolute high-voltage * * \param i index of first point * \param j index of second point * \return true if gains of given points are strictly increasing\n * as function of the absolute high-voltage; else false */ const bool areIncreasing(const Int_t i, const Int_t j) const { TGraphErrors* g1 = getData(); if ((i >= 0 && i < g1->GetN()) && (j >= 0 && i < g1->GetN())) { return ((g1->GetPointX(i) > g1->GetPointX(j) && g1->GetPointY(i) > g1->GetPointY(j)) || (g1->GetPointX(i) < g1->GetPointX(j) && g1->GetPointY(i) < g1->GetPointY(j))); } else { return false; } } /** * Checks whether two points are valid for inter-/extrapolation. * * \param i index of first point * \param j index of second point * \return true if valid for inter-/extrapolation; else false */ const bool areValid(const Int_t i, const Int_t j) const { TGraphErrors* g1 = getData(); if ((i >= 0 && i < g1->GetN()) && (j >= 0 && j < g1->GetN())) { const double HV_i = -fabs(g1->GetPointX(i)); const double HV_j = -fabs(g1->GetPointX(j)); const double gain_i = g1->GetPointY(i); const double gain_j = g1->GetPointY(j); return (areIncreasing( i, j) && checkHV(HV_i) && checkGain(gain_i) && checkHV (HV_i, HV_j) && checkHV(HV_j) && checkGain(gain_j)); } else { return false; } } /** * Inter-/Extrapolate the high-voltage value corresponding to the target gain value. * * \param gainTarget target gain value for inter-/extrapolation * \return true if inter-/extrapolation successful; else false */ bool interpolateHV(const double gainTarget) { TGraphErrors* g1 = getData(); if (g1->GetN() < 2 || !checkGain(gainTarget)) { return false; } // Search for valid inter-/extrapolation points Int_t i = 0; Int_t j = 1; for (Int_t k = 2; k < g1->GetN(); ++k) { const double dGain_i = fabs(g1->GetPointY(i) - gainTarget); const double dGain_j = fabs(g1->GetPointY(j) - gainTarget); const double dGain_k = fabs(g1->GetPointY(k) - gainTarget); if (dGain_k < dGain_j) { i = (areValid(j,k) ? j : i); j = k; } else if ((dGain_k < dGain_i || !areValid(i,j)) && areValid(j,k)) { i = j; j = k; } } // Inter-/Extrapolate high-voltage corresponding to given gain if (areValid(i, j)) { const double logHV0 = log(fabs(g1->GetPointX(i))); const double logHV1 = log(fabs(g1->GetPointX(j))); const double logG0 = log(g1->GetPointY(i)); const double logG1 = log(g1->GetPointY(j)); const double elogG0 = g1->GetErrorY(i) / g1->GetPointY(i); const double elogG1 = g1->GetErrorY(j) / g1->GetPointY(j); const double dlogG0 = log(gainTarget) - logG0; const double dlogG1 = log(gainTarget) - logG1; const double slope = (logG1 - logG0) / (logHV1 - logHV0); const double HVnom = exp(dlogG0 / slope + logHV0); const double eHVnom = HVnom * sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 + dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0)); const double distance = fabs((log(HVnom) - logHV0) / (logHV0 - logHV1)); static const double maxDist = 2.0; if (checkHV(-HVnom) && distance < maxDist) { TGraphErrors* g0 = getResult(); g0->SetPoint (0, HVnom, gainTarget); g0->SetPointError(0, eHVnom, 0.0); return true; } } return false; } /** * Get interpolated high-voltage. * * \return inter-/extrapolated high-voltage\n * corresponding to target gain value [V] */ double getHV() const { TGraphErrors* g0 = getResult(); if (g0->GetN() > 0) { return g0->GetPointX(0); } else { THROW(JNoValue, "JHVInterpolator::getHV(): Missing HV inter-/extrapolation point. Please call JHVInterpolator::interpolateHV() first."); } } /** * Get error estimate on interpolated high-voltage. * * \return error on inter-/extrapolated high-voltage\n * corresponding to the target gain value [V] */ double getHVError() const { TGraphErrors* g0 = getResult(); if (g0->GetN() > 0) { return g0->GetErrorX(0); } else { THROW(JNoValue, "JHVInterpolator::getHVError(): Missing HV inter-/extrpolation point. Please call JHVInterpolator::interpolateHV() first."); } } /** * Get graph with the input data for the interpolation. * * \return pointer to graph with input data */ TGraphErrors* getData() const { TList* list = data->GetListOfGraphs(); if (list != NULL && list->FindObject(HVINTERPOLATOR_DATA) != NULL) { return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_DATA); } else if (list == NULL) { THROW(JNullPointerException, "JHVInterpolator()::getData(): Emtpy list of graphs."); } else { THROW(JNullPointerException, "JHVInterpolator()::getData(): No " << HVINTERPOLATOR_DATA << " graph."); } } /** * Get graph with the interpolation result. * * \return pointer to graph with the interpolation result */ TGraphErrors* getResult() const { TList* list = data->GetListOfGraphs(); if (list != NULL && list->FindObject(HVINTERPOLATOR_RESULT) != NULL) { return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_RESULT); } else if (list == NULL) { THROW(JNullPointerException, "JHVInterpolator::getResult(): Empty list of graphs."); } else { THROW(JNullPointerException, "JHVInterpolator::getResult(): No " << HVINTERPOLATOR_RESULT << " graph."); } } /** * Set minimal separation distance for high-voltage. * * \param minDist minimal separation distance for high-voltage */ static void setMinHVDistance(const double minDist) { dHVmin = minDist; } /** * Set valid gain range. * * \param range valid high-voltage range [V] */ static void setHVRange(const JRange range) { hvRange = range; } /** * Set valid gain range. * * \param range valid gain range */ static void setGainRange(const JRange range) { gainRange = range; } private: TMultiGraph* data; //!< HV-versus-gain data static double dHVmin; //!< Minimal high-voltage difference between two points [V] static JRange hvRange; //!< Allowed high-voltage range [V] static JRange gainRange; //!< Allowed gain range }; /** * Default values. */ double JHVInterpolator::dHVmin = 2 * 3.14; JRange JHVInterpolator::hvRange = JRange(-1500, -80); JRange JHVInterpolator::gainRange = JRange(FITTOT_GAIN_MIN, FITTOT_GAIN_MAX); } #endif