/********************************************************************\

  Name:         mwpc.c
  Created by:   Nikolai Khomutov
  Contents:     PiBeta Analyzer module for MWPC calculations

  $Log: mwpc.c,v $
  Revision 1.32  2000/06/06 13:31:29  midas
  Fixed compiler warning

  Revision 1.31  2000/06/06 10:23:41  midas
  Removed efficiency calculation (now lives in histo.c)

  Revision 1.30  2000/05/05 13:12:51  midas
  Added run mode check to chamber efficiency

  Revision 1.29  2000/03/31 15:00:47  midas
  Fixed bugs with debug messages, added 12th trigger everywhere

  Revision 1.28  2000/03/30 14:03:47  midas
  Adde global debug message flag

  Revision 1.27  2000/03/17 10:41:31  midas
  Fixed typo in comment

  Revision 1.26  1999/12/13 11:35:24  midas
  Switch automatically between old and new chamber 1 depending on run number

  Revision 1.25  1999/10/12 13:04:23  midas
  Added Theta/Phi 2D histos

  Revision 1.24  1999/08/24 08:00:20  midas
  Calculate efficiencies only for pienu:Hi trigger

  Revision 1.23  1999/08/20 14:23:12  midas
  Added errors to gauss_fit()

  Revision 1.22  1999/08/19 15:27:34  midas
  Clear histos independent of global histo clear flag

  Revision 1.21  1999/08/19 15:17:58  midas
  Added "guess initial param" flag to gauss_fit

  Revision 1.20  1999/08/18 08:03:39  midas
  Don't clear raw ADC histos after each run

  Revision 1.19  1999/08/17 07:14:13  midas
  Write average pedestal width to ODB

  Revision 1.18  1999/08/16 10:48:40  midas
  Introduced 2D histos by CLS

  Revision 1.17  1999/08/12 10:49:22  midas
  Changed histo binning

  Revision 1.16  1999/07/14 08:09:24  midas
  Re-arranged histo ID's

  Revision 1.15  1999/07/06 09:40:56  midas
  Fixed crash caused by wire_fired being only 192 entries long

  Revision 1.14  1999/06/30 11:58:28  midas
  Added statistics check on pedestal evaluation

  Revision 1.13  1999/05/29 10:58:05  midas
  Fixed compiler warning

  Revision 1.12  1999/05/29 10:56:44  midas
  Fixed compiler warnings

  Revision 1.11  1999/02/16 10:33:03  midas
  Defined INVALID as -1000

  Revision 1.10  1999/02/05 08:43:58  midas
  Fixed minor bug reported by Kolya

  Revision 1.9  1999/02/03 13:52:08  midas
  Fixed minor bug in MWPC bank filling

  Revision 1.8  1999/02/01 14:39:55  midas
  Modified MWPC bank into arrays

  Revision 1.7  1999/02/01 14:03:08  midas
  Added 4th hit im MWPC bank

  Revision 1.6  1999/02/01 11:28:05  midas
  - Clear mwpchits array for each event
  - Defined histo 4803

  Revision 1.5  1999/01/22 10:21:55  midas
  Added bug fixes received from Kolya on 22-jan-99

  Revision 1.3  1998/12/16 12:54:05  midas
  Initial code for cahtode analysis

  Revision 1.2  1998/12/03 11:22:49  midas
  Added Log in header

\********************************************************************/

/*-- Include files -------------------------------------------------*/

/* standard includes */
#include <stdio.h>
#include <time.h>
#include <math.h>

/* midas includes */
#include "midas.h"
#include "hardware.h"
#include "experim.h"
#include "analyzer.h"      

#include "util.h"

/* cernlib includes */
#ifdef OS_WINNT
#define VISUAL_CPLUSPLUS
#endif

#ifdef __linux__
#define f2cFortran
#endif

#include <cfortran.h>
#include <hbook.h>

#ifndef PI
#define PI 3.14159265359
#endif

/*-- Parameters ----------------------------------------------------*/

MWPC_DECODE_PARAM       mwpc_decode_param;
extern TRIGGER_SETTINGS trigger_settings;
extern GLOBAL_PARAM     global_param;
extern EXP_PARAM        exp_param;
extern RUNINFO          runinfo;
extern RTES_BANK        rates;
extern char             *analyzer_name;

/*-- Static parameters ---------------------------------------------*/

#define CLUSTER_GAP 1 /* allow one wire missing in a cluster */

#define N_GROUP_MEMBER 16 /* # of channels in each group for 2nd ped. corr. */

#define WRITE_EVENT_SAMPLE 0

#define NZ_MAX 1000

const int i_sfc[] = { 0, 64, 128, 192, 256, 448};
const int n_sfc[] = {64, 64,  64,  64, 192, 192};

