draw_differential_binding.C

#include "tkmanager.h"
#include "TGraph.h"
#include "TLine.h"
#include "TF1.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TROOT.h"
#include "TLatex.h"
 
using namespace tkn;
using namespace std;
 
int mn[]={2,8,20,28,50,82};
 
Double_t fitfunc(Double_t *x,Double_t *par);
TLine* get_line(unsigned int MagicNumber);
TLatex* get_shell(unsigned int MagicNumber, const char* shell);
TLegend* get_legend();
void set_style(TGraph* g1, TGraph* g2);
 
void draw_differential_binding()
{
    gROOT->SetStyle("tkn-histo");
 
    TGraph* nuc_neqz = new TGraph;
    nuc_neqz->SetTitle("N = Z");
    TGraph* nuc_other = new TGraph;
    nuc_other->SetTitle("N #neq Z");
 
    set_style(nuc_neqz, nuc_other);
 
    tkmanager dtm;
    for(const auto &nuc : dtm.get_nuclei([&](auto nuc) {return (nuc->get_z()%2==0 &&
                                                               nuc->get_n()%2==0 &&
                                                               nuc->get_z()<50&&
                                                               dtm.known_nucleus(nuc->get_z()-2, nuc->get_a()-2)&&
                                                               dtm.known_nucleus(nuc->get_z(), nuc->get_a()-2)&&
                                                               dtm.known_nucleus(nuc->get_z()-2, nuc->get_a()-4));}))
    {
        int Z = nuc->get_z();
        int A = nuc->get_a();
 
        double binding_diff = -0.25*(nuc->get_binding_energy_over_a()*A
                                   - dtm.get_nucleus(Z-2, A-2)->get_binding_energy_over_a()*(A-2)
                                   - dtm.get_nucleus(Z,A-2)->get_binding_energy_over_a()*(A-2)
                                   + dtm.get_nucleus(Z-2,A-4)->get_binding_energy_over_a()*(A-4));
 
        if(Z==A-Z) nuc_neqz->AddPoint(Z, binding_diff);
        else nuc_other->AddPoint(Z, binding_diff);
    }
 
    TF1* func = new TF1("fit",fitfunc,2,50,5);
    func->SetParameters(-8,-4.1,-2,-1.5,-1);
    func->SetMarkerStyle(1);
    func->SetMarkerColor(func->GetLineColor());
    nuc_neqz->Fit(func,"QN","",2,50);
 
    TH2* hdum = new TH2F("dum","dum",10,2,50,10,-8,0);
    hdum->GetYaxis()->SetTitle("#Delta^{2}B (MeV)");
    hdum->GetXaxis()->SetTitle("Proton Number");
    hdum->Draw();
 
    for(auto z: {1,2,3}) get_line(mn[z])->Draw();
 
    const char* shells[] = {"1p","2s/1d","1f","1g,1f,2p"};
    for(auto z: {1,2,3,4}) get_shell(mn[z],shells[z-1])->Draw();
 
    func->Draw("same");
    nuc_neqz->Draw("P");
    nuc_other->Draw("P");
    get_legend()->Draw();
}
 
Double_t fitfunc(Double_t *x,Double_t *par)
{
    if(x[0]<mn[0]||x[0]>mn[4]) return std::nan("1");
    int idx = 0;
    for(; idx<6; idx++) if(x[0]<=mn[idx]) break;
    return (par[idx]-par[idx-1])/(mn[idx]-mn[idx-1])*(x[0]-mn[idx-1]) + par[idx-1];
}
 
TLine* get_line(unsigned int MagicNumber)
{
    TLine *ll = new TLine(MagicNumber,-8+0.01,MagicNumber,-0.01);
    ll->SetLineColor(kGray+1);
    return ll;
}
 
TLatex* get_shell(unsigned int MagicNumber, const char* shell)
{
    TLatex* tex = new TLatex(MagicNumber-1,-7.5,shell);
    tex->SetTextAlign(32);
    tex->SetTextFont(132);
    tex->SetTextColor(kGray+2);
    tex->SetTextSize(0.0630252);
    tex->SetLineWidth(2);
    return tex;
}
 
TLegend *get_legend()
{
    TLegend *leg = new TLegend(0.747492,0.432773,0.956522,0.64916,NULL,"brNDC");
    leg->SetBorderSize(0);
    leg->SetTextFont(132);
    leg->SetLineColor(0);
    leg->SetLineStyle(0);
    leg->SetLineWidth(0);
    leg->SetFillColor(0);
    leg->SetFillStyle(1001);
 
    TLegendEntry *entry = 0;
    entry=leg->AddEntry("","N #neq Z","lpf");
    entry->SetFillStyle(1000);
    entry->SetLineStyle(0);
    entry->SetLineColor(kWhite);
    entry->SetLineWidth(0);
    entry->SetMarkerColor(kAzure-3);
    entry->SetMarkerStyle(20);
    entry->SetMarkerSize(1);
    entry->SetTextFont(132);
 
    entry=leg->AddEntry("","N = Z","lpf");
    entry->SetFillStyle(1000);
    entry->SetLineStyle(0);
    entry->SetLineWidth(0);
    entry->SetLineColor(kWhite);
    entry->SetMarkerColor(kCyan-2);
    entry->SetMarkerStyle(21);
    entry->SetMarkerSize(1.2);
    entry->SetTextFont(132);
 
    return leg;
}
 
void set_style(TGraph* g1, TGraph* g2)
{
    g1->SetMarkerStyle(21);
    g1->SetMarkerSize(1.2);
    g1->SetMarkerColor(kCyan-2);
    g1->SetLineColor(kWhite);
    g1->SetTitle("N = Z");
 
    g2->SetMarkerStyle(20);
    g2->SetMarkerSize(1);
    g2->SetMarkerColor(kAzure-3);
    g2->SetLineColor(kWhite);
    g2->SetTitle("N #neq Z");
}