cluster.cc

Go to the documentation of this file.

/*
A.Bazilevsky (IHEP, Protvino) - Dec. 1998
(from March 99 - RIKEN-BNL Research Center)

shura@bnl.gov

Cluster-Shower reconstruction routines based on the A.Lednev's 
method developed for GAMS (IHEP-CERN).

The Classes and Basic Functions description is in cluster.h

The method is based on the parameterized shower profile (energy deposited 
in the tower as a function of the distance between tower center and shower 
Center of Gravity, shower energy and impact angle). No difference between 
the shower produced by electrons or gammas is assumed.

************ The shower reconstruction follows three steps: *********

1.  Cluster search. A cluster is one or several neighboring cells 
separated from the other clusters with zero cells. The cluster are analyzed 
independently. 
Implementation: Class "cluster". The list of clusters is obtained via 
"Find_Clusters" function that has the list of hits as an input (look at 
"hit" structure description in cluster.h). For the brief description of 
the "cluster" methods look in cluster.h.

2.  Each cluster may contain several peaks. A peak is located in a cell whose 
deposited energy is higher then that of any adjoined cell. Peak regions are 
calculated by sharing the energy deposited in each cell according the  
energies expected from the gammas located in each peak region. This is 
iterative procedure. 
Implementation: Class "peakarea"; most of the methods are inherited from 
"cluster". The array of peakareas for the given cluster is obtained via 
"cluster::GetPeaks( peakarea*, hit* )" method. For the brief description of 
the "peakarea" methods look in cluster.h.

3.  Gamma reconstruction within the peak region. One or two gammas within a 
peak region. Minimizes Hi**2 function in two-dimensional space ( X,Y shower 
position ). If it is more then some limit, two-gamma hypothesis is taken. 
The two-gamma separation is done by minimizing Hi**2 function with the 
maximum slope method in three-dimensional space: a=(E1-E2)/E0, dX=X1-X2, 
dY=Y1-Y2. 
Implementation: Class "emshower". The array of emshower's for the given 
peakarea is obtained via "peakarea::GetGammas( emshower* )" method. 
For the brief description of the "emshower" methods look in cluster.h.

************* Main tips for using these library ***************

1. Parameters setting

SetThreshold(TowerThresh) - sets the threshold (GeV) in each tower; 
needed for the calculation of sigmas of the energy fluctuation in the 
towers; default value is 0.02 GeV;

SetChi2Limit(i) - sets Chi2 limit for the splitting of peakarea onto 
two showers; default value is 1;

2. Program flow

All sectors are processed independently. Everything starts with 
Find_Clusters call which returns the list of clusters. It has 
the list of hits as an input. Towers are counted from 0. 
Before every Find_Cluster call user should set the sector geometry 
and vertex position via SetGeometry( SecGeom Geom, float* pVert) call. 

As soon as one gets the list of clusters, he can do whatever he wants: 
get cluster's characteristics, split each of them onto peakarea's, 
get peakarea's characteristics, split each of them onto emshower's, 
get emshower's characteristics. 

3. Shower Profile

The parameterized Shower Profile is used when splitting cluster onto 
peakareas and peakarea onto emshowers. The parameters for the Shower 
Profile calculation are set via SetProfileParameters(...) function for 
each peakarea object. The energy deposited in the particular tower is 
obtained via EMCcell(...) function (based on the Shower Profile)

*/

#include <math.h>
#include <stdio.h>
#include "cluster.h"

// If it's defined the cluster with adjacent corner are glued;
// otherwise the cluster with adjacent border are glued
//#define GLUE_CLUSTER_CORNER

// Minimal Peak energy when looking for peakarea's; used in cluster::GetPeaks
#define MinPeakEnergy 0.08

// Minimal shower energy when splitting peakarea onto showers; used in gamma(...) function
#define MinShowerEnergy 0.1

// Max number of clusters in sector; used in Find_Clusters(...)
#define MaxLen 1000

// Max number of peaks in cluster; used in cluster::GetPeaks(...)
#define MaxNofPeaks 10

// Used in cluster::GetPeaks(...)
#define PeakIter 6

// Emin cuts cluster energy: Ecl >= Epk1+Epk2 +...+Epkn !!!!!!!!!!!!!
#define Emin 0.002

// chisq=3 devides about 1% of single showers (now it isn't used)
#define chisq 3.

// Parameters for sigma in Hi2 calculations (p.36-37 v.3)
float epar00 = 0.005;
float epar0 = 0.0014;
float epar1 = 0.03;
float epar2 = -0.03;
float epar3 = 0.;
float epar4 = 4.0;

// define meaningless values for (x,y)
#define XABSURD -999999.
#define YABSURD -999999.

int EMCzSize = 72;
int EMCySize = 36;
float Tower_zSize = 5.535; // Tower size in z dir (cm)
float Tower_ySize = 5.535; // Tower size in y dir (cm)

float xVertex=0, yVertex=0, zVertex=0;
float xSector=500, ySector=-97, zSector=-196;
float normx=1, normy=0, normz=0;

// Default level, now for the Conf Level: 1% for GEANT, 2%-5% for TestBeam
float Chi2Level[50]={
    6.634899, 4.605171, 3.780564, 3.318915, 3.017103,    
    2.801872, 2.639259, 2.511249, 2.407341, 2.320967,    
    2.247720, 2.184744, 2.129863, 2.081515, 2.038526,    
    1.999994, 1.965214, 1.933627, 1.904781, 1.878311,    
    1.853912, 1.831334, 1.810365, 1.790825, 1.772564, 
    1.755449, 1.739367, 1.724222, 1.709926, 1.696406,    
    1.683593, 1.671430, 1.659864, 1.648850, 1.638344,    
    1.628311, 1.618716, 1.609528, 1.600721, 1.592268,    
    1.584148, 1.576338, 1.568822, 1.561579, 1.554596,    
    1.547856, 1.541346, 1.535055, 1.528968, 1.523077 };   

// For the Conf Level: 1% for GEANT, 2%-5% for TestBeam
float Chi2Level1[50]={
    6.634899, 4.605171, 3.780564, 3.318915, 3.017103,    
    2.801872, 2.639259, 2.511249, 2.407341, 2.320967,    
    2.247720, 2.184744, 2.129863, 2.081515, 2.038526,    
    1.999994, 1.965214, 1.933627, 1.904781, 1.878311,    
    1.853912, 1.831334, 1.810365, 1.790825, 1.772564, 
    1.755449, 1.739367, 1.724222, 1.709926, 1.696406,    
    1.683593, 1.671430, 1.659864, 1.648850, 1.638344,    
    1.628311, 1.618716, 1.609528, 1.600721, 1.592268,    
    1.584148, 1.576338, 1.568822, 1.561579, 1.554596,    
    1.547856, 1.541346, 1.535055, 1.528968, 1.523077 };   

// For the Conf Level: 2% for GEANT, 4%-7% for TestBeam
float Chi2Level2[50]={
    5.411895, 3.912024, 3.278443, 2.916812, 2.677547,    
    2.505458, 2.374582, 2.271008, 2.186567, 2.116065,    
    2.056169, 2.004491, 1.959343, 1.919481, 1.883964,    
    1.852072, 1.823237, 1.797008, 1.773021, 1.750981,    
    1.730640, 1.711795, 1.694274, 1.677931, 1.662643,    
    1.648301, 1.634814, 1.622101, 1.610093, 1.598727,    
    1.587948, 1.577709, 1.567968, 1.558684, 1.549824,    
    1.541357, 1.533256, 1.525494, 1.518051, 1.510903,    
    1.504033, 1.497424, 1.491059, 1.484924, 1.479006,    
    1.473292, 1.467771, 1.462433, 1.457267, 1.452265 };

float sin4a, sinax, sinay;
float ppar1, ppar2, ppar3, ppar4, pShiftx, pShifty;

#define max(a,b)  ((a) > (b) ? (a) : (b))
#define min(a,b)  ((a) < (b) ? (a) : (b))
#define ABS(x)  ((x) < 0 ? -(x) : (x))
#define lowint(x) ((x) < 0 ? int(x-1) : int(x))
#define asinht(x) (log((x)+sqrt((x)*(x)+1)))
#define sinht(x) ((exp(x)-exp(-(x)))/2.)

