Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2310) +++ /trunk/MagicSoft/Mars/Changelog (revision 2311) @@ -4,5 +4,6 @@ * manalysis/MCT1FindSupercuts.[h,cc], manalysis/MCT1Supercuts.[h,cc], - manalysis/MCT1SupercutsCalc.[h,cc], manalysis/MMinuitInterface.[h,cc]: + manalysis/MCT1SupercutsCalc.[h,cc], manalysis/MMinuitInterface.[h,cc], + mhist/MHFindSignificance.[h,cc]: - changed some variables and member functions with respect to an upcoming Minimization Class @@ -11,4 +12,7 @@ - added some sanity checks - simplified some variable usage + + * mhist/MHCT1Supercuts.[h,cc]: + - removed obsolete SetupFill Index: /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2310) +++ /trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2311) @@ -31,5 +31,5 @@ Bool_t SetParameters(const TArrayD &d); - const TArrayD &GetParameters() const; + const TArrayD &GetParameters() const { return fParameters; } const Double_t *GetLengthUp() const { return fLengthUp; } Index: /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2310) +++ /trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2311) @@ -61,4 +61,5 @@ #include +#include #include #include @@ -849,182 +850,162 @@ while (1) { - - //-------------------------------------------------- - // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...) - - TString funcname = "Poly"; - Double_t xmin = 0.0; - Double_t xmax = 90.0; - - TString formula = "[0]"; - TString bra1 = "+["; - TString bra2 = "]"; - TString xpower = "*x"; - TString newpower = "*x"; - for (Int_t i=1; i<=fDegree; i++) - { - formula += bra1; - formula += i; - formula += bra2; - formula += xpower; - - xpower += newpower; - } - - //*fLog << "FitPolynomial : formula = " << formula << endl; - - fPoly = new TF1(funcname, formula, xmin, xmax); - TList *funclist = fHist->GetListOfFunctions(); - funclist->Add(fPoly); - - //------------------------ - // attention : the dimensions must agree with those in CallMinuit() - - char parname[80][100]; - Double_t vinit[80]; - Double_t step[80]; - Double_t limlo[80]; - Double_t limup[80]; - Int_t fix[80]; - - UInt_t npar = fDegree+1; - - for (UInt_t j=0; jGetEntries(); - - - // use the subsequernt loop if you want to apply the - // constraint : uneven derivatives (at alpha=0) = zero - for (UInt_t j=1; jmnstat(fmin, fedm, errdef, npari, nparx, fIstat); - - //*fLog << "MHFindSignificance::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = " - // << fmin << ", " << fedm << ", " << errdef << ", " << npari - // << ", " << nparx << ", " << fIstat << endl; - - - //------------------- - // store the results - - Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0; - fChisq = fmin; - fNdf = fMbins - nparfree; - fProb = TMath::Prob(fChisq, fNdf); - - - // get fitted parameter values and errors - fValues.Set(fDegree+1); - fErrors.Set(fDegree+1); - - for (Int_t j=0; j<=fDegree; j++) - { - Double_t val, err; - if (gMinuit) - gMinuit->GetParameter(j, val, err); - - fPoly->SetParameter(j, val); - fPoly->SetParError(j, err); - - fValues[j] = val; - fErrors[j] = err; - } - - //-------------------------------------------------- - // if the highest coefficient (j0) of the polynomial - // is consistent with zero reduce the degree of the polynomial - - Int_t j0 = 0; - for (Int_t j=fDegree; j>1; j--) - { - // ignore fixed parameters - if (fErrors[j] == 0.0) - continue; - - // this is the highest coefficient - j0 = j; - break; - } - - if ( !fReduceDegree || j0 == 0 || fabs(fValues[j0]) > fErrors[j0] ) - break; - - // reduce the degree of the polynomial - *fLog << "MHFindSignificance::FitPolynomial; reduce the degree of the polynomial from " - << fDegree << " to " << (j0-2) << endl; - fDegree = j0 - 2; - - funclist->Remove(fPoly); - //if (fPoly) - delete fPoly; - fPoly = NULL; - - // delete the Minuit object in order to have independent starting - // conditions for the next minimization - //if (gMinuit) - delete gMinuit; - gMinuit = NULL; + //-------------------------------------------------- + // prepare fit of a polynomial : (a0 + a1*x + a2*x**2 + a3*x**3 + ...) + + TString funcname = "Poly"; + Double_t xmin = 0.0; + Double_t xmax = 90.0; + + TString formula = "[0]"; + TString bra1 = "+["; + TString bra2 = "]"; + TString xpower = "*x"; + TString newpower = "*x"; + for (Int_t i=1; i<=fDegree; i++) + { + formula += bra1; + formula += i; + formula += bra2; + formula += xpower; + + xpower += newpower; + } + + //*fLog << "FitPolynomial : formula = " << formula << endl; + + fPoly = new TF1(funcname, formula, xmin, xmax); + TList *funclist = fHist->GetListOfFunctions(); + funclist->Add(fPoly); + + //------------------------ + // attention : the dimensions must agree with those in CallMinuit() + const UInt_t npar = fDegree+1; + + TString parname[npar]; + TArrayD vinit(npar); + TArrayD step(npar); + TArrayD limlo(npar); + TArrayD limup(npar); + TArrayI fix(npar); + + vinit[0] = mean; + vinit[2] = a2init; + + for (UInt_t j=0; jGetEntries(); + + // use the subsequernt loop if you want to apply the + // constraint : uneven derivatives (at alpha=0) = zero + for (UInt_t j=1; jmnstat(fmin, fedm, errdef, npari, nparx, fIstat); + + //*fLog << "MHFindSignificance::FitPolynomial; fmin, fedm, errdef, npari, nparx, fIstat = " + // << fmin << ", " << fedm << ", " << errdef << ", " << npari + // << ", " << nparx << ", " << fIstat << endl; + + + //------------------- + // store the results + + Int_t nparfree = gMinuit!=NULL ? gMinuit->GetNumFreePars() : 0; + fChisq = fmin; + fNdf = fMbins - nparfree; + fProb = TMath::Prob(fChisq, fNdf); + + + // get fitted parameter values and errors + fValues.Set(npar); + fErrors.Set(npar); + + for (Int_t j=0; j<=fDegree; j++) + { + Double_t val, err; + if (gMinuit) + gMinuit->GetParameter(j, val, err); + + fPoly->SetParameter(j, val); + fPoly->SetParError(j, err); + + fValues[j] = val; + fErrors[j] = err; + } + + //-------------------------------------------------- + // if the highest coefficient (j0) of the polynomial + // is consistent with zero reduce the degree of the polynomial + + Int_t j0 = 0; + for (Int_t j=fDegree; j>1; j--) + { + // ignore fixed parameters + if (fErrors[j] == 0) + continue; + + // this is the highest coefficient + j0 = j; + break; + } + + if (!fReduceDegree || j0==0 || TMath::Abs(fValues[j0]) > fErrors[j0]) + break; + + // reduce the degree of the polynomial + *fLog << "MHFindSignificance::FitPolynomial; reduce the degree of the polynomial from " + << fDegree << " to " << (j0-2) << endl; + fDegree = j0 - 2; + + funclist->Remove(fPoly); + //if (fPoly) + delete fPoly; + fPoly = NULL; + + // delete the Minuit object in order to have independent starting + // conditions for the next minimization + //if (gMinuit) + delete gMinuit; + gMinuit = NULL; } //=========== end of loop for reducing the degree ================== // of the polynomial - + //-------------------------------------------------- @@ -1034,42 +1015,40 @@ if (fIstat >= 1) { - // error matrix was calculated - if (gMinuit) - gMinuit->mnemat(&fEmat[0][0], fNdim); - - // copy covariance matrix into a matrix which includes also the fixed - // parameters - TString name; - Double_t bnd1, bnd2, val, err; - Int_t jvarbl; - Int_t kvarbl; - for (Int_t j=0; j<=fDegree; j++) - { + // error matrix was calculated if (gMinuit) - gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl); - for (Int_t k=0; k<=fDegree; k++) - { - if (gMinuit) - gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl); - if (jvarbl == 0 || kvarbl == 0) - fEma[j][k] = 0.0; - else - fEma[j][k] = fEmat[jvarbl-1][kvarbl-1]; + gMinuit->mnemat(&fEmat[0][0], fNdim); + + // copy covariance matrix into a matrix which includes also the fixed + // parameters + TString name; + Double_t bnd1, bnd2, val, err; + Int_t jvarbl; + Int_t kvarbl; + for (Int_t j=0; j<=fDegree; j++) + { + if (gMinuit) + gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl); + + for (Int_t k=0; k<=fDegree; k++) + { + if (gMinuit) + gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl); + + fEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fEmat[jvarbl-1][kvarbl-1]; + } } - } } else { - // error matrix was not calculated, construct it - *fLog << "MHFindSignificance::FitPolynomial; error matrix not defined" + // error matrix was not calculated, construct it + *fLog << "MHFindSignificance::FitPolynomial; error matrix not defined" << endl; - for (Int_t j=0; j<=fDegree; j++) - { - for (Int_t k=0; k<=fDegree; k++) - { - fEma[j][k] = 0.0; + for (Int_t j=0; j<=fDegree; j++) + { + for (Int_t k=0; k<=fDegree; k++) + fEma[j][k] = 0; + + fEma[j][j] = fErrors[j]*fErrors[j]; } - fEma[j][j] = fErrors[j]*fErrors[j]; - } } @@ -1078,13 +1057,9 @@ // calculate correlation matrix for (Int_t j=0; j<=fDegree; j++) - { - for (Int_t k=0; k<=fDegree; k++) - { - if (fEma[j][j]*fEma[k][k] != 0.0) - fCorr[j][k] = fEma[j][k] / sqrt( fEma[j][j]*fEma[k][k] ); - else - fCorr[j][k] = 0.0; - } - } + for (Int_t k=0; k<=fDegree; k++) + { + const Double_t sq = fEma[j][j]*fEma[k][k]; + fCorr[j][k] = sq==0 ? 0 : fEma[j][k] / TMath::Sqrt(fEma[j][j]*fEma[k][k]); + } @@ -1092,8 +1067,5 @@ // reset the errors of the points in the histogram for (Int_t i=1; i<=nbins; i++) - { - fHist->SetBinError(i, saveError[i-1]); - } - + fHist->SetBinError(i, saveError[i-1]); return kTRUE; @@ -1114,13 +1086,13 @@ // search bin i0 which has x0 as lower edge - Int_t i0 = -1; - Int_t nbold = fHistOrig->GetNbinsX(); + Int_t i0 = -1; + Int_t nbold = fHistOrig->GetNbinsX(); for (Int_t i=1; i<=nbold; i++) { - if ( fabs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 ) - { - i0 = i; - break; - } + if (TMath::Abs(fHistOrig->GetBinLowEdge(i) - x0) < 1.e-4 ) + { + i0 = i; + break; + } } @@ -1137,8 +1109,7 @@ // get new bin edges - Int_t nbnew = (nbold-istart+1) / nrebin; - Double_t xmin = fHistOrig->GetBinLowEdge(istart); - Double_t xmax = xmin + - ((Double_t)nbnew) * ((Double_t)nrebin) * fHistOrig->GetBinWidth(1); + const Int_t nbnew = (nbold-istart+1) / nrebin; + const Double_t xmin = fHistOrig->GetBinLowEdge(istart); + const Double_t xmax = xmin + (Double_t)nbnew * nrebin * fHistOrig->GetBinWidth(1); fHist->SetBins(nbnew, xmin, xmax); @@ -1152,17 +1123,17 @@ for (Int_t i=1; i<=nbnew; i++) { - Int_t j = nrebin*(i-1) + istart; - - Double_t content = 0.0; - Double_t error2 = 0.0; - for (Int_t k=0; kGetBinContent(j+k); - error2 += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k); - } - fHist->SetBinContent(i, content); - fHist->SetBinError (i, sqrt(error2)); - } - fHist->SetEntries( fHistOrig->GetEntries() ); + Int_t j = nrebin*(i-1) + istart; + + Double_t content = 0; + Double_t error2 = 0; + for (Int_t k=0; kGetBinContent(j+k); + error2 += fHistOrig->GetBinError(j+k) * fHistOrig->GetBinError(j+k); + } + fHist->SetBinContent(i, content); + fHist->SetBinError (i, sqrt(error2)); + } + fHist->SetEntries( fHistOrig->GetEntries() ); return kTRUE; @@ -1286,6 +1257,4 @@ Double_t xmax = 90.0; - TString bra1 = "["; - TString bra2 = "]"; TString xpower = "*x"; TString newpower = "*x"; @@ -1293,30 +1262,8 @@ TString formulaBackg = "[0]"; for (Int_t i=1; i<=fDegree; i++) - { - formulaBackg += "+"; - formulaBackg += bra1; - formulaBackg += i; - formulaBackg += bra2; - formulaBackg += xpower; - - xpower += newpower; - } - - TString formulaGauss = bra1; - formulaGauss += fDegree+1; - formulaGauss += bra2; - formulaGauss += "/"; - formulaGauss += bra1; - formulaGauss += fDegree+3; - formulaGauss += bra2; - formulaGauss += "*exp(-0.5*((x-"; - formulaGauss += bra1; - formulaGauss += fDegree+2; - formulaGauss += bra2; - formulaGauss += ")/"; - formulaGauss += bra1; - formulaGauss += fDegree+3; - formulaGauss += bra2; - formulaGauss += ")^2)"; + formulaBackg += Form("+[%d]^%d", i, i); + + const TString formulaGauss = Form("[%d]/[%d]*exp(-0.5*((x-[%d])/[%d])^2)", + fDegree+1, fDegree+3, fDegree+2, fDegree+3); TString formula = formulaBackg; @@ -1337,27 +1284,22 @@ //------------------------ // attention : the dimensions must agree with those in CallMinuit() - - char parname[80][100]; - Double_t vinit[80]; - Double_t step[80]; - Double_t limlo[80]; - Double_t limup[80]; - Int_t fix[80]; - Int_t npar = fDegree+1 + 3; + + TString parname[npar]; + TArrayD vinit(npar); + TArrayD step(npar); + TArrayD limlo(npar); + TArrayD limup(npar); + TArrayI fix(npar); + // take as initial values for the polynomial // the result from the polynomial fit for (Int_t j=0; j<=fDegree; j++) - { vinit[j] = fPoly->GetParameter(j); - } - - - - Double_t sigma = 8.0; - vinit[fDegree+1] = 2.0 * fNex * fHist->GetBinWidth(1) / - sqrt( 2.0*TMath::Pi() ); - vinit[fDegree+2] = 0.0; + + Double_t sigma = 8; + vinit[fDegree+1] = 2.0 * fNex * fHist->GetBinWidth(1) / TMath::Sqrt(TMath::Pi()*2); + vinit[fDegree+2] = 0; vinit[fDegree+3] = sigma; @@ -1367,40 +1309,31 @@ for (Int_t j=0; jGetEntries(); + limup[0] = fHist->GetEntries()*10; // limit the sigma of the Gauss function - limlo[fDegree+3] = 0.0; - limup[fDegree+3] = 20.0; + limup[fDegree+3] = 20; // use the subsequernt loop if you want to apply the // constraint : uneven derivatives (at alpha=0) = zero - for (Int_t j=1; j<=fDegree; j += 2) - { - vinit[j] = 0.0; - step[j] = 0.0; - fix[j] = 1; + for (Int_t j=1; j<=fDegree; j+=2) + { + vinit[j] = 0; + step[j] = 0; + fix[j] = 1; } // fix position of Gauss function - vinit[fDegree+2] = 0.0; - step[fDegree+2] = 0.0; - fix[fDegree+2] = 1; + vinit[fDegree+2] = 0; + step[fDegree+2] = 0; + fix[fDegree+2] = 1; // if a constant background has been assumed (due to low statistics) @@ -1408,16 +1341,12 @@ if (fConstantBackg) { - step[0] = 0.0; - fix[0] = 1; - } - - TString method = "Migrad"; - TObject *objectfit = fHist; - Bool_t nulloutput = kFALSE; + step[0] = 0; + fix[0] = 1; + } MMinuitInterface inter; - Bool_t rc = inter.CallMinuit( fcnpolygauss, npar, parname, vinit, step, - limlo, limup, fix, objectfit, method, - nulloutput); + const Bool_t rc = inter.CallMinuit(fcnpolygauss, parname, vinit, step, + limlo, limup, fix, fHist, "Migrad", + kFALSE); if (rc != 0) @@ -1480,6 +1409,6 @@ fdSigmaGauss = fGErrors[fDegree+3]; // fitted total number of excess events - fNexGauss = fGValues[fDegree+1] * sqrt(2.0*TMath::Pi()) / - ( 2.0 * fHist->GetBinWidth(1) ); + fNexGauss = fGValues[fDegree+1] * TMath::Sqrt(TMath::Pi()*2) / + (fHist->GetBinWidth(1)*2 ); fdNexGauss = fNexGauss * fGErrors[fDegree+1]/fGValues[fDegree+1]; @@ -1502,15 +1431,14 @@ for (Int_t j=0; jmnpout(j, name, val, err, bnd1, bnd2, jvarbl); - for (Int_t k=0; kmnpout(k, name, val, err, bnd1, bnd2, kvarbl); - if (jvarbl == 0 || kvarbl == 0) - fGEma[j][k] = 0.0; - else - fGEma[j][k] = fGEmat[jvarbl-1][kvarbl-1]; - } + gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl); + + for (Int_t k=0; kmnpout(k, name, val, err, bnd1, bnd2, kvarbl); + + fGEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fGEmat[jvarbl-1][kvarbl-1]; + } } } @@ -1522,9 +1450,8 @@ for (Int_t j=0; jSetBinError(i, saveError[i-1]); - } return kTRUE; Index: /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h (revision 2310) +++ /trunk/MagicSoft/Mars/mhist/MHFindSignificance.h (revision 2311) @@ -1,9 +1,4 @@ #ifndef MARS_MHFindSignificance #define MARS_MHFindSignificance - -#ifdef MARS_MLogManip -#error Please make ensure that MLogManip.h are included _after_ MHFindSignificance.h -#endif - #ifndef MARS_MH @@ -11,11 +6,11 @@ #endif +#ifndef ROOT_TArrayD #include -#include -#include +#endif -class TArrayD; class TF1; - +class TH1; +class TCanvas; class MHFindSignificance : public MH