/*
 * hitFDC - registers hits for forward drift chambers
 *
 *        This is a part of the hits package for the
 *        HDGeant simulation program for Hall D.
 *
 *        version 1.0         -Richard Jones July 16, 2001
 *
 * changes: Wed Jun 20 13:19:56 EDT 2007 B. Zihlmann 
 *          add ipart to the function hitForwardDC
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <HDDM/hddm_s.h>
#include <geant3.h>
#include <bintree.h>

#include "calibDB.h"
extern s_HDDM_t* thisInputEvent;
extern double asic_response(double t);
extern double Ei(double x);

typedef struct{
  int writeenohits;
  int showersincol;
  int driftclusters;
}controlparams_t;

extern controlparams_t controlparams_;


const float wire_dead_zone_radius[4]={3.0,3.0,3.9,3.9};
const float strip_dead_zone_radius[4]={1.3,1.3,1.3,1.3};

#define CATHODE_ROT_ANGLE 1.309 // 75 degrees

// Drift speed 2.2cm/us is appropriate for a 90/10 Argon/Methane mixture
static float DRIFT_SPEED           =.0055;
static float ACTIVE_AREA_OUTER_RADIUS =48.5;
static float ANODE_CATHODE_SPACING =0.5;
static float TWO_HIT_RESOL         =25.;
static int   WIRES_PER_PLANE       =96;
static float WIRE_SPACING          =1.0;
static float U_OF_WIRE_ZERO        =0;//(-((WIRES_PER_PLANE-1.)*WIRE_SPACING)/2)
static float STRIPS_PER_PLANE      =192;
static float STRIP_SPACING         =0.5;
static float U_OF_STRIP_ZERO	   =0;//  (-((STRIPS_PER_PLANE-1.)*STRIP_SPACING)/2)
static float STRIP_GAP             =0.1;
static int   MAX_HITS             =100;
static float K2                  =1.15;
static float STRIP_NODES          = 3;
static float THRESH_KEV           =1. ;
static float THRESH_ANODE         = 1.;
static float THRESH_STRIPS        =5. ;  /* pC */
static float ELECTRON_CHARGE =1.6022e-4; /* fC */
static float DIFFUSION_COEFF  =   1.1e-6; // cm^2/s --> 200 microns at 1 cm
static float FDC_TIME_WINDOW = 1000.0; //time window for accepting FDC hits, ns
static float GAS_GAIN = 8e4;


#if 0
static float wire_dx_offset[2304];
static float wire_dz_offset[2304];
#endif

binTree_t* forwardDCTree = 0;
static int stripCount = 0;
static int wireCount = 0;
static int pointCount = 0;
static int initializedx=0;

void gpoiss_(float*,int*,const int*); // avoid solaris compiler warnings
void rnorml_(float*,int*);

typedef int (*compfn)(const void*, const void*);

// Sort functions for sorting clusters
int fdc_anode_cluster_sort(const void *a,const void *b){
  const s_FdcAnodeTruthHit_t *ca=a;
  const s_FdcAnodeTruthHit_t *cb=b;
  if (ca->t<cb->t) return -1;
  else if (ca->t>cb->t) return 1;
  else return 0;
}
int fdc_cathode_cluster_sort(const void *a,const void *b){
  const s_FdcCathodeTruthHit_t *ca=a;
  const s_FdcCathodeTruthHit_t *cb=b;
  if (ca->t<cb->t) return -1;
  else if (ca->t>cb->t) return 1;
  else return 0;
}



// Locate a position in array xx given x
void locate(float *xx,int n,float x,int *j){
  int ju,jm,jl;
  int ascnd;
  
  jl=-1;
  ju=n;
  ascnd=(xx[n-1]>=xx[0]);
  while(ju-jl>1){
    jm=(ju+jl)>>1;
    if (x>=xx[jm]==ascnd)
      jl=jm;
    else
      ju=jm;
  }
  if (x==xx[0]) *j=0;
  else if (x==xx[n-1]) *j=n-2;
  else *j=jl; 
}

// Polynomial interpolation on a grid.
// Adapted from Numerical Recipes in C (2nd Edition), pp. 121-122.
void polint(float *xa, float *ya,int n,float x, float *y,float *dy){
  int i,m,ns=0;
  float den,dif,dift,ho,hp,w;

  float *c=(float *)calloc(n,sizeof(float));
  float *d=(float *)calloc(n,sizeof(float));

  dif=fabs(x-xa[0]);
  for (i=0;i<n;i++){
    if ((dift=fabs(x-xa[i]))<dif){
      ns=i;
      dif=dift;
    }
    c[i]=ya[i];
    d[i]=ya[i];
  }
  *y=ya[ns--];

  for (m=1;m<n;m++){
    for (i=1;i<=n-m;i++){
      ho=xa[i-1]-x;
      hp=xa[i+m-1]-x;
      w=c[i+1-1]-d[i-1];
      if ((den=ho-hp)==0.0) {
        free(c);
        free(d);
        return;
      }
      
      den=w/den;
      d[i-1]=hp*den;
      c[i-1]=ho*den;
      
    }
    
    *y+=(*dy=(2*ns<(n-m) ?c[ns+1]:d[ns--]));
  }
  free(c);
  free(d);
}

