File indexing completed on 2021-02-14 14:30:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
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
0042
0043
0044 ClassImp(TMultiDimFet);
0045
0046
0047
0048
0049
0050
0051 TMultiDimFet::TMultiDimFet() {
0052
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
0080
0081 }
0082
0083 const TMultiDimFet &TMultiDimFet::operator=(const TMultiDimFet &in) {
0084 if (this == &in) {
0085 return in;
0086 }
0087
0088 fMeanQuantity = in.fMeanQuantity;
0089
0090 fMaxQuantity = 0.0;
0091 fMinQuantity = 0.0;
0092 fSumSqQuantity = 0.0;
0093 fSumSqAvgQuantity = 0.0;
0094
0095 fNVariables = in.fNVariables;
0096
0097 fMaxVariables.ResizeTo(in.fMaxVariables.GetLwb(), in.fMaxVariables.GetUpb());
0098 fMaxVariables = in.fMaxVariables;
0099
0100 fMinVariables.ResizeTo(in.fMinVariables.GetLwb(), in.fMinVariables.GetUpb());
0101 fMinVariables = in.fMinVariables;
0102
0103 fSampleSize = 0;
0104 fTestSampleSize = 0;
0105 fMinAngle = 1;
0106 fMaxAngle = 0.0;
0107
0108 fMaxTerms = in.fMaxTerms;
0109
0110 fMinRelativeError = 0.0;
0111
0112 fMaxPowers.clear();
0113 fPowerLimit = 1;
0114
0115 fMaxFunctions = in.fMaxFunctions;
0116
0117 fFunctionCodes.clear();
0118 fMaxStudy = 0;
0119
0120 fMaxPowersFinal.clear();
0121
0122 fMaxFunctionsTimesNVariables = in.fMaxFunctionsTimesNVariables;
0123 fPowers = in.fPowers;
0124
0125 fPowerIndex = in.fPowerIndex;
0126
0127 fMaxResidual = 0.0;
0128 fMinResidual = 0.0;
0129 fMaxResidualRow = 0;
0130 fMinResidualRow = 0;
0131 fSumSqResidual = 0.0;
0132
0133 fNCoefficients = in.fNCoefficients;
0134
0135 fCoefficients.ResizeTo(in.fCoefficients.GetLwb(), in.fCoefficients.GetUpb());
0136 fCoefficients = in.fCoefficients;
0137
0138 fRMS = 0.0;
0139 fChi2 = 0.0;
0140 fParameterisationCode = 0;
0141 fError = 0.0;
0142 fTestError = 0.0;
0143 fPrecision = 0.0;
0144 fTestPrecision = 0.0;
0145 fCorrelationCoeff = 0.0;
0146 fTestCorrelationCoeff = 0.0;
0147 fHistograms = nullptr;
0148 fHistogramMask = 0;
0149
0150
0151 fPolyType = in.fPolyType;
0152 fShowCorrelation = in.fShowCorrelation;
0153 fIsUserFunction = in.fIsUserFunction;
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
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
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
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
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
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
0267
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
0275 fQuantity(fSampleSize - 1) = D;
0276 fSqError(fSampleSize - 1) = (E == 0 ? D : E);
0277
0278
0279
0280
0281 size = fVariables.GetNrows();
0282 if (fSampleSize * fNVariables > size)
0283 fVariables.ResizeTo(size + size / 2);
0284
0285
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
0297 if (x[i] >= fMaxVariables(i))
0298 fMaxVariables(i) = x[i];
0299
0300
0301 if (x[i] <= fMinVariables(i))
0302 fMinVariables(i) = x[i];
0303 }
0304
0305
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
0314
0315
0316
0317
0318
0319
0320
0321
0322 if (fTestSampleSize++ == 0) {
0323 fTestQuantity.ResizeTo(fNVariables);
0324 fTestSqError.ResizeTo(fNVariables);
0325 fTestVariables.ResizeTo(fNVariables * 100);
0326 }
0327
0328
0329
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
0337 fTestQuantity(fTestSampleSize - 1) = D;
0338 fTestSqError(fTestSampleSize - 1) = (E == 0 ? D : E);
0339
0340
0341
0342
0343 size = fTestVariables.GetNrows();
0344 if (fTestSampleSize * fNVariables > size)
0345 fTestVariables.ResizeTo(size + size / 2);
0346
0347
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
0363 Int_t i, j, n = fNVariables, m = fMaxFunctions;
0364
0365
0366 fQuantity.Zero();
0367 fSqError.Zero();
0368 fMeanQuantity = 0;
0369 fMaxQuantity = 0;
0370 fMinQuantity = 0;
0371 fSumSqQuantity = 0;
0372 fSumSqAvgQuantity = 0;
0373
0374
0375 fVariables.Zero();
0376 fNVariables = 0;
0377 fSampleSize = 0;
0378 fMeanVariables.Zero();
0379 fMaxVariables.Zero();
0380 fMinVariables.Zero();
0381
0382
0383 fTestQuantity.Zero();
0384 fTestSqError.Zero();
0385 fTestVariables.Zero();
0386 fTestSampleSize = 0;
0387
0388
0389 fFunctions.Zero();
0390 fMaxFunctions = 0;
0391 fMaxStudy = 0;
0392 fMaxFunctionsTimesNVariables = 0;
0393 fOrthFunctions.Zero();
0394 fOrthFunctionNorms.Zero();
0395
0396
0397 fMinRelativeError = 0;
0398 fMinAngle = 0;
0399 fMaxAngle = 0;
0400 fMaxTerms = 0;
0401
0402
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
0412 fMaxResidual = 0;
0413 fMinResidual = 0;
0414 fMaxResidualRow = 0;
0415 fMinResidualRow = 0;
0416 fSumSqResidual = 0;
0417
0418
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
0430 fCoefficients.Zero();
0431 fCoefficientsRMS.Zero();
0432 fResiduals.Zero();
0433 fHistograms->Clear(option);
0434
0435
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
0444
0445
0446
0447 Double_t returnValue = fMeanQuantity;
0448 Double_t term = 0;
0449 Int_t i, j;
0450
0451 for (i = 0; i < fNCoefficients; i++) {
0452
0453 term = (coeff ? coeff[i] : fCoefficients(i));
0454 for (j = 0; j < fNVariables; j++) {
0455
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
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
0517
0518 Double_t s = 0;
0519 Double_t epsilon = 1e-6;
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
0530
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
0564
0565
0566
0567
0568
0569
0570 MakeNormalized();
0571 MakeCandidates();
0572 MakeParameterization();
0573 MakeCoefficients();
0574 ReducePolynomial(precision);
0575 }
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656
0657
0658 void TMultiDimFet::MakeCandidates() {
0659
0660
0661
0662
0663 Int_t i = 0;
0664 Int_t j = 0;
0665 Int_t k = 0;
0666
0667
0668
0669 Int_t *powers = new Int_t[fNVariables * fMaxFunctions];
0670
0671
0672 Double_t *control = new Double_t[fMaxFunctions];
0673
0674
0675 Int_t *iv = new Int_t[fNVariables];
0676 for (i = 0; i < fNVariables; i++)
0677 iv[i] = 1;
0678
0679 if (!fIsUserFunction) {
0680
0681 Int_t numberFunctions = 0;
0682
0683
0684 Int_t maxNumberFunctions = 1;
0685 for (i = 0; i < fNVariables; i++)
0686 maxNumberFunctions *= fMaxPowers[i];
0687
0688 while (kTRUE) {
0689
0690 Double_t s = EvalControl(iv);
0691
0692 if (s <= fPowerLimit) {
0693
0694
0695 if (Select(iv)) {
0696 numberFunctions++;
0697
0698
0699
0700 if (numberFunctions > fMaxFunctions)
0701 break;
0702
0703
0704
0705 control[numberFunctions - 1] = Int_t(1.0e+6 * s);
0706
0707
0708 for (i = 0; i < fNVariables; i++) {
0709 j = (numberFunctions - 1) * fNVariables + i;
0710 powers[j] = iv[i];
0711 }
0712 }
0713 }
0714
0715 for (i = 0; i < fNVariables; i++)
0716 if (iv[i] < fMaxPowers[i])
0717 break;
0718
0719
0720
0721 if (i == fNVariables) {
0722 fMaxFunctions = numberFunctions;
0723 fMaxFunctionsTimesNVariables = fMaxFunctions * fNVariables;
0724 break;
0725 }
0726
0727
0728 iv[i]++;
0729
0730 for (j = 0; j < i; j++)
0731 iv[j] = 1;
0732 }
0733 } else {
0734
0735 for (i = 0; i < fMaxFunctions; i++) {
0736
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
0747
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
0786
0787
0788
0789
0790 fChi2 = 0;
0791 Int_t i, j;
0792 Double_t *x = new Double_t[fNVariables];
0793 for (i = 0; i < fTestSampleSize; i++) {
0794
0795 for (j = 0; j < fNVariables; j++)
0796 x[j] = fTestVariables(i * fNVariables + j);
0797
0798
0799 Double_t f = Eval(x, coeff);
0800
0801
0802 fChi2 += 1. / TMath::Max(fTestSqError(i), 1e-20) * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
0803 }
0804
0805
0806 delete[] x;
0807
0808 return fChi2;
0809 }
0810
0811
0812 void TMultiDimFet::MakeCode(const char *filename, Option_t *option) {
0813
0814
0815
0816
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
0832
0833
0834
0835
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
0846
0847
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
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
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
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
0895
0896
0897
0898
0899
0900
0901
0902 Int_t i = 0, j = 0;
0903 Int_t col = 0, row = 0;
0904
0905
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
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
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
0934
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
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 }
0966 }
0967
0968
0969 void TMultiDimFet::MakeCorrelation() {
0970
0971
0972 if (!fShowCorrelation)
0973 return;
0974
0975 fCorrelationMatrix.ResizeTo(fNVariables, fNVariables + 1);
0976
0977 Double_t d2 = 0;
0978 Double_t ddotXi = 0;
0979 Double_t xiNorm = 0;
0980 Double_t xidotXj = 0;
0981 Double_t xjNorm = 0;
0982
0983 Int_t i, j, k, l, m;
0984 for (i = 0; i < fSampleSize; i++)
0985 d2 += fQuantity(i) * fQuantity(i);
0986
0987 for (i = 0; i < fNVariables; i++) {
0988 ddotXi = 0.;
0989 xiNorm = 0.;
0990 for (j = 0; j < fSampleSize; j++) {
0991
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.;
1000 xjNorm = 0.;
1001 for (k = 0; k < fSampleSize; k++) {
1002
1003
1004 l = k * fNVariables + j;
1005 m = k * fNVariables + i;
1006
1007
1008 xidotXj +=
1009 (fVariables(m) - fMeanVariables(i)) * (fVariables(l) - fMeanVariables(j));
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
1020
1021
1022
1023
1024
1025
1026
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
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
1044 f2 += fFunctions(fNCoefficients, j) * fFunctions(fNCoefficients, j);
1045
1046 fOrthFunctions(fNCoefficients, j) = fFunctions(fNCoefficients, j);
1047 }
1048
1049
1050 for (j = 0; j < fNCoefficients; j++) {
1051 Double_t fdw = 0;
1052
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
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
1065 fOrthFunctionNorms(fNCoefficients) += fOrthFunctions(fNCoefficients, j) * fOrthFunctions(fNCoefficients, j);
1066
1067
1068 fOrthCoefficients(fNCoefficients) += fQuantity(j) * fOrthFunctions(fNCoefficients, j);
1069 }
1070
1071
1072 if (!fIsUserFunction)
1073 if (TMath::Sqrt(fOrthFunctionNorms(fNCoefficients) / (f2 + 1e-10)) < TMath::Sin(fMinAngle * DEGRAD))
1074 return 0;
1075
1076
1077
1078
1079
1080 fOrthCurvatureMatrix(fNCoefficients, fNCoefficients) = 1;
1081 Double_t b = fOrthCoefficients(fNCoefficients);
1082 fOrthCoefficients(fNCoefficients) /= fOrthFunctionNorms(fNCoefficients);
1083
1084
1085 Double_t dResidur = fOrthCoefficients(fNCoefficients) * b;
1086
1087 return dResidur;
1088 }
1089
1090
1091 void TMultiDimFet::MakeHistograms(Option_t *option) {
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
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
1119 Int_t i = 0;
1120
1121
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
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
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
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
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
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
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
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249 MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
1250 }
1251
1252
1253 void TMultiDimFet::MakeNormalized() {
1254
1255
1256
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
1277 if (TESTBIT(fHistogramMask, HIST_XORIG))
1278 ((TH1D *)fHistograms->FindObject(Form("x_%d_orig", j)))->Fill(fVariables(k));
1279
1280
1281 fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1282
1283
1284 if (TESTBIT(fHistogramMask, HIST_XNORM))
1285 ((TH1D *)fHistograms->FindObject(Form("x_%d_norm", j)))->Fill(fVariables(k));
1286 }
1287 }
1288
1289 fMaxQuantity -= fMeanQuantity;
1290 fMinQuantity -= fMeanQuantity;
1291
1292
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
1302
1303
1304
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
1335
1336 while (kTRUE) {
1337
1338 if (studied++ >= fMaxStudy) {
1339 fParameterisationCode = PARAM_MAXSTUDY;
1340 break;
1341 }
1342
1343
1344 if (k >= maxPass) {
1345 fParameterisationCode = PARAM_SEVERAL;
1346 break;
1347 }
1348
1349
1350 i++;
1351
1352
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
1367 if (fIsVerbose && studied == 1)
1368 edm::LogInfo("TMultiDimFet") << "Coeff SumSqRes Contrib Angle QM Func"
1369 << " Value W^2 Powers"
1370 << "\n";
1371
1372
1373 Double_t dResidur = MakeGramSchmidt(i);
1374
1375 if (dResidur == 0) {
1376
1377
1378 fFunctionCodes[i] = 1;
1379 continue;
1380 }
1381
1382
1383 if (!fIsUserFunction) {
1384
1385 fFunctionCodes[i] = 2;
1386
1387
1388 if (!TestFunction(squareResidual, dResidur)) {
1389 fFunctionCodes[i] = 1;
1390 continue;
1391 }
1392 }
1393
1394
1395
1396
1397
1398 fFunctionCodes[i] = 3;
1399 fPowerIndex[fNCoefficients] = i;
1400 fNCoefficients++;
1401
1402
1403
1404 squareResidual -= dResidur;
1405
1406
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
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 ) {
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
1447
1448
1449
1450
1451
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
1468
1469
1470 outFile << "// -*- mode: c++ -*-"
1471 << "\n";
1472
1473 outFile << "// "
1474 << "\n"
1475 << "// File " << filename << " generated by TMultiDimFet::MakeRealCode"
1476 << "\n";
1477
1478 TDatime date;
1479 outFile << "// on " << date.AsString() << "\n";
1480
1481 outFile << "// ROOT version " << gROOT->GetVersion() << "\n"
1482 << "//"
1483 << "\n";
1484
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
1505 if (isMethod)
1506
1507 outFile << "#include \"" << classname << ".h\""
1508 << "\n";
1509
1510
1511
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
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
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
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
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
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
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
1677 outFile << "// EOF for " << filename << "\n";
1678
1679
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
1690
1691
1692
1693
1694
1695
1696
1697
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
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
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
1836
1837
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
1923
1924 Int_t i = 0;
1925
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
1986
1987
1988
1989
1990
1991
1992
1993
1994 return kTRUE;
1995 }
1996
1997
1998 void TMultiDimFet::SetMaxAngle(Double_t ang) {
1999
2000
2001
2002
2003
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
2015
2016
2017
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
2029
2030
2031
2032
2033
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
2049
2050
2051
2052 fPowerLimit = limit;
2053 }
2054
2055
2056 void TMultiDimFet::SetMaxPowers(const Int_t *powers) {
2057
2058
2059
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
2070
2071
2072
2073 fMinRelativeError = error;
2074 }
2075
2076
2077 Bool_t TMultiDimFet::TestFunction(Double_t squareResidual, Double_t dResidur) {
2078
2079
2080
2081
2082
2083 if (fNCoefficients != 0) {
2084
2085 if (fMaxAngle == 0) {
2086
2087 if (dResidur < squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2088 return kFALSE;
2089 }
2090 } else {
2091
2092
2093 if (TMath::Sqrt(dResidur / fSumSqAvgQuantity) < TMath::Cos(fMaxAngle * DEGRAD)) {
2094 return kFALSE;
2095 }
2096 }
2097 }
2098
2099 return kTRUE;
2100 }
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112