From d3cd2749ffd852bbd54d7c8539ce05586e4e298f Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 20 Oct 2014 14:02:10 -0400 Subject: [PATCH 01/19] implement histogram class --- b3d2/src/hist.cc | 20 ++++++++++++++++++++ b3d2/src/hist.h | 14 ++++++++++++++ b3d2/src/test_hist.cc | 38 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 b3d2/src/hist.cc create mode 100644 b3d2/src/hist.h create mode 100644 b3d2/src/test_hist.cc diff --git a/b3d2/src/hist.cc b/b3d2/src/hist.cc new file mode 100644 index 0000000..17b875e --- /dev/null +++ b/b3d2/src/hist.cc @@ -0,0 +1,20 @@ +#include +#include "hist.h" + +Histogram::Histogram(size_t nbins, double x_min, double x_max) { + this->x_min = x_min; + this->x_max = x_max; + this->nbins = nbins; + this->histogram.resize(nbins); +} + +void Histogram::add(double x) { + const double bin_float = (x - x_min) / ((x_max - x_min) / nbins); + const size_t bin_int = (size_t) bin_float; + //printf("x=%10g, bin=%10i\n", x, bin); + if ( (bin_float >= 0) && (bin_int < nbins) ) + histogram[bin_int]++; + else + printf("warning: %g out bounds of histogram!\n", x); +} + diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h new file mode 100644 index 0000000..cdc2d29 --- /dev/null +++ b/b3d2/src/hist.h @@ -0,0 +1,14 @@ +#include + +class Histogram { + public: + double x_min; + double x_max; + size_t nbins; + std::vector histogram; + + /// Initialize histogram. + Histogram(size_t nbins, double x_min, double x_max); + /// Add value x to given histogram. + void add(double x); +}; diff --git a/b3d2/src/test_hist.cc b/b3d2/src/test_hist.cc new file mode 100644 index 0000000..ab0ee29 --- /dev/null +++ b/b3d2/src/test_hist.cc @@ -0,0 +1,38 @@ +#include +#include +#include +#include "hist.h" + +double gaussrand() { + static double V1, V2, S; + static int phase = 0; + double X; + + if (phase == 0) { + do { + double U1 = (double) rand() / RAND_MAX; + double U2 = (double) rand() / RAND_MAX; + + V1 = 2 * U1 - 1; + V2 = 2 * U2 - 1; + S = V1 * V1 + V2 * V2; + } while (S >= 1 || S == 0); + + X = V1 * sqrt(-2 * log(S) / S); + } else + X = V2 * sqrt(-2 * log(S) / S); + + phase = 1 - phase; + + return X; +} + +int main() { + const size_t nbins = 100; + auto h = Histogram(nbins, -5, 5); + for (size_t i = 0; i < 100000000; i++) + h.add(gaussrand()); + for (size_t i = 0; i < nbins; i++) + std::cout << h.histogram[i] << " "; + std::cout << std::endl; +} From 254688abdf3c72b3a2fc6cf55f21d4e906d72af0 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 20 Oct 2014 18:27:45 -0400 Subject: [PATCH 02/19] make histogram more generic --- b3d2/src/hist.cc | 20 ----------------- b3d2/src/hist.h | 50 ++++++++++++++++++++++++++++++++++++++----- b3d2/src/test_hist.cc | 3 ++- 3 files changed, 47 insertions(+), 26 deletions(-) delete mode 100644 b3d2/src/hist.cc diff --git a/b3d2/src/hist.cc b/b3d2/src/hist.cc deleted file mode 100644 index 17b875e..0000000 --- a/b3d2/src/hist.cc +++ /dev/null @@ -1,20 +0,0 @@ -#include -#include "hist.h" - -Histogram::Histogram(size_t nbins, double x_min, double x_max) { - this->x_min = x_min; - this->x_max = x_max; - this->nbins = nbins; - this->histogram.resize(nbins); -} - -void Histogram::add(double x) { - const double bin_float = (x - x_min) / ((x_max - x_min) / nbins); - const size_t bin_int = (size_t) bin_float; - //printf("x=%10g, bin=%10i\n", x, bin); - if ( (bin_float >= 0) && (bin_int < nbins) ) - histogram[bin_int]++; - else - printf("warning: %g out bounds of histogram!\n", x); -} - diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h index cdc2d29..5baedd4 100644 --- a/b3d2/src/hist.h +++ b/b3d2/src/hist.h @@ -1,14 +1,54 @@ #include +#include +#include +template class Histogram { public: - double x_min; - double x_max; + T x_min; + T x_max; size_t nbins; - std::vector histogram; + std::vector histogram; /// Initialize histogram. - Histogram(size_t nbins, double x_min, double x_max); + Histogram(size_t nbins, T x_min, T x_max); /// Add value x to given histogram. - void add(double x); + void add(T x, int64_t count=1); + /// Return the count of the bin corresponding to x. + int64_t get_count(T x); + private: + double get_bin(T x); }; + +template +Histogram::Histogram(size_t nbins, T x_min, T x_max) { + this->x_min = x_min; + this->x_max = x_max; + this->nbins = nbins; + this->histogram.resize(nbins); // zero by default +} + +template +double Histogram::get_bin(T x) { + return (x - x_min) / ((x_max - x_min) / nbins); +} + +template +void Histogram::add(T x, int64_t count) { + const double bin_float = get_bin(x); + const size_t bin_int = (size_t) bin_float; + //printf("x=%10g, bin=%10i\n", x, bin_float); + if ( (bin_float >= 0) && (bin_int < nbins) ) + histogram[bin_int] += count; + else + printf("warning: %g out bounds of histogram!\n", x); +} + +template +int64_t Histogram::get_count(T x) { + const double bin_float = get_bin(x); + assert(bin_float >= 0); + const size_t bin_int = (size_t) bin_float; + assert(bin_int < nbins); + return histogram[bin_int]; +} diff --git a/b3d2/src/test_hist.cc b/b3d2/src/test_hist.cc index ab0ee29..81285a7 100644 --- a/b3d2/src/test_hist.cc +++ b/b3d2/src/test_hist.cc @@ -29,10 +29,11 @@ double gaussrand() { int main() { const size_t nbins = 100; - auto h = Histogram(nbins, -5, 5); + auto h = Histogram(nbins, -5, 5); for (size_t i = 0; i < 100000000; i++) h.add(gaussrand()); for (size_t i = 0; i < nbins; i++) std::cout << h.histogram[i] << " "; std::cout << std::endl; + std::cout << "\n" << h.get_count(0) << std::endl; } From 8f29ec896dbbee9f1d3870a65b8ff8412f83b8c8 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 21 Oct 2014 11:08:11 -0400 Subject: [PATCH 03/19] implement calculation of balance function --- b3d2/makefile | 4 ++ b3d2/src/anal_balance.cc | 105 +++++++++++++++++++++++++++++++++++++++ b3d2/src/b3d.h | 3 +- 3 files changed, 111 insertions(+), 1 deletion(-) create mode 100644 b3d2/src/anal_balance.cc diff --git a/b3d2/makefile b/b3d2/makefile index 5d84910..a4ccb34 100644 --- a/b3d2/makefile +++ b/b3d2/makefile @@ -79,6 +79,7 @@ build/part.o\ build/inelastic.o\ build/resonances.o\ build/anal_spectra.o\ +build/anal_balance.o\ build/anal_hbt.o\ build/anal_v2.o\ build/consolidate.o\ @@ -215,6 +216,9 @@ build/inelastic.o : src/inelastic.cc ${B3D_HFILES} build/anal_spectra.o : src/anal_spectra.cc ${B3D_HFILES} ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/anal_spectra.o src/anal_spectra.cc + +build/anal_balance.o : src/anal_balance.cc ${B3D_HFILES} + ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/anal_balance.o src/anal_balance.cc build/anal_hbt.o : src/anal_hbt.cc ${B3D_HFILES} ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/anal_hbt.o src/anal_hbt.cc diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc new file mode 100644 index 0000000..df5667f --- /dev/null +++ b/b3d2/src/anal_balance.cc @@ -0,0 +1,105 @@ +#include +#include +#include + +#include "b3d.h" +#include "hist.h" + +/// Represent all balance functions. +/// { (PID, PID): B(0, eta), ... } +//typedef unordered_map, Histogram> BalanceMap; + +/// Calculate the balance functions. +void CB3D::CalcBalance() { + NSAMPLE = parameter::getI(parmap, "B3D_NSAMPLE", 1); + const int neventsmax = parameter::getI(parmap, "B3D_NEVENTSMAX", 10); + + const string dirname = "analysis/" + run_name + "/" + qualifier; + const string command = "mkdir -p " + dirname; + system(command.c_str()); + const string anal_output_filename = dirname + "/results_balance.dat"; + auto anal_output = fopen(anal_output_filename.c_str(), "w"); + + parameter::set(parmap, string("B3D_RESONANCES_DECAYS_FILE"), string("progdata/madai/resinfo/decays_pdg_weak.dat")); + parameter::set(parmap, string("B3D_RESONANCES_INFO_FILE"), string("progdata/madai/resinfo/resonances_pdg_weak.dat")); + reslist = new CResList(); + reslist->ReadResInfo(); + + const double DELRAPIDITY = 2.0 * parameter::getD(parmap, "B3D_ETAMAX", 1.0); + COLLISIONS = false; + + /// Balance functions to calculate. + vector> balance_pairs = { {211, -211}, {321, -321} }; + /// Relevant species. + unordered_set balance_species; + for (const auto& pair : balance_pairs) { + balance_species.insert(pair.first); + balance_species.insert(-pair.first); + balance_species.insert(pair.second); + balance_species.insert(-pair.second); + } + /// A histogram for each species + unordered_map> histograms; + histograms.reserve(balance_species.size()); + // TODO: figure out what the following constants should be. + const size_t nbins = 100; + const double min_eta = -2; + const double max_eta = 2; + for (const auto& PID : balance_species) { + histograms.emplace(make_pair(PID, Histogram(nbins, min_eta, max_eta))); + } + + int nevents = 0; + do { + KillAllParts(); + const int nparts = ReadOSCAR(nevents+1); + if (nparts > 0){ + nevents += 1; + //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); + PerformAllActions(); // Decays unstable particles + //printf("after decays, nparts=%d\n",int(FinalPartMap.size())); + + // We need a histogram in eta for each species and antispecies. + for (const auto& ppos : FinalPartMap) { + const auto part = ppos.second; + const int PID = part->resinfo->code; + const auto search = histograms.find(PID); + if (search != histograms.end()) { + const double eta = part->eta; + const int weight = part->weight; + const int charge = abs(part->resinfo->charge); + search->second.add(eta, weight*charge); + } + } + } + } while (!feof(oscarfile) && nevents < neventsmax); + if (nevents != neventsmax) { + printf("EVENT SHORTAGE?\n"); + } + fclose(oscarfile); + oscarfile = NULL; + + // Calculate balance function. + // B_ab(eta) = (N_a(0) - N_-a(0))(N_b(eta) - N_-b(eta)) / (N_b(eta) + N_-b(eta)) + for (const auto& mypair : balance_pairs) { + vector B(nbins); + const int a = mypair.first; + const int b = mypair.second; + for (size_t i = 0; i < B.size(); i++) { + const int Na_pos = histograms.at(a).get_count(0); + const int Na_neg = histograms.at(-a).get_count(0); + const int Nb_pos = histograms.at(b).histogram[i]; + const int Nb_neg = histograms.at(-b).histogram[i]; + B[i] = (Na_pos - Na_neg) * (Nb_pos - Nb_neg) / (Nb_pos + Nb_neg); + } + fprintf(anal_output, "B(%i, %i)\n", a, b); + fprintf(anal_output, "eta B\n"); + for (size_t i = 0; i < B.size(); i++) { + const double eta = min_eta + (max_eta - min_eta) / nbins * (i + 0.5); + fprintf(anal_output, "%f %f", eta, B[i]); + } + } + + delete reslist; + fclose(anal_output); +} diff --git a/b3d2/src/b3d.h b/b3d2/src/b3d.h index 9c71ec5..fb9638d 100644 --- a/b3d2/src/b3d.h +++ b/b3d2/src/b3d.h @@ -175,6 +175,7 @@ class CB3D{ void CalcSpectra_PHENIXppbar(); double CalcSpectra_ALICE(); void Calc3DSpectra(); + void CalcBalance(); void CalcV2_STAR(); void CalcV2_ALICE(); void CalcHBT_STAR(); @@ -238,4 +239,4 @@ class CAction{ -#endif \ No newline at end of file +#endif From 8030ed6be8f5218c3551b836d08a9ec1e313f9a7 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 23 Oct 2014 18:01:17 -0400 Subject: [PATCH 04/19] calculate balance function correctly --- b3d2/src/anal_balance.cc | 87 ++++++++++++++++++++++++++-------------- 1 file changed, 58 insertions(+), 29 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index df5667f..512943e 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -2,6 +2,8 @@ #include #include +#include // needed for hashing pairs + #include "b3d.h" #include "hist.h" @@ -30,47 +32,78 @@ void CB3D::CalcBalance() { /// Balance functions to calculate. vector> balance_pairs = { {211, -211}, {321, -321} }; - /// Relevant species. + /// Relevant species (ignore anti particles). unordered_set balance_species; for (const auto& pair : balance_pairs) { - balance_species.insert(pair.first); - balance_species.insert(-pair.first); - balance_species.insert(pair.second); - balance_species.insert(-pair.second); + balance_species.insert(abs(pair.first)); + balance_species.insert(abs(pair.second)); } - /// A histogram for each species - unordered_map> histograms; - histograms.reserve(balance_species.size()); + /// A rapidity histogram for each balance function. + unordered_map, Histogram, boost::hash>> histograms; + histograms.reserve(balance_pairs.size()); // TODO: figure out what the following constants should be. const size_t nbins = 100; - const double min_eta = -2; - const double max_eta = 2; + const double min_y = -1; + const double max_y = 1; + const double max_dy = max_y - min_y; + const double min_dy = 0; + for (const auto& pair : balance_pairs) { + histograms.emplace(make_pair(pair, Histogram(nbins, min_dy, max_dy))); + } + /// Number for each species. + /// (Needed for denominator of balance function.) + unordered_map total_number(balance_species.size()); for (const auto& PID : balance_species) { - histograms.emplace(make_pair(PID, Histogram(nbins, min_eta, max_eta))); + total_number[PID] = 0; } int nevents = 0; do { KillAllParts(); - const int nparts = ReadOSCAR(nevents+1); - if (nparts > 0){ + const int nparts = ReadOSCAR(nevents + 1); + if (nparts > 0) { nevents += 1; //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); PerformAllActions(); // Decays unstable particles //printf("after decays, nparts=%d\n",int(FinalPartMap.size())); - // We need a histogram in eta for each species and antispecies. + // We are only interested in particles occuring in the balance function + // and with y_min <= y <= y_max. + vector relevant_particles; // TODO: only remeber weight, charge, rapidity and momentum for (const auto& ppos : FinalPartMap) { const auto part = ppos.second; const int PID = part->resinfo->code; - const auto search = histograms.find(PID); + const int charge = part->resinfo->charge; + const int weight = part->weight; + const double y = part->y; + // Skip particles we do not care about. + if (balance_species.count(abs(PID)) == 0 || y < min_y || y > max_y) + continue; + relevant_particles.push_back(*part); + total_number[abs(PID)] += weight * abs(charge); + } + + // Iterate over all pairs to calculate dy and fill histograms. + for (size_t i = 0; i < relevant_particles.size(); i++) for (size_t j = 0; j < i; j++) { + const auto& p_i = relevant_particles[i]; + const auto& p_j = relevant_particles[j]; + const double dy = fabs(p_i.y - p_j.y); + const auto mypair = make_pair(p_i.resinfo->code, p_j.resinfo->code); + const auto mypair_reversed = make_pair(mypair.second, mypair.first); + + auto search = histograms.find(mypair); + if (search != histograms.end()) { + const int charge = p_j.resinfo->charge; + const int weight = p_j.weight; + search->second.add(dy, weight * charge); + } + search = histograms.find(mypair_reversed); if (search != histograms.end()) { - const double eta = part->eta; - const int weight = part->weight; - const int charge = abs(part->resinfo->charge); - search->second.add(eta, weight*charge); + const int charge = p_i.resinfo->charge; + const int weight = p_i.weight; + search->second.add(dy, weight * charge); } - } + } } } while (!feof(oscarfile) && nevents < neventsmax); if (nevents != neventsmax) { @@ -80,23 +113,19 @@ void CB3D::CalcBalance() { oscarfile = NULL; // Calculate balance function. - // B_ab(eta) = (N_a(0) - N_-a(0))(N_b(eta) - N_-b(eta)) / (N_b(eta) + N_-b(eta)) + // B_ab(y) = (N_a(0) - N_-a(0))(N_b(y) - N_-b(y)) / (N_b(y) + N_-b(y)) for (const auto& mypair : balance_pairs) { vector B(nbins); const int a = mypair.first; const int b = mypair.second; for (size_t i = 0; i < B.size(); i++) { - const int Na_pos = histograms.at(a).get_count(0); - const int Na_neg = histograms.at(-a).get_count(0); - const int Nb_pos = histograms.at(b).histogram[i]; - const int Nb_neg = histograms.at(-b).histogram[i]; - B[i] = (Na_pos - Na_neg) * (Nb_pos - Nb_neg) / (Nb_pos + Nb_neg); + B[i] = (double) histograms.find(mypair)->second.histogram[i] / total_number[abs(b)]; } fprintf(anal_output, "B(%i, %i)\n", a, b); - fprintf(anal_output, "eta B\n"); + fprintf(anal_output, "y B\n"); for (size_t i = 0; i < B.size(); i++) { - const double eta = min_eta + (max_eta - min_eta) / nbins * (i + 0.5); - fprintf(anal_output, "%f %f", eta, B[i]); + const double y = min_y + (max_y - min_y) / nbins * (i + 0.5); + fprintf(anal_output, "%f %f\n", y, B[i]); } } From 876b40bfd9cfb3b3da47b1f4c316eee3a867df61 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Fri, 24 Oct 2014 10:16:16 -0400 Subject: [PATCH 05/19] calculate balance function ignoring species --- b3d2/src/anal_balance.cc | 40 ++++++++++++++++++++++++++++++++++++---- b3d2/src/hist.h | 2 ++ b3d2/src/pow.h | 15 +++++++++++++++ 3 files changed, 53 insertions(+), 4 deletions(-) create mode 100644 b3d2/src/pow.h diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 512943e..ff01de2 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -6,10 +6,11 @@ #include "b3d.h" #include "hist.h" +#include "pow.h" -/// Represent all balance functions. -/// { (PID, PID): B(0, eta), ... } -//typedef unordered_map, Histogram> BalanceMap; +int sign(double x) { + return (x > 0) - (x < 0); +} /// Calculate the balance functions. void CB3D::CalcBalance() { @@ -50,12 +51,18 @@ void CB3D::CalcBalance() { for (const auto& pair : balance_pairs) { histograms.emplace(make_pair(pair, Histogram(nbins, min_dy, max_dy))); } + /// A pseudo-rapidity histogram for the charge-only balance function. + // TODO: use something other than min_dy/max_dy? + auto charge_histogram = Histogram(nbins, min_dy, max_dy); + /// Number for each species. /// (Needed for denominator of balance function.) unordered_map total_number(balance_species.size()); for (const auto& PID : balance_species) { total_number[PID] = 0; } + /// Total number including all species. + int64_t total_number_all = 0; int nevents = 0; do { @@ -79,8 +86,10 @@ void CB3D::CalcBalance() { // Skip particles we do not care about. if (balance_species.count(abs(PID)) == 0 || y < min_y || y > max_y) continue; + // TODO: p_T cuts and efficiency relevant_particles.push_back(*part); total_number[abs(PID)] += weight * abs(charge); + total_number_all += weight * abs(charge); } // Iterate over all pairs to calculate dy and fill histograms. @@ -91,6 +100,7 @@ void CB3D::CalcBalance() { const auto mypair = make_pair(p_i.resinfo->code, p_j.resinfo->code); const auto mypair_reversed = make_pair(mypair.second, mypair.first); + // For B(ab). auto search = histograms.find(mypair); if (search != histograms.end()) { const int charge = p_j.resinfo->charge; @@ -103,6 +113,19 @@ void CB3D::CalcBalance() { const int weight = p_i.weight; search->second.add(dy, weight * charge); } + + // For B(+-). + if (sign(p_i.resinfo->charge) == -sign(p_j.resinfo->charge)) { + const double dpseudorapidity = fabs( + atanh(p_i.p[3]/sqrt(square(p_i.p[1]) + square(p_i.p[2]) + square(p_i.p[3]))) - + atanh(p_i.p[3]/sqrt(square(p_j.p[1]) + square(p_j.p[2]) + square(p_j.p[3]))) + ); + // We want to calculate B(+-), so we want the negative charge. + const auto& p = p_i.resinfo->charge < 0 ? p_i : p_j; + const int charge = p.resinfo->charge; + const int weight = p.weight; + charge_histogram.add(dpseudorapidity, weight * charge); + } } } } while (!feof(oscarfile) && nevents < neventsmax); @@ -112,7 +135,7 @@ void CB3D::CalcBalance() { fclose(oscarfile); oscarfile = NULL; - // Calculate balance function. + // Calculate balance functions. // B_ab(y) = (N_a(0) - N_-a(0))(N_b(y) - N_-b(y)) / (N_b(y) + N_-b(y)) for (const auto& mypair : balance_pairs) { vector B(nbins); @@ -127,6 +150,15 @@ void CB3D::CalcBalance() { const double y = min_y + (max_y - min_y) / nbins * (i + 0.5); fprintf(anal_output, "%f %f\n", y, B[i]); } + fprintf(anal_output, "\n"); + } + + vector B_charge(nbins); + fprintf(anal_output, "B(+, -)\n"); + for (size_t i = 0; i < B_charge.size(); i++) { + B_charge[i] = (double) charge_histogram.histogram[i] / total_number_all; + const double pseudorapidity = min_y + (max_y - min_y) / nbins * (i + 0.5); + fprintf(anal_output, "%f %f\n", pseudorapidity, B_charge[i]); } delete reslist; diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h index 5baedd4..1fedfb2 100644 --- a/b3d2/src/hist.h +++ b/b3d2/src/hist.h @@ -1,3 +1,5 @@ +#pragma once + #include #include #include diff --git a/b3d2/src/pow.h b/b3d2/src/pow.h new file mode 100644 index 0000000..0d56b61 --- /dev/null +++ b/b3d2/src/pow.h @@ -0,0 +1,15 @@ +/// Efficient templates for calculating integer powers. +#pragma once + +/// Calculate integer powers using squaring. +template +inline constexpr T pow(const T base, unsigned const exponent) { + return (exponent == 0) ? 1 : + (exponent % 2 == 0) ? pow(base, exponent/2) * pow(base, exponent/2) : + base * pow(base, exponent - 1); +} + +template +inline constexpr T square(const T base) { + return pow(base, 2); +} From 822b2d708d5e872dd4a6fa2cae1b4456d67d114f Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 28 Oct 2014 15:34:12 -0400 Subject: [PATCH 06/19] check for eta bounds only if BJORKEN is set --- b3d2/src/part.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/b3d2/src/part.cc b/b3d2/src/part.cc index 220d376..1b13bbb 100644 --- a/b3d2/src/part.cc +++ b/b3d2/src/part.cc @@ -102,7 +102,7 @@ void CPart::Init(int IDset,double rxset,double ryset,double tauset,double etaset nscatt=0; weight=weightset; reality=realset; - if(fabs(eta)>b3d->ETAMAX){ + if(fabs(eta)>b3d->ETAMAX && b3d->BJORKEN){ printf("in part->init, eta out of bounds, =%g\n",eta); } resinfoptr=b3d->reslist->GetResInfoPtr(ID); From acedc9408e22b22ca0824d9db93815d40c4eec7a Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 28 Oct 2014 16:01:13 -0400 Subject: [PATCH 07/19] anal_balance: fix calculation of pseudorapidity A wrong index was used, which resulted in nans. Also made sure the code works for $p = 0$. --- b3d2/src/anal_balance.cc | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index ff01de2..a1dc774 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -14,6 +14,7 @@ int sign(double x) { /// Calculate the balance functions. void CB3D::CalcBalance() { + std::cout << "*** Calculating balance functions..." << std::endl; NSAMPLE = parameter::getI(parmap, "B3D_NSAMPLE", 1); const int neventsmax = parameter::getI(parmap, "B3D_NEVENTSMAX", 10); @@ -30,6 +31,8 @@ void CB3D::CalcBalance() { const double DELRAPIDITY = 2.0 * parameter::getD(parmap, "B3D_ETAMAX", 1.0); COLLISIONS = false; + oscarfile = NULL; + //^ We need to set this, so the header gets read. Otherwise `feof(oscarfile) == true`. /// Balance functions to calculate. vector> balance_pairs = { {211, -211}, {321, -321} }; @@ -68,6 +71,7 @@ void CB3D::CalcBalance() { do { KillAllParts(); const int nparts = ReadOSCAR(nevents + 1); + std::cout << nparts << std::endl; if (nparts > 0) { nevents += 1; //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); @@ -116,10 +120,15 @@ void CB3D::CalcBalance() { // For B(+-). if (sign(p_i.resinfo->charge) == -sign(p_j.resinfo->charge)) { - const double dpseudorapidity = fabs( - atanh(p_i.p[3]/sqrt(square(p_i.p[1]) + square(p_i.p[2]) + square(p_i.p[3]))) - - atanh(p_i.p[3]/sqrt(square(p_j.p[1]) + square(p_j.p[2]) + square(p_j.p[3]))) - ); + double p_i_mag = sqrt(square(p_i.p[1]) + square(p_i.p[2]) + square(p_i.p[3])); + double p_j_mag = sqrt(square(p_j.p[1]) + square(p_j.p[2]) + square(p_j.p[3])); + // Avoid dividing by zero (which would result in nans). + if (p_i_mag == 0) + p_i_mag = 1; + if (p_j_mag == 0) + p_j_mag = 1; + const double dpseudorapidity = fabs( atanh(p_i.p[3] / p_i_mag) - atanh(p_j.p[3] / p_j_mag ) ); + assert(dpseudorapidity == dpseudorapidity); // We want to calculate B(+-), so we want the negative charge. const auto& p = p_i.resinfo->charge < 0 ? p_i : p_j; const int charge = p.resinfo->charge; @@ -130,7 +139,8 @@ void CB3D::CalcBalance() { } } while (!feof(oscarfile) && nevents < neventsmax); if (nevents != neventsmax) { - printf("EVENT SHORTAGE?\n"); + printf("EVENT SHORTAGE? Got %i, expected %i.\n", nevents, neventsmax); + std::cout << oscarfile << " vs. " << feof(oscarfile) << std::endl; } fclose(oscarfile); oscarfile = NULL; From 6e4c0ed7991c48abc07053ac9c8ccd55672fd369 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 28 Oct 2014 16:08:31 -0400 Subject: [PATCH 08/19] replace tabs by 4 spaces --- b3d2/src/anal_balance.cc | 68 ++++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index a1dc774..6c55422 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -16,21 +16,21 @@ int sign(double x) { void CB3D::CalcBalance() { std::cout << "*** Calculating balance functions..." << std::endl; NSAMPLE = parameter::getI(parmap, "B3D_NSAMPLE", 1); - const int neventsmax = parameter::getI(parmap, "B3D_NEVENTSMAX", 10); + const int neventsmax = parameter::getI(parmap, "B3D_NEVENTSMAX", 10); const string dirname = "analysis/" + run_name + "/" + qualifier; - const string command = "mkdir -p " + dirname; - system(command.c_str()); - const string anal_output_filename = dirname + "/results_balance.dat"; - auto anal_output = fopen(anal_output_filename.c_str(), "w"); - - parameter::set(parmap, string("B3D_RESONANCES_DECAYS_FILE"), string("progdata/madai/resinfo/decays_pdg_weak.dat")); - parameter::set(parmap, string("B3D_RESONANCES_INFO_FILE"), string("progdata/madai/resinfo/resonances_pdg_weak.dat")); - reslist = new CResList(); - reslist->ReadResInfo(); - - const double DELRAPIDITY = 2.0 * parameter::getD(parmap, "B3D_ETAMAX", 1.0); - COLLISIONS = false; + const string command = "mkdir -p " + dirname; + system(command.c_str()); + const string anal_output_filename = dirname + "/results_balance.dat"; + auto anal_output = fopen(anal_output_filename.c_str(), "w"); + + parameter::set(parmap, string("B3D_RESONANCES_DECAYS_FILE"), string("progdata/madai/resinfo/decays_pdg_weak.dat")); + parameter::set(parmap, string("B3D_RESONANCES_INFO_FILE"), string("progdata/madai/resinfo/resonances_pdg_weak.dat")); + reslist = new CResList(); + reslist->ReadResInfo(); + + const double DELRAPIDITY = 2.0 * parameter::getD(parmap, "B3D_ETAMAX", 1.0); + COLLISIONS = false; oscarfile = NULL; //^ We need to set this, so the header gets read. Otherwise `feof(oscarfile) == true`. @@ -67,23 +67,23 @@ void CB3D::CalcBalance() { /// Total number including all species. int64_t total_number_all = 0; - int nevents = 0; - do { - KillAllParts(); - const int nparts = ReadOSCAR(nevents + 1); + int nevents = 0; + do { + KillAllParts(); + const int nparts = ReadOSCAR(nevents + 1); std::cout << nparts << std::endl; - if (nparts > 0) { - nevents += 1; - //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); - PerformAllActions(); // Decays unstable particles - //printf("after decays, nparts=%d\n",int(FinalPartMap.size())); + if (nparts > 0) { + nevents += 1; + //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); + PerformAllActions(); // Decays unstable particles + //printf("after decays, nparts=%d\n",int(FinalPartMap.size())); // We are only interested in particles occuring in the balance function // and with y_min <= y <= y_max. vector relevant_particles; // TODO: only remeber weight, charge, rapidity and momentum - for (const auto& ppos : FinalPartMap) { - const auto part = ppos.second; - const int PID = part->resinfo->code; + for (const auto& ppos : FinalPartMap) { + const auto part = ppos.second; + const int PID = part->resinfo->code; const int charge = part->resinfo->charge; const int weight = part->weight; const double y = part->y; @@ -136,14 +136,14 @@ void CB3D::CalcBalance() { charge_histogram.add(dpseudorapidity, weight * charge); } } - } - } while (!feof(oscarfile) && nevents < neventsmax); - if (nevents != neventsmax) { - printf("EVENT SHORTAGE? Got %i, expected %i.\n", nevents, neventsmax); + } + } while (!feof(oscarfile) && nevents < neventsmax); + if (nevents != neventsmax) { + printf("EVENT SHORTAGE? Got %i, expected %i.\n", nevents, neventsmax); std::cout << oscarfile << " vs. " << feof(oscarfile) << std::endl; - } - fclose(oscarfile); - oscarfile = NULL; + } + fclose(oscarfile); + oscarfile = NULL; // Calculate balance functions. // B_ab(y) = (N_a(0) - N_-a(0))(N_b(y) - N_-b(y)) / (N_b(y) + N_-b(y)) @@ -171,6 +171,6 @@ void CB3D::CalcBalance() { fprintf(anal_output, "%f %f\n", pseudorapidity, B_charge[i]); } - delete reslist; - fclose(anal_output); + delete reslist; + fclose(anal_output); } From c4fdd4a0beff411b421abc4621d909c80d188153 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 28 Oct 2014 16:16:22 -0400 Subject: [PATCH 09/19] improve output; more reasonable eta histogram cuts --- b3d2/src/anal_balance.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 6c55422..e926770 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -51,12 +51,13 @@ void CB3D::CalcBalance() { const double max_y = 1; const double max_dy = max_y - min_y; const double min_dy = 0; + const double min_deta = 0; + const double max_deta = 7; for (const auto& pair : balance_pairs) { histograms.emplace(make_pair(pair, Histogram(nbins, min_dy, max_dy))); } /// A pseudo-rapidity histogram for the charge-only balance function. - // TODO: use something other than min_dy/max_dy? - auto charge_histogram = Histogram(nbins, min_dy, max_dy); + auto charge_histogram = Histogram(nbins, min_deta, max_deta); /// Number for each species. /// (Needed for denominator of balance function.) @@ -71,7 +72,7 @@ void CB3D::CalcBalance() { do { KillAllParts(); const int nparts = ReadOSCAR(nevents + 1); - std::cout << nparts << std::endl; + std::cout << "Processing " << nparts << " particles..." << std::endl; if (nparts > 0) { nevents += 1; //printf("before, nparts=%d =? %d\n",nparts,int(FinalPartMap.size())); From a208a4534866d4e81513ecc1bf2e9158b41885d4 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Tue, 28 Oct 2014 18:14:57 -0400 Subject: [PATCH 10/19] refactor histogram to reduce code duplication --- b3d2/src/anal_balance.cc | 9 +++++---- b3d2/src/hist.h | 17 +++++++++++++---- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index e926770..7212002 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -81,7 +81,7 @@ void CB3D::CalcBalance() { // We are only interested in particles occuring in the balance function // and with y_min <= y <= y_max. - vector relevant_particles; // TODO: only remeber weight, charge, rapidity and momentum + vector relevant_particles; // TODO: only remember weight, charge, rapidity and momentum for (const auto& ppos : FinalPartMap) { const auto part = ppos.second; const int PID = part->resinfo->code; @@ -152,13 +152,14 @@ void CB3D::CalcBalance() { vector B(nbins); const int a = mypair.first; const int b = mypair.second; + const auto& histogram = histograms.find(mypair)->second; for (size_t i = 0; i < B.size(); i++) { - B[i] = (double) histograms.find(mypair)->second.histogram[i] / total_number[abs(b)]; + B[i] = (double) histogram.histogram[i] / total_number[abs(b)]; } fprintf(anal_output, "B(%i, %i)\n", a, b); fprintf(anal_output, "y B\n"); for (size_t i = 0; i < B.size(); i++) { - const double y = min_y + (max_y - min_y) / nbins * (i + 0.5); + const double y = histogram.get_value(i); fprintf(anal_output, "%f %f\n", y, B[i]); } fprintf(anal_output, "\n"); @@ -168,7 +169,7 @@ void CB3D::CalcBalance() { fprintf(anal_output, "B(+, -)\n"); for (size_t i = 0; i < B_charge.size(); i++) { B_charge[i] = (double) charge_histogram.histogram[i] / total_number_all; - const double pseudorapidity = min_y + (max_y - min_y) / nbins * (i + 0.5); + const double pseudorapidity = charge_histogram.get_value(i); fprintf(anal_output, "%f %f\n", pseudorapidity, B_charge[i]); } diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h index 1fedfb2..17332c2 100644 --- a/b3d2/src/hist.h +++ b/b3d2/src/hist.h @@ -17,9 +17,12 @@ class Histogram { /// Add value x to given histogram. void add(T x, int64_t count=1); /// Return the count of the bin corresponding to x. - int64_t get_count(T x); + int64_t get_count(T x) const; + /// Return the value corresponding to the given bin. + T get_value(size_t bin) const; private: - double get_bin(T x); + /// Get the bin correponding to a value. + double get_bin(T x) const; }; template @@ -31,7 +34,7 @@ Histogram::Histogram(size_t nbins, T x_min, T x_max) { } template -double Histogram::get_bin(T x) { +double Histogram::get_bin(T x) const { return (x - x_min) / ((x_max - x_min) / nbins); } @@ -47,10 +50,16 @@ void Histogram::add(T x, int64_t count) { } template -int64_t Histogram::get_count(T x) { +int64_t Histogram::get_count(T x) const { const double bin_float = get_bin(x); assert(bin_float >= 0); const size_t bin_int = (size_t) bin_float; assert(bin_int < nbins); return histogram[bin_int]; } + +template +T Histogram::get_value(size_t bin) const { + assert(bin < histogram.size()); + return x_min + (x_max - x_min) / nbins * (bin + 0.5); +} From d2013d8e06675c6fb4ca4bff401dfd6e36419daa Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 30 Oct 2014 18:23:14 -0400 Subject: [PATCH 11/19] fix balance functions They have been calculated incorrectly. The charge-only balance function still does not look right. --- b3d2/src/anal_balance.cc | 36 ++++++++++++++++++++++++++---------- b3d2/src/hist.h | 10 ++++++++++ 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 7212002..9bcabad 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -36,6 +36,15 @@ void CB3D::CalcBalance() { /// Balance functions to calculate. vector> balance_pairs = { {211, -211}, {321, -321} }; + /// We need to calculate the histograms for (a,-b), (-a,b) and (-a,-b) too. + const size_t old_size = balance_pairs.size(); + for (size_t i = 0; i < old_size; i++) { + const auto a = balance_pairs[i].first; + const auto b = balance_pairs[i].second; + balance_pairs.emplace_back(make_pair(-a, b)); + balance_pairs.emplace_back(make_pair(a, -b)); + balance_pairs.emplace_back(make_pair(-a, -b)); + } /// Relevant species (ignore anti particles). unordered_set balance_species; for (const auto& pair : balance_pairs) { @@ -44,7 +53,6 @@ void CB3D::CalcBalance() { } /// A rapidity histogram for each balance function. unordered_map, Histogram, boost::hash>> histograms; - histograms.reserve(balance_pairs.size()); // TODO: figure out what the following constants should be. const size_t nbins = 100; const double min_y = -1; @@ -87,6 +95,7 @@ void CB3D::CalcBalance() { const int PID = part->resinfo->code; const int charge = part->resinfo->charge; const int weight = part->weight; + assert(weight == 1); const double y = part->y; // Skip particles we do not care about. if (balance_species.count(abs(PID)) == 0 || y < min_y || y > max_y) @@ -96,6 +105,7 @@ void CB3D::CalcBalance() { total_number[abs(PID)] += weight * abs(charge); total_number_all += weight * abs(charge); } + std::cout << "Relevant particles: " << relevant_particles.size() << std::endl; // Iterate over all pairs to calculate dy and fill histograms. for (size_t i = 0; i < relevant_particles.size(); i++) for (size_t j = 0; j < i; j++) { @@ -129,12 +139,10 @@ void CB3D::CalcBalance() { if (p_j_mag == 0) p_j_mag = 1; const double dpseudorapidity = fabs( atanh(p_i.p[3] / p_i_mag) - atanh(p_j.p[3] / p_j_mag ) ); - assert(dpseudorapidity == dpseudorapidity); - // We want to calculate B(+-), so we want the negative charge. - const auto& p = p_i.resinfo->charge < 0 ? p_i : p_j; - const int charge = p.resinfo->charge; - const int weight = p.weight; - charge_histogram.add(dpseudorapidity, weight * charge); + //assert(dpseudorapidity == dpseudorapidity); + // B_+- ~ N_+ N_- - N_+ N_+ - N_- N_- + N_- N_+ + charge_histogram.add(dpseudorapidity, + -p_i.resinfo->charge * p_i.weight * p_j.resinfo->charge * p_j.weight); } } } @@ -149,17 +157,24 @@ void CB3D::CalcBalance() { // Calculate balance functions. // B_ab(y) = (N_a(0) - N_-a(0))(N_b(y) - N_-b(y)) / (N_b(y) + N_-b(y)) for (const auto& mypair : balance_pairs) { - vector B(nbins); + const auto& histogram = histograms.find(mypair)->second; + vector B(nbins, 0); const int a = mypair.first; const int b = mypair.second; - const auto& histogram = histograms.find(mypair)->second; + // We need to calculate N_a N_b - N_a N_-b - N_-a N_b + N_-a N_-b. + const vector> pairs = { {a, b}, {a, -b}, {-a, b}, {-a, -b} }; for (size_t i = 0; i < B.size(); i++) { - B[i] = (double) histogram.histogram[i] / total_number[abs(b)]; + B[i] += histogram.histogram[i]; + B[i] -= histograms.find(pairs[1])->second.histogram[i]; + B[i] -= histograms.find(pairs[2])->second.histogram[i]; + B[i] += histograms.find(pairs[3])->second.histogram[i]; + B[i] /= total_number[abs(b)]; } fprintf(anal_output, "B(%i, %i)\n", a, b); fprintf(anal_output, "y B\n"); for (size_t i = 0; i < B.size(); i++) { const double y = histogram.get_value(i); + //^ Does not matter which histogram we use. fprintf(anal_output, "%f %f\n", y, B[i]); } fprintf(anal_output, "\n"); @@ -172,6 +187,7 @@ void CB3D::CalcBalance() { const double pseudorapidity = charge_histogram.get_value(i); fprintf(anal_output, "%f %f\n", pseudorapidity, B_charge[i]); } + std::cout << "B(+, -): " << (double) charge_histogram.sum()/charge_histogram.nbins << " counts per bin." << std::endl; delete reslist; fclose(anal_output); diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h index 17332c2..4639dd5 100644 --- a/b3d2/src/hist.h +++ b/b3d2/src/hist.h @@ -20,6 +20,8 @@ class Histogram { int64_t get_count(T x) const; /// Return the value corresponding to the given bin. T get_value(size_t bin) const; + /// Sum all counts. + int64_t sum() const; private: /// Get the bin correponding to a value. double get_bin(T x) const; @@ -63,3 +65,11 @@ T Histogram::get_value(size_t bin) const { assert(bin < histogram.size()); return x_min + (x_max - x_min) / nbins * (bin + 0.5); } + +template +int64_t Histogram::sum() const { + int64_t sum = 0; + for (const auto s : histogram) + sum += s; + return sum; +} From 93a0ff9e61f22888deb0c53b6dbc506b9b2af660 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Fri, 31 Oct 2014 14:52:01 -0400 Subject: [PATCH 12/19] choose slightly small bin size for balance functions --- b3d2/src/anal_balance.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 9bcabad..262a6a6 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -54,7 +54,7 @@ void CB3D::CalcBalance() { /// A rapidity histogram for each balance function. unordered_map, Histogram, boost::hash>> histograms; // TODO: figure out what the following constants should be. - const size_t nbins = 100; + const size_t nbins = 80; const double min_y = -1; const double max_y = 1; const double max_dy = max_y - min_y; From c629c3303c43a6e4a41225ea1306a5e43e87c104 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 3 Nov 2014 17:11:42 -0500 Subject: [PATCH 13/19] make path to boost configurable --- b3d2/makefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/b3d2/makefile b/b3d2/makefile index a4ccb34..c6c2bea 100644 --- a/b3d2/makefile +++ b/b3d2/makefile @@ -14,6 +14,7 @@ GSLPATH=${MADAI_GSLPATH} #e.g. /opt/local HDF5_HOME=${MADAI_HDF5_HOME} #e.g. /opt/local +BOOST_INCLUDE=${MADAI_BOOST_HOME}/include ######################################################################### #HOPEFULLY, everything below need not be touched @@ -36,7 +37,7 @@ installdirs : clean : rm -f build/*.o include/*.h lib/*.a -INC=-I include -I ${GSLPATH}/include -I ${CORAL_INCLUDE} -I ${HDF5_HOME}/include +INC=-I include -I ${GSLPATH}/include -I ${CORAL_INCLUDE} -I ${BOOST_INCLUDE} -I ${HDF5_HOME}/include uninstall : rm -f ${INSTALLDIR}/lib/madai/libb3d.a ${INSTALLDIR}/include/madai/b3d.h ${INSTALLDIR}/include/madai/sampler.h ${INSTALLDIR}/include/madai/resonances.h ${INSTALLDIR}/include/madai/inelastic.h ${INSTALLDIR}/include/madai/joshconvert.h ${INSTALLDIR}/include/madai/hydrotob3d.h From 99b10679ba10024fa7d360b2f782cbe8ad6acd21 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 3 Nov 2014 18:08:03 -0500 Subject: [PATCH 14/19] don't normalice balance functions Instead output the numerator and denominator separately. This has the advantage that the results of several separate runs can be easily combined, because the division can be delayed until the last step, right before plotting the balance function. --- b3d2/src/anal_balance.cc | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 262a6a6..6e39da1 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -158,7 +158,7 @@ void CB3D::CalcBalance() { // B_ab(y) = (N_a(0) - N_-a(0))(N_b(y) - N_-b(y)) / (N_b(y) + N_-b(y)) for (const auto& mypair : balance_pairs) { const auto& histogram = histograms.find(mypair)->second; - vector B(nbins, 0); + vector B(nbins, 0); const int a = mypair.first; const int b = mypair.second; // We need to calculate N_a N_b - N_a N_-b - N_-a N_b + N_-a N_-b. @@ -168,25 +168,26 @@ void CB3D::CalcBalance() { B[i] -= histograms.find(pairs[1])->second.histogram[i]; B[i] -= histograms.find(pairs[2])->second.histogram[i]; B[i] += histograms.find(pairs[3])->second.histogram[i]; - B[i] /= total_number[abs(b)]; } fprintf(anal_output, "B(%i, %i)\n", a, b); fprintf(anal_output, "y B\n"); for (size_t i = 0; i < B.size(); i++) { const double y = histogram.get_value(i); //^ Does not matter which histogram we use. - fprintf(anal_output, "%f %f\n", y, B[i]); + fprintf(anal_output, "%f %lli\n", y, B[i]); } + fprintf(anal_output, "total number of %i\n%lli\n", b, total_number[abs(b)]); fprintf(anal_output, "\n"); } - vector B_charge(nbins); + vector B_charge(nbins); fprintf(anal_output, "B(+, -)\n"); for (size_t i = 0; i < B_charge.size(); i++) { - B_charge[i] = (double) charge_histogram.histogram[i] / total_number_all; + B_charge[i] = charge_histogram.histogram[i]; const double pseudorapidity = charge_histogram.get_value(i); - fprintf(anal_output, "%f %f\n", pseudorapidity, B_charge[i]); + fprintf(anal_output, "%f %lli\n", pseudorapidity, B_charge[i]); } + fprintf(anal_output, "total number\n%lli\n", total_number_all); std::cout << "B(+, -): " << (double) charge_histogram.sum()/charge_histogram.nbins << " counts per bin." << std::endl; delete reslist; From c8cdcb3ef7ba6cd01947d4cb881735794db7df2f Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 6 Nov 2014 11:25:06 -0500 Subject: [PATCH 15/19] seed RNG with hash of run name This ensures deterministic results while giving different results for parallel runs. --- b3d2/src/anal_balance.cc | 2 +- b3d2/src/b3d.cc | 3 ++- b3d2/src/joshconvert.cc | 3 ++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 6e39da1..3854d8e 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -14,7 +14,7 @@ int sign(double x) { /// Calculate the balance functions. void CB3D::CalcBalance() { - std::cout << "*** Calculating balance functions..." << std::endl; + std::cout << "*** Calculating balance functions for run " << run_name << "..." << std::endl; NSAMPLE = parameter::getI(parmap, "B3D_NSAMPLE", 1); const int neventsmax = parameter::getI(parmap, "B3D_NEVENTSMAX", 10); diff --git a/b3d2/src/b3d.cc b/b3d2/src/b3d.cc index 93a06de..9eece16 100644 --- a/b3d2/src/b3d.cc +++ b/b3d2/src/b3d.cc @@ -1,5 +1,6 @@ #ifndef __B3D_CC__ #define __B3D_CC__ +#include #include "b3d.h" CB3D::CB3D(){ @@ -74,7 +75,7 @@ CB3D::CB3D(string run_name_set){ FinalPartMap.clear(); ActionMap.clear(); DeadActionMap.clear(); - randy=new CRandom(-1234); + randy=new CRandom(std::hash()(run_name)); CAction::b3d=this; oscarfile=NULL; sampler=NULL; diff --git a/b3d2/src/joshconvert.cc b/b3d2/src/joshconvert.cc index 441bd69..59cf38a 100644 --- a/b3d2/src/joshconvert.cc +++ b/b3d2/src/joshconvert.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include "b3d.h" #define __EVALUATE_SURFACE_VELOCITIES__ @@ -26,7 +27,7 @@ void CjoshConvert::Init(string run_name_set){ NRINGSMAX=parameter::getD(parmap,"B3D_NRINGSMAX",100); NBOSE=parameter::getI(parmap,"B3D_NBOSE",1); boseweight.resize(NBOSE); - randy=new CRandom(-1234); + randy=new CRandom(std::hash()(run_name)); ring.resize(NRINGSMAX+1); reslist->CalcEoS(Tf,epsilonf,Pf,nhadronsf,cs2,densityf); reslist->CalcEoS(Tf,epsilonf,Pf,nhadronsf,densityf,boseweight); From fb6365f11942a2ae9997cae88c695422fe315ac8 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 6 Nov 2014 11:47:26 -0500 Subject: [PATCH 16/19] also calculate proton balance function --- b3d2/src/anal_balance.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index 3854d8e..a7b57a2 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -35,7 +35,8 @@ void CB3D::CalcBalance() { //^ We need to set this, so the header gets read. Otherwise `feof(oscarfile) == true`. /// Balance functions to calculate. - vector> balance_pairs = { {211, -211}, {321, -321} }; + // (pi+, pi-), (K+, K-), (p, pbar) + vector> balance_pairs = { {211, -211}, {321, -321}, {2212, -2212} }; /// We need to calculate the histograms for (a,-b), (-a,b) and (-a,-b) too. const size_t old_size = balance_pairs.size(); for (size_t i = 0; i < old_size; i++) { From bdf534f5d7ee35fbeef7c659021125a2f2056966 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 6 Nov 2014 15:46:15 -0500 Subject: [PATCH 17/19] calculate balance functions correctly The wrong signs were used, so B(pi+, pi-) was negative instead of positive. This is fixed here. Note that B(p, pbar) still seems to have the wrong sign for some reason. --- b3d2/src/anal_balance.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index a7b57a2..eac851b 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -116,18 +116,18 @@ void CB3D::CalcBalance() { const auto mypair = make_pair(p_i.resinfo->code, p_j.resinfo->code); const auto mypair_reversed = make_pair(mypair.second, mypair.first); - // For B(ab). + // For B(ab): Calculate N_a N_b. auto search = histograms.find(mypair); if (search != histograms.end()) { const int charge = p_j.resinfo->charge; const int weight = p_j.weight; - search->second.add(dy, weight * charge); + search->second.add(dy, weight * abs(charge)); } search = histograms.find(mypair_reversed); if (search != histograms.end()) { const int charge = p_i.resinfo->charge; const int weight = p_i.weight; - search->second.add(dy, weight * charge); + search->second.add(dy, weight * abs(charge)); } // For B(+-). From bf7858efa22a1e0c300f918d8c696b40b81d2148 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 10 Nov 2014 14:57:25 -0500 Subject: [PATCH 18/19] implement calculation of STAR acceptances This is a slight refactoring of Scott's and Garry's code. Use switch statements where appropriate (in one case this should improve performance) and return structs instead of modifying parameters. --- b3d2/makefile | 16 +++- b3d2/src/accept.cc | 208 +++++++++++++++++++++++++++++++++++++++++++++ b3d2/src/accept.h | 28 ++++++ 3 files changed, 249 insertions(+), 3 deletions(-) create mode 100644 b3d2/src/accept.cc create mode 100644 b3d2/src/accept.h diff --git a/b3d2/makefile b/b3d2/makefile index c6c2bea..0d5d7a7 100644 --- a/b3d2/makefile +++ b/b3d2/makefile @@ -52,7 +52,8 @@ include/hydrotob3d.h\ include/mutinfo.h\ include/seinfo.h\ include/regen.h\ -include/cell.h +include/cell.h\ +include/accept.h B3D_OBJFILES = build/b3d.o\ build/b3dsubs.o\ @@ -86,7 +87,8 @@ build/anal_v2.o\ build/consolidate.o\ build/anal_realitydiff.o\ build/joshconvert.o\ -build/hydrotob3d.o +build/hydrotob3d.o\ +build/accept.o ####################### @@ -134,6 +136,9 @@ include/cell.h : src/cell.h include/regen.h : src/regen.h cp src/regen.h include/ +include/accept.h : src/accept.h + cp src/accept.h include/ + ###################### lib/libb3d.a : ${B3D_HFILES} ${B3D_OBJFILES} @@ -230,7 +235,7 @@ build/anal_v2.o : src/anal_v2.cc ${B3D_HFILES} build/anal_realitydiff.o : src/anal_realitydiff.cc ${B3D_HFILES} ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/anal_realitydiff.o src/anal_realitydiff.cc -build/consolidate.o : src/consolidate.cc ${B3D_HFIES} +build/consolidate.o : src/consolidate.cc ${B3D_HFILES} ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/consolidate.o src/consolidate.cc build/joshconvert.o : src/joshconvert.cc ${B3D_HFILES} @@ -239,6 +244,8 @@ build/joshconvert.o : src/joshconvert.cc ${B3D_HFILES} build/hydrotob3d.o : src/hydrotob3d.cc ${B3D_HFILES} ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/hydrotob3d.o src/hydrotob3d.cc +build/accept.o : src/accept.cc ${B3D_HFILES} + ${CPP} -c ${OPT} ${INC} ${MADAI_COMPILERDEFS} -o build/accept.o src/accept.cc ############ ${INSTALLDIR}/lib/madai/libb3d.a : lib/libb3d.a @@ -277,6 +284,9 @@ ${INSTALLDIR}/include/madai/regen.h : include/regen.h ${INSTALLDIR}/include/madai/cell.h : include/cell.h cp -f include/cell.h ${INSTALLDIR}/include/madai/ +${INSTALLDIR}/include/madai/accept.h : include/accept.h + cp -f include/accept.h ${INSTALLDIR}/include/madai/ + ${INSTALLDIR}/progdata/madai/resinfo/resonances_pdg_new.dat : resinfo/resonances_pdg_new.dat cp -f resinfo/resonances_pdg_new.dat ${INSTALLDIR}/progdata/madai/resinfo/ diff --git a/b3d2/src/accept.cc b/b3d2/src/accept.cc new file mode 100644 index 0000000..07d4ac3 --- /dev/null +++ b/b3d2/src/accept.cc @@ -0,0 +1,208 @@ +#include +#include + +#include "accept.h" +#include "pow.h" + +/// Translate the general particle ID to the specific STAR PID. +/// Returns -1 if there is not corresponding STAR PID. +/// Exits if the PID is not in the STAR list. +int get_starpid(int pid) { + switch (pid) { + case 211: case -211: + return 1; + case 321: case -321: + return 2; + case -2212: + return 3; + case 2212: + return 4; + default: + if (abs(pid) != 2112 && abs(pid) != 311 && abs(pid) != 111 && abs(pid) != 22) { + printf("pid=%d isn't in STAR list\n", pid); + exit(1); + } + return -1; + } +} + +/// Routine to return acceptance and efficiency +/// Gary D. Westfall +/// January, 2014 +/// Uses results published by STAR PRC 79 034909 2009, page 034909-13 for functional forms +/// Uses embedding results from Hui Wang for Run 10, 200 GeV +/// Determined TOF matching efficiency looking at pion spectra +/// and double checked with kaons and protons +/// Input +/// pid=1 pion (either sign) +/// pid=2 kaon (with sign) +/// pid=3 anti-proton +/// pid=4 proton +/// cen=0 0-5% +/// cen=1 5-10% +/// cen=2 10-20% +/// cen=3 20-30% +/// cen=4 30-40% +/// cen=5 40-50% +/// cen=6 50-60% +/// cen=7 60-70% +/// cen=8 70-80% +/// basic acceptance cuts are +/// |eta| < 1.0 for TPC only +/// |eta| < 0.9 for TPC+TOF +/// full phi acceptance +/// pt dependence is calculated in detail +Acceptance star_acc_eff(int pid,double pt,double eta,double phi,int cen){ + int i; + float f1,f2,f3,f4,f5; + float P0,P1,P2,P3,P4,P5,P6; + float p,theta; + // Pion parameters + float P0_pion[9]={0.714,0.745,0.779,0.810,0.830,0.842,0.851,0.862,0.874}; + float P1_pion[9]={0.083,0.082,0.079,0.087,0.105,0.085,0.121,0.109,0.121}; + float P2_pion[9]={1.906,2.021,2.031,2.331,2.915,2.498,3.813,3.052,3.287}; + // Kaon parameters + float P0_kaon[9]={0.634,0.646,0.691,0.711,0.748,0.742,0.783,0.758,0.748}; + float P1_kaon[9]={0.226,0.220,0.215,0.205,0.204,0.194,0.202,0.198,0.218}; + float P2_kaon[9]={1.387,1.490,1.447,1.533,1.458,1.485,1.372,1.559,1.665}; + float P3_kaon[9]={0.021,0.026,0.021,0.025,0.019,0.027,0.016,0.029,0.042}; + // pbar parameters + float P0_pbar[9]={0.707,0.735,0.770,0.803,0.821,0.836,0.849,0.891,0.889}; + float P1_pbar[9]={0.191,0.185,0.196,0.194,0.191,0.200,0.199,0.156,0.197}; + float P2_pbar[9]={2.440,2.695,3.336,3.684,4.021,5.037,4.850,2.701,4.538}; + // p parameters + float P0_p[9]={0.720,0.750,0.785,0.818,0.837,0.852,0.865,0.885,0.921}; + float P1_p[9]={0.201,0.197,0.201,0.193,0.202,0.202,0.192,0.192,0.178}; + float P2_p[9]={3.247,3.467,4.124,4.267,5.414,5.621,5.080,4.950,3.957}; + + + // Check the input data + + if ( (pid < 1) || (pid > 4) + || (pt < 0.2) || (pt > 3.0) + || (fabs(eta) > 1.0) + || (cen < 0) || (cen > 8) ) { + return Acceptance(false, 0); + } + + // The basic parameters are OK, check the details + + bool accept=false; + double eff=0.0; + + switch (pid) { + case 1: // Pions + theta=2.0*atan(exp(-eta)); + p=pt/sin(theta); + if((pt > 0.2) && (p < 1.6)){ + accept=true; + P0=P0_pion[cen]; + P1=P1_pion[cen]; + P2=P2_pion[cen]; + f1=P1/pt; + eff=P0*exp(-pow(f1,P2)); + if(pt > 0.60){ + if(fabs(eta) < 0.9){ + eff=0.59*eff; + } else { + eff=0.0; + } + } + } + break; + + case 2: // Kaons + theta=2.0*atan(exp(-eta)); + p=pt/sin(theta); + if((pt > 0.20) && (p < 1.6)){ + accept=true; + P0=P0_kaon[cen]; + P1=P1_kaon[cen]; + P2=P2_kaon[cen]; + P3=P3_kaon[cen]; + f1=P1/pt; + eff=P0*exp(-pow(f1,P2))+P3*pt; + if(pt > 0.6){ + if(fabs(eta) < 0.9){ + eff=0.59*eff; + } else { + eff=0.0; + } + } + } + break; + + case 3: // pbar + theta=2.0*atan(exp(-eta)); + p=pt/sin(theta); + if((pt > 0.4) && (p < 3.0)){ + accept=true; + P0=P0_pbar[cen]; + P1=P1_pbar[cen]; + P2=P2_pbar[cen]; + f1=P1/pt; + eff=P0*exp(-pow(f1,P2)); + if(pt > 1.0){ + if(fabs(eta) < 0.9){ + eff=0.59*eff; + } else { + eff=0.0; + } + } + } + break; + + case 4: // p + theta=2.0*atan(exp(-eta)); + p=pt/sin(theta); + if((pt > 0.4) && (p < 3.0)){ + accept=true; + P0=P0_p[cen]; + P1=P1_p[cen]; + P2=P2_p[cen]; + f1=P1/pt; + eff=P0*exp(-pow(f1,P2)); + if(pt > 1.0){ + if(fabs(eta) < 0.9){ + eff=0.59*eff; + } else { + eff=0.0; + } + } + } + break; + default: + assert(0); + // unreachable + } + + return Acceptance(accept, eff); +} + +Acceptance calc_acceptance_STAR(const CPart& part, int centrality, double* dca) { + // Those are the cuts independent of the STAR acceptance. + // They might be more restrictive. + const double ETAMIN = -0.9; + const double ETAMAX = 0.9; + const double PTMIN = 200; + const double PTMAX = 1600; + const int CENTRALITY = 4; + + const int starpid = get_starpid(part.resinfo->code); + if (starpid == -1) + return Acceptance(false, 0); + const auto& p = part.p; + const double pt = sqrt(square(p[1]) + square(p[2])); + const double pmag = sqrt(square(pt) + square(p[3])); + const double eta = atanh(p[3]/pmag); + + if (ETAMIN < eta && eta < ETAMAX + && dca[0] < 2.0 + && PTMIN < pt && pt < PTMAX) { + const double phi = atan2(p[2], p[1]); + return star_acc_eff(starpid, 0.001*pt, eta, phi, CENTRALITY); + } + + return Acceptance(false, 0); +} + diff --git a/b3d2/src/accept.h b/b3d2/src/accept.h new file mode 100644 index 0000000..125c029 --- /dev/null +++ b/b3d2/src/accept.h @@ -0,0 +1,28 @@ +#pragma once + +#include "b3d.h" + +/// Reprecent whether a particle has been accepted, and with which efficiency. +struct Acceptance { + bool accept; + double efficiency; + + Acceptance(bool accept=false, double efficiency=0) { + this->accept = accept; + this->efficiency = efficiency; + } +}; + +/// Calculate the acceptance at STAR. +/// The possible values for `centrality` are the following: +/// 0 => 0- 5% +/// 1 => 5-10% +/// 2 => 10-20% +/// 3 => 20-30% +/// 4 => 30-40% +/// 5 => 40-50% +/// 6 => 50-60% +/// 7 => 60-70% +/// 8 => 70-80% +// TODO: figure out WTH dca is. +Acceptance calc_acceptance_STAR(const CPart& part, int centrality, double* dca); From 229e6a90da08f8f2f09206138c45fd64a2eb35fd Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Mon, 10 Nov 2014 17:08:21 -0500 Subject: [PATCH 19/19] use STAR cuts when calculating balance functions --- b3d2/src/accept.cc | 6 +++--- b3d2/src/accept.h | 2 +- b3d2/src/anal_balance.cc | 9 ++++++++- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/b3d2/src/accept.cc b/b3d2/src/accept.cc index 07d4ac3..c7c677f 100644 --- a/b3d2/src/accept.cc +++ b/b3d2/src/accept.cc @@ -179,14 +179,14 @@ Acceptance star_acc_eff(int pid,double pt,double eta,double phi,int cen){ return Acceptance(accept, eff); } -Acceptance calc_acceptance_STAR(const CPart& part, int centrality, double* dca) { +Acceptance calc_acceptance_STAR(const CPart& part, int centrality, const double* dca) { // Those are the cuts independent of the STAR acceptance. // They might be more restrictive. const double ETAMIN = -0.9; const double ETAMAX = 0.9; const double PTMIN = 200; const double PTMAX = 1600; - const int CENTRALITY = 4; + //const int CENTRALITY = 4; const int starpid = get_starpid(part.resinfo->code); if (starpid == -1) @@ -200,7 +200,7 @@ Acceptance calc_acceptance_STAR(const CPart& part, int centrality, double* dca) && dca[0] < 2.0 && PTMIN < pt && pt < PTMAX) { const double phi = atan2(p[2], p[1]); - return star_acc_eff(starpid, 0.001*pt, eta, phi, CENTRALITY); + return star_acc_eff(starpid, 0.001*pt, eta, phi, centrality); } return Acceptance(false, 0); diff --git a/b3d2/src/accept.h b/b3d2/src/accept.h index 125c029..6091603 100644 --- a/b3d2/src/accept.h +++ b/b3d2/src/accept.h @@ -25,4 +25,4 @@ struct Acceptance { /// 7 => 60-70% /// 8 => 70-80% // TODO: figure out WTH dca is. -Acceptance calc_acceptance_STAR(const CPart& part, int centrality, double* dca); +Acceptance calc_acceptance_STAR(const CPart& part, int centrality, const double* dca); diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc index eac851b..843d8bb 100644 --- a/b3d2/src/anal_balance.cc +++ b/b3d2/src/anal_balance.cc @@ -5,6 +5,7 @@ #include // needed for hashing pairs #include "b3d.h" +#include "accept.h" #include "hist.h" #include "pow.h" @@ -101,7 +102,13 @@ void CB3D::CalcBalance() { // Skip particles we do not care about. if (balance_species.count(abs(PID)) == 0 || y < min_y || y > max_y) continue; - // TODO: p_T cuts and efficiency + // Apply STAR cuts and determine efficiency. + // (Assuming 0-5% centrality.) + const double dca[1] = { 1.0 }; // WTH?! + auto acceptance = calc_acceptance_STAR(*part, 0, &dca[0]); + if (!acceptance.accept) + continue; + // TODO: efficiency relevant_particles.push_back(*part); total_number[abs(PID)] += weight * abs(charge); total_number_all += weight * abs(charge);