misc/plotPostNuisance_combine.C at master · AndrewLevin/misc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
// originally from https://github.com/GuillelmoGomezCeballos/Combination/blob/master/plotPostNuisance_combine.C
//
// Script that visualizes the post-fit nuisance parameters
//
// Usage :
//  root -b plotPostNuisance_combine.C'("mlfit.root", true)'
//
// Arguments :
//  * char* mlfit : name of the root file returned by combine
//                  This can be generated by the following command
//                  $ combine -M MaxLikelihoodFit --saveNorm [card]
//  * bool skipmu :      flag to decide whether to show mu or not
//

void printaline(int index, char* nuiname, Double_t centralval, Double_t uncert){

    cout.width(5);  cout << index << " ::: ";
    cout.width(40); cout << nuiname << " : ";
    cout.width(10);  cout << Form("%1.3f", centralval);
    cout.width(5);  cout << " +/- ";
    cout.width(10);  cout << Form("%1.3f", uncert) << endl;

}

void plotPostNuisance_combine(char* mlfit="mlfit.root", double MaxUncertainty = 0.9,
                              double MaxValue = 1.0, bool skipmu = true) {

    //
    // Read a mlfit.root file and the list of post-fit nuisance parameters
    // from S+B fit (fit_s)
    //
    TFile *File = TFile::Open(mlfit, "READ");
    RooFitResult *fit_s = (RooFitResult*)File->Get("fit_s");
    if(!fit_s)  return;
    //if(fit_s->status() != 3) { cout << " MINUIT status is not 3 "  << endl; continue; }
    RooArgList parlist = fit_s->floatParsFinal();

    //
    // Containers for post-fit nuisances
    //
    std::vector<TString> nuisance;  // nuisance name
    std::vector<Double_t> central;  // post-fit central value  : normalized to the input uncertainty
    std::vector<Double_t> uncert;   // post-fit uncertainty    : normalized to the input uncertainty

    //
    // Loop over RooArgList and store the fit results
    // in the above containers
    //
    cout << "MaxUncertainty: " << MaxUncertainty << endl;
    cout << "-------------------------------------------------------------------------------" << endl;
    cout << "Index :::               Nuisance name              :       Pull +/- uncertainty" << endl;
    cout << "-------------------------------------------------------------------------------" << endl;
    for(int i=0; i<parlist.getSize(); i++) {

        TString name_tstr = parlist[i].GetName();

        if( skipmu && name_tstr == "r" ) continue;
        if(((RooRealVar&)parlist[i]).getError() < MaxUncertainty)
        printaline( i+1, parlist[i].GetName(),  // index here is same as plot (i+1), not parlist (i)
                    ((RooRealVar&)parlist[i]).getVal(),
                    ((RooRealVar&)parlist[i]).getError());

        central.push_back( ((RooRealVar&)parlist[i]).getVal() );
        uncert.push_back( ((RooRealVar&)parlist[i]).getError() );
        nuisance.push_back( parlist[i].GetName() );

    }
    cout << "-------------------------------------------------------------------------------" << endl;

    //
    // Visualize the result
    //
    // * gray area corresponds to post-fit cenral value and pre-fit uncertainty
    // * solid lines correspond to post-fit central value and post-fit uncertainty
    // Purpose of gray area is to compare pre/post-fit uncertainty and
    // check how better data constrains a nuisance than prediction
    TCanvas *c = new TCanvas("c", "c", 800, 500);
    //c->Divide(1,2);

    // store nuisances with pull >= 1
    // value is set in line 76 : if(TMath::Abs(central[i])>=1)
    int nchangednuisance = 0;
    std::vector<TString> changednuisance;

    c->cd(1);
    int nbins = central.size()+2;
    TH1D* h1        =   new TH1D("h1",      "h1",       nbins, -0.5, nbins-0.5); // post-fit central value and post-fit uncertainty
    TH1D* h1_ref    =   new TH1D("h1_ref",  "h1_ref",   nbins, -0.5, nbins-0.5); // post-fit central value and pre-fit uncertainty

    cout << endl;
    cout << " ****************************** LARGE VARIATION ****************************** " << endl;
    for(int i=0; i<nbins-2; i++) {
        h1->SetBinContent(i+2, central[i]);
        h1->SetBinError(i+2, uncert[i]);
        h1_ref->SetBinContent(i+2, central[i]);
        h1_ref->SetBinError(i+2, 1.);
        if(TMath::Abs(central[i]) >= MaxValue) {
            nchangednuisance++;
            changednuisance.push_back(Form("%2i ::: %s : %1.3f +/- %1.3f", i+1, nuisance[i].Data(), central[i], uncert[i]));

            printaline( i+1, nuisance[i], central[i], uncert[i]);
        }
    }
    cout << " ****************************************************************************** " << endl;
    cout << endl;
    h1_ref->SetStats(0);
    h1_ref->SetTitle("");
    h1_ref->SetYTitle("Pull");
    h1_ref->SetFillColor(18);
    h1_ref->SetMinimum(-2.);
    h1_ref->SetMaximum(2.);
    //h1_ref->SetFillStyle(3004);
    h1->SetLineColor(kBlack);

    // Draw a horizontal line at 0
    TH1D* h1_zero    =   new TH1D("h1_zero", "h1_zero", 1, -0.5, nbins-0.5);
    h1_zero->SetBinContent(1,0);
    h1_zero->SetLineStyle(2);

    h1_ref->Draw("E2");
    h1_zero->Draw("histo same");
    h1->Draw("same E");

    TString plotname = mlfit;
    plotname.ReplaceAll(".root", "_Postnuisance.pdf");
    //c->SaveAs(Form("/afs/cern.ch/user/j/jaehyeok/www/tmp/%s", plotname.Data()));
    c->SaveAs(Form("%s", plotname.Data()));

/*
    // nuisances with large change
    c->cd(2);
    TLatex l;
    l.DrawLatex(0.1, 0.9, "Nuicance : pull +/- postfit uncertainty");
    l.SetTextSize(0.04);

    for(int i=0; i<nchangednuisance; i++) {
        l.DrawLatex(0.1, 0.7-0.08*i, changednuisance[i]);
    }
*/

}