// $Id$
//
//    File: DBCALShower_factory_KLOE.cc
// Created: Tue Jul  3 18:25:12 EDT 2007
// Creator: Matthew Shepherd
//

#include <cassert>
#include <math.h>
#include <map>

#include "BCAL/DBCALHit.h"
#include "BCAL/DBCALGeometry.h"
#include "BCAL/DBCALShower_factory_KLOE.h"

#include "units.h"

using namespace std;

//------------------
// DBCALShower_factory_KLOE
//------------------
DBCALShower_factory_KLOE::DBCALShower_factory_KLOE()
{
    // this should be lower than cut in mcsmear
	ethr_cell=0.0001;     // MIN ENERGY THRESD OF cell in GeV

	CLUST_THRESH = 0.03;  // MIN ENERGY THRESD OF CLUSTER IN GEV (make this match the value used in the other algorithm)
    
	elyr = 1;
	xlyr = 2; 
	ylyr = 3; 
	zlyr = 4; 
	tlyr = 5;
    
  // these four parameters are used in merging clusters
	MERGE_THRESH_DIST   = 40.0;  // CENTROID DISTANCE THRESHOLD
	MERGE_THRESH_TIME   =  2.5;  // CENTROID TIME THRESHOLD
	MERGE_THRESH_ZDIST  = 30.0;  // FIBER DISTANCE THRESHOLD
	MERGE_THRESH_XYDIST = 40.0;  // CENTROID TRANSVERSE DISTANCE THRESHOLD

  // this parameter is used to break clusters based on rms time
  BREAK_THRESH_TRMS= 5.0;   // T RMS THRESHOLD
    
  gPARMS->SetDefaultParameter( "BCALRECON:CLUST_THRESH", CLUST_THRESH );
  gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_DIST", MERGE_THRESH_DIST );
  gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_TIME", MERGE_THRESH_TIME );
  gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_ZDIST", MERGE_THRESH_ZDIST );
  gPARMS->SetDefaultParameter( "BCALRECON:MERGE_THRESH_XYDIST", MERGE_THRESH_XYDIST );
  gPARMS->SetDefaultParameter( "BCALRECON:BREAK_THRESH_TRMS", BREAK_THRESH_TRMS );
  
  if( !DBCALGeometry::summingOn() ){

    // these are energy calibration parameters -- no summing of cells

    m_scaleZ_p0 =  0.9597;
    m_scaleZ_p1 =  0.000454875;
    m_scaleZ_p2 =  -2.29912e-06;
    m_scaleZ_p3 =  1.49757e-09;
  
    m_nonlinZ_p0 =  -0.00154122;
    m_nonlinZ_p1 =  6.73594e-05;
    m_nonlinZ_p2 =  0;
    m_nonlinZ_p3 =  0;
  
  }
  else{
    
    // these are energy calibration parameters -- 1.2.3.4 summing
  
    //last updated for svn revision 9233
    m_scaleZ_p0 = 0.99798;
    m_scaleZ_p1 = 0.000361096;
    m_scaleZ_p2 = -2.17338e-06;
    m_scaleZ_p3 = 1.32201e-09;
  
    m_nonlinZ_p0 = -0.0201272;
    m_nonlinZ_p1 = 0.000103649;
    m_nonlinZ_p2 = 0;
    m_nonlinZ_p3 = 0;
  }
  
}

//------------------
// brun
//------------------
jerror_t DBCALShower_factory_KLOE::brun(JEventLoop *loop, int runnumber)
{
    
    vector<const DBCALGeometry*> bcalGeomVect;
    loop->Get( bcalGeomVect );
    const DBCALGeometry& bcalGeom = *(bcalGeomVect[0]);
    
    //////////////////////////////////////////////////////////////////
    // Calculate Cell Position
    //////////////////////////////////////////////////////////////////
    
    ATTEN_LENGTH = bcalGeom.ATTEN_LENGTH;
    C_EFFECTIVE = bcalGeom.C_EFFECTIVE;
    
    fiberLength = bcalGeom.BCALFIBERLENGTH; // fiber length in cm
    zOffset = bcalGeom.GLOBAL_CENTER;

    //the following uses some sad notation in which modmin=0 and modmax=48, when in fact there are 48 modules labelled either 0-47 or 1-48 depending on one's whim, although if we are using the methods from DBCALGeometry (e.g. cellId()), we must start counting from 1 and if we are accessing arrays we must of course start from 0
    int   modmin = 0;
    int   modmax = bcalGeom.NBCALMODS;
    int   rowmin1=0;
    int   rowmax1= bcalGeom.NBCALLAYSIN;
    int   rowmin2= rowmax1;
    int   rowmax2= bcalGeom.NBCALLAYSOUT+rowmin2; 
    int   colmin1=0;
    int   colmax1=bcalGeom.NBCALSECSIN;
    int   colmin2=0;
    int   colmax2=bcalGeom.NBCALSECSOUT;
    
    float r_inner= bcalGeom.BCALINNERRAD;
    
    for (int i = (rowmin1+1); i < (rowmax1+1); i++){
        //this loop starts from 1, so we can use i in cellId with no adjustment
        int cellId = bcalGeom.cellId(1,i,1); //this gives us a cellId that we can use to get the radius. the module and sector numbers are irrelevant
        //rt is radius of center of layer - BCAL inner radius
        rt[i]=bcalGeom.r(cellId)-r_inner;
    }
    
    for (int i = (rowmin2+1); i < (rowmax2+1); i++){
        int cellId = bcalGeom.cellId(1,i,1);
        rt[i]=bcalGeom.r(cellId)-r_inner;
    }
    
    //these are r and phi positions of readout cells
    float r[modulemax_bcal][layermax_bcal][colmax_bcal];
    float phi[modulemax_bcal][layermax_bcal][colmax_bcal];
    
    // Now start to extract cell position information from Geometry class
    for (int k = modmin; k < modmax; k++){
        for (int i = rowmin1; i < rowmax1; i++){
            for (int j = colmin1; j < colmax1; j++){
                //in this case the loops start at 0, so we have to add 1 to the indices when calling cellId(). hooray!
                //use DBCALGeometry to get r/phi position of each cell
                int cellId = bcalGeom.cellId(k+1,i+1,j+1);
                r[k][i][j]=bcalGeom.r(cellId);
                phi[k][i][j]=bcalGeom.phi(cellId);
                //set x and y positions
                xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]);
                yy[k][i][j]=r[k][i][j]*sin(phi[k][i][j]);          
            }
        }
        
        for (int i = rowmin2; i < rowmax2; i++){
            for (int j = colmin2; j < colmax2; j++){
                int cellId = bcalGeom.cellId(k+1,i+1,j+1);
                r[k][i][j]=bcalGeom.r(cellId);
                phi[k][i][j]=bcalGeom.phi(cellId);

                xx[k][i][j]=r[k][i][j]*cos(phi[k][i][j]);
                yy[k][i][j]=r[k][i][j]*sin(phi[k][i][j]);          
            }
        }         
    }
    
    ////////////////////////////////////////////////////////////////////////////
    // Now the cell information are already contained in xx and yy arrays.
    // xx and yy arrays are private members of this class
    ////////////////////////////////////////////////////////////////////////////
    
    // these are timing offsets -- a global side offset and individual cell offsets
    ta0=0.0;
    tb0=0.0;    
    
    for (int i = 0; i < modulemax_bcal; i++){
        for (int j = 0; j < layermax_bcal; j++){
            for (int k = 0; k < colmax_bcal; k++){                       
                
                ta_offset[i][j][k]=0.0;
                tb_offset[i][j][k]=0.0;
            }
        }
    }
	 
	 return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DBCALShower_factory_KLOE::evnt(JEventLoop *loop, int eventnumber)
{
    // Needed for associated objects later
    vector<const DBCALHit*> bcalhits;
    loop->Get(bcalhits);

    // Call core KLOE reconstruction routines
    CellRecon(loop);
    CeleToArray();   
    PreCluster(loop); 
    ClusNorm();
    ClusAnalysis();
    Trakfit();
    
    // Loop over reconstructed clusters and make DBCALShower objects out of them
    vector<DBCALShower*> clusters;
    
    int id = 0;
    for (int i = 1; i < (clstot+1); i++){
        
        int  j=clspoi[i];
  
        if( e_cls[j] < CLUST_THRESH ) continue;
        
        // Time to cook a final shower
        DBCALShower *shower = new DBCALShower;
  
        shower->id                  = id++;
        shower->E_raw               = e_cls[j];
        shower->x                   = x_cls[j];
        shower->y                   = y_cls[j];
        shower->z                   = z_cls[j] + zOffset;   
        shower->t                   = t_cls[j];
        shower->N_cell              = ncltot[j];
      
        shower->xErr                = eapx[1][j];
        shower->yErr                = eapx[2][j];
        shower->zErr                = eapx[3][j];

        shower->tErr                = 0.5 * sqrt( trms_a[j] * trms_a[j] +
                                                  trms_b[j] * trms_b[j] );
      
        // calibrate energy:
        // Energy calibration has a z dependence -- the
        // calibration comes from fitting E_rec / E_gen to scale * E_gen^nonlin
        // for slices of z.  These fit parameters (scale and nonlin) are then plotted 
        // as a function of z and fit.
  
        // center of target should be a geometry lookup
        float zTarget = 65*k_cm;
        float r = sqrt( shower->x * shower->x + shower->y * shower->y );
      
        float zEntry = ( shower->z - zTarget ) * ( DBCALGeometry::BCALINNERRAD / r );
      
        float scale = m_scaleZ_p0  + m_scaleZ_p1*zEntry + 
                      m_scaleZ_p2*(zEntry*zEntry) + m_scaleZ_p3*(zEntry*zEntry*zEntry);
        float nonlin = m_nonlinZ_p0  + m_nonlinZ_p1*zEntry + 
                       m_nonlinZ_p2*(zEntry*zEntry) + m_nonlinZ_p3*(zEntry*zEntry*zEntry);

        shower->E = pow( (shower->E_raw ) / scale, 1 / ( 1 + nonlin ) );

        //copy xyz errors into covariance matrix
        shower->xyzCovariance.ResizeTo(3,3);
        shower->xyzCovariance[0][0] = shower->xErr*shower->xErr;
        shower->xyzCovariance[1][1] = shower->yErr*shower->yErr;
        shower->xyzCovariance[2][2] = shower->zErr*shower->zErr;

	// Trace back to the DBCALHit objects used in this shower and
	// add them as associated objects.
	vector<const DBCALHit*> hitsInShower;
	FindHitsInShower(j, bcalhits, hitsInShower);
	for(unsigned int j=0; j<hitsInShower.size(); j++){
		shower->AddAssociatedObject(hitsInShower[j]);
	}

      _data.push_back(shower);  
    }
    
	return NOERROR;
}


