scatt_stats.C Source File

00001 
00002 // scatt_stats.C: Determine all cross-sections for three-body scattering.
00003 
00004 // Starlab application:  get_sigma3.
00005 
00006 #include "sigma3.h"
00007 
00008 // Global statistics data:
00009 
00010 #define NA 10
00011 #define NE 10
00012 
00013 static int n = 0;
00014 static int counter[NA][NE][N_RHO_ZONE_MAX];
00015 
00016 static real sigma[NA][NE];
00017 static real sigma_err_sq[NA][NE];
00018 
00019 local void initialize_stats()
00020 {
00021     for (int ia = 0; ia < NA; ia++)
00022         for (int je = 0; je < NE; je++)
00023             for (int kr = 0; kr < N_RHO_ZONE_MAX; kr++)
00024                 counter[ia][je][kr] = 0;
00025 }
00026 
00027 // accumulate_stats: Gather statistics supplied by get_sigma3. Customize here
00028 //                   to the specific application at hand.
00029 
00030 local void accumulate_stats(scatter_profile& prof,
00031                             initial_state3& init,
00032                             intermediate_state3& inter,
00033                             final_state3& final,
00034                             int rho_bin,
00035                             sigma_out& out)
00036 {
00037     if (final.descriptor == exchange_1 || final.descriptor == exchange_2) {
00038 
00039 //      cerr << "acc_stats: " << state_string(inter.descriptor)
00040 //           << " " << state_string(final.descriptor) << endl;
00041 
00042         if (prof.r3 > 0) {
00043 
00044             real ratio = final.sma / prof.r3;
00045 
00046             int ia = (int) (log(ratio) / log(2.0));
00047             if (ia >= NA) ia = NA - 1;
00048 
00049             int je = (int) (10 * final.ecc);
00050             if (je >= NE) je = NE - 1;
00051 
00052             n++;
00053             counter[ia][je][rho_bin]++;
00054 
00055 //          cerr << "           a/R = " << ratio
00056 //               << "  e = " << final.ecc << endl;
00057         }
00058     }
00059 }
00060 
00061 // normalize_counts: convert raw counts into cross-sections and errors
00062 
00063 local void normalize_counts(sigma_out& out)
00064 {
00065     for (int ia = 0; ia < NA; ia++)
00066         for (int je = 0; je < NE; je++) {
00067 
00068             sigma[ia][je] = 0;
00069             sigma_err_sq[ia][je] = 0;
00070 
00071             real rho_sq_min, rho_sq_max;
00072             int kr;
00073             for (kr = 0, rho_sq_min = 0, rho_sq_max = out.rho_sq_init;
00074                  kr <= out.i_max;
00075                  kr++, rho_sq_min = rho_sq_max, rho_sq_max *= RHO_SQ_FACTOR) {
00076 
00077                 real factor = (rho_sq_max - rho_sq_min) / out.trials_per_zone;
00078                 sigma[ia][je] += factor * counter[ia][je][kr];
00079                 sigma_err_sq[ia][je] += factor * factor * counter[ia][je][kr];
00080             }
00081         }
00082 } 
00083 
00084 local void print_stats(scatter_profile& prof, sigma_out& out)
00085 {
00086     normalize_counts(out);
00087     real v2 = prof.v_inf * prof.v_inf;
00088 
00089     int ia, je;
00090 
00091     cerr << "\nDifferential cross sections:\n";
00092 
00093     cerr << "\nRaw counts (column = log2(sma/r3), row = 10*ecc, total = "
00094          << n << "):\n\n";
00095     for (ia = 0; ia < NA; ia++) {
00096         for (je = 0; je < NE; je++) {
00097             int total = 0;
00098             for (int kr = 0; kr <= out.i_max; kr++)
00099                 total += counter[ia][je][kr];
00100             fprintf(stderr, " %7d", total);
00101         }
00102         cerr << endl;
00103     }
00104 
00105     cerr << "\nv^2 sigma / (pi a^2)\n\n";
00106     for (ia = 0; ia < NA; ia++) {
00107         for (je = 0; je < NE; je++)
00108             fprintf(stderr, " %7.3f", v2 * sigma[ia][je]);
00109         cerr << endl;
00110     }
00111 
00112     cerr << "\nv^2 (sigma error) / (pi a^2)\n\n";
00113     for (ia = 0; ia < NA; ia++) {
00114         for (je = 0; je < NE; je++)
00115             fprintf(stderr, " %7.3f", v2 * sqrt(sigma_err_sq[ia][je]));
00116         cerr << endl;
00117     }
00118 }
00119 
00120 main(int argc, char **argv)
00121     {
00122     int  debug  = 0;
00123     int  seed   = 0;
00124     int  n_rand = 0;
00125     real max_trial_density = 1.0;
00126 
00127     real cpu_time_check = 3600; // One check per CPU hour!
00128     real dt_snap = VERY_LARGE_NUMBER;
00129     real snap_cube_size = 10;
00130     int scatter_summary_level = 0;
00131 
00132     bool print_counts = FALSE;
00133 
00134     scatter_profile prof;
00135     make_standard_profile(prof);
00136 
00137     extern char *poptarg;
00138     int c;
00139     char* param_string = "A:c:C:d:D:e:m:M:N:pqQs:v:V:x:y:z:";
00140 
00141     while ((c = pgetopt(argc, argv, param_string)) != -1)
00142         switch(c)
00143             {
00144             case 'A': prof.eta = atof(poptarg);
00145                       break;
00146             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00147                       break;
00148             case 'C': snap_cube_size = atof(poptarg);
00149                       break;
00150             case 'd': max_trial_density = atof(poptarg);
00151                       break;
00152             case 'D': dt_snap = atof(poptarg);
00153                       scatter_summary_level = 2;  // Undo with later "-q/Q"
00154                       break;
00155             case 'e': prof.ecc = atof(poptarg);
00156                       prof.ecc_flag = 1;
00157                       break;
00158             case 'm': prof.m2 = atof(poptarg);
00159                       break;
00160             case 'M': prof.m3 = atof(poptarg);
00161                       break;
00162             case 'N': n_rand = atoi(poptarg);
00163                       break;
00164             case 'p': print_counts = 1 - print_counts;
00165                       break;
00166             case 'q': if (scatter_summary_level > 0)
00167                           scatter_summary_level = 0;
00168                       else
00169                           scatter_summary_level = 1;
00170                       break;
00171             case 'Q': if (scatter_summary_level > 0)
00172                           scatter_summary_level = 0;
00173                       else
00174                           scatter_summary_level = 2;
00175                       break;
00176             case 's': seed = atoi(poptarg);
00177                       break;
00178             case 'v': prof.v_inf = atof(poptarg);
00179                       break;
00180             case 'V': debug = atoi(poptarg);
00181                       break;
00182             case 'x': prof.r1 = atof(poptarg);
00183                       break;
00184             case 'y': prof.r2 = atof(poptarg);
00185                       break;
00186             case 'z': prof.r3 = atof(poptarg);
00187                       break;
00188             case '?': params_to_usage(cerr, argv[0], param_string);
00189                       exit(1);
00190             }
00191 
00192     // Debugging:    |debug| = 0 ==> no debugging output
00193     //               |debug| = 1 ==> output only at end (default)
00194     //               |debug| = 2 ==> output at end of each top-level iteration
00195     //               |debug| = 3 ==> output after each sub-trial
00196     //
00197     //                debug  < 0 ==> simple statistics also (default: off)
00198 
00199     cpu_init();
00200 
00201     sigma_out out;
00202     int first_seed = srandinter(seed, n_rand);
00203 
00204     cerr << "Random seed = " << first_seed << endl;
00205     print_profile(cerr, prof);
00206 
00207     initialize_stats();
00208 
00209     get_sigma3(max_trial_density, prof, out,
00210                debug, cpu_time_check,
00211                dt_snap, snap_cube_size,
00212                scatter_summary_level,
00213                accumulate_stats);
00214 
00215     print_sigma3(out, prof.v_inf * prof.v_inf);
00216     if (print_counts) print_all_sigma3_counts(out, cerr);
00217 
00218     print_stats(prof, out);
00219 }