energyhistogram.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <string>
#include <TCanvas.h>
#include <TH2F.h>
#include <TMath.h>
#include <TFile.h>
#include <TGraph.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TColor.h>
#include <TTimeStamp.h>
#define AUGER_RADIO_ENABLED
#include "RecEventFile.h"
#include "RecEvent.h"
#include "RdEvent.h"
#include "DetectorGeometry.h"
#include "SDEvent.h"
#include "DetectorGeometry.h"
#include "SdRecStation.h"
#include "ShowerRRecDataQuantities.h"
#include "StationRRecDataQuantities.h"
#include "PolarHist.h"
#include <TEllipse.h>
#include <TGraphPolar.h>
#include <TPaletteAxis.h>
#include <TText.h>
//Added
#include <TVector3.h>
#include <TProfile.h>
#include <RdTrace.h>
#include <vector>
#include <TArrow.h>

using namespace std;

void reverse(double arr[], int count)
{
   double temp;
   for (int i = 0; i < count/2; ++i)
   {
      temp = arr[i];
      arr[i] = arr[count-i-1];
      arr[count-i-1] = temp;
   }
}

//from https://root-forum.cern.ch/t/set-real-aspect-ratio-of-th2-or-tgraph/20734
bool SetRealAspectRatio(TCanvas* const c, const Int_t axis = 1)
{
    if(!c)
    {
        cout << "Error in SetRealAspectRatio: canvas is NULL";
        return false;
    }
    
    {
        //Get the current min-max values if SetUserRange has been called
        c->Update();
        const Double_t xmin = c->GetUxmin();
        const Double_t xmax = c->GetUxmax();
        const Double_t ymin = c->GetUymin();
        const Double_t ymax = c->GetUymax();

        //Get the length of zoomed x and y axes
        const Double_t xlength = xmax - xmin;
        const Double_t ylength = ymax - ymin;
        const Double_t ratio = xlength/ylength;

        //Get how many pixels are occupied by the canvas
        const Int_t npx = c->GetWw();
        const Int_t npy = c->GetWh();
        
        //Get x-y coordinates at the edges of the canvas
        //(extrapolating outside the axes, NOT at the edges of the histogram)
        const Double_t x1 = c->GetX1();
        const Double_t y1 = c->GetY1();
        const Double_t x2 = c->GetX2();
        const Double_t y2 = c->GetY2();
        
        //Get the length of extrapolated x and y axes
        const Double_t xlength2 = x2 - x1;
        const Double_t ylength2 = y2 - y1;
        const Double_t ratio2 = xlength2/ylength2;
        
        //Get now number of pixels including window borders
        const Int_t bnpx = c->GetWindowWidth();
        const Int_t bnpy = c->GetWindowHeight();
        
        /*
        cout << "WindX\tWindY\tCanvX\tCanvY\tx1\ty1\tx2\ty2\tratiox/y\tCanvX/CanvY" << endl;
        cout << bnpx << "\t" << bnpy << "\t" << npx << "\t" << npy 
             << "\t" << x1 << "\t" << y1 << "\t" << x2 << "\t" << y2 << "\t" 
             << ratio2 << "\t" << (double)npx/npy << "\tOriginal Canvas" << endl;
        */
        
        if(axis==1)
        {
            c->SetCanvasSize(TMath::Nint(npy*ratio2), npy);
            c->SetWindowSize((bnpx-npx)+TMath::Nint(npy*ratio2), bnpy);
        }
        else if(axis==2)
        {
            c->SetCanvasSize(npx, TMath::Nint(npx/ratio2));
            c->SetWindowSize(bnpx, (bnpy-npy)+TMath::Nint(npx/ratio2));
        }
        else
        {
            cout << "Error in SetRealAspectRatio: axis value " << axis 
                 << " is neither 1 (resize along x-axis) nor 2 (resize along y-axis).";
            return false;
        }
    }
    
    //Check now that resizing has worked
    {
        //Get the current min-max values if SetUserRange has been called
        c->Update();
        const Double_t xmin = c->GetUxmin();
        const Double_t xmax = c->GetUxmax();
        const Double_t ymin = c->GetUymin();
        const Double_t ymax = c->GetUymax();

        //Get the length of zoomed x and y axes
        const Double_t xlength = xmax - xmin;
        const Double_t ylength = ymax - ymin;
        const Double_t ratio = xlength/ylength;

        //Get how many pixels are occupied by the canvas
        const Int_t npx = c->GetWw();
        const Int_t npy = c->GetWh();
        
        //Get x-y coordinates at the edges of the canvas 
        //(extrapolating outside the axes, NOT at the edges of the histogram)
        const Double_t x1 = c->GetX1();
        const Double_t y1 = c->GetY1();
        const Double_t x2 = c->GetX2();
        const Double_t y2 = c->GetY2();
        
        //Get the length of extrapolated x and y axes
        const Double_t xlength2 = x2 - x1;
        const Double_t ylength2 = y2 - y1;
        const Double_t ratio2 = xlength2/ylength2;
        
        //Get now number of pixels including window borders
        const Int_t bnpx = c->GetWindowWidth();
        const Int_t bnpy = c->GetWindowHeight();

        /*
        cout << bnpx << "\t" << bnpy << "\t" << npx << "\t" << npy << "\t" 
             << x1 << "\t" << y1 << "\t" << x2 << "\t" << y2 << "\t" << ratio2 
             << "\t" << (double)npx/npy << "\tModified Canvas" << endl;*/
        
        //Check accuracy +/-1 pixel due to rounding
        if(abs(TMath::Nint(npy*ratio2) - npx)<2)
        {
            cout << "Resizing finished successfully." << endl;
            return true;
        }
        else
        {
            cout << "Resizing failed." << endl;
            return false;
        }
    }
    // References:
    // https://root.cern.ch/root/roottalk/roottalk01/3676.html
    // https://root-forum.cern.ch/t/making-the-both-axes-square-on-the-pad/4325/1
}