//------------------
// FindHitsInShower()
//------------------
void DBCALShower_factory_KLOE::FindHitsInShower(int indx, vector<const DBCALHit*> &bcalhits, vector<const DBCALHit*> &hitsInShower)
{
	/// This is called after the clusters have been completely formed. Our
	/// job is simply to find the DBCALHit objects used to form a given cluster.
	/// This is so the DBCALHit objects can be added to the DBCALShower object
	/// as AssociatedObjects.

	// The variable indx indexes the next[] array as the starting cell for the
	// cluster. It also indexes the narr[][] array which holds the module, layer,
	// sector(column) values for the hits.
	//
	// Here, we need to loop over next[] elements starting at next[indx] until
	// we find the one pointing to element "indx" (i.e. the start of the list of
	// cells in the cluster.) For each of these, we must find the LAST 
	// member in bcalhits to have the same module, layer, number indicating 
	// that is a DBCALHit to be added.
	
	int start_indx = indx;
	do{
		int module = narr[1][indx];
		int layer  = narr[2][indx];
		int sector = narr[3][indx];
		
		// Loop over BCAL hits, trying to find this one
		const DBCALHit *uphit=NULL;
		const DBCALHit *downhit=NULL;
		for(unsigned int i=0; i<bcalhits.size(); i++){
			if(bcalhits[i]->module !=module)continue;
			if(bcalhits[i]->layer  !=layer)continue;
			if(bcalhits[i]->sector !=sector)continue;
			
			if(bcalhits[i]->end == DBCALGeometry::kUpstream) uphit = bcalhits[i];
			if(bcalhits[i]->end == DBCALGeometry::kDownstream) downhit = bcalhits[i];
		}
		
		if(uphit)hitsInShower.push_back(uphit);
		if(downhit)hitsInShower.push_back(downhit);

		indx = next[indx];
	}while(indx != start_indx);

}


