#include <iostream>
#include "PAM.h"
#include "PHNode.h"
#include "PHCompositeNode.h"
#include "PHIODataNode.h"
#include "PHNodeIterator.h"
#include "PHTable.hh"
#include "mMutXYtoFEE.hh"

#include "mumhitsWrapper.h"

#include "mumgeoWrapper.h"

#include "mumdigiparWrapper.h"

#include "mMutCathodePulseParams.hh"

#include "dMutDeadChannelsWrapper.h"

#include "dMutCathodeChargeWrapper.h"

#include "dMutCathodesFEMWrapper.h"

#include "dMutCathodesFEMMCWrapper.h"

#include "dMutAnodesFEMWrapper.h"

#include "dMutAnodesFEMMCWrapper.h"

#include "calibData.h"

#include "MutGeom.h"

#include <math.h>
#include <string.h>
#include "gconst.h"
#include "PSTORE.h"
#include "sublink.h"
#include "fpmlink.h"
#include "fpnlink.h"
#include "MUMPAR.h"
#include "COUNTER.h"
#include "dMutArrayParams.h"
#include "dMutFEMPar.h"
#include "TRandom.h"

// Prototypes for functions called:
#include "mMutDeadChannelP.h"
#include "mMutFillCathodeFEMTablesP.h"
#include "mMutCathStripP.h"
#include "bfieldP.h"
#include "gufld.h"
#include "mMutCalcADCSamplesP.h"

extern "C" {

        int nint(float x);
	float ranf_(void); 
	float ranlan_(float *x);
	void vzero(long ar[],int no_element);
}

extern calibData *calib;

#ifdef CXX_MISSING_BOOL
#define bool short
#define TRUE 1
#define true TRUE
#define FALSE 0
#define false FALSE
#endif
typedef PHIODataNode<PHTable> TableNode_t;
typedef PHIODataNode<TObject> ObjectNode_t;

// Default constructor and destructor to pacify CINT
mMutXYtoFEE::mMutXYtoFEE(){}

mMutXYtoFEE::~mMutXYtoFEE(){}

PHBoolean
mMutXYtoFEE::event(PHCompositeNode *root) {
 PHPointerList<PHNode> nodes;
 PHNodeIterator iiter(root), *jiter;
 ObjectNode_t *nn;
 TableNode_t *d;
 PHCompositeNode *dstNode, *outNode;

  mumhitsWrapper*            tmumhitsWrapper;
  TABLE_HEAD_ST              mumhits_h;
  MUMHITS_ST*                mumhits;

  mumgeoWrapper*             tmumgeoWrapper;
  TABLE_HEAD_ST              mumgeo_h;
  MUMGEO_ST*                 mumgeo;

  mumdigiparWrapper*         tmumdigiparWrapper;
  TABLE_HEAD_ST              mumdigipar_h;
  MUMDIGIPAR_ST*             mumdigipar;

  dMutDeadChannelsWrapper*   tdMutDeadChannelsWrapper;
  TABLE_HEAD_ST              dMutDeadChannels_h;
  DMUTDEADCHANNELS_ST*       dMutDeadChannels;

  dMutCathodeChargeWrapper*  tdMutCathodeChargeWrapper;
  TABLE_HEAD_ST              dMutCathodeCharge_h;
  DMUTCATHODECHARGE_ST*      dMutCathodeCharge;

  dMutCathodesFEMWrapper*    tdMutCathodesFEMWrapper;
  TABLE_HEAD_ST              dMutCathodesFEM_h;
  DMUTCATHODESFEM_ST*        dMutCathodesFEM;

  dMutCathodesFEMMCWrapper*  tdMutCathodesFEMMCWrapper;
  TABLE_HEAD_ST              dMutCathodesFEMMC_h;
  DMUTCATHODESFEMMC_ST*      dMutCathodesFEMMC;

  dMutAnodesFEMWrapper*      tdMutAnodesFEMWrapper;
  TABLE_HEAD_ST              dMutAnodesFEM_h;
  DMUTANODESFEM_ST*          dMutAnodesFEM;

  dMutAnodesFEMMCWrapper*    tdMutAnodesFEMMCWrapper;
  TABLE_HEAD_ST              dMutAnodesFEMMC_h;
  DMUTANODESFEMMC_ST*        dMutAnodesFEMMC;

 dstNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "DST"));
 if (!dstNode) {
   dstNode = new PHCompositeNode("DST");
   root->addNode(dstNode);
 }

