include/AsymmetryCalc/ShuffleBunches.h

0001 #ifndef __SHUFFLE_BUNCHES_H__
0002 #define __SHUFFLE_BUNCHES_H__
0003
0004 #include "AsymmetryCalc/Constants.h"
0005 #include "AsymmetryCalc/WrapperArrays.h"
0006 #include "AsymmetryCalc/WrapperAsymArrays.h"
0007 #include "AsymmetryCalc/WrapperCorrAsymArrays.h"
0008 #include "AsymmetryCalc/WrapperFitAsymArrays.h"
0009
0010
0011 #include <map>
0012 #include <vector>
0013
0014 #include <TFile.h>
0015 #include <TTree.h>
0016 #include <TH1F.h>
0017 #include <TF1.h>
0018 #include <TGraphErrors.h>
0019
0020 namespace AsymmetryCalc
0021 {
0022   class ShuffleBunches {
0023
0024   public:
0025
0026     ShuffleBunches(const int iMin,
0027                    const int iMax,
0028                    const std::string& inputfilename_fill_run,
0029                    const std::string& inputfilename_spin_pattern,
0030                    const std::string& inputfile_template,
0031                    const std::string& input_csv_pt_pi0_name,
0032                    const std::string& input_csv_eta_pi0_name,
0033                    const std::string& input_csv_xf_pi0_name,
0034                    const std::string& input_csv_pt_eta_name,
0035                    const std::string& input_csv_eta_eta_name,
0036                    const std::string& input_csv_xf_eta_name,
0037                    const std::string& input_csv_pt_ratios_pi0_name,
0038                    const std::string& input_csv_eta_ratios_pi0_name,
0039                    const std::string& input_csv_xf_ratios_pi0_name,
0040                    const std::string& input_csv_pt_ratios_eta_name,
0041                    const std::string& input_csv_eta_ratios_eta_name,
0042                    const std::string& input_csv_xf_ratios_eta_name);
0043
0044     void get_average_bin_pt_pi0(const std::string &input_csv_name);
0045     void get_average_bin_eta_pi0(const std::string &input_csv_name);
0046     void get_average_bin_xf_pi0(const std::string &input_csv_name);
0047
0048     void get_average_bin_pt_eta(const std::string &input_csv_name);
0049     void get_average_bin_eta_eta(const std::string &input_csv_name);
0050     void get_average_bin_xf_eta(const std::string &input_csv_name);
0051
0052     void get_pt_ratios_pi0(const std::string &input_csv_name);
0053     void get_eta_ratios_pi0(const std::string &input_csv_name);
0054     void get_xf_ratios_pi0(const std::string &input_csv_name);
0055
0056     void get_pt_ratios_eta(const std::string &input_csv_name);
0057     void get_eta_ratios_eta(const std::string &input_csv_name);
0058     void get_xf_ratios_eta(const std::string &input_csv_name);
0059
0060     ~ShuffleBunches();
0061
0062     void reserve_fill_yields();
0063     void reserve_fill_raw_asyms();
0064     void reserve_average_raw_asyms();
0065     void reserve_corr_asyms();
0066     void reserve_fit_asyms();
0067
0068     void clear_fill_yields();
0069     void clear_fill_raw_asyms();
0070     void clear_average_raw_asyms();
0071     void clear_corr_asyms();
0072     void clear_fit_asyms();
0073
0074     void reset();
0075
0076     //! Shuffle bunches and get raw asymmetries fill by fill
0077     void run_1();
0078
0079     //! Average raw asymmetries and get a single value of corrected asymmetry per iteration
0080     void run_2();
0081
0082     void compute_fill(const int ifill);
0083
0084     void read_binary_array(const std::string& inputfilename);
0085
0086     void check_array();
0087
0088     void shuffle_bunch_indices(const int seednumber);
0089
0090     void sum_yields(const int iIter);
0091
0092     void set_shuffle(bool val) { do_shuffle = val; }
0093
0094     void set_store_fill_histos(bool val, const std::string& outputfilename_template);
0095
0096     void set_store_iter_histos(bool val, const std::string& outputfilename_template);
0097
0098     void get_yield_histograms();
0099
0100     void get_fill_raw_asym_histograms();
0101
0102     void book_outfile_fill_histograms(const std::string& outputfilename);
0103
0104     void save_outfile_fill_histograms();
0105
0106     void get_iter_raw_asym_histograms();
0107
0108     void get_iter_corr_asym_histograms();
0109
0110     void get_iter_fit_asym_histograms();
0111
0112     void book_outfile_iter_histograms(const std::string& outputfilename);
0113
0114     void save_outfile_iter_histograms();
0115
0116     void compute_raw_asyms(const int iIter);
0117
0118     void average_asyms(const int iter);
0119
0120     void fit_raw_asyms(const int iter);
0121
0122     void fit_asym(double& asym_amplitude,
0123                   double& asym_amplitude_err,
0124                   const double (&asym_azimuth_values)[ASYM_CONSTANTS::nPhiBins],
0125                   const double (&asym_azimuth_uncertainties)[ASYM_CONSTANTS::nPhiBins],
0126                   bool fit_geom = true);
0127
0128     void compute_corrected_asyms(const int iter);
0129
0130     void compute_corrected_asym(double &corrected_asym,
0131                                 double &corrected_asym_err,
0132                                 const double &raw_asym,
0133                                 const double &bkg_asym,
0134                                 const double &raw_asym_err,
0135                                 const double &bkg_asym_err,
0136                                 const double &R,
0137                                 const double &R_err);
0138
0139     void compute_rellum_asym(const uint32_t (&array_yield)[ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins],
0140                              const double pol,
0141                              double (&array_asym)[ASYM_CONSTANTS::nPhiBins],
0142                              double (&array_asym_err)[ASYM_CONSTANTS::nPhiBins]);
0143
0144     void compute_sqrt_asym(const uint32_t (&array_yield)[ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins],
0145                            const double pol,
0146                            double (&array_asym)[ASYM_CONSTANTS::nPhiBins],
0147                            double (&array_asym_err)[ASYM_CONSTANTS::nPhiBins]);
0148
0149     const std::vector<pt_corr_array>& get_corrected_mean_pt() { return corrected_mean_pt; }
0150     const std::vector<pt_corr_array>& get_corrected_unc_pt() { return corrected_unc_pt; }
0151     const std::vector<pt_nodir_corr_array>& get_corrected_mean_pt_nodir() { return corrected_mean_pt_nodir; }
0152     const std::vector<pt_nodir_corr_array>& get_corrected_unc_pt_nodir() { return corrected_unc_pt_nodir; }
0153     const std::vector<eta_corr_array>& get_corrected_mean_eta() { return corrected_mean_eta; }
0154     const std::vector<eta_corr_array>& get_corrected_unc_eta() { return corrected_unc_eta; }
0155     const std::vector<xf_corr_array>& get_corrected_mean_xf() { return corrected_mean_xf; }
0156     const std::vector<xf_corr_array>& get_corrected_unc_xf() { return corrected_unc_xf; }
0157
0158   private:
0159
0160     int iterMin = 0;
0161     int iterMax = 0;
0162     int nIterations = 0;
0163     int nFills = 0;
0164
0165     std::string infile_template = "";
0166
0167     std::vector<std::pair<int, std::vector<int>>> fill_to_runs;
0168
0169     // Read the Spin Pattern
0170     TFile *inputfile_spin = nullptr;
0171     TTree *spin_patterns = nullptr;
0172     int fill_spin;
0173     double blue_polarization;
0174     double yellow_polarization;
0175     double polarization[2];
0176     int8_t yellow_spin_pattern[ASYM_CONSTANTS::nBunches];
0177     int8_t blue_spin_pattern[ASYM_CONSTANTS::nBunches];
0178     int8_t spin_pattern[2][ASYM_CONSTANTS::nBunches];
0179     int8_t shuffled_spin_pattern[2][ASYM_CONSTANTS::nBunches];
0180     double relative_luminosity = 1;
0181
0182     bool do_shuffle = false;
0183     int shuffled_indices[ASYM_CONSTANTS::nBunches] = {0};
0184     bool spin_flip = false;
0185     int global_irun = 0;
0186
0187     uint32_t array_yield_pt_bunches[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins][120];
0188     uint32_t array_yield_eta_bunches[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins][120];
0189     uint32_t array_yield_xf_bunches[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins][120];
0190
0191     uint32_t array_yield_pt[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins] = {0};
0192     uint32_t array_yield_eta[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins] = {0};
0193     uint32_t array_yield_xf[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nSpins][ASYM_CONSTANTS::nPhiBins] = {0};
0194
0195     std::vector<pt_yield_array> vec_yield_pt;
0196     std::vector<eta_yield_array> vec_yield_eta;
0197     std::vector<xf_yield_array> vec_yield_xf;
0198
0199     std::vector<pt_asym_array> fill_pt_asyms;
0200     std::vector<pt_nodir_asym_array> fill_pt_nodir_asyms;
0201     std::vector<eta_asym_array> fill_eta_asyms;
0202     std::vector<xf_asym_array> fill_xf_asyms;
0203
0204     std::vector<pt_asym_array> fill_pt_asym_errs;
0205     std::vector<pt_nodir_asym_array> fill_pt_nodir_asym_errs;
0206     std::vector<eta_asym_array> fill_eta_asym_errs;
0207     std::vector<xf_asym_array> fill_xf_asym_errs;
0208
0209     std::vector<pt_asym_array> mean_pt;
0210     std::vector<pt_nodir_asym_array> mean_pt_nodir;
0211     std::vector<eta_asym_array> mean_eta;
0212     std::vector<xf_asym_array> mean_xf;
0213
0214     std::vector<pt_asym_array> unc_pt;
0215     std::vector<pt_nodir_asym_array> unc_pt_nodir;
0216     std::vector<eta_asym_array> unc_eta;
0217     std::vector<xf_asym_array> unc_xf;
0218
0219     std::vector<pt_fit_array> fit_mean_pt;
0220     std::vector<pt_nodir_fit_array> fit_mean_pt_nodir;
0221     std::vector<eta_fit_array> fit_mean_eta;
0222     std::vector<xf_fit_array> fit_mean_xf;
0223
0224     std::vector<pt_fit_array> fit_unc_pt;
0225     std::vector<pt_nodir_fit_array> fit_unc_pt_nodir;
0226     std::vector<eta_fit_array> fit_unc_eta;
0227     std::vector<xf_fit_array> fit_unc_xf;
0228
0229     std::vector<pt_corr_array> corrected_mean_pt;
0230     std::vector<pt_nodir_corr_array> corrected_mean_pt_nodir;
0231     std::vector<eta_corr_array> corrected_mean_eta;
0232     std::vector<xf_corr_array> corrected_mean_xf;
0233
0234     std::vector<pt_corr_array> corrected_unc_pt;
0235     std::vector<pt_nodir_corr_array> corrected_unc_pt_nodir;
0236     std::vector<eta_corr_array> corrected_unc_eta;
0237     std::vector<xf_corr_array> corrected_unc_xf;
0238
0239     // Mean kinematic values
0240     double pTMeans[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins] = {0};
0241     double etaMeans[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins] = {0};
0242     double xfMeans[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins] = {0};
0243
0244     // background ratios
0245     double bkg_ratio_pt[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins] = {0};
0246     double bkg_ratio_err_pt[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nPtBins] = {0};
0247     double bkg_ratio_eta[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins] = {0};
0248     double bkg_ratio_err_eta[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nEtaBins] = {0};
0249     double bkg_ratio_xf[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins] = {0};
0250     double bkg_ratio_err_xf[ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nXfBins] = {0};
0251
0252     // Store fill-dependent histograms (NOT for bunch shuffling)
0253     bool store_fill_histos = false;
0254     std::string outfilename_fill_template = "";
0255     TFile *outfile_fill_histograms = nullptr;
0256
0257     bool store_iter_histos = false;
0258     std::string outfilename_iter_template = "";
0259     TFile *outfile_iter_histograms = nullptr;
0260
0261     // Verification histograms
0262     TH1F *h_yield_pt[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nPtBins][ASYM_CONSTANTS::nDirections][ASYM_CONSTANTS::nSpins] = {nullptr}; // forward vs backward
0263     TH1F *h_yield_eta[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nEtaBins][ASYM_CONSTANTS::nSpins] = {nullptr};
0264     TH1F *h_yield_xf[ASYM_CONSTANTS::nBeams][ASYM_CONSTANTS::nParticles][ASYM_CONSTANTS::nRegions][ASYM_CONSTANTS::nXfBins][ASYM_CONSTANTS::nSpins] = {nullptr};
0265
0266   };
0267 };
0268
0269 #endif