#ifndef __CINT__
#include "Riostream.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMatrixDLazy.h"
#include "TVectorD.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"
#endif

void energy2()
{
#ifdef __CINT__
  gSystem->Load("libMatrix");
#endif
     /*
      * Provided Courtesy of Jon Habif!
      */
    const Int_t NPoints = 100;

    /*
     * Unit of quantum flux.
     */
    Double_t phi0 = TMath::H()/(2.0*TMath::Qe());
    Double_t df   = 5.0E-3;  // Total flux sweep
    Double_t diff = 0.0;
    Double_t t1   = 1.9;     // GHz of energy splitting. 
    Double_t t2   = 0.45;   
    Double_t ec   = 0.5;     // 500 MHz  coupling energy.
    Double_t StartFlux = 0.5 - df;
    Double_t EndFlux   = 0.5 + df;
    //    Double_t Step = 1.001E-6;
    Double_t Step = (2.0*df)/((Double_t) NPoints);

    Double_t f1, f2;
    Double_t *flux1, *flux2, *energy1, *energy2; 
    Double_t eps1, eps2;
    Int_t    i;
    TMatrixD H(4,4);
    TVectorD ev(4);

    // root invasion ---------------------------------------------------
    TMultiGraph *tmg = new TMultiGraph();
    TGraph *E1       = new TGraph(NPoints);
    flux1            = E1->GetX();
    energy1          = E1->GetY();
    tmg->Add(E1);

    TGraph *E2       = new TGraph(NPoints);
    flux2            = E2->GetX();
    energy2          = E2->GetY();
    tmg->Add(E2);

    TCanvas *Hobbes = new TCanvas("Dist","A test",5,5,1200,600);
    Hobbes->cd();
    TPad    *Calvin = new TPad("Calvin","Silly",0.01,0.01,0.99,0.99,10);
    Calvin->Draw();
    Calvin->SetGrid();
    Calvin->cd();

    // ------------------------------------------------------------------
    H(0,1) = t1;
    H(1,0) = t1;
    H(0,2) = t2;
    H(2,0) = t2;
    H(1,3) = t2;
    H(2,3) = t1;
    H(3,1) = t2;
    H(3,2) = t1;

    for(i=0;i<NPoints;i++)
    {
	f1 = StartFlux + (Double_t)i * Step;
	f2 = f1 - diff;
	flux1[i] = f1+f2;
	flux2[i] = f1+f2;
	eps1  = 1.0E-9 * (392.0E-9 * phi0 * (f1-0.5))/TMath::H();
	eps2  = 1.0E-9 * (512.0E-9 * phi0 * (f2-0.5))/TMath::H();
	H( 0,0)  =  eps1 + eps2 + ec;
	H( 1,1)  = -eps1 + eps2 - ec;
	H( 2,2)  =  eps1 - eps2 - ec;
	H( 3,3)  = -eps1 - eps2 + ec;

	TMatrixDEigen EG(H);
	ev = EG.GetEigenValuesRe();
	energy1[i] = ev[0];
	energy2[i] = ev[1];
    }

    tmg->Draw("ALP");
    TH1F *h = tmg->GetHistogram();
    if (h)
    {
	h->SetXTitle("Flux");
	h->SetYTitle("Energy");
	h->SetLabelSize(0.03,"X");
	h->SetLabelSize(0.03,"Y");
    }
}