GridUtilityScripts/CentralityDependence.C at master · ariedel-cern/GridUtilityScripts · 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
/**
 * File              : CentralityDependence.C
 * Author            : Anton Riedel <anton.riedel@tum.de>
 * Date              : 20.02.2022
 * Last Modified Date: 21.03.2022
 * Last Modified By  : Anton Riedel <anton.riedel@tum.de>
 */

#include "GridHelperMacros.H"
#include <TGraphErrors.h>
#include <boost/algorithm/string/predicate.hpp>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <nlohmann/json.hpp>
#include <string>
#include <utility>
#include <vector>

Int_t CentralityDependence(const char *FileName) {

  // load config file
  std::fstream ConfigFile("config.json");
  nlohmann::json Jconfig = nlohmann::json::parse(ConfigFile);

  // get list of all tasks
  std::vector<std::string> Tasks;
  std::vector<Double_t> CentralityBinEdges =
      Jconfig["task"]["CentralityBinEdges"].get<std::vector<Double_t>>();
  for (std::size_t i = 0; i < CentralityBinEdges.size() - 1; i++) {
    Float_t lowCentralityBinEdge = CentralityBinEdges.at(i);
    Float_t highCentralityBinEdge = CentralityBinEdges.at(i + 1);
    Tasks.push_back(std::string(Form(
        "%s_%.1f-%.1f", Jconfig["task"]["BaseName"].get<std::string>().c_str(),
        lowCentralityBinEdge, highCentralityBinEdge)));
  }

  // get list of all observables
  std::vector<std::string> Observables;
  TFile *Mean = new TFile("Bootstrap/Mean_ReTerminated.root", "READ");
  TDirectoryFile *tdirFile = dynamic_cast<TDirectoryFile *>(Mean->Get(
      Jconfig["task"]["OutputTDirectory"].get<std::string>().c_str()));
  TList *FinalResultsList = dynamic_cast<TList *>(
      dynamic_cast<TList *>(
          tdirFile->Get(tdirFile->GetListOfKeys()->First()->GetName()))
          ->FindObject("FinalResults"));

  for (auto cor : *dynamic_cast<TList *>(
           FinalResultsList->FindObject("FinalResultCorrelator"))) {
    for (auto dep : *dynamic_cast<TList *>(cor)) {
      if (boost::contains(std::string(dep->GetName()),
                          std::string("kINTEGRATED"))) {
        Observables.push_back(std::string(dep->GetName()));
      }
    }
  }
  for (auto sc : *dynamic_cast<TList *>(
           FinalResultsList->FindObject("FinalResultSymmetricCumulant"))) {
    for (auto dep : *dynamic_cast<TList *>(sc)) {
      if (boost::contains(std::string(dep->GetName()),
                          std::string("kINTEGRATED"))) {
        Observables.push_back(std::string(dep->GetName()));
      }
    }
  }

  // setup for final loop
  TFile *bootstrap = new TFile(FileName, "READ");
  TFile *out = new TFile("CentralityDependence.root", "RECREATE");
  out->cd();
  TGraphErrors *gstat, *gsys, *hstat, *hsys;
  std::vector<Double_t> x, ex, y, estat, esys;
  Int_t i;

  // create a plot for each observable
  for (auto observable : Observables) {

    std::cout << "Working on observable " << observable << std::endl;

    x.clear();
    ex.clear();
    y.clear();
    estat.clear();
    esys.clear();

    i = 0;
    for (auto task : Tasks) {

      std::cout << "in task " << task << std::endl;

      hstat = dynamic_cast<TGraphErrors *>(bootstrap->Get(
          (std::string("STAT_") + task + std::string("_") + observable)
              .c_str()));
      hsys = dynamic_cast<TGraphErrors *>(bootstrap->Get(
          (std::string("SYS_") + task + std::string("_") + observable)
              .c_str()));

      x.push_back((CentralityBinEdges[i + 1] + CentralityBinEdges[i]) / 2.);
      ex.push_back((CentralityBinEdges[i + 1] - CentralityBinEdges[i]) / 2.);
      i++;
      y.push_back(hstat->GetPointY(0));
      estat.push_back(hstat->GetErrorY(0));
      esys.push_back(hsys->GetErrorY(0));
    }

    gstat =
        new TGraphErrors(x.size(), x.data(), y.data(), ex.data(), estat.data());
    gstat->Write((std::string("STAT_") + observable).c_str());
    delete gstat;

    gsys =
        new TGraphErrors(x.size(), x.data(), y.data(), ex.data(), esys.data());
    gsys->Write((std::string("SYS_") + observable).c_str());
    delete gsys;
  }

  bootstrap->Close();
  out->Close();

  return 0;
}