// Insert code here to navigate node hierarchy and find
// or create specific nodes to pass to physics module...

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "GEA"));
  if (!outNode) {
    outNode = new PHCompositeNode("GEA");
    root->addNode(outNode);
  }

  d = static_cast<TableNode_t*>(iiter.findFirst("PHIODataNode", "mumhits"));
  if (d) {
    tmumhitsWrapper = static_cast<mumhitsWrapper*>(d->getData());
  }
  else{
    cout << "ERROR:  'in' parameter mumhits not found" << endl;
    return False;
  }
  mumhits_h = tmumhitsWrapper->TableHeader();
  mumhits   = tmumhitsWrapper->TableData();


  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "PAR"));
  if (!outNode) {
    outNode = new PHCompositeNode("PAR");
    root->addNode(outNode);
  }

  d = static_cast<TableNode_t*>(iiter.findFirst("PHIODataNode", "mumgeo"));
  if (d) {
    tmumgeoWrapper = static_cast<mumgeoWrapper*>(d->getData());
    mumgeo_h = tmumgeoWrapper->TableHeader();
    mumgeo   = tmumgeoWrapper->TableData();
  }
  else{
    cout << "ERROR:  'in' parameter mumgeo not found" << endl;
    mumgeo = 0;
  }


  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "PAR"));
  if (!outNode) {
    outNode = new PHCompositeNode("PAR");
    root->addNode(outNode);
  }

  d = static_cast<TableNode_t*>(iiter.findFirst("PHIODataNode", "mumdigipar"));
  if (d) {
    tmumdigiparWrapper = static_cast<mumdigiparWrapper*>(d->getData());
  }
  else{
    cout << "ERROR:  'in' parameter mumdigipar not found" << endl;
    return False;
  }
  mumdigipar_h = tmumdigiparWrapper->TableHeader();
  mumdigipar   = tmumdigiparWrapper->TableData();

  mMutCathodePulseParams* dMutCathodePulseParams;
  nn = static_cast<ObjectNode_t*>(iiter.findFirst("PHIODataNode","dMutCathodePulseParams"));
  if (nn){
    dMutCathodePulseParams = static_cast<mMutCathodePulseParams*>(nn->getData());
  }  
  else{
    cout << "ERROR:  'in' parameter dMutCathodePulseParams not found" << endl;
    return False;
  }

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "PAR"));
  if (!outNode) {
    outNode = new PHCompositeNode("PAR");
    root->addNode(outNode);
  }

  d = static_cast<TableNode_t*>(iiter.findFirst("PHIODataNode", "dMutDeadChannels"));
  if (d) {
    tdMutDeadChannelsWrapper = static_cast<dMutDeadChannelsWrapper*>(d->getData());
  }
  else{
    cout << "ERROR:  'in' parameter dMutDeadChannels not found" << endl;
    return False;
  }
  dMutDeadChannels_h = tdMutDeadChannelsWrapper->TableHeader();
  dMutDeadChannels   = tdMutDeadChannelsWrapper->TableData();

  outNode = dstNode;
  jiter = new PHNodeIterator(outNode);
  d = static_cast<TableNode_t*>(jiter->findFirst("PHIODataNode","dMutCathodeCharge"));
  if (d){
    tdMutCathodeChargeWrapper = static_cast<dMutCathodeChargeWrapper*>(d->getData());
  }
  else{
     tdMutCathodeChargeWrapper = new dMutCathodeChargeWrapper("dMutCathodeCharge", 5000);
     if (!tdMutCathodeChargeWrapper) {
       return 1;
     }
     d = new TableNode_t(tdMutCathodeChargeWrapper,"dMutCathodeCharge");
     outNode->addNode(d);
  }
  dMutCathodeCharge_h = tdMutCathodeChargeWrapper->TableHeader();
  dMutCathodeCharge   = tdMutCathodeChargeWrapper->TableData();
  delete jiter;

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "MUT"));
  if (!outNode) {
    outNode = new PHCompositeNode("MUT");
    root->addNode(outNode);
  }

  jiter = new PHNodeIterator(outNode);
  d = static_cast<TableNode_t*>(jiter->findFirst("PHIODataNode","dMutCathodesFEM"));
  if (d){
    tdMutCathodesFEMWrapper = static_cast<dMutCathodesFEMWrapper*>(d->getData());
  }
  else{
     tdMutCathodesFEMWrapper = new dMutCathodesFEMWrapper("dMutCathodesFEM", 510);
     if (!tdMutCathodesFEMWrapper) {
       return 1;
     }
     d = new TableNode_t(tdMutCathodesFEMWrapper,"dMutCathodesFEM");
     outNode->addNode(d);
  }
  dMutCathodesFEM_h = tdMutCathodesFEMWrapper->TableHeader();
  dMutCathodesFEM   = tdMutCathodesFEMWrapper->TableData();
  delete jiter;

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "EVA"));
  if (!outNode) {
    outNode = new PHCompositeNode("EVA");
    root->addNode(outNode);
  }

  jiter = new PHNodeIterator(outNode);
  d = static_cast<TableNode_t*>(jiter->findFirst("PHIODataNode","dMutCathodesFEMMC")) ;
  if (d){
    tdMutCathodesFEMMCWrapper = static_cast<dMutCathodesFEMMCWrapper*>(d->getData());
  }
  else{
     tdMutCathodesFEMMCWrapper = new dMutCathodesFEMMCWrapper("dMutCathodesFEMMC", 510);
     if (!tdMutCathodesFEMMCWrapper) {
       return 1;
     }
     d = new TableNode_t(tdMutCathodesFEMMCWrapper,"dMutCathodesFEMMC");
     outNode->addNode(d);
  }
  dMutCathodesFEMMC_h = tdMutCathodesFEMMCWrapper->TableHeader();
  dMutCathodesFEMMC   = tdMutCathodesFEMMCWrapper->TableData();
  delete jiter;

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "MUT"));
  if (!outNode) {
    outNode = new PHCompositeNode("MUT");
    root->addNode(outNode);
  }

  jiter = new PHNodeIterator(outNode);
  d = static_cast<TableNode_t*>(jiter->findFirst("PHIODataNode","dMutAnodesFEM"));
  if (d){
    tdMutAnodesFEMWrapper = static_cast<dMutAnodesFEMWrapper*>(d->getData());
  }
  else{
     tdMutAnodesFEMWrapper = new dMutAnodesFEMWrapper("dMutAnodesFEM", 20);
     if (!tdMutAnodesFEMWrapper) {
       return 1;
     }
     d = new TableNode_t(tdMutAnodesFEMWrapper,"dMutAnodesFEM");
     outNode->addNode(d);
  }
  dMutAnodesFEM_h = tdMutAnodesFEMWrapper->TableHeader();
  dMutAnodesFEM   = tdMutAnodesFEMWrapper->TableData();
  delete jiter;

  outNode = static_cast<PHCompositeNode*>(iiter.findFirst("PHCompositeNode", "EVA"));
  if (!outNode) {
    outNode = new PHCompositeNode("EVA");
    root->addNode(outNode);
  }

  jiter = new PHNodeIterator(outNode);
  d = static_cast<TableNode_t*>(jiter->findFirst("PHIODataNode","dMutAnodesFEMMC"));
  if (d){
    tdMutAnodesFEMMCWrapper = static_cast<dMutAnodesFEMMCWrapper*>(d->getData());
  }
  else{
     tdMutAnodesFEMMCWrapper = new dMutAnodesFEMMCWrapper("dMutAnodesFEMMC", 20);
     if (!tdMutAnodesFEMMCWrapper) {
       return 1;
     }
     d = new TableNode_t(tdMutAnodesFEMMCWrapper,"dMutAnodesFEMMC");
     outNode->addNode(d);
  }
  dMutAnodesFEMMC_h = tdMutAnodesFEMMCWrapper->TableHeader();
  dMutAnodesFEMMC   = tdMutAnodesFEMMCWrapper->TableData();
  delete jiter;