// Simulation of signal on a wire
double wire_signal(double t,s_FdcAnodeTruthHits_t* ahits){
  double t0=1.0; // ns; rough order of magnitude
  int m;
  double asic_gain=0.76; // mV/fC
  double func=0;
  for (m=0;m<ahits->mult;m++){
    if (t>ahits->in[m].t){
      double my_time=t-ahits->in[m].t;
      func+=asic_gain*ahits->in[m].dE*asic_response(my_time);
    }
  }
  return func;
}

// Simulation of signal on a cathode strip (ASIC output)
double cathode_signal(double t,s_FdcCathodeTruthHits_t* chits){
  double t0=1.0; // ns; rough order of magnitude
  int m;
  double asic_gain=2.3;
  double func=0;
  for (m=0;m<chits->mult;m++){
    if (t>chits->in[m].t){
      double my_time=t-chits->in[m].t;
      func+=asic_gain*chits->in[m].q*asic_response(my_time);
    }
  }
  return func;
}

// Generate hits in two cathode planes flanking the wire plane  
void AddFDCCathodeHits(int PackNo,float xwire,float avalanche_y,float tdrift,
		       int n_p,int track,int ipart,int chamber,int module,
		       int layer, int global_wire_number){

  s_FdcCathodeTruthHits_t* chits;    	  

  // Anode charge
  float q_anode;
  int n_t;
  // Average number of secondary ion pairs for 40/60 Ar/CO2 mixture
  float n_s_per_p=1.89; 
  if (controlparams_.driftclusters==0){    
    /* Total number of ion pairs.  On average for each primary ion 
       pair produced there are n_s secondary ion pairs produced.  The
       probability distribution is a compound poisson distribution
       that requires generating two Poisson variables.
    */
    int n_s,one=1;  
    float n_s_mean = ((float)n_p)*n_s_per_p;
    gpoiss_(&n_s_mean,&n_s,&one);
    n_t = n_s+n_p;
    q_anode=((float)n_t)*GAS_GAIN*ELECTRON_CHARGE;
  }
  else{
    // Distribute the number of secondary ionizations for this primary
    // ionization according to a Poisson distribution with mean n_s_over_p.
    // For simplicity we assume these secondary electrons and the primary
    // electron stay together as a cluster.
    int n_s;
    int one=1;
    gpoiss_(&n_s_per_p,&n_s,&one);
    // Anode charge in units of fC
    n_t=1+n_s;
    q_anode=GAS_GAIN*ELECTRON_CHARGE*((float)n_t);
  }

  /* Mock-up of cathode strip charge distribution */ 
  int plane, node;
  for (plane=1; plane<4; plane+=2){
    float theta = (plane == 1)? -CATHODE_ROT_ANGLE : +CATHODE_ROT_ANGLE;
    float cathode_u = xwire*cos(theta)+avalanche_y*sin(theta);
    int strip1 = ceil((cathode_u-U_OF_STRIP_ZERO)/STRIP_SPACING +0.5);
    float cathode_u1 = (strip1-1)*STRIP_SPACING + U_OF_STRIP_ZERO;
    float delta = cathode_u-cathode_u1;
    float half_gap=ANODE_CATHODE_SPACING;

#if 0
    half_gap+=(plane==1)?+wire_dz_offset[global_wire_number]:
	       -wire_dz_offset[global_wire_number];
#endif

    for (node=-STRIP_NODES; node<=STRIP_NODES; node++){
      /* Induce charge on the strips according to the Mathieson 
	 function tuned to results from FDC prototype
      */
      float lambda1=(((float)node-0.5)*STRIP_SPACING+STRIP_GAP/2.
		     -delta)/half_gap;
      float lambda2=(((float)node+0.5)*STRIP_SPACING-STRIP_GAP/2.
		     -delta)/half_gap;
      float factor=0.25*M_PI*K2;
      float q = 0.25*q_anode*(tanh(factor*lambda2)-tanh(factor*lambda1));
      
      int strip = strip1+node;
      /* Throw away hits on strips falling within a certain dead-zone
	 radius */
      float strip_outer_u=cathode_u1
	+(STRIP_SPACING+STRIP_GAP/2.)*(int)node;
      float cathode_v=-xwire*sin(theta)+avalanche_y*cos(theta);
      float check_radius=sqrt(strip_outer_u*strip_outer_u
			      +cathode_v*cathode_v);
      
      if ((strip > 0) 
	  && (check_radius>strip_dead_zone_radius[PackNo]) 
	  && (strip <= STRIPS_PER_PLANE)){
	int mark = (chamber<<20) + (plane<<10) + strip;
	void** cathodeTwig = getTwig(&forwardDCTree, mark);
	if (*cathodeTwig == 0){
	  s_ForwardDC_t* fdc = *cathodeTwig = make_s_ForwardDC();
	  s_FdcChambers_t* chambers = make_s_FdcChambers(1);
	  s_FdcCathodeStrips_t* strips = make_s_FdcCathodeStrips(1);
	  strips->mult = 1;
	  strips->in[0].plane = plane;
	  strips->in[0].strip = strip;
	  strips->in[0].fdcCathodeTruthHits = chits
	    = make_s_FdcCathodeTruthHits(MAX_HITS);
	  chambers->mult = 1;
	  chambers->in[0].module = module;
	  chambers->in[0].layer = layer;
	  chambers->in[0].fdcCathodeStrips = strips;
	  fdc->fdcChambers = chambers;
	  stripCount++;
	}
	else{
	  s_ForwardDC_t* fdc = *cathodeTwig;
	  chits = fdc->fdcChambers->in[0].fdcCathodeStrips
	    ->in[0].fdcCathodeTruthHits;
	}
	
	int nhit;
	for (nhit = 0; nhit < chits->mult; nhit++){
	  // To cut down on the number of output clusters, combine 
	  // those that would be indistiguishable in time given the 
	  // expected timing resolution
	  if (fabs(chits->in[nhit].t - tdrift) <TWO_HIT_RESOL)
	    {
	      break;
	    }
	}
	if (nhit < chits->mult)		/* merge with former hit */
	  {
	    /* Use the time from the earlier hit but add the charge */
	    chits->in[nhit].q += q;
	    if(chits->in[nhit].t>tdrift){
	      chits->in[nhit].t = tdrift;
	      chits->in[nhit].itrack = gidGetId(track);
	      chits->in[nhit].ptype = ipart;
	    }
	  }
	else if (nhit < MAX_HITS){        /* create new hit */
	  chits->in[nhit].t = tdrift;
	  chits->in[nhit].q = q;
	  chits->in[nhit].itrack = gidGetId(track);
	  chits->in[nhit].ptype = ipart;
	  chits->mult++;
	}
	else{
	  // supress warning
	  /*
	    fprintf(stderr,"HDGeant error in hitForwardDC: ");
	    fprintf(stderr,"max hit count %d exceeded, truncating!\n",
	    MAX_HITS);
	  */
	}
	
      }
    } // loop over cathode strips
  } // loop over cathode views
}