//------------------
// CellRecon()
//------------------
void DBCALShower_factory_KLOE::CellRecon(JEventLoop *loop)
{
    //======================================================================
    // This code is used to reconstruct cell information cell by cell.
    // This code is adapted from kloe code by Chuncheng Xu. June 29,2005
    
    //**********************************************************************
    // The main purpose of this function is extracting information
    // from DBCALHit
    // objects to form the arrays: ecel_a,tcel_a,ecel_b,tcel_b 
    // and xcel,ycel,zcel,tcel,ecel,tcell_anor,tcell_bnor;
    // Among these arrays,
    // ecel_a,tcel_a,ecel_b,tcel_b, xcel,ycel,zcel,tcel,ecel,tcell_anor and
    // tcell_bnor are the input of function CeleToArray()
    // These arrys are 3-D arrays and are rather logically are indexed by
    // [module][layer][column]
    //********************************************************************** 
    
    /////////////////////////////////////////////////////////////////////
    // Now start to take take out DBCALHit to form our private 
    // data member ecel_a,tcel_a,ecel_b,tcel_b
    /////////////////////////////////////////////////////////////////////
    
    memset( ecel_a, 0, modulemax_bcal * layermax_bcal *
            colmax_bcal * sizeof( float ) );
    memset( tcel_a, 0, modulemax_bcal * layermax_bcal *
            colmax_bcal * sizeof( float ) );
    memset( ecel_b, 0, modulemax_bcal * layermax_bcal *
            colmax_bcal * sizeof( float ) );
    memset( tcel_b, 0, modulemax_bcal * layermax_bcal *
            colmax_bcal * sizeof( float ) );
    
    // extract the BCAL hits
    
    vector<const DBCALHit*> hits;
    loop->Get(hits);
    if(hits.size() <= 0) return;
    
        
    for (unsigned int i = 0; i < hits.size(); i++) {
        
        const  DBCALHit *hit = hits[i];     
        int module = hit->module;
        int layer = hit->layer;
        int sector = hit->sector;
        int end = hit->end;
        float E = hit->E;
        float t = hit->t;
        
        // cache hits indexed by module/layer/sector
        // this allows end to end matching in next step
        
        if(end==0) {                         // upstream
            ecel_a[module-1][layer-1][sector-1] =  E;
            tcel_a[module-1][layer-1][sector-1] =  t;
        }  
        else if(end==1) {                //   downstream
            ecel_b[module-1][layer-1][sector-1] =  E;
            tcel_b[module-1][layer-1][sector-1] =  t;	      
        }
        else
        {
            cout<<"no such end, it is neither side A Nor B \n"; 
        }
    }

    //////////////////////////////////////////////////////////////////////// 
    // Now data member ecel_a,tcel_a,ecel_b,tcel_b are readily filled.
    ////////////////////////////////////////////////////////////////////////
    
    //////////////////////////////////////////////////////////////////////// 
    // Now start to reconstruct cell information from data member ecel_a,tcel_a
    // ecel_b,tcel_b. The reconstructed cell information are stored in 
    // data member arrays xcel,ycel,zcel,tcel,ecel,tcell_anor,tcell_bnor.
    ////////////////////////////////////////////////////////////////////////
    
    //
    // ************ Loop over all CELE elements ********************
    //
    // K module,I layer, J sector
    
    for (int k = 0; k < modulemax_bcal; k++){
        for (int i = 0; i < layermax_bcal; i++){
            for (int j = 0; j <  colmax_bcal; j++){    
                
                // get times
                float  ta = tcel_a[k][i][j];
                float  tb = tcel_b[k][i][j];  
                
                if( ( ta == 0 ) || ( tb == 0 ) || 
                    ( fabs( ta - tb ) > 80. ) ) {  
                    
                    xcel[k][i][j] = 0.0;
                    ycel[k][i][j] = 0.0;
                    zcel[k][i][j] = 0.0;
                    tcel[k][i][j] = 0.0;
                    ecel[k][i][j] = 0.0;

                    continue;
                }
                
                float ea = ecel_a[k][i][j];
                float eb = ecel_b[k][i][j];
                
                // algorithm requires both ends be hit
                // improve energy resolution with single hits?
                if( ( ea < ethr_cell ) ||
                    ( eb < ethr_cell ) ){

                    xcel[k][i][j] = 0.0;
                    ycel[k][i][j] = 0.0;
                    zcel[k][i][j] = 0.0;
                    tcel[k][i][j] = 0.0;
                    ecel[k][i][j] = 0.0;
                    continue;
                }
                
                
                float  ta_cal = ta - ta0 - ta_offset[k][i][j];
                float  tb_cal = tb - tb0 - tb_offset[k][i][j];
                
                float x    = xx[k][i][j];
                float y    = yy[k][i][j];
                    
                float  z1 = -0.5 * fiberLength;    // cm
                float  z2 =  0.5 * fiberLength;    // cm      
                float  z0 =  0.5 * ( z1 + z2 );
                    
                float z = 0.5 * C_EFFECTIVE * ( ta_cal - tb_cal ) + z0;
                    
                if( z > ( 0.5 * fiberLength ) ){
                        
                    z = 0.5*fiberLength;
                }
                else if( z < ( -0.5 * fiberLength ) ){
                        
                    z = -0.5*fiberLength;
                }
                    
                float t = 0.5 * ( ta_cal + tb_cal - fiberLength / C_EFFECTIVE );
                
                float dpma = min( ( ta_cal - t ) * C_EFFECTIVE, fiberLength );
                float dpmb = min( ( tb_cal - t ) * C_EFFECTIVE, fiberLength );
                float atta = exp( -dpma / ATTEN_LENGTH );
                float attb = exp( -dpmb / ATTEN_LENGTH );

                // could compute energy weighted average instead of straight average
                float e = ( ( ea / atta ) + ( eb / attb ) ) / 2;
                
                xcel[k][i][j] = x;
                ycel[k][i][j] = y;
                zcel[k][i][j] = z;
                tcel[k][i][j] = t;
                ecel[k][i][j] = e;
                tcell_anor[k][i][j] = ta_cal;
                tcell_bnor[k][i][j] = tb_cal;
            }
        }
    }    
}


//------------------
// CeleToArray()
//------------------
void DBCALShower_factory_KLOE::CeleToArray(void)
{
    //  THis code is adpapted from kloe code by Chuncheng Xu on June29,2005
    //  The following part is taken from kaloe's clurec_lib.f's subroutine
    //  cele_to_arr.
    //
    //  It's original purpose is to extract the data information from array RW
    //  and IW into array NARR(j,i), CELDATA(j,i), nclus(i),next(i),
    //  and e_cel(i), x_cel(i),y_cel(i), z_cel(i), t_cel(i),ta_cel(i)
    //  and tb_cel(i)
    //  In the above explanationn, i is the index for cell, and different j
    //  is for different component for such cell.
    //  for example, NARR(1,i), NARR(2,i),NARR(3,i) are for module number,
    //  layer number and sector number in halld bcal simulation 
    
    
    //  Now  we assume that the information is stored in arrays ECEL_A(k,i,j)
    //  , ECEL_B(k,i,j), TCEL_A(k,i,j),  TCEL_B(k,i,j), XCEL(k,I,J),YCEL(k,I,J), 
    // ZCEL(K,I,J), TCEL(K,I,J), ECEL(K,I,J),TCELL_ANOR(K,I,J),TCELL_BNOR(K,I,J),
    //  which are passed into here by event.inc through common block.  
    
    //   these values e_cel(i), x_cel(i),y_cel(i), z_cel(i), t_cel(i),ta_cel(i)
    //  and tb_cel(i)) are extremly useful in later's clusterization
    //  and they are passed by common block to precluster subroutine too. (through
    //  clurec_cal.inc)
    // 
    //  we hope this is a good "bridge" from raw data to precluster. 
    //  This way we can keep good use of their strategy of clusterization
    //  as much as possible, although we know that Fortran has it's bad fame 
    //  of not so easy to be reused.


    // Essentially this function takes the cell information (e.g. E,x,y,z,t)
    // from the 3D arrays filled in CellRecon() and puts it into 1D arrays.
    // These 1D arrays are indexed by an identifier. The module/layer/sector
    // associated with this identifier can be found using the narr[][] array,
    // as described above.
    
    celtot=0;
    
    for (int k = 0; k < modulemax_bcal; k++){
        for (int i = 0; i < layermax_bcal; i++){
            for (int j = 0; j < colmax_bcal; j++){     
                
                float   ea  = ecel_a[k][i][j];
                float   eb  = ecel_b[k][i][j];
                float   ta  = tcel_a[k][i][j];
                float   tb  = tcel_b[k][i][j];
                
                if( (min(ea,eb)>ethr_cell) & (fabs(ta-tb)<35.) & (ta!=0.) & (tb!=0.)) { 
		  celtot=celtot+1;             
		} else {
		  continue;
                }
                
                
                if(celtot>cellmax_bcal) {
                    break;
                }
                
                narr[1][celtot]=k+1;    // these numbers will
                narr[2][celtot]=i+1;    //  be used by preclusters
                narr[3][celtot]=j+1;    //  which will start from index of 1
                                        // rather than from 0.
                                           
                // why 0.145? -- these variables are used as weights
                celdata[1][celtot]=ea/0.145;
                celdata[2][celtot]=eb/0.145;
                
                nclus[celtot] = celtot;
                next[celtot]  = celtot;
                
                e_cel[celtot] = ecel[k][i][j];
                x_cel[celtot] = xcel[k][i][j];
                y_cel[celtot] = ycel[k][i][j];
                z_cel[celtot] = zcel[k][i][j];
                t_cel[celtot] = tcel[k][i][j];
                
                ta_cel[celtot]=tcell_anor[k][i][j];
                tb_cel[celtot]=tcell_bnor[k][i][j];
            }
        }
    }    
}