/*:>--------------------------------------------------------------------
**: FILE:       mmutxytofee.c.template
**: HISTORY:
**:             00jan93-v000a-hpl- Created by stic Version
**:  Id: idl.y,v 1.1.1.1 1997/09/26 20:43:47 dave Exp  
**:<------------------------------------------------------------------*/

/**
    DESCRIPTION:  Turn input PISA (x,y) hits into 4 adc samples from
      each cathode strip and charge on each anode wire in the muon 
      tracking system.  Noise should be added to the strips and wires
      and pedestals and gains put into the ADC values.  Landau
      distribution is used to determine the total charge deposited for
      each track and the Mathieson distribution is used to determine
      the charge distribution on the cathode strips.   <p>
      
    @param   <p>
             IN:   <p>
             mumhits     - Input PISA hits table for muon tracking   <p>
           mumhits_h     -    <p>
              mumgeo     - Input PISA geometry parameters for muon tr.   <p>
            mumgeo_h     -    <p>
          mumdigipar     - Input muon tracking digitization parameters   <p>
        mumdigipar_h     -    <p>
          INOUT:   <p>
            OUT:   <p>
     dMutCathodesFEM     - cathode info, in FEM packet format   <p>
     dMutCathodesFEM_h   -    <p>
     dMutAnodesFEM       - anode info, in FEM packet format   <p>
     dMutAnodesFEM_h     -    <p>

    @return     STAF Condition Value   <p>

    @author     Melynda Brooks, LANL, mbrooks@lanl.gov
*/


