Skip to content
Snippets Groups Projects
Select Git revision
  • shuo
  • master default protected
  • derp
3 results

plot_hms_singles.cxx

  • plot_hms_singles.cxx 4.71 KiB
    #ifdef __cpp_lib_filesystem
    #include <filesystem>
    namespace fs = std::filesystem;
    #else
    #include <experimental/filesystem>
    namespace fs = std::experimental::filesystem;
    #endif
    
    #include "nlohmann/json.hpp"
    
    #include "TObject.h"
    #include "TGraph.h"
    
    
    void plot_hms_singles(int start_run = 0){
    
      using nlohmann::json;
      json j;
      {
        std::ifstream json_input_file("db2/run_count_list.json");
        try {
          json_input_file >> j;
        } catch(json::parse_error) {
          std::cerr << "error: json file, db2/run_count_list.json, is incomplete or has broken syntax.\n";
          std::quick_exit(-127);
        }
      }
      json j2;
      {
        std::ifstream json_input_file("db2/run_list.json");
        try {
          json_input_file >> j2;
        } catch(json::parse_error)  {
          std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
          std::quick_exit(-127);
        }
      }
    
    
      std::vector<double> runs;
      std::vector<double> hms_yields;
      std::vector<double> hms_yield_uncs;
    
      std::vector<double> runs2;
      std::vector<double> hms_yields2;
      std::vector<double> hms_yield_uncs2;
    
      std::vector<double> runs3;
      std::vector<double> yield_ratios;
      std::vector<double> yield_ratio_uncs;
    
      std::vector<double> target_ids;
    
      TGraphErrors* gr_hms_count  = new TGraphErrors();
      TGraphErrors* gr_hms_scaler = new TGraphErrors();
      TGraphErrors* gr_yield_ratio = new TGraphErrors();
    
    
      double scaler_reduce_factor = 0.4;
    
      for (json::iterator it = j.begin(); it != j.end(); ++it) {
    
        if (j2.find(it.key()) == j2.end()) {
          continue;
        }
        if (j2[it.key()]["target"]["target_id"].get<int>() != 3) {
          continue;
        }
    
        int run_number = std::stoi(it.key());
        int run_number2 = 0;
        double yield      = 0.0;
        double yield2      =0.0;
    
        if(run_number < start_run ) {
          continue;
        }
    
        if (it.value().find("hms e yield") != it.value().end()) {
          try {
          yield      = it.value()["hms e yield"].get<double>();
          runs.push_back(double(run_number));
          hms_yields.push_back(yield);
          if (it.value().find("hms e yield unc.") != it.value().end()) {
            hms_yield_uncs.push_back(it.value()["hms e yield unc."].get<double>());
          } else {
            hms_yield_uncs.push_back(0.0);
          }
          } catch(std::domain_error ) {
            continue;
            //you suck
          }
        }
        if (it.value().find("hms scaler yield") != it.value().end()) {
          run_number2 = std::stoi(it.key());
          yield2      = it.value()["hms scaler yield"].get<double>();
          runs2.push_back(double(run_number2));
          hms_yields2.push_back(scaler_reduce_factor*yield2);
          if (it.value().find("hms scaler yield unc.") != it.value().end()) {
            hms_yield_uncs2.push_back(scaler_reduce_factor*it.value()["hms scaler yield unc."].get<double>());
          } else {
            hms_yield_uncs2.push_back(0.0);
          }
        }
        if (it.value().find("hms scaler yield") != it.value().end()) {
          if (it.value().find("hms e yield") != it.value().end()) {
            runs3.push_back(double(run_number2));
            double yield_ratio = yield / yield2;
            yield_ratios.push_back(yield_ratio);
            if ((it.value().find("hms e yield unc.") != it.value().end()) &&
                (it.value().find("hms scaler yield unc.") != it.value().end())) {
              double y1 = it.value()["hms e yield unc."].get<double>();
              double y2 = it.value()["hms scaler yield unc."].get<double>();
              yield_ratio_uncs.push_back(yield_ratio*std::sqrt(y1 * y1/(yield*yield) + y2 * y2/(yield2*yield2)));
            } else {
              yield_ratio_uncs.push_back(0.0);
            }
          }
        }
      }
    
      // ------------------------------
      //
      TGraphErrors* gr = nullptr;
      TMultiGraph*  mg = nullptr;
      TCanvas*      c  = nullptr;
    
      // ------------------------------
      //
      c  = new TCanvas("c1","c1",1600,1200);
      c->Divide(1, 2);
    
      c->cd(1);
      mg = new TMultiGraph();
      mg->SetTitle("; run number ; counts/mC");
    
      std::vector<double> zeros(runs.size());
      gr = new TGraphErrors(runs.size(), runs.data(), hms_yields.data(),zeros.data(),hms_yield_uncs.data());
      gr->SetMarkerStyle(20);
      mg->Add(gr,"p");
    
      std::vector<double> zeros2(runs2.size());
      gr = new TGraphErrors(runs2.size(), runs2.data(), hms_yields2.data(),zeros2.data(),hms_yield_uncs2.data());
      gr->SetMarkerStyle(23);
      gr->SetMarkerColor(2);
      mg->Add(gr,"p");
    
      mg->Draw("a");
    
      c->cd(2);
      std::vector<double> zeros3(runs3.size());
      gr = new TGraphErrors(runs3.size(), runs3.data(), yield_ratios.data(),zeros3.data(),yield_ratio_uncs.data());
      gr->SetTitle("; run number ; ps4*T2 + T6 counts / scaler");
      gr->SetMarkerStyle(20);
      //gr->Draw("alp");
      //gr->GetYaxis()->SetRangeUser(0.3, 0.45);
      gr->Draw("alp");
    
      c->SaveAs("results/hms_singles_vs_run_number.png");
      c->SaveAs("results/hms_singles_vs_run_number.pdf");
    }