int Hit_NCompare( const void*, const void* );
int Hit_ACompare( const void*, const void* );
void CopyVector( int*, int*, int );
void CopyVector( hit*, hit*, int );
void ZeroVector( int*, int );
void ZeroVector( float*, int );
void ZeroVector( hit*, int );
void ResizeVector( int* , int, int );
void CorrectPosition( float, float, float, float*, float* );
void CalculateErrors( float e, float x, float y, float* pde, float* pdx, float* pdy, float* pdz);
void GlobalToSector( float, float, float, float*, float*, float* );
void SectorToGlobal( float xsec, float ysec, float zsec, float* px, float* py, float* pz );
void SectorToGlobalErr( float dxsec, float dysec, float dzsec, float* pdx, float* pdy, float* pdz );
void gamma(int, hit*, float*, float*, float*, float*, float*, float*, float*, float*);
void twogamma(int, hit*, float*, float*, float*, float*, float*, float*, float*);
void mom1(int, hit*, float*, float*, float*);
void momenta(int, hit*, float*, float*, float*, float*, float*, float* );
void c3to5(float, float, float, float, float, float, float*, float*, float*,  float*, float*, float*);
void SetProfileParameters(int, float, float, float);
float ClusterChisq(int, hit*, float, float, float);
float Chi2Limit( int ND );

void emshower::GetCorrPos( float* px, float* py )
// Returns emshower corrected position in Sector (SM) frame
{
   CorrectPosition( energy, xcg, ycg, px, py );
}

void emshower::GetGlobalPos( float* px, float* py, float* pz )
// Returns emshower position in PHENIX global coord system
{
   float xc, yc;

   GetCorrPos( &xc, &yc );
// X in Sector coord is Z in Global coord !!
   SectorToGlobal( 0, yc, xc, px, py, pz );
}

void cluster::GetCorrPos( float* px, float* py )
// Returns the cluster corrected position in Sector (SM) frame
{
   float e, x, y, xx, yy, xy;

   e = GetTotalEnergy();
   GetMoments( &x, &y, &xx, &xy, &yy );
   CorrectPosition( e, x, y, px, py );
}

void cluster::GetGlobalPos( float* px, float* py, float* pz )
// Returns the cluster position in PHENIX global coord system
{
   float xc, yc;

   GetCorrPos( &xc, &yc );
// X in Sector coord is Z in Global coord !!
   SectorToGlobal( 0, yc, xc, px, py, pz );
}

float cluster::GetTowerEnergy( int ich )
// Returns the energy of the ich-tower (0 if ich not found in the HitList)
{
     list<hit>::iterator ph;
     if( HitList.empty() ) return 0;
     ph = HitList.begin();
     while( ph != HitList.end() ) {
       if( (*ph).ich == ich ) return (*ph).amp;
        ph++;
     }
     return 0;
}

float cluster::GetTowerToF( int ich )
// Returns the ToF of the ich-tower (0 if ich not found in the HitList)
{
     list<hit>::iterator ph;
     if( HitList.empty() ) return 0;
     ph = HitList.begin();
     while( ph != HitList.end() ) {
       if( (*ph).ich == ich ) return (*ph).tof;
        ph++;
     }
     return 0;
}

float cluster::GetTotalEnergy()
// Returns the cluster total energy
{
     list<hit>::iterator ph;
     float et=0;
     if( HitList.empty() ) return 0;
     ph = HitList.begin();
     while( ph != HitList.end() ) { 
        et += (*ph).amp;
        ph++;
     }
     return et;
}

float cluster::GetE9( int ich )
// Returns the energy in 3x3 towers around the tower ich
{
     int i, ic, ixc, ix0, ichmax, ichmin, nhit;
     float e;
     hit *hlist, *vv;
     list<hit>::iterator ph;

     e=0;

     nhit = HitList.size();

     if( nhit <= 0 ) return e;

     hlist = new hit[nhit];

     ph = HitList.begin();
     vv = hlist;
     while( ph != HitList.end() ) *vv++ = *ph++;

     qsort( hlist, nhit, sizeof(hit), Hit_NCompare );

     ichmax = ich + EMCzSize + 1;
     ichmin = ich - EMCzSize - 1;
     ix0 = ich - ich/EMCzSize*EMCzSize;

     for( i=0; i<nhit; i++ ) 
        if( hlist[i].ich >= ichmin ) break;
     while( (i < nhit) && ( (ic=hlist[i].ich) <= ichmax) ) {
        ixc = ic - ic/EMCzSize*EMCzSize;
        if( abs(ixc-ix0) <= 1 ) e += hlist[i].amp;
        i++;
     }
        
     delete hlist;
     return e;
}

hit cluster::GetImpactTower()
// Returns the hit corresponding to the reconstructed impact tower
{
     float x, y;
     int ix, iy, ich;
     hit ht;

     GetCorrPos( &x, &y );
     ix = lowint(x/Tower_zSize + 0.5);
     iy = lowint(y/Tower_ySize + 0.5);
     if( ix < 0 || ix > EMCzSize-1 || iy < 0 || iy > EMCySize-1 ) {
        printf("????? EmcClusterChi2: Something wrong in GetImpactTower: (x,y)=(%f,%f)  (ix,iy)=(%d,%d) \n", x, y, ix, iy);
        ht.ich = -1;
        ht.amp = 0;
        ht.tof = 0;
        return ht;
     }
     else {
        ich = iy*EMCzSize + ix;
        ht.ich = ich;
        ht.amp = GetTowerEnergy( ich );
        ht.tof = GetTowerToF( ich );
//      printf("Energy=%f  ich=%f\n",*pEnergy,*pich);
        return ht;
     }
}

hit cluster::GetMaxTower()
// Returns the hit with the maximum energy
{
     list<hit>::iterator ph;
     float emax=0;
     hit ht;

     ht.ich = -1;
     ht.amp = 0;
     ht.tof = 0;
     if( HitList.empty() ) return ht;

     ph = HitList.begin();
     while( ph != HitList.end() ) { 
        if( (*ph).amp > emax ) { emax = (*ph).amp; ht = *ph; }
        ph++;
     }
     return ht;
}

void cluster::GetHits(hit* phit, int n)
// Returns n hits (sorted) with the maximum energy
{
     int nhit;
     hit *hlist, *vv;
     list<hit>::iterator ph;
     list<hit> hl;

     ZeroVector( phit, n );
     nhit = HitList.size();

     if( nhit <= 0 ) return;

     hlist = new hit[nhit];

     ph = HitList.begin();
     vv = hlist;
     while( ph != HitList.end() ) *vv++ = *ph++;

     qsort( hlist, nhit, sizeof(hit), Hit_ACompare );
     for( int i=0; i<min(nhit,n); i++ ) phit[i]=hlist[i];
     delete hlist;
}

void cluster::GetMoments( float* px, float* py, float* pxx, float* pxy, float* pyy )
//  Returns cluster 1-st (px,py) and 2-d momenta (pxx,pxy,pyy)
{
     list<hit>::iterator ph;
     float e, x, y, xx, yy, xy;
     int nhit;
     hit *phit, *p;

     *px = XABSURD; *py = YABSURD;
     *pxx = 0; *pxy = 0; *pyy = 0;
     nhit=HitList.size();
     if( nhit <= 0 ) return;

     phit = new hit[nhit];
     ph = HitList.begin();
     p=phit;
     while( ph != HitList.end() ) { 
        p->ich = (*ph).ich;
        p->amp = (*ph).amp;
        ph++;
        p++;
     }
     momenta( nhit, phit, &e, &x, &y, &xx, &yy, &xy );
     *px = x*Tower_zSize;
     *py = y*Tower_ySize;
     *pxx = xx*Tower_zSize*Tower_zSize;
     *pxy = xy*Tower_zSize*Tower_ySize;
     *pyy = yy*Tower_ySize*Tower_ySize;
//      cout<<*pxx<<"  "<<*pyy<<"  "<<*pxy<<'\n';
     delete phit;
}

void cluster::GetErrors( float* pde, float* pdx, float* pdy, float* pdz)
//  Returns the errors for the reconstructed energy and position
{
    float e, x, y;

    e = GetTotalEnergy();
    GetCorrPos( &x, &y );
    CalculateErrors( e, x, y, pde, pdx, pdy, pdz);
}