const int    wirenum[]  = { 192,    384   };
const double stripnum[] = {  64.0,  192.0 };
const double anodrad[]  = {  60.15, 120.1 };
const double mwpclen[]  = { 352.4,  540   };

const double cathrad1[] = {  57.6,   62.5 };
const double cathrad2[] = { 117.6,  122.6 };

double phi_corr_in_1;
double phi_corr_in_2;
// WL
double z_offset_1;
double z_offset_2;
double phi_corr_wire[2];


double phi_max;

double ctg1i, ctg1o, ctg2i, ctg2o;

/*-- Module declaration --------------------------------------------*/

INT mwpc_decode(EVENT_HEADER*,void*);
INT mwpc_init(void);
INT mwpc_bor(INT run_number);  
INT mwpc_eor(INT run_number);  

MWPC_DECODE_PARAM_STR(mwpc_decode_param_str);
ANA_MODULE mwpc_decode_module = {
  "MWPC Decode",                 /* module name           */  
  "Stefan Ritt",                 /* author                */
  mwpc_decode,                   /* event routine         */
  mwpc_bor,                      /* Begin of run routine  */
  mwpc_eor,                      /* End of run routine    */
  mwpc_init,                     /* init routine          */
  NULL,                          /* exit routine          */
  &mwpc_decode_param,            /* parameter structure   */
  sizeof(mwpc_decode_param),     /* structure size        */
  mwpc_decode_param_str,         /* initial parameters    */
};

/*-- globals -------------------------------------------------------*/

static int chamber1;       /* chamber1 can be OLD or NEW */
#define OLD 1
#define NEW 2

#define ADC_N_BINS   500
#define ADC_X_LOW      0
#define ADC_X_HIGH  5000

#define PED_N_BINS   200   /* normal pedestals */
#define PED_X_LOW      0
#define PED_X_HIGH  1000

#define PED2_N_BINS  120   /* corrected pedestals */
#define PED2_X_LOW (-300)
#define PED2_X_HIGH  300

/*-- initialize constants ------------------------------------------*/

void init_constants(INT run_number)
{
  /* chamber 1 switchover between runs 21795 and 21796 */
  if (run_number < 21796)
    {
      chamber1 = OLD;
      printf("run%d  with old chamber.\n", run_number);
    } 
  else
    chamber1 = NEW;
  
  if (chamber1 == OLD)
    {
    phi_corr_in_1 = -0.05;
    
    ctg1i = (1.f/tan(PI/180 * 36.96))*(anodrad[0]/cathrad1[0]);
    ctg1o = (1.f/tan(PI/180 * 33.65))*(anodrad[0]/cathrad1[1]);
    }
  else
    {
    phi_corr_in_1 = -0.0;
    
    ctg1i = (1.f/tan(PI/180 * 45.2958))*(anodrad[0]/cathrad1[0]);
    ctg1o = (1.f/tan(PI/180 * 40.4306))*(anodrad[0]/cathrad1[1]);
    }

  /* common constants */
  phi_max       =   4.0;

  /* mwpc 2 */
  phi_corr_in_2 =   -0.362;

  /* parameters added by WL for allignment purposes   */
  phi_corr_wire[0] = 0.;         //offsets bw wire and cathods
  phi_corr_wire[1] = 0.6130; 
  
  z_offset_1 = 0.931;
  z_offset_2 = -0.931;
  /* end of parameters added by WL*/

  ctg2i = (1.f/tan(PI/180 * 44.55))*(anodrad[1]/cathrad2[0]);
  ctg2o = (1.f/tan(PI/180 * 42.30))*(anodrad[1]/cathrad2[1]);
}

/*-- init routine --------------------------------------------------*/

INT mwpc_init(void)
{
  init_constants(runinfo.run_number);

  /* cathode histo arrays */
  HBOOK2(MWPC_ADC_ID_BASE, "Cath", ADC_N_BINS, (float)ADC_X_LOW, (float)ADC_X_HIGH, N_CATH, (float)0, (float)N_CATH, (float)0);
  HBOOK2(MWPC_PED_ID_BASE, "Cath Ped", PED_N_BINS, (float)PED_X_LOW, (float)PED_X_HIGH, N_CATH, (float)0, (float)N_CATH, (float)0);
  HBOOK2(MWPC_PED2_ID_BASE, "Cath Corr Ped", PED2_N_BINS, (float)PED2_X_LOW, (float)PED2_X_HIGH, N_CATH, (float)0, (float)N_CATH, (float)0);

  /* Wire profiles */
  HBOOK1(MWPC_ID_BASE+0, "N_Wires_1", 576, 0.0f, 576.0f, 0.0f); 
  HBOOK1(MWPC_ID_BASE+1, "MWPC_profile_1", 192, 0.0f, 192.0f, 0.0f); 
  HBOOK1(MWPC_ID_BASE+2, "N_Wires_2", 576, 0.0f, 576.0f, 0.0f); 
  HBOOK1(MWPC_ID_BASE+3, "MWPC_profile_2", 384, 0.0f, 384.0f, 0.0f); 

  /* Theta/Phi map */
  HBOOK2(MWPC_ID_BASE+4, "MWPC1_theta_phi", 360, 0.0f, 360.0f, 180, 0.0f, 180.0f, 0.0f);
  HBOOK2(MWPC_ID_BASE+5, "MWPC2_theta_phi", 360, 0.0f, 360.0f, 180, 0.0f, 180.0f, 0.0f);

  /* check alignment of wires and cathode strips WL */
  HBOOK2(MWPC_ID_BASE+10, "ch1 strip phi vs diff phi", 18, 0.0, 360., 32, -4., 4.,0.f);
  HBOOK2(MWPC_ID_BASE+11, "ch2 strip phi vs diff phi", 18, 0.0, 360., 32, -4., 4.,0.f);
  
  return SUCCESS;
}

