Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:30:52

0001 /****************************************************************************
0002 *
0003 * This is a part of TOTEM ofFileline software.
0004 *
0005 * Based on TMultiDimFet.cc (from ROOT) by Christian Holm Christense.
0006 *
0007 * Authors:
0008 *   Hubert Niewiadomski
0009 *   Jan Kašpar (jan.kaspar@gmail.com)
0010 *
0011 ****************************************************************************/
0012 
0013 #include "Riostream.h"
0014 #include "SimTransport/TotemRPProtonTransportParametrization/interface/TMultiDimFet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "TMath.h"
0017 #include "TH1.h"
0018 #include "TH2.h"
0019 #include "TROOT.h"
0020 #include "TDecompChol.h"
0021 #include "TDatime.h"
0022 #include <iostream>
0023 #include <map>
0024 
0025 #define RADDEG (180. / TMath::Pi())
0026 #define DEGRAD (TMath::Pi() / 180.)
0027 #define HIST_XORIG 0
0028 #define HIST_DORIG 1
0029 #define HIST_XNORM 2
0030 #define HIST_DSHIF 3
0031 #define HIST_RX 4
0032 #define HIST_RD 5
0033 #define HIST_RTRAI 6
0034 #define HIST_RTEST 7
0035 #define PARAM_MAXSTUDY 1
0036 #define PARAM_SEVERAL 2
0037 #define PARAM_RELERR 3
0038 #define PARAM_MAXTERMS 4
0039 
0040 //____________________________________________________________________
0041 //static void mdfHelper(int&, double*, double&, double*, int);
0042 
0043 //____________________________________________________________________
0044 ClassImp(TMultiDimFet);
0045 
0046 //____________________________________________________________________
0047 // Static instance. Used with mdfHelper and TMinuit
0048 //TMultiDimFet* TMultiDimFet::fgInstance = nullptr;
0049 
0050 //____________________________________________________________________
0051 TMultiDimFet::TMultiDimFet() {
0052   // Empty CTOR. Do not use
0053   fMeanQuantity = 0;
0054   fMaxQuantity = 0;
0055   fMinQuantity = 0;
0056   fSumSqQuantity = 0;
0057   fSumSqAvgQuantity = 0;
0058   fPowerLimit = 1;
0059 
0060   fMaxAngle = 0;
0061   fMinAngle = 1;
0062 
0063   fNVariables = 0;
0064   fMaxVariables = 0;
0065   fMinVariables = 0;
0066   fSampleSize = 0;
0067 
0068   fMaxAngle = 0;
0069   fMinAngle = 0;
0070 
0071   fPolyType = kMonomials;
0072   fShowCorrelation = kFALSE;
0073 
0074   fIsUserFunction = kFALSE;
0075 
0076   fHistograms = nullptr;
0077   fHistogramMask = 0;
0078 
0079   //fFitter                  = nullptr;
0080   //fgInstance               = nullptr;
0081 }
0082 
0083 const TMultiDimFet &TMultiDimFet::operator=(const TMultiDimFet &in) {
0084   if (this == &in) {
0085     return in;
0086   }
0087 
0088   fMeanQuantity = in.fMeanQuantity;  // Mean of dependent quantity
0089 
0090   fMaxQuantity = 0.0;       //! Max value of dependent quantity
0091   fMinQuantity = 0.0;       //! Min value of dependent quantity
0092   fSumSqQuantity = 0.0;     //! SumSquare of dependent quantity
0093   fSumSqAvgQuantity = 0.0;  //! Sum of squares away from mean
0094 
0095   fNVariables = in.fNVariables;  // Number of independent variables
0096 
0097   fMaxVariables.ResizeTo(in.fMaxVariables.GetLwb(), in.fMaxVariables.GetUpb());
0098   fMaxVariables = in.fMaxVariables;  // max value of independent variables
0099 
0100   fMinVariables.ResizeTo(in.fMinVariables.GetLwb(), in.fMinVariables.GetUpb());
0101   fMinVariables = in.fMinVariables;  // min value of independent variables
0102 
0103   fSampleSize = 0;      //! Size of training sample
0104   fTestSampleSize = 0;  //! Size of test sample
0105   fMinAngle = 1;        //! Min angle for acepting new function
0106   fMaxAngle = 0.0;      //! Max angle for acepting new function
0107 
0108   fMaxTerms = in.fMaxTerms;  // Max terms expected in final expr.
0109 
0110   fMinRelativeError = 0.0;  //! Min relative error accepted
0111 
0112   fMaxPowers.clear();  //! [fNVariables] maximum powers
0113   fPowerLimit = 1;     //! Control parameter
0114 
0115   fMaxFunctions = in.fMaxFunctions;  // max number of functions
0116 
0117   fFunctionCodes.clear();  //! [fMaxFunctions] acceptance code
0118   fMaxStudy = 0;           //! max functions to study
0119 
0120   fMaxPowersFinal.clear();  //! [fNVariables] maximum powers from fit;
0121 
0122   fMaxFunctionsTimesNVariables = in.fMaxFunctionsTimesNVariables;  // fMaxFunctionsTimesNVariables
0123   fPowers = in.fPowers;
0124 
0125   fPowerIndex = in.fPowerIndex;  // [fMaxTerms] Index of accepted powers
0126 
0127   fMaxResidual = 0.0;    //! Max redsidual value
0128   fMinResidual = 0.0;    //! Min redsidual value
0129   fMaxResidualRow = 0;   //! Row giving max residual
0130   fMinResidualRow = 0;   //! Row giving min residual
0131   fSumSqResidual = 0.0;  //! Sum of Square residuals
0132 
0133   fNCoefficients = in.fNCoefficients;  // Dimension of model coefficients
0134 
0135   fCoefficients.ResizeTo(in.fCoefficients.GetLwb(), in.fCoefficients.GetUpb());
0136   fCoefficients = in.fCoefficients;  // Vector of the final coefficients
0137 
0138   fRMS = 0.0;                   //! Root mean square of fit
0139   fChi2 = 0.0;                  //! Chi square of fit
0140   fParameterisationCode = 0;    //! Exit code of parameterisation
0141   fError = 0.0;                 //! Error from parameterization
0142   fTestError = 0.0;             //! Error from test
0143   fPrecision = 0.0;             //! Relative precision of param
0144   fTestPrecision = 0.0;         //! Relative precision of test
0145   fCorrelationCoeff = 0.0;      //! Multi Correlation coefficient
0146   fTestCorrelationCoeff = 0.0;  //! Multi Correlation coefficient
0147   fHistograms = nullptr;        //! List of histograms
0148   fHistogramMask = 0;           //! Bit pattern of hisograms used
0149   //fFitter = nullptr;            //! Fit object (MINUIT)
0150 
0151   fPolyType = in.fPolyType;                // Type of polynomials to use
0152   fShowCorrelation = in.fShowCorrelation;  // print correlation matrix
0153   fIsUserFunction = in.fIsUserFunction;    // Flag for user defined function
0154   fIsVerbose = in.fIsVerbose;              //
0155   return in;
0156 }
0157 
0158 //____________________________________________________________________
0159 TMultiDimFet::TMultiDimFet(Int_t dimension, EMDFPolyType type, Option_t *option)
0160     : TNamed("multidimfit", "Multi-dimensional fit object"),
0161       fQuantity(dimension),
0162       fSqError(dimension),
0163       fVariables(dimension * 100),
0164       fMeanVariables(dimension),
0165       fMaxVariables(dimension),
0166       fMinVariables(dimension) {
0167   // Constructor
0168   // Second argument is the type of polynomials to use in
0169   // parameterisation, one of:
0170   //      TMultiDimFet::kMonomials
0171   //      TMultiDimFet::kChebyshev
0172   //      TMultiDimFet::kLegendre
0173   //
0174   // Options:
0175   //   K      Compute (k)correlation matrix
0176   //   V      Be verbose
0177   //
0178   // Default is no options.
0179   //
0180 
0181   //fgInstance = this;
0182 
0183   fMeanQuantity = 0;
0184   fMaxQuantity = 0;
0185   fMinQuantity = 0;
0186   fSumSqQuantity = 0;
0187   fSumSqAvgQuantity = 0;
0188   fPowerLimit = 1;
0189 
0190   fMaxAngle = 0;
0191   fMinAngle = 1;
0192 
0193   fNVariables = dimension;
0194   fMaxVariables = 0;
0195   fMaxFunctionsTimesNVariables = fMaxFunctions * fNVariables;
0196   fMinVariables = 0;
0197   fSampleSize = 0;
0198   fTestSampleSize = 0;
0199   fMinRelativeError = 0.01;
0200   fError = 0;
0201   fTestError = 0;
0202   fPrecision = 0;
0203   fTestPrecision = 0;
0204   fParameterisationCode = 0;
0205 
0206   fPolyType = type;
0207   fShowCorrelation = kFALSE;
0208   fIsVerbose = kFALSE;
0209 
0210   TString opt = option;
0211   opt.ToLower();
0212 
0213   if (opt.Contains("k"))
0214     fShowCorrelation = kTRUE;
0215   if (opt.Contains("v"))
0216     fIsVerbose = kTRUE;
0217 
0218   fIsUserFunction = kFALSE;
0219 
0220   fHistograms = nullptr;
0221   fHistogramMask = 0;
0222 
0223   fMaxPowers.resize(dimension);
0224   fMaxPowersFinal.resize(dimension);
0225   //fFitter                 = nullptr;
0226 }
0227 
0228 //____________________________________________________________________
0229 TMultiDimFet::~TMultiDimFet() {
0230   if (fHistograms)
0231     fHistograms->Clear("nodelete");
0232   delete fHistograms;
0233 }
0234 
0235 //____________________________________________________________________
0236 void TMultiDimFet::AddRow(const Double_t *x, Double_t D, Double_t E) {
0237   // Add a row consisting of fNVariables independent variables, the
0238   // known, dependent quantity, and optionally, the square error in
0239   // the dependent quantity, to the training sample to be used for the
0240   // parameterization.
0241   // The mean of the variables and quantity is calculated on the fly,
0242   // as outlined in TPrincipal::AddRow.
0243   // This sample should be representive of the problem at hand.
0244   // Please note, that if no error is given Poisson statistics is
0245   // assumed and the square error is set to the value of dependent
0246   // quantity.  See also the
0247   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0248   if (!x)
0249     return;
0250 
0251   if (++fSampleSize == 1) {
0252     fMeanQuantity = D;
0253     fMaxQuantity = D;
0254     fMinQuantity = D;
0255   } else {
0256     fMeanQuantity *= 1 - 1. / Double_t(fSampleSize);
0257     fMeanQuantity += D / Double_t(fSampleSize);
0258     fSumSqQuantity += D * D;
0259 
0260     if (D >= fMaxQuantity)
0261       fMaxQuantity = D;
0262     if (D <= fMinQuantity)
0263       fMinQuantity = D;
0264   }
0265 
0266   // If the vector isn't big enough to hold the new data, then
0267   // expand the vector by half it's size.
0268   Int_t size = fQuantity.GetNrows();
0269   if (fSampleSize > size) {
0270     fQuantity.ResizeTo(size + size / 2);
0271     fSqError.ResizeTo(size + size / 2);
0272   }
0273 
0274   // Store the value
0275   fQuantity(fSampleSize - 1) = D;
0276   fSqError(fSampleSize - 1) = (E == 0 ? D : E);
0277 
0278   // Store data point in internal vector
0279   // If the vector isn't big enough to hold the new data, then
0280   // expand the vector by half it's size
0281   size = fVariables.GetNrows();
0282   if (fSampleSize * fNVariables > size)
0283     fVariables.ResizeTo(size + size / 2);
0284 
0285   // Increment the data point counter
0286   Int_t i, j;
0287   for (i = 0; i < fNVariables; i++) {
0288     if (fSampleSize == 1) {
0289       fMeanVariables(i) = x[i];
0290       fMaxVariables(i) = x[i];
0291       fMinVariables(i) = x[i];
0292     } else {
0293       fMeanVariables(i) *= 1 - 1. / Double_t(fSampleSize);
0294       fMeanVariables(i) += x[i] / Double_t(fSampleSize);
0295 
0296       // Update the maximum value for this component
0297       if (x[i] >= fMaxVariables(i))
0298         fMaxVariables(i) = x[i];
0299 
0300       // Update the minimum value for this component
0301       if (x[i] <= fMinVariables(i))
0302         fMinVariables(i) = x[i];
0303     }
0304 
0305     // Store the data.
0306     j = (fSampleSize - 1) * fNVariables + i;
0307     fVariables(j) = x[i];
0308   }
0309 }
0310 
0311 //____________________________________________________________________
0312 void TMultiDimFet::AddTestRow(const Double_t *x, Double_t D, Double_t E) {
0313   // Add a row consisting of fNVariables independent variables, the
0314   // known, dependent quantity, and optionally, the square error in
0315   // the dependent quantity, to the test sample to be used for the
0316   // test of the parameterization.
0317   // This sample needn't be representive of the problem at hand.
0318   // Please note, that if no error is given Poisson statistics is
0319   // assumed and the square error is set to the value of dependent
0320   // quantity.  See also the
0321   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0322   if (fTestSampleSize++ == 0) {
0323     fTestQuantity.ResizeTo(fNVariables);
0324     fTestSqError.ResizeTo(fNVariables);
0325     fTestVariables.ResizeTo(fNVariables * 100);
0326   }
0327 
0328   // If the vector isn't big enough to hold the new data, then
0329   // expand the vector by half it's size.
0330   Int_t size = fTestQuantity.GetNrows();
0331   if (fTestSampleSize > size) {
0332     fTestQuantity.ResizeTo(size + size / 2);
0333     fTestSqError.ResizeTo(size + size / 2);
0334   }
0335 
0336   // Store the value
0337   fTestQuantity(fTestSampleSize - 1) = D;
0338   fTestSqError(fTestSampleSize - 1) = (E == 0 ? D : E);
0339 
0340   // Store data point in internal vector
0341   // If the vector isn't big enough to hold the new data, then
0342   // expand the vector by half it's size
0343   size = fTestVariables.GetNrows();
0344   if (fTestSampleSize * fNVariables > size)
0345     fTestVariables.ResizeTo(size + size / 2);
0346 
0347   // Increment the data point counter
0348   Int_t i, j;
0349   for (i = 0; i < fNVariables; i++) {
0350     j = fNVariables * (fTestSampleSize - 1) + i;
0351     fTestVariables(j) = x[i];
0352 
0353     if (x[i] > fMaxVariables(i))
0354       Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f", i, fTestSampleSize, x[i], fMaxVariables(i));
0355     if (x[i] < fMinVariables(i))
0356       Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f", i, fTestSampleSize, x[i], fMinVariables(i));
0357   }
0358 }
0359 
0360 //____________________________________________________________________
0361 void TMultiDimFet::Clear(Option_t *option) {
0362   // Clear internal structures and variables
0363   Int_t i, j, n = fNVariables, m = fMaxFunctions;
0364 
0365   // Training sample, dependent quantity
0366   fQuantity.Zero();
0367   fSqError.Zero();
0368   fMeanQuantity = 0;
0369   fMaxQuantity = 0;
0370   fMinQuantity = 0;
0371   fSumSqQuantity = 0;
0372   fSumSqAvgQuantity = 0;
0373 
0374   // Training sample, independent variables
0375   fVariables.Zero();
0376   fNVariables = 0;
0377   fSampleSize = 0;
0378   fMeanVariables.Zero();
0379   fMaxVariables.Zero();
0380   fMinVariables.Zero();
0381 
0382   // Test sample
0383   fTestQuantity.Zero();
0384   fTestSqError.Zero();
0385   fTestVariables.Zero();
0386   fTestSampleSize = 0;
0387 
0388   // Functions
0389   fFunctions.Zero();
0390   fMaxFunctions = 0;
0391   fMaxStudy = 0;
0392   fMaxFunctionsTimesNVariables = 0;
0393   fOrthFunctions.Zero();
0394   fOrthFunctionNorms.Zero();
0395 
0396   // Control parameters
0397   fMinRelativeError = 0;
0398   fMinAngle = 0;
0399   fMaxAngle = 0;
0400   fMaxTerms = 0;
0401 
0402   // Powers
0403   for (i = 0; i < n; i++) {
0404     fMaxPowers[i] = 0;
0405     fMaxPowersFinal[i] = 0;
0406     for (j = 0; j < m; j++)
0407       fPowers[i * n + j] = 0;
0408   }
0409   fPowerLimit = 0;
0410 
0411   // Residuals
0412   fMaxResidual = 0;
0413   fMinResidual = 0;
0414   fMaxResidualRow = 0;
0415   fMinResidualRow = 0;
0416   fSumSqResidual = 0;
0417 
0418   // Fit
0419   fNCoefficients = 0;
0420   fOrthCoefficients = 0;
0421   fOrthCurvatureMatrix = 0;
0422   fRMS = 0;
0423   fCorrelationMatrix.Zero();
0424   fError = 0;
0425   fTestError = 0;
0426   fPrecision = 0;
0427   fTestPrecision = 0;
0428 
0429   // Coefficients
0430   fCoefficients.Zero();
0431   fCoefficientsRMS.Zero();
0432   fResiduals.Zero();
0433   fHistograms->Clear(option);
0434 
0435   // Options
0436   fPolyType = kMonomials;
0437   fShowCorrelation = kFALSE;
0438   fIsUserFunction = kFALSE;
0439 }
0440 
0441 //____________________________________________________________________
0442 Double_t TMultiDimFet::Eval(const Double_t *x, const Double_t *coeff) const {
0443   // Evaluate parameterization at point x. Optional argument coeff is
0444   // a vector of coefficients for the parameterisation, fNCoefficients
0445   // elements long.
0446   //   int fMaxFunctionsTimesNVariables = fMaxFunctions * fNVariables;
0447   Double_t returnValue = fMeanQuantity;
0448   Double_t term = 0;
0449   Int_t i, j;
0450 
0451   for (i = 0; i < fNCoefficients; i++) {
0452     // Evaluate the ith term in the expansion
0453     term = (coeff ? coeff[i] : fCoefficients(i));
0454     for (j = 0; j < fNVariables; j++) {
0455       // Evaluate the factor (polynomial) in the j-th variable.
0456       Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
0457       Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j)) * (x[j] - fMaxVariables(j));
0458       term *= EvalFactor(p, y);
0459     }
0460     // Add this term to the final result
0461     returnValue += term;
0462   }
0463   return returnValue;
0464 }
0465 
0466 void TMultiDimFet::ReducePolynomial(double error) {
0467   if (error == 0.0)
0468     return;
0469   else {
0470     ZeroDoubiousCoefficients(error);
0471   }
0472 }
0473 
0474 void TMultiDimFet::ZeroDoubiousCoefficients(double error) {
0475   typedef std::multimap<double, int> cmt;
0476   cmt m;
0477 
0478   for (int i = 0; i < fNCoefficients; i++) {
0479     m.insert(std::pair<double, int>(TMath::Abs(fCoefficients(i)), i));
0480   }
0481 
0482   double del_error_abs = 0;
0483   int deleted_terms_count = 0;
0484 
0485   for (cmt::iterator it = m.begin(); it != m.end() && del_error_abs < error; ++it) {
0486     if (TMath::Abs(it->first) + del_error_abs < error) {
0487       fCoefficients(it->second) = 0.0;
0488       del_error_abs = TMath::Abs(it->first) + del_error_abs;
0489       deleted_terms_count++;
0490     } else
0491       break;
0492   }
0493 
0494   int fNCoefficients_new = fNCoefficients - deleted_terms_count;
0495   TVectorD fCoefficients_new(fNCoefficients_new);
0496   std::vector<Int_t> fPowerIndex_new;
0497 
0498   int ind = 0;
0499   for (int i = 0; i < fNCoefficients; i++) {
0500     if (fCoefficients(i) != 0.0) {
0501       fCoefficients_new(ind) = fCoefficients(i);
0502       fPowerIndex_new.push_back(fPowerIndex[i]);
0503       ind++;
0504     }
0505   }
0506   fNCoefficients = fNCoefficients_new;
0507   fCoefficients.ResizeTo(fNCoefficients);
0508   fCoefficients = fCoefficients_new;
0509   fPowerIndex = fPowerIndex_new;
0510   edm::LogInfo("TMultiDimFet") << deleted_terms_count << " terms removed"
0511                                << "\n";
0512 }
0513 
0514 //____________________________________________________________________
0515 Double_t TMultiDimFet::EvalControl(const Int_t *iv) {
0516   // PRIVATE METHOD:
0517   // Calculate the control parameter from the passed powers
0518   Double_t s = 0;
0519   Double_t epsilon = 1e-6;  // a small number
0520   for (Int_t i = 0; i < fNVariables; i++) {
0521     if (fMaxPowers[i] != 1)
0522       s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1);
0523   }
0524   return s;
0525 }
0526 
0527 //____________________________________________________________________
0528 Double_t TMultiDimFet::EvalFactor(Int_t p, Double_t x) const {
0529   // PRIVATE METHOD:
0530   // Evaluate function with power p at variable value x
0531   Int_t i = 0;
0532   Double_t p1 = 1;
0533   Double_t p2 = 0;
0534   Double_t p3 = 0;
0535   Double_t r = 0;
0536 
0537   switch (p) {
0538     case 1:
0539       r = 1;
0540       break;
0541     case 2:
0542       r = x;
0543       break;
0544     default:
0545       p2 = x;
0546       for (i = 3; i <= p; i++) {
0547         p3 = p2 * x;
0548         if (fPolyType == kLegendre)
0549           p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1);
0550         else if (fPolyType == kChebyshev)
0551           p3 = 2 * x * p2 - p1;
0552         p1 = p2;
0553         p2 = p3;
0554       }
0555       r = p3;
0556   }
0557 
0558   return r;
0559 }
0560 
0561 //____________________________________________________________________
0562 void TMultiDimFet::FindParameterization(double precision) {
0563   // Find the parameterization
0564   //
0565   // Options:
0566   //     None so far
0567   //
0568   // For detailed description of what this entails, please refer to the
0569   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0570   MakeNormalized();
0571   MakeCandidates();
0572   MakeParameterization();
0573   MakeCoefficients();
0574   ReducePolynomial(precision);
0575 }
0576 
0577 //____________________________________________________________________
0578 /*
0579 void TMultiDimFet::Fit(Option_t *option)
0580 {
0581    // Try to fit the found parameterisation to the test sample.
0582    //
0583    // Options
0584    //     M     use Minuit to improve coefficients
0585    //
0586    // Also, refer to
0587    // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0588 
0589    //fgInstance = this;
0590 
0591    Int_t i, j;
0592    Double_t*      x    = new Double_t[fNVariables];
0593    Double_t  sumSqD    = 0;
0594    Double_t    sumD    = 0;
0595    Double_t  sumSqR    = 0;
0596    Double_t    sumR    = 0;
0597 
0598    // Calculate the residuals over the test sample
0599    for (i = 0; i < fTestSampleSize; i++) {
0600       for (j = 0; j < fNVariables; j++)
0601          x[j] = fTestVariables(i * fNVariables + j);
0602       Double_t res =  fTestQuantity(i) - Eval(x);
0603       sumD         += fTestQuantity(i);
0604       sumSqD       += fTestQuantity(i) * fTestQuantity(i);
0605       sumR         += res;
0606       sumSqR       += res * res;
0607       if (TESTBIT(fHistogramMask,HIST_RTEST))
0608          ((TH1D*)fHistograms->FindObject("res_test"))->Fill(res);
0609    }
0610    Double_t dAvg         = sumSqD - (sumD * sumD) / fTestSampleSize;
0611    Double_t rAvg         = sumSqR - (sumR * sumR) / fTestSampleSize;
0612    fTestCorrelationCoeff = (dAvg - rAvg) / dAvg;
0613    fTestError            = sumSqR;
0614    fTestPrecision        = sumSqR / sumSqD;
0615 
0616    TString opt(option);
0617    opt.ToLower();
0618 
0619    if (!opt.Contains("m"))
0620       MakeChi2();
0621 
0622    if (fNCoefficients * 50 > fTestSampleSize)
0623       Warning("Fit", "test sample is very small");
0624 
0625    if (!opt.Contains("m"))
0626       return;
0627 
0628    fFitter = TVirtualFitter::Fitter(nullptr,fNCoefficients);
0629    fFitter->SetFCN(mdfHelper);
0630 
0631    const Int_t  maxArgs = 16;
0632    Int_t           args = 1;
0633    Double_t*   arglist  = new Double_t[maxArgs];
0634    arglist[0]           = -1;
0635    fFitter->ExecuteCommand("SET PRINT",arglist,args);
0636 
0637    for (i = 0; i < fNCoefficients; i++) {
0638       Double_t startVal = fCoefficients(i);
0639       Double_t startErr = fCoefficientsRMS(i);
0640       fFitter->SetParameter(i, Form("coeff%02d",i),
0641          startVal, startErr, 0, 0);
0642    }
0643 
0644    args                 = 1;
0645    fFitter->ExecuteCommand("MIGRAD",arglist,args);
0646 
0647    for (i = 0; i < fNCoefficients; i++) {
0648       Double_t val = 0, err = 0, low = 0, high = 0;
0649       fFitter->GetParameter(i, Form("coeff%02d",i),
0650          val, err, low, high);
0651       fCoefficients(i)    = val;
0652       fCoefficientsRMS(i) = err;
0653    }
0654 }
0655 */
0656 
0657 //____________________________________________________________________
0658 void TMultiDimFet::MakeCandidates() {
0659   // PRIVATE METHOD:
0660   // Create list of candidate functions for the parameterisation. See
0661   // also
0662   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0663   Int_t i = 0;
0664   Int_t j = 0;
0665   Int_t k = 0;
0666 
0667   // The temporary array to store the powers in. We don't need to
0668   // initialize this array however.
0669   Int_t *powers = new Int_t[fNVariables * fMaxFunctions];
0670 
0671   // store of `control variables'
0672   Double_t *control = new Double_t[fMaxFunctions];
0673 
0674   // We've better initialize the variables
0675   Int_t *iv = new Int_t[fNVariables];
0676   for (i = 0; i < fNVariables; i++)
0677     iv[i] = 1;
0678 
0679   if (!fIsUserFunction) {
0680     // Number of funcs selected
0681     Int_t numberFunctions = 0;
0682 
0683     // Absolute max number of functions
0684     Int_t maxNumberFunctions = 1;
0685     for (i = 0; i < fNVariables; i++)
0686       maxNumberFunctions *= fMaxPowers[i];
0687 
0688     while (kTRUE) {
0689       // Get the control value for this function
0690       Double_t s = EvalControl(iv);
0691 
0692       if (s <= fPowerLimit) {
0693         // Call over-loadable method Select, as to allow the user to
0694         // interfere with the selection of functions.
0695         if (Select(iv)) {
0696           numberFunctions++;
0697 
0698           // If we've reached the user defined limit of how many
0699           // functions we can consider, break out of the loop
0700           if (numberFunctions > fMaxFunctions)
0701             break;
0702 
0703           // Store the control value, so we can sort array of powers
0704           // later on
0705           control[numberFunctions - 1] = Int_t(1.0e+6 * s);
0706 
0707           // Store the powers in powers array.
0708           for (i = 0; i < fNVariables; i++) {
0709             j = (numberFunctions - 1) * fNVariables + i;
0710             powers[j] = iv[i];
0711           }
0712         }  // if (Select())
0713       }    // if (s <= fPowerLimit)
0714 
0715       for (i = 0; i < fNVariables; i++)
0716         if (iv[i] < fMaxPowers[i])
0717           break;
0718 
0719       // If all variables have reached their maximum power, then we
0720       // break out of the loop
0721       if (i == fNVariables) {
0722         fMaxFunctions = numberFunctions;
0723         fMaxFunctionsTimesNVariables = fMaxFunctions * fNVariables;
0724         break;
0725       }
0726 
0727       // Next power in variable i
0728       iv[i]++;
0729 
0730       for (j = 0; j < i; j++)
0731         iv[j] = 1;
0732     }  // while (kTRUE)
0733   } else {
0734     // In case the user gave an explicit function
0735     for (i = 0; i < fMaxFunctions; i++) {
0736       // Copy the powers to working arrays
0737       for (j = 0; j < fNVariables; j++) {
0738         powers[i * fNVariables + j] = fPowers[i * fNVariables + j];
0739         iv[j] = fPowers[i * fNVariables + j];
0740       }
0741 
0742       control[i] = Int_t(1.0e+6 * EvalControl(iv));
0743     }
0744   }
0745 
0746   // Now we need to sort the powers according to least `control
0747   // variable'
0748   Int_t *order = new Int_t[fMaxFunctions];
0749   for (i = 0; i < fMaxFunctions; i++)
0750     order[i] = i;
0751   fPowers.resize(fMaxFunctions * fNVariables);
0752 
0753   for (i = 0; i < fMaxFunctions; i++) {
0754     Double_t x = control[i];
0755     Int_t l = order[i];
0756     k = i;
0757 
0758     for (j = i; j < fMaxFunctions; j++) {
0759       if (control[j] <= x) {
0760         x = control[j];
0761         l = order[j];
0762         k = j;
0763       }
0764     }
0765 
0766     if (k != i) {
0767       control[k] = control[i];
0768       control[i] = x;
0769       order[k] = order[i];
0770       order[i] = l;
0771     }
0772   }
0773 
0774   for (i = 0; i < fMaxFunctions; i++)
0775     for (j = 0; j < fNVariables; j++)
0776       fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j];
0777 
0778   delete[] control;
0779   delete[] powers;
0780   delete[] order;
0781   delete[] iv;
0782 }
0783 
0784 Double_t TMultiDimFet::MakeChi2(const Double_t *coeff) {
0785   // Calculate Chi square over either the test sample. The optional
0786   // argument coeff is a vector of coefficients to use in the
0787   // evaluation of the parameterisation. If coeff == 0, then the found
0788   // coefficients is used.
0789   // Used my MINUIT for fit (see TMultDimFit::Fit)
0790   fChi2 = 0;
0791   Int_t i, j;
0792   Double_t *x = new Double_t[fNVariables];
0793   for (i = 0; i < fTestSampleSize; i++) {
0794     // Get the stored point
0795     for (j = 0; j < fNVariables; j++)
0796       x[j] = fTestVariables(i * fNVariables + j);
0797 
0798     // Evaluate function. Scale to shifted values
0799     Double_t f = Eval(x, coeff);
0800 
0801     // Calculate contribution to Chic square
0802     fChi2 += 1. / TMath::Max(fTestSqError(i), 1e-20) * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
0803   }
0804 
0805   // Clean up
0806   delete[] x;
0807 
0808   return fChi2;
0809 }
0810 
0811 //____________________________________________________________________
0812 void TMultiDimFet::MakeCode(const char *filename, Option_t *option) {
0813   // Generate the file <filename> with .C appended if argument doesn't
0814   // end in .cxx or .C. The contains the implementation of the
0815   // function:
0816   //
0817   //   Double_t <funcname>(Double_t *x)
0818   //
0819   // which does the same as TMultiDimFet::Eval. Please refer to this
0820   // method.
0821   //
0822   // Further, the static variables:
0823   //
0824   //     Int_t    gNVariables
0825   //     Int_t    gNCoefficients
0826   //     Double_t gDMean
0827   //     Double_t gXMean[]
0828   //     Double_t gXMin[]
0829   //     Double_t gXMax[]
0830   //     Double_t gCoefficient[]
0831   //     Int_t    gPower[]
0832   //
0833   // are initialized. The only ROOT header file needed is Rtypes.h
0834   //
0835   // See TMultiDimFet::MakeRealCode for a list of options
0836 
0837   TString outName(filename);
0838   if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
0839     outName += ".C";
0840 
0841   MakeRealCode(outName.Data(), "", option);
0842 }
0843 
0844 void TMultiDimFet::MakeCoefficientErrors() {
0845   // PRIVATE METHOD:
0846   // Compute the errors on the coefficients. For this to be done, the
0847   // curvature matrix of the non-orthogonal functions, is computed.
0848   Int_t i = 0;
0849   Int_t j = 0;
0850   Int_t k = 0;
0851   TVectorD iF(fSampleSize);
0852   TVectorD jF(fSampleSize);
0853   fCoefficientsRMS.ResizeTo(fNCoefficients);
0854 
0855   TMatrixDSym curvatureMatrix(fNCoefficients);
0856 
0857   // Build the curvature matrix
0858   for (i = 0; i < fNCoefficients; i++) {
0859     iF = TMatrixDRow(fFunctions, i);
0860     for (j = 0; j <= i; j++) {
0861       jF = TMatrixDRow(fFunctions, j);
0862       for (k = 0; k < fSampleSize; k++)
0863         curvatureMatrix(i, j) += 1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k);
0864       curvatureMatrix(j, i) = curvatureMatrix(i, j);
0865     }
0866   }
0867 
0868   // Calculate Chi Square
0869   fChi2 = 0;
0870   for (i = 0; i < fSampleSize; i++) {
0871     Double_t f = 0;
0872     for (j = 0; j < fNCoefficients; j++)
0873       f += fCoefficients(j) * fFunctions(j, i);
0874     fChi2 += 1. / TMath::Max(fSqError(i), 1e-20) * (fQuantity(i) - f) * (fQuantity(i) - f);
0875   }
0876 
0877   // Invert the curvature matrix
0878   const TVectorD diag = TMatrixDDiag_const(curvatureMatrix);
0879   curvatureMatrix.NormByDiag(diag);
0880 
0881   TDecompChol chol(curvatureMatrix);
0882   if (!chol.Decompose())
0883     Error("MakeCoefficientErrors", "curvature matrix is singular");
0884   chol.Invert(curvatureMatrix);
0885 
0886   curvatureMatrix.NormByDiag(diag);
0887 
0888   for (i = 0; i < fNCoefficients; i++)
0889     fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i, i));
0890 }
0891 
0892 //____________________________________________________________________
0893 void TMultiDimFet::MakeCoefficients() {
0894   // PRIVATE METHOD:
0895   // Invert the model matrix B, and compute final coefficients. For a
0896   // more thorough discussion of what this means, please refer to the
0897   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
0898   //
0899   // First we invert the lower triangle matrix fOrthCurvatureMatrix
0900   // and store the inverted matrix in the upper triangle.
0901 
0902   Int_t i = 0, j = 0;
0903   Int_t col = 0, row = 0;
0904 
0905   // Invert the B matrix
0906   for (col = 1; col < fNCoefficients; col++) {
0907     for (row = col - 1; row > -1; row--) {
0908       fOrthCurvatureMatrix(row, col) = 0;
0909       for (i = row; i <= col; i++)
0910         fOrthCurvatureMatrix(row, col) -= fOrthCurvatureMatrix(i, row) * fOrthCurvatureMatrix(i, col);
0911     }
0912   }
0913 
0914   // Compute the final coefficients
0915   fCoefficients.ResizeTo(fNCoefficients);
0916 
0917   for (i = 0; i < fNCoefficients; i++) {
0918     Double_t sum = 0;
0919     for (j = i; j < fNCoefficients; j++)
0920       sum += fOrthCurvatureMatrix(i, j) * fOrthCoefficients(j);
0921     fCoefficients(i) = sum;
0922   }
0923 
0924   // Compute the final residuals
0925   fResiduals.ResizeTo(fSampleSize);
0926   for (i = 0; i < fSampleSize; i++)
0927     fResiduals(i) = fQuantity(i);
0928 
0929   for (i = 0; i < fNCoefficients; i++)
0930     for (j = 0; j < fSampleSize; j++)
0931       fResiduals(j) -= fCoefficients(i) * fFunctions(i, j);
0932 
0933   // Compute the max and minimum, and squared sum of the evaluated
0934   // residuals
0935   fMinResidual = 10e10;
0936   fMaxResidual = -10e10;
0937   Double_t sqRes = 0;
0938   for (i = 0; i < fSampleSize; i++) {
0939     sqRes += fResiduals(i) * fResiduals(i);
0940     if (fResiduals(i) <= fMinResidual) {
0941       fMinResidual = fResiduals(i);
0942       fMinResidualRow = i;
0943     }
0944     if (fResiduals(i) >= fMaxResidual) {
0945       fMaxResidual = fResiduals(i);
0946       fMaxResidualRow = i;
0947     }
0948   }
0949 
0950   fCorrelationCoeff = fSumSqResidual / fSumSqAvgQuantity;
0951   fPrecision = TMath::Sqrt(sqRes / fSumSqQuantity);
0952 
0953   // If we use histograms, fill some more
0954   if (TESTBIT(fHistogramMask, HIST_RD) || TESTBIT(fHistogramMask, HIST_RTRAI) || TESTBIT(fHistogramMask, HIST_RX)) {
0955     for (i = 0; i < fSampleSize; i++) {
0956       if (TESTBIT(fHistogramMask, HIST_RD))
0957         ((TH2D *)fHistograms->FindObject("res_d"))->Fill(fQuantity(i), fResiduals(i));
0958       if (TESTBIT(fHistogramMask, HIST_RTRAI))
0959         ((TH1D *)fHistograms->FindObject("res_train"))->Fill(fResiduals(i));
0960 
0961       if (TESTBIT(fHistogramMask, HIST_RX))
0962         for (j = 0; j < fNVariables; j++)
0963           ((TH2D *)fHistograms->FindObject(Form("res_x_%d", j)))->Fill(fVariables(i * fNVariables + j), fResiduals(i));
0964     }
0965   }  // If histograms
0966 }
0967 
0968 //____________________________________________________________________
0969 void TMultiDimFet::MakeCorrelation() {
0970   // PRIVATE METHOD:
0971   // Compute the correlation matrix
0972   if (!fShowCorrelation)
0973     return;
0974 
0975   fCorrelationMatrix.ResizeTo(fNVariables, fNVariables + 1);
0976 
0977   Double_t d2 = 0;
0978   Double_t ddotXi = 0;   // G.Q. needs to be reinitialized in the loop over i fNVariables
0979   Double_t xiNorm = 0;   // G.Q. needs to be reinitialized in the loop over i fNVariables
0980   Double_t xidotXj = 0;  // G.Q. needs to be reinitialized in the loop over j fNVariables
0981   Double_t xjNorm = 0;   // G.Q. needs to be reinitialized in the loop over j fNVariables
0982 
0983   Int_t i, j, k, l, m;  // G.Q. added m variable
0984   for (i = 0; i < fSampleSize; i++)
0985     d2 += fQuantity(i) * fQuantity(i);
0986 
0987   for (i = 0; i < fNVariables; i++) {
0988     ddotXi = 0.;  // G.Q. reinitialisation
0989     xiNorm = 0.;  // G.Q. reinitialisation
0990     for (j = 0; j < fSampleSize; j++) {
0991       // Index of sample j of variable i
0992       k = j * fNVariables + i;
0993       ddotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i));
0994       xiNorm += (fVariables(k) - fMeanVariables(i)) * (fVariables(k) - fMeanVariables(i));
0995     }
0996     fCorrelationMatrix(i, 0) = ddotXi / TMath::Sqrt(d2 * xiNorm);
0997 
0998     for (j = 0; j < i; j++) {
0999       xidotXj = 0.;  // G.Q. reinitialisation
1000       xjNorm = 0.;   // G.Q. reinitialisation
1001       for (k = 0; k < fSampleSize; k++) {
1002         // Index of sample j of variable i
1003         // l =  j * fNVariables + k;  // G.Q.
1004         l = k * fNVariables + j;  // G.Q.
1005         m = k * fNVariables + i;  // G.Q.
1006         // G.Q.        xidotXj += (fVariables(i) - fMeanVariables(i))
1007         // G.Q.          * (fVariables(l) - fMeanVariables(j));
1008         xidotXj +=
1009             (fVariables(m) - fMeanVariables(i)) * (fVariables(l) - fMeanVariables(j));  // G.Q. modified index for Xi
1010         xjNorm += (fVariables(l) - fMeanVariables(j)) * (fVariables(l) - fMeanVariables(j));
1011       }
1012       fCorrelationMatrix(i, j + 1) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1013     }
1014   }
1015 }
1016 
1017 //____________________________________________________________________
1018 Double_t TMultiDimFet::MakeGramSchmidt(Int_t function) {
1019   // PRIVATE METHOD:
1020   // Make Gram-Schmidt orthogonalisation. The class description gives
1021   // a thorough account of this algorithm, as well as
1022   // references. Please refer to the
1023   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
1024 
1025   // calculate w_i, that is, evaluate the current function at data
1026   // point i
1027   Double_t f2 = 0;
1028   fOrthCoefficients(fNCoefficients) = 0;
1029   fOrthFunctionNorms(fNCoefficients) = 0;
1030   Int_t j = 0;
1031   Int_t k = 0;
1032 
1033   for (j = 0; j < fSampleSize; j++) {
1034     fFunctions(fNCoefficients, j) = 1;
1035     fOrthFunctions(fNCoefficients, j) = 0;
1036     // First, however, we need to calculate f_fNCoefficients
1037     for (k = 0; k < fNVariables; k++) {
1038       Int_t p = fPowers[function * fNVariables + k];
1039       Double_t x = fVariables(j * fNVariables + k);
1040       fFunctions(fNCoefficients, j) *= EvalFactor(p, x);
1041     }
1042 
1043     // Calculate f dot f in f2
1044     f2 += fFunctions(fNCoefficients, j) * fFunctions(fNCoefficients, j);
1045     // Assign to w_fNCoefficients f_fNCoefficients
1046     fOrthFunctions(fNCoefficients, j) = fFunctions(fNCoefficients, j);
1047   }
1048 
1049   // the first column of w is equal to f
1050   for (j = 0; j < fNCoefficients; j++) {
1051     Double_t fdw = 0;
1052     // Calculate (f_fNCoefficients dot w_j) / w_j^2
1053     for (k = 0; k < fSampleSize; k++) {
1054       fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j, k) / fOrthFunctionNorms(j);
1055     }
1056 
1057     fOrthCurvatureMatrix(fNCoefficients, j) = fdw;
1058     // and subtract it from the current value of w_ij
1059     for (k = 0; k < fSampleSize; k++)
1060       fOrthFunctions(fNCoefficients, k) -= fdw * fOrthFunctions(j, k);
1061   }
1062 
1063   for (j = 0; j < fSampleSize; j++) {
1064     // calculate squared length of w_fNCoefficients
1065     fOrthFunctionNorms(fNCoefficients) += fOrthFunctions(fNCoefficients, j) * fOrthFunctions(fNCoefficients, j);
1066 
1067     // calculate D dot w_fNCoefficients in A
1068     fOrthCoefficients(fNCoefficients) += fQuantity(j) * fOrthFunctions(fNCoefficients, j);
1069   }
1070 
1071   // First test, but only if didn't user specify
1072   if (!fIsUserFunction)
1073     if (TMath::Sqrt(fOrthFunctionNorms(fNCoefficients) / (f2 + 1e-10)) < TMath::Sin(fMinAngle * DEGRAD))
1074       return 0;
1075 
1076   // The result found by this code for the first residual is always
1077   // much less then the one found be MUDIFI. That's because it's
1078   // supposed to be. The cause is the improved precision of Double_t
1079   // over DOUBLE PRECISION!
1080   fOrthCurvatureMatrix(fNCoefficients, fNCoefficients) = 1;
1081   Double_t b = fOrthCoefficients(fNCoefficients);
1082   fOrthCoefficients(fNCoefficients) /= fOrthFunctionNorms(fNCoefficients);
1083 
1084   // Calculate the residual from including this fNCoefficients.
1085   Double_t dResidur = fOrthCoefficients(fNCoefficients) * b;
1086 
1087   return dResidur;
1088 }
1089 
1090 //____________________________________________________________________
1091 void TMultiDimFet::MakeHistograms(Option_t *option) {
1092   // Make histograms of the result of the analysis. This message
1093   // should be sent after having read all data points, but before
1094   // finding the parameterization
1095   //
1096   // Options:
1097   //     A         All the below
1098   //     X         Original independent variables
1099   //     D         Original dependent variables
1100   //     N         Normalised independent variables
1101   //     S         Shifted dependent variables
1102   //     R1        Residuals versus normalised independent variables
1103   //     R2        Residuals versus dependent variable
1104   //     R3        Residuals computed on training sample
1105   //     R4        Residuals computed on test sample
1106   //
1107   // For a description of these quantities, refer to
1108   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
1109   TString opt(option);
1110   opt.ToLower();
1111 
1112   if (opt.Length() < 1)
1113     return;
1114 
1115   if (!fHistograms)
1116     fHistograms = new TList;
1117 
1118   // Counter variable
1119   Int_t i = 0;
1120 
1121   // Histogram of original variables
1122   if (opt.Contains("x") || opt.Contains("a")) {
1123     SETBIT(fHistogramMask, HIST_XORIG);
1124     for (i = 0; i < fNVariables; i++)
1125       if (!fHistograms->FindObject(Form("x_%d_orig", i)))
1126         fHistograms->Add(
1127             new TH1D(Form("x_%d_orig", i), Form("Original variable # %d", i), 100, fMinVariables(i), fMaxVariables(i)));
1128   }
1129 
1130   // Histogram of original dependent variable
1131   if (opt.Contains("d") || opt.Contains("a")) {
1132     SETBIT(fHistogramMask, HIST_DORIG);
1133     if (!fHistograms->FindObject("d_orig"))
1134       fHistograms->Add(new TH1D("d_orig", "Original Quantity", 100, fMinQuantity, fMaxQuantity));
1135   }
1136 
1137   // Histograms of normalized variables
1138   if (opt.Contains("n") || opt.Contains("a")) {
1139     SETBIT(fHistogramMask, HIST_XNORM);
1140     for (i = 0; i < fNVariables; i++)
1141       if (!fHistograms->FindObject(Form("x_%d_norm", i)))
1142         fHistograms->Add(new TH1D(Form("x_%d_norm", i), Form("Normalized variable # %d", i), 100, -1, 1));
1143   }
1144 
1145   // Histogram of shifted dependent variable
1146   if (opt.Contains("s") || opt.Contains("a")) {
1147     SETBIT(fHistogramMask, HIST_DSHIF);
1148     if (!fHistograms->FindObject("d_shifted"))
1149       fHistograms->Add(
1150           new TH1D("d_shifted", "Shifted Quantity", 100, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity));
1151   }
1152 
1153   // Residual from training sample versus independent variables
1154   if (opt.Contains("r1") || opt.Contains("a")) {
1155     SETBIT(fHistogramMask, HIST_RX);
1156     for (i = 0; i < fNVariables; i++)
1157       if (!fHistograms->FindObject(Form("res_x_%d", i)))
1158         fHistograms->Add(new TH2D(Form("res_x_%d", i),
1159                                   Form("Computed residual versus x_%d", i),
1160                                   100,
1161                                   -1,
1162                                   1,
1163                                   35,
1164                                   fMinQuantity - fMeanQuantity,
1165                                   fMaxQuantity - fMeanQuantity));
1166   }
1167 
1168   // Residual from training sample versus. dependent variable
1169   if (opt.Contains("r2") || opt.Contains("a")) {
1170     SETBIT(fHistogramMask, HIST_RD);
1171     if (!fHistograms->FindObject("res_d"))
1172       fHistograms->Add(new TH2D("res_d",
1173                                 "Computed residuals vs Quantity",
1174                                 100,
1175                                 fMinQuantity - fMeanQuantity,
1176                                 fMaxQuantity - fMeanQuantity,
1177                                 35,
1178                                 fMinQuantity - fMeanQuantity,
1179                                 fMaxQuantity - fMeanQuantity));
1180   }
1181 
1182   // Residual from training sample
1183   if (opt.Contains("r3") || opt.Contains("a")) {
1184     SETBIT(fHistogramMask, HIST_RTRAI);
1185     if (!fHistograms->FindObject("res_train"))
1186       fHistograms->Add(new TH1D("res_train",
1187                                 "Computed residuals over training sample",
1188                                 100,
1189                                 fMinQuantity - fMeanQuantity,
1190                                 fMaxQuantity - fMeanQuantity));
1191   }
1192   if (opt.Contains("r4") || opt.Contains("a")) {
1193     SETBIT(fHistogramMask, HIST_RTEST);
1194     if (!fHistograms->FindObject("res_test"))
1195       fHistograms->Add(new TH1D("res_test",
1196                                 "Distribution of residuals from test",
1197                                 100,
1198                                 fMinQuantity - fMeanQuantity,
1199                                 fMaxQuantity - fMeanQuantity));
1200   }
1201 }
1202 
1203 //____________________________________________________________________
1204 void TMultiDimFet::MakeMethod(const Char_t *classname, Option_t *option) {
1205   // Generate the file <classname>MDF.cxx which contains the
1206   // implementation of the method:
1207   //
1208   //   Double_t <classname>::MDF(Double_t *x)
1209   //
1210   // which does the same as  TMultiDimFet::Eval. Please refer to this
1211   // method.
1212   //
1213   // Further, the public static members:
1214   //
1215   //   Int_t    <classname>::fgNVariables
1216   //   Int_t    <classname>::fgNCoefficients
1217   //   Double_t <classname>::fgDMean
1218   //   Double_t <classname>::fgXMean[]       //[fgNVariables]
1219   //   Double_t <classname>::fgXMin[]        //[fgNVariables]
1220   //   Double_t <classname>::fgXMax[]        //[fgNVariables]
1221   //   Double_t <classname>::fgCoefficient[] //[fgNCoeffficents]
1222   //   Int_t    <classname>::fgPower[]       //[fgNCoeffficents*fgNVariables]
1223   //
1224   // are initialized, and assumed to exist. The class declaration is
1225   // assumed to be in <classname>.h and assumed to be provided by the
1226   // user.
1227   //
1228   // See TMultiDimFet::MakeRealCode for a list of options
1229   //
1230   // The minimal class definition is:
1231   //
1232   //   class <classname> {
1233   //   public:
1234   //     Int_t    <classname>::fgNVariables;     // Number of variables
1235   //     Int_t    <classname>::fgNCoefficients;  // Number of terms
1236   //     Double_t <classname>::fgDMean;          // Mean from training sample
1237   //     Double_t <classname>::fgXMean[];        // Mean from training sample
1238   //     Double_t <classname>::fgXMin[];         // Min from training sample
1239   //     Double_t <classname>::fgXMax[];         // Max from training sample
1240   //     Double_t <classname>::fgCoefficient[];  // Coefficients
1241   //     Int_t    <classname>::fgPower[];        // Function powers
1242   //
1243   //     Double_t Eval(Double_t *x);
1244   //   };
1245   //
1246   // Whether the method <classname>::Eval should be static or not, is
1247   // up to the user.
1248 
1249   MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
1250 }
1251 
1252 //____________________________________________________________________
1253 void TMultiDimFet::MakeNormalized() {
1254   // PRIVATE METHOD:
1255   // Normalize data to the interval [-1;1]. This is needed for the
1256   // classes method to work.
1257 
1258   Int_t i = 0;
1259   Int_t j = 0;
1260   Int_t k = 0;
1261 
1262   for (i = 0; i < fSampleSize; i++) {
1263     if (TESTBIT(fHistogramMask, HIST_DORIG))
1264       ((TH1D *)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));
1265 
1266     fQuantity(i) -= fMeanQuantity;
1267     fSumSqAvgQuantity += fQuantity(i) * fQuantity(i);
1268 
1269     if (TESTBIT(fHistogramMask, HIST_DSHIF))
1270       ((TH1D *)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));
1271 
1272     for (j = 0; j < fNVariables; j++) {
1273       Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
1274       k = i * fNVariables + j;
1275 
1276       // Fill histograms of original independent variables
1277       if (TESTBIT(fHistogramMask, HIST_XORIG))
1278         ((TH1D *)fHistograms->FindObject(Form("x_%d_orig", j)))->Fill(fVariables(k));
1279 
1280       // Normalise independent variables
1281       fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1282 
1283       // Fill histograms of normalised independent variables
1284       if (TESTBIT(fHistogramMask, HIST_XNORM))
1285         ((TH1D *)fHistograms->FindObject(Form("x_%d_norm", j)))->Fill(fVariables(k));
1286     }
1287   }
1288   // Shift min and max of dependent variable
1289   fMaxQuantity -= fMeanQuantity;
1290   fMinQuantity -= fMeanQuantity;
1291 
1292   // Shift mean of independent variables
1293   for (i = 0; i < fNVariables; i++) {
1294     Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
1295     fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i) - fMaxVariables(i));
1296   }
1297 }
1298 
1299 //____________________________________________________________________
1300 void TMultiDimFet::MakeParameterization() {
1301   // PRIVATE METHOD:
1302   // Find the parameterization over the training sample. A full account
1303   // of the algorithm is given in the
1304   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
1305 
1306   Int_t i = -1;
1307   Int_t j = 0;
1308   Int_t k = 0;
1309   Int_t maxPass = 3;
1310   Int_t studied = 0;
1311   Double_t squareResidual = fSumSqAvgQuantity;
1312   fNCoefficients = 0;
1313   fSumSqResidual = fSumSqAvgQuantity;
1314   fFunctions.ResizeTo(fMaxTerms, fSampleSize);
1315   fOrthFunctions.ResizeTo(fMaxTerms, fSampleSize);
1316   fOrthFunctionNorms.ResizeTo(fMaxTerms);
1317   fOrthCoefficients.ResizeTo(fMaxTerms);
1318   fOrthCurvatureMatrix.ResizeTo(fMaxTerms, fMaxTerms);
1319   fFunctions = 1;
1320 
1321   fFunctionCodes.resize(fMaxFunctions);
1322   fPowerIndex.resize(fMaxTerms);
1323   Int_t l;
1324   for (l = 0; l < fMaxFunctions; l++)
1325     fFunctionCodes[l] = 0;
1326   for (l = 0; l < fMaxTerms; l++)
1327     fPowerIndex[l] = 0;
1328 
1329   if (fMaxAngle != 0)
1330     maxPass = 100;
1331   if (fIsUserFunction)
1332     maxPass = 1;
1333 
1334   // Loop over the number of functions we want to study.
1335   // increment inspection counter
1336   while (kTRUE) {
1337     // Reach user defined limit of studies
1338     if (studied++ >= fMaxStudy) {
1339       fParameterisationCode = PARAM_MAXSTUDY;
1340       break;
1341     }
1342 
1343     // Considered all functions several times
1344     if (k >= maxPass) {
1345       fParameterisationCode = PARAM_SEVERAL;
1346       break;
1347     }
1348 
1349     // increment function counter
1350     i++;
1351 
1352     // If we've reached the end of the functions, restart pass
1353     if (i == fMaxFunctions) {
1354       if (fMaxAngle != 0)
1355         fMaxAngle += (90 - fMaxAngle) / 2;
1356       i = 0;
1357       studied--;
1358       k++;
1359       continue;
1360     }
1361     if (studied == 1)
1362       fFunctionCodes[i] = 0;
1363     else if (fFunctionCodes[i] >= 2)
1364       continue;
1365 
1366     // Print a happy message
1367     if (fIsVerbose && studied == 1)
1368       edm::LogInfo("TMultiDimFet") << "Coeff   SumSqRes    Contrib   Angle      QM   Func"
1369                                    << "     Value        W^2  Powers"
1370                                    << "\n";
1371 
1372     // Make the Gram-Schmidt
1373     Double_t dResidur = MakeGramSchmidt(i);
1374 
1375     if (dResidur == 0) {
1376       // This function is no good!
1377       // First test is in MakeGramSchmidt
1378       fFunctionCodes[i] = 1;
1379       continue;
1380     }
1381 
1382     // If user specified function, assume she/he knows what he's doing
1383     if (!fIsUserFunction) {
1384       // Flag this function as considered
1385       fFunctionCodes[i] = 2;
1386 
1387       // Test if this function contributes to the fit
1388       if (!TestFunction(squareResidual, dResidur)) {
1389         fFunctionCodes[i] = 1;
1390         continue;
1391       }
1392     }
1393 
1394     // If we get to here, the function currently considered is
1395     // fNCoefficients, so we increment the counter
1396     // Flag this function as OK, and store and the number in the
1397     // index.
1398     fFunctionCodes[i] = 3;
1399     fPowerIndex[fNCoefficients] = i;
1400     fNCoefficients++;
1401 
1402     // We add the current contribution to the sum of square of
1403     // residuals;
1404     squareResidual -= dResidur;
1405 
1406     // Calculate control parameter from this function
1407     for (j = 0; j < fNVariables; j++) {
1408       if (fNCoefficients == 1 || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
1409         fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
1410     }
1411     Double_t s = EvalControl(&fPowers[i * fNVariables]);
1412 
1413     // Print the statistics about this function
1414     if (fIsVerbose) {
1415       edm::LogVerbatim("TMultiDimFet") << std::setw(5) << fNCoefficients << " " << std::setw(10) << std::setprecision(4)
1416                                        << squareResidual << " " << std::setw(10) << std::setprecision(4) << dResidur
1417                                        << " " << std::setw(7) << std::setprecision(3) << fMaxAngle << " "
1418                                        << std::setw(7) << std::setprecision(3) << s << " " << std::setw(5) << i << " "
1419                                        << std::setw(10) << std::setprecision(4) << fOrthCoefficients(fNCoefficients - 1)
1420                                        << " " << std::setw(10) << std::setprecision(4)
1421                                        << fOrthFunctionNorms(fNCoefficients - 1) << " " << std::flush;
1422       for (j = 0; j < fNVariables; j++)
1423         edm::LogInfo("TMultiDimFet") << " " << fPowers[i * fNVariables + j] - 1 << std::flush;
1424       edm::LogInfo("TMultiDimFet") << "\n";
1425     }
1426 
1427     if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
1428       fParameterisationCode = PARAM_MAXTERMS;
1429       break;
1430     }
1431 
1432     Double_t err = TMath::Sqrt(TMath::Max(1e-20, squareResidual) / fSumSqAvgQuantity);
1433     if (err < fMinRelativeError) {
1434       fParameterisationCode = PARAM_RELERR;
1435       break;
1436     }
1437   }
1438 
1439   fError = TMath::Max(1e-20, squareResidual);
1440   fSumSqResidual -= fError;
1441   fRMS = TMath::Sqrt(fError / fSampleSize);
1442 }
1443 
1444 //____________________________________________________________________
1445 void TMultiDimFet::MakeRealCode(const char *filename, const char *classname, Option_t *) {
1446   // PRIVATE METHOD:
1447   // This is the method that actually generates the code for the
1448   // evaluation the parameterization on some point.
1449   // It's called by TMultiDimFet::MakeCode and TMultiDimFet::MakeMethod.
1450   //
1451   // The options are: NONE so far
1452   Int_t i, j;
1453 
1454   Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
1455   const char *prefix = (isMethod ? Form("%s::", classname) : "");
1456   const char *cv_qual = (isMethod ? "" : "static ");
1457 
1458   std::ofstream outFile(filename, std::ios::out | std::ios::trunc);
1459   if (!outFile) {
1460     Error("MakeRealCode", "couldn't open output file '%s'", filename);
1461     return;
1462   }
1463 
1464   if (fIsVerbose)
1465     edm::LogInfo("TMultiDimFet") << "Writing on file \"" << filename << "\" ... " << std::flush;
1466   //
1467   // Write header of file
1468   //
1469   // Emacs mode line ;-)
1470   outFile << "// -*- mode: c++ -*-"
1471           << "\n";
1472   // Info about creator
1473   outFile << "// "
1474           << "\n"
1475           << "// File " << filename << " generated by TMultiDimFet::MakeRealCode"
1476           << "\n";
1477   // Time stamp
1478   TDatime date;
1479   outFile << "// on " << date.AsString() << "\n";
1480   // ROOT version info
1481   outFile << "// ROOT version " << gROOT->GetVersion() << "\n"
1482           << "//"
1483           << "\n";
1484   // General information on the code
1485   outFile << "// This file contains the function "
1486           << "\n"
1487           << "//"
1488           << "\n"
1489           << "//    double  " << prefix << "MDF(double *x); "
1490           << "\n"
1491           << "//"
1492           << "\n"
1493           << "// For evaluating the parameterization obtained"
1494           << "\n"
1495           << "// from TMultiDimFet and the point x"
1496           << "\n"
1497           << "// "
1498           << "\n"
1499           << "// See TMultiDimFet class documentation for more "
1500           << "information "
1501           << "\n"
1502           << "// "
1503           << "\n";
1504   // Header files
1505   if (isMethod)
1506     // If these are methods, we need the class header
1507     outFile << "#include \"" << classname << ".h\""
1508             << "\n";
1509 
1510   //
1511   // Now for the data
1512   //
1513   outFile << "//"
1514           << "\n"
1515           << "// Static data variables"
1516           << "\n"
1517           << "//"
1518           << "\n";
1519   outFile << cv_qual << "int    " << prefix << "gNVariables    = " << fNVariables << ";"
1520           << "\n";
1521   outFile << cv_qual << "int    " << prefix << "gNCoefficients = " << fNCoefficients << ";"
1522           << "\n";
1523   outFile << cv_qual << "double " << prefix << "gDMean         = " << fMeanQuantity << ";"
1524           << "\n";
1525 
1526   // Assignment to mean vector.
1527   outFile << "// Assignment to mean vector."
1528           << "\n";
1529   outFile << cv_qual << "double " << prefix << "gXMean[] = {"
1530           << "\n";
1531   for (i = 0; i < fNVariables; i++)
1532     outFile << (i != 0 ? ", " : "  ") << fMeanVariables(i) << std::flush;
1533   outFile << " };"
1534           << "\n"
1535           << "\n";
1536 
1537   // Assignment to minimum vector.
1538   outFile << "// Assignment to minimum vector."
1539           << "\n";
1540   outFile << cv_qual << "double " << prefix << "gXMin[] = {"
1541           << "\n";
1542   for (i = 0; i < fNVariables; i++)
1543     outFile << (i != 0 ? ", " : "  ") << fMinVariables(i) << std::flush;
1544   outFile << " };"
1545           << "\n"
1546           << "\n";
1547 
1548   // Assignment to maximum vector.
1549   outFile << "// Assignment to maximum vector."
1550           << "\n";
1551   outFile << cv_qual << "double " << prefix << "gXMax[] = {"
1552           << "\n";
1553   for (i = 0; i < fNVariables; i++)
1554     outFile << (i != 0 ? ", " : "  ") << fMaxVariables(i) << std::flush;
1555   outFile << " };"
1556           << "\n"
1557           << "\n";
1558 
1559   // Assignment to coefficients vector.
1560   outFile << "// Assignment to coefficients vector."
1561           << "\n";
1562   outFile << cv_qual << "double " << prefix << "gCoefficient[] = {" << std::flush;
1563   for (i = 0; i < fNCoefficients; i++)
1564     outFile << (i != 0 ? "," : "") << "\n"
1565             << "  " << fCoefficients(i) << std::flush;
1566   outFile << "\n"
1567           << " };"
1568           << "\n"
1569           << "\n";
1570 
1571   // Assignment to powers vector.
1572   outFile << "// Assignment to powers vector."
1573           << "\n"
1574           << "// The powers are stored row-wise, that is"
1575           << "\n"
1576           << "//  p_ij = " << prefix << "gPower[i * NVariables + j];"
1577           << "\n";
1578   outFile << cv_qual << "int    " << prefix << "gPower[] = {" << std::flush;
1579   for (i = 0; i < fNCoefficients; i++) {
1580     for (j = 0; j < fNVariables; j++) {
1581       if (j != 0)
1582         outFile << std::flush << "  ";
1583       else
1584         outFile << "\n"
1585                 << "  ";
1586       outFile << fPowers[fPowerIndex[i] * fNVariables + j]
1587               << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",") << std::flush;
1588     }
1589   }
1590   outFile << "\n"
1591           << "};"
1592           << "\n"
1593           << "\n";
1594 
1595   //
1596   // Finally we reach the function itself
1597   //
1598   outFile << "// "
1599           << "\n"
1600           << "// The " << (isMethod ? "method " : "function ") << "  double " << prefix << "MDF(double *x)"
1601           << "\n"
1602           << "// "
1603           << "\n";
1604   outFile << "double " << prefix << "MDF(double *x) {"
1605           << "\n"
1606           << "  double returnValue = " << prefix << "gDMean;"
1607           << "\n"
1608           << "  int    i = 0, j = 0, k = 0;"
1609           << "\n"
1610           << "  for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
1611           << "\n"
1612           << "    // Evaluate the ith term in the expansion"
1613           << "\n"
1614           << "    double term = " << prefix << "gCoefficient[i];"
1615           << "\n"
1616           << "    for (j = 0; j < " << prefix << "gNVariables; j++) {"
1617           << "\n"
1618           << "      // Evaluate the polynomial in the jth variable."
1619           << "\n"
1620           << "      int power = " << prefix << "gPower[" << prefix << "gNVariables * i + j]; "
1621           << "\n"
1622           << "      double p1 = 1, p2 = 0, p3 = 0, r = 0;"
1623           << "\n"
1624           << "      double v =  1 + 2. / (" << prefix << "gXMax[j] - " << prefix << "gXMin[j]) * (x[j] - " << prefix
1625           << "gXMax[j]);"
1626           << "\n"
1627           << "      // what is the power to use!"
1628           << "\n"
1629           << "      switch(power) {"
1630           << "\n"
1631           << "      case 1: r = 1; break; "
1632           << "\n"
1633           << "      case 2: r = v; break; "
1634           << "\n"
1635           << "      default: "
1636           << "\n"
1637           << "        p2 = v; "
1638           << "\n"
1639           << "        for (k = 3; k <= power; k++) { "
1640           << "\n"
1641           << "          p3 = p2 * v;"
1642           << "\n";
1643   if (fPolyType == kLegendre)
1644     outFile << "          p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)"
1645             << " / (i - 1);"
1646             << "\n";
1647   if (fPolyType == kChebyshev)
1648     outFile << "          p3 = 2 * v * p2 - p1; "
1649             << "\n";
1650   outFile << "          p1 = p2; p2 = p3; "
1651           << "\n"
1652           << "        }"
1653           << "\n"
1654           << "        r = p3;"
1655           << "\n"
1656           << "      }"
1657           << "\n"
1658           << "      // multiply this term by the poly in the jth var"
1659           << "\n"
1660           << "      term *= r; "
1661           << "\n"
1662           << "    }"
1663           << "\n"
1664           << "    // Add this term to the final result"
1665           << "\n"
1666           << "    returnValue += term;"
1667           << "\n"
1668           << "  }"
1669           << "\n"
1670           << "  return returnValue;"
1671           << "\n"
1672           << "}"
1673           << "\n"
1674           << "\n";
1675 
1676   // EOF
1677   outFile << "// EOF for " << filename << "\n";
1678 
1679   // Close the file
1680   outFile.close();
1681 
1682   if (fIsVerbose)
1683     edm::LogInfo("TMultiDimFet") << "done"
1684                                  << "\n";
1685 }
1686 
1687 //____________________________________________________________________
1688 void TMultiDimFet::Print(Option_t *option) const {
1689   // Print statistics etc.
1690   // Options are
1691   //   P        Parameters
1692   //   S        Statistics
1693   //   C        Coefficients
1694   //   R        Result of parameterisation
1695   //   F        Result of fit
1696   //   K        Correlation Matrix
1697   //   M        Pretty print formula
1698   //
1699   Int_t i = 0;
1700   Int_t j = 0;
1701 
1702   TString opt(option);
1703   opt.ToLower();
1704 
1705   if (opt.Contains("p")) {
1706     // Print basic parameters for this object
1707     edm::LogInfo("TMultiDimFet") << "User parameters:"
1708                                  << "\n"
1709                                  << "----------------"
1710                                  << "\n"
1711                                  << " Variables:                    " << fNVariables << "\n"
1712                                  << " Data points:                  " << fSampleSize << "\n"
1713                                  << " Max Terms:                    " << fMaxTerms << "\n"
1714                                  << " Power Limit Parameter:        " << fPowerLimit << "\n"
1715                                  << " Max functions:                " << fMaxFunctions << "\n"
1716                                  << " Max functions to study:       " << fMaxStudy << "\n"
1717                                  << " Max angle (optional):         " << fMaxAngle << "\n"
1718                                  << " Min angle:                    " << fMinAngle << "\n"
1719                                  << " Relative Error accepted:      " << fMinRelativeError << "\n"
1720                                  << " Maximum Powers:               " << std::flush;
1721     for (i = 0; i < fNVariables; i++)
1722       edm::LogInfo("TMultiDimFet") << " " << fMaxPowers[i] - 1 << std::flush;
1723     edm::LogInfo("TMultiDimFet") << "\n"
1724                                  << "\n"
1725                                  << " Parameterisation will be done using " << std::flush;
1726     if (fPolyType == kChebyshev)
1727       edm::LogInfo("TMultiDimFet") << "Chebyshev polynomials"
1728                                    << "\n";
1729     else if (fPolyType == kLegendre)
1730       edm::LogInfo("TMultiDimFet") << "Legendre polynomials"
1731                                    << "\n";
1732     else
1733       edm::LogInfo("TMultiDimFet") << "Monomials"
1734                                    << "\n";
1735     edm::LogInfo("TMultiDimFet") << "\n";
1736   }
1737 
1738   if (opt.Contains("s")) {
1739     // Print statistics for read data
1740     edm::LogInfo("TMultiDimFet") << "Sample statistics:"
1741                                  << "\n"
1742                                  << "------------------"
1743                                  << "\n"
1744                                  << "                 D" << std::flush;
1745     for (i = 0; i < fNVariables; i++)
1746       edm::LogInfo("TMultiDimFet") << " " << std::setw(10) << i + 1 << std::flush;
1747     edm::LogInfo("TMultiDimFet") << "\n"
1748                                  << " Max:   " << std::setw(10) << std::setprecision(7) << fMaxQuantity << std::flush;
1749     for (i = 0; i < fNVariables; i++)
1750       edm::LogInfo("TMultiDimFet") << " " << std::setw(10) << std::setprecision(4) << fMaxVariables(i) << std::flush;
1751     edm::LogInfo("TMultiDimFet") << "\n"
1752                                  << " Min:   " << std::setw(10) << std::setprecision(7) << fMinQuantity << std::flush;
1753     for (i = 0; i < fNVariables; i++)
1754       edm::LogInfo("TMultiDimFet") << " " << std::setw(10) << std::setprecision(4) << fMinVariables(i) << std::flush;
1755     edm::LogInfo("TMultiDimFet") << "\n"
1756                                  << " Mean:  " << std::setw(10) << std::setprecision(7) << fMeanQuantity << std::flush;
1757     for (i = 0; i < fNVariables; i++)
1758       edm::LogInfo("TMultiDimFet") << " " << std::setw(10) << std::setprecision(4) << fMeanVariables(i) << std::flush;
1759     edm::LogInfo("TMultiDimFet") << "\n"
1760                                  << " Function Sum Squares:         " << fSumSqQuantity << "\n"
1761                                  << "\n";
1762   }
1763 
1764   if (opt.Contains("r")) {
1765     edm::LogInfo("TMultiDimFet") << "Results of Parameterisation:"
1766                                  << "\n"
1767                                  << "----------------------------"
1768                                  << "\n"
1769                                  << " Total reduction of square residuals    " << fSumSqResidual << "\n"
1770                                  << " Relative precision obtained:           " << fPrecision << "\n"
1771                                  << " Error obtained:                        " << fError << "\n"
1772                                  << " Multiple correlation coefficient:      " << fCorrelationCoeff << "\n"
1773                                  << " Reduced Chi square over sample:        " << fChi2 / (fSampleSize - fNCoefficients)
1774                                  << "\n"
1775                                  << " Maximum residual value:                " << fMaxResidual << "\n"
1776                                  << " Minimum residual value:                " << fMinResidual << "\n"
1777                                  << " Estimated root mean square:            " << fRMS << "\n"
1778                                  << " Maximum powers used:                   " << std::flush;
1779     for (j = 0; j < fNVariables; j++)
1780       edm::LogInfo("TMultiDimFet") << fMaxPowersFinal[j] << " " << std::flush;
1781     edm::LogInfo("TMultiDimFet") << "\n"
1782                                  << " Function codes of candidate functions."
1783                                  << "\n"
1784                                  << "  1: considered,"
1785                                  << "  2: too little contribution,"
1786                                  << "  3: accepted." << std::flush;
1787     for (i = 0; i < fMaxFunctions; i++) {
1788       if (i % 60 == 0)
1789         edm::LogInfo("TMultiDimFet") << "\n"
1790                                      << " " << std::flush;
1791       else if (i % 10 == 0)
1792         edm::LogInfo("TMultiDimFet") << " " << std::flush;
1793       edm::LogInfo("TMultiDimFet") << fFunctionCodes[i];
1794     }
1795     edm::LogInfo("TMultiDimFet") << "\n"
1796                                  << " Loop over candidates stopped because " << std::flush;
1797     switch (fParameterisationCode) {
1798       case PARAM_MAXSTUDY:
1799         edm::LogInfo("TMultiDimFet") << "max allowed studies reached"
1800                                      << "\n";
1801         break;
1802       case PARAM_SEVERAL:
1803         edm::LogInfo("TMultiDimFet") << "all candidates considered several times"
1804                                      << "\n";
1805         break;
1806       case PARAM_RELERR:
1807         edm::LogInfo("TMultiDimFet") << "wanted relative error obtained"
1808                                      << "\n";
1809         break;
1810       case PARAM_MAXTERMS:
1811         edm::LogInfo("TMultiDimFet") << "max number of terms reached"
1812                                      << "\n";
1813         break;
1814       default:
1815         edm::LogInfo("TMultiDimFet") << "some unknown reason"
1816                                      << "\n";
1817         break;
1818     }
1819     edm::LogInfo("TMultiDimFet") << "\n";
1820   }
1821 
1822   if (opt.Contains("f")) {
1823     edm::LogInfo("TMultiDimFet") << "Results of Fit:"
1824                                  << "\n"
1825                                  << "---------------"
1826                                  << "\n"
1827                                  << " Test sample size:                      " << fTestSampleSize << "\n"
1828                                  << " Multiple correlation coefficient:      " << fTestCorrelationCoeff << "\n"
1829                                  << " Relative precision obtained:           " << fTestPrecision << "\n"
1830                                  << " Error obtained:                        " << fTestError << "\n"
1831                                  << " Reduced Chi square over sample:        " << fChi2 / (fSampleSize - fNCoefficients)
1832                                  << "\n"
1833                                  << "\n";
1834     /*
1835       if (fFitter) {
1836          fFitter->PrintResults(1,1);
1837          edm::LogInfo("TMultiDimFet") << "\n";
1838       }
1839 */
1840   }
1841 
1842   if (opt.Contains("c")) {
1843     edm::LogInfo("TMultiDimFet") << "Coefficients:"
1844                                  << "\n"
1845                                  << "-------------"
1846                                  << "\n"
1847                                  << "   #         Value        Error   Powers"
1848                                  << "\n"
1849                                  << " ---------------------------------------"
1850                                  << "\n";
1851     for (i = 0; i < fNCoefficients; i++) {
1852       edm::LogInfo("TMultiDimFet") << " " << std::setw(3) << i << "  " << std::setw(12) << fCoefficients(i) << "  "
1853                                    << std::setw(12) << fCoefficientsRMS(i) << "  " << std::flush;
1854       for (j = 0; j < fNVariables; j++)
1855         edm::LogInfo("TMultiDimFet") << " " << std::setw(3) << fPowers[fPowerIndex[i] * fNVariables + j] - 1
1856                                      << std::flush;
1857       edm::LogInfo("TMultiDimFet") << "\n";
1858     }
1859     edm::LogInfo("TMultiDimFet") << "\n";
1860   }
1861   if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
1862     edm::LogInfo("TMultiDimFet") << "Correlation Matrix:"
1863                                  << "\n"
1864                                  << "-------------------";
1865     fCorrelationMatrix.Print();
1866   }
1867 
1868   if (opt.Contains("m")) {
1869     edm::LogInfo("TMultiDimFet") << std::setprecision(25);
1870     edm::LogInfo("TMultiDimFet") << "Parameterization:"
1871                                  << "\n"
1872                                  << "-----------------"
1873                                  << "\n"
1874                                  << "  Normalised variables: "
1875                                  << "\n";
1876     for (i = 0; i < fNVariables; i++)
1877       edm::LogInfo("TMultiDimFet") << "\ty" << i << "\t:= 1 + 2 * (x" << i << " - " << fMaxVariables(i) << ") / ("
1878                                    << fMaxVariables(i) << " - " << fMinVariables(i) << ")"
1879                                    << "\n";
1880     edm::LogInfo("TMultiDimFet") << "\n"
1881                                  << "  f[";
1882     for (i = 0; i < fNVariables; i++) {
1883       edm::LogInfo("TMultiDimFet") << "y" << i;
1884       if (i != fNVariables - 1)
1885         edm::LogInfo("TMultiDimFet") << ", ";
1886     }
1887     edm::LogInfo("TMultiDimFet") << "] := ";
1888     for (Int_t i = 0; i < fNCoefficients; i++) {
1889       if (i != 0)
1890         edm::LogInfo("TMultiDimFet") << " " << (fCoefficients(i) < 0 ? "- " : "+ ") << TMath::Abs(fCoefficients(i));
1891       else
1892         edm::LogInfo("TMultiDimFet") << fCoefficients(i);
1893       for (Int_t j = 0; j < fNVariables; j++) {
1894         Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
1895         switch (p) {
1896           case 1:
1897             break;
1898           case 2:
1899             edm::LogInfo("TMultiDimFet") << " * y" << j;
1900             break;
1901           default:
1902             switch (fPolyType) {
1903               case kLegendre:
1904                 edm::LogInfo("TMultiDimFet") << " * L" << p - 1 << "(y" << j << ")";
1905                 break;
1906               case kChebyshev:
1907                 edm::LogInfo("TMultiDimFet") << " * C" << p - 1 << "(y" << j << ")";
1908                 break;
1909               default:
1910                 edm::LogInfo("TMultiDimFet") << " * y" << j << "^" << p - 1;
1911                 break;
1912             }
1913         }
1914       }
1915     }
1916     edm::LogInfo("TMultiDimFet") << "\n";
1917   }
1918 }
1919 
1920 //____________________________________________________________________
1921 void TMultiDimFet::PrintPolynomialsSpecial(Option_t *option) const {
1922   //   M        Pretty print formula
1923   //
1924   Int_t i = 0;
1925   //   Int_t j = 0;
1926 
1927   TString opt(option);
1928   opt.ToLower();
1929 
1930   if (opt.Contains("m")) {
1931     edm::LogInfo("TMultiDimFet") << std::setprecision(25);
1932     edm::LogInfo("TMultiDimFet") << "Parameterization:"
1933                                  << "\n"
1934                                  << "-----------------"
1935                                  << "\n"
1936                                  << "  Normalised variables: "
1937                                  << "\n";
1938     for (i = 0; i < fNVariables; i++)
1939       edm::LogInfo("TMultiDimFet") << "\tdouble y" << i << "\t=1+2*(x" << i << "-" << fMaxVariables(i) << ")/("
1940                                    << fMaxVariables(i) << "-" << fMinVariables(i) << ");"
1941                                    << "\n";
1942     edm::LogInfo("TMultiDimFet") << "\n"
1943                                  << "  f[";
1944     for (i = 0; i < fNVariables; i++) {
1945       edm::LogInfo("TMultiDimFet") << "y" << i;
1946       if (i != fNVariables - 1)
1947         edm::LogInfo("TMultiDimFet") << ", ";
1948     }
1949     edm::LogInfo("TMultiDimFet") << "] := " << fMeanQuantity << " + ";
1950     for (Int_t i = 0; i < fNCoefficients; i++) {
1951       if (i != 0)
1952         edm::LogInfo("TMultiDimFet") << " " << (fCoefficients(i) < 0 ? "-" : "+") << TMath::Abs(fCoefficients(i));
1953       else
1954         edm::LogInfo("TMultiDimFet") << fCoefficients(i);
1955       for (Int_t j = 0; j < fNVariables; j++) {
1956         Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
1957         switch (p) {
1958           case 1:
1959             break;
1960           case 2:
1961             edm::LogInfo("TMultiDimFet") << "*y" << j;
1962             break;
1963           default:
1964             switch (fPolyType) {
1965               case kLegendre:
1966                 edm::LogInfo("TMultiDimFet") << "*Leg(" << p - 1 << ",y" << j << ")";
1967                 break;
1968               case kChebyshev:
1969                 edm::LogInfo("TMultiDimFet") << "*C" << p - 1 << "(y" << j << ")";
1970                 break;
1971               default:
1972                 edm::LogInfo("TMultiDimFet") << "*y" << j << "**" << p - 1;
1973                 break;
1974             }
1975         }
1976       }
1977       edm::LogInfo("TMultiDimFet") << "\n";
1978     }
1979     edm::LogInfo("TMultiDimFet") << "\n";
1980   }
1981 }
1982 
1983 //____________________________________________________________________
1984 Bool_t TMultiDimFet::Select(const Int_t *) {
1985   // Selection method. User can override this method for specialized
1986   // selection of acceptable functions in fit. Default is to select
1987   // all. This message is sent during the build-up of the function
1988   // candidates table once for each set of powers in
1989   // variables. Notice, that the argument array contains the powers
1990   // PLUS ONE. For example, to De select the function
1991   //     f = x1^2 * x2^4 * x3^5,
1992   // this method should return kFALSE if given the argument
1993   //     { 3, 4, 6 }
1994   return kTRUE;
1995 }
1996 
1997 //____________________________________________________________________
1998 void TMultiDimFet::SetMaxAngle(Double_t ang) {
1999   // Set the max angle (in degrees) between the initial data vector to
2000   // be fitted, and the new candidate function to be included in the
2001   // fit.  By default it is 0, which automatically chooses another
2002   // selection criteria. See also
2003   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2004   if (ang >= 90 || ang < 0) {
2005     Warning("SetMaxAngle", "angle must be in [0,90)");
2006     return;
2007   }
2008 
2009   fMaxAngle = ang;
2010 }
2011 
2012 //____________________________________________________________________
2013 void TMultiDimFet::SetMinAngle(Double_t ang) {
2014   // Set the min angle (in degrees) between a new candidate function
2015   // and the subspace spanned by the previously accepted
2016   // functions. See also
2017   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2018   if (ang > 90 || ang <= 0) {
2019     Warning("SetMinAngle", "angle must be in [0,90)");
2020     return;
2021   }
2022 
2023   fMinAngle = ang;
2024 }
2025 
2026 //____________________________________________________________________
2027 void TMultiDimFet::SetPowers(const Int_t *powers, Int_t terms) {
2028   // Define a user function. The input array must be of the form
2029   // (p11, ..., p1N, ... ,pL1, ..., pLN)
2030   // Where N is the dimension of the data sample, L is the number of
2031   // terms (given in terms) and the first number, labels the term, the
2032   // second the variable.  More information is given in the
2033   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2034   fIsUserFunction = kTRUE;
2035   fMaxFunctions = terms;
2036   fMaxTerms = terms;
2037   fMaxStudy = terms;
2038   fMaxFunctionsTimesNVariables = fMaxFunctions * fNVariables;
2039   fPowers.resize(fMaxFunctions * fNVariables);
2040   Int_t i, j;
2041   for (i = 0; i < fMaxFunctions; i++)
2042     for (j = 0; j < fNVariables; j++)
2043       fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1;
2044 }
2045 
2046 //____________________________________________________________________
2047 void TMultiDimFet::SetPowerLimit(Double_t limit) {
2048   // Set the user parameter for the function selection. The bigger the
2049   // limit, the more functions are used. The meaning of this variable
2050   // is defined in the
2051   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2052   fPowerLimit = limit;
2053 }
2054 
2055 //____________________________________________________________________
2056 void TMultiDimFet::SetMaxPowers(const Int_t *powers) {
2057   // Set the maximum power to be considered in the fit for each
2058   // variable. See also
2059   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2060   if (!powers)
2061     return;
2062 
2063   for (Int_t i = 0; i < fNVariables; i++)
2064     fMaxPowers[i] = powers[i] + 1;
2065 }
2066 
2067 //____________________________________________________________________
2068 void TMultiDimFet::SetMinRelativeError(Double_t error) {
2069   // Set the acceptable relative error for when sum of square
2070   // residuals is considered minimized. For a full account, refer to
2071   // the
2072   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2073   fMinRelativeError = error;
2074 }
2075 
2076 //____________________________________________________________________
2077 Bool_t TMultiDimFet::TestFunction(Double_t squareResidual, Double_t dResidur) {
2078   // PRIVATE METHOD:
2079   // Test whether the currently considered function contributes to the
2080   // fit. See also
2081   // Begin_Html<a href="#TMultiDimFet:description">class description</a>End_Html
2082 
2083   if (fNCoefficients != 0) {
2084     // Now for the second test:
2085     if (fMaxAngle == 0) {
2086       // If the user hasn't supplied a max angle do the test as,
2087       if (dResidur < squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2088         return kFALSE;
2089       }
2090     } else {
2091       // If the user has provided a max angle, test if the calculated
2092       // angle is less then the max angle.
2093       if (TMath::Sqrt(dResidur / fSumSqAvgQuantity) < TMath::Cos(fMaxAngle * DEGRAD)) {
2094         return kFALSE;
2095       }
2096     }
2097   }
2098   // If we get here, the function is OK
2099   return kTRUE;
2100 }
2101 
2102 //____________________________________________________________________
2103 //void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
2104 //               double* coeffs, int /*flag*/)
2105 /*{
2106    // Helper function for doing the minimisation of Chi2 using Minuit
2107 
2108    // Get pointer  to current TMultiDimFet object.
2109   // TMultiDimFet* mdf = TMultiDimFet::Instance();
2110    //chi2     = mdf->MakeChi2(coeffs);
2111 }
2112 */