Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-05-23 08:11:26 

0001 #include <TFile.h>
0002 #include <TH1D.h>
0003 #include <TH2D.h>
0004 #include <TKey.h>
0005 #include <TList.h>
0006 #include <fstream>
0007 #include <iostream>
0008 #include <string>
0009 #include <vector>
0010 
0011 #include "BootstrapGenerator/TH1DBootstrap.h"
0012 #include "BootstrapGenerator/TH2DBootstrap.h"
0013 
0014 R__LOAD_LIBRARY(libBootstrapGenerator.so)
0015 
0016 void read_luminosity(const std::string& filename, double& lumi_all, double& lumi_zvertex60) {
0017   std::ifstream infile(filename);
0018   if (!infile.is_open()) {
0019     std::cerr << "Error: Cannot open " << filename << std::endl;
0020     lumi_all = 0;
0021     lumi_zvertex60 = 0;
0022     return;
0023   }
0024   double dummy1;
0025   infile >> dummy1 >> lumi_zvertex60 >> lumi_all;
0026   infile.close();
0027   std::cout << "Read from " << filename << ": lumi_all=" << lumi_all << ", lumi_zvertex60=" << lumi_zvertex60 << std::endl;
0028 }
0029 
0030 void combine_one(const std::string& pattern, int mode, int jet_radius_index, double frac_0mrad_all, double frac_1p5mrad_all, double frac_0mrad_zvertex60, double frac_1p5mrad_zvertex60) {
0031   
0032   std::vector<std::string> filenames;
0033   std::vector<double> weights_all;
0034   std::vector<double> weights_zvertex60;
0035   std::string outname;
0036   
0037   if (mode == 1) {
0038     outname = Form("output_sim_hadd/output_sim_r0%d_%s.root", jet_radius_index, pattern.c_str());
0039     filenames.push_back(Form("output_sim_hadd/output_sim_0mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0040     filenames.push_back(Form("output_sim_hadd/output_sim_1p5mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0041     weights_all = {frac_0mrad_all, frac_1p5mrad_all};
0042     weights_zvertex60 = {frac_0mrad_zvertex60, frac_1p5mrad_zvertex60};
0043   }
0044   else if (mode == 2) {
0045     outname = Form("output_sim_hadd/output_sim_r0%d_%s.root", jet_radius_index, pattern.c_str());
0046     filenames.push_back(Form("output_sim_hadd/output_sim_0mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0047     filenames.push_back(Form("output_sim_hadd/output_sim_0mrad_double_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0048     weights_all = {0.78, 0.22};
0049     weights_zvertex60 = {0.78, 0.22};
0050   }
0051   else if (mode == 3) {
0052     outname = Form("output_sim_hadd/output_sim_r0%d_%s.root", jet_radius_index, pattern.c_str());
0053     filenames.push_back(Form("output_sim_hadd/output_sim_1p5mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0054     filenames.push_back(Form("output_sim_hadd/output_sim_1p5mrad_double_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0055     weights_all = {0.921, 0.079};
0056     weights_zvertex60 = {0.921, 0.079};
0057   }
0058   else if (mode == 4) {
0059     outname = Form("output_sim_hadd/output_sim_r0%d_%s.root", jet_radius_index, pattern.c_str());
0060     filenames.push_back(Form("output_sim_hadd/output_sim_0mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0061     filenames.push_back(Form("output_sim_hadd/output_sim_0mrad_double_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0062     filenames.push_back(Form("output_sim_hadd/output_sim_1p5mrad_single_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0063     filenames.push_back(Form("output_sim_hadd/output_sim_1p5mrad_double_r0%d_%s.root", jet_radius_index, pattern.c_str()));
0064     weights_all = {0.78*frac_0mrad_all, 0.22*frac_0mrad_all, 0.921*frac_1p5mrad_all, 0.079*frac_1p5mrad_all};
0065     weights_zvertex60 = {0.78*frac_0mrad_zvertex60, 0.22*frac_0mrad_zvertex60, 0.921*frac_1p5mrad_zvertex60, 0.079*frac_1p5mrad_zvertex60};
0066   }
0067   
0068   // Open files
0069   std::vector<TFile*> files;
0070   for (const auto& fname : filenames) {
0071     TFile* f = new TFile(fname.c_str(), "READ");
0072     if (!f || f->IsZombie()) {
0073       std::cerr << "Error: Cannot open " << fname << std::endl;
0074       for (auto ff : files) { ff->Close(); delete ff; }
0075       return;
0076     }
0077     files.push_back(f);
0078   }
0079   
0080   TFile* f_out = new TFile(outname.c_str(), "RECREATE");
0081   
0082   // Loop over histograms
0083   TList* keys = files[0]->GetListOfKeys();
0084   TIter next(keys);
0085   TKey* key;
0086   
0087   while ((key = (TKey*)next())) {
0088     std::string histname = key->GetName();
0089     std::string classname = key->GetClassName();
0090     
0091     // Choose weights based on histogram name
0092     bool use_zvertex60 = (histname.find("zvertex60") != std::string::npos);
0093     const std::vector<double>& weights = use_zvertex60 ? weights_zvertex60 : weights_all;
0094     
0095     if (classname == "TH1D") {
0096       TH1D* h_combined = nullptr;
0097       for (int i = 0; i < files.size(); i++) {
0098         TH1D* h = (TH1D*)files[i]->Get(histname.c_str());
0099         if (!h) continue;
0100         if (!h_combined) {
0101           f_out->cd();
0102           h_combined = new TH1D(*h);
0103           h_combined->SetName(histname.c_str());
0104           h_combined->Scale(weights[i]);
0105         } else {
0106           h_combined->Add(h, weights[i]);
0107         }
0108       }
0109       if (h_combined) h_combined->Write();
0110     }
0111     else if (classname == "TH2D") {
0112       TH2D* h_combined = nullptr;
0113       for (int i = 0; i < files.size(); i++) {
0114         TH2D* h = (TH2D*)files[i]->Get(histname.c_str());
0115         if (!h) continue;
0116         if (!h_combined) {
0117           f_out->cd();
0118           h_combined = new TH2D(*h);
0119           h_combined->SetName(histname.c_str());
0120           h_combined->Scale(weights[i]);
0121         } else {
0122           h_combined->Add(h, weights[i]);
0123         }
0124       }
0125       if (h_combined) h_combined->Write();
0126     }
0127     else if (classname == "TH1DBootstrap") {
0128       TH1DBootstrap* h_combined = nullptr;
0129       for (int i = 0; i < files.size(); i++) {
0130         TH1DBootstrap* h = (TH1DBootstrap*)files[i]->Get(histname.c_str());
0131         if (!h) continue;
0132         if (!h_combined) {
0133           f_out->cd();
0134           h_combined = new TH1DBootstrap(*h);
0135           h_combined->SetName(histname.c_str());
0136           h_combined->Scale(weights[i]);
0137         } else {
0138           TH1DBootstrap h_temp(*h);
0139           h_temp.Scale(weights[i]);
0140           h_combined->Add(&h_temp);
0141         }
0142       }
0143       if (h_combined) h_combined->Write();
0144     }
0145     else if (classname == "TH2DBootstrap") {
0146       TH2DBootstrap* h_combined = nullptr;
0147       for (int i = 0; i < files.size(); i++) {
0148         TH2DBootstrap* h = (TH2DBootstrap*)files[i]->Get(histname.c_str());
0149         if (!h) continue;
0150         if (!h_combined) {
0151           f_out->cd();
0152           h_combined = new TH2DBootstrap(*h);
0153           h_combined->SetName(histname.c_str());
0154           h_combined->Scale(weights[i]);
0155         } else {
0156           TH2DBootstrap h_temp(*h);
0157           h_temp.Scale(weights[i]);
0158           h_combined->Add(&h_temp);
0159         }
0160       }
0161       if (h_combined) h_combined->Write();
0162     }
0163   }
0164   
0165   f_out->Close();
0166   for (auto f : files) { f->Close(); delete f; }
0167   delete f_out;
0168   
0169   std::cout << "Output: " << outname << std::endl;
0170 }
0171 
0172 void get_mixoutput(int mode, int jet_radius_index = 4) {
0173   // mode 1: 0mrad single + 1p5mrad single
0174   // mode 2: 0mrad single + 0mrad double
0175   // mode 3: 1p5mrad single + 1p5mrad double
0176   // mode 4: 0mrad single + 0mrad double + 1p5mrad single + 1p5mrad double
0177   
0178   double lumi_0mrad_all, lumi_1p5mrad_all;
0179   double lumi_0mrad_zvertex60, lumi_1p5mrad_zvertex60;
0180   read_luminosity("luminosity_0mrad.txt", lumi_0mrad_all, lumi_0mrad_zvertex60);
0181   read_luminosity("luminosity_1p5mrad.txt", lumi_1p5mrad_all, lumi_1p5mrad_zvertex60);
0182   
0183   double frac_0mrad_all = lumi_0mrad_all / (lumi_0mrad_all + lumi_1p5mrad_all);
0184   double frac_1p5mrad_all = lumi_1p5mrad_all / (lumi_0mrad_all + lumi_1p5mrad_all);
0185 
0186   double frac_0mrad_zvertex60 = lumi_0mrad_zvertex60 / (lumi_0mrad_zvertex60 + lumi_1p5mrad_zvertex60);
0187   double frac_1p5mrad_zvertex60 = lumi_1p5mrad_zvertex60 / (lumi_0mrad_zvertex60 + lumi_1p5mrad_zvertex60);
0188   
0189   std::cout << "Lumi fractions: 0mrad=" << frac_0mrad_all << ", 1p5mrad=" << frac_1p5mrad_all << std::endl;
0190   std::cout << "Lumi fractions (zvertex60): 0mrad=" << frac_0mrad_zvertex60 << ", 1p5mrad=" << frac_1p5mrad_zvertex60 << std::endl;
0191   
0192   std::vector<std::string> patterns = {"Jet12GeV", "Jet20GeV", "Jet30GeV", "Jet40GeV", "Jet50GeV"};
0193   
0194   for (const auto& pattern : patterns) {
0195     std::cout << "\nCombining " << pattern << "..." << std::endl;
0196     combine_one(pattern, mode, jet_radius_index, frac_0mrad_all, frac_1p5mrad_all, frac_0mrad_zvertex60, frac_1p5mrad_zvertex60);
0197   }
0198   
0199   std::cout << "\nAll done!" << std::endl;
0200 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!