//------------------
// PreCluster()
//------------------        
void DBCALShower_factory_KLOE::PreCluster(JEventLoop *loop)
{
  //what this function does: for each cell with a hit it finds the maximum
  //energy neighbor and Connect()'s the two. Two cells are neighbors if they
  //are within one column of each other and within one row. The situation is
  //slightly more complicated for two cells on opposite sides of the boundary
  //between inner cells and outer and is described in more detail below, but
  //essentially works out the same way. The purpose of Connect() is described
  //in that function itself.

  int k=1;     // NUMBER OF NEARBY ROWS &/OR TO LOOK FOR MAX E CELL
    
  // extract the BCAL Geometry
  vector<const DBCALGeometry*> bcalGeomVect;
  loop->Get( bcalGeomVect );
  const DBCALGeometry& bcalGeom = *(bcalGeomVect[0]);
    
  // calculate cell position
    
  int   modmin = 0;
  int   modmax = bcalGeom.NBCALMODS;


  //these values make sense as actually minima/maxima if it is implied that rowmin1=1,colmin1=1,colmin2=1
  int   rowmax1= bcalGeom.NBCALLAYSIN;
  int   rowmin2= rowmax1+1;
  //int   rowmax2= bcalGeom.NBCALLAYSOUT+rowmin2-1;
  int   colmax1=bcalGeom.NBCALSECSIN;
  int   colmax2=bcalGeom.NBCALSECSOUT;

  float r_middle= bcalGeom.BCALMIDRAD;

  //radial size of the outermost inner layer
  float thick_inner=bcalGeom.rSize(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN,1));
  //radial size of the innermost outer layer
  float thick_outer=bcalGeom.rSize(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN+1,1));

  // this is the radial distance between the center of the innermost outer layer and the outermost inner layer
  float dis_in_out=bcalGeom.r(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN+1,1))-bcalGeom.r(bcalGeom.cellId(1,bcalGeom.NBCALLAYSIN,1));
    
  float degree_permodule=360.0/(modmax-modmin);
  float half_degree_permodule=degree_permodule/2.0;

  //roughly the width of a single cell in the outermost inner layer
  float width_1=2.0*(r_middle-thick_inner/2.0)*
    sin(half_degree_permodule*3.141593/180)/colmax1;
  //roughly the width of a single cell in the innermost outer layer
  float width_2=2.0*(r_middle+thick_outer/2.0)*
    sin(half_degree_permodule*3.141593/180)/colmax2;

  //disthres is roughly the azimuthal distance between the center of an outer cell and the center of the most distant inner cell bordering an adjacent outer cell (a picture would be nice wouldn't it)
  //this value is used for determining if two cells that straddle the boundary between inner and outer layers should be considered as neighboring
  float disthres=width_2*1.5-width_1*0.5+0.0001;
    
  for (int i = 1; i < (celtot+1); i++){
        
    int maxnn=0; //cell index of the maximum energy neighbor (if one is found)
    float emin=0.; //energy of maximum energy neighbor
        
        
    for (int j = 1; j < (celtot+1); j++){
      if ( (j!=i) & (nclus[j]!=nclus[i]) & (e_cel[j]>emin)) {
                
        int k1= narr[1][i];
        int k2= narr[1][j];
        int i1= narr[2][i];
        int i2= narr[2][j];


        int  modiff = k1-k2;
        int amodif = abs(modiff);
                
        //  the following if is to check module and row distance.         
        if ( (abs(i1-i2)<=k) & ((amodif<=1) || (amodif==47)) ) { 
          //   further check col distance 
          int   j1= narr[3][i];
          int   j2= narr[3][j];

          if(amodif==0) {   // same module       
            //   further check col distance if both are inner layers
            if ( (i1<=rowmax1) & (i2<=rowmax1) & (abs(j2-j1)<=k) ) {
              emin=e_cel[j];
              maxnn=j;
            }

            //   further check col distance if both are outer layers
 
            if ( (i1>=rowmin2) & (i2>=rowmin2) & (abs(j2-j1)<=k) ) {
              emin=e_cel[j];
              maxnn=j;
            }
          }

          if(amodif>0) {  // different module          
            if( (modiff==1) || (modiff==-47) ) {      
              if ( (i1<=rowmax1) & (i2<=rowmax1) ){ 
                if(abs((j1+colmax1)-j2)<=k){
                  emin=e_cel[j];
                  maxnn=j;
                }
              }
                        
              if ( (i1>=rowmin2) & (i2>=rowmin2) ) {
                if(abs((j1+colmax2)-j2)<=k){
                  emin=e_cel[j];
                  maxnn=j;
                }
              }              
            }

            if ( (modiff==-1) || (modiff==47) ) {      

              if ( (i1<=rowmax1) & (i2<=rowmax1) ){
                if(abs((j2+colmax1)-j1)<=k){
                  emin=e_cel[j];
                  maxnn=j;
                }
              } 
                        
              if ( (i1>=rowmin2) & (i2>=rowmin2) ){
                if(abs((j2+colmax2)-j1)<=k){
                  emin=e_cel[j];
                  maxnn=j;
                }
              }              
            }
          }
                    
          // further check col distance if one is inner layer, another is outer
          // so that the two may be between the boundary of two different size
          // of cells.
          if( ( (i1 == rowmax1) & (i2 == rowmin2) ) || 
              ( (i1 == rowmin2) & (i2 == rowmax1) ) ) {

            float delta_xx=xx[k1-1][i1-1][j1-1]-xx[k2-1][i2-1][j2-1];
            float delta_yy=yy[k1-1][i1-1][j1-1]-yy[k2-1][i2-1][j2-1];

            //distance between centers of two cells
            float dis = sqrt( delta_xx * delta_xx + delta_yy * delta_yy );
            //dis_in_out is the distance in radial direction, so we now isolate distance in direction perpendicular to radius
            dis = sqrt( dis*dis - dis_in_out * dis_in_out );
            //disthres is described above
            if( dis < disthres ){
              emin = e_cel[j];
              maxnn = j;
            }
          }                    
        }             
      }
    }        // finish second loop

    if(maxnn>0){

      Connect(maxnn,i);
    }
  }       // finish first loop
}