void emshower::GetErrors( float* pde, float* pdx, float* pdy, float* pdz)
// Returns the errors for the reconstructed energy and position
{
    float e, x, y;

    e = GetTotalEnergy();
    GetCorrPos( &x, &y );
    CalculateErrors( e, x, y, pde, pdx, pdy, pdz);
}

void cluster::GetChar(  float* pe, 
                        float* pxcg, float* pycg, 
                        float* pxc, float* pyc, 
                        float* pxg, float* pyg, float* pzg,
                        float* pxx, float* pxy, float* pyy,
                        float* pde, float* pdx, float* pdy, float* pdz  )
// This method replaces methods GetTotalEnergy, GetMoments, 
// GetCorrPos, GetGlobalPos, GetErrors
{
    *pe = GetTotalEnergy();
    GetMoments( pxcg, pycg, pxx, pxy, pyy );
    CorrectPosition( *pe, *pxcg, *pycg, pxc, pyc );
    SectorToGlobal( 0, *pyc, *pxc, pxg, pyg, pzg );
    CalculateErrors( *pe, *pxc, *pyc, pde, pdx, pdy, pdz);
}

void peakarea::GetChar( float* pe, 
                        float* pxcg, float* pycg, 
                        float* pxcgmin, float* pycgmin, 
                        float* pxc, float* pyc, 
                        float* pxg, float* pyg, float* pzg,
                        float* pxx, float* pxy, float* pyy,
                        float* pchi,
                        float* pde, float* pdx, float* pdy, float* pdz  )
// This method replaces "cluster" methods GetTotalEnergy, GetMoments,
// GetCorrPos, GetGlobalPos, GetErrors and "peakarea" methods GetCGmin, GetChi2
{
    float chi, chi0, chicorr;
    float e1, x1, y1, e2, x2, y2;
    int nh;
    hit *phit, *vv;
    list<hit>::iterator ph;
    list<hit> hl;

    *pe = 0;
    hl = HitList;
    nh = hl.size();
    if( nh <= 0 ) return;

    phit = new hit[nh];

    ph = hl.begin();
    vv = phit;
    while( ph != hl.end() ) *vv++ = *ph++;

    chi = chisq*1000;
    gamma(nh,phit,&chi,&chi0,&e1,&x1,&y1,&e2,&x2,&y2);
    chicorr = Chi2Correct(chi0,nh);
    //    chi0 = chicorr;
    *pchi = chi0;
    // Shower Center of Gravity after shower profile fit
    *pxcgmin = x1*Tower_zSize;
    *pycgmin = y1*Tower_ySize;

    *pe = GetTotalEnergy();
    GetMoments( pxcg, pycg, pxx, pxy, pyy );
//    CorrectPosition( *pe, *pxcg, *pycg, pxc, pyc );
    CorrectPosition( *pe, *pxcgmin, *pycgmin, pxc, pyc );
    SectorToGlobal( 0, *pyc, *pxc, pxg, pyg, pzg );
    CalculateErrors( *pe, *pxc, *pyc, pde, pdx, pdy, pdz);

    delete phit;
}


void emshower::GetChar( float* pe, 
                        float* pxcg, float* pycg, 
                        float* pxc, float* pyc, 
                        float* pxg, float* pyg, float* pzg, 
                        float* pchi,
                        float* pde, float* pdx, float* pdy, float* pdz  )
// This method replaces methods GetTotalEnergy, GetCG, 
// GetCorrPos, GetGlobalPos, GetChi2, GetErrors
{
    *pe = GetTotalEnergy();
    GetCG( pxcg, pycg );
    CorrectPosition( *pe, *pxcg, *pycg, pxc, pyc );
    SectorToGlobal( 0, *pyc, *pxc, pxg, pyg, pzg );
    CalculateErrors( *pe, *pxc, *pyc, pde, pdx, pdy, pdz);
    *pchi = GetChi2();
}


float peakarea::GetChi2()
// Returns Hi2 after its minimization fluctuating CG position 
// (i.e. after shower profile fit)
{
    float chi, chi0, chicorr;
    float e1, x1, y1, e2, x2, y2;
    int nh;
    hit *phit, *vv;
    list<hit>::iterator ph;
    list<hit> hl;

    hl = HitList;
    nh = hl.size();
    if( nh <= 0 ) return 0;

    phit = new hit[nh];

    ph = hl.begin();
    vv = phit;
    while( ph != hl.end() ) *vv++ = *ph++;

    chi = chisq*1000;
    gamma(nh,phit,&chi,&chi0,&e1,&x1,&y1,&e2,&x2,&y2);

    chicorr = Chi2Correct(chi0,nh);
    //    chi0 = chicorr;

    delete phit;
    return chi0;
}

void peakarea::GetCGmin( float* px, float* py )
// Gets CG coordinates corresponding to min Hi2 (after shower shape fit)
{
    float chi, chi0;
    float e1, x1, y1, e2, x2, y2;
    int nh;
    hit *phit, *vv;
    list<hit>::iterator ph;
    list<hit> hl;

    *px = XABSURD;
    *py = YABSURD;
    hl = HitList;
    nh = hl.size();
    if( nh <= 0 ) return;

    phit = new hit[nh];

    ph = hl.begin();
    vv = phit;
    while( ph != hl.end() ) *vv++ = *ph++;

    chi = chisq*1000;
    gamma(nh,phit,&chi,&chi0,&e1,&x1,&y1,&e2,&x2,&y2);
    *px = x1*Tower_zSize;
    *py = y1*Tower_ySize;

    delete phit;
}

int peakarea::GetGammas( emshower* ShList)
//
// Splits the peakarea onto 1 or 2 emshower's
// Returns Number of Showers (0, 1 or 2)
//
// If the Chi2 of the peakarea is less then the value determined 
// in Chi2Limit(Number_of_Hits) function, 1-photon hypothesis is 
// accepted, otherwise the 2-shower hypothesis is checked
//
{
     float chi, chi0, e1, x1, y1, e2, x2, y2, chicorr;
     int nh, ig;
     hit *phit, *vv;
     emshower sh1, sh2;
     list<hit>::iterator ph;
     list<hit> hl;

//     (*ShList).erase( (*ShList).begin(), (*ShList).end() );

     hl = HitList;
     nh = hl.size();
     if( nh <= 0 ) return 0;

     phit = new hit[nh];

     ph = hl.begin();
     vv = phit;
     while( ph != hl.end() ) *vv++ = *ph++;

//     chi=chisq;
     chi=Chi2Limit(nh);
     gamma(nh,phit,&chi,&chi0,&e1,&x1,&y1,&e2,&x2,&y2);
//      printf("GetGammas: e,x,y=(%f,%f,%f) ,y=(%f,%f,%f) \n",e1,x1,y1,e2,x2,y2);
     chicorr = Chi2Correct(chi,nh);
     //     chi = chicorr;
     if( e1 > 0 )
         sh1.ReInitialize( e1, x1*Tower_zSize, y1*Tower_ySize, chi );
     else
         sh1.ReInitialize( 0, XABSURD, YABSURD, 0 );

     if( e2 > 0 )
         sh2.ReInitialize( e2, x2*Tower_zSize, y2*Tower_ySize, chi );
     else
         sh2.ReInitialize( 0, XABSURD, YABSURD, 0 );

     if( e2 > e1 ) {
       sh1.ReInitialize( e2, x2*Tower_zSize, y2*Tower_ySize, chi );
       sh2.ReInitialize( e1, x1*Tower_zSize, y1*Tower_ySize, chi );
     }

     ig = 0;
     if( e1>0 ) {
       ShList[ig++]=sh1;
       if( e2>0 ) {
         ShList[ig++]=sh2;
       }
     }

     delete phit;
     return ig;
}

