#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "gdecomp.h" #define MAX_SEGS 8 /* max. number of segments to take in events */ #define MAX_PARS (8*MAX_SEGS+1) /* max. number of parameters to fit */ #define MAX_INTPTS (2*MAX_SEGS) /* max. number of interaction points */ #define MAXDATASIZE 3000 #define PMODE 0644 #define GEB_TYPE_DECOMP 1 #define MAXL 500 #define SPLEN 2000 #define NPRINT 200 #define STRLEN 1024 float minr, emin, eSmear, pSmear, maxdist, vc; int maxevents, minNoInteractions, maxNoInteractions; int minNoCrystals, maxNoCrystals; int nThresholds; float e_th[120], e_0[120]; int nResolution; float enRes_0[120], enRes_1[120]; /*-----------------------------------------------------------------------------*/ int get_event_AGT (FILE * fp, int *nhits, float *egamma, int *evno, int id[MAXL], float eg[MAXL], float xx[MAXL], float yy[MAXL], float zz[MAXL]) { int i, i1, nret, type; static int firsttime = 1, nn = 0;; static char str[STRLEN], header[STRLEN]; char *p, *pc; float r1, r2, r3; /* skip the long header the first time we are called */ if (firsttime) { firsttime = 0; pc = fgets (str, STRLEN, fp); while ((p = strstr (str, "$")) == NULL) { printf ("%s", str); pc = fgets (str, STRLEN, fp); nn++; }; printf ("get_event_AGT: skipped %i lines\n", nn); /* now read the header */ pc = fgets (header, STRLEN, fp); printf ("get_event_AGT: first header: %s\n", header); }; /* interpret header */ nret = sscanf (header, "%i %f %f %f %f %i ", &type, egamma, &r1, &r2, &r3, &i1); // printf ("%i %f %f %f %f %i ", type, *egamma, r1, r2, r3, i1); assert (nret == 6); /* read this event (until we see a -1) */ i = 0; while (1) { /*read next line */ // printf ("get_event_AGT: attempting to read next line, fp=%p, str=%p\n", fp, str); // fflush (stdout); pc = fgets (str, STRLEN, fp); nn++; // printf ("get_event_AGT: %i, %p: %s", nn, pc, str); // fflush (stdout); /* done */ if (pc == NULL) { printf ("get_event_AGT: failed to read more data...\n"); fflush (stdout); return (1); }; /* interpret */ // printf ("get_event_AGT: i=%i\n", i); // fflush (stdout); nret = sscanf (str, "%i %f %f %f %f %i", &id[i], &eg[i], &xx[i], &yy[i], &zz[i], &type); // printf ("get_event_AGT: -> %i %f %f %f %f %i\n", id[i], eg[i], xx[i], yy[i], zz[i], type); // fflush (stdout); if (id[i] == -1) { /* we hit the next header */ // printf ("get_event_AGT: done reading interaction points\n"); // fflush (stdout); strcpy (header, str); *nhits = i; // printf ("get_event_AGT: return 0\n"); // fflush (stdout); return (0); } else { /* this was just an interaction point */ i++; // printf ("get_event_AGT: i now %i\n", i); // fflush (stdout); }; }; return (0); } /*-----------------------------------------------------------------------------*/ int get_event_GT (FILE * fp, int *nhits, float *egamma, int *evno, int id[MAXL], float eg[MAXL], float xx[MAXL], float yy[MAXL], float zz[MAXL]) { /* declarations */ int nret, i; /* get header */ nret = fscanf (fp, "$ %i %f %i\n", nhits, egamma, evno); // printf ("nret=%i, nhits=%3i, e=%9.2f, evno= %i\n", nret, *nhits, *egamma, *evno); if (nret != 3) return (nret); /* read event */ for (i = 0; i < *nhits; i++) { nret = fscanf (fp, "%i %f %f %f %f\n", &id[i], &eg[i], &xx[i], &yy[i], &zz[i]); if (nret != 5) { printf ("data read error\n"); exit (1); }; // printf ("id=%i e=%9.2f x=%8.0f y=%8.0f z=%8.0f\n", id[i], eg[i], xx[i], yy[i], zz[i]); }; return (0); } /*-----------------------------------------------------------------------------*/ float ranGauss (void) { double x1, x2, w, val; static int nn = 0; static float y1, y2; if (nn == 0) { do { // x1 = 2.0 * (double) rand () / RAND_MAX - 1.0; // x2 = 2.0 * (double) rand () / RAND_MAX - 1.0; x1 = 2.0 * drand48 () - 1.0; x2 = 2.0 * drand48 () - 1.0; w = x1 * x1 + x2 * x2; } while (w >= (double) 1.0); w = sqrt ((-2.0 * log (w)) / w); y1 = x1 * w; y2 = x2 * w; // printf("%f %f, %f %f; %f\n", x1, x2, y1, y2, w); nn = 1; return ((float) y1); } else { nn = 0; return ((float) y2); } } /*--------------------------------------------------------*/ void CheckNoArgs (int required, int actual, char *str) { if (required < actual) { printf ("argument problem with chat option\n"); printf ("--> %s\n", str); printf ("required # arguments: %i\n", required - 1); printf ("actual # arguments: %i\n", actual - 1); printf ("Please fix and try again, quitting...\n"); exit (1); }; } /*--------------------------------------------------------*/ int readChatFile (char *name) { /* declarations */ FILE *fp, *fp1; char *p, *pc, str[STRLEN] = { '0' }, str1[STRLEN] = { '0'}; char str2[STRLEN] = { '0' }; int echo = 0, nret = 0, nni = 0, nn = 0; int i1, i; float r2, r3; int str_decomp (char *, int, int *); int detID; float et, e0; float res0, res1; /* open chat file */ if ((fp = fopen (name, "r")) == NULL) { printf ("error: could not open chat file: <%s>\n", name); exit (0); }; printf ("chat file: <%s> open\n", name); printf ("\n"); fflush (stdout); /* read content and act */ pc = fgets (str, STRLEN, fp); while (pc != NULL) { if (echo) printf ("chat->%s", str); fflush (stdout); /* attemp to interpret the line */ if ((p = strstr (str, "echo")) != NULL) { echo = 1; if (echo) printf ("will echo command lines\n"); fflush (stdout); } else if (str[0] == 35) { /* '#' comment line, do nothing */ nni--; /* don't count as instruction */ } else if (str[0] == 59) { /* ';' comment line, do nothing */ nni--; /* don't count as instruction */ } else if (str[0] == 10) { /* empty line, do nothing */ nni--; /* don't count as instruction */ } else if ((p = strstr (str, "minr")) != NULL) { nret = sscanf (str, "%s %f", str1, &minr); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "vc")) != NULL) { nret = sscanf (str, "%s %f", str1, &vc); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "emin")) != NULL) { nret = sscanf (str, "%s %f", str1, &emin); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "eSmear")) != NULL) { nret = sscanf (str, "%s %f", str1, &eSmear); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "pSmear")) != NULL) { nret = sscanf (str, "%s %f", str1, &pSmear); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "maxdist")) != NULL) { nret = sscanf (str, "%s %f", str1, &maxdist); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "maxevents")) != NULL) { nret = sscanf (str, "%s %i", str1, &maxevents); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "minNoInteractions")) != NULL) { nret = sscanf (str, "%s %i", str1, &minNoInteractions); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "maxNoInteractions")) != NULL) { nret = sscanf (str, "%s %i", str1, &maxNoInteractions); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "minNoCrystals")) != NULL) { nret = sscanf (str, "%s %i", str1, &minNoCrystals); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "maxNoCrystals")) != NULL) { nret = sscanf (str, "%s %i", str1, &maxNoCrystals); CheckNoArgs (nret, 2, str); } else if ((p = strstr (str, "thresholds")) != NULL) { nret = sscanf (str, "%s %i", str1, &nThresholds); CheckNoArgs (nret, 2, str); for (i = 0; i < nThresholds; i++) { nn++; /* line counter */ pc = fgets (str, STRLEN, fp); /* read thresholds */ nret = sscanf (str, "%i %f %f", &detID, &et, &e0); CheckNoArgs (nret, 3, str); e_th[detID] = et; e_0[detID] = e0; } } else if ((p = strstr (str, "crystalResPar")) != NULL) { nret = sscanf (str, "%s %i", str1, &nResolution); CheckNoArgs (nret, 2, str); for (i = 0; i < nResolution; i++) { nn++; /* line counter */ pc = fgets (str, STRLEN, fp); /* read thresholds */ nret = sscanf (str, "%i %f %f", &detID, &res0, &res1); CheckNoArgs (nret, 3, str); enRes_0[detID] = res0; enRes_1[detID] = res1; } } else { /*-----------------------------*/ /* chatscript read error point */ /*-----------------------------*/ printf ("line %2.2i in chat script, option :%s \n__not understood\n", nn, str); printf ("%i\n", str[0]); printf ("aborting\n"); fflush (stdout); exit (0); }; /* read next line in chat script */ nn++; /* line counter */ nni++; /* instruction counter */ pc = fgets (str, STRLEN, fp); }; /* done */ fclose (fp); printf ("\n"); printf ("chat file: <%s> closed\n", name); printf ("__processed %i sort instructions and %i lines\n", nni, nn); printf ("\n"); fflush (stdout); return (0); } /*-----------------------------------------------------------------------------*/ int main (int argc, char **argv) { FILE *fp; int nhits, evno, nret, evcount = 0, i1, ndet, tooFar; int i, j, ng = 0, readcount = 0; float egamma, sp_obs[SPLEN], sp_emit[SPLEN], spg[SPLEN], scatter_l[SPLEN], r1, r2, r3; float sp_obs_nosmear[SPLEN], sp_emit_nosmear[SPLEN]; float esum, esum_smeared, egnow; float eg[MAXL], xx[MAXL], yy[MAXL], zz[MAXL], rr[MAXL]; int id[MAXL], pos[MAXL]; GEBDATA *gd; CRYS_INTPTS *payload; long long int ts = 1000; off_t outData; float lx, ly, lz; unsigned int seed; int siz = 0, ilok = 0, ilnotok = 0; int minNoCrystals, maxNoCrystals; double nphoto, ncompton; float polang, doppler_factor, dp, radius, emax; FILE *fp1; char str[512]; float thresh; int wr_spe (char *, int *, float *); int printCRYS_INTPTS (FILE *, CRYS_INTPTS *, GEBDATA *); if (argc != 4) { printf ("use: %s input output chatfile\n", argv[0]); exit (0); }; /* init */ for (i = 0; i < SPLEN; i++) { sp_obs[i] = 0; sp_emit[i] = 0; sp_obs_nosmear[i] = 0; sp_emit_nosmear[i] = 0; spg[i] = 0; scatter_l[i] = 0; }; gd = (GEBDATA *) calloc (1, sizeof (GEBDATA)); i1 = sizeof (CRYS_INTPTS); i1 = 10000; /* make sure there is enough room */ // payload = (CRYS_INTPTS *) calloc (1, 1); //LR payload = (CRYS_INTPTS *) calloc (1, i1); printf ("allocated %i bytes for payload\n", i1); bzero ((void *) payload, sizeof (CRYS_INTPTS)); /* initialize random number generator */ get_a_seed (&seed); srand (seed); vc = 0; minr = 5.0; emin = 30; eSmear = 1.0; pSmear = 1.0; maxdist = 300; maxevents = 2000000000; minNoInteractions = 1; maxNoInteractions = 8; minNoCrystals = 1; maxNoCrystals = 10; nThresholds = 0; for (i = 0; i < 120; i++) { e_th[i] = 0; e_0[i] = 0; } /* open chat file and read */ readChatFile (argv[3]); /* log parameters */ printf ("vc =%f\n", vc); printf ("minr =%f\n", minr); printf ("emin =%f\n", emin); printf ("eSmear =%f\n", eSmear); printf ("pSmear =%f\n", pSmear); printf ("maxdist =%f\n", maxdist); printf ("maxevents =%i\n", maxevents); printf ("minNoInteractions=%i\n", minNoInteractions); printf ("maxNoInteractions=%i\n", maxNoInteractions); printf ("minNoCrystals =%i\n", minNoCrystals); printf ("maxNoCrystals =%i\n", maxNoCrystals); printf ("nThresholds =%i\n", nThresholds); if (nThresholds > 0) { for (i = 0; i < 120; i++) { if (e_th[i] > 0) printf (" %i %f %f\n", i, e_th[i], e_0[i]); } } printf ("crystalResPar =%i\n", nResolution); if (nResolution > 0) { for (i = 0; i < 120; i++) { if (enRes_0[i] > 0) printf (" %i %f %f\n", i, enRes_0[i], enRes_1[i]); } } #if(0) /* test gaussian random number generator */ for (i = 0; i < 10000000; i++) { r1 = 300.0 * ranGauss (); i1 = (int) (r1 + 0.5); // if (i1 ) // printf ("%i, %f\n", i1,spg[i1]); i1 += SPLEN / 2; if (i1 > 0 && i1 < SPLEN) spg[i1]++; }; i1 = SPLEN; wr_spe ("gauss.spe", &i1, spg); printf ("wrote \"gauss.spe\"\n"); if (1) exit (0); #endif /* open input file from G4 */ if ((fp = fopen (argv[1], "r")) == NULL) { printf ("could not open \"%s\" as input file\n", argv[1]); exit (1); } else printf ("\"%s\" is open as input file\n", argv[1]); /* open output file */ outData = 0; outData = open ((char *) argv[2], O_WRONLY | O_CREAT | O_TRUNC, PMODE); if (outData == 0) { printf ("could not open output data file %s, quit!\n", argv[2]); exit (1); } else printf ("\"%s\" is open as out file\n", argv[2]); #if(0) nret = 0; nn = 0; while (nret == 0) { nret = get_event_AGT (fp, &nhits, &egamma, &evno, id, eg, xx, yy, zz); nn++; printf ("nhits=%3i, e=%9.2f, evno= %i\n", nhits, egamma, evno); }; if (1) exit (0); #endif /* read until we drop */ while (1) { if (evcount < NPRINT) printf ("\n-------------\n"); /* get the next event */ #if(USGEANTFORMAT==1) nret = get_event_GT (fp, &nhits, &egamma, &evno, id, eg, xx, yy, zz); #else nret = get_event_AGT (fp, &nhits, &egamma, &evno, id, eg, xx, yy, zz); // printf ("main: nret=%i\n", nret); // fflush (stdout); #endif // if (evcount >= 2) // exit (0); /* are we done? */ if (nret != 0 || evno > maxevents) { printf ("header read error at event no %i\n", evcount); printf ("made %i attempts to read events from data file\n", readcount); printf ("nphoto =%10.0f\n", (float) nphoto); printf ("ncompton=%10.0f\n", (float) ncompton); printf ("total =%10.0f\n", (float) (ncompton + nphoto)); r1 = nphoto / (nphoto + ncompton); printf ("true P/T= %4.2f\n", r1); printf ("but not all gamma rays emitted are full energy, so...\n"); fflush (stdout); /* write out sepctrum */ i1 = SPLEN; wr_spe ("G4toMode2_obs.spe", &i1, sp_obs); printf ("wrote \"G4toMode2_obs.spe\"\n"); wr_spe ("G4toMode2_emit.spe", &i1, sp_emit); printf ("wrote \"G4toMode2_emit.spe\"\n"); wr_spe ("G4toMode2_obs_nosmear.spe", &i1, sp_obs_nosmear); printf ("wrote \"G4toMode2_obs_nosmear.spe\"\n"); wr_spe ("G4toMode2_emit_nosmear.spe", &i1, sp_emit_nosmear); printf ("wrote \"G4toMode2_emit_nosmear.spe\"\n"); i1 = SPLEN; wr_spe ("scatter_l.spe", &i1, scatter_l); printf ("wrote \"scatter_l.spe\"\n"); r1 = (float) ilok / ((float) ilok + (float) ilnotok); printf ("fraction within %9.2f is %9.3f%%\n", maxdist, 100 * r1); printf ("we binned %i gamma rays\n", ng); /* full stop */ exit (0); }; readcount++; if (nhits == 0) { if (evcount < NPRINT) printf ("skip as there were no hits\n"); } else { /* open detailed listing file */ if (evcount < NPRINT) { /* open event file */ sprintf (str, "G4toMode2_%3.3i.txt", evcount + 1); fp1 = fopen (str, "w"); if (fp1 != NULL) printf ("%s is open\n", str); else { printf ("could not open %s\n", str); exit (1); }; fprintf (fp1, "----\n----file: %s\n\n", str); }; if (evcount < NPRINT) { printf ("nhits=%3i, e=%9.2f, evno= %i\n", nhits, egamma, evno); fprintf (fp1, "nhits=%3i, e=%9.2f, evno= %i\n", nhits, egamma, evno); esum = 0; for (i = 0; i < nhits; i++) { if (i == 0) fprintf (fp1, "*"); else fprintf (fp1, " "); fprintf (fp1, "read[%2i]: id=%i e=%9.2f (%8.1f %8.1f %8.1f)\n", i, id[i], eg[i], xx[i], yy[i], zz[i]); }; } if (nhits > 0) { /* find the dopler correction factor base on */ /* first interaction point reported */ #if(1) j = 0; #else emax = 0; j = 0; for (k = 0; k < nhits; k++) if (eg[k] > emax) { emax = eg[k]; j = k; }; #endif radius = xx[j] * xx[j] + yy[j] * yy[j] + zz[j] * zz[j]; radius = sqrtf (radius); dp = (zz[j]) / radius; if (dp < -1.0) dp = -1.0; if (dp > 1.0) dp = 1.0; polang = acosf (dp); radius = 1.0 - vc * vc; doppler_factor = sqrtf (radius) / (1.0 - vc * cosf (polang)); if (evcount < NPRINT) { fprintf (fp1, "true doppler factor= %9.6f, %9.2f-->%9.2f\n", doppler_factor, egamma, egamma / doppler_factor); }; #if(0) /* list all Dopler correction options */ for (j = 0; j < nhits; j++) { // printf("j=%i\n", j); radius = xx[j] * xx[j] + yy[j] * yy[j] + zz[j] * zz[j]; radius = sqrtf (radius); dp = (zz[j]) / radius; if (dp < -1.0) dp = -1.0; if (dp > 1.0) dp = 1.0; polang = acosf (dp); radius = 1.0 - vc * vc; doppler_factor = sqrtf (radius) / (1.0 - vc * cosf (polang)); if (evcount < NPRINT) printf ("%2i: doppler factor= %9.6f, %9.2f-->%9.2f\n", j, doppler_factor, egamma, egamma / doppler_factor); }; #endif #if(1) /* write out the G4 raw data */ /* do this before we sort, so first postion */ /* is truly the first interaction point as best we know */ /* notice: no smearing of energy or positions in outout data */ bzero ((void *) payload, sizeof (CRYS_INTPTS)); gd->length = 0; gd->type = GEB_TYPE_G4SIM; gd->timestamp = ts; memcpy ((char *) payload + gd->length, (char *) &nhits, sizeof (int)); gd->length += sizeof (int); memcpy ((char *) payload + gd->length, (char *) &egamma, sizeof (float)); gd->length += sizeof (float); memcpy ((char *) payload + gd->length, (char *) &evno, sizeof (int)); gd->length += sizeof (int); memcpy ((char *) payload + gd->length, (char *) &doppler_factor, sizeof (float)); gd->length += sizeof (float); for (i = 0; i < nhits; i++) { memcpy ((char *) payload + gd->length, &id[i], sizeof (int)); gd->length += sizeof (int); memcpy ((char *) payload + gd->length, &eg[i], sizeof (float)); gd->length += sizeof (float); memcpy ((char *) payload + gd->length, &xx[i], sizeof (float)); gd->length += sizeof (float); memcpy ((char *) payload + gd->length, &yy[i], sizeof (float)); gd->length += sizeof (float); memcpy ((char *) payload + gd->length, &zz[i], sizeof (float)); gd->length += sizeof (float); }; siz = write (outData, (char *) gd, sizeof (GEBDATA)); assert (siz == sizeof (GEBDATA)); siz = write (outData, (char *) payload, gd->length); assert (siz == gd->length); if (evcount < NPRINT) fprintf (fp1, "^^^^^wrote RAW G4 simulated data out with TS=%lli\n", gd->timestamp); #endif /* resort according to detector number, */ /* which is what I need to simulate decomposed data */ for (i = 0; i < nhits; i++) for (j = i + 1; j < nhits; j++) if (id[i] > id[j]) { /*swap them */ i1 = id[i]; id[i] = id[j]; id[j] = i1; r1 = eg[i]; eg[i] = eg[j]; eg[j] = r1; r1 = xx[i]; xx[i] = xx[j]; xx[j] = r1; r1 = yy[i]; yy[i] = yy[j]; yy[j] = r1; r1 = zz[i]; zz[i] = zz[j]; zz[j] = r1; }; /* check radius */ r1 = xx[i] * xx[i]; r1 += yy[i] * yy[i]; r1 += zz[i] * zz[i]; r1 = sqrtf (r1); assert (r1 < 280.0); // printf("radius= %9.2f", r1); // printf("\n"); /* debug print */ if (evcount < NPRINT) { esum = 0; fprintf (fp1, "after reordering\n"); for (i = 0; i < nhits; i++) { fprintf (fp1, "read[%2i]: id=%i e=%9.2f x,y,z=(%8.1f %8.1f %8.1f), ", i, id[i], eg[i], xx[i], yy[i], zz[i]); esum += eg[i]; fprintf (fp1, "esum=%9.2f; ", esum); r1 = xx[i] * xx[i]; r1 += yy[i] * yy[i]; r1 += zz[i] * zz[i]; r1 = sqrtf (r1); // printf ("radius= %9.2f", r1); if (i > 0) { r1 = (xx[i] - xx[i - 1]); r2 = (yy[i] - yy[i - 1]); r3 = (zz[i] - zz[i - 1]); r1 = sqrtf (r1 * r1 + r2 * r2 + r3 * r3); fprintf (fp1, "scatter len= %9.2f", r1); }; fprintf (fp1, "\n"); }; }; /* keep record of photo to compton */ esum = 0; for (i = 0; i < nhits; i++) esum += eg[i]; if ((egamma - esum) > 3.0) ncompton++; else nphoto++; /* bin the observed energy */ i1 = (int) (esum / doppler_factor); if (i1 >= 0 && i1 < SPLEN) { sp_obs_nosmear[i1]++; }; /* bin the emitted energy */ i1 = (int) (egamma / doppler_factor); if (i1 >= 0 && i1 < SPLEN) { sp_emit_nosmear[i1]++; }; /* adapted AGATA energy smearing */ r1 = sqrtf (1.0 + 3.7 * (egamma / 1000.0)); r1 /= 2.3548; if (evcount < NPRINT) fprintf (fp1, "AGATA e mearing val=%f\n", r1); r1 *= ranGauss (); if (evcount < NPRINT) fprintf (fp1, "smearing egamma=%f with %f\n", egamma, (eSmear * r1)); egamma += (eSmear * r1); i1 = (int) (egamma / doppler_factor); if (i1 >= 0 && i1 < SPLEN) { sp_emit[i1]++; }; /* find new detector points in the list */ ndet = 0; i1 = -1; for (i = 0; i < nhits; i++) if (id[i] != i1) { /* new detector in play */ pos[ndet] = i; ndet++; i1 = id[i]; }; pos[ndet] = nhits; /* only process if within requested limits */ if (ndet >= minNoCrystals && ndet <= maxNoCrystals) { /* add interaction length information */ for (i = 0; i < ndet; i++) { lx = xx[pos[i]]; ly = yy[pos[i]]; lz = zz[pos[i]]; for (j = pos[i]; j < pos[i + 1]; j++) { r1 = (lx - xx[j]) * (lx - xx[j]); r1 += (ly - yy[j]) * (ly - yy[j]); r1 += (lz - zz[j]) * (lz - zz[j]); rr[j] = sqrtf (r1); } }; if (evcount < NPRINT) { fprintf (fp1, "we have %i detectors in play\n", ndet); for (i = 0; i < ndet; i++) { fprintf (fp1, "------"); fprintf (fp1, "range %i to %i\n", pos[i], pos[i + 1] - 1); esum = 0; for (j = pos[i]; j < pos[i + 1]; j++) { fprintf (fp1, "id=%i e=%9.2f (%8.1f,%8.1f,%8.1f) dr%8.3f ", id[j], eg[j], xx[j], yy[j], zz[j], rr[j]); esum += eg[j]; fprintf (fp1, "esum=%9.2f; ", esum); if (rr[j] > minr && eg[j] > emin) fprintf (fp1, "*\n"); else fprintf (fp1, "\n"); } }; }; /* write the interactions out in decomp format */ esum_smeared = 0; for (i = 0; i < ndet; i++) { if (evcount < NPRINT) fprintf (fp1, "processing det # %i\n", id[pos[i]]); /* prepare geb header */ bzero ((void *) payload, sizeof (CRYS_INTPTS)); gd->length = sizeof (CRYS_INTPTS); gd->type = GEB_TYPE_DECOMP; gd->timestamp = ts; /* put first interaction point in */ payload->type = (int) 0xabcd5678; payload->crystal_id = id[pos[i]]; payload->timestamp = ts; payload->tot_e = eg[pos[i]]; payload->intpts[payload->num].x = xx[pos[i]]; payload->intpts[payload->num].y = yy[pos[i]]; payload->intpts[payload->num].z = zz[pos[i]]; payload->intpts[payload->num].e = eg[pos[i]]; egnow = eg[pos[i]]; if (evcount < NPRINT) fprintf (fp1, "pos[i]=%i, started payload->num=%i\n", pos[i], payload->num); payload->num++; /* add the rest if they qualify */ for (j = pos[i] + 1; j < pos[i + 1]; j++) { if (rr[j] > minr && eg[j] > emin) { if (evcount < NPRINT) fprintf (fp1, "pos[i]=%i, add as new interaction point, payload->num=%i\n", j, payload->num); /* add as new interaction point */ payload->crystal_id = id[j]; payload->timestamp = ts; payload->tot_e += eg[j]; payload->intpts[payload->num].x = xx[j]; payload->intpts[payload->num].y = yy[j]; payload->intpts[payload->num].z = zz[j]; payload->intpts[payload->num].e = eg[j]; egnow = eg[j]; payload->num++; // assert (payload->num < MAX_INTPTS); } else { /* just add the interaction energy */ if (evcount < NPRINT) { fprintf (fp1, "pos[i]=%i, just add the interaction energy: ", j); fprintf (fp1, "minr=%f, emin=%f\n", minr, emin); }; payload->tot_e += eg[j]; payload->intpts[payload->num - 1].e += eg[j]; if (eg[j] > egnow) { if (evcount < NPRINT) fprintf (fp1, " update position %i\n", j); egnow = eg[j]; #if(1) /* I think this is the best to do; but... */ payload->intpts[payload->num - 1].x = xx[j]; payload->intpts[payload->num - 1].y = yy[j]; payload->intpts[payload->num - 1].z = zz[j]; #endif }; }; }; /* add some gaussian width to energies and postitions */ payload->tot_e = 0; for (j = 0; j < payload->num; j++) { if (payload->intpts[j].e > 10) { /* Simulate energy resolution */ if (enRes_0[payload->crystal_id] > 0) { /* Use crystal resolution parameters, if defined in the chat file ... */ r1 = enRes_0[payload->crystal_id] * sqrtf (1.0 + payload->intpts[j].e * enRes_1[payload->crystal_id]); if (evcount < NPRINT) fprintf (fp1, "Crystal %i resolution=%f\n", payload->crystal_id, r1); } else { /* ... otherwise: adapted AGATA energy smearing */ r1 = sqrtf (1.0 + 3.7 * (payload->intpts[j].e / 1000.0)); r1 /= 2.3548; if (evcount < NPRINT) fprintf (fp1, "AGATA e mearing val=%f\n", r1); } /* scale by eSmear factor from the chat file */ r1 *= ranGauss (); payload->intpts[j].e += (eSmear * r1); esum_smeared += payload->intpts[j].e; /* adapted AGATA position smearing */ r1 = 0.5 * sqrtf (0.1 / (payload->intpts[j].e / 1000.0)); r1 /= 2.3548; r1 *= pSmear; payload->intpts[j].x += (r1 * ranGauss ()); payload->intpts[j].y += (r1 * ranGauss ()); payload->intpts[j].z += (r1 * ranGauss ()); }; payload->tot_e += payload->intpts[j].e; }; /* restrict on how far the interaction points can be */ /* (we cannot track across the array at the moment) */ tooFar = 0; for (j = 1; j < payload->num; j++) { r1 = payload->intpts[j].x - payload->intpts[j - 1].x; r2 = payload->intpts[j].y - payload->intpts[j - 1].y; r3 = payload->intpts[j].z - payload->intpts[j - 1].z; r1 = r1 * r1 + r2 * r2 + r3 * r3; r1 = sqrtf (r1); if (r1 > maxdist) tooFar = 1; i1 = (int) (r1 + 0.5); if (i1 >= 0 && i1 < SPLEN) scatter_l[i1]++; }; /* debug print */ if (evcount < NPRINT) printCRYS_INTPTS (fp1, payload, gd); /* write this hit out. */ if (tooFar == 0) { if (payload->num >= minNoInteractions && payload->num <= maxNoInteractions) { ilok++; /* simulate central-contact thresholds, if defined */ if (e_th[payload->crystal_id] > 0) thresh = 0.5 * (1.0 + tanh ((payload->tot_e - e_th[payload->crystal_id]) / e_0[payload->crystal_id])); else thresh = 1.0; if (drand48 () < thresh) { siz = write (outData, (char *) gd, sizeof (GEBDATA)); assert (siz == sizeof (GEBDATA)); siz = write (outData, (char *) payload, gd->length); assert (siz == gd->length); if (evcount < NPRINT) fprintf (fp1, "^^^^^wrote %i interaction points out for this detector, ts=%lli\n", payload->num, gd->timestamp); } else { if (evcount < NPRINT) fprintf (fp1, "^^^^^threshold cut out %i interaction points out for this detector, ts=%lli\n", payload->num, gd->timestamp); } } else ilnotok++; }; /* smeared esum */ r1 = sqrtf (1.0 + 3.7 * (esum_smeared / 1000.0)); if (evcount < NPRINT) fprintf (fp1, "e smearing val=%f\n", r1); r1 /= 2.3548; r1 *= ranGauss (); if (evcount < NPRINT) fprintf (fp1, "smearing esum=%f with %f\n", esum, (eSmear * r1)); esum_smeared += (eSmear * r1); i1 = (int) (esum_smeared / doppler_factor); if (i1 >= 0 && i1 < SPLEN) { ng++; sp_obs[i1]++; }; }; }; /* next timestamp */ ts += 500; evcount++; if (evcount < NPRINT) fclose (fp1); #if(0) /* stop for debugging */ if ((int) ts > 33500) if (1) exit (0); #endif }; }; }; /* done */ return (0); };