From b90bba9aa13f594ea65f90bc29334c05c2bfbda7 Mon Sep 17 00:00:00 2001
From: Hernan Asorey <asoreyh@gmail.com>
Date: Wed, 31 Dec 2014 00:12:32 -0500
Subject: [PATCH] sol.cc now fit muon hump and other features to identify muon
and other bands. At this commit, waiting for Yunior to test on one month of
data
---
sol.cc | 144 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 142 insertions(+), 2 deletions(-)
diff --git a/sol.cc b/sol.cc
index 842d0f7..bd3eed8 100644
--- a/sol.cc
+++ b/sol.cc
@@ -68,6 +68,32 @@ double itpl_chr(double *pos_avg, double *chr_avg, int i, double x) {
return (y0 + (y1 - y0) * (x - x0) / (x1 - x0));
}
+double fitHisto(int xi, int xf, int* h, double* r) {
+ double y1, x1, x2, x3, x4, xy, x2y, x, y, v1, v2, v;
+ v = v1 = v2 = x = y = y1 = x1 = x2 = x3 = x4 = xy = x2y = 0.0;
+ int n = 0;
+ for (int i = xi; i <= xf; i++) {
+ x = i;
+ y = (double) h[i];
+ x1 += x;
+ y1 += y;
+ xy += x * y;
+ x *= i;
+ x2 += x;
+ x2y += x * y;
+ x *= i;
+ x3 += x;
+ x *= i;
+ x4 += x;
+ n++;
+ }
+ v1 = x1 * x2 * x2y - n * x2y * x3 - x2 * x2 * xy + n * x4 * xy + x2 * x3 * y1 - x1 * x4 * y1;
+ v2 = x1 * x1 * x2y - n * x2 * x2y + n * x3 * xy + x2 * x2 * y1 - x1 * (x2 * xy + x3 * y1);
+ v = v1 / (2.0 * v2); // note v1 = -a1 and v2 = a2, so we have v = -b / 2a !
+ *r = y1;
+ return v;
+}
+
void Usage(char *prog) {
cout << endl << PROJECT << " " << VERSION << endl;
cout << endl << "This is the Space Weather dedicated analysis (L1 -> L2SW)" << endl<< endl;
@@ -180,7 +206,7 @@ int main (int argc, char *argv[]) {
/* preparing outputs */
char *ifile2;
- char ifile[256];
+ char ifile[256], ifile0[256];
ifile2=ifiname;
if (strrchr(ifile2,'/')!=NULL)
ifile2=strrchr(ifile2,'/')+1; // remove dirs if present
@@ -191,6 +217,7 @@ int main (int argc, char *argv[]) {
if (strrchr(ifile,'.')!=NULL)
if (strncmp(strrchr(ifile,'.'),".sol",4)==0) // remove .sol if present
*(strrchr(ifile,'.'))='\0';
+ strcpy(ifile0,ifile);
if (strncmp(ifile,"l1_",3)==0) {
*(strchr(ifile,'1'))='2'; // change l1 by l2
}
@@ -199,7 +226,6 @@ int main (int argc, char *argv[]) {
Usage(argv[0]);
}
-
/* end of preparing outputs */
/* preparing for work */
@@ -275,7 +301,119 @@ int main (int argc, char *argv[]) {
char line[16384];
double fluxband[nbands];
int currentpos=0, pos=0, i=0, is=0, j=0,k=0;
+ /// New analysis done by Y. Perez and H. Asorey
+ ///// Obtain muon charge and other features from cal histogram
+ ///// and use them for determining bands of interest
+ ///// underway: mar dic 30 23:38:56 COT 2014
+ /* now, we have to read calibration (1h) histogram from ifile.cal, so:
+ * trying to open histogram file
+ */
+ snprintf(nfi,256,"%s.cal",ifile0);
+ if ((cal = fopen(nfi,"r"))==NULL) {
+ fprintf(stderr,"\nError: Failed to open calibration histogram file (%s). Abort.\n",nfi);
+ exit(1);
+ }
+ /* read the histogram */
+ fprintf(stderr,"Reading calibration histogram %s\n",nfi);
+ int p[3], c[3], charge[maxcharge];
+ int cline=0, ncharge=0;
+ while(fgets(line, 16384,cal)) {
+ cline++;
+ if (line[0] == '#') // comments? meta-data? for now, just skipping
+ continue;
+ if (cline < 1028) {
+ do { // read
+ check=sscanf(line, "%d %d %d %d %d %d\n", &c[0], &c[1], &c[2], &p[0], &p[1], &p[2]);
+ } while (check != 6);
+ }
+ else {
+ do { // read
+ check=sscanf(line, "%d %d %d\n", &c[0], &c[1], &c[2]);
+ } while (check != 3);
+ }
+ if (ncharge>maxcharge)
+ fprintf(stderr,"Warning: Too long charge histogram. It has %d lines (expected lines = %d)\n", ncharge, maxcharge);
+ charge[ncharge] = c[channel];
+ ncharge++;
+ }
+ fclose(cal);
+ /* now, charge[] has the charge histogram for the analyzed channel
+ * average it!, but, graded
+ */
+ int l1 = 200; /* Yunior used this cut . If charge < l1, ds=6; else ds = 25
+ * I put ds=5 as a multiple of l1
+ */
+ int sum = 0, cnt_avg=0;
+ double chr_avg[maxcharge], pos_avg[maxcharge];
+ for (i=0; i<maxcharge; i++)
+ chr_avg[i] = pos_avg[i] = 0.;
+ int step=0;
+ /* low energy */
+ step=5;
+ for (i=0; i<l1; i+=step) {
+ sum = 0;
+ for (j=0; j<step; j++)
+ sum += charge[i+j];
+ chr_avg[cnt_avg]=sum/step;
+ pos_avg[cnt_avg]=i+step/2.;
+ cnt_avg++;
+ }
+ /* high energy */
+ step=50;
+ for (i=l1; i<maxcharge; i+=step) {
+ sum = 0;
+ for (j=0; j<step; j++)
+ sum += charge[i+j];
+ chr_avg[cnt_avg]=sum/step;
+ pos_avg[cnt_avg]=i+step/2.;
+ cnt_avg++;
+ }
+ /* Histogram are averaged. Now, derive it! */
+ cnt_avg-=2; //avoid last average point because it's usally too noise
+ double chr_der[cnt_avg];
+ for (i=0; i<cnt_avg; i++) //taking logs to magnify the effect (thanks Mauricio Suarez)
+ chr_der[i] = ((log10(chr_avg[i+1])-log10(chr_avg[i]))/(pos_avg[i+1]-pos_avg[i]));
+ // looking for 0s on histogram derivative...
+ double zeros[cnt_avg], zeros_pos[cnt_avg];
+ int cnt_zeros=0;
+ for (i=0; i<cnt_avg; i++)
+ zeros_pos[i]=zeros[i]=0.;
+ for (i=0; i<cnt_avg-1; i++) {
+ if (!sign(chr_der[i])) { // bingo... too good to be true
+ zeros_pos[cnt_zeros]=(double) i;
+ zeros[cnt_zeros]=chr_avg[i]; // charge
+ cnt_zeros++;
+ }
+ else if (sign(chr_der[i+1]) != sign(chr_der[i])) { // Bolzano's theorem!!
+ zeros_pos[cnt_zeros]=itpl_pos(pos_avg, chr_der,i);
+ zeros[cnt_zeros]=itpl_chr(pos_avg, chr_avg, i, zeros_pos[cnt_zeros]);
+ cnt_zeros++;
+ }
+ }
+ /* So, we have the histogram seeds, let's fit them */
+ double tmpk=0.;
+ double feat[3];
+ // TODO erase this block, only for debbuging of fitting procedure
+ snprintf(nfi,256,"%s.pos",ifile);
+ FILE *pof;
+ if ((pof = fopen(nfi,"w"))==NULL) {
+ fprintf(stderr,"\nError: Failed to open pos file. Abort.\n");
+ exit(1);
+ }
+ //
+
+ for (i=0; i<3; i++) {
+ feat[i] = fitHisto(int(zeros_pos[i]*0.9), int(zeros_pos[i]*1.3), charge, &tmpk);
+ // TODO erase this block, only for debbuging of fitting procedure
+ do {
+ check=fprintf(pof, "%d %.3f %.3f\n", i, zeros_pos[i], feat[i]);
+ } while (check<=0);
+ }
+ /* and done. We found the EM peak, the transition to hump, and the muon hump
+ * of the calibration histogram, in charge units =D
+ */
+ /////// HERE ENDS THE NEW BLOCK - WAITING FOR YUNIOR TESTS ON DATA TO INCLUDE
/* analyze the .sol histograms */
while(fgets(line, 16384,fi)) {
@@ -341,4 +479,6 @@ int main (int argc, char *argv[]) {
pclose(out);
else
fclose(out);
+ // TODO erase this block, only for debbuging of fitting procedure
+ fclose(pof);
}
--
GitLab