//------------------
// Connect(n,m);
//------------------   
void DBCALShower_factory_KLOE::Connect(int n,int m)
{
    //----------------------------------------------------------------------
    //   Purpose and Methods :CONNECTS CELL M TO THE NEAREST MAX NEIGHBOR N
    //   Created  31-JUL-1992   WON KIM
    //----------------------------------------------------------------------
    
    // This little piece of code wasn't so easy to decipher, but I *think* I
    // understand what it's doing. The idea is the following:
    //
    // (Prior to entering this routine)
    // The sparsified list of hits is copied into some 1-D arrays with
    // dimension cellmax_bcal+1. This includes the nclus[] and next[]
    // arrays. The nclus[] array keeps the cluster number which is initalized
    // to the sparsified hit number. Thus, every (double-ended) hit cell
    // is it's own cluster. A loop over pairs of hits is done above to find
    // nearest neighbors that should be merged into the same cluster.
    //
    // The next[] array contains an index to the "next" cell in the cluster
    // for the given hit cell. This is a circular list such that if one follows
    // the "next[next[next[...]]]" values, the last element points back to 
    // the first element of the cluster. As such, these are also initialized
    // to index themselves as all single hits are considered a cluster of
    // one element prior to calling this Connect() routine.
    //
    // When we are called, the value of "m" will be less than than value
    // of "n". The cluster number (kept in nclus[]) is therefore updated
    // for all members of nclus[m] to be the same as nclus[n]. Furthermore,
    // the hit cell "m" is appended to the cluster nclus[n]. This also means
    // that the cluster numbers become non-sequential since all elements of
    // nclus whose value is "m" will be changed such that their values are
    // "n".
    //
    // 5/10/2011 DL
    
    if(nclus[n]!=nclus[m]){
        int j=m;
        nclus[j]=nclus[n];
        while(next[j]!=m){
            j=next[j];
            nclus[j]=nclus[n];
	    }        
        next[j]=next[n];
        next[n]=m;
    }
}