int cluster::GetPeaks( peakarea* PkList, hit* ppeaks )
//
// Splits the cluster onto peakarea's
// The number of peakarea's is equal to the number of Local Maxima 
// in the cluster. Local Maxima can have the energy not less then 
// defined in MinPeakEnergy
//
// Output: PkList - the array of peakarea's
//         ppeaks - the array of peak hits (one for each peakarea)
//
// Returns: >= 0 Number of Peaks;
//            -1 The number of Peaks is greater then MaxNofPeaks;
//               (just increase parameter MaxNofPeaks)
//
{
     int npk, ipk, nhit;
     int ixypk, ixpk, iypk, in, nh, ic;
     int ixy, ix, iy, nn;
     int ig, ng, igmpk1[MaxNofPeaks], igmpk2[MaxNofPeaks];
     int PeakCh[MaxNofPeaks];
     float epk[MaxNofPeaks*2], xpk[MaxNofPeaks*2], ypk[MaxNofPeaks*2];
     float ratio, eg, dx, dy, a, chi, chi0;
     float *Energy[MaxNofPeaks], *totEnergy, *tmpEnergy;
     hit *ip;
     hit *phit, *hlist, *vv;
     peakarea peak;
     list<hit>::iterator ph;
     list<hit> hl;

//     (*PkList).erase( (*PkList).begin(), (*PkList).end() );
     nhit = HitList.size();

     if( nhit <= 0 ) return 0;

     hlist = new hit[nhit];

     ph = HitList.begin();
     vv = hlist;
     while( ph != HitList.end() ) *vv++ = *ph++;

     qsort( hlist, nhit, sizeof(hit), Hit_NCompare );
//
//  Find peak (maximum) position (towers with local maximum amp)
//
     npk=0;
     for( ic=0; ic<nhit; ic++ ) {
        float amp = hlist[ic].amp;
        if( amp > MinPeakEnergy ) {
           int ich = hlist[ic].ich;
           int ichmax = ich + EMCzSize + 1;
           int ichmin = ich - EMCzSize - 1;
           int ixc = ich - ich/EMCzSize*EMCzSize;
           int in = ic + 1;
           while( (in < nhit) && (hlist[in].ich <= ichmax) ) {
              int ix = hlist[in].ich - hlist[in].ich/EMCzSize*EMCzSize;
              if( (abs(ixc-ix) <= 1) && (hlist[in].amp >= amp) ) goto new_ic;
              in++;
           }
           in = ic - 1;
           while( (in >= 0) && (hlist[in].ich >= ichmin) ) {
              int ix = hlist[in].ich - hlist[in].ich/EMCzSize*EMCzSize;
              if( (abs(ixc-ix) <= 1) && (hlist[in].amp > amp) ) goto new_ic;
              in--;
           }
           if( npk >= MaxNofPeaks ) { delete hlist; return -1; }
           PeakCh[npk]=ic;
//      cout << npk << "   "<<PeakCh[npk]<<"\n";
           npk++;
        }
new_ic: continue;
     }

     if( npk <= 0 ) { delete hlist; return 0; }
     if( npk == 1 ) {
        *ppeaks = hlist[PeakCh[0]];
        hl.erase( hl.begin(), hl.end() );
        for( int ich=0; ich<nhit; ich++ ) hl.push_back(hlist[ich]);
        peak.ReInitialize(hl);
//      PkList->push_back(peak);
        PkList[0]=peak;
        delete hlist;
        return 1;
     }

     for( ipk=0; ipk<npk; ipk++ ) { Energy[ipk] = new float[nhit]; }
     totEnergy = new float[nhit];
     tmpEnergy = new float[nhit];
//
// Divide energy in towers among photons positioned in every peak
//
     ratio = 1.;
     for( int iter=0; iter<PeakIter; iter++ ) {
        ZeroVector( tmpEnergy, nhit );
        for( ipk=0; ipk<npk; ipk++ ) {

           ic = PeakCh[ipk];
           if( iter > 0 ) ratio = Energy[ipk][ic]/totEnergy[ic];
           eg = hlist[ic].amp * ratio;
           ixypk = hlist[ic].ich;
           iypk = ixypk/EMCzSize;
           ixpk = ixypk - iypk*EMCzSize;
           epk[ipk] = eg;
           xpk[ipk] = eg * ixpk;
           ypk[ipk] = eg * iypk;

           if( ic < nhit-1 ){
             for( in=ic+1; in<nhit; in++ ) {
                ixy = hlist[in].ich;
                iy = ixy/EMCzSize;
                ix = ixy - iy*EMCzSize;
                if( (ixy-ixypk) > EMCzSize+1 ) break;
                if( abs(ix-ixpk) <= 1 ) {
                   if( iter > 0 ) ratio = Energy[ipk][in]/totEnergy[in];
                   eg = hlist[in].amp * ratio;
                   epk[ipk] += eg;
                   xpk[ipk] += eg*ix;
                   ypk[ipk] += eg*iy;
                }
             }
           } // if ic

           if( ic >= 1 ){
             for( in=ic-1; in >= 0; in-- ) {
                ixy = hlist[in].ich;
                iy = ixy/EMCzSize;
                ix = ixy - iy*EMCzSize;
                if( (ixypk-ixy) > EMCzSize+1 ) break;
                if( abs(ix-ixpk) <= 1 ) {
                   if( iter > 0 ) ratio = Energy[ipk][in]/totEnergy[in];
                   eg = hlist[in].amp * ratio;
                   epk[ipk] += eg;
                   xpk[ipk] += eg*ix;
                   ypk[ipk] += eg*iy;
                }
             }
           } // if ic

           xpk[ipk] = xpk[ipk]/epk[ipk];
           ypk[ipk] = ypk[ipk]/epk[ipk];
           SetProfileParameters(0, epk[ipk], xpk[ipk], ypk[ipk]);
           for( in=0; in<nhit; in++ ) {
              ixy = hlist[in].ich;
              iy = ixy/EMCzSize;
              ix = ixy - iy*EMCzSize;
              dx = xpk[ipk]-ix;
              dy = ypk[ipk]-iy;
              a = 0;
              if( ABS(dx) < 2.5 && ABS(dy) < 2.5 ) a = epk[ipk]*EMCcell(dx, dy, -1);
              Energy[ipk][in] = a;
              tmpEnergy[in] += a;
           }

        } // for ipk

        for( in=0; in<nhit; in++ ) totEnergy[in] = max(0.000001,tmpEnergy[in]);

     } // for iter

     phit = new hit[nhit];

     ng=0;
     for( ipk=0; ipk<npk; ipk++ ) {
        nh=0;
        for( in=0; in<nhit; in++ ) {
           if( tmpEnergy[in] > 0 ) {
              ixy = hlist[in].ich;
              a = hlist[in].amp * Energy[ipk][in]/tmpEnergy[in];
              if( a > Emin ) {
                phit[nh].ich=ixy;
                phit[nh].amp=a;
                phit[nh].tof= hlist[in].tof;  // Not necessary here
                nh++;
              }
           }
        }
        if( nh>0 ) {
//         chi=chisq*10;
           chi=Chi2Limit(nh)*10;
           gamma(nh,phit,&chi,&chi0,&epk[ng],&xpk[ng],&ypk[ng],&epk[ng+1],&xpk[ng+1],&ypk[ng+1]);
           igmpk1[ipk]=ng;
           igmpk2[ipk]=ng;
           if( epk[ng+1]>0 ) { ng++; igmpk2[ipk]=ng;}
           ng++;
        }
     }

     ZeroVector( tmpEnergy, nhit );
     for( ipk=0; ipk<npk; ipk++ ) {
        ig=igmpk1[ipk];
        SetProfileParameters(0, epk[ig], xpk[ig], ypk[ig]);
        for( in=0; in<nhit; in++ ) {
           Energy[ipk][in]=0;
           for(ig=igmpk1[ipk]; ig<=igmpk2[ipk]; ig++){
              ixy = hlist[in].ich;
              iy = ixy/EMCzSize;
              ix = ixy - iy*EMCzSize;
              dx = xpk[ig]-ix;
              dy = ypk[ig]-iy;
              a = epk[ig]*EMCcell(dx, dy, epk[ig]);
              Energy[ipk][in] += a;
              tmpEnergy[in] += a;
           }
        }
     }

     ip = ppeaks;
     nn=0;
     for( ipk=0; ipk<npk; ipk++ ) {
        nh=0;
        for( in=0; in<nhit; in++ ) {
           if( tmpEnergy[in] > 0 ) {
              ixy = hlist[in].ich;
              a = hlist[in].amp * Energy[ipk][in]/tmpEnergy[in];
              if( a > Emin ) {
                phit[nh].ich=ixy;
                phit[nh].amp=a;
                phit[nh].tof=hlist[in].tof;
                nh++;
              }
           }
        }
        if( nh>0 ) {
           *ip++ = hlist[PeakCh[ipk]];
           hl.erase( hl.begin(), hl.end() );
           for( in=0; in<nh; in++ ) hl.push_back(phit[in]);
           peak.ReInitialize(hl);
//         PkList->push_back(peak);
           PkList[nn++]=peak;
        }
     }

     delete phit;
     delete hlist;
     for( ipk=0; ipk<npk; ipk++ ) delete Energy[ipk];
     delete totEnergy;
     delete tmpEnergy;

//     return npk;
     return nn;
}

