]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/icepack/iceconvert/uwi.c
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
12-oct-2005 NvE CleanTasks() invoked before execution of sub-tasks in IceConvert...
[u/mrichter/AliRoot.git] / RALICE / icepack / iceconvert / uwi.c
CommitLineData
f67e2651 1/*
2 * Read/write files in the format used at UWI Madison (see J.Jacobsen's PhD
3 * thesis)
4 */
5
6char *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
42int uwi_rd_EH(char *s, mevt *ev); /* read the several lines for UWI files */
43int uwi_rd_DH(char *s, mevt *ev);
44int uwi_rd_MU(char *s, mtrack **tr, int *ntrack);
45int uwi_rd_SH(char *s, mtrack **tr, int *ntrack);
46int uwi_rd_HT(char *s, mhit **h, int *nhits);
47int uwi_rd_FT(char *s, mtrack **fitr, mevt_special_t **fitres, int *nfit);
48int 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/****************************************************************************/
55int 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
199int 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
269int 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
303int 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
366int 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
388int 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
537int rdmc_rhd_uwi(mcfile *fp)
538{
539 return 0;
540}
541
542int 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
587int 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
610int 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
653int 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
727int 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
787int 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
873int 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
932int 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*/