void vis()
{
    
    // Initialize variables
    TString graph_title;
    
    Int_t N_Points;
    Int_t Exp_Points = 0;
    Int_t counter_offset = 0;
    
    Float_t x;
    Float_t Bx_theory;
    Float_t Bx_exp;
    Float_t sigma_Bx;
    Float_t x_offset = -0.25;
    
    Float_t exp_p0;
    Float_t exp_p1;
    Float_t exp_p2;
    Float_t th_p0;
    Float_t th_p1;
    Float_t th_p2;
    
    // Fiducal volume
    Float_t mid = 2.25;
    Float_t half_w = 1.05;
    
    // Input data
    TFile *input = new TFile("x.root");
    input_data = (TTree*)input->Get("x");
    
    // Retrieve ntuple data
    input_data->SetBranchAddress("x", &x);
    input_data->SetBranchAddress("Bx_theory", &Bx_theory);
    input_data->SetBranchAddress("Bx_exp", &Bx_exp);
    input_data->SetBranchAddress("sigma_Bx", &sigma_Bx);
    N_Points = (Int_t)input_data->GetEntries();
    
    for(Int_t i = 0; i < N_Points; i++)
    {
        input_data->GetEntry(i);
        if(Bx_exp != 0)
        {
            Exp_Points++;
        }
    }
    
    // Prepare for graphing
    TCanvas *c = new TCanvas("c","Xtal Behavior", 100, 200, 500, 400);
    c->SetFrameFillStyle(4000);
    TGraph *theory = new TGraph(N_Points);  
    TGraphErrors *experiment = new TGraphErrors(Exp_Points);
    
    // Loop through
    for(Int_t i = 0; i < N_Points; i++)
    {     
                
        input_data->GetEntry(i);

        theory->SetPoint(i, x, Bx_theory);
        if(Bx_exp != 0)
        {
            experiment->SetPoint(i - counter_offset, x + x_offset, Bx_exp);
            experiment->SetPointError(i - counter_offset, 0, sigma_Bx);
        }
        else
        {
            counter_offset++;
        }
        
    }

    // Fit data    
    TF1 *th_fit = new TF1("th_fit", "[0] * (x - [1]) * (x - [1]) + [2]");
    th_fit->SetLineColor(2);
    th_fit->SetLineWidth(0.5);
    theory->Fit("th_fit", "Q", "", -3.5, 3.5);
    th_p0 = th_fit->GetParameter(0);
    th_p1 = th_fit->GetParameter(1);
    th_p2 = th_fit->GetParameter(2);

    TF1 *exp_fit = new TF1("exp_fit", "[0] * (x - [1]) * (x - [1]) + [2]");
    exp_fit->SetLineWidth(0.5);
    experiment->Fit("exp_fit", "Q", "", -3.5, 3.5);
    exp_p0 = exp_fit->GetParameter(0);
    exp_p1 = exp_fit->GetParameter(1);
    exp_p2 = exp_fit->GetParameter(2);

    cout << th_p0 << "\t" << exp_p0 << endl;

    //Draw a frame and display the boxes
    c->DrawFrame(-3.14, -0.025,3.14,0.14);
    f_volume = new TBox(-(mid - half_w), -0.03, mid - half_w, 0.14);
    f_volume->SetFillStyle(3013);
    f_volume->SetFillColor(17);
    f_volume->Draw();
    
    f_left = new TBox(-3.5, -.03, -(mid + half_w), 0.14);
    f_left->SetFillStyle(3013);
    f_left->SetFillColor(17);
    f_left->Draw();
    
    f_right = new TBox(mid + half_w, -.03, 3.5, 0.14);
    f_right->SetFillStyle(3013);
    f_right->SetFillColor(17);
    f_right->Draw(); 
    
    // Display data
    graph_title = "Cos #theta Coil Uniformity";
    
    theory->SetMarkerStyle(1);
    theory->SetMarkerColor(0);
    
    experiment->SetMarkerStyle(7);
    experiment->SetMarkerColor(1);
    
    TMultiGraph *visual = new TMultiGraph();
    visual->SetTitle(graph_title);
    visual->Add(theory);
    visual->Add(experiment);
    visual->Draw("P");  //was "AP"
    c->Modified();
    c->Update();
    visual->GetXaxis()->SetTitle("x (cm)");
    visual->GetXaxis()->SetRangeUser(-3.5, 3.5);
    visual->GetYaxis()->SetTitle("Deviation (%)");
    visual->GetYaxis()->SetTitleOffset(1.2);
    visual->GetYaxis()->SetRangeUser(-0.03, 0.14);
    
    // Create legend
    TLegend *legend = new TLegend(.77, 0.86, 0.99, 0.97);
    legend->SetFillColor(0);
    legend->AddEntry(th_fit, "Simulation", "l");
    legend->AddEntry(experiment, "Experiment", "p");
    legend->Draw();
    


}