diff --git a/b3d2/makefile b/b3d2/makefile index 5d84910..0d5d7a7 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 @@ -51,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\ @@ -79,12 +81,14 @@ 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\ build/anal_realitydiff.o\ build/joshconvert.o\ -build/hydrotob3d.o +build/hydrotob3d.o\ +build/accept.o ####################### @@ -132,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} @@ -215,6 +222,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 @@ -225,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} @@ -234,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 @@ -272,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..c7c677f --- /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, 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 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..6091603 --- /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, const double* dca); diff --git a/b3d2/src/anal_balance.cc b/b3d2/src/anal_balance.cc new file mode 100644 index 0000000..843d8bb --- /dev/null +++ b/b3d2/src/anal_balance.cc @@ -0,0 +1,203 @@ +#include +#include +#include + +#include // needed for hashing pairs + +#include "b3d.h" +#include "accept.h" +#include "hist.h" +#include "pow.h" + +int sign(double x) { + return (x > 0) - (x < 0); +} + +/// Calculate the balance functions. +void CB3D::CalcBalance() { + 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); + + 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; + oscarfile = NULL; + //^ We need to set this, so the header gets read. Otherwise `feof(oscarfile) == true`. + + /// Balance functions to calculate. + // (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++) { + 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) { + balance_species.insert(abs(pair.first)); + balance_species.insert(abs(pair.second)); + } + /// 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 = 80; + const double min_y = -1; + 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. + auto charge_histogram = Histogram(nbins, min_deta, max_deta); + + /// 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 { + KillAllParts(); + const int nparts = ReadOSCAR(nevents + 1); + std::cout << "Processing " << nparts << " particles..." << 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())); + + // We are only interested in particles occuring in the balance function + // and with y_min <= y <= y_max. + 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; + 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) + continue; + // 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); + } + 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++) { + 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); + + // 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 * 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 * abs(charge)); + } + + // For B(+-). + if (sign(p_i.resinfo->charge) == -sign(p_j.resinfo->charge)) { + 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); + // 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); + } + } + } + } 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; + + // 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) { + const auto& histogram = histograms.find(mypair)->second; + 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. + const vector> pairs = { {a, b}, {a, -b}, {-a, b}, {-a, -b} }; + for (size_t i = 0; i < B.size(); i++) { + 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]; + } + 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 %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); + fprintf(anal_output, "B(+, -)\n"); + for (size_t i = 0; i < B_charge.size(); i++) { + B_charge[i] = charge_histogram.histogram[i]; + const double pseudorapidity = charge_histogram.get_value(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; + fclose(anal_output); +} 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/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 diff --git a/b3d2/src/hist.h b/b3d2/src/hist.h new file mode 100644 index 0000000..4639dd5 --- /dev/null +++ b/b3d2/src/hist.h @@ -0,0 +1,75 @@ +#pragma once + +#include +#include +#include + +template +class Histogram { + public: + T x_min; + T x_max; + size_t nbins; + std::vector histogram; + + /// Initialize histogram. + Histogram(size_t nbins, T x_min, T x_max); + /// 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) 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; +}; + +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) const { + 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 { + 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); +} + +template +int64_t Histogram::sum() const { + int64_t sum = 0; + for (const auto s : histogram) + sum += s; + return sum; +} 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); 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); 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); +} diff --git a/b3d2/src/test_hist.cc b/b3d2/src/test_hist.cc new file mode 100644 index 0000000..81285a7 --- /dev/null +++ b/b3d2/src/test_hist.cc @@ -0,0 +1,39 @@ +#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; + std::cout << "\n" << h.get_count(0) << std::endl; +}