In order to operate with rdmc data structures one needs a
basic set of functions, which initialize, reset, copy, add and delete
certain information. These functions are described in the following.
Please note that functions may be added. Please refer
to rdmc.h and the references therein for more details.
(The source code is the reference manual :-).
The following functions initialize a rdmc structure
void rdmc_init_array(array *ar); void rdmc_init_array_hdef(array_hdef_t *def); void rdmc_init_mevt(mevt *ev); void rdmc_init_mtrack(mtrack *tr); void rdmc_init_mhit(mhit *h); void rdmc_init_mhit_stat(mhit_stat_t *s); void rdmc_init_mevt_uses(mevt_uses_t *u); void rdmc_init_mevt_special(mevt_special_t *s, int maxtoken);
The above functions generally set all values to zero or defaults.
Counters are set to 0, dynamic fields to NULL.
The memory of nested dynamic structures is not freed!
If a object of higher hierarchy (e.g. rdmc_init_mhit()) is initialized
the lower level structures are also initialized (thus no need
to call e.g. rdmc_init_mhit_stat()).
The structure mevt_special_t contains static fields val
a counter nval. Using only e.g. maxtoken=0 in
rdmc_init_mevt_special() will only reset the counter to 0.
If maxtoken >0 the values of val[0]..val[maxtoken-1]
are set to RDMC_NA.
The following functions free previously allocated memory before initializing the structures and can thus be used only for already initialized objects.
void rdmc_clear_array(array *ar); void rdmc_clear_mevt(mevt *ev); void rdmc_clear_mhit(mhit *h); void rdmc_clear_mtrack(mtrack *t); void rdmc_clear_mevt_special(mevt_special_t *s, int maxtoken);
If a structure contains no dynamically allocated substructures
an object could be easily copied by the standard C function
memcpy().
However nested substructures are subject of changes and it is safer
to use the supported rdmc functions:
void rdmc_cp_array(array *out, const array *in); void rdmc_cp_hdef(array_hdef_t *out, array_hdef_t *in); void rdmc_channel_cp(array *out, int out_i, const array *in , int in_i); int rdmc_cp_mevt(mevt *out, mevt *in); /* return 0 on success */ void rdmc_cp_mtrack(mtrack *out, mtrack *in); void rdmc_cp_mhit(mhit *out, mhit *in); void rdmc_cp_mevt_special(mevt_special_t *out, mevt_special_t *in); void rdmc_cp_mevt_nuses(mevt_uses_t *out, mevt_uses_t *in, int nuses);
The above functions copy the information from in to out.
The target objects have to exist (declared or allocated) and initilized,
since old memory is freed. E.g.
mevt *out = malloc(sizeof(mevt)); rdmc_init_mevt(out);
Nested substructures are automatically allocated inside the functions.
Since the array structure contains static memory for each channel
no initialization is necessary for rdmc_channel_cp, which copies
the geometry and calibration information from channel in_i
of in to channel out_i of out.
The function rdmc_cp_mevt_nuses copies nuses elements.
The target field has to be allocated to sufficient size.
A special function exists to merge two (e.g. time overlapping) hits of the same channel.
void rdmc_merge_hits(mhit *h1, mhit *h2, int ch);
The result is filled into h1. This operation
is not well defined (how to merge e.g. TOT values ?).
Therefore, before using this function the user should
look at the source code (rdmc.h) to understand what is
exactly done.
In case of merging two input streams, it is necessary to compare the two file headers, e.g. if they contain the same geometry. This is done by
int rdmc_comp_array(const array *a1, const array *a2);
The function returns 0 if the file header agree or a bit-mask
of the following codes
#define RDMC_ARRAY_COMP_OK 0x0000 #define RDMC_ARRAY_COMP_HEADER 0x01 #define RDMC_ARRAY_COMP_GEO 0x02 #define RDMC_ARRAY_COMP_OMS 0x04 #define RDMC_ARRAY_COMP_CALIB 0x08 #define RDMC_ARRAY_COMP_TRIGGER 0x10 #define RDMC_ARRAY_COMP_USER 0x20 #define RDMC_ARRAY_COMP_FIT 0x40 #define RDMC_ARRAY_COMP_STATUS 0x80 #define RDMC_ARRAY_COMP_MC 0x100
The following functions compare header definitions and return the number of differences
int rdmc_comp_array_hdef(const array_hdef_t *d1, const array_hdef_t *d2);
The structures array and mevt contain
dynamically allocated fields of substructures. In order to avoid
a complex pointer arithmetic (using malloc(), memmcpy()...)
when adding or removing objects from these fields rdmc,
provides functions to do this.
int rdmc_add_user_def(array *ar, array_hdef_t *user, int iuse); int rdmc_del_user_def(array *ar , int iuse); int rdmc_add_stat_def(array *ar, array_hdef_t *sta, int ista); int rdmc_del_stat_def(array *ar, int ista); int rdmc_add_fit_def(array *ar, array_hdef_t *fd, int ifd); int rdmc_del_fit_def(array *ar, int ifd); int rdmc_add_mcinfo_def(array *ar, array_hdef_t *mc, int imc); int rdmc_del_mcinfo_def(array *ar, int imc); int rdmc_add_trigger_def(array *ar, array_hdef_t *t, int it); int rdmc_del_trigger_def(array *ar, int it);
The above functions insert or remove a definition structure in the
data header ar
as ixxxth's element (0 <= ixxx <= ar.n_xxx).
These functions return 0 on success, 1 in case
of an error. (The above functions call
internally the functions rdmc_add_array_hdef()
and rdmc_del_array_hdef(), which are defined in
rdmc.h.
The following functions add or remove objects from an event e.g. hits or tracks:
int rdmc_add_gen(mevt *ev, mtrack *tr, int itrack); /* add MC tracks */ int rdmc_del_gen(mevt *ev , int itrack); int rdmc_add_fit(mevt *ev, mtrack *ft, mevt_special_t *fr, int ifit); int rdmc_del_fit(mevt *ev , int ifit); int rdmc_add_mhit(mevt *ev, mhit *h, int ihit); int rdmc_del_mhit(mevt *ev , int ihit); int rdmc_add_user(mevt *ev, mevt_special_t *us, int iu); int rdmc_del_user(mevt *ev , int iu); int rdmc_add_trigger(mevt *ev, mevt_special_t *tr, int it, int id); int rdmc_del_trigger(mevt *ev , int it, int id); int rdmc_add_status(mevt *ev, mevt_special_t *fr, int ifr); int rdmc_del_status(mevt *ev , int ifr); int rdmc_add_mcinfo(mevt *ev, mevt_special_t *mci, int imci); int rdmc_del_mcinfo(mevt *ev , int imci); int rdmc_add_trig_uses(mevt *ev, mevt_uses_t *tu, int ituse); int rdmc_del_trig_uses(mevt *ev , int ituse); int rdmc_add_fit_uses(mevt *ev, mevt_uses_t *fu, int ifuse); int rdmc_del_fit_uses(mevt *ev , int ifuse);
They return 0 on success, 1 on failure.
The functions for the trigger automaticaly update the two bitmask entries
mevt.trigger and mevt.trigid. Therfore the functions need
as additional parameter the trigger id defined in the header.
e.g.:
array ar; int idef = rdmc_get_hdef_tag(&(ar.def_trig), ar.n_trigger, tag_s); id = ar.def_trig[idef].id;
The following functions store all currently used hits as
mevt_uses_t
elements with the given fit/trig number
int rdmc_create_fit_uses(mevt *e, int ifit); int rdmc_create_trig_uses(mevt *e, int itrig);
In order to remove all mevt_uses_t elements
which refer to a certain fit, trigger or hit following functions
can be used
int rdmc_del_ifit_uses(mevt *ev, int ifit); int rdmc_del_itrig_uses(mevt *ev, int itrig); int rdmc_del_ihit_uses(mevt *ev, int hitid);
All above functions refer to generic driver functions,
e.g. rdmc_add_mevt_mtrack(), which can be found in rdmc.h.
However we suggest to use the above (higher level) functions.
The following functions remove a hit or a fit in the mhit
structure of an event.
In addition to int rdmc_del_mhit() it also deletes
corresponding hit uses entries, if they exist
void rdmc_remove_hit(mevt *e, int ihit); void rdmc_remove_fit(mevt *e, int ifit);
The following functions removes all hits, which are not from the leading muon
int rdmc_remove_groups(mevt *e, int invert);
If invert=0, only those hits survive, which track reference
mhit.mt points to the muon, which creates most hits.
If invert=1 these hits are deleted and all others survive.
It has to be noted, that all functions that remove or add hits
do NOT update the global event informations like mevt.nch.
This has to be done by hand, see section Updating event information.
Before adding an object it has to be properly initialized (section Initialization of objects), and filled with information, which is up to the user. However certain objects need a unique id or tag. Unique values for a hit and a MC track are provided by the following functions
int rdmc_unique_mhit_id(const mevt *ev); /*return value for mevt.h[].id */ int rdmc_unique_gen_id(const mevt *ev); /*return value for mevt.gen[].tag */
Note that fits are indexed by the number in the header definitions and thus do not need a unique id. Unique tags and ids for header definitions can be generated via
int rdmc_unique_trigger_id(char *tag, const array *ar, char *tag_root); int rdmc_unique_stat_id(char *tag, const array *ar, char *tag_root); int rdmc_unique_user_id(char *tag, const array *ar, char *tag_root);
these functions return a unique id and fill a unique string
into tag on basis of tag_root.
First tag_root itself is attempted, if this is already existing
a number is appended until tag is unique.
tag has to be large enough (RDMC_MAXTOKENLENGTH)
tag_root should be sufficiently smaller.
The first unique id attempted is the number of existing objects +1.
This number is then incremented until it is unique.
The functions return RDMC_NA in case of an error.
All functions themselves call more generic driver functions
(rdmc_unique_hdef_id, rdmc_unique_hdef_tag) which are
documented in rdmc.h.
Again fits do not need a unique definition. Several fit
definitions may have the same tag (type), e.g. rdmc-jk
(section jk fit results)
In order to test, if a certain definition exists in the header
and to obtain an index to the definition field (e.g. array.def_user)
following functions search the header definitions for an id
or tag
int rdmc_get_hdef_tag(const array_hdef_t *hd, int ndef, char *tag); int rdmc_get_hdef_id(const array_hdef_t *hd, int ndef, int id);
If the id or tag was found, the index to the definition
(0..(ndef-1)) or else RDMC_NA is returned.
For example the following example tests if a trigger
"spase-1" is defined in a file header ar.
if (rdmc_get_hdef_tag(&(ar.def_trig), ar.n_trigger, "spase-1") != RDMC_NA)
return 1; /* found */
else
return 0; /* not found */
In order to get the index or to test if a certain tag is defined
in a array_hdef_t structure the following function can be used.
It returns the index if token is a defined
word or RDMC_NA if not found.
int rdmc_token_in_hdef(const array_hdef_t *p, const char *token);
Besides the implementation of data structures and I/O, rdmc
also implements various physics functions which operate
on rdmc objects.
The task of hit cleaning and determination of direct hit
observables are discussed in separate sections
(section Hit cleaning,section Direct Hit observables)
Following functions can be used to obtain observables related to tracks, coordinates and arrival time of Cerenkov-light
int rdmc_track_closest(double *dist, double poca_xyz[3],
const mtrack *it, const double x,
const double y, const double z);
/* calculates the distance of closest approach from track it
/* to the point with the coordinates x,y,z */
/* it returns 0 and fills the values: */
/* absolute perpendicular distance *dist, and a vector xyz[3] */
/* with the three coordinates of the vector to */
/* the clostest point on the track */
/* both dist or xyz may be NULL, if you do not want to calculate */
int rdmc_track_delta(float *delta, float *dist,
const mtrack *tr1, const mtrack *tr2);
/* calculates the angular difference (solid angle) between two tracks */
/* tr1 and tr2 and the closest distance between them */
double rdmc_art_xyz(const mtrack *tr, double x, double y, double z);
/* returns Cerenkov light arrival time to point x,y,z*/
double rdmc_art(const array *ar, const mtrack *tr, int i);
/* returns Cerenkov light arrival time to i-th channel (0...)*/
double rdmc_pnt_art(const array *ar, const mtrack *tr, int i);
/* returns arrival time at channel i from a spherical wave,*/
/* assuming tr a point source, e.g. a shower */
double rdmc_om_dist(const array *a, int ich1, int ich2);
/* returns the distance between 2 channels */
double rdmc_dist(double x1, double y1, double z1,
double x2, double y2, double z2);
/* returns distance between 2 points */
double rdmc_angtr(const mtrack *t1, const mtrack *t2);
/* rdmc_angtr() calculates the cosinus of the direction angles between two */
/* tracks: returns between 1. (parrallel) and -1. (antiparallel) */
double rdmc_tdist(const array *ar, const mtrack *it, int iom);
/* calcs the min distance between a track and an om */
double rdmc_tdist_orig(const mtrack *it);
/* calcs the min distance between a track and the origin (0,0,0) */
double rdmc_vdist(const array *ar, const mtrack *it, int iom);
/* calcs the distance between a vertex of track it and a pmt */
void rdmc_tang(const array *ar, const mtrack *it, int iom,
double *dist, double *csori);
/* calculates the perpendicular distance and the */
/* cos(angle) of the incident light for a track */
/* if dist or ang are NULL, they are not calculated */
/* distance is the same as rdmc_tdist */
void rdmc_vang(const array *ar, const mtrack *it, int iom,
double *dist, double *csaxis, double *csori);
/* calculates the distance, the angle between dist and the */
/* track axis and the cos(angle) of incident light */
/* if dist or csori or csaxis are NULL, they are not calculated */
rdmc also defines following global constants:
extern const double c_vacuum; /* m/nsec; light speed in vacuum */ extern const double rez_c_vacuum; /* 1/c [ ns/m ] */ extern const double c_water; /* m/nsec; light speed in water */ extern const double rez_c_water; /* 1/c_water [ns/m ] */ extern const double n_water; /* refraction index of water */ extern const double cer; /* cerenkov angle */ extern const double tg_cer; /* tangens of cerenkov angle */ extern const double sn_cer; /* sin cerencov */ extern const double rez_sn_cer; /* 1/sin cerencov */ extern const double cs_cer; /* cos cerencov */
rdmc provides support for sorting of hits mevt.h,
MC tracks mevt.gen and hit uses objects
mevt.trig_uses, mevt.fit_uses. The functions
themselves use the standard C library function
qsort().
Many algorithms can be programmed more efficiently, if the hits are sorted in an appropriate way. This is done by the function
int rdmc_sort_hits(mevt *ev, const array *ar, enum RDMC_SORT_T method);
Presently implemented algorithms can be selected by method
which is define in rdmc.h
enum RDMC_SORT_T{
RDMC_NO_SORT=0, /* no hit sorting */
RDMC_CH_SORT=1, /* sort hits by channels */
RDMC_TM_SORT=2, /* sort hits by hit time */
RDMC_CZ_SORT=3, /* sort hits by hit channel depth */
RDMC_ID_SORT=4 /* sort hits by mhit-id */
RDMC_MT_SORT=5 /* sort hits by mt-id */
};
The above methods allow to sort hits by their time (then channel number),
the channel number (then time), the depth (z-coordinate)
of the channel which was hit (then by channel number, then by time),
by the value of their unique (unless corrupted)
id mhit.id (then by channel number, then by time) and
by the track which produced this hit mhit.mt
(then by channel number, then by time )
A special entry in the structure mevt stores
the current sort status mevt.sort_status.
Initially upon reading it is initialized to RDMC_NO_SORT
but is being updated by rdmc_sort_hits(). In order to gain efficiency
(especially during hit cleaning, section Hit cleaning), this variable
is tested upon later calls to rdmc_sort_hits().
However, if a user modifies the hits structure by adding hits
or changing the ordering, she has to assure to reset
mevt.sort_status=RDMC_NO_SORT.
The above sorting algorithms can be applied to a data stream using
the program
cpfeil (section The utility program cpfeil).
Sorting and changing of the order of MC tracks of type mtrack
in the data stream
is less important and can be dangerous!
The following function sorts the MC-track mevt.gen[]
according to the largest
number of hits they produced, but should be used very cautious.
void rdmc_sort_generated_tracks(mevt *e, int keep_tracks, int keep_all);
The provided options keep_tracks, keep_all allow
to control, if tracks without hits should be implicitly removed
from the data stream !
If keep_all=1 all tracks are kept. If keep_all=0
and keep_tracks=1 all tracks which are of the type
secondaries
(BREMS,DELTAE,PAIRPROD,NUCL_INT,MU_PAIR,HADRONS
section Particle Ids)
and did not produce any hit are removed.
If keep_all=0 and keep_tracks=0 ALL generated
mtrack objects without any hit are removed.
Note that the last option also removes MC information
like e.g. an initial neutrino track, which produces no hits.
Furthermore the number of hits is obtained by evaluating the
value mhit.mt. This algorithm is not strictly correct and
has to be used with care because several tracks may contribute to one
hit (but only one is stored) or the information in the data stream may
be corrupted.
Removal of mtrack objects is also dangerous if later
unique track ids are assigned to new mtrack objects.
Thus the safest usage is with keep_all=1.
The above function can be applied by using the SiEGMuND
program geff.
The sorting of hit-uses information is straight forward
and can be performed by the following function, which
sorts a field mevt_uses_t use[nuse] first in respect
to the value use.useid and secondly
to use.hitid
void rdmc_sort_uses(mevt_uses_t *use, int nuse);
Normally the hit uses information is sorted unless modified by the user.
This function is automatically executed by the function
rdmc_wevt().
The time of each event is stored as "Modified Julian Day" which is the number of days, seconds and nanoseconds after the 1.1.1950, 12:00 UTC. To convert this from and to standard Unix time, you may use the functions
time_t rdmc_mjd_to_unixtime(int mjd, int secs) ; /* mjd to unix secs */ int rdmc_unixtime_to_mjd(time_t unix_time); /* unix secs to mjd days*/ int rdmc_unixtime_to_mjdsecs(time_t unix_time); /* unix secs to mjd secs */
The nanosecond part of the Modified Julian Day can be used unchanged. Note that these functions are more a quick hack and do NOT include leap seconds inserted into UTC.
Old HEP programs and libraries often use a six-digit number for
showing the date. This number is in the format YYMMDD with
YY being the last two digits of the year, MM the month and
DD the day of the month (Beware of the year 2000 :-).
To convert it from or to standard Unix time
(seconds after 1.1.1970 UTC), use the functions
time_t rdmc_o_dateconv(int date); /* a YYMMDD date to unix time */ int rdmc_o_rdateconv(time_t time); /* unix time to YYMMDD format */
The conversion back to Unix time gives the time corresponding to
12:00 of the given date.
Following functions can be used to obtain information about hits and their amplitudes
int rdmc_count_nch(const mevt *e); /* returns counts number of hit channels in an event */ int rdmc_count_nstr(const mevt *e); /* returns number of hit strings in an event */ double rdmc_count_pe(const mevt *e); /* returns total pe in this event (summing the adc of all hits) */ int rdmc_omhit(float *pesum, int iom, const mevt *event); /* calcs the sum of all ADCs of hits in OM iom ([0..nom-1]) into pesum*/ /* fills pesum and returns 0 if the om was not hit or else 1 */ double rdmc_count_ompe(const mevt *e, const array *a ,float pe_oms[], int hit_oms[]); /* calculates for each om if it was hit, */ /* and how many PE (ADCs summed over all hits) this are */ /* the arrays pe_oms[] and hit_oms[] are initlized to 0 */ /* but have to be provided by the user (size: ar->nch) */ /* it fills 1 into hit_om[i] if the channel i was hit */ /* and fills total pe of all hits in this OM into pe_oms[i] */ /* returns the total number of pe of all OM */ /* the function can be called with pe_oms or hit_oms set to NULL */ /* then the corresponding values are not calculated (efficiency) */ int rdmc_centre(const mevt *ev, const array *geom, double rc[3] , double amp_pwr); /* calculate the vector rc to the centre of gravity of the event */ /* amp_pwr is the power by which the amplitude is to be weighted */
The following functions can be used to check certain MC information
(they are used by rdmc_sort_generated_tracks(), section Sorting)
int rdmc_get_leading_muon(mevt *e); /* get the leading muon from the MC tracks */ /* This is the track which produced most hits */ int rdmc_is_secondary(mtrack *t); /* returns 1 if this is a secondary shower else 0 */ int rdmc_is_point_like(mtrack *t); int rdmc_is_track_like(mtrack *t); /* returns 1 if this is a point-like/track-like light source */
For several investigations it is important to know, if a certain hit channel has hits in neighbouring channels.
int rdmc_is_coinc(mevt *e, const array *ar, int ihit, int fold, float twin); /* coincidence */ int rdmc_is_scoinc(mevt *e, const array *ar, int ihit, int fold, float twin); / scip coincidence */
The above functions return 1 if a hit e.h[ihit]
has hits in neighbouring channels (is_coinc) or
next to nearest neighbours (is_scoinc) (scip coincidence)
or return 0, else.
The times of hits have to occur within the time window twin.
The value fold give the required value of the coincidence,
e.g. fold=3
requires a 3-er coincidence of neighbour channels.
In order to be considered as a neighbouring, the channels have to be
on the same string and their channel numbers have to be close.
A skip coincidence allows that a module within the row
is not hit. E.g. a hit pattern 10010101101
is a 2-fold coincidence but a 5-fold skip coincidence.
The following function calculates if coincident hits (time window
twin) within a geometrical radius radius occur.
int rdmc_neighbours(array *a, mevt *e, int ihit, double radius, double twin);
It returns the number of hits or -1 in case of an error.
If the hits in an event mevt ev
have been modified, we recommend to add following code to your
application program
ev.sort_status = RDMC_NOSORT, ev.nch = rdmc_count_nch(&ev); ev.nstr = rdmc_count_nstr(&ev);
However certain applications, e.g. the changing of the geometry information in the header, require more serious actions, which can be done by the following functions
int rdmc_fill_mhit_str(mevt *ev, const array *ar); /* for each hit fill the string number into mhit.str */ int rdmc_repair_mhit_id(mevt *ev); int rdmc_repair_gen_id(mevt *ev); /* In case of corrupted hit/track ids (<0 or double counts) */ /* this functions assign new ids to all hits (mhit.id) */ /* or all MC tracks (mtrack.tag) */ /* if inconsitiencies are detected: corresponding references */ /* e.g. uses, mhit.mt , mtrack.parent are deleted or set to RDMC_NA */ /* returns the number of id's which were inconsistent */ void rdmc_tau_tr(mtrack *tr); /* calcs the direction cosine from the entries tr.costh, tr.phi */ /* The entries tr.px, tr.py, tr.pz are updated */ void rdmc_tauinv_tr(mtrack *tr); /* calcs the track angles costh, phi angle from tr.px, tr.py, tr.pz */ /* the entries tr.costh, tr.phi are updated */ int rdmc_count_geo(const array *a); /* returns number of strings in an geometry-header *a */ void rdmc_make_str(array *a); /* cleans up the geometry data structure *a */ /* the number of strings (a->nstr) is recalculated */ /* and the string number for each chanel (a->str[i]) are recalculated */ void rdmc_make_clust(array *a); /* calculates the iclust number for each channel in *a */ /* the old value a->clust[i] is deleted (Baikal specific) */ int rdmc_give_clust(const array *a, int ich); /* this function calcs the value iclust for channel ich in array *a */ /* it returns the result iclust or -1 on error */ /* (Baikal specific) */
Older versions of rdmc contained a special fixed
structure for fit results. This structure is used by the reconstruction
program recoos and is booked when reading the old
SiEGMuND format (DUMAND_ASCII_F).
Note that the old SiEGMuND format contained no header information
on the number of fits. In order to read these files and
create the array.def_fit entries correctly, rdmc
scans the history information in the data stream.
If thus the history information is corrupted fits cannot be
read successfully by rdmc.
For compatibility reasons rdmc supports a standard definition for
fit results (array.def_fit and mevt.fresult)
called "rdmc-jk" (array.def_fit[i].tag="rdmc-jk").
This definition is created when reading the old format
and presently also by the SiEGMuND reconstruction program
recoos. It is supported by the SiEGMuND N-Tuple programs,
e.g. munt.
The rdmc-jk convention defines 8 words and indices to
mevt.fresult according to the old definition of
the jk structure of rdmc-1.x:
#define JK_N_TOKEN 8 /* number of FRESULT tokens */ #define JK_FITID 0 /* an id of the fit algorithm */ #define JK_RCHI2 1 /* chi2/NDF or Likelihood/nch */ #define JK_PROB 2 /* some reco result parameter */ #define JK_SIGTH 3 /* some reco result parameter */ #define JK_COVMIN 4 /* some reco result parameter */ #define JK_COVMAX 5 /* some reco result parameter */ #define JK_CUTFLAG 6 /* some reco result parameter */ #define JK_CHI2 7 /* chi2 or likelihood */
The following exported variable holds a copy of the definition
extern const array_hdef_t rdmc_jk_fit_def;
Following functions support this fits of this type:
void rdmc_jk_to_fitdef(array_hdef_t *jk_def, int id);
/* initilizes a header definition as "rdmc-jk" */
void rdmc_init_fit_jk(mevt_special_t *result, int id);
/* init rdmc-jk" default values to fresult */
int rdmc_is_fresult_jk(const array *ar, const mevt *ev, int ifit);
/* checks if ev.fresult[ifit] is a "rdmc-jk" type fresult */
/* returns 0 if yes else 0*/
int rdmc_is_fitdef_jk(const array *ar, int idef);
/* returns 1 if ar.def_fit[idef] is of jk type */
int rdmc_is_this_jk(const array_hdef_t *jk_def,
const mevt_special_t *result);
/* checks if result agees with "rdmc-jk" and if
/* the definition jk_def is of type "rdmc-jk" */
/* returns 1 if yes else 0*/
It has to be noted that the naming
of the jk structure elements is historically and should not be
regarded as a definition but rather as 8 float
place-holders into which a reconstruction
program may store information.
The conversion between human readable f2000 tokens
and (magic) rdmc numbers/ids is done via lookup tables
(section Id tables). Not only for the file output but
also for user specific debug and printing purpose (e.g. in the SiEGMuND
program rtt) it is desirable to access these tables.
The tables are exported in global structures of the type:
typedef struct {
int id;
char *name;
} rdmc_idtable_t;
Each table is a linear field of elements of the type rdmc_idtable_t.
The end of the list can be tagged by a zero entry (id=RDMC_NA
and name=NULL).
Following tables can be accessed:
/* some tables defined in rdmc_local.c */ extern const rdmc_idtable_t rdmc_pmt_idtable[]; extern const rdmc_idtable_t rdmc_sphere_idtable[]; extern const rdmc_idtable_t rdmc_datatrans_idtable[]; extern const rdmc_idtable_t rdmc_detector_idtable[]; extern const rdmc_idtable_t rdmc_particle_idtable[];
Following functions allow the conversion forth and back:
/* convert OMs and PMT types to f2000 names and vice versa */ int rdmc_amanda_iomid(char *str, int ar_id); int rdmc_amanda_ipmtid(char *str); char *rdmc_amanda_spmtid(int type); char *rdmc_amanda_somid(int id); /* convert rdmc particle ids to f2000 names and vice versa */ int rdmc_amanda_ipartid(char *name); char *rdmc_amanda_spartid(int id); /* convert rdmc detector ids to f2000 names and vice versa */ char *rdmc_amanda_sdet(int geo_id); int rdmc_amanda_idet(char *detnam);
The functions return the rdmc id or a string with
the f2000 name, respectively. In case the element was not
found, the first element in the table is returned as default.
The following function do not fit into any of the above sections they are mentioned here for completeness
void rdmc_add_timeoffset(mevt *e, float offset);
/* adds a time offset to all time related values of an event */
/* this are presently all the hit times, and the fit/gen track times */
int rdmc_is_f2000_supported(int major, int minor);
/* returns 1 if the f2000.major.minor version is supported */
char *rdmc_which_format(mcfile *fp);
/* returns string with the format */
long rdmc_nint(double a);
/* converts a double to the nearest integer */
double rdmc_lgamma(double xx);
/* return the log of the gamma function */
char *rdmc_poem(void);
/* some fortune like story-teller */
The hit cleaning module in rdmc supports a various set
of algorithms. It also provides functions for:
a standard set of command line option
parsing, online help and the correct processing of cleaning tasks.
The functions are divided into a high level interface of
functions which allow to book, parse and print hit cleaning
tasks and low level functions which do the work.
We suggest the use of the high level interface, which
usual scheme is (see the SiEGMuND program soff.c as example)
is illustrated by following meta-code:
#include <unistd.h> /* getopt */
const char *program_name "myprog"
main(int argc, char **argv){
....
/* print all options as help */
fprintf(stderr," -y opt Hit reduction/ event cleaning options \n"
" opt:\n");
rdmc_usage_cleaning();
.....
/* parse options from the command line */
while ((c = getopt(argc, argv, COMMAND_SWITCHES)) != EOF)
switch (c) {
...
case 'y': /* select -y for Hit cleaning options */
if (rdmc_parse_cleaning(optarg, program_name)) /* parse the options */
errflg++;
break;
...
}
....
/* check if cleaning is to be excecuted and print debug on what is done */
if (rdmc_clean_hit.modify_mhit) /* global varibale defined */
rdmc_print_cleaning(program_name); /* in rdmc_clean_hit*/
....
/* loop all event */ {
....
/* now do cleaning before you work with mevt */
rdmc_exec_cleaning(&e,&a); /* take from global rdmc variable */
....
}
}
You may save the original event using rdmc_cp_mevt()
... and recover it later (e.g. before writing the event).
The above example would implement the following help output
-y opt
<opt>:
r=<first>:<last> reduce hits before fit
<first>= number of early hits to be removed
<last>= number of late hits to be removed
if a value is <0, take percentage of total hits
R=<time1>:<time2> reduce hits before <time1>
and after <time2> (time window)
c=<fold>:<win> use only <val>-fold coincident hits
within a window of <win> ns
C=<fold>:<win> use only <val>-fold SKIP-coincident hits
within a window of <win> ns
a=<ampl>:<amph>:<ch1>:<ch2> use only hits with
<ampl> < amp < amph> for channels
<ch1> <= ch <= <ch2> (default = all)
b=<totl>:<toth>:<ch1>:<ch2> use only hits with
<totl> < tot < toth> for channels
<ch1> <= ch <= <ch2> (default = all)
i=<time> remove hits isolated more than <time>ns
I=<radius>:<time>:<num> remove space isolated hits
more than <num> coincident (<time> [ns] window) hits
should occur within the <radius> [m]
A=<n> remove additional (subsequent) hits
after the <n>th hit in a channel
(e.g. <n>=1 keeps only first LE)
N=<n> remove hits in channel <i>
u=<i> remove hits not used by trigger <i>
<i> is the ith TRIG_DEF (1..ntrigger_def)
This evaluetes the USES lines.
WARNING: If no trigger and/or no USES lines are present,
all hits are deleted !
U=<i> remove hits not used by fit <i>
<i> is the ith FIT_DEF (1..ntrigger_def)
This evaluates the USES lines.
WARNING: If no trigger and/or no USES lines are present,
all hits are deleted !
WARNING: The cuts are processed in their order on the command line
(i.e. they are booked by rdmc_parse_cleaning)
All SiEGMuND programs, which are able to do hit cleaning (presently:
soff, deff, muff, munt, recoos), use these
functions and are thus based on exactly the same physics algorithms.
rdmc provides following interface functions:
void rdmc_usage_cleaning(void);
stderr, which options and algorithms
are implemented. (see e.g. soff -h for a printout).
rdmc_parse_cleaning(char *optarg, const char *progname);
rdmc_usage_cleaning() and books the internal data
structures (task list). Most of the algorithms depend on the ordering
of tasks (e.g. it makes a difference if you first
remove hits of small TOT and later remove channels with
no neighbouring hit or if you do it vice versa).
The hit cleaning tasks will be performed later exactly in the same
order as rdmc_parse_cleaning() was called
to book the different tasks.
void rdmc_print_cleaning(const char *progname);
stderr what cleaning
tasks have been selected.
int rdmc_exec_cleaning(mevt *e, array *a);
mevt.nch and mevt.nstr)
void rdmc_clear_cleaning_list();
int rdmc_clean_backup_hits(mevt *e);
rdmc_exec_cleaning().
rdmc_clean_restore_hits(mevt *e);
In order to check, if any task list has been booked at all,
rdmc creates a global variable rdmc_clean_hit
which is set to 1, if any cleaning tasks are selected (0 else):
typedef struct {
int modify_mhit; /* flag to indicate modification of mhit during fit*/
} rdmc_clean_hit_t;
extern rdmc_clean_hit_t rdmc_clean_hit;
Following routines do the dirty work and are given for completeness
(They do NOT update mevt.nch and mevt.nstr)
int rdmc_rm_h(array *a, mevt *e, int hit_nr, int nhits);/*remove hit: hit_nr*/
int rdmc_rm_rnd_h(array *a, mevt *e); /* remove one random hit */
int rdmc_rm_rnd_nh(array *a, mevt *e, int n); /* remove n random hits */
int rdmc_rm_earliest_h(array *a, mevt *e, int num);/*remove num earliest hit*/
int rdmc_rm_latest_h(array *a, mevt *e, int num); /* remove num latest hit */
int rdmc_rm_interval_h(array *a, mevt *e, float before, float after);
/* remove hits outside time window*/
int rdmc_rm_amp_h(array *a, mevt *e, float amp_low, float amp_high);
/*remove hits: amp_low > amp > amp_high */
int rdmc_rm_tot_h(array *a, mevt *e, float tot_low, float tot_high);
/*remove hits tot_low > tot > tot_high */
int rdmc_rm_nocoinc_h(array *a, mevt *e, int coinc, float twin);
/*remove hits with less <fold> neighbour hit channel (time window: twin) */
int rdmc_rm_snocoinc_h(array *a, mevt *e , int coinc, float twin);
/*remove hits with less <fold> SKIP-neighbour hit channels within twin */
int rdmc_rm_isolate_h(array *a, mevt *e, float window);
/* remove hits if no other hit occurs within time: window*/
int rdmc_rm_local_h(array *a, mevt *e, double radius, double twin, int nhits);
/* remove hits with less then nhits other(neighbour) hits */
/* within distance radius and time twin */
int rdmc_rm_additional_h(array *a, mevt *e, int from);
/* remove hits more than the <from> hits in each channel */
int rdmc_rm_channel_h(array *a, mevt *e, int ich);
/*remove all hits in channel ich*/
int rdmc_rm_uses_h(array *a, mevt *e, int itrig);
/* remove all hits which are not marked in a USES line by trigger itrig*/
int rdmc_rm_fuses_h(array *a, mevt *e, int ifit);
/* remove all hits which are not marked in a USES line by fit iffit*/
The functions return the number of hits, which were removed.
Similar to the hit cleaning the rdmc provides a higher
level interface which calculates the time residuals in respect
to the Cerenkov model and then calculates various hit observables
based on these times.
The basic idea is to provide a function similar
to the standard C function scanf(), which
receives a variable number of arguments and fills the values
according to a format string.
The function is:
double rdmc_get_direct_hits(const array *a , const mtrack *tr , mevt *e , char *fmt , ...);
It receives the file header a, the track tr
for which the hits are to be calculated and the event e,
which contains the hits. What is to be done is passed
by a format string fmt, which is followed by pointers to the
variables, which are to be filled.
The syntax of one format string element is:
t-:t+;f
The values t- and t+ (with a colon in between)
give the allowed time range
for the time residual. A hit which time residual is between
is by convention in the following called 'direct-hit', even
if the time window is large.
Separated by a ; f is a
character code
(similar to scanf()) which selects a specific observable.
The time residual is defined by the difference of the hit time
minus the time expected for the hit, calculated from the track
parameters, the position of the channel and the Cerenkov model.
Possible values for f and the type of the parameter are:
c
s
x
y
z
d
D
R
L
M
T B N n
Note that by changing the time window arbitrarily any time window can be selected. This is illustrated by the following example:
#define DIRECT_PARSE "-1000.:-5.;c -5.:20.;c -5.:20.;s 75.:1000.;c -5.:20.;L"
int n_early, n_direct, s_direct, n_late;
float z_direct;
float gdir = get_direct_hits( &a, &tr, &e, DIRECT_PARSE
, &n_early, &n_direct, &s_direct
, &n_late, &z_direct );
which calculates the number of "early" hits n_early
with residuals between -1000ns..-5ns,
the number of "direct" hits n_direct
with residuals between -5ns..20ns,
the number of "direct" strings s_direct
with residuals between -5ns..20ns,
the number of "late" hits n_late
with residuals between 75ns..1000ns
and the project length
of "direct" hits z_direct
with residuals between -5ns..20ns.
The function returns a value called gdir which is
calculated by summing over all hits
gdir += 1./sqrt(2*pi*(3*8ns)^2) * exp( -0.5* (residual/(3*8ns))^2)
The SiEGMuND program muff accepts codes for get_direct_hits()
from the command line and filters events according to the results.
Note that the ; has to be escaped on a shell command line.
Pandel functions are especially important for AMANDA. Though not
directly connected to rdmc datastructures, it is convenient to
modularize them and keep them in a library.
There are a lot of functions to be used, the following
pt
pp
td
ph, pnh
lg
norm
int, diff
mpe
psa
patched
Following functions are implemented:
/* now the functions */
double rdmc_pt_td(double delay, double perp_dist, double cs_ori);
double rdmc_pt_lgtd(double delay, double perp_dist, double cs_ori);
double rdmc_pt_td_norm(double perp_dist, double cs_ori);
double rdmc_pt_lgtd_norm(double perp_dist, double cs_ori);
double rdmc_pt_td_int(double delay, double perp_dist, double cs_ori);
double rdmc_pt_lgtd_int(double delay, double perp_dist, double cs_ori);
double rdmc_pt_td_diff(double delay, double perp_dist, double cs_ori);
double rdmc_pt_td_mpe(double delay, double perp_dist, double cs_ori
,double pe);
double rdmc_pt_lgtd_mpe(double delay, double perp_dist
,double cs_ori, double pe);
double rdmc_pt_td_mpe_int(double delay,double perp_dist
,double cs_ori,double pe);
double rdmc_pt_lgtd_mpe_int(double delay,double perp_dist
,double cs_ori,double pe);
double rdmc_pt_td_mpe_diff(double delay,double perp_dist
,double cs_ori,double pe);
double rdmc_pt_td_psa(double delay, double perp_dist
,double cs_ori, double mean_pe);
double rdmc_pt_lgtd_psa(double delay, double perp_dist
,double cs_ori, double mean_pe);
double rdmc_pt_td_psa_int(double delay, double perp_dist
,double cs_ori, double mean_pe);
double rdmc_pt_lgtd_psa_int(double delay, double perp_dist
,double cs_ori, double mean_pe);
double rdmc_pt_td_psa_diff(double delay, double perp_dist
,double cs_ori, double mean_pe);
double rdmc_pt_td_patched(double delay, double perp_dist, double cs_ori);
double rdmc_pt_lgtd_patched(double delay, double perp_dist, double cs_ori);
double rdmc_pt_td_mpe_patched(double delay, double perp_dist
, double cs_ori, double pe);
double rdmc_pt_lgtd_mpe_patched(double delay, double perp_dist
, double cs_ori, double pe);
double rdmc_pt_td_psa_patched(double delay, double perp_dist
, double cs_ori, double mean_pe);
double rdmc_pt_lgtd_psa_patched(double delay, double perp_dist
, double cs_ori, double mean_pe);
double rdmc_pt_ph(double perp_dist, double cs_ori, double energy
, double sensit);
double rdmc_pt_pnh(double perp_dist, double cs_ori, double energy
, double sensit);
double rdmc_pp_td(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_lgtd(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_td_norm(double perp_dist, double cs_ori, double cs_axis);
double rdmc_pp_lgtd_norm(double perp_dist, double cs_ori, double cs_axis);
double rdmc_pp_td_int(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_lgtd_int(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_td_diff(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_td_mpe(double delay,double perp_dist,double cs_ori
, double cs_axis,double pe);
double rdmc_pp_lgtd_mpe(double delay,double perp_dist,double cs_ori
, double cs_axis,double pe);
double rdmc_pp_td_mpe_int(double delay,double perp_dist,double cs_ori
, double cs_axis,double pe);
double rdmc_pp_lgtd_mpe_int(double delay,double perp_dist,double cs_ori
, double cs_axis,double pe);
double rdmc_pp_td_mpe_diff(double delay,double perp_dist,double cs_ori
, double cs_axis,double pe);
double rdmc_pp_td_psa(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_lgtd_psa(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_td_psa_int(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_lgtd_psa_int(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_td_psa_diff(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_td_patched(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_lgtd_patched(double delay, double perp_dist, double cs_ori
, double cs_axis);
double rdmc_pp_td_mpe_patched(double delay, double perp_dist, double cs_ori
, double cs_axis, double pe);
double rdmc_pp_lgtd_mpe_patched(double delay, double perp_dist, double cs_ori
, double cs_axis, double pe);
double rdmc_pp_td_psa_patched(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
double rdmc_pp_lgtd_psa_patched(double delay, double perp_dist, double cs_ori
, double cs_axis, double mean_pe);
Upon start the phyysics parameters are initilized to default values. The can be overwritten by externally calling
void rdmc_td_init(rdmc_pandel_par_t *pandel_par);
The passed structure is defined:
typedef struct {
double td_tau; /* tau for tracks */
double td_lam; /* lambda for tracks */
double td_att; /* Absorption lenght */
double ps_tau; /* tau for tracks */
double ps_lam; /* lambda for tracks */
double td_sigma; /* pmt jitter -> patched functions */
double td_dist_p1; /* scale for the distance */
double td_dist_p0_cs0; /* const for distance ped (P0) */
double td_dist_p0_cs1; /* const for distance ped (P1*cs_ori) */
double td_dist_p0_cs2; /* const for distance ped (P2*cs_ori^2) */
double ph_lam; /* lambda for Phit function track */
double ph_eps_pe;
double ph_eps_ori;
} rdmc_pandel_par_t;
The structure may contain specific values, for certain parameters.
If a value is set to RDMC_SPECIAL_NA (section Special values),
the current default is not changed. The library als defines
several data sets of predefined optics:
extern const rdmc_pandel_par_t rdmc_pandel_par_h0, /* No irregulaarities in the hole */ rdmc_pandel_par_h1, /* hole irregularities scat=100cm */ rdmc_pandel_par_h2, /* hole irregularities scat=50cm */ rdmc_pandel_par_h3, /* hole irregularities scat=30cm */ rdmc_pandel_par_h4; /* hole irregularities scat=10cm */
In order unify the handling of message output, rdmc provides functions,
which do message printing to stderr
void rdmc_msgprintf(char *fmt, ...);
void rdmc_errorprintf(char *fmt, ...);
void rdmc_warnprintf(char *fmt, ...);
The format string fmt is similar to the standard C
function printf().
The functions automatically append a '\n' to the output
and print as standard message the program name followed by the
UNIX process id (this is important if the same program occurs several times
with in a long pipe) and the message. In case of errors or warnings,
the string "Error:" or "Warning:" is appended
after the process id.
In order to use these functions the global char *program_name
variable has to be defined, e.g.
char *program_name = "myname"
In case rdmc returns an error of the type RDMC_ERROR_T
(see section Return values)
the following function may be called to get a more detailed error
output (For an example output see section The utility program cpfeil)
void rdmc_err_print(mcfile *fp, enum RDMC_ERROR_T ierr);
Go to the first, previous, next, last section, table of contents.