#define	NOISE 1
#define DEGRAD 0.01745329
#define PIBY2 1.570796
#define TO_RAD(dgree)   ( (dgree)* M_PI/180.0 )
#define TO_DGREE(rad)   ( (rad) *180.0 /M_PI )

#define acosd(x)      	TO_DGREE(acos( (x) ) )
#define atan2d(x,y)     TO_DGREE(atan2((x), (y)))


//	digitization counters
  long int iarm, istation, iplane, isegment, ihalf;
  long int igap, icath;
  Int_t NumberOfGaps, NumberOfStrips;
  Int_t WireNoMin, WireNoMax;

//	values retrieved from pisa bank
  long int ipid, iposition, itrack, mmul;
  float ptof, energy, xpos, ypos, zpos, xmom, ymom, zmom;
  float rpos, theta, phi;

//	variables used for digitization
  long int j, iplane1;
  long int nhits, ihit;
  long int truewireno;
  long int ntracks, datastart;

  float trise, tfall, tspace, t0, tsamples[4], vpeak;
  float pedestal, gain, rms;
  float chargesamp[4];
  long int adcsamp[4];
  int nsamples, goodadc, isamp;
  

// variables to include error from cathode strip stereo angle
  float Stereo[2]; /* angle of Cathode plane 1 strip w.r.t. perpendicular */
  float stererr; /* shift of charge along strip from stereo angle */

// variables to give position error to non-normal tracks
  float octang;
  float anmom;  /* component of momentum parallel to anode wire */
  float phian;  /* angle from normal */
  float phierr; /* error in hit position from non-normal incidence */
#define DEG45 0.785398163
#define RADDEG 57.29577951

// variables to calculate position error due to Lorentz angle:
  float rad;
  float langle; /* Lorentz angle for this position */  
  float lorerr; /* Equivalent error in position */
  float br, bz;  /* variables for calling bfield */
  float xyz[3], bxyz[3];
  long order, mcode;

#define	sqrt12	3.464101
  float	ranx, r;
  float	localx,dedx;
  float TrueWireHit[2];
  float	padwidth;

  static long int nmmul[MAX_FEM_TRACKS][700][2][3][8][2][3][3];
  static long int ntrks[700][2][3][8][2][3][3];
  static float truehit[MAX_FEM_TRACKS][700][2][3][8][2][3][3];

// Electronics information required for packing into FEM packet:

   static short CathFemTot;    // Total number of cathodes FEMs in system
   short AnFemMin;      // Min. FEM channel to zero
   short CathFemMin;    // Min. FEM channel to zero
   short data;          // Variable to use for zeroing FEM data
   int goodfill;

   TRandom* RanGen = new TRandom;	// Random number generator class

static int fitv0=1;    // Argument passed to mMutCalcADCSamples which says: determine
                       // the pulse-shape parameter V0 from the peak charge of the strip

// Geometry class variables:

   MutArm *ArmPointer;
   MutStrip *ptrStrip[2];
   MutWire *ptrWire[1];
   double stripIP[2];
   double wireIP;
   PHPoint hit;

   float CathodeAngle[2], AnodeAngle;
   float ReadoutSpacing;

   Int_t octant;
   int PacketID;
   int DCMChannel;
   int CathodePlane;

   int GoodConvert;

// ***********************************************************************  

