44double probability(
const vector<double>& p,
const vector<int>& i)
46 if (p.size() < i.size())
47 throw runtime_error(
"p not big enough");
52 const int n = accumulate(i.begin(), i.end(), 0);
56 vector<int> p_count = i, C_count = i;
61 for (
int count=0; count<3*n; )
63 if (n_count && result<=1)
70 while ((result>=1 || n_count==0) && accumulate(C_count.begin(), C_count.end(), 0))
72 vector<int>::iterator it = find_if(C_count.begin(), C_count.end(),
nonzero);
73 if (it == C_count.end())
throw runtime_error(
"blech");
79 while ((result>=1 || n_count==0) && accumulate(p_count.begin(), p_count.end(), 0))
81 vector<int>::iterator it = find_if(p_count.begin(), p_count.end(),
nonzero);
82 if (it == p_count.end())
throw runtime_error(
"blech2");
83 size_t index = it - p_count.begin();
112 if (
os_) *
os_ <<
"C distribution: " << md << endl;
115 for (MassDistribution::const_iterator it=md.begin(); it!=md.end(); ++it)
116 p.push_back(it->abundance);
118 vector<int> neutron0; neutron0.push_back(100);
119 vector<int> neutron1; neutron1.push_back(99); neutron1.push_back(1);
120 vector<int> neutron2; neutron2.push_back(98); neutron2.push_back(2);
121 vector<int> neutron3; neutron3.push_back(97); neutron3.push_back(3);
122 vector<int> neutron4; neutron4.push_back(96); neutron4.push_back(4);
124 vector<double> abundances;
133 *
os_ <<
"C100 abundances (manually calculated):\n";
134 copy(abundances.begin(), abundances.end(), ostream_iterator<double>(*
os_,
"\n"));
139 if (
os_) *
os_ <<
"C100 distribution (from IsotopeCalculator):\n"
143 for (
unsigned int i=0; i<abundances.size(); i++)