int main() {
    //TH1F* Hist_Energy_sd  = new TH1F("Hist_Energy_sd", "Hist_Energy_sd", 40, 17, 20);
    //TH1F* Hist_Energy_sd  = new TH1F("Hist_Energy_sd", "Hist_Energy_sd", 40, 17.5, 20);
    TH1F* Hist_Energy_sd  = new TH1F("Hist_Energy_sd", "Hist_Energy_sd", 15, 18.5, 20);
    TH1F* Hist_Zenith_sd  = new TH1F("Hist_Zenith_sd", "Hist_Zenithy_sd", 25, 60, 90);
    TH1F* Hist_Azimuth_sd = new TH1F("Hist_Azimuth_sd", "Hist_Azimuth_sd", 10, 0, 360);
    //0.9891 = sin^2(84 degree)
    TH1F* Hist_sin2theta  = new TH1F("foo", "all E", 11, 0.75, 0.9891);
    TH1F* Hist_sin2theta2 = new TH1F("foo2", "lg E > 18.2", 12, 0.75, 0.9891);
    TH1F* Hist_sin2theta3 = new TH1F("foo3", "lg E > 18.5", 12, 0.75, 0.9891);

    vector<string> v_filename;
    ifstream fileList ("./filelist.txt");
    char filename[10000] ;
    while (fileList >> filename) {
        v_filename.push_back(filename);
    }

    const TVector3 vmag(-884.62,-19699.87,+14154.43);
    RecEventFile *file = new RecEventFile(v_filename);
    RecEvent* theevent= new RecEvent();
    int nevents=file->GetNEvents();
    file->SetBuffers(&theevent);
    DetectorGeometry *geo =new DetectorGeometry();
    file->ReadDetectorGeometry(*geo);

    std::vector<double> pphi;
    std::vector<double> ttheta;

    //Loop over events
    while(file->ReadNextEvent()==RecEventFile::eSuccess) {
        const RdEvent& theradio= theevent->GetRdEvent();
        const RdRecShower& rshow=theradio.GetRdRecShower();

        if (!theradio.HasRdStations()) continue ;
        const std::vector<RdRecStation>& vecRstat = theradio.GetRdStationVector();
        const SDEvent& sdev = theevent->GetSDEvent();
        const SdRecShower& sshow = sdev.GetSdRecShower();

        const TVector3& score = sshow.GetCoreSiteCS();
        const TVector3& rcore = rshow.GetCoreSiteCS();

        const TVector3& saxis = sshow.GetAxisSiteCS();

        double rtheta = TMath::RadToDeg()*rshow.GetZenith();
        double rphi = TMath::RadToDeg()*rshow.GetAzimuth();
        double stheta = TMath::RadToDeg()*sshow.GetZenith();
        double sphi = TMath::RadToDeg()*sshow.GetAzimuth();

        if (stheta<60 || stheta>84)
            continue;

        //SD quality check
        const int minRecLevel=3; // has reconstructed energy
        const int minCandidateStations=4;
        const int T4Trigger=2; // has T4 and SD triggered independently of FD
        const int T5Trigger=3; //has T5 Posterior, which is the T5 HAS in this case

        const float minZenithSD=60;
        const float maxZenithSD=84;

//        if (sdev.GetRecLevel()<minRecLevel) continue;// alle events
//        if (sdev.GetNumberOfCandidates()<minCandidateStations) continue; //alle events
//        if (sdev.GetT4Trigger()<T4Trigger) continue;
//        if (sdev.GetT5Trigger()<T5Trigger) continue;
        if (stheta>maxZenithSD) continue;
        if (stheta<minZenithSD) continue;
//        if (log10(sshow.GetEnergy()) < 18.5) continue;

        //Hist SD energy
        double energy_sd = sshow.GetEnergy();
        double energy_err_sd = sshow.GetEnergyError();
        const double power_sd = energy_sd>0.?(int) log10(energy_sd):0.;
        const double norm_sd =pow(10, power_sd);
        energy_err_sd /= norm_sd;
        Hist_Energy_sd->Fill(log10(energy_sd));

        //Hist Zenith angle
        double zenith_sd = sshow.GetZenith()*TMath::RadToDeg();
        double zenith_err_sd = sshow.GetZenithError()*TMath::RadToDeg();
        Hist_Zenith_sd->Fill(zenith_sd);

        // Hist Azimuth angle
        double azimuth_sd = sshow.GetAzimuth()*TMath::RadToDeg();
        double azimuth_err_sd = sshow.GetAzimuthError()*TMath::RadToDeg();
        Hist_Azimuth_sd->Fill(azimuth_sd);

        // Hist sin2theta
        double sintheta = sin(sshow.GetZenith());
        Hist_sin2theta->Fill(sintheta*sintheta);
        if(log10(sshow.GetEnergy()) > 18.2)
           Hist_sin2theta2->Fill(sintheta*sintheta);
        if(log10(sshow.GetEnergy()) > 18.5)
           Hist_sin2theta3->Fill(sintheta*sintheta);
        
        pphi.push_back(sphi);
        ttheta.push_back(stheta);
    }
    //Saving histogramms
    TFile *f = new TFile("energyhistogram.root","RECREATE");
    Hist_Energy_sd->Write();
    f->Close();

    //Making plots
    TCanvas Histo("Histo","Histo",1800,1000);
    Histo.SetLeftMargin(0.1);
    Histo.SetRightMargin(0.05);
    Histo.SetTopMargin(0.06);
    Histo.SetBottomMargin(0.14);

    Hist_Energy_sd->SetStats(kFALSE);
    Hist_Energy_sd->GetXaxis()->SetTitle("log_{10}(E_{SD}/eV)");
    Hist_Energy_sd->GetXaxis()->SetLabelFont(42);
    Hist_Energy_sd->GetXaxis()->SetLabelOffset(0.01);
    Hist_Energy_sd->GetXaxis()->SetLabelSize(0.05);
    Hist_Energy_sd->GetXaxis()->SetTitleSize(0.055);
    Hist_Energy_sd->GetXaxis()->SetTitleOffset(1.15);
    Hist_Energy_sd->GetXaxis()->SetNdivisions(804);
    Hist_Energy_sd->GetXaxis()->SetRangeUser(18.4, 20.0);
    Hist_Energy_sd->GetXaxis()->SetTitleFont(42);
    Hist_Energy_sd->GetYaxis()->SetTitle("Entries");
    Hist_Energy_sd->GetYaxis()->SetLabelFont(42);
    Hist_Energy_sd->GetYaxis()->SetLabelOffset(0.01);
    Hist_Energy_sd->GetYaxis()->SetLabelSize(0.05);
    Hist_Energy_sd->GetYaxis()->SetTitleSize(0.055);
    Hist_Energy_sd->GetYaxis()->SetTitleOffset(0.8);
    Hist_Energy_sd->GetYaxis()->SetTitleFont(42);
    //Hist_Energy_sd->SetMarkerStyle(20);
    //Hist_Energy_sd->SetMarkerSize(0.4);
    //Hist_Energy_sd->SetMarkerColor(kBlack);
    //Hist_Energy_sd->GetXaxis()->SetTitleOffset(1.2);
    //Hist_Energy_sd->GetYaxis()->SetTitleOffset(1.2);
    Hist_Energy_sd->SetTitle("");
    //Hist_Energy_sd->SetFillColor(4);
    Hist_Energy_sd->SetLineWidth(3);
    Hist_Energy_sd->SetLineColor(kBlack);
    Hist_Energy_sd->Draw("HIST E");

    TH1F* test = (TH1F*)Hist_Energy_sd->Clone();
//    test->SetLineWidth(3);
//    test->Draw("E SAME");

    Histo.SetLogy();
    Histo.SaveAs("Hist_Energy.pdf");

    Hist_Zenith_sd->SetStats(kFALSE);
    Hist_Zenith_sd->GetXaxis()->SetTitle("Zenith angle / #circ");
    Hist_Zenith_sd->GetYaxis()->SetTitle("Entries");
    Hist_Zenith_sd->GetXaxis()->SetRangeUser(0.7, 1.);
    Hist_Zenith_sd->GetXaxis()->SetLabelFont(42);
    Hist_Zenith_sd->GetXaxis()->SetLabelOffset(0.01);
    Hist_Zenith_sd->GetXaxis()->SetLabelSize(0.05);
    Hist_Zenith_sd->GetXaxis()->SetTitleSize(0.055);
    Hist_Zenith_sd->GetXaxis()->SetTitleOffset(1.15);
    Hist_Zenith_sd->GetXaxis()->SetTitleFont(42);
    Hist_Zenith_sd->GetYaxis()->SetLabelFont(42);
    Hist_Zenith_sd->GetYaxis()->SetLabelOffset(0.01);
    Hist_Zenith_sd->GetYaxis()->SetLabelSize(0.05);
    Hist_Zenith_sd->GetYaxis()->SetTitleSize(0.055);
    Hist_Zenith_sd->GetYaxis()->SetTitleOffset(0.8);
    Hist_Zenith_sd->GetYaxis()->SetTitleFont(42);
    Hist_Zenith_sd->SetTitle("");
    Hist_Zenith_sd->SetFillColor(4);
    Hist_Zenith_sd->Draw("HIST");
    Histo.SaveAs("Hist_Zenith.pdf");

    /*
    gStyle->SetStatY(0.9);
    gStyle->SetStatX(0.5);                
    gStyle->SetStatW(0.3);                
    gStyle->SetStatH(0.2);                
    */

    TCanvas Histo2("Histo2","Histo2",1800,1000);
    Histo2.SetLeftMargin(0.1);
    Histo2.SetRightMargin(0.05);
    Histo2.SetTopMargin(0.06);
    Histo2.SetBottomMargin(0.14);
//    Histo2.Divide(3);
//    Histo2.cd(1);

    Hist_sin2theta->GetXaxis()->SetTitle("sin^{2}(#theta_{SD})");
    Hist_sin2theta->GetYaxis()->SetTitle("Entries");
    Hist_sin2theta->GetXaxis()->SetRangeUser(0.7, 1.);
    Hist_sin2theta->GetXaxis()->SetLabelFont(42);
    Hist_sin2theta->GetXaxis()->SetLabelOffset(0.01);
    Hist_sin2theta->GetXaxis()->SetLabelSize(0.05);
    Hist_sin2theta->GetXaxis()->SetTitleSize(0.055);
    Hist_sin2theta->GetXaxis()->SetTitleOffset(1.15);
    Hist_sin2theta->GetXaxis()->SetTitleFont(42);
    Hist_sin2theta->GetYaxis()->SetLabelFont(42);
    Hist_sin2theta->GetYaxis()->SetLabelOffset(0.01);
    Hist_sin2theta->GetYaxis()->SetLabelSize(0.05);
    Hist_sin2theta->GetYaxis()->SetTitleSize(0.055);
    Hist_sin2theta->GetYaxis()->SetTitleOffset(0.8);
    Hist_sin2theta->GetYaxis()->SetTitleFont(42);
    Hist_sin2theta->SetStats(kFALSE);
    Hist_sin2theta->SetTitle("");
    Hist_sin2theta->SetLineColor(kBlack);
    //Hist_sin2theta->SetFillColor(4);
    Hist_sin2theta->SetLineWidth(3);
    Hist_sin2theta->Draw("HIST E");

    //test = (TH1F*)Hist_sin2theta->Clone();
    //test->SetLineWidth(3);
    //test->Draw("E SAME");

/*
    //Histo2.cd(2);
    Hist_sin2theta2->SetLineColor(kBlack);
    Hist_sin2theta2->SetFillColor(5);
    Hist_sin2theta2->Draw("HIST SAME");

    //Histo2.cd(3);
    Hist_sin2theta3->SetLineColor(kBlack);
    Hist_sin2theta3->SetFillColor(6);
    Hist_sin2theta3->Draw("HIST SAME");
*/
    Histo2.SetLogy(0);
    Histo2.SaveAs("Hist_sin2theta.pdf");

    TCanvas Histo3("Histo2","Histo2",1800,1000);
    Histo3.SetLeftMargin(0.1);
    Histo3.SetRightMargin(0.05);
    Histo3.SetTopMargin(0.06);
    Histo3.SetBottomMargin(0.14);
    Hist_Azimuth_sd->SetStats(kFALSE);
    Hist_Azimuth_sd->GetXaxis()->SetTitle("#phi_{SD} / #circ");
    Hist_Azimuth_sd->GetYaxis()->SetTitle("Entries");
    Hist_Azimuth_sd->GetXaxis()->SetLabelFont(42);
    Hist_Azimuth_sd->GetXaxis()->SetLabelOffset(0.01);
    Hist_Azimuth_sd->GetXaxis()->SetLabelSize(0.05);
    Hist_Azimuth_sd->GetXaxis()->SetTitleSize(0.055);
    Hist_Azimuth_sd->GetXaxis()->SetTitleOffset(1.15);
    Hist_Azimuth_sd->GetXaxis()->SetTitleFont(42);
    Hist_Azimuth_sd->GetYaxis()->SetLabelFont(42);
    Hist_Azimuth_sd->GetYaxis()->SetLabelOffset(0.01);
    Hist_Azimuth_sd->GetYaxis()->SetLabelSize(0.05);
    Hist_Azimuth_sd->GetYaxis()->SetTitleSize(0.055);
    Hist_Azimuth_sd->GetYaxis()->SetTitleOffset(0.8);
    Hist_Azimuth_sd->GetYaxis()->SetTitleFont(42);
    Hist_Azimuth_sd->SetTitle("");
    //Hist_Azimuth_sd->GetYaxis()->SetNdivisions(810);
    Hist_Azimuth_sd->GetYaxis()->SetRangeUser(0, 100);
    Hist_Azimuth_sd->SetLineColor(kBlack);
    //Hist_Azimuth_sd->SetFillColor(4);
    Hist_Azimuth_sd->SetLineWidth(3);
    Hist_Azimuth_sd->Draw("HIST E");

    //test = (TH1F*)Hist_Azimuth_sd->Clone();
    //test->SetLineWidth(3);
    //test->Draw("E SAME");

    Histo3.SaveAs("Hist_Azimuth.pdf");

    TCanvas canvas("Histo","Histo",1800,1000);
    canvas.SetLeftMargin(0.07);
    canvas.SetRightMargin(0.28);
    canvas.SetTopMargin(0.06);
    canvas.SetBottomMargin(0.06);
    //TH2F skymap("skymap", "", 72, 0, 360, 30, 0, 90); // 5x3 Grad bins
    TH2F skymap("skymap", "", 36, 0, 360, 18, 0, 90); // 10x5 Grad bins
   
    skymap.SetStats(kFALSE);
    for (int k=0; k<pphi.size(); k++){
        skymap.Fill(pphi[k], ttheta[k]);
    }
    cout << pphi.size() << " entries in map" << endl;

    double maxBinContent = skymap.GetMaximum();
    for (int i=0; i<skymap.GetSize(); i++){
        skymap.SetBinContent(i, skymap.GetBinContent(i)/maxBinContent);
    }


    Double_t stops[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 
                          0.6250, 0.7500, 0.8750, 1.0000};
    double alpha = 1.;

    //Viridis    
/*    Double_t red[9] = { 26./255., 51./255.,  43./255.,  33./255.,  28./255.,  
                        35./255.,  74./255., 144./255., 246./255.};
    Double_t green[9] = {  9./255., 24./255.,  55./255.,  87./255., 118./255.,
                        150./255., 180./255., 200./255., 222./255.};
    Double_t blue[9] = { 30./255., 96./255., 112./255., 114./255., 112./255.,
                        101./255.,  72./255.,  35./255.,   0./255.};
*/

    //inverted kCherry
/*    Double_t red[9] = { 251./255., 235./255., 223./255., 214./255., 196./255., 
                        188./255., 157./255., 102./255., 37./255.};
    Double_t green[9] = { 251./255., 185./255., 132./255., 91./255., 67./255.,
                        37./255., 25./255., 29./255., 37./255.};
    Double_t blue[9] = { 251./255., 137./255., 98./255., 66./255., 45./255.,
                        33./255., 32./255., 37./255.};
*/
    //ocean
/*    Double_t red[9] = { 14./255.,  7./255.,  2./255.,  0./255.,  5./255.,
                            11./255.,  55./255., 131./255., 229./255.};
    Double_t green[9] = {105./255., 56./255., 26./255.,  1./255., 42./255.,
                        74./255., 131./255., 171./255., 229./255.};
    Double_t blue[9] = { 2./255., 21./255., 35./255., 60./255., 92./255.,
                         113./255., 160./255., 185./255., 229./255.};
*/

    //avocado
/*    Double_t red[9] = { 0./255.,  4./255., 12./255.,  30./255.,  52./255.,
                            101./255., 142./255., 190./255., 237./255.};
    Double_t green[9] = { 0./255., 40./255., 86./255., 121./255., 140./255.,
                            172./255., 187./255., 213./255., 240./255.};
    Double_t blue[9] = { 0./255., 9./255., 14./255., 18./255., 21./255.,
                         23./255., 27./255., 35./255., 101./255.};
*/

    //lake
/*    Double_t red[9] = {  57./255.,  72./255.,  94./255., 117./255., 136./255.,
                             154./255., 174./255., 192./255., 215./255.};
    Double_t green[9] = {   0./255.,  33./255.,  68./255., 109./255., 140./255., 
                            171./255., 192./255., 196./255., 209./255.};
    Double_t blue[9] = { 116./255., 137./255., 173./255., 201./255., 200./255.,
                             201./255., 203./255., 190./255., 187./255.};
*/

    //sunset
    Double_t red[9]   = { 0./255., 48./255., 119./255., 173./255., 212./255.,
                         224./255., 228./255., 228./255., 245./255.};
    Double_t green[9] = { 0./255., 13./255.,  30./255.,  47./255.,  79./255.,
                         127./255., 167./255., 205./255., 245./255.};
    Double_t blue[9] = { 0./255., 68./255., 75./255., 43./255., 16./255., 
                        22./255., 55./255., 128./255., 245./255.};

    reverse(red, 9);
    reverse(green, 9);
    reverse(blue, 9);

    TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, alpha);

    skymap.Draw("SURF1Z POL");
    gPad->SetTheta(90.0);
    gPad->SetPhi(0.0);
    //gPad->SetLogz(1);
    gStyle->SetOptStat(10);

    SetRealAspectRatio(&canvas, 2);
    TText *t = new TText(-.02,0.6,"N");
    t->Draw();
    t = new TText(-.02,-.65,"S");
    t->Draw();
    t = new TText(-.68,-.02,"W");
    t->Draw();
    t = new TText(.65,-.02,"E");
    t->Draw();
    t = new TText(.007,.29,"*");
    t->SetTextSize(0.08);
    t->Draw();

