/*==> rate.c ===================================================
  calculate the rates,.. measured by scalers
  !! -> measures the TaCoS-system time !!!

  ----------------------------------------------------------- */

/*=== includes ============================================== */
#include <unistd.h> /* -> usleep*/
#include <math.h>   /* -> fabs and simu*/
#include <stdio.h>  /* -> rand */
#include <minmax.h> /* -> MAX*/

#include "rate.h"
#include "rsimu.h"
#include "dev/scalers.h"
#include "target/main/target.h"  /* for rate simu */
#include "target/main/threads.h" /* for rate simu */

/*=== defines =============================================== */
#define RATE_BG        7
#define RATE_REF_RATE  0
#define RATE_REF_FREQ  10.416667L /*MHz*/

/*=== structs =============================================== */
struct sc_data {
  struct timespec time;
  scaler_t        value[NETRATE_NUM_RATES];
};
struct sc_data_diff {
  int            time;    /* time difference in usec */
  scaler_t       value[NETRATE_NUM_RATES];
};
typedef struct {
  struct sc_data      old;
  struct sc_data      new;
  struct sc_data_diff diff;
} scaler_val_t;


/* a buffer to remember the last n readouts of the charge Integrator rates */
/* 10 Hz readout is too much for this statistics */
#define R_NUM_WRRATES    20
typedef struct {
  rate_t  rate[NP_NUM_WIRES][R_NUM_WRRATES]; /* at each position ill have the 
						(rate at this wire)/R_NUM_CIRATES in MHz*/
  int     pos;                               /* actual position in the buffer*/
  rate_t  pedestal[NP_NUM_WIRES]; /* for calculating the real IA rate at the wire*/
  double  calf[NP_NUM_WIRES];     /* calibration factor */
  int     rpos[NP_NUM_WIRES];     /* positions, where the rates can be found */
} wrate_t;

/* this structure will contain the information need to compute 
   the first derivatives of the rate-"function" */
typedef struct {
  int    pos;
  rate_t rate[NETRATE_NUM_DRSLICES];
  rate_t srateq[NETRATE_NUM_DRSLICES];
  rate_t s1, s2, s3, s4, s5;
  drate_t dr;
} drateh_t;

/*=== variables ============================================= */
/* general options for rate.c */
rate_t          *Acc;
int              Testrun;
int              Debugmode;
int              Use_ref_rate;
the_rates_t      Rate;

/* Scaler values for calculating the rates */
scaler_val_t     Scaler;
/* single wire rate calculation structure */
wrate_t          WR;
drateh_t         DR;

/*=== missing system prototypes ============================= */


/*--- init-rates() ----------------------------------------
  initializes the scalers and
  fills the rate structures with the first values
  returns: void
  ------------------------------------------------------- */
