2 * Read/write files in the format used at UWI Madison (see J.Jacobsen's PhD
7 "$Header: /net/local/cvsroot/siegmund/rdmc/uwi.c,v 1.33 1999/12/23 12:23:15 wiebusch Exp $";
10 #define _USE_MATH_DEFINES
16 ///////////////#include <unistd.h>
22 #if defined(OSF1) || defined(AIX) || defined (SunOS)
28 #include "rdmc_local.h"
33 #define DEBUG 0 /* DEBUG 1 writes some info to stderr */
36 #define ERRLINE (-__LINE__)
38 #define PID180 (M_PI/180.)
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);
52 /****************************************************************************/
53 /* The function rarr_uwi() reads the geometry of a UWI like file */
54 /****************************************************************************/
55 int rdmc_rarr_uwi(char *geofile, array *ar)
59 char omfile[RDMC_MAXLINE];
65 int nchstr; /* temp. x, y, z coor, number of channels per string */
67 rdmc_init_array(ar); /* reset the array */
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) */
74 fp = fopen(geofile,"r");
76 if (fp == NULL) return ERRLINE; /* no geometry file found */
79 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
81 form = sscanf(s,"%f",&(ar->depth)); /* 1. line: detector depth */
84 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
86 form = sscanf(s,"%i",&(ar->nstr)); /* 2. line: Number of strings */
88 for (istr = 0; istr < ar->nstr; istr++) { /* read all strings */
90 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
92 form = sscanf(s,"%f %f",&x,&y); /* 3. line: x,y of 1st string */
95 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
97 form = sscanf(s,"%i",&nchstr); /* 4. line: nch of 1st string */
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 */
102 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
104 form = sscanf(s,"%f",&z); /* 5. line: z of 1st channel */
107 if (fgets(s,RDMC_MAXLINE,fp) == NULL)
109 form = sscanf(s,"%s",omfile); /* 6. line: om file of 1st ch */
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;
127 ar->costh[ich] = 0.0;
128 ar->sensit[ich] = 0.0;
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 */
138 array_hdef_t def_trig;
140 rdmc_init_array_hdef(&def_trig);
142 strcpy(def_trig.tag, "amanda-a");
144 strcpy(def_trig.pars[0], "type=external");
146 rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
148 rdmc_init_array_hdef(&def_trig);
150 strcpy(def_trig.tag, "amanda-b");
152 strcpy(def_trig.pars[0], "type=majority");
154 rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
156 rdmc_init_array_hdef(&def_trig);
158 strcpy(def_trig.tag, "spase-1");
160 strcpy(def_trig.pars[0], "type=external");
162 rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
164 rdmc_init_array_hdef(&def_trig);
166 strcpy(def_trig.tag, "spase-2");
168 strcpy(def_trig.pars[0], "type=external");
170 rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
172 rdmc_init_array_hdef(&def_trig);
174 strcpy(def_trig.tag, "gasp");
176 strcpy(def_trig.pars[0], "type=external");
178 rdmc_add_trigger_def(ar,&def_trig,ar->n_trigger);
181 /* now set the geo cal flag */
184 /* now add one FRESULT line */
187 rdmc_init_array_hdef(&fd);
188 rdmc_add_fit_def(ar,&fd,0);
189 rdmc_jk_to_fitdef(ar->def_fit, 1);
195 /****************************************************************************/
196 /* revt_uwi() reads the next UWI format event */
197 /****************************************************************************/
199 int rdmc_revt_uwi(mcfile *fp, mevt *ev, const array *ar)
202 char s[RDMC_MAXLINE]; /* input line */
203 int c; /* first char */
204 int header_read = 0; /* detect if the event header was read */
206 if (fp->fmajor < 2) return ERRLINE; /* only newer formats are supported */
208 rdmc_clear_mevt(ev); /* clear old event info */
210 if (feof(fp->fp)) /* if file end already reached */
213 ev->nrun = ar->nrun; /* copy the run number from fp */
215 for (*s = c = getc(fp->fp); c != EOF; *s=c=getc(fp->fp)) {/*for all lines */
217 if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL) { /* try to read a line */
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)))
234 } else if (strncmp(s,"SP",2) == 0) { /* Spase Fit results */
235 if (uwi_rd_SP(s+2, &(ev->rec), &(ev->fresult), &(ev->nfit)))
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 */
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);
250 } else /* unknown id */
251 continue; /* do nothing */
253 } /* for all lines */
258 ev->nch=rdmc_count_nch(ev); /* calc the number of hit channels */
265 /****************************************************************************/
266 /* skipevt_uwi() skips the next event from a UWI-like file */
267 /****************************************************************************/
269 int rdmc_skipevt_uwi(mcfile *fp)
271 char s[RDMC_MAXLINE]; /* input line */
272 int c; /* first char */
275 if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL) /* try to read first line */
278 for (*s = c = getc(fp->fp); c != EOF; *s=c=getc(fp->fp)) {/*for all lines */
280 if (fgets(s+1,RDMC_MAXLINE-1,fp->fp) == NULL) { /* try to read a line */
284 if (strncmp(s,"EN",2) == 0) /* Event end */
287 } /* for all lines */
292 } /* function skipevt_ascii() */
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 /****************************************************************************/
303 int rdmc_warr_uwi(char *geofile,const array *geo)
308 int istr, ich, nchstr;
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) */
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 */
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 */
332 xstr = 0.0; ystr = 0.0;
333 for (ich = 0; ich < geo->nch; ich++) /* calc nchstr */
334 if (geo->str[ich] == istr) {
338 } /* if str[ich] == istr */
340 xstr /= nchstr; /* mean x, y coordinate */
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);
345 for (ich = 0; ich < geo->nch; ich++)
346 if (geo->str[ich] == istr) {
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 */
358 } /* if fp == NULL */
362 /****************************************************************************/
363 /* whd_uwi() writes the header to UWI file */
364 /****************************************************************************/
366 int rdmc_whd_uwi(const mcfile *fp)
370 day = gmtime(&(fp->time))->tm_yday;
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);
378 fprintf(fp->fp, "FH %i MC unknown\n", UWI_ASCII_VERSION);
384 /****************************************************************************/
385 /* function wevt_uwi() writes an event to an UWI file */
386 /****************************************************************************/
388 int rdmc_wevt_uwi(const mcfile *fp,const mevt *event, const array *ar)
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)) {
403 if (imuon < 0) imuon = itrack; /* pointer to first muon */
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);
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 */
426 fprintf(fp->fp, "EH %i %i %i\n", event->enr, 0, event->nch);
429 /* write (gen) muons */
431 for (itrack = 0; itrack < event->ntrack; itrack++) {
434 switch (event->gen[itrack].id) {
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))
443 fprintf(fp->fp, "MU %i %.0f %i %.2f %.2f ",
444 nmuon, event->gen[itrack].e * 1e-3,
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 */
466 /* write (gen) showers */
469 for (itrack = 0; itrack < event->ntrack; itrack++) {
471 switch (event->gen[itrack].id) {
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 */
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);
496 for (ihit = 0; ihit < event->nhits; ihit++)
497 fprintf(fp->fp, "HT %i %i %i %.1f %.1f %.1f\n",
499 0, /* nr of p.e. generated */
500 1, /* nr. of output pulses generated */
501 event->h[ihit].amp * UWI_CHANNELS_PER_PE,
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);
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]);
521 fprintf(fp->fp, " %f",(float) RDMC_NA);
523 fprintf(fp->fp, " %f %i\n" , event->rec[itrack].t, 0);
526 /* write end of event */
527 fprintf(fp->fp,"EN\n");
532 /****************************************************************************/
533 /* rhd_uwi() reads the format relevant informations for UWI like */
535 /****************************************************************************/
537 int rdmc_rhd_uwi(mcfile *fp)
542 int rdmc_uwi_FH(const char *s, mcfile *fp)
545 int format, mcversion;
546 char filename[RDMC_MAXLINE];
547 char histline[RDMC_MAXLINE];
549 form = sscanf(s,"FH %i MC %i", &format, &mcversion);
550 if (form >= 2) { /* MC file */
552 fp->info.uwi.mc_id = 'M';
553 sprintf(fp->info.uwi.mc_vers,"%i",mcversion);
556 form = sscanf(s, "FH %i DATA %s", &format, filename);
557 if (form >= 1) { /* data file */
560 fp->info.uwi.mc_id = 'D';
561 strcpy(fp->info.uwi.mc_vers,"0.0");
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);
569 fp->creator = realloc(fp->creator,
570 strlen(fp->creator) + strlen(histline)+1);
571 strcat(fp->creator, "!");
572 strcat(fp->creator, histline);
581 /****************************************************************************/
582 /* Functions for reading UWI format files */
583 /****************************************************************************/
585 /* Read the "EH " line (just ignore the primary) */
587 int uwi_rd_EH(char *s, mevt *ev)
591 float prim_theta, prim_phi;
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);
608 /* read a Data event header */
610 int uwi_rd_DH(char *s, mevt *ev)
613 int reg; /* Data register */
615 int year, day, sec, nsec; /* time */
618 char nsecstr[RDMC_MAXLINE];
619 char nsecfloat[RDMC_MAXLINE];
621 form = sscanf(s, "%i %i %i.%s %f %i %i %i",
622 &year, &day, &sec, nsecstr, &gmt, &enr, ®, &nhits);
628 case 3: strcpy(nsecstr, "0");
635 sprintf(nsecfloat, "0.%s",nsecstr);
636 sscanf(nsecfloat, "%f", &fnsecs);
640 ev->mjd = rdmc_gps_to_mjd(year, day); /* convert gps date into mjd */
643 ev->trigger = reg; /* I hope! that reg contains the trigger bitfield */
644 /* (but where can I get the trigger info?) */
651 /* read a muon track and reserve the memory for all shower tracks */
653 int uwi_rd_MU(char *s, mtrack **tr, int *ntrack)
657 int muon_number, num_showers;
659 float theta_mu, phi_mu;
660 float xstart, ystart, zstart;
661 float xend, yend, zend;
665 form = sscanf(s, "%i %f %i %f %f %f %f %f %f %f %f",
666 &muon_number, &muon_energy, &num_showers,
668 &xstart, &ystart, &zstart,
669 &xend, ¥d, &zend);
678 default: return ERRLINE;
681 if (num_showers < 1) return ERRLINE;
683 *tr = (mtrack *)realloc(*tr, sizeof(mtrack) * (*ntrack+num_showers));
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);
699 muon.length = l * 1e-2; /* file: cm, rdmc: m */
701 if (l < 1e-12) /* if the track is too short for a good calculation */
703 (*tr); /* calc direction cosinuus from phi and costh */
705 muon.px = (xend - xstart) / l;/* else calc them from end-start */
706 muon.py = (yend - ystart) / l;
707 muon.pz = (zend - zstart) / l;
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) */
716 for (i = 1; i < num_showers; i++)
717 (*tr)[*ntrack+i].id = -1; /* dummies for showers */
719 (*ntrack) += num_showers;
725 /* read a shower and fill it to the (in uwi_rd_MU allocated) track */
727 int uwi_rd_SH(char *s, mtrack **tr, int *ntrack)
730 int which_muon, which_shower, showr_type;
731 float shwr_energy, shwr_photons;
736 form = sscanf(s, "%i %i %i %f %f %f %f %f",
737 &which_muon, &which_shower,
739 &shwr_energy, &shwr_photons,
742 if (form < 8) return ERRLINE;
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;
753 shower.e = shwr_energy * 1e3;
758 /* lets look for the muon number to fill into the track direction */
760 for (i = 0; (which_muon > 0) && (i < *ntrack); i++)
761 if (((*tr)[i].id == MUON_PLUS) || ((*tr)[i].id == MUON_MINUS))
764 if (i >= *ntrack) return ERRLINE;
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;
772 shower.length = 0; /* a shower has zero length (lets assume) */
775 shower.t = 0; /* normally, we should calculate the time which
776 * corresponds to the track I do not do this.
779 /* insert the shower at the right place (after the corresponding muon) */
781 memcpy((*tr) + i + which_shower - 2, &shower, sizeof(mtrack));
787 int uwi_rd_HT(char *s, mhit **h, int *nhits)
790 int om_num, n_pulses;
792 float tdc0, tdc1, tdc2, tdc3, tdc4;
793 float tot0, tot1, tot2, tot3, tot4;
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,
805 * We will skip all hits without hits!!! (That means, without TDC)
810 if (n_pulses < 0) return ERRLINE;
812 if ((n_pulses == 0) && (form == 6))
815 if (n_pulses > 5) n_pulses = 5;
816 if (form < (4 + 2*n_pulses) )
817 n_pulses = (form - 4 )/2;
819 *h = (mhit *)realloc(*h, sizeof(mhit) * (*nhits + n_pulses));
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.;
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.;
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.;
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.;
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;
865 default: return ERRLINE;
868 (*nhits) += n_pulses;
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 */
878 float theta_fit, phi_fit;
883 form = sscanf(s, "%lf %f %f %f %f %f %f",
884 &gmt,&theta_fit, &phi_fit,
885 &x_fit, &y_fit, &s30,
892 default: return ERRLINE;
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);
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;
912 rdmc_tau_tr(*fit + *nfit); /* calc direction cosinus from phi and costh */
914 (*fit)[*nfit].nmuon = 1; /* we dont know this, lets assume one muon */
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;
932 int uwi_rd_FT(char *s, mtrack **fit, mevt_special_t **fitres, int *nfit)
935 float theta_fit, phi_fit;
936 float x_fit, y_fit, z_fit;
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);
949 case 7: n_in_time = 0;
951 default: return ERRLINE;
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);
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;
973 rdmc_tau_tr(*fit + *nfit); /* calc direction cosinuus from phi and costh */
975 (*fit)[*nfit].nmuon = 1; /* we dont know this, lets assume one muon */
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;
993 #endif /* UWI_ASCII_F */
995 /****************************************************************************/
996 /********************************** E O F ***********************************/
997 /****************************************************************************/
999 This is just for EMACS:
1001 compile-command: "cd .. ; make -k rdmc"