/*
    TH2F bfield("fuu", "", 36, 0, 360, 18, 0, 90);
    bfield.Fill(87., 55.);
    bfield.Fill(0., 0.);
    bfield.SetMarkerStyle(29);
    bfield.SetMarkerSize(20);
    bfield.SetMarkerColor(kBlack);
    bfield.Draw("SAMES P");
*/

    gPad->Update();
    TPaletteAxis *palette = (TPaletteAxis*)skymap.GetListOfFunctions()->FindObject("palette");

    palette->SetX1NDC(0.85);
    palette->SetX2NDC(0.80);
    palette->SetY1NDC(0.05);
    palette->SetY2NDC(0.95);

    skymap.GetZaxis()->SetNdivisions(5, 0, 0, kFALSE);
    skymap.GetZaxis()->SetTitle("normalized density");
    skymap.GetXaxis()->SetTitle("foo");
    skymap.GetYaxis()->SetTitle("bar");
    skymap.GetZaxis()->SetLabelFont(42);
    skymap.GetZaxis()->SetTickLength(-0.07);
    skymap.GetZaxis()->SetLabelOffset(0.07);
    skymap.GetZaxis()->SetLabelSize(0.05);
    skymap.GetZaxis()->SetTitleSize(0.055);
    skymap.GetZaxis()->SetTitleOffset(1.7);

    gPad->Modified();
    gPad->Update();

    canvas.SaveAs("skymap_hist_pol.pdf");
}