/*	Initialization	*/
	

  dMutCathodePulseParams->getTspace(tspace);
  dMutCathodePulseParams->getTrise(trise);
  dMutCathodePulseParams->getTfall(tfall);
  dMutCathodePulseParams->getNsamples(nsamples);
  dMutCathodePulseParams->getTsamples(tsamples);

   nhits = 0;
   vzero(mtdigi.nmuidhits, 50);
   for ( iplane = 0; iplane < 3; iplane++)
   {
     for ( igap = 0; igap < 3; igap++)
     {
       for ( ihalf = 0; ihalf < 2 ; ihalf ++)
       {
         for ( isegment = 0; isegment < mum_segments_max; isegment++)
         {
   	    for ( istation = 0; istation < mum_stations_max; istation++)
  	    {
  	      for ( iarm = 0; iarm < mum_arms_max; iarm++)
	      {
	         for ( j = 0; j < 700; j++)
	         {
	  	    mtdigi.readout_strip[j][iarm][istation][isegment][ihalf][igap][iplane]=0.0;
	  	    if (iplane < 9  &&  isegment < 8  &&  istation < 3 &&  iarm < 2)
		    {
		     ntrks[j][iarm][istation][isegment][ihalf][igap][iplane]=0;
		    }
		    if (j < 11)mtdigi.strip[j][iarm][istation][isegment][ihalf][igap][iplane]=0.0;
	         }
	       } // arms
	    }  // stations
	 }  // octants
       }  // half octants
     }  // gaps
   } // planes within a gap

