KaliVeda  1.13/01
Heavy-Ion Analysis Toolkit
KVPartitionFunction.cpp
Go to the documentation of this file.
1 //Created by KVClassFactory on Thu Sep 30 16:48:24 2010
2 //Author: John Frankland,,,,
3 
4 #include "KVPartitionFunction.h"
5 #include "TMath.h"
6 
8 
9 
10 
11 
12 
16  : fTable(10000000, 5)
17 {
18  // Default constructor
19 }
20 
21 
22 
25 
27 {
28  // Destructor
29 }
30 
31 
32 
47 
49 {
50  // recursive calculation of number of partitions of A nucleons
51  // into M fragments of size 1, 2, ..., A
52  // from Eq. (2.6) of Bondorf et al., Nucl. Phys. A443, 321 (1985):
53  /* "We now consider the division of the partition space into the subsets characterized
54  by having a common multiplicity. Let P(A,, M) be the number of partitions in one
55  of such subsets. To compute P(A,, M) it is sufficient to observe that it should be
56  equal to the number of partitions having at least one fragment composed of just
57  one nucleon (which number is immediately shown to equal P(A0 - 1, M - 1)) plus
58  those not having any unitary fragment. The number of this last portion can be
59  computed by observing that, if we take one nucleon from each of its fragments one
60  has the set of all partitions of A0 - M nucleons into M fragments. We thus arrive at
61  P(A0,M) = P(A0-1,M-l) + P(A0-M, M) Eq. (2.5)"
62  */
63 // The expression below follows from recursive use of this relation.
64 
65  if (!A) return 0;
66  if (A == M) return 1;
67  if (M == 1) return 1;
68  Double_t p = 0;
69  int klim = A / M - 1;
70  for (int k = 0; k <= klim; k++) {
71  p += PartFunc(A - k * M - 1, M - 1);
72  }
73  return p;
74 }
75 
76 
77 
82 
84 {
85  // calculate total number of partitions of A nucleons into
86  // 1, 2, ..., A fragments.
87  // from Eq. (2.7) of Bondorf et al., Nucl. Phys. A443, 321 (1985)
88 
89  Double_t p = 0.;
90  for (int m = 1; m <= A; m++) p += PartFunc(A, m);
91  return p;
92 }
93 
94 
95 
96 
98 
100 {
101  if (A > 0 && Z >= 0 && M > 0 && B >= 0) {
102  Double_t snc = get_value(A, Z, M, B);
103  if (snc < 0) {
104  snc = calc_sneppen_Nclass(A, Z, M, B);
105  store_value(snc, A, Z, M, B);
107  }
108  return snc;
109  }
110  return 0;
111 }
112 
113 
114 
116 
118 {
119  if (A > 0 && Z >= 0 && M > 0) {
120  if (Z > A - Z) Z = A - Z; // symmetry
121  Double_t snc = get_value(A, Z, M);
122  if (snc < 0) {
123  snc = calc_sneppen_Np(A, Z, M);
124  store_value(snc, A, Z, M);
126  }
127  return snc;
128  }
129  return 0;
130 }
131 
132 
133 
136 
138 {
139  // calculates Eqs. 5a-d
140 
141  if (M <= 0) return 0;
142  if (A <= 0) return 0;
143  if (Z < 0) return 0;
144  if (B < 0) return 0;
145  if (B > A - Z) return 0;
146  if (Z == 0 && B == 0) return 0;
147  if ((M - 1) > (A - B)) return 0;
148  if (M == A) {
149  if (Z == A && B > 0) return 0;
150  if (B != 1) return 0;
151  }
152  if (Z == 0) {
153  if (A - B * M > 0) {
154  Double_t sn1 = sneppen_Np(A - (B - 1) * M, 0, M);
155  Double_t sn2 = sneppen_Np(A - B * M, 0, M);
156  return sn1 - sn2;
157  }
158  return 0;
159  }
160  if (B == 0) {
161  if (A - M >= M && Z - M >= 0) {
162  Double_t snp, snc;
163  snp = sneppen_Np(A - M, Z - M, M);
164  snc = sneppen_Nclass(A - 1, Z - 1, M - 1, 0);
165  return snp + snc;
166  }
167  return 0;
168  }
169  Double_t ncl = 0;
170  if (A - B > 0 && M - 1 > 0 && B > 1) {
171  for (int i = 1; i <= B - 1; i++) {
172  ncl += sneppen_Nclass(A - B, Z, M - 1, i);
173  }
174  }
175  ncl = sneppen_Np(A - B, Z, M - 1) - ncl;
176  return ncl;
177 }
178 
179 
180 
183 
185 {
186  // calculates Eq. 4
187 
188  if (A <= 0) return 0;
189  if (M <= 0) return 0;
190  if (Z < 0) return 0;
191  if (Z == 0) {
192  // one-component system
193  return PartFunc(A, M);
194  }
195  if (M > A) return 0;
196  if (M == 1) return 1;
197  if (M == A) return 1;
198  Double_t Np = 0;
199  int N = A - Z;
200  for (int b = 0; b <= N; b++) Np += sneppen_Nclass(A, Z, M, b);
201  return Np;
202 }
203 
204 
205 
210 
212 {
213  // Returns the total number of partitions of Z protons and (A-Z) neutrons
214  // into M fragments, using the method given by K. Sneppen
215  // in Nucl. Phys. A470, 213 (1987), Eqs. (4)-(6).
216 
217  Double_t p = sneppen_Np(A, Z, M);
218  /*Info("PartFunc(A,Z,M)", "p=%f array use=%f max value=%f\n",
219  p, NvalsNcl/pow(SNEPPENMAXTAB,4),GetMaxValueNclass());*/
220  return p;
221 }
222 
223 
224 
229 
231 {
232  // Returns the total number of partitions of Z protons and (A-Z) neutrons
233  // summed over all multiplicities, using the method given by K. Sneppen
234  // in Nucl. Phys. A470, 213 (1987), Eqs. (4)-(6).
235 
236  Double_t p = 0;
237  for (int m = 1; m <= A; m++) p += sneppen_Np(A, Z, m);
238  /*Info("PartSum(A,Z)", "p=%f array use=%f max value=%f\n",
239  p, NvalsNcl/pow(SNEPPENMAXTAB,4),GetMaxValueNclass());*/
240  return p;
241 }
242 
243 
244 
253 
255 {
256  // Calculate the mean number of clusters of size A when a system of size A0
257  // fragments in all possible ways with equal probability
258  // Using Eq. (3) of K. Sneppen Nucl. Phys. A470, 213 (1987)
259  //
260  // Correction for A=A0:
261  // there is always exactly 1 partition of M=1 composed of one fragment with A=A0
262  // the mean multiplicity of A=A0 is therefore 1/PartSum(A0)
263 
264  if (A == A0) return 1. / PartSum(A0);
265  Double_t NA = 0;
266  Int_t imax = (Int_t)(1.*A0 / (1.*A));
267  for (int i = 1; i <= imax; i++) {
268  Int_t A1 = A0 - i * A;
269  if (A1 > 0) NA += PartSum(A1);
270  }
271  NA /= PartSum(A0);
272  return NA;
273 }
274 
275 
276 
294 
296 {
297  // Calculate the mean number of clusters of size A when a system of size A0
298  // breaks up into M fragments, each partition having equal probability
299  // Using Eqs. (2b) & (3) of K. Sneppen Nucl. Phys. A470, 213 (1987)
300  //
301  // Correction for M=1:
302  // if A==A0, mean multiplicity is 1; otherwise 0.
303  // Correction for M=2 and even values of A0:
304  // when an odd-A0 splits into M=2, there are PartFunc(A0,M) partitions
305  // in which each A between 1 and A0-1 occurs once only, therefore
306  // the mean multiplicity for any A is 1/PartFunc(A0,M)
307  // when an even-A0 splits into M=2, one of the PartFunc(A0,M) partitions
308  // is the symmetric split (A0/2, A0/2). Therefore the mean multiplicity
309  // of A=A0/2 is twice that of the other A, i.e. 2/PartFunc(A0,M)
310  // Correction for M=A0:
311  // in this case only one partition exists, made of M=A0 monomers A=1
312  // therefore mean multiplicity of A=1 is M, for all other A it is 0.
313 
314  if (M == 1) {
315  if (A == A0) return 1.;
316  else return 0.;
317  }
318  if (M == 2) {
319  if (!(A0 % 2) && A == A0 / 2) return 2. / PartFunc(A0, M);
320  else if (A >= 1 && A < A0) return 1. / PartFunc(A0, M);
321  else return 0.;
322  }
323  if (M == A0) {
324  if (A == 1) return M;
325  else return 0.;
326  }
327  Double_t NA = 0;
328  Int_t imax = (Int_t)(1.*A0 / (1.*A));
329  imax = TMath::Min(imax, M - 1);
330  for (int i = 1; i <= imax; i++) {
331  Int_t A1 = A0 - i * A;
332  if (A1 > 0)NA += PartFunc(A1, M - i);
333  }
334  NA /= PartFunc(A0, M);
335  return NA;
336 }
337 
338 /*
339  Double_t MeanNA(int A0, int Z0, int A);
340  Double_t MeanNZ(int A0, int Z0, int Z);*/
341 
342 
347 
348 Double_t KVPartitionFunction::MeanNAZ(int A0, int Z0, int A, int Z)
349 {
350  // Calculate the mean number of clusters of mass A and charge Z
351  // when a two-component system (A0,Z0) fragments in all possible ways with equal probability
352  // Using Eq. (7) of K. Sneppen Nucl. Phys. A470, 213 (1987)
353 
354  Double_t NAZ = 0;
355  Int_t imax = (Int_t)(1.*A0 / (1.*A));
356  for (int i = 1; i <= imax; i++) {
357  Int_t A1 = A0 - i * A;
358  Int_t Z1 = Z0 - i * Z;
359  if (A1 > 0) NAZ += PartSum(A1, Z1);
360  }
361  NAZ /= PartSum(A0, Z0);
362  return NAZ;
363 }
364 
365 
366 
371 
372 Double_t KVPartitionFunction::MeanNA(int A0, int Z0, int A)
373 {
374  // Calculate the mean number of clusters of mass A
375  // when a two-component system (A0,Z0) fragments in all possible ways with equal probability
376  // This is just the sum of MeanNAZ(A0,Z0,A,Z) with 0<=Z<= min(A,Z0)
377 
378  Double_t NA = 0;
379  for (Int_t Z = 0; Z <= TMath::Min(A, Z0); Z++) NA += MeanNAZ(A0, Z0, A, Z);
380  return NA;
381 }
382 
383 
384 
389 
390 Double_t KVPartitionFunction::MeanNZ(int A0, int Z0, int Z)
391 {
392  // Calculate the mean number of clusters of charge Z
393  // when a two-component system (A0,Z0) fragments in all possible ways with equal probability
394  // This is just the sum of MeanNAZ(A0,Z0,A,Z) with Z<=A<=(Z+A0-Z0)
395 
396  Double_t NZ = 0;
397  for (Int_t A = Z; A <= (Z + A0 - Z0); A++) NZ += MeanNAZ(A0, Z0, A, Z);
398  return NZ;
399 }
400 
401 
402 
406 
408 {
409  // Calculate the mean multiplicity of all partitions of a single-component
410  // system of mass A0.
411  Double_t mult = 0.;
412  Double_t part_sum = 0;
413  Double_t part_func = 0;
414  for (int i = 1; i <= A0; i++) {
415  part_func = PartFunc(A0, i);
416  mult += i * part_func;
417  part_sum += part_func;
418  }
419  mult /= part_sum;
420  return mult;
421 }
422 
423 
424 
428 
430 {
431  // Calculate the mean size of clusters in all partitions of a single-component
432  // system of mass A0.
433  Double_t mult, moy = 0.;
434  Double_t sum = 0;
435  for (int i = 1; i <= A0; i++) {
436  mult = MeanNA(A0, i);
437  moy += i * mult;
438  sum += mult;
439  }
440  moy /= sum;
441  return moy;
442 }
443 
444 
int Int_t
KVIonRangeTableMaterial * M
ClassImp(KVPartitionList) void KVPartitionList
Initialisation.
#define b(i)
double Double_t
#define N
Calculates number of partitions of (A,Z,M)
Double_t get_value(Int_t A, Int_t Z, Int_t M, Int_t B=-1)
Double_t sneppen_Np(int A, int Z, int M)
Double_t PartFunc(int A, int M)
Double_t MeanNAZ(int A0, int Z0, int A, int Z)
Double_t MeanNA(int A0, int A)
Double_t calc_sneppen_Np(int A, int Z, int M)
calculates Eq. 4
Double_t calc_sneppen_Nclass(int A, int Z, int M, int B)
calculates Eqs. 5a-d
Double_t sneppen_Nclass(int A, int Z, int M, int B)
Double_t MeanNZ(int A0, int Z0, int Z)
void store_value(Double_t val, Int_t A, Int_t Z, Int_t M, Int_t B=-1)
virtual ~KVPartitionFunction()
Destructor.
Double_t MeanNA_M(int A0, int A, int M)
const long double m
Definition: KVUnits.h:70
Double_t Min(Double_t a, Double_t b)
Double_t Max(Double_t a, Double_t b)