NuSystFitter Class Reference

#include <NuSystFitter.h>

Inheritance diagram for NuSystFitter:

List of all members.

Public Member Functions
	NuSystFitter ()
virtual	~NuSystFitter ()
virtual double	Up () const
virtual void	SetErrorDef (const Double_t up)
virtual double	operator() (const std::vector< double > &pars) const
virtual double	Likelihood (const NuMMParameters &mmPars) const
virtual void	Setup (const std::vector< std::string > helperInputs, const std::vector< NuHistInterpolator * > fdDataInput, const std::vector< NuHistInterpolator * > ndDataInput, const std::string xsecFile, const std::vector< std::string >, const std::string, Int_t)
virtual void	push_back (NuMMRun *run)
virtual const NuMMParameters	Minimise (const NuMMParameters &mmPars)
virtual const NuMMParameters	Minimise (const NuMMParameters &mmPars, bool &success)
virtual const NuMMParameters	FCMinimise (const NuMMParameters &mmPars)
	More robustly minimizes by doing a wide grid search first - same as the FC fitting.
virtual const NuMMParameters	FCMinimiseNu (const NuMMParameters &mmPars)
virtual const NuMMParameters	FCMinimiseNSI (const NuMMParameters &mmPars)
virtual const TGraph *	PlotContours (const NuMMParameters &mmPars)
virtual void	NoOscPrediction ()
virtual void	NumberOfContourPoints (const Int_t numContourPoints)
virtual const std::pair< Double_t, Double_t >	SingleParDm2Errors (const NuMMParameters &mmPars, const bool quiet=false)
virtual const std::pair< Double_t, Double_t >	SingleParDm2BarErrors (const NuMMParameters &mmPars, const bool quiet=false)
virtual const std::pair< Double_t, Double_t >	SingleParSn2BarErrors (const NuMMParameters &mmPars, const bool quiet=false)
virtual const std::pair< Double_t, Double_t >	SingleParErrors (const NuMMParameters &mmPars, const int idx, const bool quiet)
virtual const std::vector< std::pair< Double_t, Double_t > >	NeutrinoContourFinder (const NuMMParameters &mmPars)
virtual const std::pair< Double_t, Double_t >	DmScanForContour (const NuMMParameters &mmPars, Double_t targetChi2)
virtual void	Sn2PointsForNeutrinoContour (const std::vector< Double_t > &sn2PointsForNeutrinoContour)
virtual const TGraph *	OneSidedNeutrinoContourFinder (const NuMMParameters &inputPars)
virtual const TGraph *	TwoSidedNeutrinoContourFinder (const NuMMParameters &inputPars)
virtual const TGraph *	Chi2SliceSin2bar (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minval=0.8, const Double_t maxval=1.2)
virtual const TGraph *	Chi2SliceDm2bar (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minval=1e-3, const Double_t maxval=3e-3)
virtual const TGraph *	Chi2SliceDm2 (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minval=1e-3, const Double_t maxval=3e-3)
virtual const TGraph *	Chi2ValleyBar (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minSin2bar=0., const Double_t maxSin2bar=2.)
virtual const TGraph *	Chi2ValleyBarDm2Bar (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minDm2bar=0., const Double_t maxDm2bar=10e-3)
virtual const TGraph *	Chi2ValleyDm2 (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minDm2=0., const Double_t maxDm2=10e-3)
virtual const TGraph *	Chi2ValleySn2 (const NuMMParameters &mmPars, const Int_t numpoints=100, const Double_t minSn2=0., const Double_t maxSn2=1.)
virtual TH2D *	LikelihoodSurfaceBar (const NuMMParameters &mmPars, Int_t sinPoints, Int_t massPoints, Double_t Sn2BarMin=0.0, Double_t Sn2BarMax=1.0, Double_t Dm2BarMin=0.0, Double_t Dm2BarMax=12e-3, Bool_t True_Minimum=true)
	Function to generate a likelihood surface.
virtual Double_t	FillLikelihoodSurface (TH2D &surface, NuMMParameters mmPars)
	Fills a histogram with likelihood values. This function accepts an arbitrarily binned histogram, and fills the values according to the likelihood. It returns the value of the minimum likelihood. (This does the CC neutirno spectrum).
virtual Double_t	FillLikelihoodSurfaceBar (TH2D &surface, NuMMParameters mmPars)
	Fills a histogram with likelihood values. This function accepts an arbitrarily binned histogram, and fills the values according to the likelihood. It returns the value of the minimum likelihood.
virtual Double_t	FillLikelihoodSurfaceNSI (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair)
virtual Double_t	FillLikelihoodSurfaceNSIBar (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair)
virtual Double_t	FillLikelihoodSurfaceNSIDm2Eps (TH2D &surface, NuMMParameters mmPars, TH2D &sinvals)
virtual Double_t	FillLikelihoodSurfaceNSISn2Eps (TH2D &surface, NuMMParameters mmPars, TH2D &massvals)
virtual Double_t	FillLikelihoodSurfaceNSIDm2Sn2 (TH2D &surface, NuMMParameters mmPars, TH2D &epsvals)
virtual Double_t	FillLikelihoodSurfaceNSI3D (TH3D &surface, NuMMParameters mmPars)
virtual Double_t	FillLikelihoodSurfaceDecay (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair=1)
virtual Double_t	FillLikelihoodSurfaceDecoherence (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair=1)
virtual Double_t	FillLikelihoodSurfaceDecoherenceBar (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair=1)
virtual Double_t	FillLikelihoodSurfaceCPT (TH2D &surface, NuMMParameters mmPars)
virtual Double_t	FillLikelihoodSurfaceSterile (TH2D &surface, NuMMParameters mmPars, Int_t ContourPair)
virtual void	BatchModeOn ()
	Supresses much of the usual output.
virtual void	QuietModeOn ()
virtual void	QuietModeOff ()
virtual void	FillCPT (TH2D &surface)
	Fills a 2D histogram that has dimensions Dm2 (Y) and DM2Bar(X). Uses the center of the passed bins.
NuMMParameters	FillAlphaBeta (TH2D &J, TH2D &K, const std::string outfilename)
Protected Member Functions
virtual Double_t	FillLikelihoodSurfaceGeneric (TH2D &surface, NuMMParameters mmPars)
virtual Double_t	FillLikelihoodSurfaceGeneric (TH3D &surface, NuMMParameters mmPars)
NuMMParameters	CoarseGridSearch (NuMMParameters mmGrid)
NuMMParameters	GridSearch (TH2D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref=0)
NuMMParameters	FineGridSearch (NuMMParameters mmStart)
NuMMParameters	CoarseGridSearchNu (NuMMParameters mmGrid)
NuMMParameters	GridSearchNu (TH2D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref=0)
NuMMParameters	FineGridSearchNu (NuMMParameters mmStart)
NuMMParameters	CoarseGridSearchNSI (NuMMParameters mmGrid)
NuMMParameters	GridSearchNSI (TH2D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref=0)
NuMMParameters	FineGridSearchNSI (NuMMParameters mmStart)
NuMMParameters	Sn2GridSearch (Double_t Dm2, Double_t Dm2Bar, Int_t resolution=0)
	Does a grid search over sn2/sn2bar, with resolution bins per parameter.
Protected Attributes
Int_t	fQuietMode
Double_t	fUp
Int_t	fNumContourPoints
std::vector< Double_t >	fsn2PointsForNeutrinoContour
std::vector< Double_t >	fsn2PointsForNeutrinoContourPhysical
std::vector< Double_t >	fsn2PointsForNeutrinoContourUnphysical
std::vector< NuMMRun * > *	fRuns

---

Detailed Description

Definition at line 18 of file NuSystFitter.h.

---

Constructor & Destructor Documentation

NuSystFitter::NuSystFitter

(

)

Definition at line 29 of file NuSystFitter.cxx.

  : fQuietMode(0),
    fUp(1.0),
    fNumContourPoints(80),
    fsn2PointsForNeutrinoContour(),
    fsn2PointsForNeutrinoContourPhysical(),
    fsn2PointsForNeutrinoContourUnphysical()
{
  fRuns = new vector<NuMMRun*>();

  for (Double_t sn2 = 1.0; sn2 > 0.871; sn2 -= 0.02){
    fsn2PointsForNeutrinoContourPhysical.push_back(sn2);
  }
  fsn2PointsForNeutrinoContourPhysical.push_back(0.87);
  fsn2PointsForNeutrinoContourPhysical.push_back(0.86);
  for (Double_t sn2 = 0.8595; sn2 > 0.7999; sn2 -= 0.0005){
    fsn2PointsForNeutrinoContourPhysical.push_back(sn2);
  }

  for (Double_t sn2 = 1.0; sn2 < 3; sn2 += 0.02){
    fsn2PointsForNeutrinoContourUnphysical.push_back(sn2);
  }

  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourPhysical;
}

NuSystFitter::~NuSystFitter

(

)

[virtual]

Definition at line 56 of file NuSystFitter.cxx.

{
}

---

Member Function Documentation

void NuSystFitter::BatchModeOn

(

)

[virtual]

Supresses much of the usual output.

Definition at line 655 of file NuSystFitter.cxx.

References fQuietMode, and QuietModeOn().

Referenced by NuFCFitter::NuFCFitter(), NuFCFitterNSI::NuFCFitterNSI(), NuFCFitterNSINu::NuFCFitterNSINu(), and NuFCFitterNSINubar::NuFCFitterNSINubar().

{
  QuietModeOn();
  fQuietMode = 2;
}

const TGraph * NuSystFitter::Chi2SliceDm2	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minval = 1e-3,
		const Double_t	maxval = 3e-3
	)			[virtual]

Definition at line 2187 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), fQuietMode, minimum, and NuMMParameters::VectorParameters().

{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;


        if (!fQuietMode) cout << "Generating Likelyhood surface for dm2" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Dm2", "#Delta#chi^{2} Surface for #Deltam^{2}");
        chisurf->GetXaxis()->SetTitle("#Deltam^{2}");

        Double_t dm2value = 0, like = 0, minimum = 0;

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                dm2value = minval + increment*i;
                mmPars.Dm2(dm2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());

                // Fill the array
//              dm2s[i] = dm2value;
//              likelyhoods[i] = like;

                chisurf->SetPoint(i, dm2value, like);

                // Is this the lowest value?
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }

        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}

const TGraph * NuSystFitter::Chi2SliceDm2bar	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minval = 1e-3,
		const Double_t	maxval = 3e-3
	)			[virtual]

Definition at line 738 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), fQuietMode, minimum, and NuMMParameters::VectorParameters().

{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;


        if (!fQuietMode) cout << "Generating Likelyhood surface for dm2bar" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Dm2Bar", "#Delta#chi^{2} Surface for #Delta#barm^{2}");
        chisurf->GetXaxis()->SetTitle("#Delta#barm^{2}");

        Double_t dm2value = 0, like = 0, minimum = 0;

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                dm2value = minval + increment*i;
                mmPars.Dm2Bar(dm2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());

                // Fill the array
//              dm2s[i] = dm2value;
//              likelyhoods[i] = like;

                chisurf->SetPoint(i, dm2value, like);

                // Is this the lowest value?
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }

        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}

const TGraph * NuSystFitter::Chi2SliceSin2bar	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minval = 0.8,
		const Double_t	maxval = 1.2
	)			[virtual]

Definition at line 684 of file NuSystFitter.cxx.

References fQuietMode, minimum, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;

        if (!fQuietMode) cout << "Generating Likelyhood surface for sin2bar" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        Double_t sn2value = 0, like = 0, minimum = 0;

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Sin2Bar", "#Delta#chi^{2} Surface for Sin^{2}2#bar#theta");
        chisurf->GetXaxis()->SetTitle("Sin^{2}2#bar#theta");

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                sn2value = minval + increment*i;
                mmPars.Sn2Bar(sn2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());

                // Fill the array
                //sin2s[i] = sn2value;
                //likelyhoods[i] = like;
                chisurf->SetPoint(i, sn2value, like);

                // Is this the lowest value?
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }

        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}

const TGraph * NuSystFitter::Chi2ValleyBar	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minSin2bar = 0.,
		const Double_t	maxSin2bar = 2.
	)			[virtual]