void gamma(int nh, hit* phit, float* pchi, float* pchi0, float* pe1, float* px1, float* py1, float* pe2, float* px2, float* py2)
//
// Tests for 1 or 2 photon hypothesis by minimizing the chi2.
// If the energy of one of two showers is less then MinShowerEnergy 
// they are merged
//
// Returns two shower parameters (pe1,px1,py1,pe2,px2,py2). 
// If one shower hypothesis accepted -> *pe2=0
// *pchi  - chi2 after splitting
// *pchi0 - chi2 before splitting
//
{
     float e1, x1, y1, e2, x2, y2;
     float chi, chi0, chi00, chisq0, chisave;
     float chir, chil, chiu, chid;
     int dof;
     float x0, y0, d2, xm2;
     float stepx, stepy, parx, pary;
     const float dxy=0.06;
     const float stepmin=0.01;
     const float zTG=100;
     const float xmcut=0.0015; // (GeV), for overlapping showers separation

     *pe1=0;
     *px1=0;
     *py1=0;
     *pe2=0;
     *px2=0;
     *py2=0;
     if( nh <= 0 ) return;
     mom1(nh,phit,&e1,&x1,&y1);
     *pe1=e1;     
     if( e1 <= 0 ) return;

     SetProfileParameters(0, e1,x1,y1);
     chisave = *pchi;
     chi = *pchi;
     chi0 = ClusterChisq(nh,phit,e1,x1,y1);
//      cout << " ChiSq=" << chi0 << "  " << e1 << "  " << x1 << "  " << y1 << '\n';
     chisq0 = chi0;
//     dof = nh-2;
     dof = nh;
     if( dof < 1 ) dof=1;
     chi = chisq0/dof;
     x0 = x1;
     y0 = y1;
     for(;;){
        chir = ClusterChisq(nh,phit,e1,x0+dxy,y0);
        chil = ClusterChisq(nh,phit,e1,x0-dxy,y0);
        chiu = ClusterChisq(nh,phit,e1,x0,y0+dxy);
        chid = ClusterChisq(nh,phit,e1,x0,y0-dxy);

        if( (chi0 > chir) || (chi0 > chil) ) {
           stepx = dxy;
           if( chir > chil ) stepx = -stepx;
        }
        else {
           stepx = 0;
           parx = chir+chil-2*chi0;
           if( parx > 0 ) stepx = -dxy*(chir-chil)/2/parx;
        }

        if( (chi0 > chiu) || (chi0 > chid) ) {
           stepy = dxy;
           if( chiu > chid ) stepy = -stepy;
        }
        else {
           stepy = 0;
           pary = chiu+chid-2*chi0;
           if( pary > 0 ) stepy = -dxy*(chiu-chid)/2/pary;
        }
        if( (ABS(stepx) < stepmin) && (ABS(stepy) < stepmin) ) break;
        chi00 = ClusterChisq(nh,phit,e1,x0+stepx,y0+stepy);
        if( chi00 >= chi0 ) break;
        chi0 = chi00;
        x0 += stepx;
        y0 += stepy;
     }
     if( chi0 < chisq0 ) {
        x1 = x0;
        y1 = y0;
        chi = chi0/dof;
     }

     *pchi0 = chi;
     *pchi = chi;
     *px1 = x1;
     *py1 = y1;

     if( chi > chisave ) {
        twogamma(nh,phit,&chi,&e1,&x1,&y1,&e2,&x2,&y2);
        if( e2 > 0 ) {
           d2 = ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))/zTG/zTG;
           xm2 = e1*e2*d2;
           if( xm2 > 0 ) xm2 = sqrt(xm2);
           if( xm2 > xmcut && e1 > MinShowerEnergy && e2 > MinShowerEnergy) {
                *pe1 = e1;
                *px1 = x1;
                *py1 = y1;
                *pe2 = e2;
                *px2 = x2;
                *py2 = y2;
                *pchi = chi;
           }
        }       
     }
//      printf("%f %f %f\n",*px1*Tower_zSize, *py1*Tower_ySize, *pchi);
}

void twogamma(int nh, hit* phit, float* pchi, float* pe1, float* px1, float* py1, float* pe2, float* px2, float* py2)
{
     float e0, x0, y0, xx, yy, yx;
     float dxy, rsg2, rsq;
     float dxc, dyc, r, epsc;
     int ix, iy, ixy, in, iter, dof;
     float step, cosi, chisq2, u;
     float e1c, x1c, y1c, e2c, x2c, y2c;
     float eps0, eps1, eps2, chisqc, ex;
     float dx1, dy1, dx2, dy2, a0, d;
     float dchi, dchi0, dd, dchida, gr, a1, a2;
     float gre, grx, gry, grec, grxc, gryc, grc, gx1, gx2, gy1, gy2;
     float scal, dx0, dy0;
     const float epsmax=0.9999;
     const float stpmin=0.025;
     const float delch=2;

     momenta(nh,phit,&e0,&x0,&y0,&xx,&yy,&yx);
     *pe2 = 0;
     *px2 = 0;
     *py2 = 0;
     if( nh <= 0 ) return;
//  choosing of the starting point
     dxy = xx-yy;
     rsg2 = dxy*dxy + 4*yx*yx;
     if( rsg2 < 1e-20 ) rsg2 = 1e-20;
     rsq = sqrt(rsg2);
     dxc = -sqrt((rsq+dxy)*2);
     dyc =  sqrt((rsq-dxy)*2);
     if( yx >= 0 ) dyc = -dyc;
     r = sqrt(dxc*dxc + dyc*dyc);
     epsc = 0;
     for( in=0; in<nh; in++ ) {
        ixy = phit[in].ich;
        iy = ixy/EMCzSize;
        ix = ixy - iy*EMCzSize;
        u = (ix-x0)*dxc/r + (iy-y0)*dyc/r;
        epsc -= phit[in].amp * u * ABS(u);
     }
     epsc /= (e0*rsq);
     if( epsc >  0.8 ) epsc = 0.8;
     if( epsc < -0.8 ) epsc =-0.8;
     dxc /= sqrt(1-epsc*epsc);
     dyc /= sqrt(1-epsc*epsc);
//  Start of iterations
     step = 0.1;
     cosi = 0;
     chisq2 = 1.e35;
     for( iter=0; iter<100; iter++)
     {
        c3to5(e0,x0,y0,epsc,dxc,dyc,&e1c,&x1c,&y1c,&e2c,&x2c,&y2c);
        eps1 = (1+epsc)/2;
        eps2 = (1-epsc)/2;
        chisqc = 0;
        for( in=0; in<nh; in++ ) {
           ex = phit[in].amp;
           ixy = phit[in].ich;
           iy = ixy/EMCzSize;
           ix = ixy - iy*EMCzSize;
           dx1 = x1c - ix;
           dy1 = y1c - iy;
           dx2 = x2c - ix;
           dy2 = y2c - iy;
           a0 = e1c*EMCcell(dx1, dy1, e1c) + e2c*EMCcell(dx2, dy2, e2c);
           d = epar00*epar00 + e0*( epar1*a0/e0 + epar2*a0*a0/e0/e0 +epar3*a0*a0*a0/e0/e0/e0 ) + e0*sqrt(e0)*epar4*a0/e0*(1-a0/e0)*sin4a + e0*e0*epar0*epar0;
           chisqc += (a0-ex)*(a0-ex)/d;
        }
        if( chisqc >= chisq2 ) {
           if( iter > 0 ) {
                dchi = chisqc-chisq2;
                dchi0 = gr*step;
                step /= (2*sqrt(1+dchi/dchi0));
           }
           step /= 2;
        }
        else {
// Calculation of gradient
           grec = 0;
           grxc = 0;
           gryc = 0;
           for( in=0; in<nh; in++ ) {
                ex = phit[in].amp;
                ixy = phit[in].ich;
                iy = ixy/EMCzSize;
                ix = ixy - iy*EMCzSize;
                dx1 = x1c - ix;
                dy1 = y1c - iy;
                dx2 = x2c - ix;
                dy2 = y2c - iy;
                a1 = e1c*EMCcell(dx1,dy1,e1c);
                a2 = e2c*EMCcell(dx2,dy2,e2c);
                a0 = a1 + a2;
                d = epar00*epar00 + e0*( epar1*a0/e0 + epar2*a0*a0/e0/e0 +epar3*a0*a0*a0/e0/e0/e0 ) + e0*sqrt(e0)*epar4*a0/e0*(1-a0/e0)*sin4a + e0*e0*epar0*epar0;
                dd = (a0-ex)/d;
                dchida = dd*( 2 - dd*(epar1 + 2*epar2*a0/e0 + 3*epar3*a0*a0/e0/e0 + e0*sqrt(e0)*epar4*sin4a*(1-2*a0/e0) + 2*epar0*epar0*a0) );
                gx1 = ( e1c*EMCcell(x1c+0.05-ix,dy1,e1c) - a1 )*20;
                gx2 = ( e2c*EMCcell(x2c+0.05-ix,dy2,e2c) - a2 )*20;
                gy1 = ( e1c*EMCcell(dx1, y1c+0.05-iy,e1c) - a1 )*20;
                gy2 = ( e2c*EMCcell(dx2, y2c+0.05-iy,e2c) - a2 )*20;
                grec += (dchida*((a1/e1c-a2/e2c)*e0 - (gx1+gx2)*dxc -(gy1+gy2)*dyc)/2);
                grxc += (dchida*(gx1*eps2-gx2*eps1));
                gryc += (dchida*(gy1*eps2-gy2*eps1));
           }
           grc = sqrt(grec*grec + grxc*grxc + gryc*gryc);
           if( grc < 1e-10 ) grc = 1e-10;
           if( iter > 0 ) {
                cosi = (gre*grec + grx*grxc + gry*gryc ) / (gr*grc);
                scal = ABS(gr/grc - cosi);
                if( scal < 0.1 ) scal = 0.1;
                step /= scal;
           }
           chisq2 = chisqc;
           eps0 = epsc;
           dx0 = dxc;
           dy0 = dyc;
           gre = grec;
           grx = grxc;
           gry = gryc;
           gr = grc;
        }
        epsc = eps0 - step*gre/gr;
        while( ABS(epsc) >= epsmax ) {
           step /= 2;
           epsc = eps0 - step*gre/gr;
        }
        dxc = dx0 - step*grx/gr;
        dyc = dy0 - step*gry/gr;
        if( step*gr < stpmin ) break;
     }
//     if( (*pchi)*(nh-3)-chisq2 < delch ) return;
     if( (*pchi)*nh-chisq2 < delch ) return;
//     dof = nh - 5;
     dof = nh;
     if( dof < 1 ) dof = 1;
     *pchi = chisq2/dof;
     c3to5(e0,x0,y0,eps0,dx0,dy0,pe1,px1,py1,pe2,px2,py2);
}

