]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/icepack/iceconvert/baikal.c
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
08-mar-2006 NvE Time offset correction in IceF2k extended to allow also a user define...
[u/mrichter/AliRoot.git] / RALICE / icepack / iceconvert / baikal.c
CommitLineData
f67e2651 1/*
2 * Read/write files in the BAIKAL formats
3 * Zeuthen: Krabi/Wischnewski
4 * Moscow: an older Nora/Belolaptikov format (only reading)
5 */
6
7char *baikal_rdmc_cvsid =
8"$Header: /net/local/cvsroot/siegmund/rdmc/baikal.c,v 1.20 2000/02/07 13:02:17 ole Exp $";
9
10#include <stdio.h>
11#define _USE_MATH_DEFINES
12#include <math.h>
13#include <string.h>
14#include <stdlib.h>
15#include <ctype.h>
16
17////////////#include <unistd.h>
18
19
20#include <fcntl.h>
21#include <sys/stat.h>
22#include <errno.h>
23
24#include "rdmc.h"
25
26#include "rdmc_local.h"
27#include "baikal.h" /* function prototyping */
28
29
30#ifndef DEBUG
31#define DEBUG 0 /* DEBUG 1 writes some info to stderr */
32#endif
33
34#define ERRLINE (-__LINE__)
35
36#define PID180 (M_PI/180.)
37
38#ifdef BAIKAL_BIN_F
39
40#define BAIKAL_PMTPERCHANNEL 2
41
42
43
44
45static void rdmc_add_fitdef_baikal(array *ar);
46
47
48/****************************************************************************/
49/* The function rarr_mc() reads the header of a baikal like file */
50/****************************************************************************/
51
52int rdmc_rarr_baikal_mc(mcfile *fp, array *ar)
53{
54
55 int nw; /* number of words to read */
56 int ptr; /* pointer to the current word in the header */
57 int i,j; /* loop variables */
58 int *iclust; /* pointer to current channel number in string */
59
60 rdmc_init_array(ar); /* reset the array */
61
62/**************************** the array info *******************************/
63 ar->nrun=fp->info.bai.nrun;
64
65 ptr = 0;
66 if (fp->info.bai.fortran_recs)
67 rdmc_mcgeti(fp); /* skip fortran record header */
68 nw = rdmc_mcgeti(fp); /* 1: number of words */
69
70 if (nw <= 2) { /* too less words: nonreconstructed data file */
71 ar->id = NT_36; /* default: baikal array */
72 ar->nch = 0; /* no telescope info aviable */
73 ar->nstr = 0;
74 rdmc_mcseek(fp, nw - ptr); /* read the unused fields */
75 if (fp->info.bai.fortran_recs)
76 rdmc_mcgeti(fp); /* unused dummy */
77 return 0; /* and return */
78 } /* if nw <= 2 */
79
80 rdmc_mcseek(fp,2); /* here unused values */
81 ar->nch = rdmc_mcgetf(fp); /* 3: number of channels */
82 ar->nstr = rdmc_mcgetf(fp); /* 4: number of strings */
83 if (ar->nch <= 18)
84 ar->id = NT_36;
85 else if (ar->nch <= 36)
86 ar->id = NT_72;
87 else if (ar->nch <= 48)
88 ar->id = NT_96;
89 else if (ar->nch <= 72)
90 ar->id = NT_144;
91 else if (ar->nch <= 96)
92 ar->id = NT_192;
93 else
94 ar->id = NT_200;
95 ptr += 4; /* 3 words read */
96
97 if (nw < ptr + ar->nch * 5) return ERRLINE; /* invalid format */
98 if (ar->nch > RDMC_MAXCHANNELS) return ERRLINE;
99 if (ar->nch <= 0) return ERRLINE;
100 if (ar->nstr <= 0) return ERRLINE;
101
102 iclust = (int *)calloc(ar->nstr,sizeof(int));
103 /* alloc mem for channel nr array */
104
105/*************************** Channel position and parameters ***************/
106
107 for (i = 0; i < ar->nch; i++) {
108 ar->x[i] = rdmc_mcgetf(fp)/100.0; /* coordinates */
109 ar->y[i] = rdmc_mcgetf(fp)/100.0; /* in meters of the */
110 ar->z[i] = rdmc_mcgetf(fp)/100.0; /* i-th channel */
111 ar->costh[i] = rdmc_mcgetf(fp); /* cos theta */
112 ar->str[i] = rdmc_mcgetf(fp); /* string number */
113 ar->type[i] = QUASAR; /* default pmt type */
114 ar->clust[i] = iclust[ar->str[i] - 1]++; /* channel number */
115 } /* for i */
116 ptr += 5 * ar->nch;
117 ar->is_calib.geo=1;
118
119 free (iclust); /* free the mem for the ch nr array */
120
121/****************************** Pmt parameters *****************************/
122
123 if (nw < ptr + ar->nch * 2 * BAIKAL_PMTPERCHANNEL) return ERRLINE;
124 for (i = 0; i < ar->nch; i++) {
125 ar->thresh[i] = rdmc_mcgetf(fp); /* threshold for each pmt */
126 ar->sensit[i] = rdmc_mcgetf(fp); /* sensitivity of each pmt */
127 for (j = 1; j < BAIKAL_PMTPERCHANNEL; j++) {
128 rdmc_mcgetf(fp); /* threshold for each pmt */
129 rdmc_mcgetf(fp); /* sensitivity of each pmt */
130 } /* for j */
131 ar->serial[i] = 0; /* no serial number */
132 } /* for i */
133 ptr += ar->nch * 2 * BAIKAL_PMTPERCHANNEL;
134
135 rdmc_mcseek(fp, nw - ptr); /* read the last unused fields */
136
137 if (fp->info.bai.fortran_recs)
138 rdmc_mcgeti(fp); /* unused dummy */
139
140 rdmc_add_fitdef_baikal(ar);
141
142 return 0;
143
144} /* function rarr_mc() */
145
146/****************************************************************************/
147/* revt_baikal() tries to read the different BAIKAL dialects */
148/****************************************************************************/
149int rdmc_revt_baikal(mcfile *fp, mevt *ev, const array *ar)
150{
151 long filepos; /* shows the current file position */
152 int r;
153
154 filepos = ftell(fp->fp); /* store the file position */
155
156 if (fp->info.bai.mc) { /* if it is (probably) a mc file */
157 r = rdmc_revt_baikal_mc(fp,ev,ar); /* try to read as an mc event */
158 if (r < EOF) {
159 fseek(fp->fp,filepos,SEEK_SET); /* set to old file position */
160 r = rdmc_revt_baikal_data(fp,ev,ar); /* try to read as data file */
161 if (r == 0) /* if it was successfull */
162 fp->info.bai.mc = 0; /* set file type to data */
163 else /* if not */
164 fseek(fp->fp,filepos,SEEK_SET); /* restore file position */
165 } /* if r == ILF */
166 } /* if mc file */
167
168 else { /* if it is (probably) a data file */
169 r = rdmc_revt_baikal_data(fp,ev,ar); /* try to read as data event */
170 if (r < EOF) { /* if it fails */
171 fseek(fp->fp,filepos,SEEK_SET); /* set to old file position */
172 r = rdmc_revt_baikal_mc(fp,ev,ar); /* try to read as an mc event */
173 if (r == 0) /* if it was successfull */
174 fp->info.bai.mc = 1; /* set file type to mc */
175 else /* if not */
176 fseek(fp->fp,filepos,SEEK_SET); /* restore file position */
177 } /* if r == ILF */
178 } /* if data file */
179
180 fp->fpos = ftell(fp->fp); /* store file position for debugging */
181
182 /* now update stuff which was not included in the data */
183 /* now uopdate the channel and srtring reference of the event */
184 ev->nch = rdmc_count_nch(ev); /* calc the number of hit channels */
185 rdmc_fill_mhit_str(ev,ar); /* fill the mhit->str fields */
186 ev->nstr = rdmc_count_nstr(ev); /* calc the number of hit channels */
187
188
189
190 return r;
191
192} /* revt_baikal() */
193
194/****************************************************************************/
195/* revt_mc() reads the next *mc*-event */
196/* it is just a switch to the correct routines (moscow/zeuthen format) */
197/****************************************************************************/
198
199int rdmc_revt_baikal_mc(mcfile *fp, mevt *ev, const array *ar)
200{
201 switch(fp->fmajor) {
202 case BAIKAL_OLD:
203 case BAIKAL_NEW:
204 return rdmc_revt_baikal_zeu(fp,ev,ar);
205 case BAIKAL_MOS:
206 return rdmc_revt_baikal_mos(fp,ev,ar);
207 default:
208 return ERRLINE;
209 }
210
211} /* function revt_mc() */
212
213/****************************************************************************/
214/* revt_zeu() reads the zeuthen formats (BAIKAL_OLD,BAIKAL_NEW) */
215/****************************************************************************/
216
217int rdmc_revt_baikal_zeu(mcfile *fp, mevt *ev, const array *ar)
218{
219
220 int nw; /* number of words to read */
221 int ptr = 0; /* pointer to the current word in the header */
222 int i; /* loop variable */
223
224 rdmc_clear_mevt(ev); /* clear old event info */
225
226 if (fp->info.bai.fortran_recs)
227 rdmc_mcseek(fp, 1); /* skip fortran record header */
228
229 nw = rdmc_mcgeti(fp) - 1; /* number of words in the record */
230 if (feof(fp->fp)) return EOF; /* file end reached? */
231
232 if (nw < 8) return ERRLINE;
233
234/***************** read the first (common) informations ********************/
235
236 ev->t_offset = 0.0; /* time offset in baikal is zero */
237
238 if (fp->fmajor == BAIKAL_NEW) { /* for the new format */
239 i = rdmc_mcgeti(fp);
240 ev->nhits = i % 1000; /* number of hits (higher part is nrun)*/
241 ev->nrun = i/1000;
242 ev->enr = rdmc_mcgeti(fp); /* event number */
243 fp->info.bai.enr = ev->enr;
244 ptr += 2;
245 } /* if new format */
246 else { /* for the old format */
247 ev->nrun = ar->nrun;
248 ev->enr = ++(fp->info.bai.enr); /* set dummies */
249 } /* if old format */
250
251/************ allocate memory for generated track and fill it **************/
252
253 {
254 mtrack gen;
255 rdmc_init_mtrack(&gen);
256 gen.costh = rdmc_mcgetf(fp); /* cos theta + 10 * nhits */
257
258 if (fp->fmajor == BAIKAL_OLD) { /* for the old format */
259 ev->nhits = (fabs(gen.costh)+2.0)/10; /* fill the nr of hits */
260 if (ev->nhits == 0) return ERRLINE; /* there must be >= 1 hit! */
261 } /* if old format */
262 else /* for the new format */
263 if ((int)((fabs(gen.costh)+2.0)/10.0) != ev->nhits) return ERRLINE;
264 /* test for nhits */
265 gen.costh = fmod(gen.costh,10.);
266 /* this is now really cos theta */
267 gen.phi = rdmc_mcgetf(fp)*PID180; /* phi in rad */
268 gen.x = rdmc_mcgetf(fp)/100.0; /* read the */
269 gen.y = rdmc_mcgetf(fp)/100.0; /* coordinates */
270 gen.z = rdmc_mcgetf(fp)/100.0; /* of the track */
271 gen.e = rdmc_mcgetf(fp)*1000.0; /* Energy in MeV */
272 gen.id = MUON_PLUS; /* we generated muons */
273 gen.tag = 1;
274 gen.length = RDMC_BIG; /* infinity track length */
275 rdmc_tau_tr(&gen); /* calc the direction cosinuus */
276 rdmc_add_gen(ev,&gen,0);
277 }
278 ptr += 5; /* increase record pointer */
279
280 if (nw < ptr + ev->nhits * 3 + 2) /* if there are not enough words */
281 return ERRLINE;
282
283/************* allocate memory for the hits and fill them ******************/
284
285 if (ev->nhits > 0)
286 ev->h = malloc(ev->nhits * sizeof(mhit));
287 /* allocate mem for the hits */
288 for (i = 0; i < ev->nhits; i++)
289 rdmc_init_mhit( &(ev->h[i]));
290 for (i = 0; i < ev->nhits; i++) {
291 ev->h[i].ch = rdmc_mcgetf(fp); /* contains ch and m */
292 ev->h[i].ma = ev->h[i].ch/10000; /* the first 2 digits */
293 ev->h[i].mt = fmod(ev->h[i].ch/100,100.); /* the digits 3 and 4 */
294 ev->h[i].ch = fmod(ev->h[i].ch,100.); /* the last 2 digits */
295 ev->h[i].ch--; /* for "C" style (index 0,...) */
296 ev->h[i].amp = rdmc_mcgetf(fp); /* amplitude */
297 ev->h[i].t = rdmc_mcgetf(fp); /* time */
298 ev->h[i].tot = -1.0; /* no TOT information */
299 } /* for i */
300
301
302 rdmc_mcseek(fp, 1); /* one dummy word */
303 ev->gen->t = rdmc_mcgetf(fp); /* time corresponding gen->xyz */
304 ev->gen->nmuon = ev->ntrack; /* number of muons */
305
306 if (ev->gen->costh < -1.5) { /* if costh = -2 (no gen. track) */
307 free(ev->gen); /* free the gen. track */
308 ev->gen = NULL; /* reset the pointer */
309 ev->ntrack = 0; /* set the number of tracks to zero */
310 } /* if costh < -1.5 */
311
312 ptr += ev->nhits * 3 + 2;
313
314 if (nw <= ptr + 2) { /* if no reconstruction */
315 ev->nfit = 0; /* number of reconstr. tracks is 0 */
316 rdmc_mcseek(fp, nw-ptr); /* read the last (unused) fields */
317 if (fp->info.bai.fortran_recs)
318 rdmc_mcseek(fp, 1); /* read the field used for fortran format */
319 return 0;
320 } /* if no reco */
321
322 if (nw < ptr + 10) return ERRLINE;
323
324/*********** allocate memory for reconstructed track and fill it ***********/
325
326 {
327 mtrack rec;
328 mevt_special_t fresult;
329 rdmc_init_mtrack(&rec);
330 rdmc_init_fit_jk(&fresult,1);
331
332 rec.costh = rdmc_mcgetf(fp); /* cos theta, or -99 if no reco */
333 rec.phi = rdmc_mcgetf(fp)*PID180; /* reconstructed phi */
334 rec.x = rdmc_mcgetf(fp)/100.0; /* reconstructed vertex */
335 rec.y = rdmc_mcgetf(fp)/100.0;
336 rec.z = rdmc_mcgetf(fp)/100.0;
337 rec.t = rdmc_mcgetf(fp); /* corresponding time */
338 rec.id = MUON_PLUS; /* we reconstructed muons */
339 rec.e = 0.0; /* no energy reconstructed */
340 rec.length = RDMC_BIG; /* infinity track length */
341 rec.nmuon = 1; /* default: one muon */
342 rec.tag = 1;
343 rdmc_tau_tr(&rec); /* calc the direction cosinuus */
344
345 /************** read jaanus' reconstruction results ************************/
346
347 fresult.val[JK_SIGTH] = rdmc_mcgetf(fp);
348 fresult.val[JK_COVMAX] = rdmc_mcgetf(fp);
349 fresult.val[JK_COVMIN] = fmod(fresult.val[JK_COVMAX],1000.)/1000.;
350 fresult.val[JK_COVMAX] = (fresult.val[JK_COVMAX] /1000.
351 - fresult.val[JK_COVMIN])/1000.;
352 fresult.val[JK_CUTFLAG] = rdmc_mcgetf(fp);
353 fresult.val[JK_FITID] =
354 fp->info.bai.rec_vers + 10000 * fp->info.bai.min_vers;
355 ptr += 9;
356
357 rdmc_mcseek(fp, nw-ptr); /* read the last (unused) fields */
358
359 fresult.val[JK_RCHI2] = -1.0; /* in the baikal format is no rchi2 */
360 fresult.val[JK_CHI2] = -1.0; /* in the baikal format is no chi2 */
361 fresult.val[JK_PROB] = -1.0; /* in the baikal format is no prob */
362
363 if (rec.costh >= -2.0) { /* if cos theta is -99 */
364 rdmc_add_fit(ev,&rec,&fresult,0);
365 }
366 } /* read reco */
367
368 if (fp->info.bai.fortran_recs)
369 rdmc_mcseek(fp, 1); /* read the field used for fortran format */
370
371 return 0;
372
373} /* function revt_zeu() */
374
375/****************************************************************************/
376/* revt_mos() reads the next moscow format event (BAIKAL_MOS) */
377/****************************************************************************/
378
379int rdmc_revt_baikal_mos(mcfile *fp, mevt *ev, const array *ar)
380{
381
382 int nw; /* number of words to read */
383 int ptr = 0; /* pointer to the current word in the header */
384 int i,j; /* loop variable */
385 int nint; /* number of interactions of a muon */
386 mtrack gentrack; /* struct to hold the first track */
387
388 rdmc_clear_mevt(ev); /* clear old event info */
389
390 if (fp->info.bai.fortran_recs)
391 rdmc_mcseek(fp, 1); /* skip fortran record header */
392
393 nw = rdmc_mcgeti(fp) - 1; /* number of words in the record */
394 if (feof(fp->fp)) return EOF; /* file end reached? */
395
396 if (nw < 9) return ERRLINE;
397
398/***************** read the first (common) informations ********************/
399
400 ev->t_offset = 0.0; /* time offset in baikal is zero */
401
402
403 ev->nhits = rdmc_mcgetf(fp); /* number of hits */
404 ev->nrun = 0; /* no run number is stored */
405 ev->enr = ++(fp->info.bai.enr); /* event number from file pointer */
406 ptr += 1;
407
408/************ read the first generated track (central muon) ****************/
409
410 rdmc_init_mtrack(&gentrack);
411 gentrack.costh = -cos(rdmc_mcgetf(fp)*PID180); /* cos theta */
412 gentrack.phi = (180.0+rdmc_mcgetf(fp))*PID180; /* phi in rad */
413 if (gentrack.phi > 360.0) gentrack.phi -= 360.0; /* phi in (0,360) */
414 gentrack.x = rdmc_mcgetf(fp); /* read the */
415 gentrack.y = rdmc_mcgetf(fp); /* coordinates */
416 gentrack.z = rdmc_mcgetf(fp); /* of the track */
417 gentrack.t = rdmc_mcgetf(fp); /* time of the vertex */
418 gentrack.id = MUON_PLUS; /* we generated muons */
419 gentrack.length = RDMC_BIG; /* infinity track length */
420 gentrack.nmuon = 1; /* one muon produces this track */
421 gentrack.tag=1;
422 rdmc_tau_tr(&gentrack); /* calc the direction cosinuus */
423 gentrack.e = rdmc_mcgetf(fp) * 1e3; /* Energy in MeV */
424 ev->ntrack = (int)rdmc_mcgetf(fp) % 100; /* number of tracks */
425 /* higher part is nmuonall */
426 rdmc_mcgetf(fp); /* one dummy ("EMMAX + 1000* ESUM") */
427 ptr += 8; /* increase record pointer */
428
429/******* allocate memory for generated track(s) and fill them **************/
430
431#if 0 /* i dont know if I really need the first track (axis) */
432 ev->ntrack++; /* one additional track for the group axis */
433#endif
434
435 if (ev->ntrack > 0)
436 ev->gen = malloc(ev->ntrack*sizeof(mtrack));
437 /* allocate mem for generated track */
438
439#if 0 /* i dont know if I really need the first track (axis) */
440 memcpy(ev->gen,&gentrack,sizeof(mtrack)); /* fill the first track (axis) */
441
442 for (i = 1; i < ev->ntrack; i++){ /* for all tracks */
443#else
444 for (i = 0; i < ev->ntrack; i++){ /* for all tracks */
445#endif
446 rdmc_init_mtrack(ev->gen+i);
447 nint = rdmc_mcgetf(fp); /* read the number of interactions */
448 ev->gen[i].costh = gentrack.costh; /* theta */
449 ev->gen[i].phi = gentrack.phi; /* and phi are the same */
450 ev->gen[i].t = gentrack.t; /* and the time, too */
451 ev->gen[i].x = rdmc_mcgetf(fp); /* koordinates */
452 ev->gen[i].y = rdmc_mcgetf(fp); /* of the */
453 ev->gen[i].z = rdmc_mcgetf(fp); /* vertex */
454 rdmc_tau_tr(ev->gen+i); /* calc the direction cosinuus */
455 ev->gen[i].nmuon = 1; /* one muon produces this track */
456 ev->gen[i].length = RDMC_BIG; /*infinity length */
457 ev->gen[i].id = MUON_PLUS;
458 ev->gen[i].tag = i+1;
459 ev->gen[i].e = 0.0; /* reset the muon energy */
460 ptr += 4; /* 4 words read */
461 for (j = 0; j < nint; j++) { /* for all interactions */
462 ev->gen[i].e += rdmc_mcgetf(fp)*1e3; /* add the interaction energy */
463 rdmc_mcseek(fp,3); /* skip the interaction point */
464 ptr += 4; /* 4 words read */
465 } /* for j */
466
467 } /* for i */
468
469 ev->gen[0].e = gentrack.e; /* the first track gets all muon energy */
470
471 if (nw < ptr + ev->nhits * 5) /* if there are not enough words */
472 return ERRLINE;
473
474/************* allocate memory for the hits and fill them ******************/
475
476 if (ev->nhits > 0)
477 ev->h = malloc(ev->nhits * sizeof(mhit));
478 for (i = 0; i < ev->nhits; i++) {
479 rdmc_init_mhit(&(ev->h[i]));
480 }
481 /* allocate mem for the hits */
482 for (i = 0; i < ev->nhits; i++) {
483 ev->h[i].ch = rdmc_mcgetf(fp); /* channel number */
484 ev->h[i].ma = ((int)rdmc_mcgetf(fp)+1000) % 1000; /* amp-defining muon */
485 ev->h[i].mt = ((int)rdmc_mcgetf(fp)+1000) % 1000; /* time-defining muon */
486 ev->h[i].ch--; /* for "C" style (index 0,...) */
487 ev->h[i].amp = rdmc_mcgetf(fp); /* amplitude */
488 ev->h[i].t = rdmc_mcgetf(fp); /* time */
489 ev->h[i].tot = -1.0; /* no TOT information */
490 } /* for i */
491
492 ptr += ev->nhits * 5;
493
494 rdmc_mcseek(fp, nw-ptr); /* read the last (unused) fields */
495
496 if (fp->info.bai.fortran_recs)
497 rdmc_mcseek(fp, 1); /* read the field used for fortran format */
498
499 return 0;
500
501} /* function revt_mos() */
502
503/****************************************************************************/
504/* revt_data() reads the next *data* event */
505/****************************************************************************/
506
507int rdmc_revt_baikal_data(mcfile *fp, mevt *ev, const array *ar)
508{
509
510 int nw; /* number of words to read */
511 int ptr = 0; /* pointer to the current word in the header */
512 int i; /* loop variable */
513
514 rdmc_clear_mevt(ev); /* clear old event info */
515
516 if (fp->info.bai.fortran_recs)
517 rdmc_mcseek(fp, 1); /* skip fortran record header */
518
519 nw = rdmc_mcgeti(fp); /* number of words in the record */
520 if (feof(fp->fp)) return EOF; /* file end reached? */
521
522 if (nw < 2) return ERRLINE;
523
524/***************** read the first (common) informations ********************/
525
526 i = rdmc_mcgeti(fp); /* number of channels and run number */
527 ev->nhits = fmod(i,1000.); /* number of channels */
528 ev->nrun = i/1000; /* run number */
529
530 ev->enr = rdmc_mcgeti(fp); /* event number */
531 fp->info.bai.enr = ev->enr;
532 ev->ntrack = 0; /* no generation */
533
534 ptr += 2;
535 if (ptr + 3*ev->nhits > nw)/* if there are not enough words for hit info */
536 return ERRLINE;
537
538/************* allocate memory for the hits and fill them ******************/
539
540 if (ev->nhits > 0)
541 ev->h = malloc(ev->nhits * sizeof(mhit));
542 for (i = 0; i < ev->nhits; i++) { /* init */
543 rdmc_init_mhit( &(ev->h[i]));
544 }
545
546 for (i = 0; i < ev->nhits; i++) { /* read the hit infos */
547 ev->h[i].ch = rdmc_mcgetf(fp); /* channel nr (1,...) */
548 ev->h[i].ch--; /* for "C" style (index 0,...) */
549 ev->h[i].amp = rdmc_mcgetf(fp); /* amplitude of this hit */
550 ev->h[i].t = rdmc_mcgetf(fp); /* time of this hit */
551 ev->h[i].mt = ev->h[i].ma = RDMC_NA; /* there is no particle info */
552 ev->h[i].tot = -1.0; /* no TOT information */
553 } /* for i */
554
555
556 ptr += 3*ev->nhits;
557 if (nw <= ptr + 2) { /* if no reconstruction */
558 ev->nfit = 0; /* reset fit flag */
559 rdmc_mcseek(fp, nw-ptr); /* read the last (unused) fields */
560 if (fp->info.bai.fortran_recs)
561 rdmc_mcseek(fp, 1); /* read the field used for fortran format */
562 return 0;
563 } /* if no reco */
564
565 if (nw < ptr + 9) return ERRLINE;
566
567/*********** allocate memory for reconstructed track and fill it ***********/
568 {
569 mtrack rec;
570 mevt_special_t fresult;
571
572 rdmc_init_mtrack(&rec);
573 rdmc_init_fit_jk(&fresult,1);
574
575 rec.costh = rdmc_mcgetf(fp); /* cos theta, or -99 if no reco */
576 rec.phi = rdmc_mcgetf(fp)*PID180; /* reconstructed phi */
577 rec.x = rdmc_mcgetf(fp)/100.0; /* reconstructed coordinates */
578 rec.y = rdmc_mcgetf(fp)/100.0; /* of track origin */
579 rec.z = rdmc_mcgetf(fp)/100.0;
580 rec.t = rdmc_mcgetf(fp); /* reconstructed time of track origin */
581 rec.id = MUON_PLUS; /* we reconstructed muons */
582 rec.e = 0.0; /* no energy reconstructed */
583 rec.length = RDMC_BIG; /* infinity track length */
584 rec.nmuon = 1; /* default: one muon */
585 rec.tag = 1;
586 rdmc_tau_tr(&rec); /* calc the direction cosinuus */
587
588
589/************** read jaanus' reconstruction results ************************/
590
591 fresult.val[JK_SIGTH] = rdmc_mcgetf(fp);
592 fresult.val[JK_COVMAX] = rdmc_mcgetf(fp);
593 fresult.val[JK_COVMIN] = fmod(fresult.val[JK_COVMAX],1000.)/1000.;
594 fresult.val[JK_COVMAX] = (fresult.val[JK_COVMAX] /1000.
595 - fresult.val[JK_COVMIN])/1000.;
596 fresult.val[JK_CUTFLAG] = rdmc_mcgetf(fp);
597 ptr += 9;
598
599 rdmc_mcseek(fp, nw-ptr); /* read the last (unused) fields */
600
601 fresult.val[JK_RCHI2] = -1.0; /* in the baikal format is no rchi2 */
602 fresult.val[JK_CHI2] = -1.0; /* in the baikal format is no chi2 */
603 fresult.val[JK_PROB] = -1.0; /* in the baikal format is no prob */
604
605 if (ev->rec->costh >= -2.0) { /* if cos theta is -99 */
606 rdmc_add_fit(ev,&rec,&fresult,0);
607 } /* if no reco */
608 }
609 if (fp->info.bai.fortran_recs)
610 rdmc_mcseek(fp, 1); /* read the field used for fortran format */
611 return 0;
612
613} /* function revt_data() */
614
615/****************************************************************************/
616/* warr_mc() writes the array info to the baikal-like file */
617/****************************************************************************/
618
619int rdmc_warr_baikal_mc(mcfile *fp, const array *ar)
620{
621 int nw;
622 int i,j;
623 double mc_ind_muon, mc_id_spec; /* Jaanus' mc flags */
624 double mc_vers ;
625
626 nw = 14 + ar->nch * (5 + 2 * BAIKAL_PMTPERCHANNEL); /* number of words */
627 if (fp->info.bai.fortran_recs) /* if it is a formatted fortran file */
628 rdmc_mcputi(4 * nw + 4, fp); /* write the record length */
629
630 rdmc_mcputi(nw, fp); /* write the number of words as first int */
631
632 rdmc_mcputf(0.0, fp); /* write something (?) */
633
634 mc_vers=atof(fp->info.bai.mc_vers);
635 mc_vers = ((int) mc_vers) % 10000 ;
636 rdmc_mcputf(mc_vers + 10000.0 * fp->fmajor, fp); /* format version */
637 rdmc_mcputf((double)ar->nch, fp);
638 rdmc_mcputf((double)ar->nstr, fp);
639
640 for (i = 0; i < ar->nch; i++) { /* for all channels */
641 rdmc_mcputf(100.0*ar->x[i], fp);
642 rdmc_mcputf(100.0*ar->y[i], fp);
643 rdmc_mcputf(100.0*ar->z[i], fp);
644 rdmc_mcputf(ar->costh[i], fp);
645 rdmc_mcputf((double)ar->str[i], fp);
646 } /* for i */
647
648 for (i = 0; i < ar->nch; i++) {
649 for (j = 0; j < BAIKAL_PMTPERCHANNEL; j++) {
650 rdmc_mcputf(ar->thresh[i], fp); /* threshold for each pmt */
651 rdmc_mcputf(ar->sensit[i], fp); /* sensitivity of each pmt */
652 } /* for j */
653 } /* for i */
654
655 rdmc_mcputf((double)rdmc_o_rdateconv(fp->info.bai.time), fp); /* write time in YYMMDD format */
656
657#if 0
658 /* igen does contain the random seed since ASCII version -6 */
659 mc_ind_muon = fp->info.bai.igen /10000;
660 mc_id_spec = fp->info.bai.igen % 10000;
661#else
662 mc_ind_muon=mc_id_spec=0.;
663#endif
664
665 rdmc_mcputf(mc_ind_muon, fp);
666 rdmc_mcputf(mc_id_spec, fp);
667 rdmc_mcputf((double)fp->info.bai.igtrack, fp);
668 rdmc_mcputf(0.0, fp);
669 rdmc_mcputf(0.0, fp);
670
671 rdmc_mcputf((double)fp->info.bai.nw_rec_arr, fp);
672 rdmc_mcputf((double)fp->info.bai.nw_rec_evt, fp);
673 rdmc_mcputf((double)fp->info.bai.rec_vers, fp);
674 rdmc_mcputf((double)fp->info.bai.min_vers, fp);
675
676 if (fp->info.bai.fortran_recs) /* if it is a formatted fortran file */
677 rdmc_mcputi(4 * nw + 4, fp); /* write the record length */
678
679#if (DEBUG == 1)
680 fprintf(stderr, "<Write><Header: nch=%i nstr=%i>\n",ar->nch,ar->nstr);
681#endif
682
683 return errno;
684
685} /* function warr_mc() */
686
687/****************************************************************************/
688/* function wevt_mc() writes an event to file */
689/****************************************************************************/
690
691int rdmc_wevt_baikal_mc(mcfile *fp, const mevt *ev, const array *ar)
692{
693 int nw; /* number of words */
694 int i; /* loop variable */
695
696 nw = 20 + 3*ev->nhits;
697
698 if (fp->info.bai.fortran_recs) /* record length for formatted files */
699 rdmc_mcputi(4*nw+4,fp);
700
701/***************** write the first (common) informations ********************/
702
703 rdmc_mcputi(nw,fp); /* first word is number of words */
704 if (fp->fmajor == BAIKAL_NEW) { /* for the new format */
705 rdmc_mcputi(ev->nhits+1000*ev->nrun,fp);
706 rdmc_mcputi(ev->enr,fp);
707 } /* if new format */
708
709/******************* write the generating track *****************************/
710
711 if (ev->ntrack > 0) { /* if there is a generating track */
712 if (ev->gen[0].costh > 0)
713 rdmc_mcputf(ev->nhits*10.0+ev->gen[0].costh, fp);
714 else
715 rdmc_mcputf(-ev->nhits*10.0+ev->gen[0].costh, fp);
716 rdmc_mcputf(ev->gen[0].phi/PID180, fp);
717 rdmc_mcputf(ev->gen[0].x*100.0, fp);
718 rdmc_mcputf(ev->gen[0].y*100.0, fp);
719 rdmc_mcputf(ev->gen[0].z*100.0, fp);
720 rdmc_mcputf(ev->gen[0].e/1000.0, fp); /* baikal: energy in GeV */
721 } /* if gen */
722 else { /* if no generating track */
723 rdmc_mcputf(-ev->nhits*10.0-2.0, fp); /* seems to be not good ole */
724/* mcputf(ev->nhits*10.0, fp); */ /* this could be an alternative */
725 rdmc_mcputf(0.0, fp); /* set all to zero */
726 rdmc_mcputf(0.0, fp); /* set all to zero */
727 rdmc_mcputf(0.0, fp); /* set all to zero */
728 rdmc_mcputf(0.0, fp); /* set all to zero */
729 rdmc_mcputf(0.0, fp); /* set all to zero */
730 } /* if not gen */
731
732/***************************** write the hits *******************************/
733
734 for (i = 0; i < ev->nhits; i++) { /* store hits */
735 rdmc_mcputf(ev->h[i].ch+1.0+ev->h[i].mt*100.0+ev->h[i].ma*10000.0, fp);
736 rdmc_mcputf((double)ev->h[i].amp, fp);
737 rdmc_mcputf((double)ev->h[i].t, fp);
738 }
739
740 rdmc_mcputf(1.0, fp); /* I dont know what it is */
741
742 if (ev->ntrack > 0) /* if there is a generating track */
743 rdmc_mcputf(ev->gen[0].t, fp); /* time for origin of generated track */
744 else
745 rdmc_mcputf(0.0, fp);
746
747/******************* write the reconstructed track **************************/
748
749 if (ev->nfit > 0) { /* if there is a reconstrunction */
750 rdmc_mcputf(ev->rec[0].costh, fp);
751 rdmc_mcputf(ev->rec[0].phi/PID180, fp);
752 rdmc_mcputf(ev->rec[0].x*100.0, fp);
753 rdmc_mcputf(ev->rec[0].y*100.0, fp);
754 rdmc_mcputf(ev->rec[0].z*100.0, fp);
755 rdmc_mcputf(ev->rec[0].t, fp);
756 if ( (ev->fresult != NULL) /* if there is an jk record */
757 && (0 <= ev->fresult->id )
758 && (ar->n_fit > ev->fresult->id ) /* there is a fit defined */
759 && (rdmc_is_this_jk(&(ar->def_fit[ev->fresult->id])
760 ,ev->fresult) ) ) {
761 rdmc_mcputf(ev->fresult->val[JK_SIGTH], fp);
762 rdmc_mcputf(ev->fresult->val[JK_COVMIN]*1e3
763 + 1e6*ev->fresult->val[JK_COVMAX], fp);
764 rdmc_mcputf(ev->fresult->val[JK_CUTFLAG], fp);
765 } /* if res */
766 else { /* if there is no jk record */
767 rdmc_mcputf(0.0, fp); /* set all to zero */
768 rdmc_mcputf(0.0, fp); /* set all to zero */
769 rdmc_mcputf(0.0, fp); /* set all to zero */
770 } /* if not res */
771 } /* if reconstruction */
772 else { /* if there is no reconstruction */
773 rdmc_mcputf(-99., fp);
774 rdmc_mcputf(0.0, fp); /* set all to zero */
775 rdmc_mcputf(0.0, fp); /* set all to zero */
776 rdmc_mcputf(0.0, fp); /* set all to zero */
777 rdmc_mcputf(0.0, fp); /* set all to zero */
778 rdmc_mcputf(0.0, fp); /* set all to zero */
779 rdmc_mcputf(0.0, fp); /* set all to zero */
780 rdmc_mcputf(0.0, fp); /* set all to zero */
781 rdmc_mcputf(0.0, fp); /* set all to zero */
782 } /* if no reconstruction */
783
784 if (fp->fmajor == BAIKAL_NEW)
785 rdmc_mcputi(0, fp); /* flag for further subrecords */
786
787 if (fp->info.bai.fortran_recs) /* record length for formatted files */
788 rdmc_mcputi(4*nw+4,fp);
789
790 return 0;
791
792} /* function wevt() */
793
794
795/****************************************************************************/
796/* function wevt_data() writes an event to file */
797/****************************************************************************/
798
799int rdmc_wevt_baikal_data(mcfile *fp, const mevt *ev, const array *ar)
800{
801
802 return rdmc_wevt_baikal_mc(fp,ev,ar); /* we use mc format instead of data format */
803
804} /* function wevt_data() */
805
806/****************************************************************************/
807/* rhd_mc() reads the format relevant informations for baikal-like formats, */
808/* and then rewinds */
809/****************************************************************************/
810
811int rdmc_rhd_baikal_mc(mcfile *fp)
812{
813 long i1,i2;
814 float f1,f2,f3,f4;
815
816 i1 = rdmc_mcgeti(fp);
817 i2 = rdmc_mcgeti(fp);
818
819 if ((i1 == i2*4+4) || (rdmc_swap(i1) == rdmc_swap(i2) * 4 + 4)){/* if the first int */
820 fp->info.bai.fortran_recs = 1; /* contains the record length -> fortran file */
821 }
822 rewind(fp->fp); /* rewind file to the begin */
823 if (fp->info.bai.fortran_recs) rdmc_mcseek(fp,1); /* skip the fortran record header */
824 i1 = rdmc_mcgeti(fp); /* contains the number of words in header */
825 if (i1 == 2) { /* if there are only two words */
826 fp->info.bai.mc = 0; /* it is a non-reconstrcted data file */
827 fp->info.bai.nrun = rdmc_mcgetf(fp); /* first word is the run number */
828 fp->info.bai.time = 0; /* no time stored */
829 fp->info.bai.rec_vers = 0; /* no reconstruction */
830 fp->info.bai.min_vers = 0; /* no minimization */
831 rewind(fp->fp); /* rewind the file */
832 return 0; /* and return */
833 } /* if i1 = 2 */
834 f1 = rdmc_mcgetf(fp); /* dummy a(1)*/
835 f2 = rdmc_mcgetf(fp); /* contains mc and format version */
836 f3 = rdmc_mcgetf(fp); /* number of channels */
837 f4 = rdmc_mcgetf(fp); /* the number of strings */
838
839 if ((f4 != (int)f4) /* should be a integer */
840 ||(f4 <= 0.0) /* and greater than 0 */
841 ||(f4 > 20.0)) { /* and less than 20 (we havn't such big telescopes) */
842 fp->info.bai.swap = 1; /* perhaps the bytes are swapped */
843 } /* try it with swapped bytes! */
844 rewind(fp->fp); /* rewind file to the begin */
845 if (fp->info.bai.fortran_recs) rdmc_mcseek(fp,1); /* skip the fortran record header */
846 i1 = rdmc_mcgeti(fp); /* contains the number of words in header */
847 if (i1 <= 2) { /* if there are too less words */
848 fp->info.bai.mc = 0; /* it is a non-reconstrcted data file */
849 fp->info.bai.nrun = rdmc_mcgetf(fp); /* first word is the run number */
850 fp->info.bai.time = 0; /* no time stored */
851 fp->info.bai.rec_vers = 0; /* no reconstruction */
852 fp->info.bai.min_vers = 0; /* no minimization */
853 rewind(fp->fp); /* rewind the file */
854 return 0; /* and return */
855 } /* if i1 = 2 */
856 f1 = rdmc_mcgetf(fp); /* dummy a(1) */
857 f2 = rdmc_mcgetf(fp); /* contains mc and format version */
858 f3 = rdmc_mcgetf(fp); /* number of channels */
859 f4 = rdmc_mcgetf(fp); /* the number of strings */
860
861 fp->fmajor = (int) f2 / 10000; /* format version in second word */
862 fp->fminor = 0;
863 sprintf(fp->info.bai.mc_vers,"%i", ((int) f2) % 10000);
864
865 if ((f1 >= 0.9) && (f1 <= 10.0)&& (fp->fmajor == BAIKAL_OLD))
866 fp->fmajor = BAIKAL_MOS; /* 1.0 as 1st word: moscow format */
867
868 rdmc_mcseek(fp, (int)f3*9); /* skip the pmt data */
869 fp->info.bai.time = rdmc_o_dateconv(rdmc_mcgetf(fp)); /* time in sec since 1.1.70 */
870 fp->info.bai.igen = 10000 * rdmc_mcgetf(fp); /* k+2 is mc model */
871 fp->info.bai.igen += rdmc_mcgetf(fp); /* k+3 is mc spectrum */
872 if (fp->info.bai.igen >= 10000 )
873 fp->info.bai.igen = 103;
874 else
875 fp->info.bai.igen = 0;
876#if 0
877 fp->igtrack = rdmc_mcgetf(fp); /* k+4 is mc trigger */
878#else
879 fp->info.bai.nrun = rdmc_mcgetf(fp); /*ignore them */
880#endif
881 rdmc_mcseek(fp, 2); /* two dummies */
882#if 0
883 fp->info.bai.nrun = fp->igtrack; /* it seems to be so */
884#endif
885
886 i1 -= 10 + f3 * 9; /* decrease the header length s a pointer */
887
888 if (i1-- > 0) fp->info.bai.nw_rec_arr = rdmc_mcgetf(fp);
889 if (i1-- > 0) fp->info.bai.nw_rec_evt = rdmc_mcgetf(fp);
890 if (i1-- > 0) fp->info.bai.rec_vers = rdmc_mcgetf(fp);
891 if (i1-- > 0) fp->info.bai.min_vers = rdmc_mcgetf(fp);
892
893 rewind(fp->fp); /* rewind the file so that one */
894 /* starts reading from the beginning */
895 fp->info.bai.mc_id = (fp->info.bai.mc) ? 'U':'D';
896
897 return 0;
898
899} /* function rhd_mc() */
900
901
902/****************************************************************************/
903/* mcgeti() reads an 4 byte int from a mc/data file */
904/****************************************************************************/
905
906long rdmc_mcgeti(mcfile *fp)
907{
908
909 long r;
910 fread(&r, sizeof(r), 1, fp->fp);
911 return (fp->info.bai.swap)? rdmc_swap(r):r;
912
913} /* function rdmc_mcgeti() */
914
915
916/****************************************************************************/
917/* mcputi() writes an 4 bytes int to a mc/data file */
918/****************************************************************************/
919
920long rdmc_mcputi(long i, mcfile *fp)
921{
922 long r;
923
924 r = (fp->info.bai.swap)? rdmc_swap(i) : i;
925 fwrite(&r, sizeof(r), 1, fp->fp);
926 return r;
927
928} /* function mcputi() */
929
930
931/****************************************************************************/
932/* mcgetf() reads a 4 byte float form a mc/data file */
933/* attention: altough mcgetf() reads a float it returns double according to */
934/* the C standard */
935/****************************************************************************/
936
937double rdmc_mcgetf(mcfile *fp)
938{
939
940 union {
941 long r;
942 float f;
943 } u;
944
945 fread(&u.r, sizeof(u.r), 1, fp->fp);
946
947 if (fp->info.bai.swap) u.r = rdmc_swap(u.r);
948 return (double)u.f;
949
950} /* function mcgetf() */
951
952
953/****************************************************************************/
954/* mcputf() writes a 4 byte float to a mc/data file */
955/* attention: altought mcputf writes a 4 byte float, it requires a double */
956/* argument because this is standard in C */
957/****************************************************************************/
958
959double rdmc_mcputf(double d, mcfile *fp)
960{
961
962 union {
963 long r;
964 float f;
965 } u;
966
967
968 u.f = d;
969 if (fp->info.bai.swap) u.r = rdmc_swap(u.r);
970 fwrite(&u.r, sizeof(u.r), 1, fp->fp);
971
972 return d;
973
974} /* function mcgetf() */
975
976
977/****************************************************************************/
978/* mcseek() seeks over a number of 4-byte-words (int, float) */
979/****************************************************************************/
980
981int rdmc_mcseek(mcfile *fp, long n)
982{
983
984 if (feof(fp->fp)) return EOF;
985
986 return fseek(fp->fp, n * sizeof(long), SEEK_CUR);
987
988} /* function mcseek() */
989
990/****************************************************************************/
991/* skipevt_baikal_mc() skips to the next record (event) in a baikal mc/data file */
992/****************************************************************************/
993
994int rdmc_skipevt_baikal_mc(mcfile *fp)
995{
996
997 int nw; /* number of words in the current record */
998 int r; /* return value */
999
1000 if (fp->info.bai.fortran_recs)
1001 rdmc_mcgeti(fp);
1002
1003 nw = rdmc_mcgeti(fp);
1004 if (nw < 0) return ERRLINE;
1005
1006 r = rdmc_mcseek(fp,nw);
1007
1008 fp->info.bai.enr++; /* increase event number */
1009
1010 if (fp->info.bai.fortran_recs)
1011 rdmc_mcgeti(fp);
1012
1013 return r;
1014
1015} /* function skipevt_mc() */
1016
1017/****************************************************************************/
1018/* swap() swaps the high and low bytes of a 4 byte int */
1019/****************************************************************************/
1020
1021long rdmc_swap(long i)
1022{
1023
1024 union {
1025 unsigned char c[sizeof(long)];
1026 long l;
1027 } u; /* help unions for swapping */
1028 unsigned char r;
1029
1030 int n;
1031
1032 u.l = i;
1033
1034 for (n = 0; n < sizeof(long)/2; n++) { /* for all bytes in the first half */
1035 r = u.c[n]; /* swap the first byte */
1036 u.c[n] = u.c[sizeof(long) - n - 1]; /* with the last */
1037 u.c[sizeof(long) - n - 1] = r;
1038 } /* for n */
1039
1040 return u.l;
1041
1042} /* function swap() */
1043
1044static void rdmc_add_fitdef_baikal(array *ar){
1045
1046 array_hdef_t def_fit;
1047
1048 rdmc_jk_to_fitdef(&def_fit, 1);
1049 rdmc_add_fit_def(ar,&def_fit, 0);
1050}
1051
1052
1053#endif /* BAIKAL_BIN_F */
1054/****************************************************************************/
1055/********************************** E O F ***********************************/
1056/****************************************************************************/
1057/*
1058 This is just for EMACS:
1059 Local Variables:
1060 compile-command: "cd .. ; make -k rdmc"
1061 End:
1062*/