Definition at line 793 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), NuMMParameters::FixDm2(), NuMMParameters::FixDm2Bar(), NuMMParameters::FixNCBackground(), NuMMParameters::FixNorm(), NuMMParameters::FixShwScale(), NuMMParameters::FixSn2(), NuMMParameters::FixSn2Bar(), fQuietMode, Minimise(), minimum, NuMMParameters::MinuitPass(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

{
    NuMMParameters mmPars(mmParsC);

    //  mmPars.Dm2Bar(2.5e-3);
    // mmPars.Sn2Bar(1.0);
    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.FixDm2();
    mmPars.FixSn2();
    mmPars.FixDm2Bar(false);
    mmPars.FixSn2Bar();

    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t sinval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot..." << endl;

    TGraph *sin2like = new TGraph(points);

    for (Int_t i = 0; i < points; i++) {

        //        cout << endl << "Sin2 " << i+1 << " / " << points << endl;
        // Set up the parameter structure
        sinval = min + interval*i;
        mmPars.Sn2Bar(sinval);

        NuMMParameters mmFit;
        if (sinval == 0) {
            mmFit = mmPars;
        }
        else {
            // Now, minimize the likelyhood for this value
            mmFit = this->Minimise(mmPars);
            int failCount = 0;
            while (!mmFit.MinuitPass()) {
                if (failCount > 10) break;
                failCount++;
                mmPars.Dm2Bar(mmPars.Dm2Bar()*1.1);
                cout << "Starting at dm2= " << mmPars.Dm2Bar() << endl;
                mmFit = this->Minimise(mmPars);
            }
            mmPars.Dm2Bar(mmParsC.Dm2Bar());
            if (failCount > 10) continue;
        }

        likeval = (*this)(mmFit.VectorParameters());
        sin2like->SetPoint(i, sinval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
        //              cout << "Minimum is now " << minimum << endl;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;
    // Offset the data points
    for (Int_t i = 0; i < points; i++)
        {
        Double_t x, y;
        sin2like->GetPoint(i, x, y);
        //              cout << "From: " << x << ", " << y << " to ";
        sin2like->SetPoint(i, x, y-minimum);
        //              cout << x << ", " << y-minimum << endl;
        }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    // Make the graphs
    return sin2like;
}

const TGraph * NuSystFitter::Chi2ValleyBarDm2Bar	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minDm2bar = 0.,
		const Double_t	maxDm2bar = 10e-3
	)			[virtual]

Definition at line 873 of file NuSystFitter.cxx.

References NuMMParameters::ConstrainPhysical(), NuMMParameters::Dm2Bar(), NuMMParameters::FixDm2(), NuMMParameters::FixDm2Bar(), NuMMParameters::FixNCBackground(), NuMMParameters::FixNorm(), NuMMParameters::FixShwScale(), NuMMParameters::FixSn2(), NuMMParameters::FixSn2Bar(), fQuietMode, Minimise(), minimum, NuMMParameters::MinuitPass(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.FixDm2();
    mmPars.FixSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar(false);
    mmPars.ConstrainPhysical();


    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t dmval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for dm2..." << endl;

    TGraph *dm2like = new TGraph(points);
    dm2like->SetName("valleyplot_mass");

    for (Int_t i = 0; i < points; i++) {
        if (i % 10 == 0) {
            cout << i << "/" << points << endl;
        }
        // Set up the parameter structure
        dmval = min + interval*i;
        mmPars.Dm2Bar(dmval);

        NuMMParameters mmFit;
        if (dmval == 0) {
            mmFit = mmPars;
        }
        else {
            // Now, minimize the likelyhood for this value
            mmFit = this->Minimise(mmPars);
            int failCount = 0;
            while (!mmFit.MinuitPass()) {
                if (failCount > 10) break;
                failCount++;
                mmPars.Sn2Bar(mmPars.Sn2Bar()*0.9);
                cout << "Starting at sn2= " << mmPars.Sn2Bar() << endl;
                mmFit = this->Minimise(mmPars);
            }
            mmPars.Sn2Bar(mmParsC.Sn2Bar());
            if (failCount > 10) continue;
        }
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        dm2like->SetPoint(i, dmval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++) {
        Double_t x, y;
        dm2like->GetPoint(i, x, y);
        dm2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return dm2like;
}

const TGraph * NuSystFitter::Chi2ValleyDm2	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minDm2 = 0.,
		const Double_t	maxDm2 = 10e-3
	)			[virtual]

Definition at line 956 of file NuSystFitter.cxx.

References NuMMParameters::ConstrainPhysical(), NuMMParameters::Dm2(), NuMMParameters::FixDm2(), NuMMParameters::FixDm2Bar(), NuMMParameters::FixNCBackground(), NuMMParameters::FixNorm(), NuMMParameters::FixShwScale(), NuMMParameters::FixSn2Bar(), fQuietMode, Minimise(), minimum, NuUtilities::ProgressBar(), NuMMParameters::ReleaseSn2(), and NuMMParameters::VectorParameters().

{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.FixDm2();
    mmPars.ReleaseSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar();
    mmPars.ConstrainPhysical();


    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t dmval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for dm2..." << endl;

    TGraph *dm2like = new TGraph(points);
    dm2like->SetName("valleyplot_mass");

    for (Int_t i = 0; i < points; i++)
    {
        // Print a progress bar
        NuUtilities::ProgressBar(i, points-1,5);

        // Set up the parameter structure
        dmval = min + interval*i;
        mmPars.Dm2(dmval);

        // Now, minimize the likelyhood for this value
        NuMMParameters mmFit = this->Minimise(mmPars);
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        dm2like->SetPoint(i, dmval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++)
    {
        Double_t x, y;
        dm2like->GetPoint(i, x, y);
        dm2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return dm2like;
}

const TGraph * NuSystFitter::Chi2ValleySn2	(	const NuMMParameters &	mmPars,
		const Int_t	numpoints = 100,
		const Double_t	minSn2 = 0.,
		const Double_t	maxSn2 = 1.
	)			[virtual]

Definition at line 1025 of file NuSystFitter.cxx.

References NuMMParameters::ConstrainPhysical(), NuMMParameters::FixDm2Bar(), NuMMParameters::FixNCBackground(), NuMMParameters::FixNorm(), NuMMParameters::FixShwScale(), NuMMParameters::FixSn2(), NuMMParameters::FixSn2Bar(), fQuietMode, Minimise(), minimum, NuUtilities::ProgressBar(), NuMMParameters::ReleaseDm2(), NuMMParameters::Sn2(), and NuMMParameters::VectorParameters().

{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.ReleaseDm2();
    mmPars.FixSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar();
    mmPars.ConstrainPhysical();


    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t snval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for sn2..." << endl;

    TGraph *sn2like = new TGraph(points);
    sn2like->SetName("valleyplot_sin");

    for (Int_t i = 0; i < points; i++)
    {
        // Print a progress bar
        NuUtilities::ProgressBar(i, points-1,5);

        // Set up the parameter structure
        snval = min + interval*i;
        mmPars.Sn2(snval);

        // Now, minimize the likelyhood for this value
        NuMMParameters mmFit = this->Minimise(mmPars);
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        sn2like->SetPoint(i, snval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++)
    {
        Double_t x, y;
        sn2like->GetPoint(i, x, y);
        sn2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return sn2like;
}

NuMMParameters NuSystFitter::CoarseGridSearch

(

NuMMParameters

mmGrid

)

[protected]

Definition at line 1601 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), MuELoss::e, GridSearch(), Msg::kInfo, MSG, and NuMMParameters::Sn2Bar().

Referenced by FCMinimise().

{
  // The width and height of the grid
  const Int_t gridWidth = 5;
  const Int_t gridHeight = 15;

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, gridHeight, 2e-3, 30e-3);
  // TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, 998*10, 2e-3, 1000e-3);
  mmGrid = GridSearch(likesurf, mmGrid);
  // cout << "Found dm2: " << mmGrid.Dm2Bar() << ", sn2: " << mmGrid.Sn2Bar() << endl;

  // Now do a log search
  vector<Double_t> bins;
  for (Double_t i = -1.5; i <= -0.; i += 0.1) {
    bins.push_back(pow(10, i));
  }
  Int_t bincount = bins.size()-1;
  TH2D logsurf("logsurf", "logsurf", gridWidth, 0.2, 1, bincount, &(bins.front()));
  mmGrid = GridSearch(logsurf, mmGrid, &mmGrid);

  // This was code to do a much sparser grid search instead of the log search
  // TH2D sparser("sparser", "sparser", gridWidth, 0.2, 1, 400, 30e-3, 1000e-3);
  //  mmGrid = GridSearch(sparser, mmGrid, &mmGrid);

  // Output the minimum
  MSG("NuSystFitter", Msg::kInfo) << "Coarse grid search found sn2: "
    << mmGrid.Sn2Bar() << ", dm: " << mmGrid.Dm2Bar() << '\n';

  return mmGrid;
}

NuMMParameters NuSystFitter::CoarseGridSearchNSI

(

NuMMParameters

mmGrid

)

[protected]

Definition at line 1926 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), MuELoss::e, NuMMParameters::Epsilon(), GridSearchNSI(), Msg::kInfo, and MSG.

Referenced by FCMinimiseNSI().

{
  // The width and height of the grid
  const Int_t gridWidth = 5;
  const Int_t gridHeight = 15;

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, 0.01, 0.5, gridHeight, 0e-3, 10e-3);
  mmGrid = GridSearchNSI(likesurf, mmGrid);
  cout << "Found dm2: " << mmGrid.Dm2() << ", epsilon: " << mmGrid.Epsilon() << endl;

  // Now do a log search
  vector<Double_t> bins;
  for (Double_t i = -1.5; i <= -0.; i += 0.1) {
    bins.push_back(pow(10, i));
  }
  Int_t bincount = bins.size()-1;
  TH2D logsurf("logsurf", "logsurf", gridWidth, 0.01, 0.5, bincount, &(bins.front()));
  mmGrid = GridSearchNSI(logsurf, mmGrid, &mmGrid);

  // Output the minimum
  MSG("NuSystFitter", Msg::kInfo) << "Coarse grid search found epsilon: "
    << mmGrid.Epsilon() << ", dm: " << mmGrid.Dm2() << '\n';

  return mmGrid;
}

NuMMParameters NuSystFitter::CoarseGridSearchNu

(

NuMMParameters

mmGrid

)

[protected]

Definition at line 1761 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), MuELoss::e, GridSearchNu(), Msg::kInfo, MSG, and NuMMParameters::Sn2().

Referenced by FCMinimiseNu().

{
  // The width and height of the grid
  const Int_t gridWidth = 5;
  const Int_t gridHeight = 15;
  
  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, gridHeight, 2e-3, 30e-3);
  // TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, 998*10, 2e-3, 1000e-3);
  mmGrid = GridSearchNu(likesurf, mmGrid);
  // cout << "Found dm2: " << mmGrid.Dm2Bar() << ", sn2: " << mmGrid.Sn2Bar() << endl;
  
  // Now do a log search
  vector<Double_t> bins;
  for (Double_t i = -1.5; i <= -0.; i += 0.1) {
    bins.push_back(pow(10, i));
  }
  Int_t bincount = bins.size()-1;
  TH2D logsurf("logsurf", "logsurf", gridWidth, 0.2, 1, bincount, &(bins.front()));
  mmGrid = GridSearchNu(logsurf, mmGrid, &mmGrid);
  
  // This was code to do a much sparser grid search instead of the log search
  // TH2D sparser("sparser", "sparser", gridWidth, 0.2, 1, 400, 30e-3, 1000e-3);
  //  mmGrid = GridSearch(sparser, mmGrid, &mmGrid);
  
  // Output the minimum
  MSG("NuSystFitter", Msg::kInfo) << "Coarse grid search found sn2: "
  << mmGrid.Sn2() << ", dm: " << mmGrid.Dm2() << endl;
  
  return mmGrid;
}

const std::pair< Double_t, Double_t > NuSystFitter::DmScanForContour	(	const NuMMParameters &	mmPars,
		Double_t	targetChi2
	)			[virtual]

Definition at line 362 of file NuSystFitter.cxx.

References NuMMParameters::AreAllParametersFixed(), NuMMParameters::Dm2(), MuELoss::e, NuMMParameters::FixDm2(), NuMMParameters::NCBackgroundScale(), NuMMParameters::Normalisation(), NuMMParameters::ShwEnScale(), NuMMParameters::Sn2(), and NuMMParameters::VectorParameters().

{
  std::pair<Double_t,Double_t> errors;

  NuMMParameters bestFitPoint = this->Minimise(mmPars);
  Double_t bestFitChi2 = (*this)(bestFitPoint.VectorParameters());
  Double_t bestFitDm2 = bestFitPoint.Dm2();

  if (bestFitChi2 > targetChi2){
    errors.first = -1.0;
    errors.second = -1.0;
    return errors;
  }

  cout << "Best chi2: " << bestFitChi2
       << " at dm2 = " << bestFitDm2
       << ", sn2 = " << bestFitPoint.Sn2()
       << ", norm = " << bestFitPoint.Normalisation()
       << ", NC = " << bestFitPoint.NCBackgroundScale()
       << ", shwEn = " << bestFitPoint.ShwEnScale()
       << endl
       << ". Target chi2: " << targetChi2
       << endl;

  for (Int_t direction = -1; direction < 2; direction += 2){

    cout << "Direction = " << direction << endl;

    NuMMParameters scanMMPars = mmPars;
    scanMMPars.Dm2(bestFitDm2);
    scanMMPars.FixDm2();

    Double_t incrementArray[9] = {0.1e-3,
                                  0.03e-3,
                                  0.01e-3,
                                  0.003e-3,
                                  0.001e-3,
                                  0.0003e-3,
                                  0.0001e-3,
                                  0.00003e-3,
                                  0.00001e-3};
    Double_t firstChi2 = bestFitChi2;
    Double_t firstDm2 = bestFitDm2;
    Double_t currentChi2 = -1.0;
    Double_t lastChi2 = -1.0;
    Double_t lastDm2 = -1.0;
    for (Int_t incCount = 0; incCount<9; ++incCount){
      Double_t increment = incrementArray[incCount];

      currentChi2 = firstChi2;
      scanMMPars.Dm2(firstDm2);

      cout << "Starting while with firstChi2 = " << firstChi2
           << ", firstDm2 = " << firstDm2
           << ", increment = " << increment
           << endl;

      while (currentChi2 < targetChi2){
        lastChi2 = currentChi2;
        lastDm2 = scanMMPars.Dm2();

        scanMMPars.Dm2(scanMMPars.Dm2() + direction*increment);
        cout << "Trying dm2 = " << scanMMPars.Dm2() << "... ";
        flush(cout);
        if (scanMMPars.AreAllParametersFixed()){
          currentChi2 = (*this)(scanMMPars.VectorParameters());
        }
        else{
          currentChi2 = (*this)(this->Minimise(scanMMPars).VectorParameters());
        }
        cout << "chi2 = " << currentChi2 << endl;
      }

      firstDm2 = lastDm2;
      firstChi2 = lastChi2;

    }

    Double_t chi2Difference = lastChi2 - currentChi2;
    Double_t dm2Difference = lastDm2 - scanMMPars.Dm2();
    Double_t value = -1.0;
    cout << "chi2Difference = " << chi2Difference << endl;
    cout << "dm2Difference = " << dm2Difference << endl;
    if ((fabs(chi2Difference)<1e-12) || (fabs(dm2Difference)<1e-12)){
      cout << "Chi2 or Dm2 difference is zero. Supplying currentChi2."
           << endl;
      value = scanMMPars.Dm2();
    }
    else{
      Double_t fractionBetween = (targetChi2 - currentChi2)/chi2Difference;
      value = fractionBetween*dm2Difference + scanMMPars.Dm2();
    }
    cout << "Error is between " << lastDm2
         << " and " << scanMMPars.Dm2()
         << ", with chi2 between " << lastChi2
         << " and " << currentChi2
         << endl;
    cout << "Value is " << value << endl;
    if (-1 == direction){errors.first = value;}
    if (1 == direction){errors.second = value;}
  }

  cout << "Final answer is dm2 = " << bestFitDm2
       << " + " << errors.second - bestFitDm2
       << " - " << errors.first - bestFitDm2
       << endl;
  return errors;
}

const NuMMParameters NuSystFitter::FCMinimise

(

const NuMMParameters &

mmPars

)

[virtual]

More robustly minimizes by doing a wide grid search first - same as the FC fitting.

Definition at line 1567 of file NuSystFitter.cxx.

References CoarseGridSearch(), NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), FineGridSearch(), Msg::kInfo, Minimise(), NuMMParameters::MinuitPass(), MSG, NuMMParameters::Sn2(), and NuMMParameters::Sn2Bar().