// Add wire information
int AddFDCAnodeHit(s_FdcAnodeTruthHits_t* ahits,int layer,int ipart,int track,
		   float xwire,float xyz[3],float dE,float t,float *tdrift){
 
  // Generate 2 random numbers from a Gaussian distribution
  // 
  float rndno[2];
  int two=2;

  // Only and always use the built-in Geant random generator,
  // otherwise debugging is a problem because sequences are not
	  // reproducible from a given pair of random seeds. [rtj]
  
  /* rnorml_(rndno,&two); */ {
    float rho,phi1;
    grndm_(rndno,&two);
    rho = sqrt(-2*log(rndno[0]));
    phi1 = rndno[1]*2*M_PI;
    rndno[0] = rho*cos(phi1);
    rndno[1] = rho*sin(phi1);
  }

   // Get the magnetic field at this cluster position	    
  float x[3],B[3];
  transformCoord(xyz,"local",x,"global");
  gufld_db_(x,B);
  
  // Find the angle between the wire direction and the direction of the
  // magnetic field in the x-y plane
  float wire_dir[2];
  float wire_theta=1.0472*(float)((layer%3)-1);
  float phi=0.;;
  float Br=sqrt(B[0]*B[0]+B[1]*B[1]);
  
  wire_dir[0]=sin(wire_theta);
  wire_dir[1]=cos(wire_theta);
  if (Br>0.) phi= acos((B[0]*wire_dir[0]+B[1]*wire_dir[1])/Br);
  
  // useful combinations of dx and dz
  float dx=xyz[0]-xwire;
  float dx2=dx*dx;
  float dx4=dx2*dx2;
  float dz2=xyz[2]*xyz[2];
  float dz4=dz2*dz2;
  
  // Next compute the avalanche position along wire.  
  // Correct avalanche position with deflection along wire due to 
  // Lorentz force.
  xyz[1]+=( 0.1458*B[2]*(1.-0.048*Br) )*dx
    +( 0.1717+0.01227*B[2] )*(Br*cos(phi))*xyz[2]
    +( -0.000176 )*dx*dx2/(dz2+0.001);
  // Add transverse diffusion
  xyz[1]+=(( 0.01 )*pow(dx2+dz2,0.125)+( 0.0061 )*dx2)*rndno[1];

  // Do not use this cluster if the Lorentz force would deflect 
  // the electrons outside the active region of the detector
  if (sqrt(xyz[1]*xyz[1]+xwire*xwire)>ACTIVE_AREA_OUTER_RADIUS) 
    return 0;

  // Model the drift time and longitudinal diffusion as a function of 
  // position of the cluster within the cell 	 	  
  float tdrift_unsmeared=( 1086.0-106.7*B[2] )*dx2+( 1068.0 )*dz2
    +dx4*(( -2.675 )/(dz2+0.001)+( 2.4e4 )*dz2);	
  float dt=(( 39.44   )*dx4/(0.5-dz2)+( 56.0  )*dz4/(0.5-dx2)
      +( 0.01566 )*dx4/(dz4+0.002)/(0.251-dx2))*rndno[1];
  // Only allow deviations to longer times if near the cell border
  double dradius=sqrt(dx2+dz2);
  if (dradius-0.5>=0. && dt<0){ 
    dt=0.;
  }

  // Minimum drift time for docas near wire (very crude approximation)
  double v_max=0.08; // guess for now based on Garfield, near wire 
  double tmin=dradius/v_max;
  double tdrift_smeared=tdrift_unsmeared+dt;
  if (tdrift_smeared<tmin){
    tdrift_smeared=tmin;
  }

  // Avalanche time
  *tdrift=t+tdrift_smeared;
	  
  // Skip cluster if the time would go beyond readout window
  if (t>FDC_TIME_WINDOW) return 0;

  int nhit;

  // Record the anode hit
  for (nhit = 0; nhit < ahits->mult; nhit++)
    {
      if (fabs(ahits->in[nhit].t - *tdrift) < TWO_HIT_RESOL)
	{
	  break;
	}
    }
  if (nhit < ahits->mult)                 /* merge with former hit */
    {
      /* use the time from the earlier hit but add the energy */
      ahits->in[nhit].dE += dE;
      if(ahits->in[nhit].t>*tdrift){
	ahits->in[nhit].t = *tdrift;
	ahits->in[nhit].t_unsmeared=tdrift_unsmeared;
	ahits->in[nhit].d = sqrt(dx2+dz2);
	
	ahits->in[nhit].itrack = gidGetId(track);
	ahits->in[nhit].ptype = ipart;
      }
    
      /*ahits->in[nhit].t = 
	(ahits->in[nhit].t * ahits->in[nhit].dE + tdrift * dE)
	/ (ahits->in[nhit].dE += dE);
      */
    }
  else if (nhit < MAX_HITS)              /* create new hit */
    {
      ahits->in[nhit].t = *tdrift;
      ahits->in[nhit].t_unsmeared=tdrift_unsmeared;
      ahits->in[nhit].dE = dE;
      ahits->in[nhit].d = sqrt(dx2+dz2);
      ahits->in[nhit].itrack = gidGetId(track);
      ahits->in[nhit].ptype = ipart;
      ahits->mult++;
    }
  else
    {
      fprintf(stderr,"HDGeant error in hitForwardDC: ");
      fprintf(stderr,"max hit count %d exceeded, truncating!\n",MAX_HITS);
    }

  return 1;
}