/*-- bor routine ---------------------------------------------------*/

INT mwpc_bor(int run_number)
{
  init_constants(run_number);

  /* clear histos independent of global clear histo flag */
  HRESET(MWPC_ADC_ID_BASE, " ");
  HRESET(MWPC_PED_ID_BASE, " ");
  HRESET(MWPC_PED2_ID_BASE, " ");

  return SUCCESS;
}

/*-- eor routine ---------------------------------------------------*/

static float yped[N_CATH][PED_N_BINS];
static float yped2[N_CATH][PED2_N_BINS];

INT mwpc_eor(INT run_number)
{
int    i, i_channel, pedestals_present;
double x[PED_N_BINS], y[PED_N_BINS];
double param[3], error[3], ymax, xmax, chisqr;
double sy, sxy, sxxy, mean, sigma;
char   str[256];
HNDLE  hDB, hkey;
FILE   *f;
time_t now;

  cm_get_experiment_database(&hDB, NULL);

  if (mwpc_decode_param.correct_pedestals)
   /* move ADC values into array */    
   HUNPAK(MWPC_PED_ID_BASE, yped[0][0], "HIST", 1);
    {
    /*---- analyze pedestal histo ------------------------------------*/

    for (i=0 ; i<PED_N_BINS ; i++)
      x[i] = (double) i/PED_N_BINS*(PED_X_HIGH-PED_X_LOW)+PED_X_LOW;

    /* loop through ADC channels */
    pedestals_present = TRUE;
    for (i_channel=0; i_channel<N_CATH; i_channel++)
      {
      ymax = 0;
      for (i=0 ; i<PED_N_BINS ; i++)
        {
        y[i] = yped[i_channel][i];
        if (y[i] > ymax && i>0)
          ymax = y[i];
        }

      if (ymax < 20)
        {
        pedestals_present = FALSE;
        break;
        }

      /* calculate mean and sigma */
      for (i=0,sy=sxy=sxxy=0.f ; i<PED_N_BINS ; i++)
        {
        sxy += x[i]*y[i];
        sxxy += x[i]*x[i]*y[i];
        sy += y[i];
        }

      if (sy > 0)
        {
        mean = sxy/sy;
        sigma = sqrt((sxxy - (sxy*sxy)/sy)/sy);

        mwpc_decode_param.pedestal[i_channel] = (int) mean;
        mwpc_decode_param.pedestal_sigma[i_channel] = (float) sigma;

        pedestals_present++;
        }
      }

    /* write new pedestal values to online database. */
    if (pedestals_present)
      {
      sprintf(str,"/%s/Parameters/MWPC decode", analyzer_name);
      db_find_key(hDB, 0, str, &hkey);
      db_set_record(hDB, hkey, &mwpc_decode_param, sizeof(mwpc_decode_param), 0);
      }

    /*---- analyze corrected pedestal histos -------------------------*/

    if (pedestals_present)
      {
      /* move ADC values into array */
      HUNPAK(MWPC_PED2_ID_BASE, yped2[0][0], "HIST", 1);

      for (i=0 ; i<PED2_N_BINS ; i++)
        x[i] = (float) i/PED2_N_BINS*(PED2_X_HIGH-PED2_X_LOW)+PED2_X_LOW;

      /* loop through ADC channels */
      for (i_channel=0; i_channel<N_CATH; i_channel++)
        {

        /* do gauss fit */
        for (i=0,ymax=0 ; i<PED2_N_BINS ; i++)
          {
          y[i] = yped2[i_channel][i];
          if (y[i] > ymax)
            ymax = y[i];
          }

        if (ymax > 0)
          {
          /* printf("Fitting channel %d\r", i_channel); */
          chisqr = gauss_fit(x, y, PED2_N_BINS, param, error, TRUE);
          if (chisqr != -1)
            {
            mwpc_decode_param.pedestal_2nd[i_channel] = (int) (param[1]+0.5);
            mwpc_decode_param.pedestal_sigma_2nd[i_channel] = (float) fabs(param[2]/2.35482);

            /* mark histo */
            i = (int) (((param[1]+3*mwpc_decode_param.pedestal_sigma_2nd[i_channel])-PED2_X_LOW) / 
                      (PED2_X_HIGH-PED2_X_LOW)*PED2_N_BINS+0.5);
            if (i>0 && i<PED2_N_BINS)
              yped2[i_channel][i] = (float) ymax;              
            }
          }
        }

      HPAK(MWPC_PED2_ID_BASE, yped2[0]);

      /* calculate average sigma */
      for (i=0,xmax=0; i<N_CATH ; i++)
        xmax += mwpc_decode_param.pedestal_sigma_2nd[i];
      mwpc_decode_param.pedestal_average = (float) (xmax / N_CATH);

      /* write new pedestal values to online database. */
      sprintf(str,"/%s/Parameters/MWPC decode", analyzer_name);
      db_find_key(hDB, 0, str, &hkey);
      db_set_record(hDB, hkey, &mwpc_decode_param, sizeof(mwpc_decode_param), 0);

      if (global_param.debug_messages)
        cm_msg(MINFO, "mwpc_eor", "Average pedestal sigma MWPC ADC: %1.2lf", 
               mwpc_decode_param.pedestal_average);

      /* write pedestals sigma to file */
      add_data_dir(str, "cath_ped2_sigma.txt");
      f = fopen(str, "a");
      if (f == NULL)
        cm_msg(MERROR, "mwpc_eor", "Cannot open file %s", str);
      else
        {
        time(&now);
        strcpy(str, ctime(&now));
        str[24] = 0;
        fprintf(f, "%s\trun %d\t", str, run_number);
        for (i=0 ; i<N_CATH ; i++)
          fprintf(f, "%1.2f\t", mwpc_decode_param.pedestal_sigma_2nd[i]);
        fprintf(f, "\n");
        fclose(f);
        }
      }
    /*
    else
      cm_msg(MINFO, "mwpc_eor", "Not enough statistics to evaluate cathode pedestals");
    */

    }

  return SUCCESS;
}

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