Referenced by NuFCFitterNuMuBar::Fit().

{
  // Start with a coarse grid search, to seed minuit
  const NuMMParameters mmStart = CoarseGridSearch(mmPars);

  // Now run the normal minuit procedure with this as the starting point
  NuMMParameters mmFit = Minimise(mmStart);

  // If minuit failed, do a finer grid search around the coarse starting point
  // Did minuit succeed?
  // if (mmFit.MinuitPass())
  // {

  // } else {
  if (!mmFit.MinuitPass())
  {
      // Minuit failed to fit. Use a more intensive grid search,
      // starting on the minuit starting point (that we get from
      // a grid search)
      MSG("NuSystFitter", Msg::kInfo) << "Minuit fit failed. Doing Finer Grid search" << endl;
      mmFit = FineGridSearch(mmStart);
  } else {
      // Tell them what minuit found
      MSG("NuSystFitter", Msg::kInfo)
        << "Minuit found: dm2: " << mmFit.Dm2()
        << " sn2: " << mmFit.Sn2()
        << " dm2bar: " << mmFit.Dm2Bar()
        << " sn2bar: " << mmFit.Sn2Bar()
        << '\n';
  }
  return mmFit;
}

const NuMMParameters NuSystFitter::FCMinimiseNSI

(

const NuMMParameters &

mmPars

)

[virtual]

Definition at line 1886 of file NuSystFitter.cxx.

References CoarseGridSearchNSI(), NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), FineGridSearchNSI(), Msg::kInfo, Minimise(), NuMMParameters::MinuitPass(), MSG, NuMMParameters::PrintStatus(), NuMMParameters::Sn2(), and NuMMParameters::Sn2Bar().

{
  // Start with a coarse grid search, to seed minuit
  mmPars.PrintStatus();
  const NuMMParameters mmStart = CoarseGridSearchNSI(mmPars);

  // Now run the normal minuit procedure with this as the starting point
  mmStart.PrintStatus();
  NuMMParameters mmFit = Minimise(mmStart);

  // Test seeding directly with no coarse search
  //NuMMParameters mmFit = Minimise(mmPars);

  // If minuit failed, do a finer grid search around the coarse starting point
  // Did minuit succeed?
  // if (mmFit.MinuitPass())
  // {

  // } else {
  if (!mmFit.MinuitPass())
  {
      // Minuit failed to fit. Use a more intensive grid search,
      // starting on the minuit starting point (that we get from
      // a grid search)
      MSG("NuSystFitter", Msg::kInfo) << "Minuit fit failed. Doing Finer Grid search" << endl;
      mmFit = FineGridSearchNSI(mmStart);
  } else {
      // Tell them what minuit found
      MSG("NuSystFitter", Msg::kInfo)
        << "Minuit found: dm2: " << mmFit.Dm2()
        << " sn2: " << mmFit.Sn2()
        << " dm2bar: " << mmFit.Dm2Bar()
        << " sn2bar: " << mmFit.Sn2Bar()
        << "epsilon: " << mmFit.Epsilon()
        << '\n';
  }
  return mmFit;
}

const NuMMParameters NuSystFitter::FCMinimiseNu

(

const NuMMParameters &

mmPars

)

[virtual]

Definition at line 1727 of file NuSystFitter.cxx.

References CoarseGridSearchNu(), NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), FineGridSearchNu(), Msg::kInfo, Minimise(), NuMMParameters::MinuitPass(), MSG, NuMMParameters::Sn2(), and NuMMParameters::Sn2Bar().

{
  // Start with a coarse grid search, to seed minuit
  const NuMMParameters mmStart = CoarseGridSearchNu(mmPars);    
  
  // Now run the normal minuit procedure with this as the starting point
  NuMMParameters mmFit = Minimise(mmStart);
  
  // If minuit failed, do a finer grid search around the coarse starting point
  // Did minuit succeed?
  // if (mmFit.MinuitPass())
  // {
  
  // } else {
  if (!mmFit.MinuitPass())
  {
    // Minuit failed to fit. Use a more intensive grid search,
    // starting on the minuit starting point (that we get from
    // a grid search)
    MSG("NuSystFitter", Msg::kInfo) << "Minuit fit failed. Doing Finer Grid search" << endl;
    mmFit = FineGridSearchNu(mmStart);
  } else {
    // Tell them what minuit found
    MSG("NuSystFitter", Msg::kInfo)
    << "Minuit found: dm2: " << mmFit.Dm2()
    << " sn2: " << mmFit.Sn2()
    << " dm2bar: " << mmFit.Dm2Bar()
    << " sn2bar: " << mmFit.Sn2Bar()
    << '\n';
  }
  return mmFit;
}

NuMMParameters NuSystFitter::FillAlphaBeta	(	TH2D &	J,
		TH2D &	K,
		const std::string	outfilename
	)

Definition at line 2381 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), minimum, NuUtilities::ProgressBar(), NuMMParameters::Sn2(), and NuMMParameters::Sn2Bar().

{
  // characteristics of each surface
  Int_t NpointsJX = J.GetNbinsX();
  Int_t NpointsJY = J.GetNbinsY();
  Int_t NpointsKX = K.GetNbinsX();
  Int_t NpointsKY = K.GetNbinsY();

  Double_t widthX = J.GetXaxis()->GetBinWidth(1);
  Double_t widthY = J.GetYaxis()->GetBinWidth(1);

  Double_t minKX = K.GetXaxis()->GetBinCenter(1);
  Double_t maxKX = K.GetXaxis()->GetBinCenter(NpointsKX);
  Double_t minKY = K.GetYaxis()->GetBinCenter(1);
  Double_t maxKY = K.GetYaxis()->GetBinCenter(NpointsKY);

  Double_t minJX = J.GetXaxis()->GetBinCenter(1);
  Double_t maxJX = J.GetXaxis()->GetBinCenter(NpointsJX);
  Double_t minJY = J.GetYaxis()->GetBinCenter(1);
  Double_t maxJY = J.GetYaxis()->GetBinCenter(NpointsJY);

  Double_t minABX = minKX - maxJX + widthX/2;
  Double_t maxABX = maxKX - minJX - widthX/2.;
  Double_t minABY = minKY - maxJY + widthY/2.;
  Double_t maxABY = maxKY - minJY - widthY/2.;

  Int_t NpointsAlpha = (int)((maxABY-minABY)/widthY) + 1;
  Int_t NpointsBeta  = (int)((maxABX-minABX)/widthX) + 1;

  // resulting (alpha, beta) surfaces
  TH2D *absurf = new TH2D("ab", "(#alpha,#beta) likelihood surface",
                          NpointsBeta, minABX, maxABX,
                          NpointsAlpha, minABY, maxABY );
  absurf->GetXaxis()->SetTitle("#beta");
  absurf->GetYaxis()->SetTitle("#alpha");

  TH2D *hDmMin = new TH2D("hDmMin", "Dm2 for Minimum Likehood",
                          NpointsBeta, minABX, maxABX,
                          NpointsAlpha, minABY, maxABY );
  hDmMin->GetXaxis()->SetTitle("#beta");
  hDmMin->GetYaxis()->SetTitle("#alpha");

  TH2D *hSnMin = new TH2D("hSnMin", "Sn2 for Minimum Likehood",
                          NpointsBeta, minABX, maxABX,
                          NpointsAlpha, minABY, maxABY );
  hSnMin->GetXaxis()->SetTitle("#beta");
  hSnMin->GetYaxis()->SetTitle("#alpha");

  TH2D *hDmMinBar = new TH2D("hDmMinBar", "Dm2Bar for Minimum Likehood",
                          NpointsBeta, minABX, maxABX,
                          NpointsAlpha, minABY, maxABY );
  hDmMinBar->GetXaxis()->SetTitle("#beta");
  hDmMinBar->GetYaxis()->SetTitle("#alpha");

  TH2D *hSnMinBar = new TH2D("hSnMinBar", "Sn2Bar for Minimum Likehood",
                          NpointsBeta, minABX, maxABX,
                          NpointsAlpha, minABY, maxABY );
  hSnMinBar->GetXaxis()->SetTitle("#beta");
  hSnMinBar->GetYaxis()->SetTitle("#alpha");

  // building the surfaces
  Double_t likeTot, minimum;
  Double_t dmMin=0, snMin=0, dmMinBar=0, snMinBar=0;

  for (Double_t beta=minABX; beta<=maxABX; beta+=widthX){
    Int_t iB = absurf->GetXaxis()->FindFixBin(beta);

    for(Double_t alpha=minABY; alpha<=maxABY; alpha+=widthY){
      Int_t iA = absurf->GetYaxis()->FindFixBin(alpha);
      minimum = -1;
      NuUtilities::ProgressBar(iA*iB, NpointsAlpha*NpointsBeta, 1);

      for(Int_t indexJx=1; indexJx<=NpointsJX; ++indexJx){
        for(Int_t indexJy=1; indexJy<=NpointsJY; ++indexJy){

          Double_t xval = J.GetXaxis()->GetBinCenter(indexJx);
          Double_t yval = J.GetYaxis()->GetBinCenter(indexJy);

          xval += beta;
          if (xval < minKX) continue;
          if (xval > maxKX) continue;

          yval += alpha;
          if (yval < minKY) continue;
          if (yval > maxKY) continue;

          Int_t indexKx = K.GetXaxis()->FindFixBin(xval);
          Int_t indexKy = K.GetYaxis()->FindFixBin(yval);
          likeTot = J.GetBinContent(indexJx, indexJy) + K.GetBinContent(indexKx, indexKy);

          if ( (minimum==-1) || (likeTot<minimum) ) {
            minimum = likeTot;
            dmMin    = J.GetYaxis()->GetBinCenter(indexJy);
            snMin    = J.GetXaxis()->GetBinCenter(indexJx);
            dmMinBar = K.GetYaxis()->GetBinCenter(indexKy);
            snMinBar = K.GetXaxis()->GetBinCenter(indexKx);
          }
        }
      }
      absurf->SetBinContent(iB, iA, minimum);
      hDmMin->SetBinContent(iB, iA, dmMin);
      hSnMin->SetBinContent(iB, iA, snMin);
      hDmMinBar->SetBinContent(iB, iA, dmMinBar);
      hSnMinBar->SetBinContent(iB, iA, snMinBar);
    }
  }

  // shift by minimum
  Int_t minbinX, minbinY, minbinZ;
  absurf->GetMinimumBin(minbinX, minbinY, minbinZ);
  Double_t minContent = absurf->GetBinContent(minbinX, minbinY);

  for (Int_t iterB=1; iterB<=NpointsBeta; ++iterB){
    for (Int_t iterA=1; iterA<=NpointsAlpha; ++iterA){
      Double_t content = absurf->GetBinContent(iterB, iterA);
      absurf->SetBinContent(iterB, iterA, content-minContent);
    }
  }

  // min of Alpha Beta surface
  Double_t minalpha = absurf->GetYaxis()->GetBinCenter(minbinY);
  Double_t minbeta  = absurf->GetXaxis()->GetBinCenter(minbinX);
  TGraph *minG = new TGraph(1, &minbeta, &minalpha);

  // save
  cout << "Plots saved in " << outfilename.c_str() << endl;
  TFile *file = new TFile(outfilename.c_str(), "RECREATE");

  J.Write("J", TObject::kOverwrite);
  K.Write("K", TObject::kOverwrite);
  absurf->Write("AlphaBeta", TObject::kOverwrite);
  hDmMin->Write("dm2Min", TObject::kOverwrite);
  hSnMin->Write("sn2Min", TObject::kOverwrite);
  hDmMinBar->Write("dm2MinBar", TObject::kOverwrite);
  hSnMinBar->Write("sn2MinBar", TObject::kOverwrite);
  minG->Write("AlphaBetaMinimum", TObject::kOverwrite);
  file->Close(0);

  // done!
  NuMMParameters pars;
  pars.Dm2( hDmMin->GetBinContent(minbinX, minbinY) );
  pars.Sn2( hSnMin->GetBinContent(minbinX, minbinY) );
  pars.Dm2Bar( hDmMinBar->GetBinContent(minbinX, minbinY) );
  pars.Sn2Bar( hSnMinBar->GetBinContent(minbinX, minbinY) );
  return pars;
}