void init_rates(int testrun, int use_ref_rate, rate_t *acc, data_t *d) {
  int      i, j;

  Acc = acc;

  init_rsimu();

  Use_ref_rate       = use_ref_rate;
  Testrun            = testrun;
  Rate.integrated_ia = 0;
  Rate.ia_mean       = 0.0;
  Rate.sigma_mean    = 0.0;

  WR.pos              = 0;
  for(i=0; i<NP_NUM_WIRES; i++) {
    Rate.wrate[i]        = 0;
    Rate.wrate_nomean[i] = 0;
    for(j=0; j< R_NUM_WRRATES; j++) {
      WR.rate[i][j] = 0.0;
    }
  }

  /* CI values up to 20.7.99
  WR.calf[0] =   87000;
  WR.calf[1] =  140000;
  WR.calf[2] =   30000;
  WR.calf[3] =       0;
  WR.calf[4] =  120000;
  WR.calf[5] =   54000;
  WR.calf[6] =       0;
  WR.calf[7] =  130000;

  WR.pedestal[0]     =    -40.5;
  WR.pedestal[1]     =    -52.7;
  WR.pedestal[2]     =    -10.9;
  WR.pedestal[3]     =      0.0;
  WR.pedestal[4]     =    -64.8;
  WR.pedestal[5]     =    -27.1;
  WR.pedestal[6]     =      0.0;
  WR.pedestal[7]     =    -56.4;
  */ 


  /* implemented 28.07.99 */
  /* CI values up to 24. Jan 2000 
  WR.calf[0] =  64683;
  WR.calf[1] =  38168;
  WR.calf[2] =  39683;
  WR.calf[3] =  69109;
  WR.calf[4] =  13563;
  WR.calf[5] =  99403;
  WR.calf[6] =  17749;
  WR.calf[7] =  60240;    

  WR.pedestal[0]     =    -27.620;
  WR.pedestal[1]     =    -16.040;
  WR.pedestal[2]     =    -23.214;
  WR.pedestal[3]     =    -30.667;
  WR.pedestal[4]     =    -7.975;
  WR.pedestal[5]     =    -49.220;
  WR.pedestal[6]     =    -10.428;
  WR.pedestal[7]     =    -33.639; */


  /* implemented 24. Jan. 2000 */

  /* CI values up to 26. Jan 2000 
  WR.calf[0] =  1000000.0/96.8;
  WR.calf[1] =  1000000.0/41.8;
  WR.calf[2] =  1000000.0/46.0;
  WR.calf[3] =  1000000.0/18.2;
  WR.calf[4] =  1000000.0/8280;
  WR.calf[5] =  1000000.0/2130;
  WR.calf[6] =  1000000.0/10900;
  WR.calf[7] =  1000000.0/9800;    

  WR.pedestal[0]     =      -159.0 /    96.8;
  WR.pedestal[1]     =      -149.0 /    41.8;
  WR.pedestal[2]     =      -153.0 /    46.0;
  WR.pedestal[3]     =      -151.0 /    18.2;
  WR.pedestal[4]     =    -18890.0 /  8280.0;
  WR.pedestal[5]     =    -19270.0 /  2130.0;
  WR.pedestal[6]     =    -19670.0 / 10900.0;
  WR.pedestal[7]     =    -19770.0 /  9800.0; */

 /* implemented 26. Jan. 2000 */

  /* CI values up to 29. Jan. 2000 
  WR.calf[0] =  12800.;
  WR.calf[1] =  24300.;
  WR.calf[2] =  18500.;
  WR.calf[3] =  54600.;
  WR.calf[4] =  158.;
  WR.calf[5] =  476.;
  WR.calf[6] =  91.7;
  WR.calf[7] =  105.;    

  WR.pedestal[0]     =      -3.5 ;
  WR.pedestal[1]     =      -4.5;
  WR.pedestal[2]     =      -4.2;
  WR.pedestal[3]     =      -8.1;
  WR.pedestal[4]     =      -2.9;
  WR.pedestal[5]     =      -9.4;
  WR.pedestal[6]     =      -1.8;
  WR.pedestal[7]     =      -2.0; */

 /* implemented 29. Jan. 2000 */

  /* values up to 11. Feb. 2000 
  WR.calf[0] =  141.;
  WR.calf[1] =  20600.;
  WR.calf[2] =  16507.;
  WR.calf[3] =  24400.;
  WR.calf[4] =  158.;
  WR.calf[5] =  476.;
  WR.calf[6] =  91.7;
  WR.calf[7] =  105.;    

  WR.pedestal[0]     =      -2.7 ;
  WR.pedestal[1]     =      -4.1;
  WR.pedestal[2]     =      -3.73;
  WR.pedestal[3]     =      -3.0;
  WR.pedestal[4]     =      -2.9;
  WR.pedestal[5]     =      -9.4;
  WR.pedestal[6]     =      -1.8;
  WR.pedestal[7]     =      -2.0; */


/* implemented 11.Feb. 2000 */

  WR.calf[0] =  141.;
  WR.calf[1] =  231.06;
  WR.calf[2] =  210.73;
  WR.calf[3] =  240.34;
  WR.calf[4] =  158.;
  WR.calf[5] =  22000.;
  WR.calf[6] =  10017;
  WR.calf[7] =  15468.;    

  WR.pedestal[0]     =      -2.7 ;
  WR.pedestal[1]     =      -4.761;
  WR.pedestal[2]     =      -4.56;
  WR.pedestal[3]     =      -4.67;
  WR.pedestal[4]     =      -2.9;
  WR.pedestal[5]     =      -9.4;
  WR.pedestal[6]     =      -2.3;
  WR.pedestal[7]     =      -3.24;



  WR.rpos[0]     = 24;
  WR.rpos[1]     = 25;
  WR.rpos[2]     = 26;
  WR.rpos[3]     = 27;
  WR.rpos[4]     = 28;
  WR.rpos[5]     = 29;
  WR.rpos[6]     = 30;
  WR.rpos[7]     = 31;


  DR.pos = 0;
  for(i=0; i<NETRATE_NUM_DRSLICES; i++) {
    DR.rate[i] = DR.srateq[i] = 0;
 }

  if(!Testrun)
    init_scaler();

  /* Set the time of the scaler readout and */
  /* Read all rates from the scaler and write it to the structure */
  getclock(TIMEOFDAY, &Scaler.old.time);

  if(!Testrun) {
    for(i=0; i<NETRATE_NUM_RATES; i++)
      Scaler.old.value[i] = read_scaler(i);
  } else {
    rate_sim_setstat(3);
  }

  update_rates(d);
  usleep(100000);
}