float ClusterChisq(int nh, hit* phit, float e, float x, float y)
{
     float chi=0;
     int ixy, ix, iy;
     float et, a, d;

     for( int in=0; in<nh; in++ ) {
        ixy = phit[in].ich;
        iy = ixy/EMCzSize;
        ix = ixy - iy*EMCzSize;
        et = phit[in].amp;
        a = EMCcell(x-ix, y-iy, -1);
        d = epar00*epar00 + e*(epar1*a + epar2*a*a + epar3*a*a*a) + 
            e*sqrt(e)*epar4*a*(1-a)*sin4a + e*e*epar0*epar0;
        a *= e;
        chi += (et-a)*(et-a)/d;
//      cout << nh << "  "<< in << "  "<< a/e << "  "<< x << "  "<< y << '\n';
     }
     return chi;
}

void mom1(int nh, hit* phit, float* pe, float* px, float* py)
// First momentum calculation
{
     float a, xw, yw, e;
     int ix, iy;
     hit* p;

     *pe=0;
     *px=0;
     *py=0;
     if( nh <= 0 ) return;
     p=phit;
     xw=0;
     yw=0;
     e=0;
     for( int i=0; i<nh; i++ ) {
        a = p->amp;
        iy = p->ich / EMCzSize;
        ix = p->ich - iy*EMCzSize;
        e += a;
        xw += ix*a;
        yw += iy*a;
        p++;
     }
     *pe = e;
     if( e <= 0 ) return;
     *px = xw/e;
     *py = yw/e;
}

void momenta(int nh, hit* phit, float* pe, float* px, float* py, float* pxx, float* pyy, float* pyx )
// First and second momenta calculation
{
     float a, x, y, e, xx, yy, yx;
     int ix, iy, i;
     hit* p;

     *pe=0;
     *px=0;
     *py=0;
     *pxx=0;
     *pyy=0;
     *pyx=0;
     if( nh <= 0 ) return;

     p=phit;
     x=0;
     y=0;
     e=0;
     xx=0;
     yy=0;
     yx=0;
     for( i=0; i<nh; i++ ) {
        a = p->amp;
        iy = p->ich / EMCzSize;
        ix = p->ich - iy*EMCzSize;
        e += a;
        x += ix*a;
        y += iy*a;
        xx += a*ix*ix;
        yy += a*iy*iy;
        yx += a*ix*iy;
        p++;
     }
     *pe = e;
     if( e <= 0 ) return;
     x /= e;
     y /= e;
     xx = xx/e - x*x;
     yy = yy/e - y*y;
     yx = yx/e - y*x;

     *px = x;
     *py = y;
     *pxx = xx;
     *pyy = yy;
     *pyx = yx;
}

void c3to5(float e0, float x0, float y0, float eps, float dx, float dy, float* pe1, float* px1, float* py1,  float* pe2, float* px2, float* py2)
// 3 to 5-dim space conversion subroutine.
// eps=(e1-e2)/e0,  (e0=e1+e2), x0*e0=x1*e1+x2*e2, dx=x1-x2
{
     *pe1 = e0*(1+eps)/2;
     *pe2 = e0 - *pe1;
     *px1 = x0 + dx*(1-eps)/2;
     *py1 = y0 + dy*(1-eps)/2;
     *px2 = x0 - dx*(1+eps)/2;
     *py2 = y0 - dy*(1+eps)/2;
}

void SetProfileParameters( int sec, float Energy, float z, float y )
//
// Axis Z here means X in two dim Sector coord system !!! 
// So one has to supply (x,y) coordinates in cell units (x instead of z)
// If sec < 0 this routine changes only Energy dependent parameters - 
// the angle dependent ones are set in the previous call
//
{
     float t;
     static float sin2ax, sin2ay, sin2a, lgE;
     float vx, vy, vz;
     float xVert, yVert, zVert;
     int sign;

     if( sec >= 0 ) {
// Vertex coordinates in Sector coordinate system
        GlobalToSector( xVertex, yVertex, zVertex, &xVert, &yVert, &zVert );
        vx = xVert;
// It is if we use the coordinates of the Sector Center (not zero-tower)
//      vy = y - (EMCySize-1)*Tower_ySize/2 - yVert;
//      vz = z - (EMCzSize-1)*Tower_zSize/2 - zVert;
        vy = y*Tower_ySize - yVert;
        vz = z*Tower_zSize - zVert;
// From this point X coord in sector frame is Z coord in Global Coord System !!!
        sinax = vz/sqrt(vz*vz+vx*vx);
        sinay = vy/sqrt(vy*vy+vx*vx);
        if( ABS(sinax) > 0.5 || ABS(sinay) > 0.5 ) 
          printf("!!!!! mEmcClusterChi2: Something strange in SetProfileParameters: sinax=%f  sinay=%f\n",sinax, sinay);
//      printf("Vertx=(%f,%f,%f)\n",xVert,yVert,zVert);
//      printf("(Ax,Ay)=(%f,%f)\n", sinax, sinay);
        t = (vz*vz+vy*vy)/(vz*vz+vy*vy+vx*vx);
        sin2a=t;
        sin4a=t*t;
        sin2ax=sinax*sinax;
        sin2ay=sinay*sinay;
     }

     if( Energy <= 1.e-10 ) lgE=0;
     else lgE=log(Energy);

     ppar1=0.59-(1.45+0.13*lgE)*sin2a;
     ppar2=0.265+(0.80+0.32*lgE)*sin2a;
     ppar3=0.25+(0.45-0.036*lgE)*sin2a;
     ppar4=0.42;

     if( sinax > 0 ) sign = 1;
     else sign = -1;
     pShiftx = (1.05+0.12*lgE) * sin2ax * sign;
//      pShiftx=0;

     if( sinay > 0 ) sign = 1;
     else sign = -1;
     pShifty = (1.05+0.12*lgE) * sin2ay * sign;
//      pShifty=0;
}


