RALICE/icepack/iceconvert/uwi.c

   1 /*
   2  * Read/write files in the format used at UWI Madison (see J.Jacobsen's PhD
   3  * thesis)
   4  */
   5
   6 char *uwi_rdmc_cvsid =
   7 "$Header: /net/local/cvsroot/siegmund/rdmc/uwi.c,v 1.33 1999/12/23 12:23:15 wiebusch Exp $";
   8
   9 #include <stdio.h>
  10 #define _USE_MATH_DEFINES
  11 #include <math.h>
  12 #include <string.h>
  13 #include <stdlib.h>
  14 #include <ctype.h>
  15
  16 ///////////////#include <unistd.h>
  17
  18
  19 #include <fcntl.h>
  20 #include <sys/stat.h>
  21
  22 #if defined(OSF1) || defined(AIX) || defined (SunOS)
  23 #include <errno.h>
  24 #endif
  25
  26 #include "rdmc.h"
  27
  28 #include "rdmc_local.h"
  29
  30 #include "uwi.h"
  31
  32 #ifndef DEBUG
  33 #define DEBUG 0                       /* DEBUG 1 writes some info to stderr */
  34 #endif
  35
  36 #define ERRLINE (-__LINE__)
  37
  38 #define PID180 (M_PI/180.)
  39
  40 #ifdef UWI_ASCII_F
  41
  42 int uwi_rd_EH(char *s, mevt *ev);   /* read the several lines for UWI files */
  43 int uwi_rd_DH(char *s, mevt *ev);
  44 int uwi_rd_MU(char *s, mtrack **tr, int *ntrack);
  45 int uwi_rd_SH(char *s, mtrack **tr, int *ntrack);
  46 int uwi_rd_HT(char *s, mhit **h, int *nhits);
  47 int uwi_rd_FT(char *s, mtrack **fitr, mevt_special_t **fitres, int *nfit);
  48 int uwi_rd_SP(char *s, mtrack **fitr, mevt_special_t **fitres, int *nfit);
  49
  50
  51
  52 /****************************************************************************/
  53 /* The function rarr_uwi() reads the geometry of a UWI like file            */
  54 /****************************************************************************/
  55 int rdmc_rarr_uwi(char *geofile, array *ar)
  56 {
  57   FILE *fp;
  58   char s[RDMC_MAXLINE];
  59   char omfile[RDMC_MAXLINE];
  60   int form;
  61   int ich = 0;
  62   int ichstr = 0;
  63   int istr;
  64   float x, y, z;
  65   int nchstr;          /* temp. x, y, z coor, number of channels per string */
  66
  67   rdmc_init_array(ar);                               /* reset the array */
  68
  69   if (geofile == NULL)      /* use the default file (environement variable) */
  70     geofile = getenv(DEFAULT_UWI_GEO_FILE_ENV);
  71   if (geofile == NULL)                              /* use the default file */
  72     geofile = DEFAULT_UWI_GEO_FILE;                            /* (builtin) */
  73
  74   fp = fopen(geofile,"r");
  75
  76   if (fp == NULL) return ERRLINE; /* no geometry file found */
  77
  78   do {
  79     if (fgets(s,RDMC_MAXLINE,fp) == NULL)
  80       return EOF;
  81     form = sscanf(s,"%f",&(ar->depth));          /* 1. line: detector depth */
  82   } while (form < 1);
  83   do {
  84     if (fgets(s,RDMC_MAXLINE,fp) == NULL)
  85       return EOF;
  86     form = sscanf(s,"%i",&(ar->nstr));        /* 2. line: Number of strings */
  87   } while (form < 1);
  88   for (istr = 0; istr < ar->nstr; istr++) {             /* read all strings */
  89     do {
  90       if (fgets(s,RDMC_MAXLINE,fp) == NULL)
  91         return EOF;
  92       form = sscanf(s,"%f %f",&x,&y);         /* 3. line: x,y of 1st string */
  93     } while (form < 2);
  94     do {
  95       if (fgets(s,RDMC_MAXLINE,fp) == NULL)
  96         return EOF;
  97       form = sscanf(s,"%i",&nchstr);          /* 4. line: nch of 1st string */
  98     } while (form < 1);
  99     ar->nch += nchstr;        /* add the number of channels to the total nr */
 100     for (ichstr = 0; ichstr < nchstr; ich++, ichstr++) {  /* read ch in str */
 101       do {
 102         if (fgets(s,RDMC_MAXLINE,fp) == NULL)
 103           return EOF;
 104         form = sscanf(s,"%f",&z);              /* 5. line: z of 1st channel */
 105       } while (form < 1);
 106       do {
 107         if (fgets(s,RDMC_MAXLINE,fp) == NULL)
 108           return EOF;
 109         form = sscanf(s,"%s",omfile);         /* 6. line: om file of 1st ch */
 110       } while (form < 1);
 111       ar->x[ich] = x;
 112       ar->y[ich] = y;
 113       ar->z[ich] = z;
 114       ar->str[ich] = istr+1;
 115       ar->type[ich] = STD_AMANDA_B_OM;
 116       ar->clust[ich] = ichstr;
 117       if (strcmp(omfile,"up.pmt") == 0){
 118         ar->costh[ich] = 1.0;
 119         ar->sensit[ich] = 1.0;
 120       }else if (strcmp(omfile,"down.pmt") == 0){
 121         ar->costh[ich] = -1.0;
 122         ar->sensit[ich] = 1.0;
 123       }else if (strcmp(omfile,"dead.pmt") == 0){
 124         ar->costh[ich] = 0.0;
 125         ar->sensit[ich] = 0.0;
 126       } else{
 127         ar->costh[ich] = 0.0;
 128         ar->sensit[ich] = 0.0;
 129       }
 130   } /* for ich */
 131   } /* for istr */
 132
 133   ar->id = AMANDA_B_4;                             /* default: amanda array */
 134   ar->longitude = -90.0;                  /* geographic longitude of AMANDA */
 135   ar->lattitude = 0.0;                              /* geographic lattitude */
 136
 137   {
 138     array_hdef_t def_trig;
 139
 140     rdmc_init_array_hdef(&def_trig);
 141     def_trig.id = 1;
 142     strcpy(def_trig.tag, "amanda-a");
 143     def_trig.npars = 1;
 144     strcpy(def_trig.pars[0], "type=external");
 145     def_trig.nwords = 0;
 146     rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
 147
 148     rdmc_init_array_hdef(&def_trig);
 149     def_trig.id = 2;
 150     strcpy(def_trig.tag, "amanda-b");
 151     def_trig.npars = 1;
 152     strcpy(def_trig.pars[0], "type=majority");
 153     def_trig.nwords = 0;
 154     rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
 155
 156     rdmc_init_array_hdef(&def_trig);
 157     def_trig.id  = 3;
 158     strcpy(def_trig.tag, "spase-1");
 159     def_trig.npars = 1;
 160     strcpy(def_trig.pars[0], "type=external");
 161     def_trig.nwords = 0;
 162     rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
 163
 164     rdmc_init_array_hdef(&def_trig);
 165     def_trig.id  = 4;
 166     strcpy(def_trig.tag, "spase-2");
 167     def_trig.npars = 1;
 168     strcpy(def_trig.pars[0], "type=external");
 169     def_trig.nwords = 0;
 170     rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
 171
 172     rdmc_init_array_hdef(&def_trig);
 173     def_trig.id = 5;
 174     strcpy(def_trig.tag, "gasp");
 175     def_trig.npars = 1;
 176     strcpy(def_trig.pars[0], "type=external");
 177     def_trig.nwords = 0;
 178     rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
 179   }
 180
 181   /* now set the geo cal flag */
 182   ar->is_calib.geo=1;
 183
 184   /* now add one FRESULT line */
 185   {
 186     array_hdef_t fd;
 187     rdmc_init_array_hdef(&fd);
 188     rdmc_add_fit_def(ar,&fd,0);
 189     rdmc_jk_to_fitdef(ar->def_fit, 1);
 190   }
 191   return 0;
 192
 193 } /* rarr_uwi() */
 194
 195 /****************************************************************************/
 196 /* revt_uwi() reads the next UWI format event                               */
 197 /****************************************************************************/
 198
 199 int rdmc_revt_uwi(mcfile *fp, mevt *ev, const array *ar)
 200 {
 201
 202   char s[RDMC_MAXLINE];                                            /* input line */
 203   int c;                                                      /* first char */
 204   int header_read = 0;               /* detect if the event header was read */
 205
 206   if (fp->fmajor < 2)  return ERRLINE; /* only newer formats are supported */
 207
 208   rdmc_clear_mevt(ev);                           /* clear old event info */
 209
 210   if (feof(fp->fp))                          /* if file end already reached */
 211     return EOF;
 212
 213   ev->nrun = ar->nrun;                       /* copy the run number from fp */
 214
 215   for (*s = c = getc(fp->fp); c != EOF; *s=c=getc(fp->fp)) {/*for all lines */
 216     fp->fpos += 1;
 217     if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL) {        /* try to read a line */
 218       return EOF;
 219     }
 220
 221     if (strncmp(s,"DH",2) == 0) {                      /* data event header */
 222       if (header_read++) return ERRLINE;
 223       if (uwi_rd_DH(s+2, ev)) return ERRLINE;
 224     } else if (strncmp(s,"EH",2) == 0) {                /* MC event header */
 225       if (header_read++) return ERRLINE;
 226       if (uwi_rd_EH(s+2, ev)) return ERRLINE;
 227     } else if (strncmp(s,"MU",2) == 0) {          /* Muon track information */
 228       if (uwi_rd_MU(s+2, &(ev->gen), &(ev->ntrack))) return ERRLINE;
 229     } else if (strncmp(s,"SH",2) == 0) {             /* shower information */
 230       if (uwi_rd_SH(s+2, &(ev->gen), &(ev->ntrack))) return ERRLINE;
 231     } else if (strncmp(s,"FT",2) == 0) {                    /* Fit results */
 232       if (uwi_rd_FT(s+2, &(ev->rec), &(ev->fresult), &(ev->nfit)))
 233         return ERRLINE;
 234     } else if (strncmp(s,"SP",2) == 0) {               /* Spase Fit results */
 235       if (uwi_rd_SP(s+2, &(ev->rec), &(ev->fresult), &(ev->nfit)))
 236         return ERRLINE;
 237     } else if (strncmp(s,"HT",2) == 0) {                /* Hit information */
 238       if (uwi_rd_HT(s+2, &(ev->h), &(ev->nhits))) return ERRLINE;
 239     } else if (strncmp(s,"EN",2) == 0) {                      /* Event end */
 240
 241    /* now fill necessary data structures */
 242       rdmc_repair_gen_id(ev);
 243       ev->nch = rdmc_count_nch(ev);     /* calc the number of hit channels */
 244       rdmc_fill_mhit_str(ev,ar);
 245       ev->nstr = rdmc_count_nstr(ev);
 246
 247
 248
 249       return 0;
 250     } else                                                   /* unknown id */
 251       continue;  /* do nothing */
 252
 253   } /* for all lines */
 254
 255   if (c == EOF)
 256     return EOF;
 257
 258   ev->nch=rdmc_count_nch(ev);                         /* calc the number of hit channels */
 259
 260   return 0;
 261
 262 } /* revt_uwi() */
 263
 264
 265 /****************************************************************************/
 266 /* skipevt_uwi() skips the next event from a UWI-like file                  */
 267 /****************************************************************************/
 268
 269 int rdmc_skipevt_uwi(mcfile *fp)
 270 {
 271   char s[RDMC_MAXLINE];                                            /* input line */
 272   int c;                                                      /* first char */
 273
 274   fp->fpos += 1;
 275   if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL)        /* try to read first line */
 276     return EOF;
 277
 278   for (*s = c = getc(fp->fp); c != EOF; *s=c=getc(fp->fp)) {/*for all lines */
 279     fp->fpos += 1;
 280     if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL) {        /* try to read a line */
 281       return EOF;
 282     }
 283
 284     if (strncmp(s,"EN",2) == 0)                                /* Event end */
 285       return 0;
 286
 287   } /* for all lines */
 288
 289
 290   return 0;
 291
 292 } /* function skipevt_ascii() */
 293
 294
 295 /****************************************************************************/
 296 /* warr_uwi() may write the geometry to a UWI like geometry file            */
 297 /*  This function is somehow funny: It first looks if there is already a    */
 298 /*  geometry file. If yes, it compares them, and returns a zero if the      */
 299 /*  geometries are the same, and and error if not.                          */
 300 /*  only if there is no geometry file, we will write it.                    */
 301 /****************************************************************************/
 302
 303 int rdmc_warr_uwi(char *geofile,const array *geo)
 304 {
 305   FILE *fp;
 306   array stored_array;
 307   int r;
 308   int istr, ich, nchstr;
 309   float xstr, ystr;
 310
 311   if (geofile == NULL)      /* use the default file (environement variable) */
 312     geofile = getenv(DEFAULT_UWI_GEO_FILE_ENV);
 313   if (geofile == NULL)                              /* use the default file */
 314     geofile = DEFAULT_UWI_GEO_FILE;                            /* (builtin) */
 315
 316   fp = fopen(geofile,"r");                  /* try to open the file to read */
 317   if (fp != NULL) {                    /* there was already a geometry file */
 318     fclose(fp);                                      /* just close it again */
 319     if ((r = rdmc_rarr_uwi(geofile,&stored_array)))    /* read the geometry */
 320       return r;                                             /* if errornous */
 321     if (rdmc_comp_array(geo, &stored_array) == 1)        /* compare arrays */
 322       return ERRLINE;
 323     else
 324       return 0;
 325   } else {                                 /* if there was no geometry file */
 326     fp = fopen(geofile,"w");          /* open the geometry file for writing */
 327     if (fp == NULL) return ERRLINE;
 328     fprintf(fp, "%.0f\t\t! Depth of detector center\n", geo->depth);
 329     fprintf(fp, "%i\t\t! Number of strings\n", geo->nstr);
 330     for (istr = 1; istr <= geo->nstr; istr++) {          /* for all strings */
 331       nchstr = 0;
 332       xstr = 0.0; ystr = 0.0;
 333       for (ich = 0; ich < geo->nch; ich++)                   /* calc nchstr */
 334         if (geo->str[ich] == istr) {
 335           nchstr++;
 336           xstr += geo->x[ich];
 337           ystr += geo->y[ich];
 338         } /* if str[ich] == istr */
 339       if (nchstr > 0) {
 340         xstr /= nchstr;                             /* mean x, y coordinate */
 341         ystr /= nchstr;
 342         fprintf(fp, "\n%.3f %.3f\t! x, y of string %i\n",xstr, ystr, istr);
 343         fprintf(fp, "%i\t\t! Nr of OMs\n\n",nchstr);
 344         nchstr = 0;
 345         for (ich = 0; ich < geo->nch; ich++)
 346           if (geo->str[ich] == istr) {
 347             nchstr++;
 348             fprintf(fp, "%.3f\t\t! z of OM %i\n",geo->z[ich],nchstr);
 349             if (geo->costh[ich] < -0.5) fprintf(fp, "down.pmt\n");
 350             else if (geo->costh[ich] > 0.5) fprintf(fp, "up.pmt\n");
 351             else fprintf(fp, "unknown.pmt\n");
 352         } /* if str[ich] == istr */
 353       } /* if nchstr */
 354
 355     } /* for istr */
 356     fclose(fp);
 357     return 0;
 358   } /* if fp == NULL */
 359
 360 } /* warr_uwi() */
 361
 362 /****************************************************************************/
 363 /* whd_uwi() writes the header to UWI file                                  */
 364 /****************************************************************************/
 365
 366 int rdmc_whd_uwi(const mcfile *fp)
 367 {
 368 #if 0
 369   int day;
 370   day = gmtime(&(fp->time))->tm_yday;
 371 #endif
 372
 373   if (fp->info.uwi.mc_id == 'D')
 374     fprintf(fp->fp, "FH %i DATA unknown\n", UWI_ASCII_VERSION);
 375   else if (fp->info.uwi.mc_id == 'M')
 376     fprintf(fp->fp, "FH %i MC %s\n",UWI_ASCII_VERSION, fp->info.uwi.mc_vers);
 377   else
 378     fprintf(fp->fp, "FH %i MC unknown\n", UWI_ASCII_VERSION);
 379
 380   return 0;
 381
 382 } /* whd_uwi() */
 383
 384 /****************************************************************************/
 385 /* function wevt_uwi() writes an event to an UWI file                       */
 386 /****************************************************************************/
 387
 388 int rdmc_wevt_uwi(const mcfile *fp,const mevt *event, const array *ar)
 389 {
 390   int nmuon = 0;
 391   int nshower = 0;
 392   int imuon = -1;
 393   int shwr_type;
 394   int itrack;
 395   int ihit;
 396   int gpsyear, gpsday;
 397
 398   /* look how many muons we have */
 399   for (itrack = 0; itrack < event->ntrack; itrack++)
 400     if ((event->gen[itrack].id == MUON_PLUS)
 401         || (event->gen[itrack].id == MUON_MINUS)) {
 402       nmuon++;
 403       if (imuon < 0) imuon = itrack; /* pointer to first muon */
 404     } else
 405       nshower++;
 406
 407   /* write event header  */
 408   if (fp->info.uwi.mc_id == 'D') {
 409     rdmc_mjd_to_gps(event->mjd, &gpsyear, &gpsday);
 410     fprintf(fp->fp, "DH %i %i %i.%06li %i.%06li %i %li %i\n",
 411             gpsyear, gpsday, event->secs,
 412             rdmc_nint(event->nsecs*1e-3),
 413             event->secs + 86400 * gpsday,
 414             rdmc_nint(event->nsecs*1e-3),
 415             event->enr, event->trigger, event->nhits);
 416   } else {
 417     if (nmuon > 0)
 418       fprintf(fp->fp, "EH %i %i %i %.0f %.2f %.2f %.0f\n",
 419               event->enr, nmuon, event->nch,
 420               -1.0, /* primary energy */
 421               acos(event->gen[imuon].costh)/PID180,
 422               event->gen[imuon].phi/PID180,
 423               -1.0 /* AMU mass of primary */
 424               );
 425     else
 426       fprintf(fp->fp, "EH %i %i %i\n", event->enr, 0, event->nch);
 427   }
 428
 429   /* write (gen) muons */
 430   nmuon = 0;
 431   for (itrack = 0; itrack < event->ntrack; itrack++) {
 432     float len;
 433     shwr_type = -1;
 434     switch (event->gen[itrack].id) {
 435     case MUON_PLUS:
 436     case MUON_MINUS:
 437       nmuon++;
 438       for (nshower = 1; itrack+nshower < event->ntrack; nshower++)
 439         if ((event->gen[itrack+nshower].id == MUON_PLUS)
 440             || (event->gen[itrack+nshower].id == MUON_MINUS))
 441           break;
 442
 443       fprintf(fp->fp, "MU %i %.0f %i %.2f %.2f ",
 444               nmuon, event->gen[itrack].e * 1e-3,
 445               nshower,
 446               acos(event->gen[itrack].costh)/PID180,
 447               event->gen[itrack].phi/PID180);
 448       len = ( event->gen[itrack].length < 0. )
 449         ? RDMC_BIG : event->gen[itrack].length ;
 450       fprintf(fp->fp, "%.1f %.1f %.1f %.1f %.1f %.1f %.1f\n",
 451               event->gen[itrack].x * 1e2,
 452               event->gen[itrack].y * 1e2,
 453               event->gen[itrack].z * 1e2,
 454               (event->gen[itrack].x +
 455                len  * event->gen[itrack].px) * 1e2,
 456               (event->gen[itrack].y +
 457                len * event->gen[itrack].py) * 1e2,
 458               (event->gen[itrack].z +
 459                len * event->gen[itrack].pz) * 1e2,
 460               0.0 /* closest distance to origin */
 461               );
 462       break;
 463     } /* switch id */
 464   } /* for itrack */
 465
 466   /* write (gen) showers */
 467   nmuon = 0;
 468   nshower = 0;
 469   for (itrack = 0; itrack < event->ntrack; itrack++) {
 470     shwr_type = -1;
 471     switch (event->gen[itrack].id) {
 472     case MUON_PLUS:
 473     case MUON_MINUS:
 474       nmuon++;
 475       nshower = 0;
 476       break;
 477     case BREMS: if (shwr_type < 0) shwr_type = 1;
 478     case PAIRPROD: if (shwr_type < 0) shwr_type = 2;
 479     case NUCL_INT: if (shwr_type < 0) shwr_type = 3;
 480     case DELTAE: if (shwr_type < 0) shwr_type = 4;
 481     case MU_PAIR: if (shwr_type < 0) shwr_type = 5;
 482     default:                                 /* store unknown particles */
 483       nshower++;
 484       fprintf(fp->fp, "SH %i %i %i %.0f %.0f %.1f %.1f %.1f\n",
 485               nmuon, /* muon responsible for this shower */
 486               nshower+1, shwr_type,
 487               event->gen[itrack].e * 1e-3,
 488               0.0, /* light output from shower */
 489               event->gen[itrack].x * 1e2,
 490               event->gen[itrack].y * 1e2,
 491               event->gen[itrack].z * 1e2);
 492     } /* switch id */
 493   } /* for itrack */
 494
 495   /* write hits */
 496   for (ihit = 0; ihit < event->nhits; ihit++)
 497     fprintf(fp->fp, "HT %i %i %i %.1f %.1f %.1f\n",
 498             event->h[ihit].ch+1,
 499             0, /* nr of p.e. generated */
 500             1, /* nr. of output pulses generated */
 501             event->h[ihit].amp * UWI_CHANNELS_PER_PE,
 502             event->h[ihit].t,
 503             event->h[ihit].tot);
 504
 505   /* write fits */
 506   for (itrack = 0; itrack < event->nfit; itrack++){
 507     fprintf(fp->fp, "FT %.1f %.1f %.1f %.1f %.1f",
 508             acos(event->rec[itrack].costh)/PID180,
 509             event->rec[itrack].phi/PID180,
 510             event->rec[itrack].x * 1e2,
 511             event->rec[itrack].y * 1e2,
 512             event->rec[itrack].z * 1e2);
 513
 514     if ((event->fresult != NULL)              /* if there is an jk record */
 515         && (0 <= event->fresult[itrack].id )
 516         && (ar->n_fit > event->fresult[itrack].id )  /* there is a fit defined  */
 517         && (rdmc_is_this_jk(&(ar->def_fit[event->fresult[itrack].id])
 518                                ,&(event->fresult[itrack])) ) ) {
 519       fprintf(fp->fp, " %f", event->fresult[itrack].val[JK_CHI2]);
 520     }else{
 521       fprintf(fp->fp, " %f",(float) RDMC_NA);
 522     }
 523     fprintf(fp->fp, " %f %i\n" , event->rec[itrack].t,  0);
 524
 525   }
 526   /* write end of event */
 527   fprintf(fp->fp,"EN\n");
 528
 529   return 0;
 530 } /* wevt_uwi() */
 531
 532 /****************************************************************************/
 533 /* rhd_uwi() reads the format relevant informations for UWI like            */
 534 /*          formats                                                         */
 535 /****************************************************************************/
 536
 537 int rdmc_rhd_uwi(mcfile *fp)
 538 {
 539   return 0;
 540 }
 541
 542 int rdmc_uwi_FH(const char *s, mcfile *fp)
 543 {
 544   int form;
 545   int format, mcversion;
 546   char filename[RDMC_MAXLINE];
 547   char histline[RDMC_MAXLINE];
 548
 549   form = sscanf(s,"FH %i MC %i", &format, &mcversion);
 550   if (form >= 2) {                                         /* MC file */
 551     fp->fmajor = format;
 552     fp->info.uwi.mc_id = 'M';
 553     sprintf(fp->info.uwi.mc_vers,"%i",mcversion);
 554     return 0;
 555   }
 556   form = sscanf(s, "FH %i DATA %s", &format, filename);
 557   if (form >= 1) {                                        /* data file */
 558     fp->fmajor = format;
 559     fp->fminor = 0;
 560     fp->info.uwi.mc_id = 'D';
 561     strcpy(fp->info.uwi.mc_vers,"0.0");
 562
 563     sprintf(histline,"uwicalib (%i) %s\n",format, (form > 1)?filename:"");
 564     if (fp->creator == NULL) {
 565       fp->creator = malloc(strlen(histline)+2);
 566       strcpy(fp->creator, "!");
 567       strcat(fp->creator, histline);
 568     } else {
 569       fp->creator = realloc(fp->creator,
 570                             strlen(fp->creator) + strlen(histline)+1);
 571       strcat(fp->creator, "!");
 572       strcat(fp->creator, histline);
 573     }
 574     return 0;
 575   }
 576
 577   return ERRLINE;
 578
 579 } /* rhd_uwi() */
 580
 581 /****************************************************************************/
 582 /* Functions for reading UWI format files                                   */
 583 /****************************************************************************/
 584
 585 /* Read the "EH " line (just ignore the primary) */
 586
 587 int uwi_rd_EH(char *s, mevt *ev)
 588 {
 589   int form;
 590   float prim_energy;
 591   float prim_theta, prim_phi;
 592   float prim_mass;
 593   int nparts;
 594   int nhits;
 595
 596   form = sscanf(s,"%i %i %i %f %f %f %f",
 597                 &(ev->enr), &(nparts), &nhits,
 598                 &prim_energy, &prim_theta, &prim_phi, &prim_mass);
 599
 600
 601   if (form == 0)
 602     ev->enr = 0;
 603
 604   return 0;
 605
 606 } /* uwi_rd_EH() */
 607
 608 /* read a Data event header */
 609
 610 int uwi_rd_DH(char *s, mevt *ev)
 611 {
 612   int form;
 613   int reg;                  /* Data register */
 614   int nhits;
 615   int year, day, sec, nsec; /* time */
 616   float gmt, fnsecs;
 617   int enr;
 618   char nsecstr[RDMC_MAXLINE];
 619   char nsecfloat[RDMC_MAXLINE];
 620
 621   form = sscanf(s, "%i %i %i.%s %f %i %i %i",
 622                 &year, &day, &sec, nsecstr, &gmt, &enr, &reg, &nhits);
 623
 624   switch(form) {
 625   case 0: year = 70;
 626   case 1: day = 0;
 627   case 2: sec = 0;
 628   case 3: strcpy(nsecstr, "0");
 629   case 4: gmt = 0.0;
 630   case 5: enr = 0;
 631   case 6: reg = 0;
 632   case 7: nhits = 0;
 633   }
 634
 635   sprintf(nsecfloat, "0.%s",nsecstr);
 636   sscanf(nsecfloat, "%f", &fnsecs);
 637   nsec = 1e6 * fnsecs;
 638   nsec *= 1000;
 639
 640   ev->mjd = rdmc_gps_to_mjd(year, day);  /* convert gps date into mjd */
 641   ev->secs = sec;
 642   ev->nsecs = nsec;
 643   ev->trigger = reg;   /* I hope! that reg contains the trigger bitfield */
 644                        /* (but where can I get the trigger info?)        */
 645   ev->enr = enr;
 646
 647   return 0;
 648
 649 } /* uwi_rd_DH() */
 650
 651 /* read a muon track and reserve the memory for all shower tracks */
 652
 653 int uwi_rd_MU(char *s, mtrack **tr, int *ntrack)
 654 {
 655   int form;
 656   float l;
 657   int muon_number, num_showers;
 658   float muon_energy;
 659   float theta_mu, phi_mu;
 660   float xstart, ystart, zstart;
 661   float xend, yend, zend;
 662   mtrack muon;
 663   int i;
 664
 665   form = sscanf(s, "%i %f %i %f %f %f %f %f %f %f %f",
 666                 &muon_number, &muon_energy, &num_showers,
 667                 &theta_mu, &phi_mu,
 668                 &xstart, &ystart, &zstart,
 669                 &xend, &yend, &zend);
 670
 671   switch (form) {
 672   case 8:
 673     xend = xstart;
 674     yend = ystart;
 675     zend = zstart;
 676   case 11:
 677   case 12: break;
 678   default: return ERRLINE;
 679   }
 680
 681   if (num_showers < 1) return ERRLINE;
 682
 683   *tr = (mtrack *)realloc(*tr, sizeof(mtrack) * (*ntrack+num_showers));
 684
 685   muon.id = MUON_PLUS;
 686   muon.e = muon_energy * 1e3; /* file: GeV, rdmc: MeV */
 687   muon.x = xstart*1e-2;
 688   muon.y = ystart*1e-2;
 689   muon.z = zstart*1e-2;
 690   muon.phi = phi_mu * PID180;
 691   muon.costh = - cos(theta_mu*PID180);
 692   l = (xstart - xend) * (xstart - xend)
 693     + (ystart - yend) * (ystart - yend)
 694     + (zstart - zend) * (zstart - zend);
 695   if (l < 0.)
 696     l = RDMC_NA;
 697   else
 698     l = sqrt(l);
 699   muon.length = l * 1e-2; /* file: cm, rdmc: m */
 700
 701   if (l < 1e-12) /* if the track is too short for a good calculation */
 702     rdmc_tau_tr
 703 (*tr);     /* calc direction cosinuus from phi and costh */
 704   else {
 705     muon.px = (xend - xstart) / l;/* else calc them from end-start */
 706     muon.py = (yend - ystart) / l;
 707     muon.pz = (zend - zstart) / l;
 708   }
 709
 710   muon.t = 0.;
 711   muon.nmuon = 1;
 712
 713   for (i = 0; i < num_showers; i++)            /* copy the muon to the event */
 714     memcpy((*tr)+*ntrack+i, &muon, sizeof(mtrack));   /* (dummy for showers) */
 715
 716   for (i = 1; i < num_showers; i++)
 717     (*tr)[*ntrack+i].id = -1;                         /* dummies for showers */
 718
 719   (*ntrack) += num_showers;
 720
 721   return 0;
 722
 723 } /* uwi_rd_MU() */
 724
 725 /* read a shower and fill it to the (in uwi_rd_MU allocated) track */
 726
 727 int uwi_rd_SH(char *s, mtrack **tr, int *ntrack)
 728 {
 729   int form;
 730   int which_muon, which_shower, showr_type;
 731   float shwr_energy, shwr_photons;
 732   float x, y, z;
 733   int i;
 734   mtrack shower;
 735
 736   form = sscanf(s, "%i %i %i %f %f %f %f %f",
 737                 &which_muon, &which_shower,
 738                 &showr_type,
 739                 &shwr_energy, &shwr_photons,
 740                 &x, &y, &z);
 741
 742   if (form < 8) return ERRLINE;
 743
 744   switch (showr_type) {
 745   case 1: shower.id = BREMS; break;
 746   case 2: shower.id = PAIRPROD; break;
 747   case 3: shower.id = NUCL_INT; break;
 748   case 4: shower.id = DELTAE; break;
 749   case 5: shower.id = MU_PAIR; break;
 750   default: shower.id = -1;
 751   }
 752
 753   shower.e = shwr_energy * 1e3;
 754   shower.x = x * 1e-2;
 755   shower.y = y * 1e-2;
 756   shower.z = z * 1e-2;
 757
 758   /* lets look for the muon number to fill into the track direction */
 759
 760   for (i = 0; (which_muon > 0) && (i < *ntrack); i++)
 761     if (((*tr)[i].id == MUON_PLUS) || ((*tr)[i].id == MUON_MINUS))
 762       which_muon--;
 763
 764   if (i >= *ntrack) return ERRLINE;
 765
 766   shower.px = (*tr)[i].px;
 767   shower.py = (*tr)[i].py;
 768   shower.pz = (*tr)[i].pz;
 769   shower.phi = (*tr)[i].phi;
 770   shower.costh = (*tr)[i].costh;
 771
 772   shower.length = 0; /* a shower has zero length (lets assume)  */
 773   shower.nmuon = 0;
 774
 775   shower.t = 0; /* normally, we should calculate the time which
 776                          * corresponds to the  track I do not do this.
 777                          */
 778
 779   /* insert the shower at the right place (after the corresponding muon) */
 780
 781   memcpy((*tr) + i + which_shower - 2, &shower, sizeof(mtrack));
 782
 783   return 0;
 784
 785 } /* uwi_rd_SH() */
 786
 787 int uwi_rd_HT(char *s, mhit **h, int *nhits)
 788 {
 789   int form;
 790   int om_num, n_pulses;
 791   float adc, n_pe;
 792   float tdc0, tdc1, tdc2, tdc3, tdc4;
 793   float tot0, tot1, tot2, tot3, tot4;
 794
 795   form = sscanf(s, "%i %f %i %f %f %f %f %f %f %f %f %f %f %f",
 796                 &om_num, &n_pe, &n_pulses,
 797                 &adc,
 798                 &tdc0, &tot0,
 799                 &tdc1, &tot1,
 800                 &tdc2, &tot2,
 801                 &tdc3, &tot3,
 802                 &tdc4, &tot4);
 803
 804   /*
 805    * We will skip all hits without hits!!! (That means, without TDC)
 806    */
 807
 808   if (form < 4)
 809     return ERRLINE;
 810   if (n_pulses < 0) return ERRLINE;
 811
 812   if ((n_pulses == 0) && (form == 6))
 813     n_pulses = 1;
 814
 815   if (n_pulses > 5) n_pulses = 5;
 816   if (form < (4 + 2*n_pulses) )
 817     n_pulses = (form - 4 )/2;
 818
 819   *h = (mhit *)realloc(*h, sizeof(mhit) * (*nhits + n_pulses));
 820
 821   switch(n_pulses) {
 822   case 5:
 823     (*h)[*nhits+4].ch = om_num-1;
 824     (*h)[*nhits+4].str = 0;  /* default: no string (will be set at the end) */
 825     (*h)[*nhits+4].mt = RDMC_NA;
 826     (*h)[*nhits+4].ma = RDMC_NA;
 827     (*h)[*nhits+4].t = tdc4; /* units are nanonseconds */
 828     (*h)[*nhits+4].tot = tot4;
 829     (*h)[*nhits+4].amp = 0.;
 830   case 4:
 831     (*h)[(*nhits)+3].ch = om_num-1;
 832     (*h)[(*nhits)+3].str = 0;
 833     (*h)[(*nhits)+3].mt = RDMC_NA;
 834     (*h)[(*nhits)+3].ma = RDMC_NA;
 835     (*h)[(*nhits)+3].t = tdc3;
 836     (*h)[(*nhits)+3].tot = tot3;
 837     (*h)[(*nhits)+3].amp = 0.;
 838   case 3:
 839     (*h)[(*nhits)+2].ch = om_num-1;
 840     (*h)[(*nhits)+2].str = 0;
 841     (*h)[(*nhits)+2].mt = RDMC_NA;
 842     (*h)[(*nhits)+2].ma = RDMC_NA;
 843     (*h)[(*nhits)+2].t = tdc2;
 844     (*h)[(*nhits)+2].tot = tot2;
 845     (*h)[(*nhits)+2].amp = 0.;
 846   case 2:
 847     (*h)[(*nhits)+1].ch = om_num-1;
 848     (*h)[(*nhits)+1].str = 0;
 849     (*h)[(*nhits)+1].mt = RDMC_NA;
 850     (*h)[(*nhits)+1].ma = RDMC_NA;
 851     (*h)[(*nhits)+1].t = tdc1;
 852     (*h)[(*nhits)+1].tot = tot1;
 853     (*h)[(*nhits)+1].amp = 0.;
 854   case 1:
 855     (*h)[*nhits].ch = om_num-1;
 856     (*h)[*nhits].str = 0;
 857     (*h)[*nhits].mt = RDMC_NA;
 858     (*h)[*nhits].ma = RDMC_NA;
 859     (*h)[*nhits].t = tdc0;
 860     (*h)[*nhits].tot = tot0;
 861     (*h)[*nhits].amp = adc / UWI_CHANNELS_PER_PE;
 862     break;
 863   case 0:
 864     break;
 865   default: return ERRLINE;
 866   }
 867
 868   (*nhits) += n_pulses;
 869
 870   return 0;
 871 } /* uwi_rd_HT() */
 872
 873 int uwi_rd_SP(char *s, mtrack **fit, mevt_special_t **fitres, int *nfit)
 874 { /* Tmiller: SP gps_time theta phi x y Nparticles dummy(currently 0) */
 875   /* serap:   SP gmt, theta, phi, x, y, s30, xxx */
 876
 877   int form;
 878   float theta_fit, phi_fit;
 879   float x_fit, y_fit;
 880   float s30,xxx;
 881   double gmt;
 882
 883   form = sscanf(s, "%lf %f %f %f %f %f %f",
 884                 &gmt,&theta_fit, &phi_fit,
 885                 &x_fit, &y_fit, &s30,
 886                 &xxx);
 887
 888   switch (form) {
 889   case 5: s30 = 0;
 890   case 6: xxx = 0.;
 891   case 7: break;
 892   default: return ERRLINE;
 893   }
 894
 895   *fit = (mtrack *)realloc(*fit, sizeof(mtrack) * (*nfit+1));
 896   *fitres = (mevt_special_t *)
 897     realloc( *fitres , sizeof(mevt_special_t)  *  (*nfit+1));
 898   rdmc_init_mtrack(&((*fit)[*nfit]));
 899   rdmc_init_fit_jk(&((*fitres)[*nfit]),RDMC_NA);
 900
 901
 902   (*fit)[*nfit].id = MUON_PLUS;
 903   (*fit)[*nfit].e = 0.0;
 904   (*fit)[*nfit].x = x_fit;
 905   (*fit)[*nfit].y = y_fit;
 906   (*fit)[*nfit].z = 0.0;
 907   (*fit)[*nfit].phi = phi_fit * PID180;
 908   (*fit)[*nfit].costh = cos(theta_fit*PID180);
 909   (*fit)[*nfit].length = RDMC_BIG;
 910   (*fit)[*nfit].tag = *nfit;
 911
 912   rdmc_tau_tr(*fit + *nfit); /* calc direction cosinus from phi and costh */
 913
 914   (*fit)[*nfit].nmuon = 1;  /* we dont know this, lets assume one muon */
 915
 916   (*fitres)[*nfit].id = 0;
 917   (*fitres)[*nfit].val[JK_FITID] = 9000;
 918   (*fitres)[*nfit].val[JK_CHI2] = 0.;
 919   (*fitres)[*nfit].val[JK_RCHI2] = 0.;
 920   (*fitres)[*nfit].val[JK_PROB] = 0.;
 921   (*fitres)[*nfit].val[JK_SIGTH] = 0.;
 922   (*fitres)[*nfit].val[JK_COVMIN] = s30;
 923   (*fitres)[*nfit].val[JK_COVMAX] = xxx;
 924   (*fitres)[*nfit].val[JK_CUTFLAG] = -1;
 925
 926   (*nfit)++;
 927
 928   return 0;
 929
 930 } /* uwi_rd_FT() */
 931
 932 int uwi_rd_FT(char *s, mtrack **fit, mevt_special_t  **fitres, int *nfit)
 933 {
 934   int form;
 935   float theta_fit, phi_fit;
 936   float x_fit, y_fit, z_fit;
 937   float chi2;
 938   int n_in_time;
 939   float t_fit;
 940
 941   form = sscanf(s, "%f %f %f %f %f %f %f %i",
 942                 &theta_fit, &phi_fit,
 943                 &x_fit, &y_fit, &z_fit,
 944                 &chi2, &t_fit, &n_in_time);
 945
 946   switch (form) {
 947   case 5: chi2 = 0.0;
 948   case 6: t_fit = 0.0;
 949   case 7: n_in_time = 0;
 950   case 8: break;
 951   default: return ERRLINE;
 952   }
 953
 954   *fit = (mtrack *)realloc(*fit, sizeof(mtrack) * (*nfit+1));
 955   *fitres = (mevt_special_t *)
 956     realloc( *fitres , sizeof(mevt_special_t)  *  (*nfit+1));
 957   rdmc_init_mtrack(&((*fit)[*nfit]));
 958   rdmc_init_fit_jk(&((*fitres)[*nfit]),RDMC_NA);
 959
 960
 961
 962   (*fit)[*nfit].id = MUON_PLUS;
 963   (*fit)[*nfit].e = 0.0;
 964   (*fit)[*nfit].x = x_fit*1e-2;
 965   (*fit)[*nfit].y = y_fit*1e-2;
 966   (*fit)[*nfit].z = z_fit*1e-2;
 967   (*fit)[*nfit].t = t_fit;
 968   (*fit)[*nfit].phi = phi_fit * PID180;
 969   (*fit)[*nfit].costh = cos(theta_fit*PID180);
 970   (*fit)[*nfit].length = RDMC_BIG;
 971   (*fit)[*nfit].tag = *nfit;
 972
 973   rdmc_tau_tr(*fit + *nfit); /* calc direction cosinuus from phi and costh */
 974
 975   (*fit)[*nfit].nmuon = 1;     /* we dont know this, lets assume one muon */
 976
 977   (*fitres)[*nfit].id = 0;
 978   (*fitres)[*nfit].val[JK_FITID] = 0;
 979   (*fitres)[*nfit].val[JK_CHI2] = chi2;
 980   (*fitres)[*nfit].val[JK_RCHI2] = -1.0;
 981   (*fitres)[*nfit].val[JK_PROB] = 0.0;
 982   (*fitres)[*nfit].val[JK_SIGTH] = 0.0;
 983   (*fitres)[*nfit].val[JK_COVMIN] = 0.0;
 984   (*fitres)[*nfit].val[JK_COVMAX] = 0.0;
 985   (*fitres)[*nfit].val[JK_CUTFLAG] = -1;
 986
 987   (*nfit)++;
 988
 989   return 0;
 990
 991 } /* uwi_rd_FT() */
 992
 993 #endif /* UWI_ASCII_F */
 994
 995 /****************************************************************************/
 996 /********************************** E O F ***********************************/
 997 /****************************************************************************/
 998 /*
 999    This is just for EMACS:
1000    Local Variables:
1001    compile-command: "cd .. ; make -k rdmc"
1002    End:
1003 */