/*----------------------------------------------------------
  ------------------------------------------------------- */
void term_rates() {
  term_rsimu();
}

/*----------------------------------------------------------
  ------------------------------------------------------- */
void calc_dr(data_t *d) {
  int i, tt;
  rate_t xds, nen;

  /* first calculate the sums */
  tt   =0;
  DR.s1=0;
  DR.s2=0;
  DR.s3=0;
  DR.s4=0;
  DR.s5=0;

  for(i=DR.pos; i>=0; i--) { 
    xds    = tt/DR.srateq[i];

    DR.s1 += 1/DR.srateq[i];
    DR.s2 += DR.rate[i]*xds;
    DR.s3 += DR.rate[i]/DR.srateq[i];
    DR.s4 += xds;
    DR.s5 += tt*xds;

    tt--;
  }  
  for(i=d->autom.dr_slices-1; i>DR.pos; i--) { 
    xds    = tt/DR.srateq[i];

    DR.s1 += 1/DR.srateq[i];
    DR.s2 += DR.rate[i]*xds;
    DR.s3 += DR.rate[i]/DR.srateq[i];
    DR.s4 += xds;
    DR.s5 += tt*xds;

    tt--;
  }  

  nen = (DR.s1*DR.s5 - DR.s4*DR.s4);
  
  DR.dr.a  = (DR.s1*DR.s2 - DR.s3*DR.s4) / nen;
  DR.dr.b  = (DR.s5*DR.s3 - DR.s4*DR.s2) / nen;
  DR.dr.sa = sqrt(DR.s1/nen);
  DR.dr.sb = sqrt(DR.s5/nen);

  /*print_scrollX("%3d>> a: %5.2g/%5.2g | b: %5.2g/%5.2g | nen: %5.2g | s: %5.2g - %5.2g - %5.2g - %5.2g - %5.2g",
		DR.pos, DR.dr.a, DR.dr.sa, DR.dr.b, DR.dr.sb, nen,
		DR.s1, DR.s2, DR.s3, DR.s4, DR.s5);*/

}

/*--- update_rates() ------------------------------------------
  measures the rates,.. defines the TaCoS-systemtime,..
  calculates some rates and rate relevant stuff
  _and_ checks for a overall emergency,..
  the actual time for the _whole_ TaCoS(target) system is the
  time measured right before reading out the scalers,..
  ---------------------------------------------------------- */