/* register hits during tracking (from gustep) */

void hitForwardDC (float xin[4], float xout[4],
                   float pin[5], float pout[5], float dEsum,
                   int track, int stack, int history, int ipart)
{
  float x[3], t;
  float dx[3], dr;
  float dEdx;
  float xlocal[3];
  float xinlocal[3];
  float xoutlocal[3];
  float dradius;
  float alpha,sinalpha,cosalpha;
  int i,j;

  if (!initializedx){
      mystr_t strings[250];
      float values[250];
      int nvalues = 250;

      // Get parameters related to the geometry and the signals 
      int status = GetConstants("FDC/fdc_parms", &nvalues,values,strings);
      if (!status) {
        int ncounter = 0;
        int i;
        for ( i=0;i<(int)nvalues;i++){
          //printf("%d %s %f\n",i,strings[i].str,values[i]);
          if (!strcmp(strings[i].str,"FDC_DRIFT_SPEED")) {
            DRIFT_SPEED  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_ACTIVE_AREA_OUTER_RADIUS")) {
            ACTIVE_AREA_OUTER_RADIUS  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_ANODE_CATHODE_SPACING")) {
            ANODE_CATHODE_SPACING  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_TWO_HIT_RESOL")) {
            TWO_HIT_RESOL  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_WIRES_PER_PLANE")) {
            WIRES_PER_PLANE  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_WIRE_SPACING")) {
            WIRE_SPACING  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_STRIPS_PER_PLANE")) {
            STRIPS_PER_PLANE  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_STRIP_SPACING")) {
            STRIP_SPACING  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_STRIP_GAP")) {
            STRIP_GAP  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_MAX_HITS")) {
            MAX_HITS  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_K2")) {
            K2  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_STRIP_NODES")) {
            STRIP_NODES  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_THRESH_KEV")) {
            THRESH_KEV  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_THRESH_STRIPS")) {
            THRESH_STRIPS  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_ELECTRON_CHARGE")) {
            ELECTRON_CHARGE  = values[i];
            ncounter++;
          }
          if (!strcmp(strings[i].str,"FDC_DIFFUSION_COEFF")) {
            DIFFUSION_COEFF  = values[i];
            ncounter++;
          }
	}
	U_OF_WIRE_ZERO     = (-((WIRES_PER_PLANE-1.)*WIRE_SPACING)/2);
	U_OF_STRIP_ZERO	   = (-((STRIPS_PER_PLANE-1.)*STRIP_SPACING)/2);
  
	if (ncounter==16){
          printf("FDC: ALL parameters loaded from Data Base\n");
        } else if (ncounter<16){
          printf("FDC: NOT ALL necessary parameters found in Data Base %d out of 16\n",ncounter);
        } else {
          printf("FDC: SOME parameters found more than once in Data Base\n");
        }       
#if 0
	{
	  int num_values=2304*2;
	  float my_values[2304*2]; 
	  mystr_t my_strings[2304*2];
	  status=GetArrayConstants("FDC/fdc_wire_offsets",&num_values,
				   my_values,my_strings);
	  if (!status){
	    int i;
	    for (i=0;i<num_values;i+=2){
	      int j=i/2;
	      wire_dx_offset[j]=my_values[i];
	      wire_dz_offset[j]=my_values[i+1];
	    }
	  }
	}
#endif
      }
      initializedx = 1 ;
  }

  /* Get chamber information */ 
  int layer = getlayer_wrapper_(); 
  if (layer==0){
    printf("hitFDC: FDC layer number evaluates to zero! THIS SHOULD NEVER HAPPEN! drop this particle.\n");
    return;
  }
  //int module = getmodule_wrapper_();
  //int chamber = (module*10)+layer;
  //int PackNo = (chamber-11)/20;
  int PackNo = getpackage_wrapper_()-1;
  if (PackNo==-1){
     printf("hitFDC: FDC package number evaluates to zero! THIS SHOULD NEVER HAPPEN! drop this particle.\n");
    return;
  }
  int glayer=6*PackNo+layer-1;
  int module = 2*(PackNo)+(layer-1)/3+1;
  int chamber = (module*10)+(layer-1)%3+1;
  int wire1,wire2;
  int wire,dwire;

  // transform layer number into Richard's scheme
  layer=(layer-1)%3+1;

  transformCoord(xin,"global",xinlocal,"local");

  wire1 = ceil((xinlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
  transformCoord(xout,"global",xoutlocal,"local");
  wire2 = ceil((xoutlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
  // Check that wire numbers are not out of range
  if ((wire1>WIRES_PER_PLANE && wire2==WIRES_PER_PLANE) ||
      (wire2>WIRES_PER_PLANE && wire1==WIRES_PER_PLANE)) 
    wire1=wire2=WIRES_PER_PLANE;  
  if ((wire1==0 && wire2 == 1) || (wire1==1 && wire2== 0)){
    wire1=wire2=1;
  }
  // Make sure at least one wire number is valid
  if (wire1>WIRES_PER_PLANE&&wire2>WIRES_PER_PLANE) return;
  if (wire1==0 && wire2==0) return;

  if (wire1>WIRES_PER_PLANE) wire1=wire2;
  else if (wire2>WIRES_PER_PLANE) wire2=wire1;
  if (wire1==0) wire1=wire2;
  else if (wire2==0) wire2=wire1;

  dwire = (wire1 < wire2)? 1 : -1;
  alpha = atan2(xoutlocal[0]-xinlocal[0],xoutlocal[2]-xinlocal[2]);
  sinalpha=sin(alpha);
  cosalpha=cos(alpha);
  xlocal[0] = (xinlocal[0] + xoutlocal[0])/2;
  xlocal[1] = (xinlocal[1] + xoutlocal[1])/2;
  xlocal[2] = (xinlocal[2] + xoutlocal[2])/2;

  wire = ceil((xlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
  x[0] = (xin[0] + xout[0])/2;
  x[1] = (xin[1] + xout[1])/2;
  x[2] = (xin[2] + xout[2])/2;
  t = (xin[3] + xout[3])/2 * 1e9;
  dx[0] = xin[0] - xout[0];
  dx[1] = xin[1] - xout[1];
  dx[2] = xin[2] - xout[2];
  dr = sqrt(dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]);
  if (dr > 1e-3)
  {
    dEdx = dEsum/dr;
  }
  else
  {
    dEdx = 0;
  }

  /* Make a fuzzy boundary around the forward dead region 
   * by killing any track segment whose midpoint is within the boundary */

  if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
      < wire_dead_zone_radius[PackNo])
  {
    return;
  }

  /* post the hit to the truth tree */
 
  if (history == 0)
  {
    int mark = (1<<16) + (chamber<<20) + pointCount;
    void** twig = getTwig(&forwardDCTree, mark);
    if (*twig == 0)
    {
      s_ForwardDC_t* fdc = *twig = make_s_ForwardDC();
      s_FdcChambers_t* chambers = make_s_FdcChambers(1);
      s_FdcTruthPoints_t* points = make_s_FdcTruthPoints(1);
      float xwire = U_OF_WIRE_ZERO + (wire-1)*WIRE_SPACING;
      float u[2];
      u[0] = xinlocal[2];
      u[1] = xinlocal[0]-xwire;
      dradius = fabs(u[1]*cosalpha-u[0]*sinalpha);
      points->mult = 1;
        int a = thisInputEvent->physicsEvents->in[0].reactions->in[0].vertices->in[0].products->mult;
      points->in[0].primary = (stack <= a);
      points->in[0].track = track;
      points->in[0].x = x[0];
      points->in[0].y = x[1];
      points->in[0].z = x[2];
      points->in[0].t = t;
      points->in[0].px = pin[0]*pin[4];
      points->in[0].py = pin[1]*pin[4];
      points->in[0].pz = pin[2]*pin[4];
      points->in[0].E = pin[3];
      points->in[0].dradius = dradius;
      points->in[0].dEdx = dEdx;
      points->in[0].ptype = ipart;
      chambers->mult = 1;
      chambers->in[0].module = module;
      chambers->in[0].layer = layer;
      chambers->in[0].fdcTruthPoints = points;
      fdc->fdcChambers = chambers;
      pointCount++;
    }
  }

  /* post the hit to the hits tree, mark cell as hit */

  if (dEsum > 0)
  {
    float sign=1.; // for dealing with the y-position for tracks crossing two cells

    for (wire=wire1; wire-dwire != wire2; wire+=dwire)
    {
      int valid_hit=1;
      float dE,dt;
      float u[2];
      float x0[3],x1[3];
      float avalanche_y;
      float xwire = U_OF_WIRE_ZERO + (wire-1)*WIRE_SPACING;
      int global_wire_number=96*glayer+wire-1;

#if 0
      xwire+=wire_dx_offset[global_wire_number];
#endif     

      if (wire1==wire2){
	dE=dEsum; 
	x0[0] = xinlocal[0];
	x0[1] = xinlocal[1];
	x0[2] = xinlocal[2]; 
	x1[0] = xoutlocal[0];
	x1[1] = xoutlocal[1];
	x1[2] = xoutlocal[2];
      }
      else{
	x0[0] = xwire-0.5*dwire*WIRE_SPACING;
	x0[1] = xinlocal[1] + (x0[0]-xinlocal[0]+1e-20)*
	  (xoutlocal[1]-xinlocal[1])/(xoutlocal[0]-xinlocal[0]+1e-20);
	x0[2] = xinlocal[2] + (x0[0]-xinlocal[0]+1e-20)*
	  (xoutlocal[2]-xinlocal[2])/(xoutlocal[0]-xinlocal[0]+1e-20);
	if (fabs(x0[2]-xoutlocal[2]) > fabs(xinlocal[2]-xoutlocal[2]))
	  {
	    x0[0] = xinlocal[0];
	    x0[1] = xinlocal[1];
	    x0[2] = xinlocal[2];
	  }
	x1[0] = xwire+0.5*dwire*WIRE_SPACING;
	x1[1] = xinlocal[1] + (x1[0]-xinlocal[0]+1e-20)*
	  (xoutlocal[1]-xinlocal[1])/(xoutlocal[0]-xinlocal[0]+1e-20);
	x1[2] = xinlocal[2] + (x1[0]-xinlocal[0]+1e-20)*
	  (xoutlocal[2]-xinlocal[2])/(xoutlocal[0]-xinlocal[0]+1e-20);
	if (fabs(x1[2]-xinlocal[2]) > fabs(xoutlocal[2]-xinlocal[2]))
	  {
	    x1[0] = xoutlocal[0];
	    x1[1] = xoutlocal[1];
	    x1[2] = xoutlocal[2];
	  }
	dE = dEsum*(x1[2]-x0[2])/(xoutlocal[2]-xinlocal[2]);
      }

      if (dE > 0){
	s_FdcAnodeTruthHits_t* ahits;    

	// Create (or grab) an entry in the tree for the anode wire
	int mark = (chamber<<20) + (2<<10) + wire;
	void** twig = getTwig(&forwardDCTree, mark);
       	
	if (*twig == 0)
	  {
	    s_ForwardDC_t* fdc = *twig = make_s_ForwardDC();
	    s_FdcChambers_t* chambers = make_s_FdcChambers(1);
	    s_FdcAnodeWires_t* wires = make_s_FdcAnodeWires(1);
	    wires->mult = 1;
	    wires->in[0].wire = wire;
	    wires->in[0].fdcAnodeTruthHits = ahits = make_s_FdcAnodeTruthHits(MAX_HITS);
	    chambers->mult = 1;
	    chambers->in[0].module = module;
	    chambers->in[0].layer = layer;
	    chambers->in[0].fdcAnodeWires = wires;
	    fdc->fdcChambers = chambers;
	    wireCount++;          
	  }
	else
	  {
	    s_ForwardDC_t* fdc = *twig;
	    ahits = fdc->fdcChambers->in[0].fdcAnodeWires->in[0].fdcAnodeTruthHits;
	  }
	

	float rndno[2];
	int two=2;
      
	// Find the number of primary ion pairs:
	/* The total number of ion pairs depends on the energy deposition 
	   and the effective average energy to produce a pair, w_eff.
	   On average for each primary ion pair produced there are n_s_per_p 
	   secondary ion pairs produced. 
	*/
	int one=1;
	// Average number of secondary ion pairs for 40/60 Ar/CO2 mixture
	float n_s_per_p=1.89; 
	//Average energy needed to produce an ion pair for 50/50 mixture
	float w_eff=30.2e-9; // GeV
	// Average number of primary ion pairs
	float n_p_mean = dE/w_eff/(1.+n_s_per_p);
	int n_p; // number of primary ion pairs
	gpoiss_(&n_p_mean,&n_p,&one);
   
	// Drift time
	float tdrift=0;

	if (controlparams_.driftclusters==0){
	  float zrange=x1[2]-x0[2];
	  float tany=(x1[1]-x0[1])/zrange;
	  float tanx=(x1[0]-x0[0])/zrange;
	  float dz=ANODE_CATHODE_SPACING-dradius*sign*sinalpha;
#if 0	  
	  dz+=wire_dz_offset[global_wire_number];
#endif
	  xlocal[0]=x0[0]+tanx*dz;
	  if (fabs(xlocal[0]-xwire)>0.5){
	    xlocal[0]=x1[0];
	    xlocal[1]=x1[1];
	    xlocal[2]=x1[2];
	  }
	  else{
	    xlocal[1]=x0[1]+tany*dz;
	    xlocal[2]=x0[2]+dz;
	  }
	  
	  /* If the cluster position is within the wire-deadened region of the 
	     detector, skip this cluster 
	  */
	  if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
	      >=wire_dead_zone_radius[PackNo]){	 
	    if (AddFDCAnodeHit(ahits,layer,ipart,track,xwire,xlocal,dE,t,
			       &tdrift)){
	      AddFDCCathodeHits(PackNo,xwire,xlocal[1],tdrift,n_p,track,ipart,
				chamber,module,layer,global_wire_number);
	    }
	  
	  }
	}
	else{
	  // Loop over the number of primary ion pairs
	  int n;
	  for (n=0;n<n_p;n++){
	    // Generate a cluster at a random position along the path with cell
	    float rndno[2];
	    grndm_(rndno,&two);
	    float dzrand=(x1[2]-x0[2])*rndno[0];
	    // Position of the cluster
	    xlocal[0]=x0[0]+(x1[0]-x0[0])*rndno[0];
	    xlocal[1]=x0[1]+(x1[1]-x0[1])*rndno[0];
	    xlocal[2]=x0[2]+dzrand;
	    /* If the cluster position is within the wire-deadened region of the 
	       detector, skip this cluster 
	    */
	    if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
		>=wire_dead_zone_radius[PackNo]){	   
	      if (AddFDCAnodeHit(ahits,layer,ipart,track,xwire,xlocal,dE,t,
				 &tdrift)){
		AddFDCCathodeHits(PackNo,xwire,xlocal[1],tdrift,n_p,track,ipart,
				  chamber,module,layer,global_wire_number);
	      }
	    }

	  } // loop over primary ion pairs 
	}
      } // Check for non-zero energy

      sign*=-1; // for dealing with the y-position for tracks crossing two cells
    } // loop over wires
  } // Check that total energy deposition is not zero
}

/* entry points from fortran */

void hitforwarddc_(float* xin, float* xout,
                   float* pin, float* pout, float* dEsum,
                   int* track, int* stack, int* history, int* ipart)
{
   hitForwardDC(xin,xout,pin,pout,*dEsum,*track,*stack,*history,*ipart);
}


/* pick and package the hits for shipping */

s_ForwardDC_t* pickForwardDC ()
{
   s_ForwardDC_t* box;
   s_ForwardDC_t* item;

   if ((stripCount == 0) && (wireCount == 0) && (pointCount == 0))
   {
      return HDDM_NULL;
   }

   box = make_s_ForwardDC();
   box->fdcChambers = make_s_FdcChambers(32);
   box->fdcChambers->mult = 0;
   while (item = (s_ForwardDC_t*) pickTwig(&forwardDCTree))
   {
      s_FdcChambers_t* chambers = item->fdcChambers;
      int module = chambers->in[0].module;
      int layer = chambers->in[0].layer;
      int m = box->fdcChambers->mult;

      /* compress out the hits below threshold */
      s_FdcAnodeWires_t* wires = chambers->in[0].fdcAnodeWires;
      int wire;
      s_FdcCathodeStrips_t* strips = chambers->in[0].fdcCathodeStrips;
      int strip;
      s_FdcTruthPoints_t* points = chambers->in[0].fdcTruthPoints;
      int point;
      int mok=0;
      for (wire=0; wire < wires->mult; wire++)
      {
         s_FdcAnodeTruthHits_t* ahits = wires->in[wire].fdcAnodeTruthHits;
	   
	 // Sort the clusters by time
	 qsort(ahits->in,ahits->mult,sizeof(s_FdcAnodeTruthHit_t),
	       (compfn)fdc_anode_cluster_sort);
	 
	 int i,iok=0;

	 if (controlparams_.driftclusters==0){
	   for (iok=i=0; i < ahits->mult; i++)
	     {
	       if (ahits->in[i].dE >= THRESH_KEV/1e6)
		 {
               if (iok < i)
		 {
		   ahits->in[iok] = ahits->in[i];
		 }
               ++iok;
               ++mok;
		 }
	     }
	 }
	 else{  // Simulate clusters within the cell
	
	   // printf("-------------\n");
	  
	   
	   // Temporary histogram in 1 ns bins to store waveform data
	   int num_samples=(int)FDC_TIME_WINDOW;
	   float *samples=(float *)malloc(num_samples*sizeof(float));
	   for (i=0;i<num_samples;i++){
	     samples[i]=wire_signal((float)i,ahits);
	     //printf("%f %f\n",(float)i,samples[i]);
	   }
	   
	   int returned_to_baseline=0;
	   float q=0;
	   for (i=0;i<num_samples;i++){
	     if (samples[i]>=THRESH_ANODE){
	       if (returned_to_baseline==0){
		 ahits->in[iok].itrack = ahits->in[0].itrack;
		 ahits->in[iok].ptype = ahits->in[0].ptype;
		 
		 // Do an interpolation to find the time at which the threshold
		 // was crossed.
		 float t_array[4];
		 int k;
		 float my_t,my_terr;
		 for (k=0;k<4;k++) t_array[k]=i-1+k;
		 polint(&samples[i-1],t_array,4,THRESH_ANODE,&my_t,&my_terr);
		 ahits->in[iok].t=my_t;
		 
		 returned_to_baseline=1;
		 iok++;
		 mok++;
	       }
	       q+=samples[i];
	     }
	     if (returned_to_baseline 
		 && (samples[i]<THRESH_ANODE)){
	       returned_to_baseline=0;   
	     }
	   }
	   free(samples);
	 } // Simulation of clusters within cell

	 if (iok)
	   {
	     ahits->mult = iok;
	   }
	 else if (ahits != HDDM_NULL)
	   {
	     FREE(ahits);
	   }
      } 
      if ((wires != HDDM_NULL) && (mok == 0))
	{
	  FREE(wires);
	  wires = HDDM_NULL;
	}
      
      
      mok = 0;
      for (strip=0; strip < strips->mult; strip++)
	{
	  s_FdcCathodeTruthHits_t* chits = strips->in[strip].fdcCathodeTruthHits;
	  
	  // Sort the clusters by time
	  qsort(chits->in,chits->mult,sizeof(s_FdcCathodeTruthHit_t),
		  (compfn)fdc_cathode_cluster_sort);
	  
	  int i,iok=0;
	  
	  if (controlparams_.driftclusters==0){
	    for (iok=i=0; i < chits->mult; i++)
	      {
		if (chits->in[i].q > 0.)
		  {
		    if (iok < i)
		      {
			chits->in[iok] = chits->in[i];
		      }
		    ++iok;
		    ++mok;
		  }
	      }	    
	    
	  }
	  else{
	   
	    // Temporary histogram in 1 ns bins to store waveform data
	    int num_samples=(int)(FDC_TIME_WINDOW);
	    float *samples=(float *)malloc(num_samples*sizeof(float));
	    for (i=0;i<num_samples;i++){
	      samples[i]=cathode_signal((float)i,chits);
		 //printf("t %f V %f\n",(float)i,samples[i]);
	    }	 
	    
	    int threshold_toggle=0;
	    int istart=0;
	    float q=0;
	    int FADC_BIN_SIZE=1;
	    for (i=0;i<num_samples;i+=FADC_BIN_SIZE){
	      if (samples[i]>=THRESH_STRIPS){
		if (threshold_toggle==0){
		  chits->in[iok].itrack = chits->in[0].itrack;
		  chits->in[iok].ptype = chits->in[0].ptype;
		  chits->in[iok].t=(float) i;
		  //chits->in[iok].q=samples[i];
		  istart=i-1;
		  threshold_toggle=1;
		  //iok++;
		  //mok++;
		}
	      }
	      if (threshold_toggle && 
		  (samples[i]<THRESH_STRIPS)){
		int j;
		// Find the first peak
		for (j=istart+1;j<i-1;j++){
		  if (samples[j]>samples[j-1] && samples[j]>samples[j+1]){
		    chits->in[iok].q=samples[j];
		    break;
		  }
		} 
		threshold_toggle=0; 
		iok++;
		mok++;
		//break;
	      }
	    }
	    i=num_samples-1;
	    if (samples[i]>=THRESH_STRIPS&&threshold_toggle){
	      int j;
	      for (j=istart+1;j<i-1;j++){
		if (samples[j]>samples[j-1] && samples[j]>samples[j+1]){
		  chits->in[iok].q=samples[j];
		  break;
		}
	      }
	    }
	    
	    free(samples);
	  }// Simulate clusters within cell
	
	  if (iok)
	    {
	      chits->mult = iok;
	      //chits->mult=1;
	    }
	  else if (chits != HDDM_NULL)
	    {
	      FREE(chits);
	    }
	  
	}
      if ((strips != HDDM_NULL) && (mok == 0))
	{
	  FREE(strips);
	  strips = HDDM_NULL;
	}
      
      if ((wires != HDDM_NULL) || 
          (strips != HDDM_NULL) ||
          (points != HDDM_NULL))
	{
	  if ((m == 0) || (module > box->fdcChambers->in[m-1].module)
                     || (layer  > box->fdcChambers->in[m-1].layer))
	    {
          box->fdcChambers->in[m] = chambers->in[0];

          box->fdcChambers->in[m].fdcCathodeStrips = 
                              make_s_FdcCathodeStrips(stripCount);
          box->fdcChambers->in[m].fdcAnodeWires =
                              make_s_FdcAnodeWires(wireCount);
          box->fdcChambers->in[m].fdcTruthPoints =
                              make_s_FdcTruthPoints(pointCount);
          box->fdcChambers->mult++;
        }
        else
        {
          m--;
        }
        for (strip=0; strip < strips->mult; ++strip)
        {
           int mm = box->fdcChambers->in[m].fdcCathodeStrips->mult++;
           box->fdcChambers->in[m].fdcCathodeStrips->in[mm] = strips->in[strip];
        }
        if (strips != HDDM_NULL)
        {
           FREE(strips);
        }
        for (wire=0; wire < wires->mult; ++wire)
        {
           int mm = box->fdcChambers->in[m].fdcAnodeWires->mult++;
           box->fdcChambers->in[m].fdcAnodeWires->in[mm] = wires->in[wire];
        }
        if (wires != HDDM_NULL)
        {
           FREE(wires);
        }
        for (point=0; point < points->mult; ++point)
        {
           int mm = box->fdcChambers->in[m].fdcTruthPoints->mult++;
           box->fdcChambers->in[m].fdcTruthPoints->in[mm] = points->in[point];
        }
        if (points != HDDM_NULL)
        {
           FREE(points);
        }
     }
     FREE(chambers);
     FREE(item);
   }

   stripCount = wireCount = pointCount = 0;

   if ((box->fdcChambers != HDDM_NULL) &&
       (box->fdcChambers->mult == 0))
   {
      FREE(box->fdcChambers);
      box->fdcChambers = HDDM_NULL;
   }
   if (box->fdcChambers->mult == 0)
   {
      FREE(box);
      box = HDDM_NULL;
   }
   return box;
}