/* data structures used for MWPC calculations */

typedef struct {
  double z;
  double phi;
  double theta;
  int phi_ok;
} MWPCPOINT;

typedef struct {
  int first;
  int last;
  int i_max;
  int v_max;
  float center;
} PEAK;

/*---- helper functions for cathode analysis -----------------------*/

void find_center(float strip[], int surface, PEAK *peak)
{
int    i, j, width, v;
long   s1, s2;
double c, delta;
static const double fcorr = 0.5;

  width = peak->last - peak->first + 1;
  if (width < 0) 
    width += n_sfc[surface];

  peak->v_max = -1;
  s1 = 0;
  s2 = 0;

  /* calculate center of gravity for peak */
  for (i = peak->first; i<(peak->first+width); i++)
    {
    if (i < n_sfc[surface])
      j = i;
    else
      j = i-n_sfc[surface];

    v = (int) strip[i_sfc[surface]+j];

    if (v > peak->v_max)
      {
      peak->v_max = v;
      peak->i_max = i;
      }

    s1 += v;
    s2 += i*v;
    }

  if (s1)
    c = (double) s2 / (double) s1;
  else
    c = 0;

  /* Correction of systematic deviations ( Dissertation by C.Grab, I-23 ) */
  delta = c - peak->i_max;
  if (fabs(delta) < 1.0)
    {
    delta += fcorr*delta*(0.5 - fabs(delta));
    c = peak->i_max + delta;
    }

  if (c > n_sfc[surface]) 
    c -= n_sfc[surface];
  peak->center = (float) c;
}

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