void NuSystFitter::FillCPT

(

TH2D &

surface

)

[virtual]

Fills a 2D histogram that has dimensions Dm2 (Y) and DM2Bar(X). Uses the center of the passed bins.

Definition at line 2310 of file NuSystFitter.cxx.

References NuMMParameters::ConstrainPhysical(), MuELoss::e, NuMMParameters::FixDm2(), NuMMParameters::FixDm2Bar(), Msg::kError, Msg::kWarning, Minimise(), minimum, NuMMParameters::MinuitPass(), MSG, NuUtilities::ProgressBar(), NuMMParameters::Sn2(), NuMMParameters::Sn2Bar(), Sn2GridSearch(), and NuMMParameters::VectorParameters().

{
  Double_t minimum = 0.0;

  const Double_t bincount = surf.GetNbinsX() * surf.GetNbinsY();

  // Loop over every bin in the histogram
  for(Int_t xBin = 1; xBin <= surf.GetNbinsX(); ++xBin) {
    for (Int_t yBin = 1; yBin <= surf.GetNbinsY(); ++yBin) {
      // Show a progress bar
      NuUtilities::ProgressBar(surf.GetNbinsY()*(xBin-1) + yBin, bincount, 5);

      const Double_t dm2 =    surf.GetYaxis()->GetBinCenter(yBin);
      const Double_t dm2bar = surf.GetXaxis()->GetBinCenter(xBin);

      // Find the sn2 minimum at this value. First do a grid search
      NuMMParameters pars = Sn2GridSearch(dm2, dm2bar);
      // Get the likelihood, so we can be certain MINUIT improved the fit
      const Double_t likeSparse = (*this)(pars.VectorParameters());

      // Now, minimise from this found point
      pars.FixDm2();
      pars.FixDm2Bar();
      pars.ConstrainPhysical();
      NuMMParameters parsFit = Minimise(pars);

      // Handle MINUIT failure
      if (!parsFit.MinuitPass()) {
        MSG("NuSystFitter",Msg::kWarning)
          << "MINUIT fit did not pass when attempting to make Dm2-Dm2Bar plot. " << '\n'
          << "Dm2: " << dm2 << ", Dm2Bar: " << dm2bar << '\n'
          << "Grid Search returned Sn2: " << pars.Sn2() << ", Sn2Bar: " << pars.Sn2Bar() << '\n'
          << "Doing Higher resolution grid search" << endl;
        // Just do a higher resolution grid search
        parsFit = Sn2GridSearch(dm2, dm2bar, 31);
        MSG("NuSystFitter",Msg::kWarning)
          << "  Found Sn2: " << parsFit.Sn2() << ", Sn2Bar: " << parsFit.Sn2Bar() << endl;
      }

      // We now have a best fit for this point. Get the likelihood
      const Double_t like = (*this)(parsFit.VectorParameters());
      // Sort out the minimums
      if (xBin == 1 && yBin == 1) minimum = like;
      if (like < minimum) minimum = like;

      // Validate this found minimum was actually better than the grid search
      if (like - likeSparse > 1e-4) {
        MSG("NuSystFitter",Msg::kError)
          << "Dm2: " << dm2 << ", Dm2Bar: " << dm2bar << '\n'
          << "Grid Search  Sn2: " << pars.Sn2() << ", Sn2Bar: " << pars.Sn2Bar() << " = " << likeSparse << '\n'
          << "MINUIT Found Sn2: " << parsFit.Sn2() << ", Sn2Bar: " << parsFit.Sn2Bar() << " = " << like << '\n'
          << "MINUIT fit was worse than seeding grid search!\n"  << '\n'
          << "Since this error is not understood, Skipping this bin." << endl;
        continue;
      }
      // Finally, fill the histogram with this likelihood
      surf.SetBinContent(xBin, yBin, like);
    }
  }

  // Now, shift the whole histogram based on the minimum
  for(Int_t xBin = 1; xBin <= surf.GetNbinsX(); ++xBin) {
    for (Int_t yBin = 1; yBin <= surf.GetNbinsY(); ++yBin) {
      Double_t val = surf.GetBinContent(xBin, yBin);
      surf.SetBinContent(xBin, yBin, val-minimum);
    }
  }
}

Double_t NuSystFitter::FillLikelihoodSurface	(	TH2D &	surface,
		NuMMParameters	mmPars
	)			[virtual]

Fills a histogram with likelihood values. This function accepts an arbitrarily binned histogram, and fills the values according to the likelihood. It returns the value of the minimum likelihood. (This does the CC neutirno spectrum).

Definition at line 1323 of file NuSystFitter.cxx.

References NuMMParameters::ContourForNuMus(), and FillLikelihoodSurfaceGeneric().

Referenced by GridSearchNu().