void update_rates(data_t *d) {
  int     i, j, tt, tc;
  rate_t  helper;
  
  /* set the new target-system time */
  getclock(TIMEOFDAY, &Scaler.new.time);
  d->time = Scaler.new.time;

  if(!Testrun) {
    /* !!!!!!!!!!!!! REAL RUNNING !!!!!!!!!!!!!!!! */
    for(i=0; i<NETRATE_NUM_RATES; i++) {
      /* read all scaler values */ 
      Scaler.new.value[i]  = read_scaler(i);
      /* calculate the difference,.. */
      if(Scaler.new.value[i]<Scaler.old.value[i]) 
	Scaler.diff.value[i] = (SCALER_OVERFL-Scaler.old.value[i])+Scaler.new.value[i];
      else
	Scaler.diff.value[i] = Scaler.new.value[i]-Scaler.old.value[i];  
    }
  } else {
    /* !!!!!!!!!!!!! TESTRUN !!!!!!!!!!!!!!!!!!!!! */
    for(i=0; i<NETRATE_NUM_RATES; i++) {
      Scaler.diff.value[i] = (scaler_t) 20000 
	  + 5000 * sin((double) Scaler.new.time.tv_sec/400 * (double)i );  
    }
    Scaler.diff.value[RATE_REF_RATE] = (scaler_t) 1000000;  
    Scaler.diff.value[18] = (scaler_t) 100; /* just for this very special setting*/  
    Scaler.diff.value[27] = (scaler_t) 13; /* just for this very special setting*/  
    Scaler.diff.value[25] = (scaler_t) 14; /* just for this very special setting*/  
    Scaler.diff.value[20] = (scaler_t) 12; /* just for this very special setting*/  
    Scaler.diff.value[28] = (scaler_t) 15; /* just for this very special setting*/  
  }
  
  /* calculate the time difference in microseconds*/
  if (!Use_ref_rate) {
    /* Calculate the rates using the time interval between two scaler readouts */
    /* but be careful, the timmer resolution is only 0.01 s !!! */
    Scaler.diff.time = (long) 
      (   (Scaler.new.time.tv_sec  - Scaler.old.time.tv_sec) * 1000000 
	+ (Scaler.new.time.tv_nsec - Scaler.old.time.tv_nsec)/ 1000 );
  } else {
    /* via reference rate */
    Scaler.diff.time = (long) 
      ((double)Scaler.diff.value[RATE_REF_RATE]/(double)RATE_REF_FREQ ) ;
  }
  

  /* Now calculate all rates */          
  for(i=0; i<NETRATE_NUM_RATES; i++) {
    Rate.misc[i] = (rate_t) (  (double)Scaler.diff.value[i]/
			       ( (double) Acc[i] * (double)Scaler.diff.time )  );
  }

  /*-------- some rate calculation follow here --------------*/
  /* calculate the interaction rate and the bg*/
  if(!Testrun) {
    /* !!!!!!!!!!!!! REAL RUNNING !!!!!!!!!!!!!!!! */
    Rate.ia = (Rate.misc[2] + Rate.misc[3] +  
	       Rate.misc[4] + Rate.misc[5])/4; 
    Rate.bg = Rate.misc[RATE_BG];
  } else {
    /* !!!!!!!!!!!!! SIMULATION !!!!!!!!!!!!!!!!!! */
    Rate.ia = getSimIA();
    Rate.bg = getSimBG();
  }

  /*---------------------------------------------------------*/
  /* intergrated ia/ means/sigmas */
  tc = d->autom.timeco;
  /* calculate the integrated interaction rate,.. -> lumi (Hz*sec)*/
  Rate.integrated_ia += (intrate_t) (Rate.ia * Scaler.diff.time);
  /* Calculate the mean of the measured rate */
  Rate.ia_mean = Rate.ia * (((double)   1)/((double)tc)) 
    + Rate.ia_mean       * (((double)tc-1)/((double)tc));
  /* Calculate the sigma of the measured rate,..  and its mean  */
  Rate.sigma =  fabs(Rate.ia - Rate.ia_mean);
  Rate.sigma_mean = Rate.sigma  * (((double)   1)/((double)tc))
    + Rate.sigma_mean           * (((double)tc-1)/((double)tc));
  /*---------------------------------------------------------*/
  /* single wire rates calculation */
  for(i=0; i<NP_NUM_WIRES; i++) {
    helper = (rate_t) MAX((double) 0.0, 
			  (double) ((double)Rate.misc[WR.rpos[i]]*(double)WR.calf[i]+(double)WR.pedestal[i]) 
			  );
    WR.rate[i][WR.pos]   = helper;
    Rate.wrate_nomean[i] = helper;

    Rate.wrate[i] = 0;
    for(j=0; j<R_NUM_WRRATES; j++) {
      Rate.wrate[i] += WR.rate[i][j];
    }
    Rate.wrate[i] /= R_NUM_WRRATES;
  }
  if(WR.pos == R_NUM_WRRATES -1) {
    WR.pos = 0;
  } else {
    WR.pos ++;
  }
  /*---------------------------------------------------------*/
  /* lin. reg for ia-rate -> filling of the array */
  DR.pos++;
  if(DR.pos>d->autom.dr_slices-1) {
    DR.pos = 0;
  }
  tt = 0;
  DR.rate[DR.pos] = Rate.ia;
  for(i=DR.pos; i<d->autom.dr_slices; i++) {
    DR.srateq[i] = (rate_t) 
      ((double)DR.rate[i]*d->autom.dr_factor*
       (double)((double)((double)d->autom.dr_offset+(double)tt)/(double)d->autom.dr_slices));
    DR.srateq[i] *= DR.srateq[i];
    tt++;
  }
  for(i=0; i<DR.pos; i++) {
    DR.srateq[i] = (rate_t) 
      ((double)DR.rate[i]*d->autom.dr_factor*
       (double)((double)((double)d->autom.dr_offset+(double)tt)/(double)d->autom.dr_slices));
    DR.srateq[i] *= DR.srateq[i];
    tt++;
  }

  /*print_scrollX("%f - %g", DR.rate[DR.pos], DR.srateq[DR.pos]);*/
  calc_dr(d);
  Rate.dr = DR.dr;

  /* return back the rates*/
  d->rate = Rate;

  /* fill in the 'old'-structure */
  Scaler.old.time = Scaler.new.time;
  for(i=0; i<NETRATE_NUM_RATES; i++) {
    Scaler.old.value[i] = Scaler.new.value[i];
  }
}