void find_z(int surface, double cogi, double cogo, int *nz, MWPCPOINT **zp)
{
double strip_step;
double cgi, cgo, r;
double z, z1, z2, phi, theta, ct1, ct2, delta_phi;
int    k, k1, k2;
 double z_offset;  //WL

  switch(surface)
    {
    case 0 :
      r = anodrad[0];
      strip_step = 2*PI/stripnum[0];
      ct1 = ctg1i;
      ct2 = ctg1o;

      if (chamber1 == OLD)
        cgi = cogi;
      else
        cgi = cogi+0.5;

      cgi += phi_corr_in_1;
      if (cgi >= stripnum[0])
        cgi -= stripnum[0];
      if (cgi < 0)
        cgi += stripnum[0];

      if (chamber1 == OLD)
        cgo = stripnum[0] - 1 - cogo;
      else
        cgo = stripnum[0] - 1 - cogo + 0.5;

      if (cgo < 0)
        cgo += stripnum[0];

      z1 = -0.5 * mwpclen[0];
      z2 = 0.f;

      k1 = -2;
      k2 = 0;
      break;

    case 2 :
      r = anodrad[0];
      strip_step = 2*PI/stripnum[0];
      ct1 = ctg1i;
      ct2 = ctg1o;

      if (chamber1 == OLD)
        cgi = stripnum[0] - 1 - cogi;
      else
        cgi = stripnum[0] - 1 - cogi + 0.5;

      cgi += phi_corr_in_1;
      if (cgi < 0)
        cgi += stripnum[0];

      if (chamber1 == OLD)
        cgo = cogo;
      else
        cgo = cogo+0.5;      

      if (cgo < 0)
        cgo += stripnum[0];

      z1 = 0.f;
      z2 = 0.5 * mwpclen[0];

      k1 = 0;
      k2 = 2;
      break;

    case 4 :
      r = anodrad[1];
      strip_step = 2*PI/stripnum[1];
      ct1 = ctg2i;
      ct2 = ctg2o;

      cgi = stripnum[1] - 1 - cogi;
      cgi += phi_corr_in_2;

      if (cgi < 0)
        cgi += stripnum[1];

      cgo = cogo;

      z1 = -0.5 * mwpclen[1];
      z2 =  0.5 * mwpclen[1];

      k1 = -1;
      k2 = 1;
      break;
    }


  cgi *= strip_step;
  cgo *= strip_step;

  // by WL
  if (surface<4)
    z_offset = z_offset_1;
  else
    z_offset = z_offset_2;
  for (k=k1; k<=k2; k++)
    {
    delta_phi = cgo - cgi + 2*PI*k;
    //z = r * delta_phi/(ct1 + ct2);  by WL
    z = r * delta_phi/(ct1 + ct2)+z_offset;
    if ((z>=z1) && (z<=z2))
      {
      (*zp)->z = z;

      phi = 180 - 180/PI*((delta_phi * ct1/(ct1 + ct2)) + cgi);
      while (phi>=360.f)
        phi -= 360.f;
      while (phi<0)
        phi += 360.f;
      (*zp)->phi = phi;

      if (z==0.f)
        theta = 90.f;
      else
        theta = 180/PI*atan2(r, z);
      (*zp)->theta = theta;
      if (*nz < (NZ_MAX-1))
        {
        (*nz)++;
        (*zp)++;
        }
      }
    }
}

/*-- event routine -------------------------------------------------*/