{
  mmPars.ContourForNuMus();
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

Double_t NuSystFitter::FillLikelihoodSurfaceBar	(	TH2D &	surface,
		NuMMParameters	mmPars
	)			[virtual]

Fills a histogram with likelihood values. This function accepts an arbitrarily binned histogram, and fills the values according to the likelihood. It returns the value of the minimum likelihood.

Definition at line 1332 of file NuSystFitter.cxx.

References NuMMParameters::ContourForNuBars(), and FillLikelihoodSurfaceGeneric().

Referenced by NuFCFitter::Archive(), GridSearch(), NuFCFitter::GridSearch(), LikelihoodSurfaceBar(), and NuFCDelegateArchiver::UseFitter().

{
  mmPars.ContourForNuBars();
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

Double_t NuSystFitter::FillLikelihoodSurfaceCPT	(	TH2D &	surface,
		NuMMParameters	mmPars
	)			[virtual]

Definition at line 1525 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), fQuietMode, minimum, NuUtilities::ProgressBar(), NuMMParameters::Sn2(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

{
  Double_t minimum = -1;
  const int Nx = surface.GetNbinsX();
  const int Ny = surface.GetNbinsY();

  for (Int_t x = 1; x <= Nx; x++)
  {
      for (Int_t y = 1; y <= Ny; y++)
      {
          mmPars.Sn2(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2(surface.GetYaxis()->GetBinCenter(y));
          mmPars.Sn2Bar(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2Bar(surface.GetYaxis()->GetBinCenter(y));
          Double_t likelihood = (*this)(mmPars.VectorParameters());

          surface.SetBinContent(x, y, likelihood);

          // Work out if this is the minimum
          if (minimum < 0 || minimum > likelihood) minimum = likelihood;

          // Now print the progress bar unless in "batch mode"
          if(fQuietMode < 2) NuUtilities::ProgressBar((x-1)*Ny+y, Nx*Ny, 1);
      }
  }

  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceDecay	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair = 1
	)			[virtual]

Definition at line 1495 of file NuSystFitter.cxx.

References NuMMParameters::ContourForDecay(), and FillLikelihoodSurfaceGeneric().

{
  mmPars.ContourForDecay(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

Double_t NuSystFitter::FillLikelihoodSurfaceDecoherence	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair = 1
	)			[virtual]

Definition at line 1505 of file NuSystFitter.cxx.

References NuMMParameters::ContourForDecoherence(), and FillLikelihoodSurfaceGeneric().

{
  mmPars.ContourForDecoherence(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

Double_t NuSystFitter::FillLikelihoodSurfaceDecoherenceBar	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair = 1
	)			[virtual]

Definition at line 1515 of file NuSystFitter.cxx.

References NuMMParameters::ContourForDecoherenceBar(), and FillLikelihoodSurfaceGeneric().

{
  mmPars.ContourForDecoherenceBar(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

Double_t NuSystFitter::FillLikelihoodSurfaceGeneric	(	TH3D &	surface,
		NuMMParameters	mmPars
	)			[protected, virtual]

Definition at line 1258 of file NuSystFitter.cxx.

References NuMMParameters::ContourPars3D(), NuMMParameters::FixParameterByIndex(), fQuietMode, Minimise(), minimum, NuUtilities::ProgressBar(), NuMMParameters::SetParameterByIndex(), and NuMMParameters::VectorParameters().

{
  const int Nx = surface.GetNbinsX();
  const int Ny = surface.GetNbinsY();
  const int Nz = surface.GetNbinsZ();

  TAxis* xax = surface.GetXaxis();
  TAxis* yax = surface.GetYaxis();
  TAxis* zax = surface.GetZaxis();

  //cout << "Before ContourPars3D()" << endl;

  const int idx1 = mmPars.ContourPars3D()[0];
  const int idx2 = mmPars.ContourPars3D()[1];
  const int idx3 = mmPars.ContourPars3D()[2];

  //cout << "After ContourPars3D(): idx1= " << idx1
  //     << ", idx2= " << idx2 << ", idx3= " << idx3 << endl;

  Double_t minimum = -1;

  // Order loops so same-dmsq go sequentially for usual osc contours. This
  // takes full advantage of the cacheing in the oscillation functions
  for(int z = 1; z <= Nz; z++){
    for(int y = 1; y <= Ny; y++){
      for(int x = 1; x <= Nx; x++){
        // Now print the progress bar unless in "batch mode"
        // Need this new mode to avoid having nothing but
        // progress bars in FC log files
        if(fQuietMode < 2){
          NuUtilities::ProgressBar((z-1)*Ny*Nx+y+x-1, Nx*Ny*Nz, 1);
        }
        //cout << "Before mmPars.SetParameterByIndex()" << endl;
        //cout << "(x, y, z): " << x << ", " << y << ", " << z << endl;
        //cout << "xax->GetBinCenter(x): " << xax->GetBinCenter(x)
        //     << ", yax->GetBinCenter(y): " << yax->GetBinCenter(y)
        //     << ", zax->GetBinCenter(z): " << zax->GetBinCenter(z) << endl;
        mmPars.SetParameterByIndex(idx3, xax->GetBinCenter(x));
        //cout << "Set parameter idx3: " << idx3 << endl;
        mmPars.FixParameterByIndex(idx3);
        mmPars.SetParameterByIndex(idx2, zax->GetBinCenter(z));
        //cout << "Set parameter idx2: " << idx2 << endl;
        mmPars.FixParameterByIndex(idx2);
        mmPars.SetParameterByIndex(idx1, yax->GetBinCenter(y));
        //cout << "Set parameter idx1: " << idx1 << endl;
        mmPars.FixParameterByIndex(idx1);
        //cout << "In three for loops about to minimize mmPars" << endl;

        mmPars = Minimise(mmPars);

        const Double_t likelihood = (*this)(mmPars.VectorParameters());
        surface.SetBinContent(x, y, z, likelihood);

        // Work out if this is the minimum
        if (minimum < 0 || minimum > likelihood) minimum = likelihood;
      }
    }
  }
  cout << "minimum is: " << minimum << endl;
  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceGeneric	(	TH2D &	surface,
		NuMMParameters	mmPars
	)			[protected, virtual]

Helper function. Fills a likelihood surface based on the parameters as defined by NuMMParameters::ContourPars()

Definition at line 1207 of file NuSystFitter.cxx.

References NuMMParameters::ContourPars(), NuMMParameters::FixParameterByIndex(), fQuietMode, Minimise(), minimum, NuUtilities::ProgressBar(), NuMMParameters::SetParameterByIndex(), and NuMMParameters::VectorParameters().

Referenced by FillLikelihoodSurface(), FillLikelihoodSurfaceBar(), FillLikelihoodSurfaceDecay(), FillLikelihoodSurfaceDecoherence(), FillLikelihoodSurfaceDecoherenceBar(), FillLikelihoodSurfaceNSI(), FillLikelihoodSurfaceNSI3D(), FillLikelihoodSurfaceNSIBar(), and FillLikelihoodSurfaceSterile().

{

  NuMMParameters surfPars;

  const int Nx = surface.GetNbinsX();
  const int Ny = surface.GetNbinsY();

  TAxis* xax = surface.GetXaxis();
  TAxis* yax = surface.GetYaxis();

  const int idx1 = mmPars.ContourPars().first;
  const int idx2 = mmPars.ContourPars().second;

  Double_t minimum = -1;

  // Order loops so same-dmsq go sequentially for usual osc contours. This
  // takes full advantage of the cacheing in the oscillation functions
  for(int y = 1; y <= Ny; y++){
    for(int x = 1; x <= Nx; x++){
      // Now print the progress bar unless in "batch mode"
      // Need this new mode to avoid having nothing but
      // progress bars in FC log files
      if(fQuietMode < 2){
        NuUtilities::ProgressBar((y-1)*Nx+x, Nx*Ny, 1);
      }

      surfPars = mmPars;

      surfPars.SetParameterByIndex(idx2, xax->GetBinCenter(x));
      surfPars.FixParameterByIndex(idx2);
      surfPars.SetParameterByIndex(idx1, yax->GetBinCenter(y));
      surfPars.FixParameterByIndex(idx1);

      surfPars = Minimise(surfPars);

      const Double_t likelihood = (*this)(surfPars.VectorParameters());

      surface.SetBinContent(x, y, likelihood);

      // Work out if this is the minimum
      if (minimum < 0 || minimum > likelihood) minimum = likelihood;
    }
  }

  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceNSI	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair
	)			[virtual]

Definition at line 1339 of file NuSystFitter.cxx.

References NuMMParameters::ContourForNSI(), and FillLikelihoodSurfaceGeneric().

Referenced by NuFCFitterNSINu::Archive(), NuFCFitterNSINu::GridSearch(), and GridSearchNSI().

{
  mmPars.ContourForNSI(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);

}

Double_t NuSystFitter::FillLikelihoodSurfaceNSI3D	(	TH3D &	surface,
		NuMMParameters	mmPars
	)			[virtual]

Definition at line 1359 of file NuSystFitter.cxx.

References NuMMParameters::ContourForNSI(), and FillLikelihoodSurfaceGeneric().

Referenced by NuFCFitterNSI::Archive(), and NuFCFitterNSI::GridSearch().

{
  mmPars.ContourForNSI();
  //  cout << "Set contours for NSI: ContourPars3D()[0]= " << mmPars.ContourPars3D()[0]
  //     << ", ContourPars3D()[1]= " << mmPars.ContourPars3D()[1]
  //     << ", ContourPars3D()[2]= " << mmPars.ContourPars3D()[2]
  //     << endl;
  return FillLikelihoodSurfaceGeneric(surface, mmPars);

}

Double_t NuSystFitter::FillLikelihoodSurfaceNSIBar	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair
	)			[virtual]

Definition at line 1349 of file NuSystFitter.cxx.

References NuMMParameters::ContourForNSIBar(), and FillLikelihoodSurfaceGeneric().

Referenced by NuFCFitterNSINubar::Archive(), and NuFCFitterNSINubar::GridSearch().

{
  mmPars.ContourForNSIBar(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);

}

Double_t NuSystFitter::FillLikelihoodSurfaceNSIDm2Eps	(	TH2D &	surface,
		NuMMParameters	mmPars,
		TH2D &	sinvals
	)			[virtual]

Definition at line 1374 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Epsilon(), NuMMParameters::FixDm2(), NuMMParameters::FixEpsilon(), fQuietMode, Likelihood(), Minimise(), minimum, NuUtilities::ProgressBar(), and NuMMParameters::Sn2().

{
 
  Int_t Neps = surface.GetNbinsX();
  Int_t Ndm2 = surface.GetNbinsY();

  Double_t minimum = 1e6;

  for (int i = 1; i <= Neps; i++){
    for (int j = 1; j <= Ndm2; j++){

      if(fQuietMode < 2){
        NuUtilities::ProgressBar((i-1)*Ndm2+j, Ndm2*Neps, 1);
      }

      Double_t eps = surface.GetXaxis()->GetBinCenter(i);
      Double_t dm2 = surface.GetYaxis()->GetBinCenter(j);
      mmPars.Dm2(dm2);
      mmPars.Epsilon(eps);
      mmPars.FixDm2();
      mmPars.FixEpsilon();

      NuMMParameters bf = this->Minimise(mmPars);
      Double_t likelihood = Likelihood(bf);
      Double_t sn2Min = bf.Sn2();

      surface.SetBinContent(i, j, likelihood);
      sinvals.SetBinContent(i, j, sn2Min);
      if(likelihood < minimum) minimum = likelihood;      
    }
  }

  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceNSIDm2Sn2	(	TH2D &	surface,
		NuMMParameters	mmPars,
		TH2D &	epsvals
	)			[virtual]

Definition at line 1454 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Epsilon(), NuMMParameters::FixDm2(), NuMMParameters::FixSn2(), fQuietMode, Likelihood(), Minimise(), minimum, NuUtilities::ProgressBar(), and NuMMParameters::Sn2().

{
 
  Int_t Nsn2 = surface.GetNbinsX();
  Int_t Ndm2 = surface.GetNbinsY();

  Double_t minimum = 1e6;

  for (int i = 1; i <= Nsn2; i++){
    for (int j = 1; j <= Ndm2; j++){

      if(fQuietMode < 2){
        NuUtilities::ProgressBar((i-1)*Ndm2+j, Ndm2*Nsn2, 1);
      }

      Double_t sn2 = surface.GetXaxis()->GetBinCenter(i);
      Double_t dm2 = surface.GetYaxis()->GetBinCenter(j);
      mmPars.Dm2(dm2);
      mmPars.Sn2(sn2);
      mmPars.FixDm2();
      mmPars.FixSn2();

      NuMMParameters bf = this->Minimise(mmPars);
      Double_t likelihood = Likelihood(bf);
      Double_t epsMin = bf.Epsilon();
      
      surface.SetBinContent(i, j, likelihood);
      epsvals.SetBinContent(i, j, epsMin);

      if(likelihood < minimum) minimum = likelihood;      
    }
  }

  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceNSISn2Eps	(	TH2D &	surface,
		NuMMParameters	mmPars,
		TH2D &	massvals
	)			[virtual]

Definition at line 1413 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Epsilon(), NuMMParameters::FixEpsilon(), NuMMParameters::FixSn2(), fQuietMode, Likelihood(), Minimise(), minimum, NuUtilities::ProgressBar(), and NuMMParameters::Sn2().

{
 
  Int_t Neps = surface.GetNbinsX();
  Int_t Nsn2 = surface.GetNbinsY();

  Double_t minimum = 1e6;

  for (int i = 1; i <= Neps; i++){
    for (int j = 1; j <= Nsn2; j++){

      if(fQuietMode < 2){
        NuUtilities::ProgressBar((i-1)*Nsn2+j, Nsn2*Neps, 1);
      }

      Double_t eps = surface.GetXaxis()->GetBinCenter(i);
      Double_t sn2 = surface.GetYaxis()->GetBinCenter(j);
      mmPars.Sn2(sn2);
      mmPars.Epsilon(eps);
      mmPars.FixSn2();
      mmPars.FixEpsilon();

      NuMMParameters bf = this->Minimise(mmPars);
      Double_t likelihood = Likelihood(bf);
      Double_t dm2Min = bf.Dm2();

      surface.SetBinContent(i, j, likelihood);
      massvals.SetBinContent(i, j, dm2Min);

      if(likelihood < minimum) minimum = likelihood;      
    }
  }

  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceSterile	(	TH2D &	surface,
		NuMMParameters	mmPars,
		Int_t	ContourPair
	)			[virtual]

Definition at line 1555 of file NuSystFitter.cxx.

References NuMMParameters::ContourForSterile(), and FillLikelihoodSurfaceGeneric().

                                                                       {
  mmPars.ContourForSterile(ContourPair);
  return FillLikelihoodSurfaceGeneric(surface, mmPars);
}

NuMMParameters NuSystFitter::FineGridSearch

(

NuMMParameters

mmStart

)

[protected]

Definition at line 1634 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), MuELoss::e, GridSearch(), Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerSn2BarLimit(), MSG, and NuMMParameters::Sn2Bar().

Referenced by FCMinimise().

{
  // Calculate the location of the fine grid search
  Double_t gridL, gridR, gridT, gridB;
  gridL = gridR = gridT = gridB = -1.0;

  // Go +- 0.2 in sin space
  gridL = mmStart.Sn2Bar() - 0.2;
  gridR = gridL + 0.4;
  // Now limit it to the physical realm
  if (gridL <= mmStart.LowerSn2BarLimit()) gridL = mmStart.LowerSn2BarLimit();
  if (gridR > 1.0) gridR = 1.0;

  // And go +-2e-3 in mass space
  gridB = mmStart.Dm2Bar() - 2e-3;
  gridT = gridB + 4e-3;
  // We only really need to worry about lower limit here
  if (gridB <= mmStart.LowerDm2BarLimit()) gridB = mmStart.LowerSn2BarLimit();


  MSG("NuSystFitter", Msg::kInfo)
    << "Doing fine grid search on sin = [ "
    << gridL << ", " << gridR << " ]\n"
    << "                              dm  = [ "
    << gridB << ", " << gridT << " ]\n";

  // Ensure that we haven't done anything stupid like set low above high
  if (gridL >= gridR || gridT <= gridB) {
//      cout << "Error: Negative range in fine grid search" << endl;
      throw runtime_error("Negative range in fine grid search");
  }

  // Make the surface

  // The width and height of the grid
  const Int_t gridWidth = 19;
  const Int_t gridHeight = 19;

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
    gridHeight, gridB, gridT);

  // And now do a grid search on this
  NuMMParameters mmGrid = GridSearch(likesurf, mmStart);

  // Output the minimum
  MSG("NusystFitter", Msg::kInfo) << "Fine grid search found sn2bar: "
    << mmGrid.Sn2Bar() << ", dm2bar: " << mmGrid.Dm2Bar() << '\n';

  return mmGrid;
}

NuMMParameters NuSystFitter::FineGridSearchNSI

(

NuMMParameters

mmStart

)

[protected]

Definition at line 1954 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), MuELoss::e, NuMMParameters::Epsilon(), GridSearchNSI(), Msg::kInfo, NuMMParameters::LowerDm2Limit(), NuMMParameters::LowerEpsilonLimit(), and MSG.

Referenced by FCMinimiseNSI().

{
  // Calculate the location of the fine grid search
  Double_t gridL, gridR, gridT, gridB;
  gridL = gridR = gridT = gridB = -1.0;

//   // Go +- 0.2 in sin space
//   gridL = mmStart.Sn2Bar() - 0.2;
//   gridR = gridL + 0.4;
//   // Now limit it to the physical realm
//   if (gridL <= mmStart.LowerSn2BarLimit()) gridL = mmStart.LowerSn2BarLimit();
//   if (gridR > 1.0) gridR = 1.0;

  // And go +-2e-3 in mass space
  gridB = mmStart.Dm2() - 2e-3;
  gridT = gridB + 4e-3;
  // We only really need to worry about lower limit here
  if (gridB <= mmStart.LowerDm2Limit()) gridB = mmStart.LowerDm2Limit();

  // Go +- 0.2 in epsilon space
  gridL = mmStart.Epsilon() - 0.2;
  gridR = gridL + 0.4;
  // Now limit it to the physical realm
  if (gridL <= mmStart.LowerEpsilonLimit()) gridL = mmStart.LowerEpsilonLimit();

  MSG("NuSystFitter", Msg::kInfo)
    << "Doing fine grid search on eps = [ "
    << gridL << ", " << gridR << " ]\n"
    << "                              dm  = [ "
    << gridB << ", " << gridT << " ]\n";

  // Ensure that we haven't done anything stupid like set low above high
  if (gridL >= gridR || gridT <= gridB) {
//      cout << "Error: Negative range in fine grid search" << endl;
      throw runtime_error("Negative range in fine grid search");
  }

  // Make the surface

  // The width and height of the grid
  const Int_t gridWidth = 19;
  const Int_t gridHeight = 19;

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
    gridHeight, gridB, gridT);

  // And now do a grid search on this
  NuMMParameters mmGrid = GridSearchNSI(likesurf, mmStart);

  // Output the minimum
  MSG("NusystFitter", Msg::kInfo) << "Fine grid search found epsilon: "
    << mmGrid.Epsilon() << ", dm2: " << mmGrid.Dm2() << '\n';

  return mmGrid;
}

NuMMParameters NuSystFitter::FineGridSearchNu

(

NuMMParameters

mmStart

)

[protected]

Definition at line 1794 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), MuELoss::e, GridSearchNu(), Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerSn2BarLimit(), MSG, NuMMParameters::Sn2(), and NuMMParameters::Sn2Bar().

Referenced by FCMinimiseNu().

{
  // Calculate the location of the fine grid search
  Double_t gridL, gridR, gridT, gridB;
  gridL = gridR = gridT = gridB = -1.0;
  
  // Go +- 0.2 in sin space
  gridL = mmStart.Sn2Bar() - 0.2;
  gridR = gridL + 0.4;
  // Now limit it to the physical realm
  if (gridL <= mmStart.LowerSn2BarLimit()) gridL = mmStart.LowerSn2BarLimit();
  if (gridR > 1.0) gridR = 1.0;
  
  // And go +-2e-3 in mass space
  gridB = mmStart.Dm2Bar() - 2e-3;
  gridT = gridB + 4e-3;
  // We only really need to worry about lower limit here
  if (gridB <= mmStart.LowerDm2BarLimit()) gridB = mmStart.LowerSn2BarLimit();
  
  
  MSG("NuSystFitter", Msg::kInfo) 
  << "Doing fine grid search on sin = [ "
  << gridL << ", " << gridR << " ]\n"
  << "                              dm  = [ "
  << gridB << ", " << gridT << " ]\n";
  
  // Ensure that we haven't done anything stupid like set low above high
  if (gridL >= gridR || gridT <= gridB) {
    //      cout << "Error: Negative range in fine grid search" << endl;
    throw runtime_error("Negative range in fine grid search");
  }
  
  // Make the surface
  
  // The width and height of the grid
  const Int_t gridWidth = 19;
  const Int_t gridHeight = 19;
  
  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
                gridHeight, gridB, gridT);
  
  // And now do a grid search on this
  NuMMParameters mmGrid = GridSearchNu(likesurf, mmStart);
  
  // Output the minimum
  MSG("NusystFitter", Msg::kInfo) << "Fine grid search found sn2: "
  << mmGrid.Sn2() << ", dm2: " << mmGrid.Dm2() << '\n';
  
  return mmGrid;
}

NuMMParameters NuSystFitter::GridSearch	(	TH2D &	likesurf,
		NuMMParameters	mmGrid,
		NuMMParameters *	Ref = 0
	)			[protected]

Definition at line 1687 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), FillLikelihoodSurfaceBar(), Msg::kDebug, MAXMSG, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by CoarseGridSearch(), and FineGridSearch().

{
    Double_t likelihood = FillLikelihoodSurfaceBar(likesurf, mmGrid);

    // Find the point on this surface that is minimum
    Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
    likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
    Double_t minsin = likesurf.GetXaxis()->GetBinCenter(minbinX);
    Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);

    // Check the zero point, to see if it is lower than anything else on the grid
    mmGrid.Sn2Bar(0);
    mmGrid.Dm2Bar(0);
    Double_t zeropoint = (*this)(mmGrid.VectorParameters());
    if (zeropoint < likelihood) {
        // Then the no-oscillation case is the minimum
        minsin = 0.0;
        minmass = 0.0;
    }

    // If we have been given a reference point, check to see if this is lower
    if (Ref) {
      Double_t refpoint = (*this)(Ref->VectorParameters());
      if (refpoint < likelihood) {
        MAXMSG("NuSystFitter", 5, Msg::kDebug) << "Reference point is lower than grid search" << endl;
        minsin = Ref->Sn2Bar();
        minmass = Ref->Dm2Bar();
      }
    }

    // Set mmPars to the newly found point, then return it
    mmGrid.Sn2Bar(minsin);
    mmGrid.Dm2Bar(minmass);

    return mmGrid;
}

NuMMParameters NuSystFitter::GridSearchNSI	(	TH2D &	likesurf,
		NuMMParameters	mmGrid,
		NuMMParameters *	Ref = 0
	)			[protected]

Definition at line 2012 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Epsilon(), FillLikelihoodSurfaceNSI(), Msg::kDebug, MAXMSG, and NuMMParameters::VectorParameters().

Referenced by CoarseGridSearchNSI(), and FineGridSearchNSI().

{
  Double_t likelihood = FillLikelihoodSurfaceNSI(likesurf, mmGrid, 1);

    // Find the point on this surface that is minimum
    Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
    likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
    Double_t mineps = likesurf.GetXaxis()->GetBinCenter(minbinX);
    Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);

    // Check the zero point, to see if it is lower than anything else on the grid
    mmGrid.Epsilon(0);
    mmGrid.Dm2(0);
    Double_t zeropoint = (*this)(mmGrid.VectorParameters());
    if (zeropoint < likelihood) {
        // Then the no-oscillation case is the minimum
        mineps = 0.0;
        minmass = 0.0;
    }

    // If we have been given a reference point, check to see if this is lower
    if (Ref) {
      Double_t refpoint = (*this)(Ref->VectorParameters());
      if (refpoint < likelihood) {
        MAXMSG("NuSystFitter", 5, Msg::kDebug) << "Reference point is lower than grid search" << endl;
        mineps = Ref->Epsilon();
        minmass = Ref->Dm2();
      }
    }

    // Set mmPars to the newly found point, then return it
    mmGrid.Epsilon(mineps);
    mmGrid.Dm2(minmass);

    return mmGrid;
}

NuMMParameters NuSystFitter::GridSearchNu	(	TH2D &	likesurf,
		NuMMParameters	mmGrid,
		NuMMParameters *	Ref = 0
	)			[protected]

Definition at line 1847 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), FillLikelihoodSurface(), Msg::kDebug, MAXMSG, NuMMParameters::Sn2(), and NuMMParameters::VectorParameters().

Referenced by CoarseGridSearchNu(), and FineGridSearchNu().

{
  Double_t likelihood = FillLikelihoodSurface(likesurf, mmGrid);
  
  // Find the point on this surface that is minimum
  Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
  likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
  Double_t minsin = likesurf.GetXaxis()->GetBinCenter(minbinX);
  Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);
  
  // Check the zero point, to see if it is lower than anything else on the grid
  mmGrid.Sn2(0);
  mmGrid.Dm2(0);
  Double_t zeropoint = (*this)(mmGrid.VectorParameters());
  if (zeropoint < likelihood) {
    // Then the no-oscillation case is the minimum
    minsin = 0.0;
    minmass = 0.0;
  }
  
  // If we have been given a reference point, check to see if this is lower
  if (Ref) {
    Double_t refpoint = (*this)(Ref->VectorParameters());
    if (refpoint < likelihood) {
      MAXMSG("NuSystFitter", 5, Msg::kDebug) << "Reference point is lower than grid search" << endl;
      minsin = Ref->Sn2();
      minmass = Ref->Dm2();
    }
  }
  
  // Set mmPars to the newly found point, then return it
  mmGrid.Sn2(minsin);
  mmGrid.Dm2(minmass);
  
  return mmGrid;
}