//  loop over all the muon hits, and extract the ones from the muon tracking
//  station detectors	


   mmul = mumhits_h.nok;

   for ( ihit = 0; ihit < mmul; ihit++)
   {
     iposition = mumhits[ihit].plane;

//  If the octant number is 0, this is not a tracking station hit
//  so go to next hit	

     isegment = (iposition%1000)/100 -1;
     if (isegment  >=  0)
     {
       iarm = iposition/1000 -1;
       istation = (iposition%100)/10 -1;

//  See if this is a muon tracking detector (not needed anymore?)

       if ((istation  >=  0)  &&  (istation  <  NumberOfStations))
       {
	 iplane1 = (int)  iposition%10 -1;

//  Option to remove all hits from certain gaps:

	 if (mumdigipar[0].keepgap[(istation)*3+iplane1]  ==  0) continue;

//  Now get all of the hit information from the bank

	 itrack = mumhits[ihit].track-1;
	 ipid = mumhits[ihit].pid;
	 ptof = mumhits[ihit].t;
	 energy = mumhits[ihit].e;
	 xpos = mumhits[ihit].x[0];
	 ypos = mumhits[ihit].x[1];
	 zpos = mumhits[ihit].x[2];
	 xmom = mumhits[ihit].p[0];
	 ymom = mumhits[ihit].p[1];
	 zmom = mumhits[ihit].p[2];
	 rpos = sqrt(xpos*xpos + ypos*ypos + zpos*zpos);
	 theta = acosd(fabs(zpos)/rpos);
	 phi = atan2d(ypos, xpos);
	 if (phi  <  0) phi = phi + 360.;

//  Skip this hit if it is outside the active area:	

	 if (theta  <  mumdigipar->st_angacc[iarm*3 + istation])  continue;

// perform digitization
// turn lab coordinate into a hit in units of wire numbers (or whatever
// is appropriate)  Put Landau in.  Check value of ranf so that floating point
// or MATH library error does not occur.	

	 r = ranf_();
	 if (r  !=  0.0  &&  r  !=  1.0) dedx = ranlan_(&r)*energy_loss + 165;
	 else
	 {
	   r = ranf_();
	   dedx = ranlan_(&r)*energy_loss + 165;
	 }

// Put a detector efficiency into the hit calculation.  When a chamber
// is inefficient, all three detector planes in the gap will be without
// a signal.

         r = ranf_();
         if (r > mumdigipar[0].hit_efficiency) dedx = 0.0;

	 dMutCathodeCharge[ihit].truecharge = dedx;

// Geometry class:

         hit.setX((double)xpos);
         hit.setY((double)ypos);
         hit.setZ((double)zpos);
         if (zpos<0){
           GoodConvert = SouthArm()->convertPisaHitToChamber(hit, 
                              ptrWire, wireIP, ptrStrip, stripIP);
         }
         else{
           GoodConvert = NorthArm()->convertPisaHitToChamber(hit, 
                              ptrWire, wireIP, ptrStrip, stripIP);
         }   

         if (ptrStrip[0]){
 
// Calculate error in the hit position for tracks that are not perpendicular 
// to the anode wire:  Note that for small angles there should actually
// be a slight IMPROVEMENT in resolution because more charge is deposited in
// the chamber giving a better signal to noise and the spread of the
// charge is still small.  For now, just add if there is a degradation.
// Parameterization was obtained from standalone CSC simulation code.
// Also note, the Lorentz and non-normal error should be the same for
// both cathodes in a gap, adjusted for stereo angles.

           octant = ptrStrip[Cathode_1]->getOctant();
           octang = (float)(octant)*DEG45;
           anmom = ymom*(float)cos((double)octang) - 
                 xmom*(float)sin((double)octang);
           if (zmom != 0) 
                 phian = (float)atan2((double)anmom, (double)fabs(zmom));
           else phian = 0.0;
           phian = fabs(phian*RADDEG);
           phierr = (-2.1*phian + 0.457*((float)pow((double)phian,2)))*1.0e-4;
           if (phierr<0) phierr = 0.0;
           ranx = RanGen->Gaus(0.0,1.0); 
           phierr*=ranx;

// Calculate Lorentz angle error:

           rad = (float)sqrt(pow((double)xpos,2) + pow((double)ypos,2));
           order = 1;
           mcode = 0;
           /* bfield_(&rad, &zpos, &br, &bz, &mcode, &order); */
	   xyz[0] = xpos;
	   xyz[1] = ypos;
	   xyz[2] = zpos;
	   gufld_(xyz, bxyz);
	   bz = bxyz[2];
	   /*
	   cout << "\n mMutXYtoFEE:";
	   cout << "  xpos = " << xpos;
	   cout << ", ypos = " << ypos;
	   cout << ", zpos = " << zpos << endl;
	   cout << "\n   bx = " << bxyz[0] << ", by = " << bxyz[1];
           cout << ",  bz = " << bz << endl;
	   exit(99);
	   */
           langle = bz*mumdigipar[0].baselor/1000.0;
           lorerr = (0.36*langle + 0.44*(float)pow((double)langle,2))*1.0E-4;
           ranx = RanGen->Gaus(0.0,1.0); 
           lorerr *= ranx;

           AnodeAngle = ptrWire[0]->getAngle();
           if (AnodeAngle < 0) AnodeAngle += 2.0*Pi;


// Loop over cathode planes, add position errors to calculated strip positions,
// turn hit positions into charge on strips:

           for (icath=0; icath<NumberOfCathodePlanes; icath++){
             if (ptrStrip[icath]){
               ReadoutSpacing = (float)(ptrStrip[icath]->getStripSpacing());
               TrueWireHit[icath] = (float)(ptrStrip[icath]->getStrip() + 
                                  stripIP[icath]/ReadoutSpacing);

               CathodeAngle[icath] = (float)( ptrStrip[icath]->getAngle() );
               if (CathodeAngle[icath] < 0) CathodeAngle[icath] += 2.0*Pi;

               Stereo[icath] = CathodeAngle[icath] + PIBY2 - AnodeAngle ;
               if (Stereo[icath] > PIBY2) Stereo[icath] -= 2.0*Pi;
               if (Stereo[icath] < -PIBY2) Stereo[icath] += 2.0*Pi;

               TrueWireHit[icath] += (float)sqrt(pow(cos((double)Stereo[icath])*
                   (double)lorerr,2) + pow(cos((double)Stereo[icath])*(double)phierr,2));

// Adjust hit position for stereo cathode strips which are NOT 
// perpendicular to the anode wires:

               if (Stereo[icath]){
                   if (istation==Station1 || istation==Station3){
                     stererr = (wireIP)*(float)sin((double)Stereo[icath])/
                             ReadoutSpacing;
                   }
                   else if (istation==Station2){
                     stererr = -(wireIP)*(float)sin((double)Stereo[icath])/
                             ReadoutSpacing;
                   }
               }
               else stererr = 0.0;

               TrueWireHit[icath] += stererr; 

//  PERFORM CATHODE STRIP CHARGE DISTRIBUTION ON BOTH CATHODE PLANES
//  BOTH CATHODE PLANES WILL BE TREATED IDENTICALLY.

// The strip width that needs to be passed to mMutCathStrip is the real
//  strip spacing, not the readout strip spacing, so padwidth = readout
//  strip spacing divided by 2 since every other strip read out:

               padwidth = ReadoutSpacing/2.0;
               localx = TrueWireHit[icath];

               truewireno = (long)(ptrStrip[icath]->getStrip()*2);
               iarm = (long)ptrStrip[icath]->getArm();
               istation = (long)ptrStrip[icath]->getStation();
               isegment = (long)ptrStrip[icath]->getOctant();
               ihalf = (long)ptrStrip[icath]->getHalfOctant();
               igap = (long)ptrStrip[icath]->getGap();
               iplane = igap*2 + icath;
               WireNoMin = 0;
               WireNoMax = ptrStrip[icath]->getNumberOfStrips();

               mMutCathStrip(localx, padwidth, truewireno, dedx,
                 iplane, icath, igap, ihalf, isegment, istation, iarm, 
                 ihit+1, nmmul, truehit, TrueWireHit[icath] ,ntrks, WireNoMin, WireNoMax);

             }  //  If a strip on this plane was hit
           }  //  Loop over cathode planes
         }  //  If this hit position converts to a good strip position
       }  //  This is a hit with good station number (still needed?)	
     }  //  This is a hit in a muon tracking chamber
   }  //  Loop over all muon hits	