INT mwpc_decode(EVENT_HEADER *pheader, void *pevent)
{
int          i, j, chamber, index, j_min;
int          n_pcos, n_cath, n_zcath;
MWPC_BANK    *mwpc;
WORD         *pcos, *trig;
LRS1882_DATA *lrs1882;
float        *ccat, *strip;
DWORD        *zcat;
float        min, median;
static FILE  *f = NULL;
char         str[256];
float        member[N_GROUP_MEMBER];

int          wire_fired[2][384];
int          n_wires_fired[2];
int          n_wire_clusters[2];
float        cluster_center[2][200];
WORD         cluster_width[2][200];

int          n_peaks[6];
PEAK         peaks[6][192];

int          n_z[2];
MWPCPOINT    zhit[2][NZ_MAX];
MWPCPOINT    mwpchits[2][384];

int          surface;
int          in_peak;
MWPCPOINT    *zp_1, *zp_2;
double       vmin, phi_delta;
int          jmin;

float        tmp;
 float        tmp_diff;
  /* look for PCOS, CATH or ZCAT bank, return if not present */
  n_pcos = bk_locate(pevent, "PCOS", &pcos);
  
  n_cath = bk_locate(pevent, "CATH", &lrs1882);
  n_zcath = bk_locate(pevent, "ZCAT", &zcat);

  if (pcos == NULL || (n_cath < N_CATH && zcat == NULL))
    {
    if (n_cath > 0 && n_cath < N_CATH)
      cm_msg(MERROR, "mwpc_decode", 
        "CATH bank only contains %d entries instead of %d", n_cath, N_CATH);
    return 1;
    }
  
  if (!bk_locate(pevent, "TRIG", &trig))
    return 1;

  /* create calibrated cathode bank */
  bk_create(pevent, "CCAT", TID_FLOAT, &ccat);

  /*---- decompress ZCAT bank --------------------------------------*/

  if (zcat)
    {
    memset(ccat, 0, sizeof(float)*N_CATH);
    for (i=0 ; i<n_zcath ; i++)
      ccat[zcat[i] >> 16] = (float) ((short int) (zcat[i] & 0xFFFF));

    bk_close(pevent, ccat+N_CATH);
    }

  /*---- decode CATH bank ------------------------------------------*/

  if (lrs1882)
    {
    memset(ccat, 0, sizeof(float)*N_CATH);
    for (i=0 ; i<N_CATH ; i++)
      {
      if (n_cath == N_CATH)
        ccat[i] = (float ) lrs1882[i].data;
      else
        {
        /* correct for 641 channels (1998 data had 65 channels on slot 3) */
        if (i>=256)
          ccat[i] = (float ) lrs1882[i+1].data;
        else
          ccat[i] = (float ) lrs1882[i].data;
        }
      }
    bk_close(pevent, ccat+N_CATH);
    }

  /*---- correct swapped cables & dead channels --------------------*/

  if (chamber1 == NEW)
    {
    for (i=0 ; i<16 ; i++)
      {
      tmp = ccat[i+16];
      ccat[i+16] = ccat[i+32];
      ccat[i+32] = tmp;
      }
    }

  ccat[124] = (float) (0.5*(ccat[123]+ccat[125]));
  ccat[468] = (float) (0.5*(ccat[467]+ccat[469]));

  /*---- pedestal run code -----------------------------------------*/

  if (exp_param.run_mode == RUN_MODE_PEDESTAL)
    {
    /* fill histos */
    for (i=0 ; i<N_CATH ; i++)
      HF2(MWPC_PED_ID_BASE, ccat[i], (float)i, 1.f);

    return 1;
    }

  /*---- pedestal correction only on raw data ----------------------*/

  if (lrs1882)
    {

    /* fill pedestal histos */
    for (i=0 ; i<N_CATH ; i++)
      HF2(MWPC_PED_ID_BASE, ccat[i], (float)i, 1.f);

    /*---- primary pedestal subtraction ----*/

    for (i=0 ; i<N_CATH ; i++)
      ccat[i] = ccat[i] - mwpc_decode_param.pedestal[i];

    /*---- secondary pedestal subtraction ----*/

    if (f == NULL && pheader->serial_number < 100 && WRITE_EVENT_SAMPLE)
      {
      add_data_dir(str, "event_sample.txt");
      f = fopen(str, "w");
      }
    if (f)
      fprintf(f, "Event #%d ---------------------------------------\n", pheader->serial_number);

    if (pheader->serial_number > 100 && f)
      {
      fclose(f);
      f = NULL;
      }

    for (index = 0 ; index < N_CATH/N_GROUP_MEMBER ; index++)
      {                                           
      /* find median */
      for (i=index*N_GROUP_MEMBER ; i<(index+1)*N_GROUP_MEMBER ; i++)
        {
        if (mwpc_decode_param.pedestal_sigma[i] == 0) /* sigmas could initially be zero */
          member[i-index*N_GROUP_MEMBER] = ccat[i] / 100;
        else
          member[i-index*N_GROUP_MEMBER] = ccat[i] / mwpc_decode_param.pedestal_sigma[i];
        }

      for (i=0 ; i<N_GROUP_MEMBER/2 ; i++)
        {
        min = 10000.f;
        for (j=0 ; j<N_GROUP_MEMBER ; j++)
          {
          if (member[j] < min)
            {
            min = member[j];
            j_min = j;
            }
          }
        median = min;
        member[j_min] = 100000.f;
        }

      if (f && index == 30)
        {
        fprintf(f, "madian %1.0f (%d)\nratio:  ", 
          median*mwpc_decode_param.pedestal_sigma[j_min+index*N_GROUP_MEMBER], j_min);

        for (i=index*N_GROUP_MEMBER ; i<(index+1)*N_GROUP_MEMBER ; i++)
          fprintf(f, "%6.2f ", mwpc_decode_param.pedestal_sigma[i]/ 
                               mwpc_decode_param.pedestal_sigma[j_min+index*N_GROUP_MEMBER]);

        fprintf(f, "\nbefore: ");
        for (i=index*N_GROUP_MEMBER ; i<(index+1)*N_GROUP_MEMBER ; i++)
          fprintf(f, "%6.0f ", ccat[i]);
        fprintf(f, "\nafter:  ");
        }

      for (i=index*N_GROUP_MEMBER ; i<(index+1)*N_GROUP_MEMBER ; i++)
        ccat[i] -= median * mwpc_decode_param.pedestal_sigma[i];
      
      if (f && index == 30)
        {
        for (i=index*N_GROUP_MEMBER ; i<(index+1)*N_GROUP_MEMBER ; i++)
          fprintf(f, "%6.0f ", ccat[i]);
        fprintf(f, "\n\n");
        }
      }

    /* fill corrected pedestal histos */
    for (i=0 ; i<N_CATH ; i++)
      HF2(MWPC_PED2_ID_BASE, ccat[i], (float)i, 1.f);
    }

  /*---- create compressed ZADC bank -------------------------------*/

  if (zcat == NULL && runinfo.online_mode == 0)
    {
    bk_create(pevent, "ZCAT", TID_DWORD, &zcat);

    /* fill ADC over threshold */
    for (i=0 ; i<N_CATH ; i++)
      if (ccat[i] > mwpc_decode_param.pedestal_2nd[i] +
                    mwpc_decode_param.pedestal_sigma_2nd[i] * 3)
        *zcat++ = (i<<16) | ((WORD) (ccat[i]+0.5));

    bk_close(pevent, zcat);
    }

  /*---- software gain calibration ---------------------------------*/

  for (i=0 ; i<N_CATH ; i++)
    ccat[i] *= mwpc_decode_param.software_gain[i];

  /*---- fill cathode histos ---------------------------------------*/

  /* fill ADC histos if above threshold and current trigger 
     pattern has bits in common with the histo mask */
  for (i=0 ; i<N_CATH ; i++)
    if (ccat[i] > mwpc_decode_param.pedestal_2nd[i] +
                  mwpc_decode_param.pedestal_sigma_2nd[i] * 3)
      if (trig[0] & global_param.histo_mask)
        HF2(MWPC_ADC_ID_BASE, ccat[i], (float)i, 1.f);

  /*---- analyze PCOS data -----------------------------------------*/

  n_wires_fired[0] = n_wires_fired[1] = 0;
  for (i=0 ; i < n_pcos ; i++)
    {
    chamber  = pcos[i]>>10 & 1;
    wire_fired[chamber][n_wires_fired[chamber]++] = pcos[i]>>1 & 0x1FF;
    }

  /* fill wire histos */
  HF1(MWPC_ID_BASE, (float) n_wires_fired[0], 1.0f);            
  for (i=0 ; i < n_wires_fired[0] ; i++) 
    HF1(MWPC_ID_BASE+1, (float) wire_fired[0][i], 1.0f); 

  HF1(MWPC_ID_BASE+2, (float) n_wires_fired[1], 1.0f);            
  for (i=0 ; i < n_wires_fired[0] ; i++) 
    HF1(MWPC_ID_BASE+3, (float) wire_fired[1][i], 1.0f); 

  /*---- combine adjacent wires into single cluster ----------------*/

  n_wire_clusters[0] = n_wire_clusters[1] = -1;
  for (chamber = 0 ; chamber < 2 ; chamber++)
    {
    if (n_wires_fired[chamber] > 0)
      {
      for (i=0 ; i<(INT)n_wires_fired[chamber] ; i++)
        {
        if (i && ((wire_fired[chamber][i]-wire_fired[chamber][i-1]) == 1))
          {  
          /* expand current cluster */
          cluster_center[chamber][n_wire_clusters[chamber]] += 0.5f;
          cluster_width[chamber][n_wire_clusters[chamber]]++;
          }
        else
          {  
          /* start new cluster */
          cluster_center[chamber][++n_wire_clusters[chamber]] = 
            (float) wire_fired[chamber][i];
          cluster_width[chamber][n_wire_clusters[chamber]] = 1;
          }
        } /* end for (i = 0; i < n_wires_fired[chamber]; i++) */

      /* concatenate first & last clusters when necessary*/
      if ( (n_wire_clusters[chamber] > 0) && 
           (wire_fired[chamber][0] == 0) &&
           (wire_fired[chamber][n_wires_fired[chamber]-1] == wirenum[chamber] - 1) )
        {
        cluster_center[chamber][0] -=
          0.5f * cluster_width[chamber][n_wire_clusters[chamber]];
        if (cluster_center[chamber][0] < 0)
          cluster_center[chamber][0] += wirenum[chamber];
        cluster_width[chamber][0] += cluster_width[chamber][n_wire_clusters[chamber]];
        n_wire_clusters[chamber]--;
        }

      /* convert wire clusters to degrees */
      for (i = 0; i <= n_wire_clusters[chamber]; i++)
        {
        cluster_center[chamber][i] =
          180.f - cluster_center[chamber][i] * 360.f / wirenum[chamber]+phi_corr_wire[chamber];
	
        if (cluster_center[chamber][i] < 0)
          cluster_center[chamber][i] += 360.f;
        }
      } /* end if (n_wires_fired[chamber]>0) */

    n_wire_clusters[chamber]++;
    } /* end for (chamber = 0 ; chamber < 2 ; chamber ++) */

  /*---- analyzer cathode data, find peaks for all surfaces --------*/

  strip = ccat;
 
  /* loop over all surfaces */
  for (surface=0; surface<6; surface++)
    {
    in_peak = 0;
    j = -1;
    n_peaks[surface] = 0;
    for (i=0; i<n_sfc[surface]; i++)
      {
      /* only consider ADC values above the noise level (3 sigma) */
      if (strip[i+i_sfc[surface]] > 
            mwpc_decode_param.pedestal_2nd[i+i_sfc[surface]] +
            mwpc_decode_param.pedestal_sigma_2nd[i+i_sfc[surface]] * 3)
        {
        if (in_peak)
          peaks[surface][j].last++;  /* continue current peak */
        else
          {  
          /* start new peak */
          in_peak = 1;
          j++;
          peaks[surface][j].first = i;
          peaks[surface][j].last  = i;
          }
        }
      else if (in_peak)
        in_peak = 0;        /* end current peak */
      } /* end for (i=0; i<n_sfc[surface]; i++) */

    if (j>=1)
      if ((peaks[surface][j].last == (n_sfc[surface]-1))
         && (peaks[surface][0].first == 0))
        {  
        /* concatenate first & last peaks when necessary*/
        peaks[surface][0].first = peaks[surface][j].first;
        j--;
        }
    n_peaks[surface] = j+1;

    /* find peak centers */
    for (i=0; i<n_peaks[surface]; i++)
      find_center(strip, surface, &(peaks[surface][i]));

    } /* end for (surface=0; surface<6; surface++) */

  /*---- find z & theta for all strip crossings --------------------*/

  n_z[0] = 0;
  n_z[1] = 0;
  zp_1 = zhit[0];
  zp_2 = zhit[1];

  for (surface=0; surface<6; surface+=2)
    {
    for (i=0; i<n_peaks[surface]; i++)
      for (j=0; j<n_peaks[surface+1]; j++)
        {
        if (surface<4)
          find_z(surface, peaks[surface][i].center, peaks[surface+1][j].center,
                &(n_z[0]), &zp_1);
        else
          find_z(surface, peaks[surface][i].center, peaks[surface+1][j].center,
                &(n_z[1]), &zp_2);
        }
    } /* end for (surface=0; surface<6; surface+=2) */

  /*---- find strip crossings corresponding to anode hits ----------*/

  memset(mwpchits, 0, sizeof(mwpchits));

  for (chamber=0; chamber<2; chamber++)
    {
    for (i=0; i<n_wire_clusters[chamber]; i++)
      {
      vmin = 1000.f;
      for (j=0; j<n_z[chamber]; j++)
        {
        phi_delta = fabs(cluster_center[chamber][i] - zhit[chamber][j].phi);
	// WL
	tmp_diff = cluster_center[chamber][i] - zhit[chamber][j].phi;
        if (phi_delta > 180) 
          phi_delta = 360.f - phi_delta;
        if (phi_delta < vmin)
          {
          jmin = j;
          vmin = phi_delta;
          }
        } /* end for (j=0; j<n_z[chamber]; j++) */

      if (vmin < phi_max)
        {
        mwpchits[chamber][i].phi    = zhit[chamber][jmin].phi;
        mwpchits[chamber][i].z    = zhit[chamber][jmin].z;
        mwpchits[chamber][i].theta  = zhit[chamber][jmin].theta;
        mwpchits[chamber][i].phi_ok = 1;
	HF2(MWPC_ID_BASE+10+chamber, mwpchits[chamber][i].phi, tmp_diff, 1.f);
        } /* end if  vmin < PHI_MAX) */
      else
        {
        mwpchits[chamber][i].phi    = cluster_center[chamber][i];
        mwpchits[chamber][i].z      = INVALID;
        mwpchits[chamber][i].theta  = INVALID;
        mwpchits[chamber][i].phi_ok = 0;
        }
      } /* end for (i=0; i<n_wire_clusters[chamber]; i++) */
    } /* end for (chamber=0; chamber<2; chamber++)*/

  /*---- create decoded MWPC bank ----------------------------------*/

  bk_create(pevent, "MWPC", TID_STRUCT, &mwpc);

  mwpc->wc1_n_hits = n_wire_clusters[0];
  mwpc->wc2_n_hits = n_wire_clusters[1];
  
  for (i=0 ; i<4 ; i++)
    if (mwpchits[0][i].phi_ok) 
      {
      mwpc->wc1_theta[i] = (float) mwpchits[0][i].theta;
      mwpc->wc1_phi[i]   = (float) mwpchits[0][i].phi;
      }
    else
      {
      mwpc->wc1_theta[i] = INVALID;
      mwpc->wc1_phi[i]   = INVALID;
      }

  for (i=0 ; i<4 ; i++)
    if (mwpchits[1][i].phi_ok) 
      {
      mwpc->wc2_theta[i] = (float) mwpchits[1][i].theta;
      mwpc->wc2_phi[i]   = (float) mwpchits[1][i].phi;
      }
    else
      {
      mwpc->wc2_theta[i] = INVALID;
      mwpc->wc2_phi[i]   = INVALID;
      }

  bk_close(pevent, mwpc+1);

  /* fill histos */
  for (i=0 ; i<4 ; i++)
    {
    if (mwpc->wc1_theta[i] != INVALID)
      HF2(MWPC_ID_BASE+4, mwpc->wc1_phi[i], mwpc->wc1_theta[i], 1.f);
    if (mwpc->wc2_theta[i] != INVALID)
      HF2(MWPC_ID_BASE+5, mwpc->wc2_phi[i], mwpc->wc2_theta[i], 1.f);
    }

  return SUCCESS;
}