double NuSystFitter::Likelihood

(

const NuMMParameters &

mmPars

)

const [virtual]

Definition at line 635 of file NuSystFitter.cxx.

References fQuietMode, fRuns, and NuMMParameters::PenaltyTerm().

Referenced by NuABFitter::DeltaChi(), NuABFitter::DoOutput(), NuABFitter::FillLikelihoodSurfaceAB(), FillLikelihoodSurfaceNSIDm2Eps(), FillLikelihoodSurfaceNSIDm2Sn2(), FillLikelihoodSurfaceNSISn2Eps(), NuFCFitterNuMuBar::Fit(), operator()(), NuABFitter::RecoverNegativeDeltaChi(), and NuFCDelegateArchiver::UseFitter().

{
    Double_t likelihood = 0.0;
    for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
         fRuns->end() != runIt;
         ++runIt){
      likelihood += (*runIt)->ComparePredWithData(mmPars);
    }
    if (!fQuietMode){
        cout << "Stats likelihood: " << likelihood;
    }
    likelihood += mmPars.PenaltyTerm();
    if (!fQuietMode){
        cout << " + penalty = " << likelihood << endl;
    }

    return likelihood;
}

TH2D * NuSystFitter::LikelihoodSurfaceBar	(	const NuMMParameters &	mmPars,
		Int_t	sinPoints,
		Int_t	massPoints,
		Double_t	Sn2BarMin = 0.0,
		Double_t	Sn2BarMax = 1.0,
		Double_t	Dm2BarMin = 0.0,
		Double_t	Dm2BarMax = 12e-3,
		Bool_t	True_Minimum = true
	)			[virtual]

Function to generate a likelihood surface.

Definition at line 1094 of file NuSystFitter.cxx.