//------------------
// ClusNorm()
//------------------   
void DBCALShower_factory_KLOE::ClusNorm(void)
{    
    // fast initialization of arrays:
    memset( e_cls,  0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( x_cls,  0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( y_cls,  0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( z_cls,  0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( t_cls,  0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( ea_cls, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( eb_cls, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( ta_cls, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( tb_cls, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( tsqr_a, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( tsqr_b, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( trms_a, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( trms_b, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( e2_a,   0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( e2_b,   0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( clspoi, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( ncltot, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    memset( ntopol, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );
    
    // Part of what is being done here is to further sparsify the
    // data into a list of clusters. This starts to fill arrays
    // with dimension clsmax_bcal+1. One important thing is that
    // the clspoi[] array is being filled with the cluster number
    // of what should be the index of the first cell hit in the
    // cluster. 
    //
    // 5/10/2011 DL

    clstot=0;
    
    for (int ix = 1; ix < (celtot+1); ix++){
        
        //----------------------------------------------------------------------
        // KEEP TALLY OF CLUSTERS
        //----------------------------------------------------------------------
        
        int n=nclus[ix];
        int j=0;
        
        for (int i = 1; i < (clstot+1); i++){
            if(n==clspoi[i]) j=i;
        }
        
        if(j==0) {
            clstot=clstot+1;
            clspoi[clstot]=n;
        }
        
        //----------------------------------------------------------------------
        // NORMALIZED QUANTITIES OF THE TOTAL CLUSTER
        //----------------------------------------------------------------------
        if(e_cel[ix]<0.000000001)continue; // just for protection
        
        x_cls[n]=(e_cls[n]*x_cls[n]+e_cel[ix]*x_cel[ix])
            /(e_cls[n]+e_cel[ix]);
        
        y_cls[n]=(e_cls[n]*y_cls[n]+e_cel[ix]*y_cel[ix])
            /(e_cls[n]+e_cel[ix]);
        
        z_cls[n]=(e_cls[n]*z_cls[n]+e_cel[ix]*z_cel[ix])
            /(e_cls[n]+e_cel[ix]);
        
        t_cls[n]=(e_cls[n]*t_cls[n]+e_cel[ix]*t_cel[ix])
            /(e_cls[n]+e_cel[ix]);
        
        e_cls[n]=e_cls[n]+e_cel[ix];
        
        //        write(*,*)' x-y-z-T-e done'
        //----------------------------------------------------------------------
        //       NORMALIZED QUANTITIES FOR EACH SIDE
        //----------------------------------------------------------------------
        
        ta_cls[n]=(ea_cls[n]*ta_cls[n]+celdata[1][ix]*ta_cel[ix])
            /(ea_cls[n]+celdata[1][ix]);        
        tsqr_a[n]=(ea_cls[n]*tsqr_a[n]+celdata[1][ix]*ta_cel[ix]*
                   ta_cel[ix])/(ea_cls[n]+celdata[1][ix]);
        ea_cls[n]=ea_cls[n]+celdata[1][ix];
        e2_a[n]=e2_a[n]+celdata[1][ix]*celdata[1][ix];
        tb_cls[n]=(eb_cls[n]*tb_cls[n]+celdata[2][ix]*tb_cel[ix])/
            (eb_cls[n]+celdata[2][ix]);
        tsqr_b[n]=(eb_cls[n]*tsqr_b[n]+celdata[2][ix]*tb_cel[ix]*
                   tb_cel[ix])/(eb_cls[n]+celdata[2][ix]);
        eb_cls[n]=eb_cls[n]+celdata[2][ix];
        e2_b[n]=e2_b[n]+celdata[2][ix]*celdata[2][ix];
        
        //----------------------------------------------------------------------
        //       TOTAL CELLS AND CELLS IN OTHER MODULES
        //----------------------------------------------------------------------
        
        ncltot[n]++;

        if( narr[1][n] != narr[1][ix] || narr[2][n] != narr[2][ix] )
            ntopol[n]++;
        
        //----------------------------------------------------------------------
        
    }
    
    for (int n = 1; n < (clstot+1); n++){
        
        int ix=clspoi[n];
        if( ncltot[ix] > 1) {
            
            float effnum = ea_cls[ix] * ea_cls[ix] / e2_a[ix];
            trms_a[ix] = effnum / ( effnum - 1 ) * 
                ( tsqr_a[ix] - ta_cls[ix] * ta_cls[ix] );

            effnum = eb_cls[ix] * eb_cls[ix] / e2_b[ix];
            trms_b[ix] = effnum / ( effnum - 1 ) * 
                ( tsqr_b[ix] - tb_cls[ix] * tb_cls[ix] );

            if( trms_a[ix] <= 0.0 ) trms_a[ix] = 0.;
            if( trms_b[ix] <= 0.0 ) trms_b[ix] = 0.;
            trms_a[ix] = sqrt( trms_a[ix] );
            trms_b[ix] = sqrt( trms_b[ix] );
        }
        else {
            trms_a[ix] = 0.;
            trms_b[ix] = 0.;
        }
    }
}

//------------------
// ClusAnalysis()
//------------------  
void DBCALShower_factory_KLOE::ClusAnalysis()
{
    // track when clusters change to cut down
    // on excess calles to expensive ClusNorm
    bool newClust = false;

    //----------------------------------------------------------------------
    //  Check for overlapping clusters
    //----------------------------------------------------------------------
    
    for (int i = 0; i < 2; i++){
        for (int j = 1; j < (clstot+1); j++){
            int ix=clspoi[j];
            if(e_cls[ix]>0.0){
                float  dist=sqrt(trms_a[ix]*trms_a[ix]+trms_b[ix]*trms_b[ix]);
                if(dist>BREAK_THRESH_TRMS) {
                    Clus_Break(ix);
                    newClust = true;
                }
            }
        }
        
        if( newClust ){
            
            ClusNorm();
            newClust = false;
        }
    }
    
    //----------------------------------------------------------------------       
    // merge clusters likely to be from the same shower
    //----------------------------------------------------------------------       

    int icls[3];
    for (int i = 1; i < clstot; i++){
        icls[1]=0;
        icls[2]=0;
        for (int j = (i+1); j < (clstot+1); j++){
            
            int ix=clspoi[i];
            int iy=clspoi[j];
            
            
            if ( (e_cls[ix]>0.0) & (e_cls[iy]>0.0) ) {
                
                float delta_x=x_cls[ix]-x_cls[iy];
                float delta_y=y_cls[ix]-y_cls[iy];          
                float delta_z=z_cls[ix]-z_cls[iy];  
                float dist=sqrt(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z);
                
                float  tdif=fabs(t_cls[ix]-t_cls[iy]);
                
                //	  float  distz=abs(z_cls[ix]-z_cls[iy]);
                //         cout<<"dist="<<dist<<" distz="<<distz<<" tdif="<<tdif<<"\n";    
                
                if ( (dist<MERGE_THRESH_DIST) & (tdif<MERGE_THRESH_TIME) ){
                    float zdif=fabs(z_cls[ix]-z_cls[iy]);
                    float distran=sqrt(delta_x*delta_x+delta_y*delta_y);
                    
                    if ( (zdif<MERGE_THRESH_ZDIST) & (distran<MERGE_THRESH_XYDIST) ){
                        if(e_cls[ix]>=e_cls[iy]) {
                            icls[1]=ix;
                            icls[2]=iy;
                        }
                        else {
                            icls[1]=iy;
                            icls[2]=ix;
                        }
                    }
                }
                
            }
            
            if(min(icls[1],icls[2])>0){
                
                Connect(icls[1],icls[2]);
                newClust = true;
            }
        }
    }
    
    if( newClust ){
        
        ClusNorm();
    }
}

//------------------
// Clus_Break(ix);
//------------------
void DBCALShower_factory_KLOE::Clus_Break(int nclust)
{
    int   nseed[5],selnum,selcel[cellmax_bcal+1];
    float tdif,tdif_a,tdif_b,tseed[5];  
    
    //----------------------------------------------------------------------
    for (int i =0; i < 5; i++){
        nseed[i]=0;
        tseed[i]=0;
    }
    //----------------------------------------------------------------------
    //   Divide cluster cells into four quadrant groups
    //----------------------------------------------------------------------
    int n=nclust;
    tdif_a=ta_cel[n]-ta_cls[nclust];
    tdif_b=tb_cel[n]-tb_cls[nclust];
    selnum=0;
    
    //----------------------------------------------------------------------
    if(tdif_a>0.0) {
        if(tdif_b>0){
            selnum=1;
        }
        else {
            selnum=2;
        }
    }
    
    else {
        if(tdif_b>0.0) {
            selnum=3;
        }
        else {
            selnum=4;
        }
    }
    
    //------------------------------------------------------------------------
    
    if(selnum>0) {
        float tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
        if(tdif>tseed[selnum]){
            nseed[selnum]=n;
            tseed[selnum]=tdif;
        }
        selcel[n]=selnum;
    }
    
    //-------------------------------------------------------------------------
    
    while(next[n]!=nclust) {
        n=next[n];
        tdif_a=ta_cel[n]-ta_cls[nclust];
        tdif_b=tb_cel[n]-tb_cls[nclust];
        selnum=0;
        
        //**************************************************************************
        
        if(tdif_a>0.0) {
            
            if(tdif_b>0.0) {
                selnum=1;
            }
            else {
                selnum=2;
            }
        }
        
        else {
            if(tdif_b>0.0) {
                selnum=3;
            } 
            else {
                selnum=4;
            }
        }
        
        //***************************************************************************
        
        //...........................................................................
        
        
        if(selnum>0){
            tdif=sqrt(tdif_a*tdif_a+tdif_b*tdif_b);
            
            if(tdif>tseed[selnum]){
                nseed[selnum]=n;
                tseed[selnum]=tdif;
            }
            
            selcel[n]=selnum;
        }
        
        //...........................................................................
        
    }             //this bracket is related to "while" above. 
    
    
    //----------------------------------------------------------------------
    // If successful, divide cluster chain into the new cluster chains
    //----------------------------------------------------------------------
    
    for (int i =1; i < 5; i++){
        
        if(nseed[i]>0) {
            
            nclus[nseed[i]]=nseed[i];
            next[nseed[i]]=nseed[i];

            for (int j =1; j < (celtot+1); j++){
	      if ( (nclus[j]==nclust) & (j!=nseed[i]) ){
                    if(selcel[j]==i) {
                        nclus[j]=j;
                        next[j]=j;
                        Connect(nseed[i],j);
                    }
                }
            }
        }
    }
}


//------------------
// Trakfit()
//------------------ 

// routine modified by MRS to do expensive filling of layer information
// contains code previously in ClusNorm which is called multiple times
// per event, but there is no need to call Trakfit multiple times per event
void DBCALShower_factory_KLOE::Trakfit( void )
{
 
    float emin=0.0001;

    memset( clslyr, 0, ( clsmax_bcal + 1 ) * 
            ( layermax_bcal + 1 ) * 6 * sizeof( float ) );
    
    memset( apx,   0, ( clsmax_bcal + 1 ) * 4 * sizeof( float ) );
    memset( eapx,  0, ( clsmax_bcal + 1 ) * 4 * sizeof( float ) );
    memset( ctrk,  0, ( clsmax_bcal + 1 ) * 4 * sizeof( float ) );
    memset( ectrk, 0, ( clsmax_bcal + 1 ) * 4 * sizeof( float ) );
    
    for (int ix = 1; ix < (celtot+1); ix++){
        
        int n = nclus[ix];
        
        //----------------------------------------------------------------------
        // NORMALIZED QUANTITIES OF THE TOTAL CLUSTER PER LAYER
        //----------------------------------------------------------------------
        
        int lyr = narr[2][ix];
        
        //	write(*,*)'X, E',E_CEL(ix),CLSLYR(XLYR,LYR,N)
        clslyr[xlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[xlyr][lyr][n]
                              +e_cel[ix]*x_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
        
        clslyr[ylyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[ylyr][lyr][n]
                              +e_cel[ix]*y_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
        //	write(*,*)' Y'
        
        clslyr[zlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[zlyr][lyr][n]
                              +e_cel[ix]*z_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
        
        //	write(*,*)' z'
        clslyr[tlyr][lyr][n]=(clslyr[elyr][lyr][n]*clslyr[tlyr][lyr][n]
                              +e_cel[ix]*t_cel[ix])/(clslyr[elyr][lyr][n]+e_cel[ix]);
        
        //	write(*,*)'E'
        clslyr[elyr][lyr][n]=clslyr[elyr][lyr][n]+e_cel[ix];
        
        //        write(*,*)' slopes done'
    }

    memset( nlrtot, 0, ( clsmax_bcal + 1 ) * sizeof( float ) );

    for (int n = 1; n < ( clstot + 1 ); n++){

        int ix=clspoi[n];
        
        for (int i = 1; i < (layermax_bcal+1); i++){
            
            if( clslyr[elyr][i][ix] > 0.0 ) nlrtot[ix]++;
        }
                    
        for (int i = 0; i < ( layermax_bcal + 1 ); i++){
        
            x[i]=0.0;
            y[i]=0.0;  
            z[i]=0.0;
            e[i]=0.0;
            sigx[i]=0.0;  // cm
            sigy[i]=0.0; // cm
            sigz[i]=0.0;
        }
        
        int nltot=0;
    
        for (int il = 1; il < (layermax_bcal+1); il++){

            if(clslyr[elyr][il][ix]>emin) {

                nltot=nltot+1;
                x[nltot]= clslyr[xlyr][il][ix];
                y[nltot]= clslyr[ylyr][il][ix];
                z[nltot]= clslyr[zlyr][il][ix];         
                e[nltot]= clslyr[elyr][il][ix]; 

                sigy[nltot] =  1.0/e[nltot];
                sigx[nltot] =  1.0/e[nltot];
                sigz[nltot] =  1.0/sqrt(e[nltot]);
            }
        }
    
        // The following error bar is the estimation of error bar
        // based on the experience of my fortran code running
        // If the structure of BCAL changes drastically, you have to make changes
        // accordingly.  Xu Chuncheng , Jan 9, 2006
        
        // **** this seems strange since it appears to overwrite what is above
    
        sigx[1]=0.5;
        sigx[2]=0.5;
        sigx[3]=0.5;
        sigx[4]=0.5;
        sigx[5]=0.5;
        sigx[6]=0.8;
        sigx[7]=0.9;
        sigx[8]=1.2;
        sigx[9]=1.3;
        
        sigy[1]=0.5;
        sigy[2]=0.5;
        sigy[3]=0.5;
        sigy[4]=0.5;
        sigy[5]=0.5;
        sigy[6]=0.8;
        sigy[7]=0.9;
        sigy[8]=1.2;
        sigy[9]=1.3;
        
        
        sigz[1]=0.5;
        sigz[2]=0.5;
        sigz[3]=0.5;
        sigz[4]=0.5;
        sigz[5]=0.5;
        sigz[6]=0.8;
        sigz[7]=0.9;
        sigz[8]=1.2;
        sigz[9]=1.3;
        
        if( nltot > 1 ){ 
            
            Fit_ls();
            
            for (int i = 1; i < 4; i++){
                
                ctrk[i][ix]=ctrk_ix[i];
                ectrk[i][ix]=ectrk_ix[i];
                apx[i][ix]=apx_ix[i];
                eapx[i][ix]=eapx_ix[i];
            }
        }    
        else{
            
            // we have 1 or less layers hit -- no fit
            
            apx[1][ix] = x[1];
            apx[2][ix] = y[1];
            apx[3][ix] = z[1]; 
            eapx[1][ix] = sigx[1];
            eapx[2][ix] = sigy[1];
            eapx[3][ix] = sigz[1];  
            ectrk[1][ix] = 0.0;
            ectrk[2][ix] = 0.0;
            ectrk[3][ix] = 0.0;
            ctrk[1][ix] = 999.0;
            ctrk[2][ix] = 999.0;
            ctrk[3][ix] = 999.0;
        }            
    }
}

//------------------
// Fit_ls()
//------------------  
void DBCALShower_factory_KLOE::Fit_ls()
{      
    float a,b,c;
    float d,e,f,chi2,q,norm;
    float siga,sigb,sigc,sigd,sige,sigf;
    float sigb2,sigd2,sigf2;
    
    //    fitting for X=a+bt
    Linefit(1,1,a,b,siga,sigb,chi2,q);
    //    fitting for Y=c+dt
    Linefit(2,1,c,d,sigc,sigd,chi2,q);
    //    fitting for Z=e+ft
    Linefit(3,1,e,f,sige,sigf,chi2,q);
    sigb2=sigb*sigb;
    sigd2=sigd*sigd;
    sigf2=sigf*sigf;
    
    apx_ix[1]=a;
    apx_ix[2]=c;
    apx_ix[3]=e;
    eapx_ix[1]=siga;
    eapx_ix[2]=sigc;
    eapx_ix[3]=sige;
    
    //      write(*,*) "b,d,f = ",b,d,f
    //		     cout<<"b= "<<b<<" d="<<d<<" f="<<f<<"\n";
    
    //      write(*,*) "sigb,sigd,sigf = ",sigb,sigd,sigf
    
    norm=sqrt(b*b+d*d+f*f);
    
    ctrk_ix[1]=b/norm;
    ctrk_ix[2]=d/norm;
    ctrk_ix[3]=f/norm;
    
    float norm3=norm*norm*norm;
    
    ectrk_ix[1]=sqrt((d*d+f*f)*(d*d+f*f)*sigb2+b*b*d*d*sigd2+b*b*f*f*sigf2)/norm3;
    ectrk_ix[2]=sqrt((b*b+f*f)*(b*b+f*f)*sigd2+d*d*b*b*sigb2+d*d*f*f*sigf2)/norm3;
    ectrk_ix[3]=sqrt((b*b+d*d)*(b*b+d*d)*sigf2+f*f*b*b*sigb2+f*f*d*d*sigd2)/norm3;
    
    return;
}


//------------------
// Linefit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
//------------------  
void DBCALShower_factory_KLOE::Linefit(int ixyz,int mwt,float &a,
                                  float &b,float &siga,float &sigb,float &chi2,float &q)
{
    
    // This programme is taken from Garth's book
    //  "Numerical Recipes in Fortran"
    //  and we tested it. It works well.  Chuncheng Xu,2005
    
    
    float sig[layermax_bcal+1],etemp;
    float xtemp[layermax_bcal+1],ytemp[layermax_bcal+1];
    // uses gammq
    
    //  Given a set of data points X(1:ndata),Y(1:ndata) with individual standard
    //  deviations sig(1:ndata), fit them to a strait line y=a+bx by minimum 
    // chisq. Returned are a,b and their respective probable uncertainties
    // siga and sigb, the chisq, and the goodness-of-fit probability q(that the
    // fit would have chisq this large or larger). if mwt=0 on input, then the 
    // the standard deviations are assumed to be unavalable:q is returned as 
    // 1.0 and the normalization of chi2 is to the unit standard deviation on all
    // points.
    
    
    float sigdat,ss,st2,sx,sxoss,sy,t,wt;
    sx=0.0;             // Initialize sums to zero
    sy=0.0;
    st2=0.0;
    b=0.0;
    
    int ndata=0;
    
    
    
    if(ixyz==1) {
        for (int i = 1; i < (layermax_bcal+1); i++){
            xtemp[i]=rt[i];
            ytemp[i]=x[i];
            sig[i]=sigx[i];
            etemp=e[i];
            if(etemp>0.0001)ndata=ndata+1;
            
        }   
    }
    else if(ixyz==2) {
        for (int i = 1; i < (layermax_bcal+1); i++){
            xtemp[i]=rt[i];
            ytemp[i]=y[i];
            sig[i]=sigy[i];
            etemp=e[i];
            if(etemp>0.000001)ndata=ndata+1;   
        }
    }
    else if(ixyz==3) {
        for (unsigned int i = 1; i < (layermax_bcal+1); i++){
            xtemp[i]=rt[i];
            ytemp[i]=z[i]; 
            sig[i]=sigz[i];  
            etemp=e[i];
            if(etemp>0.000001)ndata=ndata+1;
        }
    }
    
    if(mwt!=0) {   // Accumulate sums
        ss=0.0;
        for (int i = 1; i < (ndata+1); i++){
            wt=1.0/(sig[i]*sig[i]);
            ss=ss+wt;
            sx=sx+xtemp[i]*wt;
            sy=sy+ytemp[i]*wt;
        }
    }
    
    else  {
        for (int i = 1; i < (ndata+1); i++){
            sx=sx+xtemp[i];
            sy=sy+ytemp[i];
        } 
        ss=float(ndata);
    }
    
    sxoss=sx/ss;
    
    if(mwt!=0) {
        for (int i = 1; i < (ndata+1); i++){
            t=(xtemp[i]-sxoss)/sig[i];
            st2=st2+t*t;
            b=b+t*ytemp[i]/sig[i];
        }
        
    }
    
    else  {
        
        for (int i = 1; i < (ndata+1); i++){
            t=xtemp[i]-sxoss;
            st2=st2+t*t;
            b=b+t*ytemp[i];
        }
    }
    
    b=b/st2;                // Solve for a,b,siga and sigb
    a=(sy-sx*b)/ss;
    siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb=sqrt(1.0/st2);
    chi2=0.0;        //   Calculate chisq
    
    if(mwt==0) {
        for (int i = 1; i < (ndata+1); i++){
            chi2=chi2+(ytemp[i]-a-b*xtemp[i])*(ytemp[i]-a-b*xtemp[i]);
        }
        q=1.0;
        sigdat=sqrt(chi2/(ndata-2));
        siga=siga*sigdat;
        sigb=sigb*sigdat;
    }
    else  {
        for (int i = 1; i < (ndata+1); i++){
            chi2=chi2+((ytemp[i]-a-b*xtemp[i])/
                       sig[i])*((ytemp[i]-a-b*xtemp[i])/sig[i]);
        }
        q=Gammq(0.5*(ndata-2),0.5*chi2); 
    }
    
}



//------------------
// Gammq(a_gammq,x_gammq)
//------------------  
float DBCALShower_factory_KLOE::Gammq(float a_gammq,float x_gammq)
{
    
    //============================================================================
    
    
    float gammq;
    
    //    uses gcf,gser
    //Returns the incomplete gamma function Q(a_gammq,x_gammq)=1-P(a_gammq,x_gammq)
    
    
    float gammcf,gamser;
    
    if(a_gammq==0.0) { 
        gammq=999.0;
        return gammq;
    }
    
    
    
    if(x_gammq<0. || a_gammq<= 0.0) {
        //  cout<<"bad arguments in gammq"<<"\n";
        return 999.0;// pause 'bad arguments in gammq'
    }
    
    if(x_gammq<(a_gammq+1.)) {     // Use the series representation
        Gser(gamser,a_gammq,x_gammq);
        gammq=1.0-gamser;    // and take its complement
    }
    else  {               // Use continued fraction representation
        Gcf(gammcf,a_gammq,x_gammq);
        gammq=gammcf;
    }
    return gammq;
}


//------------------
// Gser(gamser,a_gser,x_gser,gln)
//------------------  
void DBCALShower_factory_KLOE::Gser(float &gamser,float a_gser,float x_gser)
{
    //============================================================================
    int itmax=100;
    float eps=3.0e-7;
    float gln;            
    
    //    uses gammln
    //    Returns the incomplete gamma function P(a,x) evaluated by its
    //    series representation as gamser. Also returns ln(Gamma(a)) as gln
    
	   float ap, del,sum;
       
       gln=Gammln(a_gser);
       
       if(x_gser<=0.0) {
           if(x_gser<0.0) cout<<"x_gser<0 in gser"<<"\n";
           gamser=0.0;
           return;
       }
       
       ap=a_gser;
       sum=1.0/a_gser;
       del=sum;
       
       
       for (int n = 1; n < (itmax+1); n++){
           ap=ap+1.0;
           del=del*x_gser/ap;
           sum=sum+del;
           
           if(fabs(del)<fabs(sum)*eps) {
               gamser=sum*exp(-x_gser+a_gser*log(x_gser)-gln);
               return;
           }
           
       }
       
       // cout<< "a too large, ITMAX too small in gser"<<"\n";
       return;
}


//------------------
//  Gcf(gammcf,a,x,gln)
//------------------  
void DBCALShower_factory_KLOE::Gcf(float &gammcf,float a_gcf,float x_gcf)
{
    //========================================================================
    
    int itmax=100;
    float eps=3.0e-7;
    float fpmin=1.0e-30;
    
    float gln;
    
    //   uses gammln
    //        Returns the incomplete gamma function Q(a,x) evaluated by 
    //   its continued fraction representation as gammcf. Also returns
    //   ln(Gamma(a)) as gln
    
    //   Parameters: ITMAX is the maximum allowed number of iterations;
    //   EPS is the relative accuracy; FPMIN is a number near the smallest
    //   representable floating-point number.
    
    float an,b,c,d,del,h;
    
    gln=Gammln(a_gcf);
    b=x_gcf+1.0-a_gcf;
    c=1.0/fpmin;
    d=1.0/b;
    h=d;
    
    
    for (int i = 1; i < (itmax+1); i++){
        an=-i*(i-a_gcf);
        b=b+2.0;
        d=an*d+b;
        if(fabs(d)<fpmin)d=fpmin;
        c=b+an/c;
        if(fabs(c)<fpmin)c=fpmin;
        d=1.0/d;
        del=d*c;
        h=h*del;
        if(fabs(del-1.0)<eps) {	 
            gammcf=exp(-x_gcf+a_gcf*log(x_gcf)-gln)*h; // Put factors in front 
            return;
        }
        // cout<< "a too large,ITMAX too small in gcf"<<"\n";    
        return;
    }
} //???

//------------------
// Gammln(xx)
//------------------  
float DBCALShower_factory_KLOE::Gammln(float xx_gln)
{
    //    Returns the value ln[Gamma(xx_gln)] for xx_gln>0.0
    
	   float ser,stp,tmp,x_gln,y_gln;
    float cof[7];
    float gammln; 
    
    //    Internal arithmetic will be done in double precision, a nicety that
    //    that you can omit if five-figure accuracy is good enoug
    //         save cof,stp
    //         data cof,stp/76.18009172947146d0,-86.50532032941677d0,
    //     * 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
    //     * -.5395239384953d-5,2.5066282746310005d0/
    
    stp=2.5066282746310005;
    cof[1]=76.18009172947146;
    cof[2]=-86.50532032941677;
    cof[3]=24.01409824083091;
    cof[4]=-1.231739572450155;
    cof[5]=.1208650973866179e-2;
    cof[6]=-.5395239384953e-5;
    
    x_gln=xx_gln;
    y_gln=x_gln;
    tmp=x_gln+5.5;
    tmp=(x_gln+0.5)*log(tmp)-tmp;
    
    ser=1.000000000190015;
    
    for (int j = 1; j < 7; j++){
        y_gln=y_gln+1.0;
        ser=ser+cof[j]/y_gln;
    }
    
    
    gammln=tmp+log(stp*ser/x_gln);
    return gammln;
}