float EMCcell( float xc, float yc, float en )
//
// Calculates the energy deposited in the tower, the distance between 
// its center and shower Center of Gravity being (xc,yc)
// en - shower energy
// If en<0 -> no Shower Profile parameters change is needed
//
{
     float dx, dy, r1, r2, r3, e;

     if( en > 0 ) SetProfileParameters(-1,en,xc,yc);
     dx=ABS(xc-pShiftx);
     dy=ABS(yc-pShifty);
     e=0;
     r2=dx*dx+dy*dy;
     r1=sqrt(r2);
     r3=r2*r1;
     e=ppar1*exp(-r3/ppar2)+ppar3*exp(-r1/ppar4);

     return e;
}

//********************************************************************
int Find_Clusters( list<hit> hlist, list<cluster>* pClList)
//********************************************************************
// Cluster search algorithm based on Lednev's one developed for GAMS.
// Returns -1 if MaxLen parameter is too small (just increase it)
{
        int nhit, nCl;
//      static int MaxLen=1000;
//      int* LenCl = new int[MaxLen];
        int LenCl[MaxLen];
        int next, ib, ie, ia, iab, iae, last, LastCl, leng, ich;
        hit *vv;
        hit *vhit, *vt;
        cluster Clt;
        list<hit>::iterator ph;
        list<hit> hl;

        (*pClList).erase(  (*pClList).begin(),  (*pClList).end() );
        nhit = hlist.size();

        if( nhit <= 0 ) return 0;
        if( nhit == 1 ) {
            Clt.ReInitialize( hlist );
            pClList->push_back(Clt);
            return 1;
        }

        vt = new hit[nhit];
        vhit = new hit[nhit];

        ph = hlist.begin();
        vv = vhit;
        while( ph != hlist.end() ) *vv++ = *ph++;

        qsort( vhit, nhit, sizeof(hit), Hit_NCompare );

        nCl=0;
        next = 0;
        for( ich=1; ich<nhit+1; ich++ ){
           if( ich<nhit ) ia=vhit[ich].ich;
           if( (ia-vhit[ich-1].ich > 1) || (ich >= nhit) || 
// Protect against gluing together edges of the different rows
              (ia-ia/EMCzSize*EMCzSize == 0) ){
              ib=next;
              ie=ich-1;
              next=ich;
              if( nCl >= MaxLen ) {
                return -1;
//              ResizeVector(LenCl, MaxLen, MaxLen*2 );
//              MaxLen *= 2;
              }
              nCl++;
              LenCl[nCl-1]=next-ib;
              if( nCl > 1 ) {
// Job to glue the subclusters
                iab=vhit[ib].ich;       // The last subcl begin
                iae=vhit[ie].ich;       // The last subcl end
                last=ib-1;              // The prelast subcl end
                LastCl=nCl-2;
                for( int iCl=LastCl; iCl>=0; iCl-- ){
                   leng=LenCl[iCl];
#ifdef GLUE_CLUSTER_CORNER
//*********** Begin: glue clusters with adjacent corner **********
//                 if( (iab-vhit[last].ich > EMCzSize+1) ) goto new_ich;
                   if( (iab-vhit[last].ich > EMCzSize+1) || 
// This is if iab-channel is on the left edge
                       ( (iab-vhit[last].ich == EMCzSize+1) &&
                         (iab-iab/EMCzSize*EMCzSize == 0) ) ) goto new_ich;
                   for( int ichc=last; ichc >= last-leng+1; ichc-- ) {
//                      if( iab-vhit[ichc].ich >  EMCzSize+1 ) goto new_icl;
                        if( (iab-vhit[ichc].ich >  EMCzSize+1) ||  
// This is if iab-channel is on the left edge
                            ( (iab-vhit[ichc].ich == EMCzSize+1) &&
                              (iab-iab/EMCzSize*EMCzSize == 0) ) ) goto new_icl;
//                      if( iae-vhit[ichc].ich >= EMCzSize-1 ) {
                        if( (iae-vhit[ichc].ich >= EMCzSize-1) && 
// This is if iae-channel is not on the right edge
                            ( (iae-vhit[ichc].ich != EMCzSize-1) || 
                              ((iae+1)-(iae+1)/EMCzSize*EMCzSize != 0) ) ) {
//*********** End: glue clusters with adjacent corner **********
#else
//*********** Begin: glue clusters with adjacent edge **********
                   if( iab-vhit[last].ich > EMCzSize ) goto new_ich;
                   for( int ichc=last; ichc >= last-leng+1; ichc-- ) {
                        if( iab-vhit[ichc].ich > EMCzSize ) goto new_icl;
                        if( iae-vhit[ichc].ich >= EMCzSize ) {
//*********** End: glue clusters with adjacent edge **********
#endif
//                        vt = new hit[leng];
                          CopyVector( &vhit[last+1-leng], vt, leng );
                          CopyVector( &vhit[last+1], &vhit[last+1-leng], ib-1-last );
                          CopyVector( vt, &vhit[ib-leng], leng );
                          for( int i=iCl; i<nCl-2; i++ ) LenCl[i]=LenCl[i+1];
                          ib -= leng;
                          LenCl[nCl-2]=LenCl[nCl-1]+leng;
                          nCl--;
                          goto new_icl;
                        }
                   }
new_icl:           last=last-leng;
                }
              }
           }
new_ich:   continue;
        }
//      cout << "NCluster=" << nCl << "\n";
//      for( ich=0; ich<nhit; ich++ ) cout << vhit[ich].ich <<"\n";
        if( nCl > 0 ) {
           ib=0;
           for( int iCl=0; iCl<nCl; iCl++ ) { 
                leng=LenCl[iCl];
                hl.erase( hl.begin(), hl.end() );
                for( ich=0; ich<leng; ich++ ) hl.push_back(vhit[ib+ich]);
                Clt.ReInitialize(hl);
                ib += LenCl[iCl];
                pClList->push_back(Clt);
           }
        }
//      delete LenCl;
        delete vhit;
        delete vt;

        return nCl;
}


void SetGeometry( SecGeom sg, float* Vertex )
//
// Sets the sector geometry and vertex position
// Should be called before Find_Clusters(...)
//
{

// The number of towers in Z and Y dir.

  EMCzSize = sg.nz;
  EMCySize = sg.ny;

// Tower size in Z and Y dir.

  Tower_zSize = sg.Tower_zSize;
  Tower_ySize = sg.Tower_ySize;

// Vertex coodinates (cm) in Global Coord System

   xVertex = Vertex[0];
   yVertex = Vertex[1];
   zVertex = Vertex[2];

// Zero-Tower center coodinates (cm) in Global Coord System

   //   xSector=500;
   //   ySector=0;
   //   zSector=0;
   xSector = sg.xyz0[0];
   ySector = sg.xyz0[1];
   zSector = sg.xyz0[2];

// Normal vector of the sector (tower?). For West Sectors it consides with Sector X' axis. We suppose, that Z global axis consides with Z' sector one

   //   normx=1;
   //   normy=0;
   //   normz=0;
   normx = sg.nxyz[0];
   normy = sg.nxyz[1];
   normz = sg.nxyz[2];

//      printf(" ************ %d   %d   %f   %f\n", EMCzSize, EMCySize, Tower_zSize, Tower_ySize);

}

void GlobalToSector( float xgl, float ygl, float zgl, float* px, float* py, float* pz )
{
    float x, y, z;

    x = xgl-xSector;
    y = ygl-ySector;
    z = zgl-zSector;
    *px =  normx*x + normy*y;
    *py = -normy*x + normx*y;
    *pz = z;
}

void SectorToGlobal( float xsec, float ysec, float zsec, float* px, float* py, float* pz ) 
{
    *px =  normx*xsec - normy*ysec + xSector;
    *py =  normy*xsec + normx*ysec + ySector;
    *pz = zsec + zSector;
}

void SectorToGlobalErr( float dxsec, float dysec, float dzsec, float* pdx, float* pdy, float* pdz ) 
{
    *pdx = sqrt( normx*normx*dxsec*dxsec + normy*normy*dysec*dysec );
    *pdy = sqrt( normy*normy*dxsec*dxsec + normx*normx*dysec*dysec );
    *pdz = dzsec;
}

void CalculateErrors( float e, float x, float y, float* pde, float* pdx, float* pdy, float* pdz)
//
// Returns the errors for the reconstructed energy and position 
// (in the hypothesis of EM shower)
// Should be called only just after CorrectPosition !!!
//
{
    float de, dy, dz, dxg, dyg, dzg;
    static float ae = 0.076, be = 0.022;        // de/E = a/sqrt(E)&b
    static float a = 0.57, b = 0.155, d = 1.6;  // dx = a/sqrt(E)+b (cm)
    // This is for hadrons
    //    static float a = 0.83, b = 0.95, d = 8.3;     // dx = a/sqrt(E)+b (cm)
    static float dx = 0.1;  // (cm)

//      printf("Pos=(%f,%f) a=(%f,%f)  \n", x, y, sinax, sinay);
    de = sqrt( ae*ae*e + be*be*e*e );
    dz = a/sqrt(e) + b;
    dy = dz;
    dz = sqrt( dz*dz + d*d*sinax*sinax );
    dy = sqrt( dy*dy + d*d*sinay*sinay );

    SectorToGlobalErr( dx, dy, dz, &dxg, &dyg, &dzg );

    *pde = de;
    *pdx = dxg;
    *pdy = dyg;
    *pdz = dzg;
}

void CorrectPosition( float Energy, float x, float y, float* pxc, float* pyc )
//
// Corrects the Shower Center of Gravity for the systematic error due to 
// the limited tower size and angle shift
//
// Everything here is in cm. 
// (x,y) - CG position, (*pxc,*pyc) - corrected position
//
{
   float xShift, yShift, xZero, yZero, bx, by;
   float t, x0, y0;
   int ix0, iy0;
   int signx, signy;

   SetProfileParameters( 0, Energy, x/Tower_zSize, y/Tower_ySize );
   if( sinax > 0 ) signx =  1;
   else            signx = -1;
   if( sinay > 0 ) signy =  1;
   else            signy = -1;
   t = 1.93+0.383*log(Energy);
   // This is for hadrons !!!
   //   t = 3.5;
   xShift = t*sinax;
   yShift = t*sinay;
   xZero=xShift-(0.417*ABS(sinax)+1.500*sinax*sinax)*signx;
   yZero=yShift-(0.417*ABS(sinay)+1.500*sinay*sinay)*signy;
   t = 0.98 + 0.98*sqrt(Energy);
   bx = 0.16 + t*sinax*sinax;
   by = 0.16 + t*sinay*sinay;
   // This is for hadrons !!!
   //   bx = 0.5;
   //   by = 0.5;
//      printf("Sin=(%f,%f) Shift=(%f,%f)  Zero=(%f,%f)\n",sinax,sinay,xShift,yShift,xZero,yZero);

   x0 = x/Tower_zSize;
   x0 = x0 - xShift + xZero;
   ix0 = lowint(x0 + 0.5);
   if( ABS(x0-ix0) <= 0.5 ) {
        x0 = (ix0-xZero)+bx*asinht( 2.*(x0-ix0)*sinht(0.5/bx) );
        *pxc = x0*Tower_zSize;
   }
   else {
        *pxc =  x - xShift*Tower_zSize;
        printf("????? mEmcClusterChi2: Something wrong in CorrectPosition: x=%f  dx=%f\n", x, x0-ix0);
   }

   y0 = y/Tower_ySize;
   y0 = y0 - yShift + yZero;
   iy0 = lowint(y0 + 0.5);
   if( ABS(y0-iy0) <= 0.5 ) {
        y0 = (iy0-yZero)+by*asinht( 2.*(y0-iy0)*sinht(0.5/by) );
        *pyc = y0*Tower_ySize;
   }
   else {
        *pyc = y - yShift*Tower_ySize;
        printf("????? mEmcClusterChi2: Something wrong in CorrectPosition: y=%f  dy=%f\n", y, y0-iy0);
   }

}



int Hit_NCompare( const void* h1, const void* h2 )
{
     return ( ((hit*)h1)->ich - ((hit*)h2)->ich );
}

int Hit_ACompare( const void* h1, const void* h2 )
{
     float amp1 = ((hit*)h1)->amp;
     float amp2 = ((hit*)h2)->amp;
     return (amp1<amp2) ? 1: (amp1>amp2) ? -1 : 0;
}

void ZeroVector( int* v, int N )
{
        int* p = v;
        for( int i=0; i < N; i++ ) *p++ = 0;
}

void ZeroVector( float* v, int N )
{
        float* p = v;
        for( int i=0; i < N; i++ ) *p++ = 0;
}

void ZeroVector( hit* v, int N )
{
        for( int i=0; i < N; i++ ) { v[i].ich=0; v[i].amp=0; v[i].tof=0; }
}

void ResizeVector( int* Vector, int OldSize, int NewSize )
{
        int* vsave;

        if( OldSize <= 0 ) { Vector = new int[NewSize]; return; }
        vsave = new int[OldSize];
        CopyVector( Vector, vsave, OldSize );
        delete Vector;
        Vector = new int[NewSize];
        if( NewSize > OldSize ) CopyVector( vsave, Vector, OldSize );
        else                    CopyVector( vsave, Vector, NewSize );
        delete vsave;
        return;
}

void CopyVector( int* from, int* to, int N )
{
        if( N <= 0 ) return;
        for( int i=0; i<N; i++ ) to[i]=from[i];
}

void CopyVector( hit* from, hit* to, int N )
{
        if( N <= 0 ) return;
        for( int i=0; i<N; i++ ) to[i]=from[i];
}

void SetChi2Limit( int limit )
//
// Sets the limit for PeakArea splitting onto 2 EMShowers:
// limit=0 -> For the Conf Level: 0%
// limit=1 -> For the Conf Level: 1% for GEANT, 2%-5% for TestBeam
// limit=2 -> For the Conf Level: 2% for GEANT, 4%-7% for TestBeam
//
{
  int i;

  switch ( limit ) {
  case 0:
    for( i=0; i<50; i++ ) Chi2Level[i]=9999.;
    break;
  case 1:
    for( i=0; i<50; i++ ) Chi2Level[i]=Chi2Level1[i];
    break;
  case 2:
    for( i=0; i<50; i++ ) Chi2Level[i]=Chi2Level2[i];
    break;
  default:
    for( i=0; i<50; i++ ) Chi2Level[i]=Chi2Level1[i];
    break;
  }
}


float Chi2Limit( int ND )
//  Here the reverse Chi2Correct function is used
{
  float rn, a, b, chi2;

  if( ND < 1 ) return 9999.;

  chi2 = Chi2Level[min(ND,50)-1];
  if( chi2 > 100 ) return 9999.;

  rn = ND;
  b = 0.072*sqrt(rn);
  a = 6.21/(sqrt(rn)+4.7);

  return chi2*a/(1.-chi2*b);
}


float Chi2Correct( float Chi2, int ND )
// Chi2 - is reduced Chi2: Chi2/ND !!
{
  float rn, a, b, c;

  if( ND < 1 ) return 9999.;

  rn = ND;
  b = 0.072*sqrt(rn);
  a = 6.21/(sqrt(rn)+4.7);
  c = a + b*Chi2;
  if( c < 1 ) c=1;

  return Chi2/c;
}

void SetThreshold( float Thresh )
{
  epar0 = Thresh*0.07;
  epar00 = max( Thresh/3, 0.005 );
}


/* Future improvements:

1. Find_Clusters: to ensure that all hits are above energy threshold 
set by SetThreshold routine (or default one)

*/