diff --git a/3j.cc b/3j.cc index 9c62913..8d68984 100644 --- a/3j.cc +++ b/3j.cc @@ -5,7 +5,7 @@ using namespace std; -#include "TFitterMinuit.h" +#include "Math/Minimizer.h" #include "Math/SpecFunc.h" #include "wave.h" @@ -223,13 +223,13 @@ decomposeMoment(int L, int M, const waveset& ws, const double* x) double decomposeMomentError(const std::pair& LM, - const waveset& ws, const TFitterMinuit* minuit) + const waveset& ws, const ROOT::Math::Minimizer* minuit) { return decomposeMomentError(LM.first, LM.second, ws, minuit); } double -decomposeMomentError(int L, int M, const waveset& ws, const TFitterMinuit* minuit) +decomposeMomentError(int L, int M, const waveset& ws, const ROOT::Math::Minimizer* minuit) { double resultSquare = 0; for (size_t iWs = 0; iWs < ws.size(); iWs++) @@ -245,27 +245,27 @@ decomposeMomentError(int L, int M, const waveset& ws, const TFitterMinuit* minui const wave& w2 = w[iW2]; double coeff = getCoefficient(eps, L, M, w1.l, w1.m, w2.l, w2.m); - double re1 = minuit->GetParameter(w1.getIndex()); - double im1 = minuit->GetParameter(w1.getIndex() + 1); - double re2 = minuit->GetParameter(w2.getIndex()); - double im2 = minuit->GetParameter(w2.getIndex() + 1); + double re1 = minuit->X()[w1.getIndex()]; + double im1 = minuit->X()[w1.getIndex() + 1]; + double re2 = minuit->X()[w2.getIndex()]; + double im2 = minuit->X()[w2.getIndex() + 1]; - double errRe1 = minuit->GetParError(w1.getIndex()); - double errIm1 = minuit->GetParError(w1.getIndex() + 1); - double errRe2 = minuit->GetParError(w2.getIndex()); - double errIm2 = minuit->GetParError(w2.getIndex() + 1); + double errRe1 = 0;//minuit->GetParError(w1.getIndex()); + double errIm1 = 0;//minuit->GetParError(w1.getIndex() + 1); + double errRe2 = 0;//minuit->GetParError(w2.getIndex()); + double errIm2 = 0;//minuit->GetParError(w2.getIndex() + 1); - double covReRe = minuit->GetCovarianceMatrixElement(w1.idxInCovariance(minuit), - w2.idxInCovariance(minuit)); + double covReRe = minuit->CovMatrix(w1.getIndex(), + w2.getIndex()); // The covariance between the imaginary parts is only // non-zero if neither is fixed. double covImIm = 0; - if (!minuit->IsFixed(w1.getIndex() + 1) - && !minuit->IsFixed(w2.getIndex() + 1)) + if (!minuit->IsFixedVariable(w1.getIndex() + 1) + && !minuit->IsFixedVariable(w2.getIndex() + 1)) { - covImIm = minuit->GetCovarianceMatrixElement(w1.idxInCovariance(minuit) + 1, - w2.idxInCovariance(minuit) + 1); + covImIm = minuit->CovMatrix(w1.getIndex() + 1, + w2.getIndex() + 1); } // value += coeff*(x[w1.getIndex()]*x[w2.getIndex()] - x[w1.getIndex()+1]*x[w1.getIndex()+1]); diff --git a/3j.h b/3j.h index d2f17ac..b26efba 100644 --- a/3j.h +++ b/3j.h @@ -33,8 +33,8 @@ double decomposeMoment(int L, int M, const waveset& ws, const vector& x) double decomposeMoment(int L, int M, const waveset& ws, const double* x); double decomposeMomentError(const std::pair& LM, - const waveset& ws, const TFitterMinuit* minuit); -double decomposeMomentError(int L, int M, const waveset& ws, const TFitterMinuit* minuit); + const waveset& ws, const ROOT::Math::Minimizer* minuit); +double decomposeMomentError(int L, int M, const waveset& ws, const ROOT::Math::Minimizer* minuit); #endif diff --git a/massDep.cc b/massDep.cc index d1d6327..4839dd3 100644 --- a/massDep.cc +++ b/massDep.cc @@ -1,12 +1,16 @@ #include #include +#include using namespace std; #include "TObject.h" #include "TStopwatch.h" #include "TFile.h" #include "TTree.h" -#include "TFitterMinuit.h" +#include "Minuit2/FCNBase.h" +#include "Math/Minimizer.h" +#include "Math/Factory.h" +#include "Math/Functor.h" #include "TVectorD.h" #include "TMatrixDSym.h" #include "TDecompChol.h" @@ -25,8 +29,13 @@ class chiSquare : public ROOT::Minuit2::FCNBase { double Up() const { return 1.; } double operator()(const vector& x) const; + + double operator()( const double* x ) const { + std::vector p(x, x+info->getNvars()); + return (*this)(p); + } - ~chiSquare() { delete model; } + //~chiSquare() { delete model; } private: TTree *tree; @@ -218,6 +227,8 @@ int main(int argc, char **argv) return 1; } + cout << info->getModelName() << endl; + TRandom1* gRandom = new TRandom1(); gRandom->SetSeed2(0); // truely random!!! @@ -257,8 +268,9 @@ int main(int argc, char **argv) chiSquare fitFunc(tree); - TFitterMinuit* minuit = new TFitterMinuit(); - minuit->SetMinuitFCN(&fitFunc); + ROOT::Math::Minimizer* minuit = ROOT::Math::Factory::CreateMinimizer("Minuit2", ""); + ROOT::Math::Functor func(fitFunc, info->getNvars()); + minuit->SetFunction(func); double vals[] = //{ 50, 50 , 1.275, 0.185, 50, 50, 1.515, 0.07, 50, 50, 2.13, .27, 200, 0, 1, .065, 50, 50, 1.5, 0.104, 0, 100, 1.730, 0.100, 50, 50, 1.37, .3, 50, 50, 0.5, 0.5, -0.5, 50, 50, 0.5, 0.5, -0.5} //WA102 //{ 50, 50, 1.275, 0.185, 50, 50, 1.32, .107, 50, 50, 1.515, 0.07, 200, 0, 1, .065, 100, 0, 1.5, 0.104, 100, 0, 1.730, 0.100, 50, 50, 0.5, -0.5, 0.5, 50, 50, 0.5, -0.5, 0.5} // f2/a2 @@ -304,90 +316,100 @@ int main(int argc, char **argv) //10, 10, 1.505, 0.109, 10, 10, 1.72, 0.135, 10, 10, 1.300, 0.350, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}// start over with pdg ; - minuit->SetParameter(0, "f_2 strength Re", vals[0], 1, 0, 0); + minuit->SetVariable(0, "f_2 strength Re", vals[0], 1); //minuit->FixParameter(0); - minuit->SetParameter(1, "f_2 strength Im", vals[1], 1, 0, 0); + minuit->SetVariable(1, "f_2 strength Im", vals[1], 1); //minuit->FixParameter(1); - minuit->SetParameter(2, "f_2 mass", vals[2], 1, 1.20, 1.35); + minuit->SetVariable(2, "f_2 mass", vals[2], 1); + minuit->SetVariableLimits(2, 1.2, 1.35); //minuit->FixParameter(2); - minuit->SetParameter(3, "f_2 width", vals[3], 1, 0.1, 0.5); + minuit->SetVariable(3, "f_2 width", vals[3], 1); + minuit->SetVariableLimits(3, 0.1, 0.5); //minuit->FixParameter(3); - minuit->SetParameter(4, "f_2' strength Re", vals[4], 1, 0, 0); + minuit->SetVariable(4, "f_2' strength Re", vals[4], 1); //minuit->FixParameter(4); - minuit->SetParameter(5, "f_2' strength Im", vals[5], 1, 0, 0); + minuit->SetVariable(5, "f_2' strength Im", vals[5], 1); //minuit->FixParameter(5); - minuit->SetParameter(6, "f_2' mass", vals[6], 1, 1.45, 1.6); + minuit->SetVariable(6, "f_2' mass", vals[6], 1); + minuit->SetVariableLimits(6, 1.45, 1.6); //minuit->FixParameter(6); - minuit->SetParameter(7, "f_2' width", vals[7], 1, 0.05, 0.5); + minuit->SetVariable(7, "f_2' width", vals[7], 1); + minuit->SetVariableLimits(7, 0.05, 0.5); //minuit->FixParameter(7); - minuit->SetParameter(8, "D BG strength Re", vals[8], 1, 0, 0); + minuit->SetVariable(8, "D BG strength Re", vals[8], 1); //minuit->FixParameter(8); - minuit->SetParameter(9, "D BG strength Im", vals[9], 1, 0, 0); + minuit->SetVariable(9, "D BG strength Im", vals[9], 1); ///minuit->FixParameter(9); - //minuit->SetParameter(10, "D BG const", vals[10], 1, 0, 0); + //minuit->SetVariable(10, "D BG const", vals[10], 1); //minuit->FixParameter(10); - minuit->SetParameter(10, "D BG linear", vals[10], .1, 0, 0); + minuit->SetVariable(10, "D BG linear", vals[10], .1); //minuit->FixParameter(10); - minuit->SetParameter(11, "D BG quadratic", vals[11], .1, 0, 0); + minuit->SetVariable(11, "D BG quadratic", vals[11], .1); //minuit->FixParameter(11); /* - minuit->SetParameter(12, "f_0(980) strength Re", vals[12], 1, 0, 0); + minuit->SetVariable(12, "f_0(980) strength Re", vals[12], 1); //minuit->FixParameter(12); - minuit->SetParameter(13, "f_0(980) strength Im", vals[13], 1, 0, 0); + minuit->SetVariable(13, "f_0(980) strength Im", vals[13], 1); //minuit->FixParameter(13); - minuit->SetParameter(14, "f_0(980) mass", vals[14], 0.01, 0, 0); + minuit->SetVariable(14, "f_0(980) mass", vals[14], 0.01); minuit->FixParameter(14); - minuit->SetParameter(15, "f_0(980) width", vals[15], 0.01, 0, 0); + minuit->SetVariable(15, "f_0(980) width", vals[15], 0.01); minuit->FixParameter(15); */ - //minuit->SetParameter(12, "S BG strength Re", vals[12], 1, 0, 0); + //minuit->SetVariable(12, "S BG strength Re", vals[12], 1); //minuit->FixParameter(12); - minuit->SetParameter(12, "S BG strength Re", vals[12], 1, 0, 0); + minuit->SetVariable(12, "S BG strength Re", vals[12], 1); //minuit->FixParameter(12); - //minuit->SetParameter(12, "S BG const", vals[12], 1, 0, 0); + //minuit->SetVariable(12, "S BG const", vals[12], 1); //minuit->FixParameter(12); - minuit->SetParameter(13, "S BG linear", vals[13], .1, 0, 0); + minuit->SetVariable(13, "S BG linear", vals[13], .1); //minuit->FixParameter(13); - minuit->SetParameter(14, "S BG quadratic", vals[14], .1, 0, 0); + minuit->SetVariable(14, "S BG quadratic", vals[14], .1); //minuit->FixParameter(14); - minuit->SetParameter(15, "f_0(1500) strength Re", vals[15], 1, 0, 0); + minuit->SetVariable(15, "f_0(1500) strength Re", vals[15], 1); //minuit->FixParameter(15); - minuit->SetParameter(16, "f_0(1500) strength Im", vals[16], 1, 0, 0); + minuit->SetVariable(16, "f_0(1500) strength Im", vals[16], 1); //minuit->FixParameter(16); - minuit->SetParameter(17, "f_0(1500) mass", vals[17], 1, 1.45, 1.55); + minuit->SetVariable(17, "f_0(1500) mass", vals[17], 1); + minuit->SetVariableLimits(17, 1.45, 1.55); //minuit->FixParameter(17); - minuit->SetParameter(18, "f_0(1500) width", vals[18], 1, 0.05, 0.15); + minuit->SetVariable(18, "f_0(1500) width", vals[18], 1); + minuit->SetVariableLimits(18, 0.05, 0.15); //minuit->FixParameter(18); - minuit->SetParameter(19, "f_0(1710) strength Re", vals[19], 1, 0, 0); + minuit->SetVariable(19, "f_0(1710) strength Re", vals[19], 1); //minuit->FixParameter(19); - minuit->SetParameter(20, "f_0(1710) strength Im", vals[20], 1, 0, 0); + minuit->SetVariable(20, "f_0(1710) strength Im", vals[20], 1); //minuit->FixParameter(20); - minuit->SetParameter(21, "f_0(1710) mass", vals[21], 1, 1.65, 1.75); + minuit->SetVariable(21, "f_0(1710) mass", vals[21], 1); + minuit->SetVariableLimits(21, 1.65, 1.75); //minuit->FixParameter(21); - minuit->SetParameter(22, "f_0(1710) width", vals[22], 1, 0.05, 0.2); + minuit->SetVariable(22, "f_0(1710) width", vals[22], 1); + minuit->SetVariableLimits(22, 0.05, 0.2); //minuit->FixParameter(22); - minuit->SetParameter(23, "f_0(1370) strength Re", vals[23], 1, 0, 0); + minuit->SetVariable(23, "f_0(1370) strength Re", vals[23], 1); //minuit->FixParameter(23); - minuit->SetParameter(24, "f_0(1370) strength Im", vals[24], 1, 0, 0); + minuit->SetVariable(24, "f_0(1370) strength Im", vals[24], 1); //minuit->FixParameter(24); - minuit->SetParameter(25, "f_0(1370) mass", vals[25], 1, 1.25, 1.40); + minuit->SetVariable(25, "f_0(1370) mass", vals[25], 1); + minuit->SetVariableLimits(25, 1.25, 1.4); //minuit->FixParameter(25); - minuit->SetParameter(26, "f_0(1370) width", vals[26], 1, 0.05, 0.40); + minuit->SetVariable(26, "f_0(1370) width", vals[26], 1); + minuit->SetVariableLimits(26, 0.05, 0.40); //minuit->FixParameter(26); /* - minuit->SetParameter(24, "f_0(1370) strength Re", vals[24], 1, 0, 0); + minuit->SetVariable(24, "f_0(1370) strength Re", vals[24], 1); //minuit->FixParameter(24); - minuit->SetParameter(25, "f_0(1370) strength Im", vals[25], 1, 0, 0); + minuit->SetVariable(25, "f_0(1370) strength Im", vals[25], 1); //minuit->FixParameter(25); - minuit->SetParameter(26, "f_0(1370) mass", vals[26], 0.01, 0, 0); + minuit->SetVariable(26, "f_0(1370) mass", vals[26], 0.01); minuit->FixParameter(26); - minuit->SetParameter(27, "f_0(1370) width", vals[27], 0.01, 0, 0); + minuit->SetVariable(27, "f_0(1370) width", vals[27], 0.01); minuit->FixParameter(27); */ @@ -398,20 +420,19 @@ int main(int argc, char **argv) TStopwatch sw; sw.Start(); - minuit->CreateMinimizer(); int iret = minuit->Minimize(); sw.Stop(); cout << "iret = " << iret << " after " << sw.CpuTime() << " s." << endl; if (iret) - minuit->PrintResults(1, 0.); + minuit->PrintResults(); //if (!iret) { fitModel* model = fitModel::getFitModelForName(info->getModelName()); vector x; - for (int i = 0; i < minuit->GetNumberTotalParameters(); i++) - x.push_back(minuit->GetParameter(i)); + for (size_t i = 0; i < info->getNvars(); i++) + x.push_back(minuit->X()[i]); for (size_t i = 0; i < info->getNbins(); i++) { @@ -428,7 +449,7 @@ int main(int argc, char **argv) complex SwaveBG = model->valueForWave("S0BG"); complex f2 = model->valueForWave("f2"); - complex f2prime = model->valueForWave("f2prime"); + //complex f2prime = model->valueForWave("f2prime"); complex f01370 = model->valueForWave("f01370"); complex f01500 = model->valueForWave("f01500"); @@ -454,7 +475,7 @@ int main(int argc, char **argv) gHist.getHist("hPhaseSD", "phase D - S", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, arg(Swave/Dwave)); gHist.getHist("hf2", "f2(1270)", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, norm(f2)); - gHist.getHist("hf2prime", "f2'(1525)", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, norm(f2prime)); + //gHist.getHist("hf2prime", "f2'(1525)", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, norm(f2prime)); gHist.getHist("hf01370", "f0 1370", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, norm(f01370)); gHist.getHist("hf01500", "f01500", info->getNbins(), info->getLower(), info->getUpper())->SetBinContent(i, norm(f01500)); diff --git a/myFit.cc b/myFit.cc index c81ec7b..c3a6166 100644 --- a/myFit.cc +++ b/myFit.cc @@ -12,7 +12,6 @@ using namespace std; #include "TH3.h" #include "TCanvas.h" #include "TRandom1.h" -#include "TFitterMinuit.h" #include "TStopwatch.h" #include "TDirectory.h" #include "TFile.h" @@ -51,14 +50,14 @@ class combinedLikelihood : public ROOT::Minuit2::FCNBase { { } - ~combinedLikelihood() { for (size_t i = 0; i < myLs.size(); i++) delete myLs[i]; } + //~combinedLikelihood() { for (size_t i = 0; i < myLs.size(); i++) delete myLs[i]; } void addChannel(vector& RDevents, vector& MCevents, vector& MCallEvents) { - size_t idxBranching = /*2*NWAVES*/ 16 + this->getNChannels(); + size_t idxBranching = 2*ws.getNwaves() + this->getNChannels(); myLs.push_back(new likelihood(ws, RDevents, MCevents, MCallEvents, nBins, threshold, binWidth, idxBranching)); } @@ -127,6 +126,11 @@ class combinedLikelihood : public ROOT::Minuit2::FCNBase { return -result; } + double operator()( const double * x) const { + std::vector p(x , x + 2*ws.getNwaves() + this->getNChannels()); + return (*this)(p); + } + size_t getNChannels() const { return myLs.size(); } size_t MCeventsInBin() const { size_t sum = 0; for (size_t i = 0; i < myLs.size(); i++) sum += myLs[i]->MCeventsInBin(); return sum; } }; @@ -388,7 +392,7 @@ myFit() double upper = threshold + nBins*binWidth; vector positive; - //positive.push_back(wave("P+", 1, 1, nBins, lower, upper, true)); + positive.push_back(wave("P+", 1, 1, nBins, lower, upper, true)); positive.push_back(wave("D+", 2, 1, nBins, lower, upper)); //positive.push_back(wave("F+", 3, 1, nBins, lower, upper)); //positive.push_back(wave("G+", 4, 1, 78, 1.7, upper)); @@ -396,11 +400,11 @@ myFit() vector negative; negative.push_back(wave("S0", 0, 0, nBins, lower, upper, true)); - //negative.push_back(wave("P0", 1, 0, nBins, lower, upper)); - //negative.push_back(wave("P-", 1, 1, nBins, lower, upper)); + negative.push_back(wave("P0", 1, 0, nBins, lower, upper)); + negative.push_back(wave("P-", 1, 1, nBins, lower, upper)); negative.push_back(wave("D0", 2, 0, nBins, lower, upper)); negative.push_back(wave("D-", 2, 1, nBins, lower, upper)); - // negative.push_back(wave("G0", 4, 0, nBins, lower, upper)); + //negative.push_back(wave("G0", 4, 0, nBins, lower, upper)); //negative.push_back(wave("G0", 4, 0, 78, 1.7, upper)); //negative.push_back(wave("G-", 4, 1, 78, 1.7, upper)); @@ -535,21 +539,16 @@ myFit() float tPr; float theta; float phi; - float likeK, likePi; tree->SetBranchAddress("mKpi", &mX); tree->SetBranchAddress("tPrime", &tPr); tree->SetBranchAddress("theta", &theta); tree->SetBranchAddress("phi", &phi); - //tree->SetBranchAddress("likeK", &likeK); - //tree->SetBranchAddress("likePi", &likePi); RDevents.reserve(tree->GetEntries()); for (Long_t i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); - //if (likeK != -1 && likePi > 2*likeK) - //continue; event e(mX, tPr, theta, phi); RDevents.push_back(e); fillRDhists(e); @@ -807,16 +806,25 @@ myFit() TStopwatch sw; const size_t nParams = lastIdx + myL.getNChannels(); - TFitterMinuit* minuit = new TFitterMinuit(); - minuit->SetMinuitFCN(&myL); + ROOT::Math::Minimizer* minuit = ROOT::Math::Factory::CreateMinimizer("Minuit2", ""); + ROOT::Math::Functor func(myL, 2*ws.getNwaves() + myL.getNChannels()); + minuit->SetFunction(func); - TFitterMinuit* fitamb = new TFitterMinuit(); - fitamb->SetMinuitFCN(&myL); - - // Do not show fit results: + // set tolerance , etc... + minuit->SetErrorDef(0.5); //statistical scale for negative log likelihood + // minuit->SetMaxFunctionCalls(1000000); + // minuit->SetMaxIterations(10000); + minuit->SetTolerance(0.1); minuit->SetPrintLevel(0); - fitamb->SetPrintLevel(0); + ROOT::Math::Minimizer* fitamb = ROOT::Math::Factory::CreateMinimizer("Minuit2", ""); + fitamb->SetFunction(func); + + // set tolerance , etc... + fitamb->SetPrintLevel(0); + fitamb->SetErrorDef(0.5); //statistical scale for negative log likelihood + fitamb->SetTolerance(0.1); + //TH3* hPredict = new TH3D("hPredict", "prediction", nBins, 0, nBins, 100, -1, 1, 100, -M_PI, M_PI); TH2* hthpre = new TH2D("hthpre", "prediction for cos(#theta)", nBins, threshold, threshold + nBins*binWidth, 100, -1, 1); @@ -880,29 +888,28 @@ myFit() if (!startingValues[j].fixed) { //if (j != 0) - minuit->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], 10/*vStartingValues[j]*0.01*/, 0, 0); + minuit->SetVariable(j, (startingValues[j].name).c_str(), + vStartingValues[j], 10); } else { - minuit->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], 1, 0, 0); - minuit->FixParameter(j); + minuit->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], 1); + minuit->FixVariable(j); } } for (size_t j = nParams - myL.getNChannels(); j < nParams; j++) { - minuit->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], .1, 0, 1); + minuit->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], .1); if (startingValues[j].fixed) - minuit->FixParameter(j); + minuit->FixVariable(j); } // Run minimizer. - minuit->CreateMinimizer(); int iret = minuit->Minimize(); - - while ( iret != 0 ){ + + while ( iret != 1 ){ // if fit did not converge, use random starting values until it does for (size_t iSV = 0; iSV < nParams; iSV++) { @@ -920,37 +927,41 @@ myFit() if (!startingValues[j].fixed) { //if (j != 0) - minuit->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], 10/*vStartingValues[j]*0.01*/, 0, 0); + minuit->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], 10); } else { - minuit->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], 1, 0, 0); - minuit->FixParameter(j); + minuit->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], 1); + minuit->FixVariable(j); } } iret = minuit->Minimize(); } + + minuit->PrintResults(); + sw.Stop(); - // cout << "iret = " << iret << " after " << sw.CpuTime() << " s." << endl; + cout << "Fit converged after " << sw.CpuTime() << " s." << endl; + - if (iret == 0) + if (iret == 1) { // vector vStartingValue(nParams); for (size_t j = 0; j < nParams; j++) { - vStartingValues[j] = minuit->GetParameter(j); + vStartingValues[j] = minuit->X()[j]; } for (size_t i = 0; i < ws.size(); i++) for (size_t j = 0; j < ws[i].waves.size(); j++) { size_t idx = ws[i].waves[j].getIndex(); - values[idx] = minuit->GetParameter(idx); - values[idx+1] = minuit->GetParameter(idx+1); + values[idx] = minuit->X()[idx]; + values[idx+1] = minuit->X()[idx+1]; } if (!ambiguous){ @@ -959,44 +970,42 @@ myFit() { const wave& w1 = ws[i1].waves[j1]; size_t idx1 = w1.getIndex(); - size_t idx1Cov = w1.idxInCovariance(minuit); for (size_t i2 = 0; i2 < ws.size(); i2++) for (size_t j2 = 0; j2 < ws[i2].waves.size(); j2++) { const wave& w2 = ws[i2].waves[j2]; size_t idx2 = w2.getIndex(); - size_t idx2Cov = w2.idxInCovariance(minuit); // Four possibilities: // 1) both free -> simply copy values to covMat if (!w1.phaseLocked && !w2.phaseLocked) { - (*covMat)(idx1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov); - (*covMat)(idx1,idx2+1) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov+1); - (*covMat)(idx1+1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov); - (*covMat)(idx1+1,idx2+1) = minuit->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov+1); + (*covMat)(idx1,idx2) = minuit->CovMatrix(idx1, idx2); + (*covMat)(idx1,idx2+1) = minuit->CovMatrix(idx1, idx2+1); + (*covMat)(idx1+1,idx2) = minuit->CovMatrix(idx1+1, idx2); + (*covMat)(idx1+1,idx2+1) = minuit->CovMatrix(idx1+1, idx2+1); } // 2) first free, second fixed -> covariances with Im(w1) are zero else if (!w1.phaseLocked && w2.phaseLocked) { - (*covMat)(idx1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov); + (*covMat)(idx1,idx2) = minuit->CovMatrix(idx1, idx2); (*covMat)(idx1,idx2+1) = 0; - (*covMat)(idx1+1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov); + (*covMat)(idx1+1,idx2) = minuit->CovMatrix(idx1+1, idx2); (*covMat)(idx1+1,idx2+1) = 0; } // 3) vice versa else if (w1.phaseLocked && !w2.phaseLocked) { - (*covMat)(idx1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov); - (*covMat)(idx1,idx2+1) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov+1); + (*covMat)(idx1,idx2) = minuit->CovMatrix(idx1, idx2); + (*covMat)(idx1,idx2+1) = minuit->CovMatrix(idx1, idx2+1); (*covMat)(idx1+1,idx2) = 0; (*covMat)(idx1+1,idx2+1) = 0; } // 4) both fixed (impossible) else { - (*covMat)(idx1,idx2) = minuit->GetCovarianceMatrixElement(idx1Cov, idx2Cov); + (*covMat)(idx1,idx2) = minuit->CovMatrix(idx1, idx2); (*covMat)(idx1,idx2+1) = 0; (*covMat)(idx1+1,idx2) = 0; (*covMat)(idx1+1,idx2+1) = 0; @@ -1007,9 +1016,9 @@ myFit() NMCevents = myL.MCeventsInBin(); outTree->Fill(); } - - for (int j = 0; j < minuit->GetNumberTotalParameters(); j++) - startingValues[j].value = minuit->GetParameter(j); + + for (size_t iPar = 0; iPar < nParams; iPar++) + startingValues[iPar].value = minuit->X()[iPar]; gHist.getHist("hMClikelihood", "MC likelihood of result", nBins, lower, upper) ->SetBinContent(iBin+1, myL.calc_mc_likelihood(vStartingValues)); @@ -1032,14 +1041,14 @@ myFit() if (myL.getNChannels() == 2) { - hBR->SetBinContent(iBin+1, minuit->GetParameter(nParams - 1)); - hBR->SetBinError(iBin+1, minuit->GetParError(nParams - 1)); + hBR->SetBinContent(iBin+1, minuit->X()[nParams - 1]); + hBR->SetBinError(iBin+1, minuit->Errors()[nParams - 1]); } vector result; - for (int iPar = 0; iPar < minuit->GetNumberTotalParameters(); iPar++) + for (size_t iPar = 0; iPar < nParams; iPar++) { - result.push_back(minuit->GetParameter(iPar)); + result.push_back(minuit->X()[iPar]); } @@ -1080,7 +1089,7 @@ myFit() { // calculate coefficients vector amplitudes; - for (int k = 0; k < lastIdx; k++) + for (unsigned int k = 0; k < lastIdx; k++) amplitudes.push_back(values[k]); vector > input; @@ -1156,26 +1165,26 @@ myFit() { if (!startingValues[j].fixed /*|| j < 4*/) { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - ambiguous[n][j], 10/*ambiguous[n][j]*0.01*/, 0, 0); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + ambiguous[n][j], 10/*ambiguous[n][j]*0.01*/); } else { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - ambiguous[n][j], 1, 0, 0); - fitamb->FixParameter(j); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + ambiguous[n][j], 1); + fitamb->FixVariable(j); } } for (size_t j = nParams - myL.getNChannels(); j < nParams; j++) { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], .1, 0, 1); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], .1); if (startingValues[j].fixed) - fitamb->FixParameter(j); + fitamb->FixVariable(j); } // Run minimizer. - fitamb->CreateMinimizer(); + //fitamb->CreateMinimizer(); int amret = fitamb->Minimize(); // if fit did not converge, use new random starting values until it does @@ -1191,22 +1200,22 @@ myFit() { if (!startingValues[j].fixed /*|| j < 4*/) { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - ambiguous[n][j], 10/*ambiguous[n][j]*0.01*/, 0, 0); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + ambiguous[n][j], 10/*ambiguous[n][j]*0.01*/); } else { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - ambiguous[n][j], 1, 0, 0); - fitamb->FixParameter(j); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + ambiguous[n][j], 1); + fitamb->FixVariable(j); } } for (size_t j = nParams - myL.getNChannels(); j < nParams; j++) { - fitamb->SetParameter(j, startingValues[j].name.c_str(), - vStartingValues[j], .1, 0, 1); + fitamb->SetVariable(j, startingValues[j].name.c_str(), + vStartingValues[j], .1); if (startingValues[j].fixed) - fitamb->FixParameter(j); + fitamb->FixVariable(j); } amret = fitamb->Minimize(); } @@ -1218,15 +1227,15 @@ myFit() for (int i = 0; i < 4; i++) { - positiveStart[i] = fitamb->GetParameter(i); + positiveStart[i] = fitamb->X()[i]; } for (size_t i = 0; i < ws.size(); i++) for (size_t j = 0; j < ws[i].waves.size(); j++) { size_t idx = ws[i].waves[j].getIndex(); - values[idx] = fitamb->GetParameter(idx); - values[idx+1] = fitamb->GetParameter(idx+1); + values[idx] = fitamb->X()[idx]; + values[idx+1] = fitamb->X()[idx+1]; } for (size_t i1 = 0; i1 < ws.size(); i1++) @@ -1234,44 +1243,42 @@ myFit() { const wave& w1 = ws[i1].waves[j1]; size_t idx1 = w1.getIndex(); - size_t idx1Cov = w1.idxInCovariance(fitamb); for (size_t i2 = 0; i2 < ws.size(); i2++) for (size_t j2 = 0; j2 < ws[i2].waves.size(); j2++) { const wave& w2 = ws[i2].waves[j2]; size_t idx2 = w2.getIndex(); - size_t idx2Cov = w2.idxInCovariance(fitamb); // Four possibilities: // 1) both free -> simply copy values to covMat if (!w1.phaseLocked && !w2.phaseLocked) { - (*covMat)(idx1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov); - (*covMat)(idx1,idx2+1) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov+1); - (*covMat)(idx1+1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov); - (*covMat)(idx1+1,idx2+1) = fitamb->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov+1); + (*covMat)(idx1,idx2) = fitamb->CovMatrix(idx1, idx2); + (*covMat)(idx1,idx2+1) = fitamb->CovMatrix(idx1, idx2+1); + (*covMat)(idx1+1,idx2) = fitamb->CovMatrix(idx1+1, idx2); + (*covMat)(idx1+1,idx2+1) = fitamb->CovMatrix(idx1+1, idx2+1); } // 2) first free, second fixed -> covariances with Im(w1) are zero else if (!w1.phaseLocked && w2.phaseLocked) { - (*covMat)(idx1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov); + (*covMat)(idx1,idx2) = fitamb->CovMatrix(idx1, idx2); (*covMat)(idx1,idx2+1) = 0; - (*covMat)(idx1+1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov+1, idx2Cov); + (*covMat)(idx1+1,idx2) = fitamb->CovMatrix(idx1+1, idx2); (*covMat)(idx1+1,idx2+1) = 0; } // 3) vice versa else if (w1.phaseLocked && !w2.phaseLocked) { - (*covMat)(idx1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov); - (*covMat)(idx1,idx2+1) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov+1); + (*covMat)(idx1,idx2) = fitamb->CovMatrix(idx1, idx2); + (*covMat)(idx1,idx2+1) = fitamb->CovMatrix(idx1, idx2+1); (*covMat)(idx1+1,idx2) = 0; (*covMat)(idx1+1,idx2+1) = 0; } // 4) both fixed (impossible) else { - (*covMat)(idx1,idx2) = fitamb->GetCovarianceMatrixElement(idx1Cov, idx2Cov); + (*covMat)(idx1,idx2) = fitamb->CovMatrix(idx1, idx2); (*covMat)(idx1,idx2+1) = 0; (*covMat)(idx1+1,idx2) = 0; (*covMat)(idx1+1,idx2+1) = 0; @@ -1320,7 +1327,7 @@ myFit() fulltime.Stop(); cout << "Took " << fulltime.CpuTime() << " s CPU time, " << fulltime.RealTime() << " s wall time." << endl; - cout << "In " << failBins << " bins TFitterMinuit::Minimization DID not converge !" << endl; + cout << "In " << failBins << " bins the minimization DID not converge !" << endl; } diff --git a/wave.cc b/wave.cc index 2c1d034..7dc6eb4 100644 --- a/wave.cc +++ b/wave.cc @@ -1,7 +1,7 @@ #include #include -#include "TFitterMinuit.h" +#include "Math/Minimizer.h" #include "TMatrixDSym.h" #include "TVectorD.h" @@ -48,40 +48,40 @@ wave::getHistPhase(const wave& other) // the correct indices into the covariance matrix. Therefore this // contortion ... which carries the assumption also made above that we // never fix real parts. -size_t -wave::idxInCovariance(const TFitterMinuit* minuit) const -{ - size_t countFixedBelow = 0; - for (size_t i = 0; i < idx; i++) - countFixedBelow += minuit->IsFixed(i); - return idx - countFixedBelow; -} +// size_t +// wave::idxInCovariance(const ROOT::Math::Minimizer* minuit) const +// { +// size_t countFixedBelow = 0; +// for (size_t i = 0; i < idx; i++) +// countFixedBelow += minuit->IsFixedVariable(i); +// return idx - countFixedBelow; +// } void -wave::fillHistIntensity(int iBin, const TFitterMinuit* minuit) +wave::fillHistIntensity(int iBin, const ROOT::Math::Minimizer* minuit) { - complex a(minuit->GetParameter(idx), minuit->GetParameter(idx+1)); + complex a(minuit->X()[idx], minuit->X()[idx+1]); histIntensity->SetBinContent(iBin+1,norm(a)); double error; - if (minuit->IsFixed(idx+1)) - error = 2*abs(a)*minuit->GetParError(idx); + if (minuit->IsFixedVariable(idx+1)) + error = 2*abs(a)*minuit->Errors()[idx]; else { - size_t idxCov = idxInCovariance(minuit); - double cov = minuit->GetCovarianceMatrixElement(idxCov, idxCov + 1); + //size_t idxCov = idxInCovariance(minuit); + double cov = minuit->CovMatrix(idx, idx + 1); - error = 2*(sqrt(pow(real(a) * minuit->GetParError(idx), 2) - + pow(imag(a) * minuit->GetParError(idx+1), 2) - + (2*real(a)*imag(a)*cov))); + error = 2*(sqrt(pow(real(a) * minuit->Errors()[idx], 2) + + pow(imag(a) * minuit->Errors()[idx+1], 2) + + (2*real(a)*imag(a)*cov))); } histIntensity->SetBinError(iBin+1, error); } void -wave::fillHistPhase(int iBin, const wave& other, const TFitterMinuit* minuit) +wave::fillHistPhase(int iBin, const wave& other, const ROOT::Math::Minimizer* minuit) { TH1* h = 0; if (mHistPhase.find(other.name) == mHistPhase.end()) @@ -100,8 +100,8 @@ wave::fillHistPhase(int iBin, const wave& other, const TFitterMinuit* minuit) h = mHistPhase[other.name]; - complex a1(minuit->GetParameter(idx), minuit->GetParameter(idx + 1)); - complex a2(minuit->GetParameter(other.getIndex()), minuit->GetParameter(other.getIndex()+1)); + complex a1(minuit->X()[idx], minuit->X()[idx + 1]); + complex a2(minuit->X()[other.getIndex()], minuit->X()[other.getIndex()+1]); double phi = arg(a1 / a2); double oldPhase = 0; @@ -128,47 +128,33 @@ wave::fillHistPhase(int iBin, const wave& other, const TFitterMinuit* minuit) // 2. Im(a2) fixed. // 3. Neither fixed. - if (minuit->IsFixed(idx + 1) || minuit->IsFixed(other.getIndex() + 1)) + if (minuit->IsFixedVariable(idx + 1) || minuit->IsFixedVariable(other.getIndex() + 1)) { - if (minuit->IsFixed(idx + 1) && minuit->IsFixed(other.getIndex() + 1)) + if (minuit->IsFixedVariable(idx + 1) && minuit->IsFixedVariable(other.getIndex() + 1)) { // Can't happen, ignore. h->SetBinError(iBin+1, 0.2); } else { - bool swapped = false; size_t idxFix = idx; size_t idxFree = other.getIndex(); - if (minuit->IsFixed(other.getIndex() + 1)) - { - swapped = true; - std::swap(idxFix, idxFree); - } + if (minuit->IsFixedVariable(other.getIndex() + 1)) + std::swap(idxFix, idxFree); + TMatrixDSym cov(3); - size_t idxCovFix, idxCovFree; - if (swapped) - { - idxCovFix = other.idxInCovariance(minuit); - idxCovFree = idxInCovariance(minuit); - } - else - { - idxCovFix = idxInCovariance(minuit); - idxCovFree = other.idxInCovariance(minuit); - } - cov(0,0) = minuit->GetCovarianceMatrixElement(idxCovFix, idxCovFix); - cov(1,1) = minuit->GetCovarianceMatrixElement(idxCovFree, idxCovFree); - cov(2,2) = minuit->GetCovarianceMatrixElement(idxCovFree + 1, idxCovFree + 1); + cov(0,0) = minuit->CovMatrix(idxFix, idxFix); + cov(1,1) = minuit->CovMatrix(idxFree, idxFree); + cov(2,2) = minuit->CovMatrix(idxFree + 1, idxFree + 1); - cov(0,1) = cov(1,0) = minuit->GetCovarianceMatrixElement(idxCovFix, idxCovFree); - cov(0,2) = cov(2,0) = minuit->GetCovarianceMatrixElement(idxCovFix, idxCovFree + 1); - cov(1,2) = cov(2,1) = minuit->GetCovarianceMatrixElement(idxCovFree, idxCovFree + 1); + cov(0,1) = cov(1,0) = minuit->CovMatrix(idxFix, idxFree); + cov(0,2) = cov(2,0) = minuit->CovMatrix(idxFix, idxFree + 1); + cov(1,2) = cov(2,1) = minuit->CovMatrix(idxFree, idxFree + 1); - double values[3] = { minuit->GetParameter(idxFix), - minuit->GetParameter(idxFree), - minuit->GetParameter(idxFree + 1) }; + double values[3] = { minuit->X()[idxFix], + minuit->X()[idxFree], + minuit->X()[idxFree + 1] }; TVectorD gradient(3); // The algorithm follows the one in TF1::Derivative() : // df(x) = (4 D(h/2) - D(h)) / 3 @@ -222,23 +208,24 @@ wave::fillHistPhase(int iBin, const wave& other, const TFitterMinuit* minuit) else { TMatrixDSym cov(4); - size_t idxCov, idxCovOther; - idxCov = idxInCovariance(minuit); - idxCovOther = other.idxInCovariance(minuit); - cov(0,0) = minuit->GetCovarianceMatrixElement(idxCov, idxCov); - cov(1,1) = minuit->GetCovarianceMatrixElement(idxCov + 1, idxCov + 1); - cov(2,2) = minuit->GetCovarianceMatrixElement(idxCovOther, idxCovOther); - cov(3,3) = minuit->GetCovarianceMatrixElement(idxCovOther + 1, idxCovOther + 1); - - cov(0,1) = cov(1,0) = minuit->GetCovarianceMatrixElement(idxCov, idxCov + 1); - cov(0,2) = cov(2,0) = minuit->GetCovarianceMatrixElement(idxCov, idxCovOther); - cov(0,3) = cov(3,0) = minuit->GetCovarianceMatrixElement(idxCov, idxCovOther + 1); - cov(1,2) = cov(2,1) = minuit->GetCovarianceMatrixElement(idxCov + 1, idxCovOther); - cov(1,3) = cov(3,1) = minuit->GetCovarianceMatrixElement(idxCov + 1, idxCovOther + 1); - cov(2,3) = cov(3,2) = minuit->GetCovarianceMatrixElement(idxCovOther, idxCovOther + 1); - - double values[4] = { minuit->GetParameter(idx), minuit->GetParameter(idx + 1), - minuit->GetParameter(other.getIndex()), minuit->GetParameter(other.getIndex()+1) }; + size_t idxOther = other.getIndex(); + // size_t idxCov, idxCovOther; + // idxCov = idxInCovariance(minuit); + // idxCovOther = other.idxInCovariance(minuit); + cov(0,0) = minuit->CovMatrix(idx, idx); + cov(1,1) = minuit->CovMatrix(idx + 1, idx + 1); + cov(2,2) = minuit->CovMatrix(idxOther, idxOther); + cov(3,3) = minuit->CovMatrix(idxOther + 1, idxOther + 1); + + cov(0,1) = cov(1,0) = minuit->CovMatrix(idx, idx + 1); + cov(0,2) = cov(2,0) = minuit->CovMatrix(idx, idxOther); + cov(0,3) = cov(3,0) = minuit->CovMatrix(idx, idxOther + 1); + cov(1,2) = cov(2,1) = minuit->CovMatrix(idx + 1, idxOther); + cov(1,3) = cov(3,1) = minuit->CovMatrix(idx + 1, idxOther + 1); + cov(2,3) = cov(3,2) = minuit->CovMatrix(idxOther, idxOther + 1); + + double values[4] = { minuit->X()[idx], minuit->X()[idx + 1], + minuit->X()[other.getIndex()], minuit->X()[other.getIndex()+1] }; TVectorD gradient(4); // The algorithm follows the one in TF1::Derivative() : // df(x) = (4 D(h/2) - D(h)) / 3 diff --git a/wave.h b/wave.h index 835b343..49cfc05 100644 --- a/wave.h +++ b/wave.h @@ -7,7 +7,7 @@ #include #include "TH1.h" -#include "TFitterMinuit.h" +#include "Math/Minimizer.h" using namespace std; @@ -35,7 +35,7 @@ struct wave { // Thanks to the brilliance of Minuit2, the covariance matrix doesn't // contain the fixed parameters, making the mapping of fit parameters // to covariance matrix elements awkward. This should help. - size_t idxInCovariance(const TFitterMinuit* minuit) const; + size_t idxInCovariance(const ROOT::Math::Minimizer* minuit) const; const string& getName() const { return name; } size_t getL() const { return l; } @@ -43,8 +43,8 @@ struct wave { void buildHists(int nBins, double lower, double upper); TH1* getHistIntensity() const { return histIntensity; } TH1* getHistPhase(const wave& other); - void fillHistIntensity(int iBin, const TFitterMinuit* minuit); - void fillHistPhase(int iBin, const wave& other, const TFitterMinuit* minuit); + void fillHistIntensity(int iBin, const ROOT::Math::Minimizer* minuit); + void fillHistPhase(int iBin, const wave& other, const ROOT::Math::Minimizer* minuit); ClassDefNV(wave, 1) };