list-particles-2.c

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>

#include <math.h>

#include "corsika-tape.h"
#include "corsika-particle.h"

#ifdef QUIET
#  define Q 0
#else
#  define Q 1
#endif

typedef struct {
  RWORD r, m, e;
} particle_extra;

int main(int argc, char **argv) {
  corsika_file   *tape;
  sub_block      *sb;
  particle_data  *pd;
  particle_extra pe;
  int            i, j;
  RWORD          desc, A, Z;
  double         w = 0;
  int            p = 0;

  struct {
    long n;
    double w, rmax;
  } rb[] = {
    {0, 0, 10.e2},
    {0, 0, 31.6e2},
    {0, 0, 100.e2},
    {0, 0, 316.e2},
    {0, 0, 1000.e2},
    {0, 0, 3160.e2},
    {0, 0, 10000.e2},
    {0, 0, 31600.e2},
    {0, 0, 0.e0}
  };
  RWORD rmin, rmax;

  if ((tape = corsika_fopen(argv[1], "r")) == NULL) {
    fprintf(stderr, "Cannot open file %s. Aborting...\n", argv[1]);
    exit(1);
  }

  while ((sb = corsika_get(tape)) != NULL) {
    if (is_control(sb->rh.id)) {
      Q && printf("%4.4s\n", sb->rh.id);
    } else {
      for(i=0; i<PARTICLES_IN_SUB_BLOCK; i++) {
        pd = sb->pb.particle + i;
        desc = pd->description;
        Q && printf("%8g: ",desc);
        if (is_particle_r(desc)) {
          Q && printf("%-10s", 
                      particle[(int)desc/1000].name);
        } else if(is_nucleus_r(desc)) {
          A = (int)desc / 100000;
          Z = (int)desc / 1000 - 100*A;
          Q && printf("N, Z=%g, A=%g", Z, A);
        } else {
          Q && printf("OTHER: %g\n", desc);
        }

        if (desc != 0) {
          pe.r = sqrt(pd->x * pd->x + pd->y* pd->y);
          for (j=0; ; j++) {
            if (pe.r < rb[j].rmax || rb[j].rmax == 0.) {
              rb[j].n++;
#ifdef THINNING
              rb[j].w += pd->weight;
#else
              rb[j].w += 1;
#endif
              break;
            }
          }
          Q && printf(", Distance %7.2f (m)", pe.r / 100.);
#ifdef THINNING
          Q && printf(", Weight: %g", pd->weight);
#endif
          Q && printf("\n");

          p++;
#ifdef THINNING
          w += pd->weight;
#endif
        }
      }
    }
  }
  printf("Total particles: %i", p);
#ifdef THINNING
  printf(", weighted: %g, average weight %g", w, w/p);
#endif
  printf(".\n\n");

  w = 0;
  rmin = 0;
  for (j=0; ; j++) {
    rmax = rb[j].rmax;

    printf("Bin %2i (r < ", j);
    if (rmax != 0) printf("%8.0f", rmax / 100.);
    else                 printf("infinity");
    printf(" m): %8e particles, %i hits,\n", rb[j].w, rb[j].n);
    printf("                         avg. weight %g", 
           rb[j].n ? rb[j].w / rb[j].n : 0.0) ;

    if (rmax != 0) 
      printf(", density %8.3e m^(-2)", 
             rb[j].w / (rmax*rmax - rmin*rmin) / 3.14159e-4);

    printf(".\n");

    w += rb[j].w;

    if (rmax == 0) break;
    rmin = rmax;

  }
  printf("\nSum of all bins: %g\n", w);
}