bunch_shuffling/run/BunchShuffling.cc

0001 #include <iostream>
0002 #include <sstream>
0003 #include <cstdlib>
0004 #include <cmath>
0005
0006 #include "AsymmetryCalc/WrapperFinalAsymArrays.h"
0007 #include "AsymmetryCalc/Constants.h"
0008 #include "AsymmetryCalc/ShuffleBunches.h"
0009
0010 // To check virtual and resident memory usage
0011 #include <Riostream.h>
0012 #include <TSystem.h>
0013 #include <malloc.h>
0014
0015 void monitorMemoryUsage(const std::string &label)
0016 {
0017   ProcInfo_t procInfo;
0018   gSystem->GetProcInfo(&procInfo);
0019   std::cout << label << " - Memory usage: "
0020             << "Resident = " << procInfo.fMemResident << " kB, "
0021             << "Virtual = " << procInfo.fMemVirtual << " kB" << std::endl;
0022 }
0023
0024 std::unordered_map<std::string, std::string> read_config(const std::string& filename)
0025 {
0026   std::unordered_map<std::string, std::string> config;
0027   std::ifstream file(filename);
0028   if (!file.is_open())
0029   {
0030     std::cerr << "Warning: Config file " << filename << " not found. All config paths will be set to the default." << std::endl;
0031     return config;
0032   }
0033
0034   std::string line;
0035   while (std::getline(file, line))
0036   {
0037     std::istringstream is_line(line);
0038     std::string key;
0039     if (std::getline(is_line, key, '='))
0040     {
0041       std::string value;
0042       if (std::getline(is_line, value, '='))
0043       {
0044         config[key] = value;
0045       }
0046     }
0047   }
0048   return config;
0049 }
0050
0051 void book_shuffle_file(TFile* (&output_file_shuffle),
0052                        TH1F* (&h_AN_pt_shuffle)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections],
0053                        TH1F* (&h_AN_pt_nodir_shuffle)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins],
0054                        TH1F* (&h_AN_eta_shuffle)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins],
0055                        TH1F* (&h_AN_xf_shuffle)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins],
0056                        TH1F* (&h_AN_pt_shuffle_details)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0057                        TH1F* (&h_AN_pt_nodir_shuffle_details)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0058                        TH1F* (&h_AN_eta_shuffle_details)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0059                        TH1F* (&h_AN_xf_shuffle_details)[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0060                        const std::string &outputfilename)
0061 {
0062   output_file_shuffle = new TFile(outputfilename.c_str(), "RECREATE");
0063   output_file_shuffle->cd();
0064   for (int iP = 0; iP < ASYM_CONSTANTS::nParticles; iP++) {
0065     for (int iPt = 0; iPt < ASYM_CONSTANTS::nPtBins; iPt++) {
0066       {
0067         std::stringstream hname;
0068         hname << "h_AN_pt_nodir_shuffle_"
0069             << ASYM_CONSTANTS::particle[iP] << "_pt_"
0070             << iPt << "_sqrt";
0071         std::stringstream htitle;
0072         htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0073         h_AN_pt_nodir_shuffle[iP][iPt] = new TH1F(
0074           hname.str().c_str(),
0075           htitle.str().c_str(),
0076           100, -4, 4);
0077         // Beam- and config- specific values
0078         for (int iB = 0; iB < ASYM_CONSTANTS::nBeams; iB++) {
0079           for (int iC = 0; iC < ASYM_CONSTANTS::nConfigs; iC++) {
0080             hname.str("");
0081             hname << "h_AN_pt_nodir_shuffle_"
0082                   << ASYM_CONSTANTS::beams[iB] << "_"
0083                   << ASYM_CONSTANTS::particle[iP] << "_pt_"
0084                   << iPt << "_"
0085                   << ASYM_CONSTANTS::configuration[iC] << "_sqrt";
0086             htitle.str("");
0087             htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0088             h_AN_pt_nodir_shuffle_details[iP][iPt][iB][iC] = new TH1F(
0089               hname.str().c_str(),
0090               htitle.str().c_str(),
0091               100, -4, 4);
0092           }
0093         }
0094       }
0095       for (int iDir = 0; iDir < ASYM_CONSTANTS::nDirections; iDir++) {
0096         std::stringstream hname;
0097         hname << "h_AN_pt_shuffle_"
0098             << ASYM_CONSTANTS::particle[iP] << "_"
0099             << ASYM_CONSTANTS::directions[iDir] << "_pt_"
0100             << iPt << "_sqrt";
0101         std::stringstream htitle;
0102         htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0103         h_AN_pt_shuffle[iP][iPt][iDir] = new TH1F(
0104           hname.str().c_str(),
0105           htitle.str().c_str(),
0106           100, -4, 4);
0107         // Beam- and config- specific values
0108         for (int iB = 0; iB < ASYM_CONSTANTS::nBeams; iB++) {
0109           for (int iC = 0; iC < ASYM_CONSTANTS::nConfigs; iC++) {
0110             hname.str("");
0111             hname << "h_AN_pt_shuffle_"
0112                   << ASYM_CONSTANTS::beams[iB] << "_"
0113                   << ASYM_CONSTANTS::particle[iP] << "_"
0114                   << ASYM_CONSTANTS::directions[iDir] << "_pt_"
0115                   << iPt << "_"
0116                   << ASYM_CONSTANTS::configuration[iC] << "_sqrt";
0117             htitle.str("");
0118             htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0119             h_AN_pt_shuffle_details[iP][iPt][iDir][iB][iC] = new TH1F(
0120               hname.str().c_str(),
0121               htitle.str().c_str(),
0122               100, -4, 4);
0123           }
0124         }
0125       }
0126     }
0127     for (int iEta = 0; iEta < ASYM_CONSTANTS::nEtaBins; iEta++) {
0128       std::stringstream hname;
0129       hname << "h_AN_eta_shuffle_"
0130             << ASYM_CONSTANTS::particle[iP] << "_eta_"
0131             << iEta << "_sqrt";
0132       std::stringstream htitle;
0133       htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0134       h_AN_eta_shuffle[iP][iEta] = new TH1F(
0135         hname.str().c_str(),
0136         htitle.str().c_str(),
0137         100, -4, 4);
0138       // Beam- and config- specific values
0139       for (int iB = 0; iB < ASYM_CONSTANTS::nBeams; iB++) {
0140         for (int iC = 0; iC < ASYM_CONSTANTS::nConfigs; iC++) {
0141           hname.str("");
0142           hname << "h_AN_eta_shuffle_"
0143                 << ASYM_CONSTANTS::beams[iB] << "_"
0144                 << ASYM_CONSTANTS::particle[iP] << "_eta_"
0145                 << iEta << "_"
0146                 << ASYM_CONSTANTS::configuration[iC] << "_sqrt";
0147           htitle.str("");
0148           htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0149           h_AN_eta_shuffle_details[iP][iEta][iB][iC] = new TH1F(
0150             hname.str().c_str(),
0151             htitle.str().c_str(),
0152             100, -4, 4);
0153         }
0154       }
0155     }
0156     for (int iXf = 0; iXf < ASYM_CONSTANTS::nXfBins; iXf++) {
0157       std::stringstream hname;
0158       hname << "h_AN_xf_shuffle_"
0159             << ASYM_CONSTANTS::particle[iP] << "_xf_"
0160             << iXf << "_sqrt";
0161       std::stringstream htitle;
0162       htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0163       h_AN_xf_shuffle[iP][iXf] = new TH1F(
0164         hname.str().c_str(),
0165         htitle.str().c_str(),
0166         100, -4, 4);
0167       // Beam- and config- specific values
0168       for (int iB = 0; iB < ASYM_CONSTANTS::nBeams; iB++) {
0169         for (int iC = 0; iC < ASYM_CONSTANTS::nConfigs; iC++) {
0170           hname.str("");
0171           hname << "h_AN_xf_shuffle_"
0172                 << ASYM_CONSTANTS::beams[iB] << "_"
0173                 << ASYM_CONSTANTS::particle[iP] << "_xf_"
0174                 << iXf << "_"
0175                 << ASYM_CONSTANTS::configuration[iC] << "_sqrt";
0176           htitle.str("");
0177           htitle << ";A_{N}/#delta A_{N}^{stat}; Counts";
0178           h_AN_xf_shuffle_details[iP][iXf][iB][iC] = new TH1F(
0179             hname.str().c_str(),
0180             htitle.str().c_str(),
0181             100, -4, 4);
0182         }
0183       }
0184     }
0185   }
0186 }
0187
0188 void fill_shuffle_detailed_histos(
0189        AsymmetryCalc::ShuffleBunches *shuffle_bunches_0mrad_mbd,
0190        AsymmetryCalc::ShuffleBunches *shuffle_bunches_0mrad_photon,
0191        AsymmetryCalc::ShuffleBunches *shuffle_bunches_15mrad_mbd,
0192        AsymmetryCalc::ShuffleBunches *shuffle_bunches_15mrad_photon,
0193        TH1F* h_AN_pt_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0194        TH1F* h_AN_pt_nodir_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0195        TH1F* h_AN_eta_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0196        TH1F* h_AN_xf_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs],
0197        const int iIter)
0198 {
0199   // Total of 4 configurations
0200   const int nConfigs = ASYM_CONSTANTS::nConfigs;
0201   AsymmetryCalc::ShuffleBunches *shuffle_bunches[nConfigs] = {
0202     shuffle_bunches_0mrad_mbd,
0203     shuffle_bunches_0mrad_photon,
0204     shuffle_bunches_15mrad_mbd,
0205     shuffle_bunches_15mrad_photon
0206   };
0207
0208   // Load corrected asymmetries from every configuration
0209   std::vector<const AsymmetryCalc::pt_corr_array*> corrected_mean_pt(nConfigs);
0210   std::vector<const AsymmetryCalc::pt_nodir_corr_array*> corrected_mean_pt_nodir(nConfigs);
0211   std::vector<const AsymmetryCalc::eta_corr_array*> corrected_mean_eta(nConfigs);
0212   std::vector<const AsymmetryCalc::xf_corr_array*> corrected_mean_xf(nConfigs);
0213   std::vector<const AsymmetryCalc::pt_corr_array*> corrected_unc_pt(nConfigs);
0214   std::vector<const AsymmetryCalc::pt_nodir_corr_array*> corrected_unc_pt_nodir(nConfigs);
0215   std::vector<const AsymmetryCalc::eta_corr_array*> corrected_unc_eta(nConfigs);
0216   std::vector<const AsymmetryCalc::xf_corr_array*> corrected_unc_xf(nConfigs);
0217
0218   for (int iConfig = 0; iConfig < nConfigs; iConfig++) {
0219     corrected_mean_pt[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_pt()[iIter];
0220     corrected_mean_pt_nodir[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_pt_nodir()[iIter];
0221     corrected_mean_eta[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_eta()[iIter];
0222     corrected_mean_xf[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_xf()[iIter];
0223     corrected_unc_pt[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_pt()[iIter];
0224     corrected_unc_pt_nodir[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_pt_nodir()[iIter];
0225     corrected_unc_eta[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_eta()[iIter];
0226     corrected_unc_xf[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_xf()[iIter];
0227   }
0228
0229   for (int iB = 0; iB < ASYM_CONSTANTS::nBeams; iB++) {
0230     for (int iC = 0; iC < ASYM_CONSTANTS::nConfigs; iC++) {
0231       for (int iP = 0; iP < ASYM_CONSTANTS::nParticles; iP++) {
0232         for (int iPt = 0; iPt < ASYM_CONSTANTS::nPtBins; iPt++) {
0233           {
0234             // Beam- and config- specific values
0235             h_AN_pt_nodir_shuffle_details[iP][iPt][iB][iC]->Fill(
0236               (*corrected_mean_pt_nodir[iC])(iB, iP, iPt) /
0237               (*corrected_unc_pt_nodir[iC])(iB, iP, iPt));
0238           }
0239           for (int iDir = 0; iDir < ASYM_CONSTANTS::nDirections; iDir++) {
0240             // Beam- and config- specific values
0241             h_AN_pt_shuffle_details[iP][iPt][iDir][iB][iC]->Fill(
0242               (*corrected_mean_pt[iC])(iB, iP, iPt, iDir) /
0243               (*corrected_unc_pt[iC])(iB, iP, iPt, iDir));
0244           }
0245         }
0246         for (int iEta = 0; iEta < ASYM_CONSTANTS::nEtaBins; iEta++) {
0247           h_AN_eta_shuffle_details[iP][iEta][iB][iC]->Fill(
0248             (*corrected_mean_eta[iC])(iB, iP, iEta) /
0249             (*corrected_unc_eta[iC])(iB, iP, iEta));
0250
0251         }
0252         for (int iXf = 0; iXf < ASYM_CONSTANTS::nXfBins; iXf++) {
0253           h_AN_xf_shuffle_details[iP][iXf][iB][iC]->Fill(
0254             (*corrected_mean_xf[iC])(iB, iP, iXf) /
0255             (*corrected_unc_xf[iC])(iB, iP, iXf));
0256         }
0257       }
0258     }
0259   }
0260 }
0261
0262 void fill_shuffle_histos(AsymmetryCalc::pt_final_array& merged_mean_pt,
0263                          AsymmetryCalc::pt_nodir_final_array& merged_mean_pt_nodir,
0264                          AsymmetryCalc::eta_final_array& merged_mean_eta,
0265                          AsymmetryCalc::xf_final_array& merged_mean_xf,
0266                          AsymmetryCalc::pt_final_array& merged_unc_pt,
0267                          AsymmetryCalc::pt_nodir_final_array& merged_unc_pt_nodir,
0268                          AsymmetryCalc::eta_final_array& merged_unc_eta,
0269                          AsymmetryCalc::xf_final_array& merged_unc_xf,
0270                          TH1F* h_AN_pt_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections],
0271                          TH1F* h_AN_pt_nodir_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins],
0272                          TH1F* h_AN_eta_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins],
0273                          TH1F* h_AN_xf_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins])
0274 {
0275   for (int iP = 0; iP < ASYM_CONSTANTS::nParticles; iP++) {
0276     for (int iPt = 0; iPt < ASYM_CONSTANTS::nPtBins; iPt++) {
0277       {
0278         h_AN_pt_nodir_shuffle[iP][iPt]->Fill(merged_mean_pt_nodir(iP, iPt) / merged_unc_pt_nodir(iP, iPt));
0279       }
0280       for (int iDir = 0; iDir < ASYM_CONSTANTS::nDirections; iDir++) {
0281         h_AN_pt_shuffle[iP][iPt][iDir]->Fill(merged_mean_pt(iP, iPt, iDir) / merged_unc_pt(iP, iPt, iDir));
0282       }
0283     }
0284     for (int iEta = 0; iEta < ASYM_CONSTANTS::nEtaBins; iEta++) {
0285       h_AN_eta_shuffle[iP][iEta]->Fill(merged_mean_eta(iP, iEta) / merged_unc_eta(iP, iEta));
0286     }
0287     for (int iXf = 0; iXf < ASYM_CONSTANTS::nXfBins; iXf++) {
0288       h_AN_xf_shuffle[iP][iXf]->Fill(merged_mean_xf(iP, iXf) / merged_unc_xf(iP, iXf));
0289     }
0290   }
0291 }
0292
0293 void save_shuffle_histos(TFile *output_file_shuffle)
0294 {
0295   output_file_shuffle->cd();
0296   output_file_shuffle->Write();
0297   output_file_shuffle->Close();
0298   delete output_file_shuffle;
0299   output_file_shuffle = nullptr;
0300 }
0301
0302 static bool valid_measurement(double val, double unc)
0303 {
0304   return std::isfinite(val) && std::isfinite(unc) && unc > 0.0 && val != 0;
0305 }
0306
0307 void compute_average_2(double& res_val, double& res_unc,
0308                        double in_val_1, double in_unc_1,
0309                        double in_val_2, double in_unc_2)
0310 {
0311     const bool ok1 = valid_measurement(in_val_1, in_unc_1);
0312     const bool ok2 = valid_measurement(in_val_2, in_unc_2);
0313
0314     if (!ok1 && !ok2) {
0315         res_val = 0.0;
0316         res_unc = 0.0;
0317         return;
0318     }
0319     if (!ok1) {
0320         res_val = in_val_2;
0321         res_unc = in_unc_2;
0322         return;
0323     }
0324     if (!ok2) {
0325         res_val = in_val_1;
0326         res_unc = in_unc_1;
0327         return;
0328     }
0329
0330     const double w1 = 1.0 / (in_unc_1 * in_unc_1);
0331     const double w2 = 1.0 / (in_unc_2 * in_unc_2);
0332
0333     res_val = (w1 * in_val_1 + w2 * in_val_2) / (w1 + w2);
0334     res_unc = 1.0 / std::sqrt(w1 + w2);
0335 }
0336
0337 void save_final_histograms(const std::string &name,
0338                            AsymmetryCalc::pt_final_array& merged_mean_pt,
0339                            AsymmetryCalc::pt_nodir_final_array& merged_mean_pt_nodir,
0340                            AsymmetryCalc::eta_final_array& merged_mean_eta,
0341                            AsymmetryCalc::xf_final_array& merged_mean_xf,
0342                            AsymmetryCalc::pt_final_array& merged_unc_pt,
0343                            AsymmetryCalc::pt_nodir_final_array& merged_unc_pt_nodir,
0344                            AsymmetryCalc::eta_final_array& merged_unc_eta,
0345                            AsymmetryCalc::xf_final_array& merged_unc_xf)
0346 {
0347   TFile *outputfile = new TFile(name.c_str(), "RECREATE");
0348   outputfile->cd();
0349   for (int iP = 0; iP < ASYM_CONSTANTS::nParticles; iP++) {
0350     for (int iDir = 0; iDir < ASYM_CONSTANTS::nDirections; iDir++) {
0351       TGraphErrors *graph_AN = new TGraphErrors();
0352       std::stringstream graph_name;
0353       graph_name << "graph_AN_pt_final_" << ASYM_CONSTANTS::particle[iP] << "_"
0354                  << ASYM_CONSTANTS::directions[iDir] << "_sqrt";
0355       std::stringstream graph_title;
0356       graph_title << "; iPt; A_N";
0357       graph_AN->SetName(graph_name.str().c_str());
0358       graph_AN->SetTitle(graph_title.str().c_str());
0359       for (int iPt = 0; iPt < ASYM_CONSTANTS::nPtBins; iPt++) {
0360         graph_AN->SetPoint(iPt, (double)(iPt+1), merged_mean_pt(iP, iPt, iDir));
0361         graph_AN->SetPointError(iPt, 0, merged_unc_pt(iP, iPt, iDir));
0362       }
0363       graph_AN->Write();
0364     }
0365     {
0366       TGraphErrors *graph_AN = new TGraphErrors();
0367       std::stringstream graph_name;
0368       graph_name << "graph_AN_eta_final_" << ASYM_CONSTANTS::particle[iP] << "_sqrt";
0369       std::stringstream graph_title;
0370       graph_title << "; iEta; A_N";
0371       graph_AN->SetName(graph_name.str().c_str());
0372       graph_AN->SetTitle(graph_title.str().c_str());
0373       for (int iEta = 0; iEta < ASYM_CONSTANTS::nEtaBins; iEta++) {
0374         graph_AN->SetPoint(iEta, (double)(iEta+1), merged_mean_eta(iP, iEta));
0375         graph_AN->SetPointError(iEta, 0, merged_unc_eta(iP, iEta));
0376       }
0377       graph_AN->Write();
0378     }
0379     {
0380       TGraphErrors *graph_AN = new TGraphErrors();
0381       std::stringstream graph_name;
0382       graph_name << "graph_AN_xf_final_" << ASYM_CONSTANTS::particle[iP] << "_sqrt";
0383       std::stringstream graph_title;
0384       graph_title << "; iXf; A_N";
0385       graph_AN->SetName(graph_name.str().c_str());
0386       graph_AN->SetTitle(graph_title.str().c_str());
0387       for (int iXf = 0; iXf < ASYM_CONSTANTS::nXfBins; iXf++) {
0388         graph_AN->SetPoint(iXf, (double)(iXf+1), merged_mean_xf(iP, iXf));
0389         graph_AN->SetPointError(iXf, 0, merged_unc_xf(iP, iXf));
0390       }
0391       graph_AN->Write();
0392     }
0393   }
0394   outputfile->Close();
0395   delete outputfile;
0396 }
0397
0398 void merge_configurations(AsymmetryCalc::ShuffleBunches *shuffle_bunches_0mrad_mbd,
0399                           AsymmetryCalc::ShuffleBunches *shuffle_bunches_0mrad_photon,
0400                           AsymmetryCalc::ShuffleBunches *shuffle_bunches_15mrad_mbd,
0401                           AsymmetryCalc::ShuffleBunches *shuffle_bunches_15mrad_photon,
0402                           AsymmetryCalc::pt_final_array& merged_mean_pt,
0403                           AsymmetryCalc::pt_nodir_final_array& merged_mean_pt_nodir,
0404                           AsymmetryCalc::eta_final_array& merged_mean_eta,
0405                           AsymmetryCalc::xf_final_array& merged_mean_xf,
0406                           AsymmetryCalc::pt_final_array& merged_unc_pt,
0407                           AsymmetryCalc::pt_nodir_final_array& merged_unc_pt_nodir,
0408                           AsymmetryCalc::eta_final_array& merged_unc_eta,
0409                           AsymmetryCalc::xf_final_array& merged_unc_xf,
0410                           const int iIter)
0411 {
0412   // Total of 4 configurations
0413   const int nConfigs = ASYM_CONSTANTS::nConfigs;
0414   AsymmetryCalc::ShuffleBunches *shuffle_bunches[nConfigs] = {
0415     shuffle_bunches_0mrad_mbd,
0416     shuffle_bunches_0mrad_photon,
0417     shuffle_bunches_15mrad_mbd,
0418     shuffle_bunches_15mrad_photon
0419   };
0420
0421   // Load corrected asymmetries from every configuration
0422   std::vector<const AsymmetryCalc::pt_corr_array*> corrected_mean_pt(nConfigs);
0423   std::vector<const AsymmetryCalc::pt_nodir_corr_array*> corrected_mean_pt_nodir(nConfigs);
0424   std::vector<const AsymmetryCalc::eta_corr_array*> corrected_mean_eta(nConfigs);
0425   std::vector<const AsymmetryCalc::xf_corr_array*> corrected_mean_xf(nConfigs);
0426   std::vector<const AsymmetryCalc::pt_corr_array*> corrected_unc_pt(nConfigs);
0427   std::vector<const AsymmetryCalc::pt_nodir_corr_array*> corrected_unc_pt_nodir(nConfigs);
0428   std::vector<const AsymmetryCalc::eta_corr_array*> corrected_unc_eta(nConfigs);
0429   std::vector<const AsymmetryCalc::xf_corr_array*> corrected_unc_xf(nConfigs);
0430
0431   for (int iConfig = 0; iConfig < nConfigs; iConfig++) {
0432     corrected_mean_pt[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_pt()[iIter];
0433     corrected_mean_pt_nodir[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_pt_nodir()[iIter];
0434     corrected_mean_eta[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_eta()[iIter];
0435     corrected_mean_xf[iConfig] = &shuffle_bunches[iConfig]->get_corrected_mean_xf()[iIter];
0436     corrected_unc_pt[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_pt()[iIter];
0437     corrected_unc_pt_nodir[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_pt_nodir()[iIter];
0438     corrected_unc_eta[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_eta()[iIter];
0439     corrected_unc_xf[iConfig] = &shuffle_bunches[iConfig]->get_corrected_unc_xf()[iIter];
0440   }
0441
0442   // Average the 4 configurations in each kinematic bin
0443   for (int iP = 0; iP < ASYM_CONSTANTS::nParticles; iP++) {
0444     for (int iPt = 0; iPt < ASYM_CONSTANTS::nPtBins; iPt++) {
0445       {
0446         // First average over trigger
0447         double mean_trigger[4] = {0};
0448         double unc_trigger[4] = {0};
0449         compute_average_2(
0450              mean_trigger[0], unc_trigger[0],
0451              (*corrected_mean_pt_nodir[0])(0, iP, iPt),
0452              (*corrected_unc_pt_nodir[0])(0, iP, iPt),
0453              (*corrected_mean_pt_nodir[1])(0, iP, iPt),
0454              (*corrected_unc_pt_nodir[1])(0, iP, iPt)); // Yellow 0 mrad
0455
0456         compute_average_2(
0457              mean_trigger[1], unc_trigger[1],
0458              (*corrected_mean_pt_nodir[2])(0, iP, iPt),
0459              (*corrected_unc_pt_nodir[2])(0, iP, iPt),
0460              (*corrected_mean_pt_nodir[3])(0, iP, iPt),
0461              (*corrected_unc_pt_nodir[3])(0, iP, iPt)); // Yellow 1.5 mrad
0462
0463         compute_average_2(
0464              mean_trigger[2], unc_trigger[2],
0465              (*corrected_mean_pt_nodir[0])(1, iP, iPt),
0466              (*corrected_unc_pt_nodir[0])(1, iP, iPt),
0467              (*corrected_mean_pt_nodir[1])(1, iP, iPt),
0468              (*corrected_unc_pt_nodir[1])(1, iP, iPt)); // Blue 0 mrad
0469
0470         compute_average_2(
0471              mean_trigger[3], unc_trigger[3],
0472              (*corrected_mean_pt_nodir[2])(1, iP, iPt),
0473              (*corrected_unc_pt_nodir[2])(1, iP, iPt),
0474              (*corrected_mean_pt_nodir[3])(1, iP, iPt),
0475              (*corrected_unc_pt_nodir[3])(1, iP, iPt)); // Blue 1.5 mrad
0476
0477         // Then average over crossing angle
0478         double mean_trigger_angle[2] = {0};
0479         double unc_trigger_angle[2] = {0};
0480         compute_average_2(
0481              mean_trigger_angle[0], unc_trigger_angle[0],
0482              mean_trigger[0], unc_trigger[0],
0483              mean_trigger[1], unc_trigger[1]); // Yellow
0484         compute_average_2(
0485              mean_trigger_angle[1], unc_trigger_angle[1],
0486              mean_trigger[2], unc_trigger[2],
0487              mean_trigger[3], unc_trigger[3]); // Blue
0488
0489         // And finally average over beam
0490         double mean_final = 0;
0491         double unc_final = 0;
0492         compute_average_2(
0493              mean_final, unc_final,
0494              mean_trigger_angle[0], unc_trigger_angle[0],
0495              mean_trigger_angle[1], unc_trigger_angle[1]);
0496         merged_mean_pt_nodir(iP, iPt) = mean_final;
0497         merged_unc_pt_nodir(iP, iPt) = unc_final;
0498       }
0499       for (int iDir = 0; iDir < ASYM_CONSTANTS::nDirections; iDir++) {
0500
0501         // First average over trigger
0502         double mean_trigger[4] = {0};
0503         double unc_trigger[4] = {0};
0504         compute_average_2(
0505              mean_trigger[0], unc_trigger[0],
0506              (*corrected_mean_pt[0])(0, iP, iPt, iDir),
0507              (*corrected_unc_pt[0])(0, iP, iPt, iDir),
0508              (*corrected_mean_pt[1])(0, iP, iPt, iDir),
0509              (*corrected_unc_pt[1])(0, iP, iPt, iDir)); // Yellow 0 mrad
0510
0511         compute_average_2(
0512              mean_trigger[1], unc_trigger[1],
0513              (*corrected_mean_pt[2])(0, iP, iPt, iDir),
0514              (*corrected_unc_pt[2])(0, iP, iPt, iDir),
0515              (*corrected_mean_pt[3])(0, iP, iPt, iDir),
0516              (*corrected_unc_pt[3])(0, iP, iPt, iDir)); // Yellow 1.5 mrad
0517
0518         compute_average_2(
0519              mean_trigger[2], unc_trigger[2],
0520              (*corrected_mean_pt[0])(1, iP, iPt, iDir),
0521              (*corrected_unc_pt[0])(1, iP, iPt, iDir),
0522              (*corrected_mean_pt[1])(1, iP, iPt, iDir),
0523              (*corrected_unc_pt[1])(1, iP, iPt, iDir)); // Blue 0 mrad
0524
0525         compute_average_2(
0526              mean_trigger[3], unc_trigger[3],
0527              (*corrected_mean_pt[2])(1, iP, iPt, iDir),
0528              (*corrected_unc_pt[2])(1, iP, iPt, iDir),
0529              (*corrected_mean_pt[3])(1, iP, iPt, iDir),
0530              (*corrected_unc_pt[3])(1, iP, iPt, iDir)); // Blue 1.5 mrad
0531
0532         // Then average over crossing angle
0533         double mean_trigger_angle[2] = {0};
0534         double unc_trigger_angle[2] = {0};
0535         compute_average_2(
0536              mean_trigger_angle[0], unc_trigger_angle[0],
0537              mean_trigger[0], unc_trigger[0],
0538              mean_trigger[1], unc_trigger[1]); // Yellow
0539         compute_average_2(
0540              mean_trigger_angle[1], unc_trigger_angle[1],
0541              mean_trigger[2], unc_trigger[2],
0542              mean_trigger[3], unc_trigger[3]); // Blue
0543
0544         // And finally average over beam
0545         double mean_final = 0;
0546         double unc_final = 0;
0547         compute_average_2(
0548              mean_final, unc_final,
0549              mean_trigger_angle[0], unc_trigger_angle[0],
0550              mean_trigger_angle[1], unc_trigger_angle[1]);
0551         merged_mean_pt(iP, iPt, iDir) = mean_final;
0552         merged_unc_pt(iP, iPt, iDir) = unc_final;
0553       }
0554     }
0555     for (int iEta = 0; iEta < ASYM_CONSTANTS::nEtaBins; iEta++) {
0556       // First average over trigger
0557       double mean_trigger[4] = {0};
0558       double unc_trigger[4] = {0};
0559
0560       compute_average_2(
0561            mean_trigger[0], unc_trigger[0],
0562            (*corrected_mean_eta[0])(0, iP, iEta),
0563            (*corrected_unc_eta[0])(0, iP, iEta),
0564            (*corrected_mean_eta[1])(0, iP, iEta),
0565            (*corrected_unc_eta[1])(0, iP, iEta)); // Yellow 0 mrad
0566       compute_average_2(
0567            mean_trigger[1], unc_trigger[1],
0568            (*corrected_mean_eta[2])(0, iP, iEta),
0569            (*corrected_unc_eta[2])(0, iP, iEta),
0570            (*corrected_mean_eta[3])(0, iP, iEta),
0571            (*corrected_unc_eta[3])(0, iP, iEta)); // Yellow 1.5 mrad
0572       compute_average_2(
0573            mean_trigger[2], unc_trigger[2],
0574            (*corrected_mean_eta[0])(1, iP, iEta),
0575            (*corrected_unc_eta[0])(1, iP, iEta),
0576            (*corrected_mean_eta[1])(1, iP, iEta),
0577            (*corrected_unc_eta[1])(1, iP, iEta)); // Blue 0 mrad
0578       compute_average_2(
0579            mean_trigger[3], unc_trigger[3],
0580            (*corrected_mean_eta[2])(1, iP, iEta),
0581            (*corrected_unc_eta[2])(1, iP, iEta),
0582            (*corrected_mean_eta[3])(1, iP, iEta),
0583            (*corrected_unc_eta[3])(1, iP, iEta)); // Blue 1.5 mrad
0584
0585       // Then average over crossing angle
0586       double mean_trigger_angle[2] = {0};
0587       double unc_trigger_angle[2] = {0};
0588       compute_average_2(
0589            mean_trigger_angle[0], unc_trigger_angle[0],
0590            mean_trigger[0], unc_trigger[0],
0591            mean_trigger[1], unc_trigger[1]); // Yellow
0592       compute_average_2(
0593            mean_trigger_angle[1], unc_trigger_angle[1],
0594            mean_trigger[2], unc_trigger[2],
0595            mean_trigger[3], unc_trigger[3]); // Blue
0596
0597       // And finally average over beam
0598       double mean_final = 0;
0599       double unc_final = 0;
0600       compute_average_2(
0601            mean_final, unc_final,
0602            mean_trigger_angle[0], unc_trigger_angle[0],
0603            mean_trigger_angle[1], unc_trigger_angle[1]);
0604       //std::cout << "final (" << mean_final << ", " << unc_final << ")\n\n" << std::endl;
0605
0606       merged_mean_eta(iP, iEta) = mean_final;
0607       merged_unc_eta(iP, iEta) = unc_final;
0608     }
0609     for (int iXf = 0; iXf < ASYM_CONSTANTS::nXfBins; iXf++) {
0610
0611       // First average over trigger
0612       double mean_trigger[4] = {0};
0613       double unc_trigger[4] = {0};
0614
0615       compute_average_2(
0616            mean_trigger[0], unc_trigger[0],
0617            (*corrected_mean_xf[0])(0, iP, iXf),
0618            (*corrected_unc_xf[0])(0, iP, iXf),
0619            (*corrected_mean_xf[1])(0, iP, iXf),
0620            (*corrected_unc_xf[1])(0, iP, iXf)); // Yellow 0 mrad
0621       compute_average_2(
0622            mean_trigger[1], unc_trigger[1],
0623            (*corrected_mean_xf[2])(0, iP, iXf),
0624            (*corrected_unc_xf[2])(0, iP, iXf),
0625            (*corrected_mean_xf[3])(0, iP, iXf),
0626            (*corrected_unc_xf[3])(0, iP, iXf)); // Yellow 1.5 mrad
0627       compute_average_2(
0628            mean_trigger[2], unc_trigger[2],
0629            (*corrected_mean_xf[0])(1, iP, iXf),
0630            (*corrected_unc_xf[0])(1, iP, iXf),
0631            (*corrected_mean_xf[1])(1, iP, iXf),
0632            (*corrected_unc_xf[1])(1, iP, iXf)); // Blue 0 mrad
0633       compute_average_2(
0634            mean_trigger[3], unc_trigger[3],
0635            (*corrected_mean_xf[2])(1, iP, iXf),
0636            (*corrected_unc_xf[2])(1, iP, iXf),
0637            (*corrected_mean_xf[3])(1, iP, iXf),
0638            (*corrected_unc_xf[3])(1, iP, iXf)); // Blue 1.5 mrad
0639
0640       // Then average over crossing angle
0641       double mean_trigger_angle[2] = {0};
0642       double unc_trigger_angle[2] = {0};
0643       compute_average_2(
0644            mean_trigger_angle[0], unc_trigger_angle[0],
0645            mean_trigger[0], unc_trigger[0],
0646            mean_trigger[1], unc_trigger[1]); // Yellow
0647       compute_average_2(
0648            mean_trigger_angle[1], unc_trigger_angle[1],
0649            mean_trigger[2], unc_trigger[2],
0650            mean_trigger[3], unc_trigger[3]); // Blue
0651
0652       // And finally average over beam
0653       double mean_final = 0;
0654       double unc_final = 0;
0655       compute_average_2(
0656            mean_final, unc_final,
0657            mean_trigger_angle[0], unc_trigger_angle[0],
0658            mean_trigger_angle[1], unc_trigger_angle[1]);
0659
0660       merged_mean_xf(iP, iXf) = mean_final;
0661       merged_unc_xf(iP, iXf) = unc_final;
0662     }
0663   }
0664 }
0665
0666 int main(int argc, char *argv[])
0667 {
0668   if (argc < 3) {
0669     std::cout << "Usage: ./BunchShuffling <iterMin> <iterMax> [<config file> (default: config.txt)]" << std::endl;
0670     return 1;
0671   }
0672
0673   int shuffleIterMin = std::stoi(argv[1]);
0674   int shuffleIterMax = std::stoi(argv[2]);
0675
0676   std::string config_filename = "config.txt";
0677   if (argc > 3)
0678   {
0679     config_filename = argv[3];
0680   }
0681
0682   const int shuffleMaxStorage = 10; // To avoid having too much resident memory
0683
0684   bool perform_shuffle = true; // Set to false to deactivate the shuffle (debug only)
0685   bool store_raw = false;
0686   bool store_final = false;
0687
0688   // Default path
0689   std::unordered_map<std::string, std::string> config = read_config(config_filename);
0690
0691   std::string path_to_bunch_dependent_yields_mbd = config["path_to_bunch_dependent_yields_mbd"];
0692   std::string path_to_bunch_dependent_yields_photon = config["path_to_bunch_dependent_yields_photon"];
0693   std::string path_to_fill_list_0mrad = config["path_to_fill_list_0mrad"];
0694   std::string path_to_fill_list_15mrad = config["path_to_fill_list_15mrad"];
0695   std::string path_to_spin_pattern_0mrad = config["path_to_spin_pattern_0mrad"];
0696   std::string path_to_spin_pattern_15mrad = config["path_to_spin_pattern_15mrad"];
0697   std::string path_to_csv_mean = config["path_to_csv_mean"];
0698   std::string path_to_csv_ratios = config["path_to_csv_ratios"];
0699   std::string suffix = config["suffix"]; // This parameter is allowed to be empty
0700   std::string stored_iter_asym_histos = config["stored_iter_asym_histos"];
0701   std::string stored_final_asym_histos = config["stored_final_asym_histos"];
0702
0703   if (path_to_bunch_dependent_yields_mbd.empty())
0704   {
0705     path_to_bunch_dependent_yields_mbd = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/analysis_complete_ana509_01312026_MBD_pt_l4_0mrad/fast_analysis_";
0706   }
0707   if (path_to_bunch_dependent_yields_photon.empty())
0708   {
0709     path_to_bunch_dependent_yields_photon = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/analysis_complete_ana509_01312026_efficiency_30_pt_g3_0mrad/fast_analysis_";
0710   }
0711   if (path_to_fill_list_0mrad.empty())
0712   {
0713     path_to_fill_list_0mrad = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/run_to_fill_0mrad.txt";
0714   }
0715   if (path_to_fill_list_15mrad.empty())
0716   {
0717     path_to_fill_list_15mrad = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/run_to_fill_15mrad.txt";
0718   }
0719   if (path_to_spin_pattern_0mrad.empty())
0720   {
0721     path_to_spin_pattern_0mrad = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/spin_pattern_fill_0mrad.root";
0722   }
0723   if (path_to_spin_pattern_15mrad.empty())
0724   {
0725     path_to_spin_pattern_15mrad = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/spin_pattern_fill_15mrad.root";
0726   }
0727   if (path_to_csv_mean.empty())
0728   {
0729     path_to_csv_mean = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/csv_inputs";
0730   }
0731   if (path_to_csv_ratios.empty())
0732   {
0733     path_to_csv_ratios = "/sphenix/u/virgilemahaut/work/analysis/AnNeutralMeson/macros/comparison_micro_DST_analysis/nano_analysis/csv_ratios";
0734   }
0735
0736   if (!stored_iter_asym_histos.empty()) {
0737     store_raw = true;
0738   }
0739   if (!stored_final_asym_histos.empty()) {
0740     store_final = true;
0741   }
0742
0743   // All pT/eta/xF means and ratios
0744   std::string input_csv_mbd_0mrad_pt_pi0_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0745   std::string input_csv_mbd_0mrad_pt_eta_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0746   std::string input_csv_mbd_0mrad_eta_pi0_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0747   std::string input_csv_mbd_0mrad_eta_eta_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0748   std::string input_csv_mbd_0mrad_xf_pi0_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0749   std::string input_csv_mbd_0mrad_xf_eta_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0750   std::string input_csv_mbd_0mrad_pt_ratios_pi0_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0751   std::string input_csv_mbd_0mrad_pt_ratios_eta_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0752   std::string input_csv_mbd_0mrad_eta_ratios_pi0_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0753   std::string input_csv_mbd_0mrad_eta_ratios_eta_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0754   std::string input_csv_mbd_0mrad_xf_ratios_pi0_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_MBD_pt_l3_0mrad.csv";
0755   std::string input_csv_mbd_0mrad_xf_ratios_eta_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_MBD_pt_l4_0mrad.csv";
0756
0757   std::string input_csv_photon_0mrad_pt_pi0_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0758   std::string input_csv_photon_0mrad_pt_eta_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0759   std::string input_csv_photon_0mrad_eta_pi0_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0760   std::string input_csv_photon_0mrad_eta_eta_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0761   std::string input_csv_photon_0mrad_xf_pi0_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0762   std::string input_csv_photon_0mrad_xf_eta_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0763   std::string input_csv_photon_0mrad_pt_ratios_pi0_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0764   std::string input_csv_photon_0mrad_pt_ratios_eta_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0765   std::string input_csv_photon_0mrad_eta_ratios_pi0_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0766   std::string input_csv_photon_0mrad_eta_ratios_eta_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0767   std::string input_csv_photon_0mrad_xf_ratios_pi0_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_efficiency_30_pt_g3_0mrad.csv";
0768   std::string input_csv_photon_0mrad_xf_ratios_eta_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_efficiency_30_pt_g4_0mrad.csv";
0769
0770   std::string input_csv_mbd_15mrad_pt_pi0_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0771   std::string input_csv_mbd_15mrad_pt_eta_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0772   std::string input_csv_mbd_15mrad_eta_pi0_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0773   std::string input_csv_mbd_15mrad_eta_eta_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0774   std::string input_csv_mbd_15mrad_xf_pi0_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0775   std::string input_csv_mbd_15mrad_xf_eta_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0776   std::string input_csv_mbd_15mrad_pt_ratios_pi0_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0777   std::string input_csv_mbd_15mrad_pt_ratios_eta_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0778   std::string input_csv_mbd_15mrad_eta_ratios_pi0_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0779   std::string input_csv_mbd_15mrad_eta_ratios_eta_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0780   std::string input_csv_mbd_15mrad_xf_ratios_pi0_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_MBD_pt_l3_15mrad.csv";
0781   std::string input_csv_mbd_15mrad_xf_ratios_eta_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_MBD_pt_l4_15mrad.csv";
0782
0783   std::string input_csv_photon_15mrad_pt_pi0_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0784   std::string input_csv_photon_15mrad_pt_eta_name = path_to_csv_mean + "/input_pt_mean_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0785   std::string input_csv_photon_15mrad_eta_pi0_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0786   std::string input_csv_photon_15mrad_eta_eta_name = path_to_csv_mean + "/input_eta_mean_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0787   std::string input_csv_photon_15mrad_xf_pi0_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0788   std::string input_csv_photon_15mrad_xf_eta_name = path_to_csv_mean + "/input_xf_mean_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0789   std::string input_csv_photon_15mrad_pt_ratios_pi0_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0790   std::string input_csv_photon_15mrad_pt_ratios_eta_name = path_to_csv_ratios + "/pt_ratios_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0791   std::string input_csv_photon_15mrad_eta_ratios_pi0_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0792   std::string input_csv_photon_15mrad_eta_ratios_eta_name = path_to_csv_ratios + "/eta_ratios_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0793   std::string input_csv_photon_15mrad_xf_ratios_pi0_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_efficiency_30_pt_g3_15mrad.csv";
0794   std::string input_csv_photon_15mrad_xf_ratios_eta_name = path_to_csv_ratios + "/xf_ratios_01312026" + suffix + "_efficiency_30_pt_g4_15mrad.csv";
0795
0796   int nIterations = (shuffleIterMax - shuffleIterMin + 1);
0797   int nStore = nIterations / shuffleMaxStorage;
0798   if (nIterations % shuffleMaxStorage != 0) {
0799     nStore++;
0800   }
0801
0802   std::stringstream shuffle_file_name;
0803   shuffle_file_name << "shuffle_" << shuffleIterMin << "_" << shuffleIterMax << ".root";
0804   TFile *output_file_shuffle = nullptr;
0805   TH1F *h_AN_pt_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections] = {nullptr};
0806   TH1F *h_AN_pt_nodir_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins] = {nullptr};
0807   TH1F *h_AN_eta_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins] = {nullptr};
0808   TH1F *h_AN_xf_shuffle[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins] = {nullptr};
0809   TH1F *h_AN_pt_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs] {nullptr};
0810   TH1F *h_AN_pt_nodir_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs] {nullptr};
0811   TH1F *h_AN_eta_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs] {nullptr};
0812   TH1F *h_AN_xf_shuffle_details[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nConfigs] {nullptr};
0813   book_shuffle_file(output_file_shuffle,
0814                     h_AN_pt_shuffle,
0815                     h_AN_pt_nodir_shuffle,
0816                     h_AN_eta_shuffle,
0817                     h_AN_xf_shuffle,
0818                     h_AN_pt_shuffle_details,
0819                     h_AN_pt_nodir_shuffle_details,
0820                     h_AN_eta_shuffle_details,
0821                     h_AN_xf_shuffle_details,
0822                     shuffle_file_name.str());
0823
0824   int shuffleMin = 0;
0825   int shuffleMax = 0;
0826   for (int iStore = 0; iStore < nStore; iStore++) {
0827
0828     shuffleMin = iStore * shuffleMaxStorage + shuffleIterMin;
0829     shuffleMax = std::min((iStore + 1) * shuffleMaxStorage - 1 + shuffleIterMin,
0830                           shuffleIterMax);
0831     std::cout << "shuffleMin, shuffleMax = " << shuffleMin << ", " << shuffleMax << std::endl;
0832
0833     AsymmetryCalc::ShuffleBunches *sb_mbd_0mrad = new AsymmetryCalc::ShuffleBunches(
0834         shuffleMin,
0835         shuffleMax,
0836         path_to_fill_list_0mrad,
0837         path_to_spin_pattern_0mrad,
0838         path_to_bunch_dependent_yields_mbd,
0839         input_csv_mbd_0mrad_pt_pi0_name,
0840         input_csv_mbd_0mrad_eta_pi0_name,
0841         input_csv_mbd_0mrad_xf_pi0_name,
0842         input_csv_mbd_0mrad_pt_eta_name,
0843         input_csv_mbd_0mrad_eta_eta_name,
0844         input_csv_mbd_0mrad_xf_eta_name,
0845         input_csv_mbd_0mrad_pt_ratios_pi0_name,
0846         input_csv_mbd_0mrad_eta_ratios_pi0_name,
0847         input_csv_mbd_0mrad_xf_ratios_pi0_name,
0848         input_csv_mbd_0mrad_pt_ratios_eta_name,
0849         input_csv_mbd_0mrad_eta_ratios_eta_name,
0850         input_csv_mbd_0mrad_xf_ratios_eta_name);
0851     sb_mbd_0mrad->set_shuffle(perform_shuffle);
0852     sb_mbd_0mrad->set_store_fill_histos(store_raw, (stored_iter_asym_histos + "/mbd_0mrad_"));
0853     sb_mbd_0mrad->run_1();
0854     sb_mbd_0mrad->run_2();
0855
0856     AsymmetryCalc::ShuffleBunches *sb_photon_0mrad = new AsymmetryCalc::ShuffleBunches(
0857         shuffleMin,
0858         shuffleMax,
0859         path_to_fill_list_0mrad,
0860         path_to_spin_pattern_0mrad,
0861         path_to_bunch_dependent_yields_photon,
0862         input_csv_photon_0mrad_pt_pi0_name,
0863         input_csv_photon_0mrad_eta_pi0_name,
0864         input_csv_photon_0mrad_xf_pi0_name,
0865         input_csv_photon_0mrad_pt_eta_name,
0866         input_csv_photon_0mrad_eta_eta_name,
0867         input_csv_photon_0mrad_xf_eta_name,
0868         input_csv_photon_0mrad_pt_ratios_pi0_name,
0869         input_csv_photon_0mrad_eta_ratios_pi0_name,
0870         input_csv_photon_0mrad_xf_ratios_pi0_name,
0871         input_csv_photon_0mrad_pt_ratios_eta_name,
0872         input_csv_photon_0mrad_eta_ratios_eta_name,
0873         input_csv_photon_0mrad_xf_ratios_eta_name);
0874     sb_photon_0mrad->set_shuffle(perform_shuffle);
0875     sb_photon_0mrad->set_store_fill_histos(store_raw, (stored_iter_asym_histos + "/photon_0mrad_"));
0876     sb_photon_0mrad->run_1();
0877     sb_photon_0mrad->run_2();
0878
0879     AsymmetryCalc::ShuffleBunches *sb_mbd_15mrad = new AsymmetryCalc::ShuffleBunches(
0880         shuffleMin,
0881         shuffleMax,
0882         path_to_fill_list_15mrad,
0883         path_to_spin_pattern_15mrad,
0884         path_to_bunch_dependent_yields_mbd,
0885         input_csv_mbd_15mrad_pt_pi0_name,
0886         input_csv_mbd_15mrad_eta_pi0_name,
0887         input_csv_mbd_15mrad_xf_pi0_name,
0888         input_csv_mbd_15mrad_pt_eta_name,
0889         input_csv_mbd_15mrad_eta_eta_name,
0890         input_csv_mbd_15mrad_xf_eta_name,
0891         input_csv_mbd_15mrad_pt_ratios_pi0_name,
0892         input_csv_mbd_15mrad_eta_ratios_pi0_name,
0893         input_csv_mbd_15mrad_xf_ratios_pi0_name,
0894         input_csv_mbd_15mrad_pt_ratios_eta_name,
0895         input_csv_mbd_15mrad_eta_ratios_eta_name,
0896         input_csv_mbd_15mrad_xf_ratios_eta_name);
0897     sb_mbd_15mrad->set_shuffle(perform_shuffle);
0898     sb_mbd_15mrad->set_store_fill_histos(store_raw, (stored_iter_asym_histos + "/mbd_15mrad_"));
0899     sb_mbd_15mrad->run_1();
0900     sb_mbd_15mrad->run_2();
0901
0902     AsymmetryCalc::ShuffleBunches *sb_photon_15mrad = new AsymmetryCalc::ShuffleBunches(
0903         shuffleMin,
0904         shuffleMax,
0905         path_to_fill_list_15mrad,
0906         path_to_spin_pattern_15mrad,
0907         path_to_bunch_dependent_yields_photon,
0908         input_csv_photon_15mrad_pt_pi0_name,
0909         input_csv_photon_15mrad_eta_pi0_name,
0910         input_csv_photon_15mrad_xf_pi0_name,
0911         input_csv_photon_15mrad_pt_eta_name,
0912         input_csv_photon_15mrad_eta_eta_name,
0913         input_csv_photon_15mrad_xf_eta_name,
0914         input_csv_photon_15mrad_pt_ratios_pi0_name,
0915         input_csv_photon_15mrad_eta_ratios_pi0_name,
0916         input_csv_photon_15mrad_xf_ratios_pi0_name,
0917         input_csv_photon_15mrad_pt_ratios_eta_name,
0918         input_csv_photon_15mrad_eta_ratios_eta_name,
0919         input_csv_photon_15mrad_xf_ratios_eta_name);
0920     sb_photon_15mrad->set_shuffle(perform_shuffle);
0921     sb_photon_15mrad->set_store_fill_histos(store_raw, (stored_iter_asym_histos + "/photon_15mrad_"));
0922     sb_photon_15mrad->run_1();
0923     sb_photon_15mrad->run_2();
0924
0925     int local_number_shuffles = (shuffleMax - shuffleMin + 1);
0926     for (int iIter = 0; iIter < local_number_shuffles; iIter++) {
0927       AsymmetryCalc::pt_final_array merged_mean_pt;
0928       AsymmetryCalc::pt_nodir_final_array merged_mean_pt_nodir;
0929       AsymmetryCalc::eta_final_array merged_mean_eta;
0930       AsymmetryCalc::xf_final_array merged_mean_xf;
0931       AsymmetryCalc::pt_final_array merged_unc_pt;
0932       AsymmetryCalc::pt_nodir_final_array merged_unc_pt_nodir;
0933       AsymmetryCalc::eta_final_array merged_unc_eta;
0934       AsymmetryCalc::xf_final_array merged_unc_xf;
0935       merge_configurations(
0936           sb_mbd_0mrad,
0937           sb_photon_0mrad,
0938           sb_mbd_15mrad,
0939           sb_photon_15mrad,
0940           merged_mean_pt,
0941           merged_mean_pt_nodir,
0942           merged_mean_eta,
0943           merged_mean_xf,
0944           merged_unc_pt,
0945           merged_unc_pt_nodir,
0946           merged_unc_eta,
0947           merged_unc_xf,
0948           iIter);
0949
0950       fill_shuffle_detailed_histos(
0951          sb_mbd_0mrad,
0952          sb_photon_0mrad,
0953          sb_mbd_15mrad,
0954          sb_photon_15mrad,
0955          h_AN_pt_shuffle_details,
0956          h_AN_pt_nodir_shuffle_details,
0957          h_AN_eta_shuffle_details,
0958          h_AN_xf_shuffle_details,
0959          iIter);
0960
0961       // Fill shuffle histogram
0962
0963       if (store_final) {
0964         std::string output_final_name = stored_final_asym_histos;
0965         save_final_histograms(
0966          output_final_name,
0967          merged_mean_pt,
0968          merged_mean_pt_nodir,
0969          merged_mean_eta,
0970          merged_mean_xf,
0971          merged_unc_pt,
0972          merged_unc_pt_nodir,
0973          merged_unc_eta,
0974          merged_unc_xf);
0975       }
0976
0977       // After filling shuffle histos, "merged" vectors are cleared automatically
0978       fill_shuffle_histos(merged_mean_pt,
0979                           merged_mean_pt_nodir,
0980                           merged_mean_eta,
0981                           merged_mean_xf,
0982                           merged_unc_pt,
0983                           merged_unc_pt_nodir,
0984                           merged_unc_eta,
0985                           merged_unc_xf,
0986                           h_AN_pt_shuffle,
0987                           h_AN_pt_nodir_shuffle,
0988                           h_AN_eta_shuffle,
0989                           h_AN_xf_shuffle);
0990     }
0991
0992     delete sb_mbd_0mrad;
0993     delete sb_photon_0mrad;
0994     delete sb_mbd_15mrad;
0995     delete sb_photon_15mrad;
0996   }
0997
0998   save_shuffle_histos(output_file_shuffle);
0999
1000   return 0;
1001 }