//  Write data to dMutCathodesFEM table and dMutAnodesFEM table	

   dMutAnodesFEM_h.nok = 0;
   dMutAnodesFEMMC_h.nok = 0;
   AnFemMin = 0;
   CathFemMin = 0;
   NumberOfGaps = NorthArm()->f_pMutStations[NumberOfStations-1]
                    ->f_pMutOctants[NumberOfOctants-1]
                    ->f_pMutHalfOctants[1]->getNumberOfGaps();
   NumberOfStrips = NorthArm()->f_pMutStations[NumberOfStations-1]
                      ->f_pMutOctants[NumberOfOctants-1]
                      ->f_pMutHalfOctants[1]
                      ->f_pMutGaps[NumberOfGaps-1]
                      ->f_pMutPlanes[NumberOfPlanes-1]->f_pMutStrips.size();
   CathFemTot = NorthArm()->f_pMutStations[NumberOfStations-1]
                  ->f_pMutOctants[NumberOfOctants-1]
                  ->f_pMutHalfOctants[1]
                  ->f_pMutGaps[NumberOfGaps - 1]
                  ->f_pMutPlanes[NumberOfPlanes-1]
                  ->f_pMutStrips[NumberOfStrips-1]->getPacket_ID() - 
                  PACKETIDSTART;

   dMutCathodesFEM_h.nok = CathFemTot;
   dMutCathodesFEMMC_h.nok = CathFemTot;

   data = (short)calib->getPedestal(0, 0, 0, 0, 0, 0, 0);
   goodfill = mMutZeroCathodeFEMTables(CathFemMin, CathFemTot, data, 
     &dMutCathodesFEM_h, dMutCathodesFEM, &dMutCathodesFEMMC_h,
     dMutCathodesFEMMC);

   for ( iarm = 0; iarm < NumberOfArms; iarm++ )
   {
     switch (iarm){
       case (South):
         ArmPointer = SouthArm();
         break;
       case (North):
         ArmPointer = NorthArm();
         break;
       default:
         cout << "ERROR: invalid arm number used in mMutXYtoFEE" << endl;
         ArmPointer = 0;
         break;

     } 
     for ( istation =0; istation < NumberOfStations; istation++ )
     {
       for ( isegment =0; isegment < NumberOfOctants; isegment++ )
       {
         for ( ihalf=0; ihalf< NumberOfHalfOctants; ihalf++ )
         {
           NumberOfGaps = ArmPointer->f_pMutStations[istation]->f_pMutOctants[isegment]->
                            f_pMutHalfOctants[ihalf]->getNumberOfGaps();

           for ( igap =0; igap < NumberOfGaps; igap++ )
           {
             for ( icath=0; icath< NumberOfCathodePlanes; icath++ )
             {
               if (icath==0) CathodePlane = Cathode1;
               else CathodePlane = Cathode2;

               iplane = igap*NumberOfCathodePlanes + icath;
               j = 0;

               while (ArmPointer->f_pMutStations[istation]->
                      f_pMutOctants[isegment]->f_pMutHalfOctants[ihalf]->
                      f_pMutGaps[igap]->f_pMutPlanes[CathodePlane]->
                      f_pMutStrips[j] && 
                      j < ArmPointer->f_pMutStations[istation]->
                      f_pMutOctants[isegment]->f_pMutHalfOctants[ihalf]->
                      f_pMutGaps[igap]->f_pMutPlanes[CathodePlane]->
                      getNumElements())
  	       {

//  See if channel is "dead":

                 ntracks = ntrks[j][iarm][istation][isegment][ihalf][igap][icath];
                 if (ntracks  >  10)
                 {
                   cout << "more than 10 tracks contributing to cathode strip";
                    ntracks = 10;
                 }

                 pedestal = calib->getPedestal(iarm, istation, isegment, 
                   ihalf, igap, icath, j);
                 gain = calib->getGain(iarm, istation, isegment, 
                   ihalf, igap, icath, j);
                 rms = calib->getRms(iarm, istation, isegment, 
                   ihalf, igap, icath, j);

                 t0 = 0.0;

                 vpeak = mtdigi.readout_strip[j][iarm][istation][isegment][ihalf][igap][icath];
                 if (vpeak<0) vpeak = 0;

                 goodadc = mMutCalcADCSamples(&trise, &tfall, &vpeak, &t0, 
                   &nsamples, tsamples, chargesamp, &fitv0);

//  Add noise to all of the readout strips, cathodes and anodes.
//  The noise is in units of energy loss	
                 if (NOISE)
                 {

//  Just put noise on hit strips for now.  Need to add noise to readout_strip for
//  anode signals (use this to determine latch); need to add noise to chargesamp for
//  cathode strips.

                   if(mtdigi.readout_strip[j][iarm][istation][isegment][ihalf][igap][icath] !=  0)
                   {
                     ranx = RanGen->Gaus(0.0,1.0); 
                     mtdigi.readout_strip[j][iarm][istation][isegment][ihalf][igap][icath] +=
                       rms * ranx;

                     for (isamp=0; isamp<nsamples; isamp++){
                       ranx = RanGen->Gaus(0.0,1.0); 
                       chargesamp[isamp] += rms * ranx;
                     }
                   }
                 }            /* If adding noise to strips	*/
  
                 for (isamp=0; isamp<nsamples; isamp++){
                   adcsamp[isamp] = int(gain*chargesamp[isamp] + pedestal + 0.5);
                   if (adcsamp[isamp] < adcmin) adcsamp[isamp]=adcmin;
                   if (adcsamp[isamp] > adcmax) adcsamp[isamp]=adcmax;
                 }

//  See what FEM Address this strip corresponds to:
                 DCMChannel = ArmPointer->f_pMutStations[istation]
                                        ->f_pMutOctants[isegment]
                                        ->f_pMutHalfOctants[ihalf]
                                        ->f_pMutGaps[igap]
                                        ->f_pMutPlanes[CathodePlane]
                                        ->f_pMutStrips[j]->getDCMChannel();
                 PacketID   = ArmPointer->f_pMutStations[istation]
                                        ->f_pMutOctants[isegment]
                                        ->f_pMutHalfOctants[ihalf]
                                        ->f_pMutGaps[igap]
                                        ->f_pMutPlanes[CathodePlane]
                                        ->f_pMutStrips[j]->getPacket_ID();

                 datastart = DCMChannel*nsamples;

                 goodfill = mMutFillCathodeFEMTables(iarm, istation, isegment, ihalf,
                   igap, icath, j, PacketID, DCMChannel,
                   datastart, ntracks, truehit, nmmul, adcsamp, nsamples, 
                   tsamples, tspace,
                   &dMutCathodesFEM_h, dMutCathodesFEM, &dMutCathodesFEMMC_h, 
                   dMutCathodesFEMMC);

               j++;
               }     /* Loop over strips */
             }       /* Loop over planes within a gap */
           }       /* Loop over gaps */
         }       /* Loop over half-octants*/
       }       /* Loop over octants */
     }       /* Loop over stations */
   }       /* Loop over arms */

// Set the row count of tables that have been filled:

   tdMutCathodesFEMWrapper->SetRowCount(dMutCathodesFEM_h.nok);
   tdMutCathodesFEMMCWrapper->SetRowCount(dMutCathodesFEMMC_h.nok);
   tdMutAnodesFEMWrapper->SetRowCount(dMutAnodesFEM_h.nok);
   tdMutAnodesFEMMCWrapper->SetRowCount(dMutAnodesFEMMC_h.nok);
   tdMutCathodeChargeWrapper->SetRowCount(dMutCathodeCharge_h.nok);

// Delete random number generator class:

   delete RanGen;

   return ( STAFCV_OK );
}