References NuMMParameters::Dm2Bar(), FillLikelihoodSurfaceBar(), fQuietMode, Minimise(), minimum, NuMMParameters::MinuitPass(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

{
    // Calculate the bin spacing. This is -1 because we want to 'add'
    // one bin outside of the specified range, so that the bin contents
    // are representative of the center of the bin.
    Double_t sininterval =  (sn2BarMax - sn2BarMin) / (sinPoints -1);
    Double_t massinterval = (dm2BarMax - dm2BarMin) / (massPoints - 1);

    // Offsets. This is how much to shift the bins min and max by
    Double_t sinoffset = sininterval / 2.0;
    Double_t massoffset = massinterval / 2.0;
    // Create the destination histogram
    TH2D *surface = new TH2D(   "chisurf", "#chi^{2} Likelihood Surface",
                                sinPoints, sn2BarMin-sinoffset, sn2BarMax+sinoffset,
                                massPoints, dm2BarMin-massoffset, dm2BarMax+massoffset );
    // Set it up
    surface->GetXaxis()->SetTitle("sin^{2}2#bar#theta");
    surface->GetYaxis()->SetTitle("#Delta#barm^{2}");

    // A variable to store the minimum likelihood encountered
    // Double_t minimum = -1;

    // Make a copy of the parameters
    NuMMParameters mmPars(mmParsIn);

    // Variables to store the current value. These are declared here
    // so that we can verify we have counted them correctly after the
    // looping ends.
    // Double_t sinVal = 0., dmVal = 0.;

    if (!fQuietMode) cout << "Generating Likelihood Surface..." << endl;

    Double_t minimum = this->FillLikelihoodSurfaceBar(*surface, mmPars);

    // // Now, start looping over the surface
    // for (Int_t sinBin = 1; sinBin <= sinPoints; sinBin++)
    // {
    //     // Calculate and set this sin (bin 0 starts at -sinoffset but
    //     // we want it to contain the value at 0)
    //     sinVal = sn2BarMin + sininterval*(sinBin-1);
    //     mmPars.Sn2Bar(sinVal);
    //
    //     for (Int_t massBin = 1; massBin <= massPoints; massBin++)
    //     {
    //         // Calculate and set the dmass value
    //         dmVal = dm2BarMin + massinterval*(massBin-1);
    //         mmPars.Dm2Bar(dmVal);
    //
    //         // Now, calculate the likelihood using this fitter
    //         Double_t likelihood = (*this)(mmPars.VectorParameters());
    //
    //         // Minimim calculations
    //         if (minimum < 0 || likelihood < minimum)
    //             minimum = likelihood;
    //
    //         // Set the bin content!
    //         surface->SetBinContent(sinBin, massBin, likelihood);
    //     }
    // }
    // Print a message to say we have finished
    if (!fQuietMode) {
        cout << "Surface Generation Done. Minimum Likelihood = "
             << minimum << endl;
    }


    // Get the minimum bin on this histogram
    Int_t minbinX, minbinY, minbinZ;
    surface->GetMinimumBin(minbinX, minbinY, minbinZ);
    // Set the fit parameters to this point as a start point
    mmPars.Sn2Bar(surface->GetXaxis()->GetBinCenter(minbinX));
    mmPars.Dm2Bar(surface->GetYaxis()->GetBinCenter(minbinY));

  if (True_Minimum) {
    // Now minimise at this point, and grab the likelihood
    if (!fQuietMode) cout << "Searching for true minimum" << endl;

    NuMMParameters mmFit = this->Minimise(mmPars);
    Double_t minlike = (*this)(mmFit.VectorParameters());

    // Check if the fit failed, if it did, then don't use it
    if (mmFit.MinuitPass())
    {
        // If this is smaller than this minimum we have found (it should be)
        if (minlike <= minimum)
          minimum = minlike;
        else
          cout << "WARNING: Fit Minimisation found minimum higher than grid search" << endl;
    }
  }
    // Now, we have generated likelihoods for all grid points.
    // Offset each one by the minimum, so that it is a delta chi
    // squared surface.
    for (Int_t sinBin = 1; sinBin <= sinPoints; sinBin++) {
        for (Int_t massBin = 1; massBin <= massPoints; massBin++) {
            Double_t content = surface->GetBinContent(sinBin, massBin);
            surface->SetBinContent(sinBin, massBin, content - minimum);
        }
    }

    // Now we have generated and offset the histogram. Return it!
    return surface;
}

const NuMMParameters NuSystFitter::Minimise	(	const NuMMParameters &	mmPars,
		bool &	success
	)			[virtual]

Definition at line 99 of file NuSystFitter.cxx.

References NuMMParameters::AreAllParametersFixed(), NuMMParameters::CopyConstraintsLimitsAndFixings(), fQuietMode, NuMMParameters::kAlpha, NuMMParameters::kDelta1, NuMMParameters::kDelta2, NuMMParameters::kDelta3, NuMMParameters::kDm2, NuMMParameters::kDm221, NuMMParameters::kDm243, NuMMParameters::kDm2Bar, NuMMParameters::kEpsilon, NuMMParameters::kMu2, NuMMParameters::kNcBack, NuMMParameters::kNorm, NuMMParameters::kNumParameters, NuMMParameters::kShwScale, NuMMParameters::kSn2, NuMMParameters::kSn2Bar, NuMMParameters::kTheta12, NuMMParameters::kTheta13, NuMMParameters::kTheta14, NuMMParameters::kTheta24, NuMMParameters::kTheta34, NuMMParameters::MinuitParameters(), NuMMParameters::MinuitPass(), and NuMMParameters::PrintStatus().

{
  success = true;

  if(mmPars.AreAllParametersFixed()) return mmPars;

  //Set up Minuit parameters
  ROOT::Minuit2::MnUserParameters mnPars = mmPars.MinuitParameters();
    
  vector<double> param_vec = mnPars.Params();
//    cout << "Fitting with " << param_vec.size() << " parameters." << endl;
//   cout << "Fitting with " << mnPars.VariableParameters() << " variable parameters." << endl;
  if (!fQuietMode) mmPars.PrintStatus();

  //Run Minuit fit
  ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();

  if (!funcMin.IsValid()) {
    success = false;
    cout << "Function minimum is not valid!" << endl;
  } else {
    success = true;
    if (!fQuietMode){
      const ROOT::Minuit2::MnUserParameters upars = funcMin.UserParameters();

      cout << "Function minimum is: " << endl
           << "Sn2: "       << upars.Parameter(NuMMParameters::kSn2).Value()
           << ", dm2: "     << upars.Parameter(NuMMParameters::kDm2).Value()
           << "; norm: "    << upars.Parameter(NuMMParameters::kNorm).Value()
           << "; shwEn: "   << upars.Parameter(NuMMParameters::kShwScale).Value()
           << "; NCBack: "  << upars.Parameter(NuMMParameters::kNcBack).Value()
           << "; Sn2bar: "  << upars.Parameter(NuMMParameters::kSn2Bar).Value()
           << "; dm2bar: "  << upars.Parameter(NuMMParameters::kDm2Bar).Value()
           << "; epsilon: " << upars.Parameter(NuMMParameters::kEpsilon).Value()
           << "; Theta12: " << upars.Parameter(NuMMParameters::kTheta12).Value()
           << "; Theta13: " << upars.Parameter(NuMMParameters::kTheta13).Value()
           << "; Theta14: " << upars.Parameter(NuMMParameters::kTheta14).Value()
           << "; Delta1: "  << upars.Parameter(NuMMParameters::kDelta1).Value()
           << "; Delta2: "  << upars.Parameter(NuMMParameters::kDelta2).Value()
           << "; Delta3: "  << upars.Parameter(NuMMParameters::kDelta3).Value()
           << "; Theta24: " << upars.Parameter(NuMMParameters::kTheta24).Value()
           << "; Theta34: " << upars.Parameter(NuMMParameters::kTheta34).Value()
           << "; Dm221: "   << upars.Parameter(NuMMParameters::kDm221).Value()
           << "; Dm243: "   << upars.Parameter(NuMMParameters::kDm243).Value()
           << "; mu2: "     << upars.Parameter(NuMMParameters::kMu2).Value()
           << "; alpha: "   << upars.Parameter(NuMMParameters::kAlpha).Value()
           << "; CPT: "     << upars.Parameter(NuMMParameters::kNumParameters).Value()
      << endl;
    }
  }

  //cout << "Diagonal covar elements:" << endl;
  //ROOT::Minuit2::MnUserCovariance covar = funcMin.UserCovariance();
  //for (unsigned int i = 0; i < covar.Nrow(); i++) {
  //  cout << "element " << i << " = " << covar(i,i) << endl;
  //}

  NuMMParameters bestFitPars(funcMin);
  // Make sure the parameters we return are configured the same way as what was passed
  bestFitPars.CopyConstraintsLimitsAndFixings(mmPars);
  bestFitPars.MinuitPass(success);
  //bestFitPars.PrintStatus();
  return bestFitPars;
}

const NuMMParameters NuSystFitter::Minimise

(

const NuMMParameters &

mmPars

)

[virtual]

Definition at line 92 of file NuSystFitter.cxx.

Referenced by Chi2ValleyBar(), Chi2ValleyBarDm2Bar(), Chi2ValleyDm2(), Chi2ValleySn2(), FCMinimise(), FCMinimiseNSI(), FCMinimiseNu(), FillCPT(), FillLikelihoodSurfaceGeneric(), FillLikelihoodSurfaceNSIDm2Eps(), FillLikelihoodSurfaceNSIDm2Sn2(), FillLikelihoodSurfaceNSISn2Eps(), NuFCFitterNuMuBar::Fit(), NuEZFitterNSI::Fit(), NuEZFitter::Fit(), LikelihoodSurfaceBar(), NuFCFitterNSINubar::MinuitFit(), NuFCFitterNSINu::MinuitFit(), NuFCFitterNSI::MinuitFit(), NuFCFitter::MinuitFit(), NuFCFitterNSINubar::RecoverNegativeDeltaChi(), NuFCFitterNSINu::RecoverNegativeDeltaChi(), NuFCFitterNSI::RecoverNegativeDeltaChi(), and NuFCFitter::RecoverNegativeDeltaChi().

{
  bool success;
  return Minimise(mmPars, success);
}

const std::vector< std::pair< Double_t, Double_t > > NuSystFitter::NeutrinoContourFinder

(

const NuMMParameters &

mmPars

)

[virtual]

Definition at line 301 of file NuSystFitter.cxx.

References NuMMParameters::FixSn2(), NuMMParameters::Sn2(), and NuMMParameters::VectorParameters().

{
  //Find the target chi2 for the contour
  NuMMParameters bestFitPars = this->Minimise(inputPars);
  Double_t bestFitChi2 = (*this)(bestFitPars.VectorParameters());
  Double_t targetChi2 = bestFitChi2 + this->Up();
  cout << "Target chi2 = " << targetChi2 << endl;

  //Fix sn2 ready for the dm2 scans:
  NuMMParameters mmPars = inputPars;
  mmPars.FixSn2();
  //Make somewhere to put the contour points
  std::vector<pair<Double_t,Double_t> > pointsUp;
  std::vector<pair<Double_t,Double_t> > pointsDown;

  //Run the contour in quiet mode
  this->QuietModeOn();

  //Get the points for the contour
  for (vector<Double_t>::iterator snIt = fsn2PointsForNeutrinoContour.begin();
       snIt != fsn2PointsForNeutrinoContour.end();
       ++snIt){
    Double_t sn2 = *snIt;
    mmPars.Sn2(sn2);
    const std::pair<Double_t, Double_t> values =
      this->DmScanForContour(mmPars,targetChi2);
    cout << "Contour points: sn2 = " << sn2
         << "; dm2 = " << values.first << " and " << values.second
         << endl;
    if ((values.first < 0.0) || (values.second < 0.0)){
      cout << "End of contour at sn2 = " << mmPars.Sn2() << endl;
      break;
    }
    else{
      pair<Double_t, Double_t> pointDown(sn2,values.first);
      pair<Double_t, Double_t> pointUp(sn2,values.second);
      pointsDown.push_back(pointDown);
      pointsUp.push_back(pointUp);
    }
  }

  //Make a single vector of points
  std::vector<pair<Double_t,Double_t> > contour;
  for (vector<pair<Double_t, Double_t> >::const_iterator it = pointsUp.begin();
       it != pointsUp.end();
       ++it){
    pair<Double_t,Double_t> currentPoint((*it).first,(*it).second);
    contour.push_back(currentPoint);
  }
  for (vector<pair<Double_t, Double_t> >::reverse_iterator
         itr = pointsDown.rbegin();
       itr != pointsDown.rend();
       ++itr){
    pair<Double_t,Double_t> currentPoint((*itr).first,(*itr).second);
    contour.push_back(currentPoint);
  }
  return contour;
}

void NuSystFitter::NoOscPrediction

(

)

[virtual]

Definition at line 86 of file NuSystFitter.cxx.

{

}

virtual void NuSystFitter::NumberOfContourPoints

(

const Int_t

numContourPoints

)

[inline, virtual]

Definition at line 46 of file NuSystFitter.h.

References fNumContourPoints.

    {fNumContourPoints = numContourPoints;}

const TGraph * NuSystFitter::OneSidedNeutrinoContourFinder

(

const NuMMParameters &

inputPars

)

[virtual]

Definition at line 230 of file NuSystFitter.cxx.

{
  std::vector<std::pair<Double_t,Double_t> > contPoints =
    this->NeutrinoContourFinder(inputPars);

  //Put the contour points in a TGraph
  TGraph* contour = new TGraph(contPoints.size());
  Int_t point=0;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contPoints.begin();
       it != contPoints.end();
       ++it){
    contour->SetPoint(point,(*it).first,(*it).second);
    cout << "Set contour point " << point
         << ", " << (*it).first
         << ", " << (*it).second
         << endl;
    ++point;
  }
  return contour;
}

double NuSystFitter::operator()

(

const std::vector< double > &

pars

)

const [virtual]

Definition at line 626 of file NuSystFitter.cxx.

References Likelihood().

{
  NuMMParameters mmPars(pars);
  //cout << "Called operator()" << endl;
  //mmPars.PrintStatus();
  return Likelihood(mmPars);
}

const TGraph * NuSystFitter::PlotContours

(

const NuMMParameters &

mmPars

)

[virtual]

Definition at line 167 of file NuSystFitter.cxx.

References NuMMParameters::ContourPars(), fNumContourPoints, NuMMParameters::kAlpha, NuMMParameters::kDm2, NuMMParameters::kDm2Bar, NuMMParameters::kMu2, NuMMParameters::kNcBack, NuMMParameters::kNorm, NuMMParameters::kNumParameters, NuMMParameters::kShwScale, NuMMParameters::kSn2, NuMMParameters::kSn2Bar, NuMMParameters::MinuitParameters(), and Munits::second.

{


  //Set up Minuit parameters
  ROOT::Minuit2::MnUserParameters mnPars = mmPars.MinuitParameters();

  //Run Minuit fit
  ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();

  if (!funcMin.IsValid()){
    cout << "Function minimum is not valid!" << endl;
  }
  else{
    cout << "Contour Minimization: " << endl;

    const ROOT::Minuit2::MnUserParameters upars = funcMin.UserParameters();

    cout << "Function minimum is: " << endl
         << "Sn2: " << upars.Parameter(NuMMParameters::kSn2).Value()
         << ", dm2: " << upars.Parameter(NuMMParameters::kDm2).Value()
         << "; norm: " << upars.Parameter(NuMMParameters::kNorm).Value()
         << "; shwEn: " << upars.Parameter(NuMMParameters::kShwScale).Value()
         << "; NCBack: " << upars.Parameter(NuMMParameters::kNcBack).Value()
         << "; Sn2bar: " << upars.Parameter(NuMMParameters::kSn2Bar).Value()
         << "; dm2bar: " << upars.Parameter(NuMMParameters::kDm2Bar).Value()
         << "; mu2: " << upars.Parameter(NuMMParameters::kMu2).Value()
         << "; alpha: " << upars.Parameter(NuMMParameters::kAlpha).Value()
         << "; CPT: " << upars.Parameter(NuMMParameters::kNumParameters).Value()
         << endl;
  }

  ROOT::Minuit2::MnContours contours(*this,funcMin);
  ROOT::Minuit2::ContoursError contErr = contours.Contour(mmPars.ContourPars().first,
                                                          mmPars.ContourPars().second,
                                                          fNumContourPoints);

  std::vector<std::pair<Double_t,Double_t> > vContPoints = contErr();
  //std::vector<std::pair<Double_t,Double_t> > vContPoints(10);// = contErr();

  Double_t* x = new Double_t[vContPoints.size()];
  for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
    x[pointCount] = vContPoints[pointCount].first;
  }

  Double_t* y = new Double_t[vContPoints.size()];
  for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
    y[pointCount] = vContPoints[pointCount].second;
  }
  TGraph* gCont = new TGraph(vContPoints.size(),y,x);
//  cout << "drawing contour" << endl;
//  gCont->SetName("gCont");
//  gDirectory->Print();
//  TFile* file=TFile::Open("contourplot.root","RECREATE");
//  gDirectory->Print();
//  gCont->Write("gCont");
//  file->Close();
  return gCont;
}

void NuSystFitter::push_back

(

NuMMRun *

run

)

[virtual]

Definition at line 61 of file NuSystFitter.cxx.

References fRuns, and run().

Referenced by NuFCFittingInterface::AddRun(), NuEZRunsFitter::BuildFitter(), NuFCFitter::NuFCFitter(), NuFCFitterNSI::NuFCFitterNSI(), NuFCFitterNSINu::NuFCFitterNSINu(), NuFCFitterNSINubar::NuFCFitterNSINubar(), NuEZFitterNSI::Rebuild(), and NuEZFitter::Rebuild().

{
  fRuns->push_back(run);
}

void NuSystFitter::QuietModeOff

(

)

[virtual]

Definition at line 673 of file NuSystFitter.cxx.

References fQuietMode, and fRuns.

{
  fQuietMode = 0;
  for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
       fRuns->end() != runIt;
       ++runIt){
    (*runIt)->QuietModeOff();
  }
}

void NuSystFitter::QuietModeOn

(

)

[virtual]

Definition at line 662 of file NuSystFitter.cxx.

References fQuietMode, and fRuns.

Referenced by BatchModeOn().

{
  fQuietMode = 1;
  for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
       fRuns->end() != runIt;
       ++runIt){
    (*runIt)->QuietModeOn();
  }
}

virtual void NuSystFitter::SetErrorDef

(

const Double_t

up

)

[inline, virtual]

