Changeset 2311 for trunk/MagicSoft/Mars/mhist

Timestamp:

08/21/03 11:55:13 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHFindSignificance.cc (modified) (18 diffs)
• MHFindSignificance.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc

@@ -61 +61 @@
 
 #include <TArrayD.h>
+#include <TArrayI.h>
 #include <TH1.h>
 #include <TF1.h>
@@ -849 +850 @@
   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; j<npar; j++)
-  {
-    vinit[j] = 0.0;
-  }
-  vinit[0] =   mean;
-  vinit[2] = a2init;
-
-
-  for (UInt_t j=0; j<npar; j++)
-  {
-    sprintf(&parname[j][0], "p%d", j+1);
-
-    if (vinit[j] != 0.0)
-      step[j] = fabs(vinit[j]) / 10.0;
-    else
-      step[j] = 0.000001;
-
-    limlo[j]  = 0.0;
-    limup[j]  = 0.0;
-
-    fix[j]    =   0;
-  }
-
-  // limit the first coefficient of the polynomial to positive values
-  // because the background must not be negative
-  limlo[0]  = 0.0;
-  limup[0]  = fHist->GetEntries();
-
-
-  // use the subsequernt loop if you want to apply the
-  // constraint : uneven derivatives (at alpha=0) = zero
-  for (UInt_t j=1; j<npar; j += 2)
-  {
-    vinit[j] = 0.0;
-    step[j]  = 0.0;
-    fix[j]   =   1;
-  }
-
-  TString method     = "Migrad";
-  TObject *objectfit = fHist;
-  Bool_t  nulloutput = kFALSE;
-
-  *fLog << "FitPolynomial : before CallMinuit()" << endl;
-
-  MMinuitInterface inter;
-  Bool_t rc = inter.CallMinuit( fcnpoly, npar, parname,     vinit, step,
-                                limlo,   limup,    fix, objectfit, method,
-                                nulloutput);
-
-  if (rc != 0)
-  {
-  //  *fLog << "MHFindSignificance::FitPolynomial; polynomial fit failed"
-  //        << endl;
-  //  return kFALSE;
-  }
-
-
-  //-------------------
-  // get status of minimization
-  Double_t fmin   = 0.0;
-  Double_t fedm   = 0.0;
-  Double_t errdef = 0.0;
-  Int_t    npari  =   0;
-  Int_t    nparx  =   0;
- 
-  if (gMinuit)
-    gMinuit->mnstat(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; j<npar; j++)
+      {
+          parname[j]  = "p";
+          parname[j] += j+1;
+
+          step[j] = vinit[j] != 0.0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
+      }
+
+      // limit the first coefficient of the polynomial to positive values
+      // because the background must not be negative
+      limup[0] = fHist->GetEntries();
+
+      // use the subsequernt loop if you want to apply the
+      // constraint : uneven derivatives (at alpha=0) = zero
+      for (UInt_t j=1; j<npar; j+=2)
+      {
+          vinit[j] = 0;
+          step[j]  = 0;
+          fix[j]   = 1;
+      }
+
+      *fLog << "FitPolynomial : before CallMinuit()" << endl;
+
+      MMinuitInterface inter;
+      const Bool_t rc = inter.CallMinuit(fcnpoly, parname, vinit, step,
+                                         limlo, limup, fix, fHist, "Migrad",
+                                         kFALSE);
+
+      if (rc != 0)
+      {
+          //  *fLog << "MHFindSignificance::FitPolynomial; polynomial fit failed"
+          //        << endl;
+          //  return kFALSE;
+      }
+
+
+      //-------------------
+      // get status of minimization
+      Double_t fmin   = 0;
+      Double_t fedm   = 0;
+      Double_t errdef = 0;
+      Int_t    npari  = 0;
+      Int_t    nparx  = 0;
+
+      if (gMinuit)
+          gMinuit->mnstat(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 +1015 @@
   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 +1057 @@
   // 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 +1067 @@
   // 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 +1086 @@
   // 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 +1109 @@
   // 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 +1123 @@
   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; k<nrebin; k++)
-    {
-      content += fHistOrig->GetBinContent(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; k<nrebin; k++)
+      {
+          content += fHistOrig->GetBinContent(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 +1257 @@
   Double_t xmax =  90.0;
 
-  TString bra1     = "[";
-  TString bra2     =    "]";
   TString xpower   = "*x";
   TString newpower = "*x";
@@ -1293 +1262 @@
   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 +1284 @@
   //------------------------
   // 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 +1309 @@
   for (Int_t j=0; j<npar; j++)
   {
-    sprintf(&parname[j][0], "p%d", j+1);
-
-    if (vinit[j] != 0.0)
-      step[j] = fabs(vinit[j]) / 10.0;
-    else
-      step[j] = 0.000001;
-
-    limlo[j]  = 0.0;
-    limup[j]  = 0.0;
-
-    fix[j]    =   0;
+      parname[j]  = "p";
+      parname[j] += j+1;
+
+      step[j] = vinit[j]!=0 ? TMath::Abs(vinit[j]) / 10.0 : 0.000001;
   }
 
   // limit the first coefficient of the polynomial to positive values
   // because the background must not be negative
-  limlo[0]  = 0.0;
-  limup[0]  = 10.0 * fHist->GetEntries();
+  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 +1341 @@
   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 +1409 @@
   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 +1431 @@
     for (Int_t j=0; j<npar; j++)
     {
-      if (gMinuit)
-        gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
-      for (Int_t k=0; k<npar; k++)
-      {
         if (gMinuit)
-          gMinuit->mnpout(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; k<npar; k++)
+        {
+            if (gMinuit)
+                gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
+
+            fGEma[j][k] = jvarbl==0 || kvarbl==0 ? 0 : fGEmat[jvarbl-1][kvarbl-1];
+        }
     }
   }
@@ -1522 +1450 @@
     for (Int_t j=0; j<npar; j++)
     {
-      for (Int_t k=0; k<npar; k++)
-      {
-        fGEma[j][k] = 0.0; 
-      }
-      fGEma[j][j] = fGErrors[j]*fGErrors[j];
+        for (Int_t k=0; k<npar; k++)
+            fGEma[j][k] = 0;
+
+        fGEma[j][j] = fGErrors[j]*fGErrors[j];
     }
   }
@@ -1537 +1464 @@
     for (Int_t k=0; k<npar; k++)
     {
-      if (fGEma[j][j]*fGEma[k][k] != 0.0)
-        fGCorr[j][k] = fGEma[j][k] / sqrt( fGEma[j][j]*fGEma[k][k] ); 
-      else
-        fGCorr[j][k] = 0.0;
+        const Double_t sq = fGEma[j][j]*fGEma[k][k];
+        fGCorr[j][k] = sq==0 ? 0 : fGEma[j][k] / sqrt( fGEma[j][j]*fGEma[k][k] );
     }
   }
@@ -1548 +1473 @@
   // reset the errors of the points in the histogram
   for (Int_t i=1; i<=nbins; i++)
-  {
     fHist->SetBinError(i, saveError[i-1]);
-  }
 
   return kTRUE;

trunk/MagicSoft/Mars/mhist/MHFindSignificance.h

@@ -1 +1 @@
 #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 +6 @@
 #endif
 
+#ifndef ROOT_TArrayD
 #include <TArrayD.h>
-#include <TH1.h>
-#include <TCanvas.h>
+#endif
 
-class TArrayD;
 class TF1;
-
+class TH1;
+class TCanvas;
 
 class MHFindSignificance : public MH