epic-analysis/macro/binning.C

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Duane Byer
0003
0004 Double_t const PI = TMath::Pi();
0005
0006 Double_t const bins_x[] = {  0.01, 0.0316, 0.1, 0.316, 1. };
0007 Double_t const bins_Q2[] = { 1., 10., 100., 1000. };
0008 Double_t const bins_z[] = { 0.2, 0.4, 0.6, 0.8, 1. };
0009 Double_t const bins_pt[] = { 0., 100000. };
0010
0011 Int_t const num_bins_x = sizeof(bins_x) / sizeof(Double_t) - 1;
0012 Int_t const num_bins_Q2 = sizeof(bins_Q2) / sizeof(Double_t) - 1;
0013 Int_t const num_bins_z = sizeof(bins_z) / sizeof(Double_t) - 1;
0014 Int_t const num_bins_pt = sizeof(bins_pt) / sizeof(Double_t) - 1;
0015 Int_t const num_bins = num_bins_x * num_bins_Q2 * num_bins_z * num_bins_pt;
0016
0017 Int_t get_bin(Double_t x, Double_t const* bins, Int_t num_bins) {
0018     Int_t idx = -1;
0019     for (Int_t bin = 0; bin < num_bins + 1; ++bin) {
0020         if (x >= bins[bin]) {
0021             idx += 1;
0022         } else {
0023             return idx;
0024         }
0025     }
0026     return idx;
0027 }
0028
0029 void binning(char const* file_name, char const* output="") {
0030     std::cout << "Opening files." << std::endl;
0031     TFile* file = new TFile(file_name, "OPEN");
0032     TTree* events = file->Get<TTree>("tree");
0033     TFile* file_plots = nullptr;
0034     if (std::string(output) != std::string("")) {
0035         file_plots = new TFile(output, "RECREATE");
0036     }
0037
0038     std::cout << "Initializing histograms." << std::endl;
0039     std::vector<TH2D*> angle_up_hists;
0040     std::vector<TH2D*> angle_down_hists;
0041     std::vector<TH2D*> asym_hists;
0042     for (Int_t bin = 0; bin < num_bins; ++bin) {
0043         std::string up_name = "AngleUpHist" + std::to_string(bin);
0044         std::string down_name = "AngleDownHist" + std::to_string(bin);
0045         angle_up_hists.push_back(new TH2D(
0046             up_name.c_str(), up_name.c_str(),
0047             8, -PI, PI,
0048             8, -PI, PI));
0049         angle_down_hists.push_back(new TH2D(
0050             down_name.c_str(), down_name.c_str(),
0051             8, -PI, PI,
0052             8, -PI, PI));
0053         asym_hists.push_back(nullptr);
0054     }
0055
0056     std::cout << "Creating branch addresses." << std::endl;
0057     Double_t x;
0058     Double_t Q2;
0059     Double_t z;
0060     Double_t pt;
0061     Double_t phi_h;
0062     Double_t phi_s;
0063     Int_t spin_idx;
0064     Double_t pol;
0065     Double_t depol_1;
0066     Double_t depol_2;
0067     Double_t depol_3;
0068     Double_t depol_4;
0069     Double_t weight;
0070     events->SetBranchAddress("X", &x);
0071     events->SetBranchAddress("Q2", &Q2);
0072     events->SetBranchAddress("Z", &z);
0073     events->SetBranchAddress("PhPerp", &pt);
0074     events->SetBranchAddress("PhiH", &phi_h);
0075     events->SetBranchAddress("PhiS", &phi_s);
0076     events->SetBranchAddress("Spin_idx", &spin_idx);
0077     events->SetBranchAddress("Pol", &pol);
0078     events->SetBranchAddress("Depol1", &depol_1);
0079     events->SetBranchAddress("Depol2", &depol_2);
0080     events->SetBranchAddress("Depol3", &depol_3);
0081     events->SetBranchAddress("Depol4", &depol_4);
0082     events->SetBranchAddress("Weight", &weight);
0083
0084     std::cout << "Looping over events." << std::endl;
0085     Int_t num_events = events->GetEntries();
0086     for (Int_t event = 0; event < num_events; ++event) {
0087         events->GetEntry(event);
0088         Int_t bin_x = get_bin(x, bins_x, num_bins_x);
0089         Int_t bin_Q2 = get_bin(Q2, bins_Q2, num_bins_Q2);
0090         Int_t bin_z = get_bin(z, bins_z, num_bins_z);
0091         Int_t bin_pt = get_bin(pt, bins_pt, num_bins_pt);
0092         if (!(bin_x >= 0 && bin_x < num_bins_x)
0093                 || !(bin_Q2 >= 0 && bin_Q2 < num_bins_Q2)
0094                 || !(bin_z >= 0 && bin_z < num_bins_z)
0095                 || !(bin_pt >= 0 && bin_pt < num_bins_pt)) {
0096             continue;
0097         }
0098         Int_t bin =
0099             bin_x + num_bins_x * (
0100             bin_Q2 + num_bins_Q2 * (
0101             bin_z + num_bins_z * (
0102             bin_pt)));
0103         Double_t S = 820.;
0104         Double_t y = Q2 / (S * x);
0105         if (y > 0.05) {
0106             if (spin_idx == 1) {
0107                 angle_up_hists[bin]->Fill(phi_h, phi_s, weight);
0108             } else if (spin_idx == -1) {
0109                 angle_down_hists[bin]->Fill(phi_h, phi_s, weight);
0110             } else {
0111                 // Shouldn't happen.
0112             }
0113         }
0114     }
0115
0116     std::cout << "Fitting histograms." << std::endl;
0117     std::vector<std::vector<Double_t> > params{ {}, {}, {}, {}, {} };
0118     std::vector<std::vector<Double_t> > param_errs{ {}, {}, {}, {}, {} };
0119     std::vector<std::string> param_names{ "Sivers", "Collins", "sin(3φ_h - φ_s)", "sin(φ_s)", "sin(2φ_h - φ_s)" };
0120     TF2* fit = new TF2(
0121         "asymmetry",
0122         "[0] * ([3] * sin(x - y) + [1] * [4] * sin(x + y) + [1] * [5] * sin(3 * x - y) + [2] * [6] * sin(y) + [2] * [7] * sin(2 * x - y))",
0123         -PI, PI,
0124         -PI, PI);
0125     fit->FixParameter(0, pol);
0126     // Deal with depolarization later.
0127     //fit->FixParameter(1, depol_1);
0128     fit->FixParameter(1, 1.);
0129     fit->FixParameter(2, 1.);
0130     for (std::size_t p = 0; p < params.size(); ++p) {
0131         fit->SetParameter(3 + p, 0.);
0132         fit->SetParLimits(3 + p, -1., 1.);
0133     }
0134     for (Int_t bin = 0; bin < num_bins; ++bin) {
0135         Int_t bin_x = bin % num_bins_x;
0136         Int_t bin_Q2 = (bin / num_bins_x) % num_bins_Q2;
0137         Int_t bin_z = (bin / (num_bins_x * num_bins_Q2)) % num_bins_z;
0138         Int_t bin_pt = (bin / (num_bins_x * num_bins_Q2 * num_bins_pt)) % num_bins_pt;
0139         Int_t count = angle_up_hists[bin]->GetEntries();
0140         asym_hists[bin] = static_cast<TH2D*>(angle_up_hists[bin]->GetAsymmetry(angle_down_hists[bin]));
0141         std::string asym_name = "AsymmetryHist" + std::to_string(bin);
0142         asym_hists[bin]->SetTitle(asym_name.c_str());
0143         asym_hists[bin]->SetName(asym_name.c_str());
0144         asym_hists[bin]->Fit(fit, "IDMQ");
0145         for (std::size_t p = 0; p < params.size(); ++p) {
0146             params[p].push_back(0.);
0147             param_errs[p].push_back(1.);
0148         }
0149         if (count != 0) {
0150             for (std::size_t p = 0; p < params.size(); ++p) {
0151                 params[p].back() = fit->GetParameter(3 + p);
0152                 param_errs[p].back() = fit->GetParError(3 + p);
0153             }
0154         }
0155         std::cout << "Bin " << bin << ":" << std::endl
0156             << "\tCoords:  " << bin_x << ", " << bin_Q2 << ", " << bin_z << ", " << bin_pt << std::endl
0157             << "\tCount:   " << count << std::endl;
0158         for (std::size_t p = 0; p < params.size(); ++p) {
0159             std::cout << param_names[p] << ": " << params[p].back() << " ± " << param_errs[p].back() << std::endl;
0160         }
0161         if (file_plots != nullptr) {
0162             angle_up_hists[bin]->Write();
0163             angle_down_hists[bin]->Write();
0164             asym_hists[bin]->Write();
0165         }
0166     }
0167     TList* axes = new TList();
0168     std::vector<std::string> axis_names;
0169     if (num_bins_x > 1) {
0170         axes->Add(new TAxis(num_bins_x, bins_x));
0171         axis_names.push_back("x");
0172     }
0173     if (num_bins_Q2 > 1) {
0174         axes->Add(new TAxis(num_bins_Q2, bins_Q2));
0175         axis_names.push_back("Q sq.");
0176     }
0177     if (num_bins_z > 1) {
0178         axes->Add(new TAxis(num_bins_z, bins_z));
0179         axis_names.push_back("z");
0180     }
0181     if (num_bins_pt > 1) {
0182         axes->Add(new TAxis(num_bins_pt, bins_pt));
0183         axis_names.push_back("pt");
0184     }
0185     if (file_plots != nullptr) {
0186         file_plots->WriteObject(axes, "Axes");
0187         file_plots->WriteObject(&axis_names, "AxisNames");
0188         file_plots->WriteObject(&params, "Params");
0189         file_plots->WriteObject(&param_names, "ParamNames");
0190         file_plots->WriteObject(&param_errs, "ParamErrs");
0191     }
0192     if (file_plots != nullptr) {
0193         file_plots->Close();
0194     }
0195     file->Close();
0196 }
0197