Definition at line 24 of file NuSystFitter.h.

References fUp.

{fUp = up;}

virtual void NuSystFitter::Setup	(	const std::vector< std::string >	helperInputs,
		const std::vector< NuHistInterpolator * >	fdDataInput,
		const std::vector< NuHistInterpolator * >	ndDataInput,
		const std::string	xsecFile,
		const std::vector< std::string >	,
		const std::string	,
		Int_t
	)			[virtual]

const std::pair< Double_t, Double_t > NuSystFitter::SingleParDm2BarErrors	(	const NuMMParameters &	mmPars,
		const bool	quiet = false
	)			[virtual]

Definition at line 481 of file NuSystFitter.cxx.

{
  return SingleParErrors(mmPars,5,quiet);
}

const std::pair< Double_t, Double_t > NuSystFitter::SingleParDm2Errors	(	const NuMMParameters &	mmPars,
		const bool	quiet = false
	)			[virtual]

Definition at line 474 of file NuSystFitter.cxx.

{
  return SingleParErrors(mmPars,0,quiet);
}

const std::pair< Double_t, Double_t > NuSystFitter::SingleParErrors	(	const NuMMParameters &	mmPars,
		const int	idx,
		const bool	quiet
	)			[virtual]

Definition at line 495 of file NuSystFitter.cxx.

References NuMMParameters::AreAllParametersFixed(), MuELoss::e, NuMMParameters::FixParameterByIndex(), NuMMParameters::GetParameterByIndex(), NuMMParameters::MinuitPass(), NuMMParameters::SetParameterByIndex(), and NuMMParameters::VectorParameters().

{
  std::pair<Double_t,Double_t> errors;

  NuMMParameters bestFitPoint = this->Minimise(mmPars);
  Double_t bestFitChi2 = (*this)(bestFitPoint.VectorParameters());
  Double_t targetChi2 = bestFitChi2 + this->Up();
  Double_t bestFitPar = bestFitPoint.GetParameterByIndex(idx);

  cout << "Best chi2: " << bestFitChi2
       << " at par = " << bestFitPar
       << ". Target chi2: " << targetChi2
       << endl;

  for (Int_t direction = -1; direction < 2; direction += 2){

    if(!quiet) cout << "Direction = " << direction << endl;

    NuMMParameters scanMMPars;

//    Double_t incrementArray[5] = {0.1e-3,
//                                0.01e-3,
//                                0.001e-3,
//                                0.0001e-3,
//                                0.00001e-3};
    Double_t firstChi2 = bestFitChi2;
    Double_t firstPar = bestFitPar;
    Double_t currentChi2 = -1.0;
    Double_t lastChi2 = -1.0;
    Double_t lastPar = -1.0;
    Double_t increment = 1e-3*bestFitPar;

    while (currentChi2<targetChi2 && increment < 1e6){
      increment *= 2.0;
      scanMMPars = mmPars;
      scanMMPars.FixParameterByIndex(idx);
      scanMMPars.SetParameterByIndex(idx,bestFitPar + direction*increment);
      lastPar = scanMMPars.GetParameterByIndex(idx);
      if(!quiet) cout << "Trying par = " << scanMMPars.GetParameterByIndex(idx) << "... ";
      if (scanMMPars.AreAllParametersFixed()){
        currentChi2 = (*this)(scanMMPars.VectorParameters());
      }
      else{
        currentChi2 = (*this)(this->Minimise(scanMMPars).VectorParameters());
      }
      if(!quiet){
        cout << "Starting while with increment = " << increment
             << ", currentChi2 = " << currentChi2
             << ", lastPar = " << lastPar
             << endl;
      }
    }

    increment /= 2.0;

    while (TMath::Abs(currentChi2-targetChi2)>1e-5 && increment > 1e-5*bestFitPar){

      scanMMPars = mmPars;
      scanMMPars.FixParameterByIndex(idx);
      scanMMPars.SetParameterByIndex(idx,firstPar);

      if(!quiet){
        cout << "Starting while with firstChi2 = " << firstChi2
             << ", firstPar = " << firstPar
             << ", increment = " << increment
             << endl;
      }

      lastChi2 = currentChi2;
      lastPar = scanMMPars.GetParameterByIndex(idx);
      Double_t passminuit = true;
      scanMMPars.SetParameterByIndex(idx,scanMMPars.GetParameterByIndex(idx) + direction*increment);
      if(!quiet) cout << "Trying par = " << scanMMPars.GetParameterByIndex(idx) << "... ";
      if (scanMMPars.AreAllParametersFixed()){
        currentChi2 = (*this)(scanMMPars.VectorParameters());
      }
      else{
        NuMMParameters minPars = Minimise(scanMMPars);
        currentChi2 = (*this)(minPars.VectorParameters());
        passminuit = minPars.MinuitPass();
      }
      if(!quiet) cout << "chi2 = " << currentChi2 << endl;

      if(passminuit){
        if(currentChi2 > targetChi2) firstPar = lastPar;
        else firstPar = lastPar + direction*increment;
        increment /= 2.0;
      }
      else{
        firstPar = lastPar;
        increment *= 1.01;
      }
    }

    Double_t chi2Difference = lastChi2 - currentChi2;
    Double_t parDifference = lastPar - scanMMPars.GetParameterByIndex(idx);
    Double_t value = -1.0;
    if(!quiet){
      cout << "chi2Difference = " << -chi2Difference << endl;
      cout << "parDifference = " << parDifference << endl;
    }
    if ((fabs(chi2Difference)<1e-12) || (fabs(parDifference)<1e-12)){
      cout << "Chi2 or Par difference is zero. Supplying currentChi2."
           << endl;
      value = scanMMPars.GetParameterByIndex(idx);
    }
    else{
      Double_t fractionBetween = (targetChi2 - currentChi2)/chi2Difference;
      value = fractionBetween*parDifference + scanMMPars.GetParameterByIndex(idx);
    }
    if(!quiet){
      cout << "Error is between " << lastPar
           << " and " << scanMMPars.GetParameterByIndex(idx)
           << ", with chi2 between " << lastChi2
           << " and " << currentChi2
           << endl;
      cout << "Value is " << value << endl;
    }
    if (-1 == direction){errors.first = value;}
    if (1 == direction){errors.second = value;}
  }

  cout << "Final answer is par = " << bestFitPar
       << " + " << errors.second - bestFitPar
       << " - " << -errors.first + bestFitPar
       << endl;
  return errors;
}

const std::pair< Double_t, Double_t > NuSystFitter::SingleParSn2BarErrors	(	const NuMMParameters &	mmPars,
		const bool	quiet = false
	)			[virtual]

Definition at line 488 of file NuSystFitter.cxx.

{
  return SingleParErrors(mmPars,6,quiet);
}

NuMMParameters NuSystFitter::Sn2GridSearch	(	Double_t	Dm2,
		Double_t	Dm2Bar,
		Int_t	resolution = 0
	)			[protected]

Does a grid search over sn2/sn2bar, with resolution bins per parameter.

Definition at line 2243 of file NuSystFitter.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), NuMMParameters::Sn2(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by FillCPT().

{
  if (resolution <= 0) resolution = 11;
  // Use a 2D histogram to define the search grid
  TH2D sinsurf("sinsurf","sinsurf",resolution,0,1,resolution,0,1);
  // Create a template parameters object
  NuMMParameters pars;
  pars.Dm2(Dm2);
  pars.Dm2Bar(Dm2Bar);

  Double_t minLike, minSn2, minSn2Bar;
  minLike = minSn2 = minSn2Bar = -1;

  // Loop over every bin
  for (Int_t xBin = 1; xBin <= sinsurf.GetNbinsX(); ++xBin) {
    for (Int_t yBin = 1; yBin <= sinsurf.GetNbinsY(); ++yBin) {
      Double_t sn2 = sinsurf.GetYaxis()->GetBinCenter(yBin);
      Double_t sn2bar = sinsurf.GetXaxis()->GetBinCenter(xBin);
      pars.Sn2(sn2);
      pars.Sn2Bar(sn2bar);

      // Evaluate the likelihood
      const Double_t likelihood = (*this)(pars.VectorParameters());

      // Check for minimum (or first entry)
      if ((xBin == 1 && yBin == 1) || likelihood < minLike) {
        minLike = likelihood;
        minSn2 = sn2;
        minSn2Bar = sn2bar;
      }
    }
  }

  // Check Zero and One
  pars.Sn2(0.0);
  pars.Sn2Bar(0.0);
  Double_t likelihood = (*this)(pars.VectorParameters());
  if (likelihood < minLike) {
    minLike = likelihood;
    minSn2 = 0.0;
    minSn2Bar = 0.0;
  }

  // Check Zero and One
  pars.Sn2(1.0);
  pars.Sn2Bar(1.0);
  likelihood = (*this)(pars.VectorParameters());
  if (likelihood < minLike) {
    minLike = likelihood;
    minSn2 = 1.0;
    minSn2Bar = 1.0;
  }

  // We should now have a minimum. Sanity check this
  assert(minLike > 0);
  assert(minSn2 > 0);
  assert(minSn2Bar > 0);

  // Now build our return
  pars.Sn2(minSn2);
  pars.Sn2Bar(minSn2Bar);

  return pars;
}

virtual void NuSystFitter::Sn2PointsForNeutrinoContour

(

const std::vector< Double_t > &

sn2PointsForNeutrinoContour

)

[inline, virtual]

Definition at line 62 of file NuSystFitter.h.

References fsn2PointsForNeutrinoContour.

    {fsn2PointsForNeutrinoContour = sn2PointsForNeutrinoContour;}

const TGraph * NuSystFitter::TwoSidedNeutrinoContourFinder

(

const NuMMParameters &

inputPars

)

[virtual]

Definition at line 254 of file NuSystFitter.cxx.

{
  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourPhysical;
  std::vector<std::pair<Double_t,Double_t> > contPointsPhys =
    this->NeutrinoContourFinder(inputPars);

  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourUnphysical;
  std::vector<std::pair<Double_t,Double_t> > contPointsUnphys =
    this->NeutrinoContourFinder(inputPars);

  //Make a single vector of points
  std::vector<pair<Double_t,Double_t> > contourPoints;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contPointsPhys.begin();
       it != contPointsPhys.end();
       ++it){
    pair<Double_t,Double_t> currentPoint((*it).first,(*it).second);
    contourPoints.push_back(currentPoint);
  }
  for (vector<pair<Double_t, Double_t> >::reverse_iterator
         itr = contPointsUnphys.rbegin();
       itr != contPointsUnphys.rend();
       ++itr){
    pair<Double_t,Double_t> currentPoint((*itr).first,(*itr).second);
    contourPoints.push_back(currentPoint);
  }


  //Put the contour points in a TGraph
  TGraph* contour = new TGraph(contourPoints.size());
  Int_t point=0;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contourPoints.begin();
       it != contourPoints.end();
       ++it){
    contour->SetPoint(point,(*it).first,(*it).second);
    cout << "Set contour point " << point
         << ", " << (*it).first
         << ", " << (*it).second
         << endl;
    ++point;
  }
  return contour;
}

virtual double NuSystFitter::Up

(

)

const [inline, virtual]

Definition at line 23 of file NuSystFitter.h.

References fUp.

{return fUp;}//Chi2 errors; 0.5 for likelihood

---

Member Data Documentation

Int_t NuSystFitter::fNumContourPoints [protected]

Definition at line 177 of file NuSystFitter.h.

Referenced by NumberOfContourPoints(), and PlotContours().

Int_t NuSystFitter::fQuietMode [protected]

Definition at line 175 of file NuSystFitter.h.

Referenced by BatchModeOn(), Chi2SliceDm2(), Chi2SliceDm2bar(), Chi2SliceSin2bar(), Chi2ValleyBar(), Chi2ValleyBarDm2Bar(), Chi2ValleyDm2(), Chi2ValleySn2(), FillLikelihoodSurfaceCPT(), FillLikelihoodSurfaceGeneric(), FillLikelihoodSurfaceNSIDm2Eps(), FillLikelihoodSurfaceNSIDm2Sn2(), FillLikelihoodSurfaceNSISn2Eps(), Likelihood(), LikelihoodSurfaceBar(), Minimise(), QuietModeOff(), and QuietModeOn().

std::vector<NuMMRun*>* NuSystFitter::fRuns [protected]

Definition at line 181 of file NuSystFitter.h.

Referenced by NuABFitter::Archive(), Likelihood(), push_back(), QuietModeOff(), and QuietModeOn().

std::vector<Double_t> NuSystFitter::fsn2PointsForNeutrinoContour [protected]

Definition at line 178 of file NuSystFitter.h.

Referenced by Sn2PointsForNeutrinoContour().

std::vector<Double_t> NuSystFitter::fsn2PointsForNeutrinoContourPhysical [protected]

Definition at line 179 of file NuSystFitter.h.

std::vector<Double_t> NuSystFitter::fsn2PointsForNeutrinoContourUnphysical [protected]

Definition at line 180 of file NuSystFitter.h.

Double_t NuSystFitter::fUp [protected]

Definition at line 176 of file NuSystFitter.h.

Referenced by SetErrorDef(), and Up().

---

The documentation for this class was generated from the following files:

• NuSystFitter.h
• NuSystFitter.cxx

---