#include "ProtoCluster.h"
#include "EVENT/LCIO.h"
#include "EVENT/LCCollection.h"
#include "MyCaloHitExtended.h"
#include "MyTrackExtended.h"
// GEAR include files
#include <marlin/Global.h>
#include <gear/GEAR.h>
#include <gear/CalorimeterParameters.h>
#include <gear/LayerLayout.h>
// io/s stream
#include "streamlog/streamlog.h"
#include <sstream>
#include <iostream>

#include <cmath>

using namespace marlin ;

const float pi = acos(-1.0);
const float twopi = 2.0*pi;

ProtoCluster::ProtoCluster(){ 

} 

ProtoCluster::ProtoCluster(MyCaloHitExtended* seedHit){ 
 
  // protocluster constructed from a calorimeter hit
  // assume it is growing like a PHOTON 

  _temporaryReclusteringAttempt=false;
  _ID=0;

  this->initialise();
  this->AddHit(seedHit);

  // initially assume cluster originates from IP
  
  _mip =false;

  float x = seedHit->getRadialUnitVec()[0];
  float y = seedHit->getRadialUnitVec()[1];
  float z = seedHit->getRadialUnitVec()[2];
  float r = sqrt(x*x+y*y+z*z);
  _initialDir[0] = x/r;
  _initialDir[1] = y/r;
  _initialDir[2] = z/r;
  _uR[0] = _initialDir[0];
  _uR[1] = _initialDir[1];
  _uR[2] = _initialDir[2];


  _currentDirectionOK = false; 
  _currentDir[0] = x/r;
  _currentDir[1] = y/r;
  _currentDir[2] = z/r;

  _stillGrowing = true;

  // defaults
  _tanConeAngleECAL   = 0.36;
  _tanConeAngleHCAL   = 1.00;
  _additionalPadsECAL = 1.5;
  _additionalPadsHCAL = 2.5;
  _sameLayerPadCutECAL = 1.8;
  _sameLayerPadCutHCAL = 1.8;
  _maximumDistanceForConeAssociationECAL = 250.;
  _maximumDistanceForConeAssociationHCAL = 500.;
  _trackingWeight = 2;

} 

ProtoCluster::ProtoCluster(MyTrackExtended* trackAR){ 
 
  // protocluster constructed from a track
  // assume it is growing like a MIP

  //  float rhat[3];
  _ID=0;
  this->initialise();

  _temporaryReclusteringAttempt=false;
  _mip = true;
  _trackSeeded = true;
  _centroidX[0]=trackAR->getSeedPosition()[0];
  _centroidY[0]=trackAR->getSeedPosition()[1];
  _centroidZ[0]=trackAR->getSeedPosition()[2];
  float r = sqrt(_centroidX[0]*_centroidX[0]+_centroidY[0]*_centroidY[0]+_centroidZ[0]*_centroidZ[0]);
  _uR[0] = _centroidX[0]/r;
  _uR[1] = _centroidY[0]/r;
  _uR[2] = _centroidZ[0]/r;

  // store region of detector for speed optimisation 
  // didn't make much difference to performance
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // WARNING !!!!! THIS CODE IS NOT WELL PROTECTED AGAINST OVERWRITE
  //               IF RESURRECTED CHECK CAREFULLY !!!!
  //         !!!!!
  //if(1==2){
  //  float r =_centroidX[0]*_centroidX[0]+
  //    _centroidY[0]*_centroidY[0]+
  //    _centroidZ[0]*_centroidZ[0];
  //  r = sqrt(r);
  //  rhat[0]    = _centroidX[0]/r;
  //  rhat[1]    = _centroidY[0]/r;
  //  rhat[2]    = _centroidZ[0]/r;
  //  float phi = atan2(rhat[1],rhat[0]);
  //  if(phi<0)phi += twopi;
  //  float theta = acos(static_cast<double>(rhat[2]));
  //  float thetasector = theta/pi*(NUMBER_OF_THETA_SECTORS-2);
  //  float phisector   = phi/twopi*NUMBER_OF_PHI_SECTORS;
  //  int itheta = static_cast<int>(thetasector)+1;
  //  int iphi   = static_cast<int>(phisector);
  //  for(int i = iphi-1;i<=iphi+1;++i){
  //    for(int j = itheta-1;j<=itheta+1;++j){
  //	int ii = i;
  //	if(i>=NUMBER_OF_PHI_SECTORS)ii=0;
  //	if(i<0)ii=NUMBER_OF_PHI_SECTORS-1;
  //	int isectorOfInfluence   = j*NUMBER_OF_PHI_SECTORS+ii;
  //	_sectorCount[isectorOfInfluence]++;
  //    }
  //  }
  //}

  _initialDir[0] = trackAR->getSeedDirection()[0];
  _initialDir[1] = trackAR->getSeedDirection()[1];
  _initialDir[2] = trackAR->getSeedDirection()[2];

  _currentDirectionOK = false; 
  _currentDir[0] = _initialDir[0];
  _currentDir[1] = _initialDir[1];
  _currentDir[2] = _initialDir[2];

  _lastLayer = 0;
  _firstLayer = 0;
  _hitLayers = 0;

  // defaults
  _tanConeAngleECAL   = 0.36;
  _tanConeAngleHCAL   = 1.00;
  _additionalPadsECAL = 1.5;
  _additionalPadsHCAL = 2.5;
  _sameLayerPadCutECAL = 1.8;
  _sameLayerPadCutHCAL = 1.8;
  _maximumDistanceForConeAssociationECAL = 250.;
  _maximumDistanceForConeAssociationHCAL = 500.;
  
  _stillGrowing = true;
  
} 


float ProtoCluster::DistanceToTrack(MyTrackExtended* trackAR, int layers, bool useEndCap){ 
 
  float closest = 99999999.;
  float dmin = 99999999.;


  float x0 =trackAR->getSeedPosition()[0];
  float y0 =trackAR->getSeedPosition()[1];
  float z0 =trackAR->getSeedPosition()[2]; 
  float ux = trackAR->getSeedDirection()[0];
  float uy = trackAR->getSeedDirection()[1];
  float uz = trackAR->getSeedDirection()[2];

  if(useEndCap){
    x0 = trackAR->getEndcapSeedPosition()[0];
    y0 = trackAR->getEndcapSeedPosition()[1];
    z0 = trackAR->getEndcapSeedPosition()[2]; 
    ux = trackAR->getEndcapSeedDirection()[0];
    uy = trackAR->getEndcapSeedDirection()[1];
    uz = trackAR->getEndcapSeedDirection()[2];
  }

  int layersToCheck = layers;
  // if in barrel endcap overlap region search all layers
  if(useEndCap)layersToCheck = _lastLayer+1;
  for(int ilayer=_firstLayer; ilayer<  layersToCheck; ++ilayer){
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      float xi = hiti->getPosition()[0];
      float yi = hiti->getPosition()[1];
      float zi = hiti->getPosition()[2];
      float dx = xi-x0;
      float dy = yi-y0;
      float dz = zi-z0;
      float d = dx*dx+dy*dy+dz*dz;
      d = sqrt(d);
      if(d<dmin)dmin=d;
      float dAlong =  dx*ux + dy*uy+ dz*uz;
      if(dAlong>-100 && dAlong < 100){
	float ddx = uy*dz - uz*dy;
	float ddy = uz*dx - ux*dz;
	float ddz = ux*dy - uy*dx;
	float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
	float dPerp  = sqrt(distSq);  
	if(useEndCap){ 
	  if(hiti->getPhysicalLayer()>layers)dPerp=99999999.;
	}
	if(dPerp<closest)closest=dPerp;
      }
    }
  }


  // deal with the awkward region in barrel-endcap overlap where
  // pseudo geometry doesn't fully work 
  // Also project to HCAL through lumical
  // take care of recursive nature - only call once  
  if(!useEndCap){
    float z = fabs(trackAR->getSeedPosition()[2]);
    if(z>_zOfBarrel-10.&&z<_zOfEndcap){
      float eapproach = this->DistanceToTrack(trackAR,10,true);    
      if(eapproach<closest)closest = eapproach;
    }
    z = fabs(trackAR->getEndcapSeedPosition()[2]);
    // lumical region
    if(z>_zOfEndcap+50.){
      float eapproach = this->DistanceToTrack(trackAR,10,true);    
      if(eapproach<closest)closest = eapproach;
    }
  }

  return closest;
} 




float ProtoCluster::CentroidDistanceToTrack(MyTrackExtended* trackAR){ 
 
  float closest = 99999999.;

  float x0 =trackAR->getSeedPosition()[0];
  float y0 =trackAR->getSeedPosition()[1];
  float z0 =trackAR->getSeedPosition()[2]; 
  float ux = trackAR->getSeedDirection()[0];
  float uy = trackAR->getSeedDirection()[1];
  float uz = trackAR->getSeedDirection()[2];

  if(_se>0){
    float xi = _sxe/_se;
    float yi = _sye/_se;
    float zi = _sze/_se;
    float dx = xi-x0;
    float dy = yi-y0;
    float dz = zi-z0;
    float d = dx*dx+dy*dy+dz*dz;
    d = sqrt(d);
    float dAlong =  dx*ux + dy*uy+ dz*uz;
    float ddx = uy*dz - uz*dy;
    float ddy = uz*dx - ux*dz;
    float ddz = ux*dy - uy*dx;
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    float dPerp  = sqrt(distSq);  
    if(dPerp<closest)closest=dPerp;
  }

  return closest;
} 




float ProtoCluster::Centroid10DistanceToTrack(MyTrackExtended* trackAR){ 
 
  float closest = 99999999.;

  float x0 =trackAR->getSeedPosition()[0];
  float y0 =trackAR->getSeedPosition()[1];
  float z0 =trackAR->getSeedPosition()[2]; 
  float ux = trackAR->getSeedDirection()[0];
  float uy = trackAR->getSeedDirection()[1];
  float uz = trackAR->getSeedDirection()[2];

  if(_se>10){
    float xi = _sxe10/_se10;
    float yi = _sye10/_se10;
    float zi = _sze10/_se10;
    float dx = xi-x0;
    float dy = yi-y0;
    float dz = zi-z0;
    float d = dx*dx+dy*dy+dz*dz;
    d = sqrt(d);
    float dAlong =  dx*ux + dy*uy+ dz*uz;
    float ddx = uy*dz - uz*dy;
    float ddy = uz*dx - ux*dz;
    float ddz = ux*dy - uy*dx;
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    float dPerp  = sqrt(distSq);  
    if(dPerp<closest)closest=dPerp;
  }

  return closest;
} 




float ProtoCluster::Centroid20DistanceToTrack(MyTrackExtended* trackAR){ 
 
  float closest = 99999999.;

  float x0 =trackAR->getSeedPosition()[0];
  float y0 =trackAR->getSeedPosition()[1];
  float z0 =trackAR->getSeedPosition()[2]; 
  float ux = trackAR->getSeedDirection()[0];
  float uy = trackAR->getSeedDirection()[1];
  float uz = trackAR->getSeedDirection()[2];

  if(_se20>0){
    float xi = _sxe20/_se20;
    float yi = _sye20/_se20;
    float zi = _sze20/_se20;
    float dx = xi-x0;
    float dy = yi-y0;
    float dz = zi-z0;
    float d = dx*dx+dy*dy+dz*dz;
    d = sqrt(d);
    float dAlong =  dx*ux + dy*uy+ dz*uz;
    float ddx = uy*dz - uz*dy;
    float ddy = uz*dx - ux*dz;
    float ddz = ux*dy - uy*dx;
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    float dPerp  = sqrt(distSq);  
    if(dPerp<closest)closest=dPerp;
  }

  return closest;
} 



float ProtoCluster::DistanceToAltTrack(MyTrackExtended* trackAR, int layers){ 
 
  float closest = 99999999.;
  float dmin = 99999999.;


  float x0 =trackAR->getAltSeedPosition()[0];
  float y0 =trackAR->getAltSeedPosition()[1];
  float z0 =trackAR->getAltSeedPosition()[2]; 
  float ux = trackAR->getAltSeedDirection()[0];
  float uy = trackAR->getAltSeedDirection()[1];
  float uz = trackAR->getAltSeedDirection()[2];

  int layersToCheck = layers;
  // if in barrel endcap overlap region search all layers

  for(int ilayer=_firstLayer; ilayer<  layersToCheck; ++ilayer){
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      float xi = hiti->getPosition()[0];
      float yi = hiti->getPosition()[1];
      float zi = hiti->getPosition()[2];
      float dx = xi-x0;
      float dy = yi-y0;
      float dz = zi-z0;
      float d = dx*dx+dy*dy+dz*dz;
      d = sqrt(d);
      if(d<dmin)dmin=d;
      float dAlong =  dx*ux + dy*uy+ dz*uz;
      if(dAlong>-100 && dAlong < 100){
	float ddx = uy*dz - uz*dy;
	float ddy = uz*dx - ux*dz;
	float ddz = ux*dy - uy*dx;
	float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
	float dPerp  = sqrt(distSq);  
	if(dPerp<closest)closest=dPerp;
      }
    }
  }

  return closest;
} 


float ProtoCluster::DistanceToLine(float x0, float y0, float z0, float ux, float uy , float uz, int layers){ 
 
  float closest = 99999999.;
  float dmin = 99999999.;


  for(int ilayer=_firstLayer; ilayer<  layers; ++ilayer){
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      float xi = hiti->getPosition()[0];
      float yi = hiti->getPosition()[1];
      float zi = hiti->getPosition()[2];
      float dx = xi-x0;
      float dy = yi-y0;
      float dz = zi-z0;
      float d = dx*dx+dy*dy+dz*dz;
      d = sqrt(d);
      if(d<dmin)dmin=d;
      float dAlong =  dx*ux + dy*uy+ dz*uz;
      if(dAlong>-100){
	float ddx = uy*dz - uz*dy;
	float ddy = uz*dx - ux*dz;
	float ddz = ux*dy - uy*dx;
	float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
	float dPerp  = sqrt(distSq);  
	if(dPerp<closest)closest=dPerp;
      }
    }
  }
  
  return closest;
} 




void ProtoCluster::ConfigureGrowth(protoClusterDesign_t design){

  _tanConeAngleECAL    = design.tanConeAngleECAL;
  _tanConeAngleHCAL    = design.tanConeAngleHCAL;
  _additionalPadsECAL  = design.additionalPadsECAL;
  _additionalPadsHCAL  = design.additionalPadsHCAL;
  _sameLayerPadCutECAL = design.sameLayerPadCutECAL;
  _sameLayerPadCutHCAL = design.sameLayerPadCutHCAL;
  _maximumDistanceForConeAssociationECAL = design.maximumDistanceForConeAssociationECAL;
  _maximumDistanceForConeAssociationHCAL = design.maximumDistanceForConeAssociationHCAL;
  _trackingWeight                        = design.trackingWeight;
  _trackingCutOff                        = design.trackingCutOff;
  _trackingTube                          = design.trackingTube;


}


float ProtoCluster::DistanceFromInitialProjection(ProtoCluster* parent){

  float closest = 99999.;

  int fLayer = parent->FirstLayer();

  float x0 = parent->GetCentroid(fLayer)[0];
  float y0 = parent->GetCentroid(fLayer)[1];
  float z0 = parent->GetCentroid(fLayer)[2];

  float ux = parent->GetInitialDir()[0];
  float uy = parent->GetInitialDir()[1];
  float uz = parent->GetInitialDir()[2];

  float dmin = 99999.;
  for(int ilayer=_firstLayer; ilayer<=_firstLayer+4; ++ilayer){
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      float xi = hiti->getPosition()[0];
      float yi = hiti->getPosition()[1];
      float zi = hiti->getPosition()[2];
      float dx = xi-x0;
      float dy = yi-y0;
      float dz = zi-z0;
      float d = dx*dx+dy*dy+dz*dz;
      d = sqrt(d);
      if(d<dmin)dmin=d;
      float ddx = uy*dz - uz*dy;
      float ddy = uz*dx - ux*dz;
      float ddz = ux*dy - uy*dx;
      float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
      float dPerp  = sqrt(distSq);  

      float dAlong = fabs(ux*dx + uy*dy + uz*dz);

      if(dPerp/dAlong<0.9){
	if(dPerp<closest)closest=dPerp;
      }
    }
  }
  
  return closest;

}

void ProtoCluster::IsNowDefunct(){
  
  this->initialise();
  _defunct=true;

}

void ProtoCluster::initialise(){
			
  _cluster = NULL;
  _badness = -1;
  _electronID = false;
  _completeness = 0.;
  _purity = 0;
  _associatedTrack = false;
  _connected = false;
  _defunct = false;
  _hotHadron = false;
  _bestEnergy =0.;
  _energyEM=0.;
  _energyBarrel=0.;
  _energyEndcap=0.;
  _energyInMips=0.;
  _energyHAD=0.;
  _clusterClassification = PC_UNKNOWN;
  _earliestRelationID=0;
  _earliestRelation=this;
  _firstLayer = 999;
  _lastLayer = -999;
  _hitLayers = 0;
  _showerLayer = -999;
  _mipFitResult.ok = false;
  _fixedPhoton = false;
  _defunct = false;
  _mipFraction  = 0.;
  _clusterRms = 0.;
  _zOfEndcap = 100000.;
  _zOfBarrel = 100000.;
  _sx  = 0.;
  _sy  = 0.;
  _sz  = 0.;
  _se  = 0.;
  _sxe = 0.;
  _sye = 0.;
  _sze = 0.;

  _se10  = 0.;
  _sxe10 = 0.;
  _sye10 = 0.;
  _sze10 = 0.;

  _se20  = 0.;
  _sxe20 = 0.;
  _sye20 = 0.;
  _sze20 = 0.;

  _sl  = 0.;
  _sxl = 0.;
  _syl = 0.;
  _szl = 0.;
  _sll = 0.;
  _sxx = 0.;
  _sxy = 0.;
  _sxz = 0.;
  _syy = 0.;
  _syz = 0.;
  _szz = 0.;
  _nHits =0;
  _nMips =0;
  _trackSeeded = false;
  _photonID_dcoscut = 0.94;
  _fracCharged = 0.;
  _fracNeutral = 0.;
  _purity = 0.;
  _completeness = 0.;
  _photonLongProfileFraction=-999.;
  _photonLongShowerStart=-999.;


  for(int i=0;i<MAX_NUMBER_OF_LAYERS;++i){
    _sumX[i] = 0.;
    _sumY[i] = 0.;
    _sumZ[i] = 0.;
    _energy[i] = 0.;
    _dotUrUnorm[i] = -999.;
    _hitsByPseudoLayer[i].clear();
    _nOnTrackPath[i] = 0;
  }

  for(int i=0;i<NUMBER_OF_SECTORS;++i){
    _sectorCount[i] = 0;
  }

  _associatedTracks.clear();

  const gear::CalorimeterParameters& pEcalBarrel = Global::GEAR->getEcalBarrelParameters();
  const gear::LayerLayout& ecalBarrelLayout = pEcalBarrel.getLayerLayout();
  _nEcalLayers = ecalBarrelLayout.getNLayers();
  for(int i=0;i<ecalBarrelLayout.getNLayers();++i){
    _absThicknessByPhysicalLayer[i] = ecalBarrelLayout.getAbsorberThickness(i);
    if(i==0)_calThicknessByPhysicalLayer[i] = _absThicknessByPhysicalLayer[i];
    if(i!=0)_calThicknessByPhysicalLayer[i] = _calThicknessByPhysicalLayer[i-1]+_absThicknessByPhysicalLayer[i];
  }


}

void ProtoCluster::PhotonProfileID(bool truncate){
  

  float X0BinSize = 0.5;
  const unsigned int X0Bins = 100;
  float X0Max     = X0Bins*X0BinSize;
  float X0Profile[X0Bins];
  float expectedX0Profile[X0Bins];
  float eProfile[MAX_NUMBER_OF_LAYERS];
  for(unsigned int i=0;i<MAX_NUMBER_OF_LAYERS;i++)eProfile[i]=0.;
  for(unsigned int i=0;i<X0Bins;i++)X0Profile[i]=0.;
  for(unsigned int i=0;i<X0Bins;i++)expectedX0Profile[i]=0.;

  //
  float cost[MAX_NUMBER_OF_LAYERS];
  float sep[MAX_NUMBER_OF_LAYERS];

  for(unsigned int i=1;i<MAX_NUMBER_OF_LAYERS;++i){
    cost[i] = _dotUrUnorm[i];
  }
  for(unsigned int i=1;i<MAX_NUMBER_OF_LAYERS;++i)sep[i] = -999.;
  for(int i=_firstLayer; i<=_nEcalLayers; i++){
    if(_hitsByPseudoLayer[i].size()>0){    
      eProfile[i] = _energy[i];    
      int physLayer = _hitsByPseudoLayer[i][0]->getPhysicalLayer();
      sep[i]=_absThicknessByPhysicalLayer[physLayer];
    } 
  }

  int ilast  = 0;
  for(int i=_firstLayer; i<=_nEcalLayers; i++){
    if(sep[i]>0){
      if(ilast==0){
	for(int j=1;j<i;j++)sep[j]=sep[i];
	ilast = i;
      }
      if(i>ilast+1){
	for(int j=ilast+1;j<i;j++)sep[j]=sep[ilast];
      }
      ilast = i;
    }
  }

  if(ilast==0)sep[0]=_absThicknessByPhysicalLayer[1];
  for(int i=ilast+1; i<=_nEcalLayers; i++)sep[i]=sep[ilast];

  ilast  = 0;
  for(int i=_firstLayer; i<=_nEcalLayers; i++){
    if(cost[i]>0){
      if(ilast==0){
	for(int j=1;j<i;j++)cost[j]=cost[i];
	ilast = i;
      }
      if(i>ilast+1){
	for(int j=ilast+1;j<i;j++)cost[j]=cost[ilast];
      }
      ilast = i;
    }
  }
  if(ilast==0)cost[0]=1.;
  for(int i=ilast+1; i<=_nEcalLayers; i++)cost[i]=cost[ilast];

  float x0 = 3.5;
  float E0 = this->EnergyEM();
  // Effective critical energy
  float EC = 0.08;
  double a = 1.25+0.5*log(E0/EC);
  double lngammaa = lgamma(a);
  double gammaa = exp(lngammaa);

  float dist=0.;
  int   X0ProfileEnd=0;
  // step through ECAL making profile in radiation lengths
  float ecalE = 0.;
  int maxLayer = _nEcalLayers;
  if(truncate)maxLayer = 35;

  for(int i=1; i<=_nEcalLayers; i++){
    ecalE += eProfile[i];
    float dr = sep[i]/cost[i];
    float dt = dr/x0;
    float distStart = dist;
    dist+= dt;
    float deltaX0Bins = dt/X0BinSize;
    float bStart = distStart/X0BinSize;
    float bEnd   = dist/X0BinSize;
    if(bEnd>100)bEnd=100.; 
    unsigned int iStart = static_cast<int>(bStart);
    unsigned int iEnd = static_cast<int>(bEnd);
    float deltaStart = bStart-iStart;
    float deltaEnd   = 1.0-bEnd+iEnd;
    if(iStart<=100){
      for(unsigned int ibin=iStart;ibin<=iEnd;++ibin){
	float delta = 1.;
	if(ibin==iStart)delta = delta - deltaStart;
	if(ibin==iEnd)delta   = delta - deltaEnd;
	float frac = delta/deltaX0Bins;
	if(ibin<100)X0Profile[ibin] += eProfile[i]*frac; 
      }
    }
    X0ProfileEnd = iEnd;
    if(X0ProfileEnd>=100)X0ProfileEnd=100;
  }
  if(truncate)E0=ecalE;
 

  float t =0.;
  for(int i=0; i<100; i++){
    t+= X0BinSize;
    float de = E0/2.0*pow(t/2,a-1)*exp(-t/2)*X0BinSize/gammaa;
    expectedX0Profile[i] = de;
  }

  float diffEmin = 9999.;
  int   ioffmin  = 0;
  bool  stillgoing = true;
  for(int ioffset=0; stillgoing && ioffset<_nEcalLayers; ioffset++){
    float diffE = 0.;
    for(int i=0; i<X0ProfileEnd; i++){
      if(i-ioffset<0){
	diffE += X0Profile[i];
      }else{
	diffE += fabs(expectedX0Profile[i-ioffset]-X0Profile[i]);
      }
    }
    if(diffE<diffEmin){
      diffEmin = diffE;
      ioffmin = ioffset;
    }
    if(diffE-diffEmin>0.1)stillgoing = false;
  }
  
  if(ecalE>0){
    _photonLongProfileFraction =  diffEmin/ecalE;
    _photonLongShowerStart     =  ioffmin*X0BinSize ;
  }

}

void ProtoCluster::AddTrack(MyTrackExtended* track){

  _associatedTracks.push_back(track);

  _centroidX[0]=track->getSeedPosition()[0];
  _centroidY[0]=track->getSeedPosition()[1];
  _centroidZ[0]=track->getSeedPosition()[2];
  _initialDir[0] = track->getSeedDirection()[0];
  _initialDir[1] = track->getSeedDirection()[1];
  _initialDir[2] = track->getSeedDirection()[2];
  _associatedTrack=true;

}


void ProtoCluster::AddMCPFO(MyMCPFO* p){

  _associatedMCPFOs.push_back(p);

}


void ProtoCluster::ClearTracks(){

  _associatedTracks.clear();
  _associatedTrack = false;
  
}


MyTrackExtendedVec & ProtoCluster::GetTrackVec(){

  return _associatedTracks;
  
}


MyMCPFOVec & ProtoCluster::GetMCPFOVec(){

  return _associatedMCPFOs;
  
}



float ProtoCluster::genericDistanceToHit(const MyCaloHitExtended* candHit, const int searchLayer){ 
 
  // returns a floating point value indicating whether the candhit is part of this cluster
  // extrapolates forward from layer ilayer
  // the smaller the value, the better the association
  // values of greater than 1.0 indicate that the hit is not associated according to
  // the algorithm used by this cluster

  // Check hit is in roughly right region of detector
  //  int isector = candHit->getSector();
  //  if(_sectorCount[isector]==0)return 9999;

  // TrackingWeight == 3: is a very agressive tracking algorithm
  //   only hit on the track path are added below _trackingCutOfF
  //   deeper in the calorimeter - the tracked cluster will mop up any hits
  //   which are associated by the cone algorithm


  // check track projection if this is a track seeded cluster
  if(searchLayer==0&&_trackSeeded){
    float retVal = genericDistanceToHit(candHit, searchLayer, EXTERNAL_DIRECTION);
    return retVal;
  }

  // quit if no hits in this layer
  if(_hitsByPseudoLayer[searchLayer].size()==0)return 9999.;

  float dx = (_centroidX[searchLayer]-candHit->getPosition()[0]);
  float dy = (_centroidY[searchLayer]-candHit->getPosition()[1]);
  float dz = (_centroidZ[searchLayer]-candHit->getPosition()[2]);
  float dr = sqrt(dx*dx+dy*dy+dz*dz);

  float retVal=999.;
  if(dr>1000)return retVal;


  bool conePropagation = true;
  if(_trackingWeight>=3&&_trackSeeded&&searchLayer<=_trackingCutOff)conePropagation=false;
  if(_trackingWeight==4&&_trackSeeded){
    if(searchLayer<=_trackingCutOff&&_hitsByPseudoLayer[searchLayer].size()>0)return 9999.;
  }

  if(conePropagation){
    // if looking in same layer as hit do something different
    if(candHit->getPseudoLayer()==searchLayer)return this->genericDistanceInSameLayer(candHit);
    // calculate
    retVal = genericDistanceToHit(candHit, searchLayer, INITIAL_DIRECTION);
    if(_trackingWeight==3&&_trackSeeded&&searchLayer>_trackingCutOff)retVal=retVal/5.;
    if(_currentDirectionOK){
      float retValLocal   = genericDistanceToHit(candHit, searchLayer, CURRENT_DIRECTION);
      if(retValLocal<1.0&&_mip)retValLocal=retValLocal/5.0;
      if(retValLocal<retVal)retVal=retValLocal;
    }
  }

  // Match to extroplated track: if consistent give larger weight to association
  // Simple linear track extrapolation
  if(_trackSeeded&&_trackingWeight>0){
    float retValLocal = genericDistanceToTrackSeed(candHit, searchLayer);
    if(retValLocal<1.0)retValLocal=retValLocal/5.0;
    if(retValLocal<retVal)retVal=retValLocal;
  }
  
  return retVal;
  
}


float ProtoCluster::genericDistanceToTrackSeed(const MyCaloHitExtended* newhit, const int searchLayer){
  float retVal=99999.;
  float distance = 0.;
  // only extrapolate a certain distance into calorimeter
  if(_trackingWeight==1)distance = 100.;  
  if(_trackingWeight==2)distance = 250.;  
  if(_trackingWeight==3)distance = 1000.;  
  // check to see is a hit in previous few layers on track path
  bool check = false;

  // Old style tracking
  if(_trackingWeight==1){
  // check to see if there is a hit in current layer on track path
    for(unsigned int ihit = 0;!check && ihit<_hitsByPseudoLayer[searchLayer].size();ihit++){
      float approach = genericDistanceToTrackSeed(_hitsByPseudoLayer[searchLayer][ihit],distance);
      if(approach<1.)check = true;
    }
  }

  // New style tracking
  if(_trackingWeight==2){
    if(searchLayer<=2)check=true;
    for(int il = searchLayer-3; il<searchLayer;il++){
      if(_nOnTrackPath[il]>0)check=true;
    }
  }

  // If forced tracking
  if(_trackingWeight==3){
    if(searchLayer<=5)check=true;
    for(int il = searchLayer-3; il<searchLayer;il++){
      if(_nOnTrackPath[il]>0&&searchLayer<=_trackingCutOff)check=true;
    }
  }


  if(check)retVal = genericDistanceToTrackSeed(newhit,distance);

  return retVal;
}


float ProtoCluster::genericDistanceToTrackSeed(const MyCaloHitExtended* hit, const float distance){

  float _trackPathWidth=0.0;
  if(_trackingWeight==1)_trackPathWidth= 0.0;
  if(_trackingWeight==2)_trackPathWidth= 2.0;
  if(_trackingWeight==3)_trackPathWidth= 0.0;

  float retVal=99999.;
  float hitXYZ[3];

  hitXYZ[0] = hit->getPosition()[0];
  hitXYZ[1] = hit->getPosition()[1];
  hitXYZ[2] = hit->getPosition()[2];

  float dx   = (hitXYZ[0]-_centroidX[0]);
  float dy   = (hitXYZ[1]-_centroidY[0]);
  float dz   = (hitXYZ[2]-_centroidZ[0]);
  float dr   = dx*dx+dy*dy+dz*dz;
  dr = sqrt(dr);

  if(dr<distance){
    // Vector product
    float ddx = _initialDir[1]*dz - _initialDir[2]*dy;
    float ddy = _initialDir[2]*dx - _initialDir[0]*dz;
    float ddz = _initialDir[0]*dy - _initialDir[1]*dx;      
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    float dPerp  = sqrt(distSq);  
    // Scalar product
    //    float dAlong = _initialDir[0]*dx + _initialDir[1]*dy + _initialDir[2]*dz;
    //float absDAlong = fabs(dAlong);
    float slop = (_trackPathWidth*dr/distance+1.0);
    float dCut = slop*_additionalPadsECAL*hit->getHitSizeZ();
    if(hit->getCalorimeterHitType()==CALHITTYPE_HCAL)dCut = slop*_additionalPadsHCAL*hit->getHitSizeZ();
    if(_trackingWeight==3)dCut = 0.7*hit->getHitSizeZ();
    retVal = dPerp/dCut;
  }

  return retVal;
}




float ProtoCluster::genericDistanceToHit(const MyCaloHitExtended* candHit, const int searchLayer, const ProtoClusterDirectionType directionType){ 

  float retVal=999.;
  float dir[3];

  switch(directionType){
  case INITIAL_DIRECTION:
    for(int i=0;i<3;++i)dir[i]=_initialDir[i];
    break;
  case CURRENT_DIRECTION:
    for(int i=0;i<3;++i)dir[i]=_currentDir[i];
    break;
  case EXTERNAL_DIRECTION:
    for(int i=0;i<3;++i)dir[i]=_initialDir[i];
    float hitInCluster[3];
    hitInCluster[0] = _centroidX[0];
    hitInCluster[1] = _centroidY[0];
    hitInCluster[2] = _centroidZ[0];
    retVal = this->coneApproach(candHit,dir,hitInCluster);
    return retVal;
    break;
  default:
    break;
  }
  
 
  for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[searchLayer].size();ihit++){
    float hitInCluster[3];
    hitInCluster[0] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[0];
    hitInCluster[1] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[1];
    hitInCluster[2] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[2];
    float approach = this->coneApproach(candHit,dir,hitInCluster);
    if(approach<retVal)retVal = approach;
  }
    
  return retVal;
}

clusterContact_t ProtoCluster::ContactLayers(ProtoCluster* pC, float distance){ 

  clusterContact_t retval;
  int contactLayers = 0;
  float compared = 0;
  float contactFraction = 0.;
  int istart = pC->FirstLayer();
  if(_firstLayer>istart)istart=_firstLayer;
  int iend   = pC->LastLayer();
  if(_lastLayer<iend)iend=_lastLayer;
  float ri[3];
  float rj[3];
  //  std::cout << " Contact : " << istart << " - " << iend << std::endl;
  for(int ilayer = istart;ilayer<=iend;ilayer++){
    float layerMin = 99999.;
    bool layerDone = false;
    float cut = 0.;
    if(_hitsByPseudoLayer[ilayer].size()>0){
      float pixel = _hitsByPseudoLayer[ilayer][0]->getHitSizeZ();
      cut = 1.5*distance*pixel;
      cut = cut*cut;
      if(pC->hitsInLayer(ilayer)>0)compared+=1.;
    }
    for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[ilayer].size() && !layerDone;ihit++){
      ri[0] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[0];
      ri[1] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[1];
      ri[2] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[2];
      int nj = pC->hitsInLayer(ilayer);
      for(int jhit=0;jhit<nj;++jhit){
	const MyCaloHitExtended* hitj = pC->hitInLayer(ilayer,jhit);
	rj[0] = hitj->getPosition()[0];
	rj[1] = hitj->getPosition()[1];
	rj[2] = hitj->getPosition()[2];
	float dx = ri[0]-rj[0];
	float dy = ri[1]-rj[1];
	float dz = ri[2]-rj[2];
	float dr2 = dx*dx+dy*dy+dz*dz;
	//      std::cout << dr2 << std::endl;
	if(dr2<layerMin)layerMin=dr2;
      }
      if(layerMin<cut){
	layerDone = true;
	contactLayers++;
      }
    }
    //    std::cout << " Layer : " << ilayer << " min : " << layerMin << std::endl;
  }
    
  if(compared>0)contactFraction = contactLayers/compared;
  retval.contactLayers = contactLayers;
  retval.contactFraction = contactFraction;
  return retval;
}



float ProtoCluster::DistanceToHit(const MyCaloHitExtended* candHit){ 

  float ri[3];
  ri[0] = candHit->getPosition()[0];
  ri[1] = candHit->getPosition()[1];
  ri[2] = candHit->getPosition()[2];
  
  float dr2min = 1E10;  
  for(int ilayer = _firstLayer;ilayer<_lastLayer;ilayer++){
    for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[ilayer].size();ihit++){
      float rj[3];
      rj[0] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[0];
      rj[1] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[1];
      rj[2] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[2];
      float dx = ri[0]-rj[0];
      float dy = ri[1]-rj[1];
      float dz = ri[2]-rj[2];
      float dr2 = dx*dx+dy*dy+dz*dz;
      if(dr2<dr2min)dr2min=dr2;
    }
  }
    
  return sqrt(dr2min);
}






float ProtoCluster::DistanceToIsolatedHit(const MyCaloHitExtended* candHit){ 

  float ri[3];
  ri[0] = candHit->getPosition()[0];
  ri[1] = candHit->getPosition()[1];
  ri[2] = candHit->getPosition()[2];
  
  float dr2min = 1E10;  
  for(unsigned int ihit = 0;ihit<_isolatedHits.size();ihit++){
    float rj[3];
    rj[0] = (_isolatedHits[ihit])->getPosition()[0];
    rj[1] = (_isolatedHits[ihit])->getPosition()[1];
    rj[2] = (_isolatedHits[ihit])->getPosition()[2];
    float dx = ri[0]-rj[0];
    float dy = ri[1]-rj[1];
    float dz = ri[2]-rj[2];
    float dr2 = dx*dx+dy*dy+dz*dz;
    if(dr2<dr2min)dr2min=dr2;
    if(candHit==_isolatedHits[ihit])dr2min = -999.;
  }
    
  float drmin = -999.;
  if(dr2min>=0.)drmin=sqrt(dr2min);
  return drmin;
}




float ProtoCluster::DistanceToPoint(float x, float y, float z, int layers){ 

  float ri[3];
  ri[0] = x;
  ri[1] = y;
  ri[2] = z;
  
  float dr2min = 1E10;  
  int maxLayer = _firstLayer+layers;
  if(maxLayer>MAX_NUMBER_OF_LAYERS)maxLayer=MAX_NUMBER_OF_LAYERS;
  for(int ilayer = _firstLayer;ilayer<maxLayer;ilayer++){
    for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[ilayer].size();ihit++){
      float rj[3];
      rj[0] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[0];
      rj[1] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[1];
      rj[2] = (_hitsByPseudoLayer[ilayer][ihit])->getPosition()[2];
      float dx = ri[0]-rj[0];
      float dy = ri[1]-rj[1];
      float dz = ri[2]-rj[2];
      float dr2 = dx*dx+dy*dy+dz*dz;
      if(dr2<dr2min)dr2min=dr2;
    }
  }
  if(dr2min<0.0)dr2min=0.0;
  return sqrt(dr2min);
}





float ProtoCluster::BarrelEncapEnergySplit(){

  float retval = 1.0;
  if(_energyBarrel+_energyEndcap>0.){
    retval = _energyBarrel/(_energyBarrel+_energyEndcap);
    if(retval<0.5)retval=1.0-retval;
  }
  return retval;
}


float ProtoCluster::genericDistanceInSameLayer(const MyCaloHitExtended* candHit){

  float retVal =999.;
  int searchLayer = candHit->getPseudoLayer();   
  if(_hitsByPseudoLayer[searchLayer].size()==0)return 9999.;
  float hitXYZ[3];

  hitXYZ[0] = candHit->getPosition()[0];
  hitXYZ[1] = candHit->getPosition()[1];
  hitXYZ[2] = candHit->getPosition()[2];
  
  for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[searchLayer].size();ihit++){
    float hitInCluster[3];
    hitInCluster[0] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[0];
    hitInCluster[1] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[1];
    hitInCluster[2] = (_hitsByPseudoLayer[searchLayer][ihit])->getPosition()[2];
    float dx   = (hitXYZ[0]-hitInCluster[0]);
    float dy   = (hitXYZ[1]-hitInCluster[1]);
    float dz   = (hitXYZ[2]-hitInCluster[2]);
    float dr   = dx*dx+dy*dy+dz*dz;
    dr = sqrt(dr);
    float dCut = _sameLayerPadCutECAL*candHit->getHitSizeZ();
    if(candHit->getCalorimeterHitType()==CALHITTYPE_HCAL)dCut = _sameLayerPadCutHCAL*candHit->getHitSizeZ();
    float approach = dr/dCut;
    if(approach<retVal)retVal = approach;
  }

  return retVal;
}
 

float ProtoCluster::coneApproach(const MyCaloHitExtended* candHit, const float dir[3], const float position[3]){

  float retVal =999.;
  float hitXYZ[3];
  hitXYZ[0] = candHit->getPosition()[0];
  hitXYZ[1] = candHit->getPosition()[1];
  hitXYZ[2] = candHit->getPosition()[2];

  float dx   = (hitXYZ[0]-position[0]);
  float dy   = (hitXYZ[1]-position[1]);
  float dz   = (hitXYZ[2]-position[2]);
  float dr   = dx*dx+dy*dy+dz*dz;
  dr = sqrt(dr);
  float drcut = _maximumDistanceForConeAssociationECAL;
  if(candHit->getCalorimeterHitType()==CALHITTYPE_HCAL)drcut =_maximumDistanceForConeAssociationHCAL;
  
  if(dr<1000.){
    // Vector product
    float ddx = dir[1]*dz - dir[2]*dy;
    float ddy = dir[2]*dx - dir[0]*dz;
    float ddz = dir[0]*dy - dir[1]*dx;      
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    float dPerp  = sqrt(distSq);  
    // Scalar product
    float dAlong = dir[0]*dx + dir[1]*dy + dir[2]*dz;
    float absDAlong = fabs(dAlong);
    //    float cosTheta = dAlong/dr;

    float dCut = absDAlong*_tanConeAngleECAL+_additionalPadsECAL*candHit->getHitSizeZ();
    if(candHit->getCalorimeterHitType()==CALHITTYPE_HCAL)dCut = absDAlong*_tanConeAngleHCAL+_additionalPadsHCAL*candHit->getHitSizeZ();
    float approach = dPerp/dCut;
    if(dAlong<200. && dAlong > -10.){
      if(approach<retVal)retVal = approach;
    }
  }

  return retVal;
}



//void ProtoCluster::AddHit(MyCaloHitExtended* newhit){ 
//  this->AddHit(newhit,1.);
//}

void ProtoCluster::AddHit(MyCaloHitExtended* newhit, float weight){ 

  // if it is an isolated hit add it as isolated
  if(newhit->isIsolatedHit()){
    //stream_log(WARNING4) << " USING AddHit for an isolated hit " << std::endl;
    this->AddIsolatedHit(newhit);
    return;
  }

  _nHits++;
  _hits.push_back(newhit);

  //int isector = newhit->getSector();
  //if((_nHits>1 || _trackSeeded) && _sectorCount[isector]==0){
    //    std::cout << " Sector Count[" << isector << "] : " << _sectorCount[isector] << std::endl;
    //std::cout << " Hits : " << _nHits << " " << _trackSeeded << std::endl;
    //float x = newhit->getPosition()[0];
    //float y = newhit->getPosition()[1];
    //float z = newhit->getPosition()[2];
    //std::cout << " HIT x,y,z : " << x << "," << y << "," << z << std::endl;
    //std::cout << " TRK x,y,z : " << _centroidX[0] << "," << _centroidY[0] << "," << _centroidZ[0] << std::endl;
  //  }
  
  //for(int i = 0;i<9;++i){
  // int isectorOfInfluence = newhit->getSectorOfInfluence(i);
  // _sectorCount[isectorOfInfluence]++;
  //}

  int ilayer = newhit->getPseudoLayer();
  if(ilayer<_firstLayer)_firstLayer=ilayer;
  if(ilayer>_lastLayer)_lastLayer = ilayer;
  if(_hitsByPseudoLayer[ilayer].size()==0){
    _hitLayers++;
    _sumX[ilayer] = 0;
    _sumY[ilayer] = 0;
    _sumZ[ilayer] = 0;
  }


  //  for(unsigned int ihit = 0;ihit<_hitsByPseudoLayer[ilayer].size();ihit++){
  //   if(newhit==_hitsByPseudoLayer[ilayer][ihit])std::cout << " REPEATED HIT " << std::endl;
  // }
  _hitsByPseudoLayer[ilayer].push_back(newhit);

  // try and determine normal to layer from data (keep things geo indep)
  if(_hitsByPseudoLayer[ilayer].size()==3){
    MyCaloHitExtended* hit0 = _hitsByPseudoLayer[ilayer][0];
    MyCaloHitExtended* hit1 = _hitsByPseudoLayer[ilayer][1];
    MyCaloHitExtended* hit2 = _hitsByPseudoLayer[ilayer][2];
    float x0 = hit0->getPosition()[0];
    float y0 = hit0->getPosition()[1];
    float z0 = hit0->getPosition()[2];
    float dx1 = (hit1->getPosition()[0] -x0);
    float dy1 = (hit1->getPosition()[1] -y0);
    float dz1 = (hit1->getPosition()[2] -z0);
    float dx2 = (hit2->getPosition()[0] -x0);
    float dy2 = (hit2->getPosition()[1] -y0);
    float dz2 = (hit2->getPosition()[2] -z0);
    // use cross-product of two vectors in plane to get normal
    float xn = dy1*dz2-dz1*dy2;
    float yn = dz1*dx2-dx1*dz2;
    float zn = dx1*dy2-dy1*dx2;
    float rn = sqrt(xn*xn+yn*yn+zn*zn);
    if(rn>10.){
      xn = xn/rn;
      yn = yn/rn;
      zn = zn/rn;
      float cost = fabs(xn*_uR[0]+yn*_uR[1]+zn*_uR[2]);
      if(cost>0.3)_dotUrUnorm[ilayer] = cost;
    }
  }

  if(newhit->isCandidateMip()){
    _nMips++;
    _mipHitsByPseudoLayer[ilayer].push_back(newhit);
  }

  float x = newhit->getPosition()[0];
  float y = newhit->getPosition()[1];
  float z = newhit->getPosition()[2];
  float l = ilayer;

  float e = newhit->getEnergyEM();
  _sxe  += x*e;
  _sye  += y*e;
  _sze  += z*e;
  _se   += e;
  if(ilayer>0&&ilayer<=10){
    _sxe10  += x*e;
    _sye10  += y*e;
    _sze10  += z*e;
    _se10   += e;
  }
  if(ilayer>10&&ilayer<=20){
    _sxe20  += x*e;
    _sye20  += y*e;
    _sze20  += z*e;
    _se20   += e;
  }

  _sx  += x;
  _sy  += y;
  _sz  += z;
  _sl  += l;
  _sxl += x*l;
  _syl += y*l;
  _szl += z*l;
  _sll += l*l;
  _sxx += x*x;
  _sxy += x*y;
  _sxz += x*z;
  _syy += y*y;
  _syz += y*z;
  _szz += z*z;

  _sumX[ilayer] += x;
  _sumY[ilayer] += y;
  _sumZ[ilayer] += z;

  float Eem = newhit->getEnergyEM();
  float Ehad = newhit->getEnergyHAD();

  _bestEnergy += newhit->getEnergy()*weight;
  _energyEM += Eem*weight;
  _energyHAD += Ehad*weight;

  if(fabs(z)>_zOfEndcap-10){
    _energyEndcap+=Eem*weight;
  }else{
    _energyBarrel+=Eem*weight;
  }

  if(_electronID)_energyHAD= _energyHAD + (Eem-Ehad)*weight;
  _energyInMips += newhit->getEnergyInMips()*weight;

  // if(newhit->getEnergyHAD()>0.5){
  //  std::cout << " HOT HIT " << newhit->getEnergyHAD() << std::endl;
  // }
  
  _energy[ilayer] = _energy[ilayer] + Eem*weight;

  float ninlayer = static_cast<float>(_hitsByPseudoLayer[ilayer].size());

  _centroidX[ilayer] = _sumX[ilayer]/ninlayer;
  _centroidY[ilayer] = _sumY[ilayer]/ninlayer;
  _centroidZ[ilayer] = _sumZ[ilayer]/ninlayer;

  // check to see if hit lies on track extroplation
  if(_trackSeeded){
    float retValLocal = genericDistanceToTrackSeed(newhit, ilayer);
    if(retValLocal<1.0)_nOnTrackPath[ilayer]++;
  }


} 






void ProtoCluster::AddIsolatedHit(MyCaloHitExtended* newhit){ 


  for(unsigned int ihit = 0;ihit<_isolatedHits.size();ihit++){
    if(newhit==_isolatedHits[ihit])std::cout << " REPEATED ISOLATED HIT " << std::endl;
  }


  _isolatedHits.push_back(newhit);
  
  float Eem = newhit->getEnergyEM();
  float Ehad = newhit->getEnergyHAD();
  
  _bestEnergy += newhit->getEnergy();
  _energyEM += Eem;
  int ilayer = newhit->getPseudoLayer();
  _isolatedHitsByPseudoLayer[ilayer].push_back(newhit);

  _energy[ilayer] = _energy[ilayer] + Eem;
  _energyHAD += Ehad;
  _energyInMips+= newhit->getEnergyInMips();

  

} 



ProtoCluster::~ProtoCluster(){ 

  for(int i = _firstLayer;i<=_lastLayer;i++){
    _hitsByPseudoLayer[i].clear();
    _isolatedHitsByPseudoLayer[i].clear();
    _mipHitsByPseudoLayer[i].clear();
  }
  _daughters.clear();
  _parents.clear();
  _hits.clear();
  _isolatedHits.clear();
}

void ProtoCluster::Dump(){ 

  std::cout << " ***** PC : hits :  " << _nHits << ":" << _bestEnergy << std::endl;
  for(int i = _firstLayer;i<=_lastLayer;i++){
    std::cout << i << " : " << _hitsByPseudoLayer[i].size() << "{" << 
      _mipHitsByPseudoLayer[i].size() << "}" << std::endl; 
  }

  
  return;
}


void ProtoCluster::Update(int ilayer){
  
  // new hits may have been added so need to 
  //               recalculate local direction
  //               fragment track
  //               change shower <-> track  


  // redefine local shower direction based on first plane and last planes
  _mipFraction = this->MipFraction();
 
  if(_lastLayer-_firstLayer>10){
    float dx = _centroidX[_lastLayer] - _centroidX[_firstLayer];
    float dy = _centroidY[_lastLayer] - _centroidY[_firstLayer];
    float dz = _centroidZ[_lastLayer] - _centroidZ[_firstLayer];
    float dirMag = dx*dx+dy*dy+dz*dz;
    dirMag = sqrt(dirMag);
    _currentDir[0] = dx/dirMag;
    _currentDir[1] = dy/dirMag;
    _currentDir[2] = dz/dirMag;
    _currentDirectionOK = true; 
  }
  
  if(ilayer-_lastLayer>3)_stillGrowing = false;


  // 
  if(_lastLayer-_firstLayer>6){
    std::vector<vec4>points;
    int istart = _lastLayer-8; 
    if(_mipFraction<0.5)istart = _lastLayer-16;
    if(istart<_firstLayer)istart = _firstLayer;
    for(int i = istart; i <= _lastLayer;++i){
      if(_hitsByPseudoLayer[i].size()>0){
	vec4 point;
	point.x = _centroidX[i];
	point.y = _centroidY[i];
	point.z = _centroidZ[i];
	point.ilayer = i - _firstLayer;
	point.padSize = _hitsByPseudoLayer[i][0]->getHitSizeZ();
	points.push_back(point);
      }
    }
    fitResult fit = this->FitPoints(points);
    if(fit.ok){
      _currentDir[0] = fit.dir[0];  // valgrind claims uninit. var
      _currentDir[1] = fit.dir[1];  // valgrind claims uninit. var
      _currentDir[2] = fit.dir[2];  // valgrind claims uninit. var
      _currentDirectionOK = true; 
      float dot =  _currentDir[0]*_initialDir[0] +
	           _currentDir[1]*_initialDir[1] +
	           _currentDir[2]*_initialDir[2];
      // allow direction to deviate from initial direction but
      // if the difference is large only use if rms is small
      if(dot<0.75 && fit.chi2 > 5.0)_currentDirectionOK = false;
      if(dot<0.5 && fit.chi2 > 2.5)_currentDirectionOK = false;
      if(fit.chi2>2.5)_mip =false;
    }else{
      _currentDirectionOK =false;
    }
  }
}




int ProtoCluster::hitsInLayer(int ilayer){
  return _hitsByPseudoLayer[ilayer].size();
}


int ProtoCluster::isolatedHitsInLayer(int ilayer){
  return _isolatedHitsByPseudoLayer[ilayer].size();
}


MyCaloHitExtended* ProtoCluster::hitInLayer(int ilayer, int ihit){
  return _hitsByPseudoLayer[ilayer][ihit];
}


MyCaloHitExtended* ProtoCluster::isolatedHitInLayer(int ilayer, int ihit){
  return _isolatedHitsByPseudoLayer[ilayer][ihit];
}


void ProtoCluster::RemoveHit(MyCaloHitExtended* hit){ 

  std::vector<MyCaloHitExtended*>::iterator literH;
  bool notFound;
  
  int layer = hit->getPseudoLayer();
  notFound = true;
  literH = _hits.begin();
  while (notFound && literH != _hits.end()){
    if(*literH==hit){
      _hits.erase(literH);
      notFound = false;
    }
    literH++;
  }
  if(notFound)std::cout << "ProtoCluster::RemoveHit failed to remove hit from _hits" << layer << std::endl;


  notFound = true;
  literH = _hitsByPseudoLayer[layer].begin();
  while (notFound && literH != _hitsByPseudoLayer[layer].end()){
    if(*literH==hit){
      _hitsByPseudoLayer[layer].erase(literH);
      notFound = false;
    }
    literH++;
  }
  if(notFound)std::cout << "ProtoCluster::RemoveHit failed to remove hit from _hitsByPseudoLayer" << layer << std::endl;



  if(hit->isCandidateMip()){
    _nMips--;
    notFound = true;
    literH = _mipHitsByPseudoLayer[layer].begin();
    while (notFound && literH != _mipHitsByPseudoLayer[layer].end()){
      if(*literH==hit){
	_mipHitsByPseudoLayer[layer].erase(literH);
	notFound = false;
      }
      literH++;
    }
    if(notFound)std::cout << "ProtoCluster::RemoveHit failed to remove hit from _mipHitsByPseudoLayer" << std::endl;
  }

  if(_hitsByPseudoLayer[layer].size()==0){
    _sumX[layer] = 0;
    _sumY[layer] = 0;
    _sumZ[layer] = 0;
    if(layer==_lastLayer){
      for(int il = _firstLayer; il<_lastLayer;++il){
	if(_hitsByPseudoLayer[il].size()>0)_lastLayer=il;
      }
    }
  }else{
    float x = hit->getPosition()[0];
    float y = hit->getPosition()[1];
    float z = hit->getPosition()[2];
    float e = hit->getEnergyEM();
    float l = layer;
    _sx  -= x;
    _sy  -= y;
    _sz  -= z;
    _se  -= e;
    _sxe -= x*e;
    _sye -= y*e;
    _sze -= z*e;
    if(layer>0&&layer<=10){
      _se10  -= e;
      _sxe10 -= x*e;
      _sye10 -= y*e;
      _sze10 -= z*e;
    }
   if(layer>10&&layer<=20){
     _se20  -= e;
     _sxe20 -= x*e;
     _sye20 -= y*e;
     _sze20 -= z*e;
   }
    _sl  -= l;
    _sxl -= x*l;
    _syl -= y*l;
    _szl -= z*l;
    _sll -= l*l;
    _sxx -= x*x;
    _sxy -= x*y;
    _sxz -= x*z;
    _syy -= y*y;
    _syz -= y*z;
    _szz -= z*z; 
    _sumX[layer] -= x;
    _sumY[layer] -= y;
    _sumZ[layer] -= z;
    float ninlayer = static_cast<float>(_hitsByPseudoLayer[layer].size());
    _centroidX[layer] = _sumX[layer]/ninlayer;
    _centroidY[layer] = _sumY[layer]/ninlayer;
    _centroidZ[layer] = _sumZ[layer]/ninlayer;
  }

  
  _bestEnergy -= hit->getEnergy();
  _energyEM -= hit->getEnergyEM();
  float z = hit->getPosition()[2];
  if(fabs(z)>_zOfEndcap-10.){
    _energyEndcap -= hit->getEnergyEM();
  }else{
    _energyBarrel -= hit->getEnergyEM();
  }

  _energyHAD -= hit->getEnergyHAD(); 
  _energyInMips -= hit->getEnergyInMips(); 
  _energy[layer] = _energy[layer] -  hit->getEnergyEM();

   
  _nHits--;
  
  return;
} 

float* ProtoCluster::GetCentroid(int ilayer){

  _centroid[0] = _centroidX[ilayer];
  _centroid[1] = _centroidY[ilayer];
  _centroid[2] = _centroidZ[ilayer];

  return _centroid;

}


void ProtoCluster::TransverseProfile(std::vector<protoClusterPeaks_t> &peaks, int maxLayers){

  peaks.clear(); // valgrind claims invalid read and write somewhere in this method (stl)!
  int nbins   = 41;
  int ioffset = nbins/2; 
  int ifirst = this->FirstLayer();
  float x0   = this->GetCentroid(ifirst)[0];
  float y0   = this->GetCentroid(ifirst)[1];
  float z0   = this->GetCentroid(ifirst)[2];
  float r0   = sqrt(x0*x0+y0*y0+z0*z0);
  float ux   = x0/r0;
  float uy   = y0/r0;
  float uz   = z0/r0;
  float utx  = uy/sqrt(ux*ux+uy*uy);
  float uty  = -ux/sqrt(ux*ux+uy*uy);
  float utz  = 0;
  float vtx=0;
  float vty=0;
  float vtz=0;
  if(utx!=0&&uz!=0){
    vtx  = -uty/utx;
    vtz  = (uty*ux-uy*utx)/utx/uz;
    vty  = 1;
    float v    = sqrt(vtx*vtx+vty*vty+vtz*vtz);
    vtx = vtx/v;
    vty = vty/v;
    vtz = vtz/v;
  }
  if(utx==0&&uz!=0){
    vtx  = 1;
    vtz  = ux/uz;
    vty  = 0;
    float v    = sqrt(vtx*vtx+vty*vty+vtz*vtz);
    vtx = vtx/v;
    vty = vty/v;
    vtz = vtz/v;
  }
  if(utx==0&&uz==0){
    vtx  = 0;
    vtz  = 1;
    vty  = 0;
  }

  int  done[nbins][nbins];
  bool assigned[nbins][nbins];
  float tprofile[nbins][nbins];
  float tlprofile[nbins][nbins][MAX_NUMBER_OF_LAYERS];
  for(int i=0; i<nbins; ++i){
    for(int j=0; j<nbins; ++j){
      tprofile[i][j]=0;
      done[i][j]=false;
      assigned[i][j]=false;
      for(int k=0;k<=maxLayers ; k++)tlprofile[i][j][k]=0.;
    }
  }

  bool first = true;
  float pixelsize = 10.; 
  for(int i=0; i<maxLayers; i++){
    for(int ihit =0; ihit < this->hitsInLayer(i);++ihit){
      MyCaloHitExtended* hiti = this->hitInLayer(i,ihit);
      if(first){
	pixelsize  = hiti->getHitSizeZ();
      }
      float dx = (hiti->getPosition()[0] -x0)/pixelsize;
      float dy = (hiti->getPosition()[1] -y0)/pixelsize;
      float dz = (hiti->getPosition()[2] -z0)/pixelsize;
      float dut = dx*utx+dy*uty+dz*utz;
      float dvt = dx*vtx+dy*vty+dz*vtz;
      int idut  = static_cast<int>(dut+0.5+ioffset);
      int idvt  = static_cast<int>(dvt+0.5+ioffset);
      if(idut>=0&&idut<nbins&&idvt>=0&&idvt<nbins){
	tprofile[idut][idvt] += hiti->getEnergyEM();  
	tlprofile[idut][idvt][i] += hiti->getEnergyEM();  
      }
    }
  }

  // mask low ph region
  float threshold = 0.025;
  for(int i=0; i<nbins; i++){
    for(int j=0; j<nbins; j++){
      if(tprofile[i][j]<threshold)assigned[i][j]=true;
    }
  }
  // Search for peaks in profile
  // 
  bool stillfindingpeaks=true;
  int  npeaks=0;
  float dmin=9999.;

  while(stillfindingpeaks){
    float peakheight = 0.;
    int   ipeak =0;
    int   jpeak =0;
    float peakenergy=0;
    float xbar=0;
    float ybar=0;
    float xxbar=0;
    float yybar=0;
    float longitudinalProfile[MAX_NUMBER_OF_LAYERS];
    for(int i=0;i<=maxLayers;i++)longitudinalProfile[i]=0.;
    
    for(int i=1; i<nbins-1; i++){
      for(int j=1; j<nbins-1; j++){
	if(!assigned[i][j]){
	  if(tprofile[i][j]>peakheight){
	    peakheight = tprofile[i][j];
	    ipeak = i;
	    jpeak = j;
	  }
	}
      }
    }
    
    if(peakheight<threshold){
      stillfindingpeaks=false;
    }else{
      int point[1000][2];
      int pstart = 0;
      int pend   = 0;
      int pcurrent = pend;
      point[0][0] = ipeak;
      point[0][1] = jpeak;
      peakenergy=tprofile[ipeak][jpeak];
      for(int i=0;i<=maxLayers;++i)longitudinalProfile[i]+=tlprofile[ipeak][jpeak][i];
      xbar=(ipeak-ioffset)*tprofile[ipeak][jpeak];
      ybar=(jpeak-ioffset)*tprofile[ipeak][jpeak];
      xxbar=(ipeak-ioffset)*(ipeak-ioffset)*tprofile[ipeak][jpeak];
      yybar=(jpeak-ioffset)*(jpeak-ioffset)*tprofile[ipeak][jpeak];
      assigned[ipeak][jpeak]=true;
      dmin = sqrt(   (ipeak-ioffset)*(ipeak-ioffset)+
                     (jpeak-ioffset)*(jpeak-ioffset));
      npeaks++;
      float stillgoing = true;
      while(stillgoing){
	for(int ip=pstart;ip<=pend;ip++){
	  int i = point[ip][0];
	  int j = point[ip][1];
	  float height = tprofile[i][j];
	  for(int ii=-1; ii<2; ii++){
	    int is = ii+i;
	    for(int jj=-1; jj<2; jj++){
	      int js = jj+j;
	      if(is>=0&&is<nbins&&js>=0&&js<nbins){
		if(!assigned[is][js]){
		  if(tprofile[is][js]<height*2.0){
		    assigned[is][js]=true;
		    peakenergy+=tprofile[is][js];
		    for(int i=0;i<=maxLayers;++i)longitudinalProfile[i]+=tlprofile[is][js][i];
		    xbar+= (is-ioffset)*tprofile[is][js];
		    ybar+= (js-ioffset)*tprofile[is][js];
		    xxbar+= (is-ioffset)*(is-ioffset)*tprofile[is][js];
		    yybar+= (js-ioffset)*(js-ioffset)*tprofile[is][js];
		    pcurrent++;
		    point[pcurrent][0] = is;
		    point[pcurrent][1] = js;
		    if(tprofile[is][js]/tprofile[ipeak][jpeak]>0.1){
		      float d = sqrt((is-ioffset)*(is-ioffset)+(js-ioffset)*(js-ioffset));
		      if(d<dmin)dmin=d;
		    }
		  }
		}
	      }
	    }
	  }
	}
	if(pcurrent==pend)stillgoing=false;
	pstart = pend+1;
	pend=pcurrent;
      }
      xbar = xbar/peakenergy;
      ybar = ybar/peakenergy;
      xxbar = xxbar/peakenergy;
      yybar = yybar/peakenergy;

      protoClusterPeaks_t peak;
      peak.du = ipeak -ioffset;
      peak.dv = jpeak -ioffset;
      peak.energy  = peakenergy;
      peak.dmin    = dmin;
      float dr2 = xxbar+yybar -xbar*xbar-ybar*ybar;
      if(dr2>0)peak.rms  = sqrt(xxbar+yybar -xbar*xbar-ybar*ybar);
      if(dr2<=0)peak.rms = -999.;
      float sum = 0;
      bool stillGoing=true;

      peak.showerDepth90    = maxLayers;
      peak.showerDepth25    = maxLayers;
      for(int i=0;i<=maxLayers && stillGoing;i++){
	sum+= longitudinalProfile[i];
	if(sum>0.25*peakenergy && peak.showerDepth25==maxLayers){
	  peak.showerDepth25= i;
	}
	if(sum>0.9*peakenergy){
	  peak.showerDepth90= i;
	  stillGoing = false;
	}
      }

      stillGoing=true;
      peak.showerStartDepth    = maxLayers;
      for(int i=0;i<=maxLayers && stillGoing;i++){
	if(longitudinalProfile[i]>0.1){
	  peak.showerStartDepth= i;
	  stillGoing = false;
	}
      }

      // ****************************************************
      // IDEALLY WOULD INCLUDE SOME EM PROFILE CHI2 HERE
      // NOT INCLUDED AT THIS STAGE AS TOO DETECTOR SPECIFIC
      // ****************************************************
      // From v2.0 the code exists to do this - but isn't used yet

      peaks.push_back(peak);
    }
    if(npeaks>3)stillfindingpeaks=false;
  }

  return;

}





ProtoCluster* ProtoCluster::TransverseProfile(int peakForProtoCluster, int maxLayers, int extraLayers){

  ProtoCluster* newCluster=NULL;


  int nbins   = 41;
  int ioffset = nbins/2; 
  int ifirst = this->FirstLayer();
  float x0   = this->GetCentroid(ifirst)[0];
  float y0   = this->GetCentroid(ifirst)[1];
  float z0   = this->GetCentroid(ifirst)[2];
  float r0   = sqrt(x0*x0+y0*y0+z0*z0);
  float ux   = x0/r0;
  float uy   = y0/r0;
  float uz   = z0/r0;
  float utx  = uy/sqrt(ux*ux+uy*uy);
  float uty  = -ux/sqrt(ux*ux+uy*uy);
  float utz  = 0;
  float vtx=0;
  float vty=0;
  float vtz=0;
  if(utx!=0&&uz!=0){
    vtx  = -uty/utx;
    vtz  = (uty*ux-uy*utx)/utx/uz;
    vty  = 1;
    float v    = sqrt(vtx*vtx+vty*vty+vtz*vtz);
    vtx = vtx/v;
    vty = vty/v;
    vtz = vtz/v;
  }
  if(utx==0&&uz!=0){
    vtx  = 1;
    vtz  = ux/uz;
    vty  = 0;
    float v    = sqrt(vtx*vtx+vty*vty+vtz*vtz);
    vtx = vtx/v;
    vty = vty/v;
    vtz = vtz/v;
  }
  if(utx==0&&uz==0){
    vtx  = 0;
    vtz  = 1;
    vty  = 0;
  }

  int  done[nbins][nbins];
  bool assigned[nbins][nbins];
  float tprofile[nbins][nbins];
  std::vector<MyCaloHitExtended*> tlprofile[nbins][nbins][MAX_NUMBER_OF_LAYERS];

  for(int i=0; i<nbins; i++){
    for(int j=0; j<nbins; j++){
      tprofile[i][j]=0;
      done[i][j]=assigned[i][j]=false;
    }
  }

  bool first = true;
  float pixelsize = 10.;
  int  endLayer = maxLayers;
  if(extraLayers>0)endLayer+=extraLayers;
  for(int i=0; i<=endLayer; i++){
    for(int ihit =0; ihit < this->hitsInLayer(i);++ihit){
      MyCaloHitExtended* hiti = this->hitInLayer(i,ihit);
      if(first){
	pixelsize  = hiti->getHitSizeZ();
      }
      float dx = (hiti->getPosition()[0] -x0)/pixelsize;
      float dy = (hiti->getPosition()[1] -y0)/pixelsize;
      float dz = (hiti->getPosition()[2] -z0)/pixelsize;
      float dut = dx*utx+dy*uty+dz*utz;
      float dvt = dx*vtx+dy*vty+dz*vtz;
      int idut  = static_cast<int>(dut+0.5+ioffset);
      int idvt  = static_cast<int>(dvt+0.5+ioffset);
      if(idut>=0&&idut<nbins&&idvt>=0&&idvt<nbins){
	if(i<=maxLayers)tprofile[idut][idvt] += hiti->getEnergyEM();  
	if(i<=maxLayers+extraLayers)tlprofile[idut][idvt][i].push_back(hiti);  
      }
    }
  }

  // mask low ph region
  float threshold = 0.025;
  for(int i=0; i<nbins; i++){
    for(int j=0; j<nbins; j++){
      if(tprofile[i][j]<threshold)assigned[i][j]=true;
    }
  }
  // Search for peaks in profile
  // 
  bool stillfindingpeaks=true;
  int  npeaks=0;
  float dmin=9999.;

  while(stillfindingpeaks){
    float peakheight = 0.;
    int   ipeak =0;
    int   jpeak =0;
    float peakenergy=0;
    float xbar=0;
    float ybar=0;
    float xxbar=0;
    float yybar=0;
    std::vector<MyCaloHitExtended*> longitudinalProfile[MAX_NUMBER_OF_LAYERS];
    
    for(int i=1; i<nbins-1; i++){
      for(int j=1; j<nbins-1; j++){
	if(!assigned[i][j]){
	  if(tprofile[i][j]>peakheight){
	    peakheight = tprofile[i][j];
	    ipeak = i;
	    jpeak = j;
	  }
	}
      }
    }
    
    if(peakheight<threshold){
      stillfindingpeaks=false;
    }else{
      int point[1000][2];
      int pstart = 0;
      int pend   = 0;
      int pcurrent = pend;
      point[0][0] = ipeak;
      point[0][1] = jpeak;
      peakenergy=tprofile[ipeak][jpeak];
      
      for(int i=0;i<=endLayer;i++){
	for(unsigned int ihit=0;ihit<tlprofile[ipeak][jpeak][i].size();ihit++)
	  longitudinalProfile[i].push_back(tlprofile[ipeak][jpeak][i][ihit]);
      }

      xbar=(ipeak-ioffset)*tprofile[ipeak][jpeak];
      ybar=(jpeak-ioffset)*tprofile[ipeak][jpeak];
      xxbar=(ipeak-ioffset)*(ipeak-ioffset)*tprofile[ipeak][jpeak];
      yybar=(jpeak-ioffset)*(jpeak-ioffset)*tprofile[ipeak][jpeak];
      assigned[ipeak][jpeak]=true;
      dmin = sqrt(   (ipeak-ioffset)*(ipeak-ioffset)+
                     (jpeak-ioffset)*(jpeak-ioffset));
      npeaks++;
      float stillgoing = true;
      while(stillgoing){
	for(int ip=pstart;ip<=pend;ip++){
	  int i = point[ip][0];
	  int j = point[ip][1];
	  float height = tprofile[i][j];
	  for(int ii=-1; ii<2; ii++){
	    int is = ii+i;
	    for(int jj=-1; jj<2; jj++){
	      int js = jj+j;
	      if(is>=0&&is<nbins&&js>=0&&js<nbins){
		if(!assigned[is][js]){
		  if(tprofile[is][js]<height*1.5){
		    assigned[is][js]=true;
		    peakenergy+=tprofile[is][js];
		    for(int i=0;i<=endLayer;i++){
		      for(unsigned int ihit=0;ihit<tlprofile[is][js][i].size();ihit++)
			longitudinalProfile[i].push_back(tlprofile[is][js][i][ihit]);
		    }
		    xbar+= (is-ioffset)*tprofile[is][js];
		    ybar+= (js-ioffset)*tprofile[is][js];
		    xxbar+= (is-ioffset)*(is-ioffset)*tprofile[is][js];
		    yybar+= (js-ioffset)*(js-ioffset)*tprofile[is][js];
		    pcurrent++;
		    point[pcurrent][0] = is;
		    point[pcurrent][1] = js;
		    float d = sqrt((is-ioffset)*(is-ioffset)+(js-ioffset)*(js-ioffset));
		    if(d<dmin)dmin=d;
		  }
		}
	      }
	    }
	  }
	}
	if(pcurrent==pend)stillgoing=false;
	pstart = pend+1;
	pend=pcurrent;
      }

      // for the peak we actually want make a protocluster
      if(npeaks==peakForProtoCluster+1){
	for(int i=0;i<=maxLayers+extraLayers;i++){
	  for(unsigned int ihit = 0; ihit<longitudinalProfile[i].size();ihit++){
	    MyCaloHitExtended* hiti = longitudinalProfile[i][ihit];
	    if(newCluster==NULL){
	      newCluster = new ProtoCluster(hiti);
	    }else{
	      newCluster->AddHit(hiti);
	    }
	  }
	}
      }
    }
    if(npeaks>3)stillfindingpeaks=false;
  }

  return newCluster;

}

float* ProtoCluster::GetInitialDir(){
  return _initialDir;
}

float* ProtoCluster::GetCurrentDir(){
  return _currentDir;
}


fitResult ProtoCluster::FitPoints(const std::vector<vec4> & points){

  fitResult result;

  double sx  = 0.;
  double sy  = 0.;
  double su  = 0.; 
  double sv  = 0.; 
  double sz  = 0.;
  double svu = 0.;
  double szu = 0.;
  double suu = 0.;
  double av,bv,bu,az,bz;
  double nc = 0;
  float z;
  float u;
  float v;

  result.ok = true;

  // perform linear regression of x vs d and y vs d and z vs d (assuming same error on all hits)
  for(unsigned int i = 0; i<points.size(); ++i){
    nc += 1.;
    sx += points[i].x;
    sy += points[i].y; 
    sz += points[i].z; 
  }

  if(nc<=1){
    result.ok = false;
    return result;
  }

  
  float meanx = sx/nc;
  float meany = sy/nc;
  float meanz = sz/nc;

  fitResult zfit;
  zfit.ok = false;
  if(fabs(meanz)>_zOfEndcap){
    if(meanz<0)zfit = this->FitPointsEndCap(points,1);    
    if(meanz>0)zfit = this->FitPointsEndCap(points,2);    
    return zfit;
  }


  float rxy = sqrt(meanx*meanx+meany*meany);
  float cost = meanx/rxy;
  float sint = meany/rxy;
  
  sz = 0.;
  for(unsigned int i = 0; i<points.size(); ++i){
    u =  cost*points[i].x +sint*points[i].y;
    v = -sint*points[i].x +cost*points[i].y;
    z = points[i].z; 
    su  += u;
    sv  += v;
    sz  += z;
    svu += v*u;
    szu += z*u;
    suu += u*u;
  }


  if(su*su-nc*suu==0){
    result.ok = false;
    return result;
  }
  av = (su*sv-nc*svu)/(su*su-nc*suu);
  bv = (sv - av*su)/nc;
  bu = -bv/av;

  az = (su*sz-nc*szu)/(su*su-nc*suu);
  bz = (sz - az*su)/nc;

  // intercept on v = 0
  float bzprime = az*bu+bz;

  float r = sqrt(av*av+az*az+1.0);
  float diru = 1.0/r;
  float dirv = av/r;
  float dirz = az/r;

  float dirx = cost*diru - sint*dirv;
  float diry = sint*diru + cost*dirv;

  result.dir[0] = dirx;
  result.dir[1] = diry;
  result.dir[2] = dirz;

  result.intercept[0] =  cost*bu;
  result.intercept[1] =  sint*bu;
  result.intercept[2] =  bzprime;
  float bx = result.intercept[0];
  float by = result.intercept[1];

  float ri = sqrt(bx*bx+by*by+bzprime*bzprime);

  result.dirDotR = ( result.dir[0]*result.intercept[0]+
                           result.dir[1]*result.intercept[1]+
                           result.dir[2]*result.intercept[2])/ri ;
  if(result.dirDotR< 0){
    for(int i=0;i<3;i++)result.dir[i] = -result.dir[i];
    result.dirDotR = -result.dirDotR;
  }
  // now calculate something like a chi2 to straight line fit

  float chi2v = 0.;
  float chi2z = 0.;
  float rms = 0.;

  nc =0;
  double sl=0.;
  double sa=0.;
  double sll=0.;
  double sal=0.;
  
    
  for(unsigned int i = 0; i<points.size(); ++i){
    nc++;
    u =  cost*points[i].x +sint*points[i].y;
    v = -sint*points[i].x +cost*points[i].y;
    z = points[i].z; 
    float error = points[i].padSize/3.46; 
    float chiv = (v-av*u-bv)/error;
    chi2v += chiv*chiv;
    float chiz = (z-az*u-bz)/error;
    chi2z += chiz*chiz;
    float dx = points[i].x - result.intercept[0];
    float dy = points[i].y - result.intercept[1];
    float dz = points[i].z - result.intercept[2];
    float ddx = result.dir[1]*dz - result.dir[2]*dy;
    float ddy = result.dir[2]*dx - result.dir[0]*dz;
    float ddz = result.dir[0]*dy - result.dir[1]*dx;      
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    rms += distSq;

    float ai = (points[i].x-result.intercept[0])*result.dir[0]+
               (points[i].y-result.intercept[1])*result.dir[1]+ 
               (points[i].z-result.intercept[2])*result.dir[2];
    float li = (float)points[i].ilayer;
    sa += ai;
    sl += li;
    sal += ai*li;
    sll += li*li;
  } 
  if((sl*sl-nc*sll)!=0){
    float a = (sl*sa-nc*sal)/(sl*sl-nc*sll);
    //std::cout << " a : " << a << std::endl;
    if(a<0){
      for(unsigned int i=0;i<3;i++)result.dir[i]=-result.dir[i];
    }
  }

  result.chi2 = (chi2v+chi2z)/nc;
  result.rms = sqrt(rms/nc);

  return result;

}




fitResult ProtoCluster::FitPointsEndCap(const std::vector<vec4> & points, int endcap){

  fitResult result;

  double sx  = 0.;
  double sy  = 0.;
  double sz  = 0.;
  double sxz = 0.;
  double syz = 0.;
  double szz = 0.;
  double ax,bx,ay,by;
  double nc = 0;
  float x;
  float y;
  float z;

  result.ok = true;

  for(unsigned int i = 0; i<points.size(); ++i){
    nc += 1.;
    x = points[i].x; 
    y = points[i].y; 
    z = points[i].z; 
    if(endcap==1)z +=  _zOfEndcap; 
    if(endcap==2)z -=  _zOfEndcap; 
    sx  += x;
    sy  += y;
    sz  += z;
    szz += z*z;
    sxz += x*z;
    syz += y*z;
  }

  if(sz*sz-nc*szz == 0){
    result.ok = false;
    return result;
  }

  ax = (sz*sx-nc*sxz)/(sz*sz-nc*szz);
  bx = (sx - ax*sz)/nc;

  ay = (sz*sy-nc*syz)/(sz*sz-nc*szz);
  by = (sy - ay*sz)/nc;

  float r = sqrt(ax*ax+ay*ay+1.0);
  float dirx = ax/r;
  float diry = ay/r;
  float dirz = 1.0/r;

  result.dir[0] = dirx;
  result.dir[1] = diry;
  result.dir[2] = dirz;

  result.intercept[0] =  bx;
  result.intercept[1] =  by;
  if(endcap==1)result.intercept[2] =  -_zOfEndcap;
  if(endcap==2)result.intercept[2] =  +_zOfEndcap;


  bx = result.intercept[0];
  by = result.intercept[1];
  float ri = sqrt(bx*bx+by*by+_zOfEndcap*_zOfEndcap);

  result.dirDotR = ( result.dir[0]*result.intercept[0]+
                     result.dir[1]*result.intercept[1]+
                     result.dir[2]*result.intercept[2])/ri ;

  if(result.dirDotR<0){
    for(int i=0;i<3;i++)result.dir[i] = -result.dir[i];
    result.dirDotR = -result.dirDotR;
  }
  // now calculate something like a chi2 to straight line fit

  float chi2x = 0.;
  float chi2y = 0.;
  float rms = 0.;

  double sa=0.;
  double sl=0.;
  double sal=0.;
  double sll=0.;
  nc =0;
  for(unsigned int i = 0; i<points.size(); ++i){
    nc++;
    x = points[i].x;
    y = points[i].y; 
    z = points[i].z; 
    float dx = points[i].x - result.intercept[0];
    float dy = points[i].y - result.intercept[1];
    float dz = points[i].z - result.intercept[2];
    float error = points[i].padSize/3.46; 
    float chix = (dx-ax*dz)/error;
    chi2x += chix*chix;
    float chiy = (dy-ay*dz)/error;
    chi2y += chiy*chiy;
    float ddx = result.dir[1]*dz - result.dir[2]*dy;
    float ddy = result.dir[2]*dx - result.dir[0]*dz;
    float ddz = result.dir[0]*dy - result.dir[1]*dx;      
    float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
    rms += distSq;
    float ai = (points[i].x-result.intercept[0])*result.dir[0]+
               (points[i].y-result.intercept[1])*result.dir[1]+ 
               (points[i].z-result.intercept[2])*result.dir[2];
    float li = (float)points[i].ilayer;
    sa += ai;
    sl += li;
    sal += ai*li;
    sll += li*li;
  } 

  if((sl*sl-nc*sll)!=0){
    float a = (sl*sa-nc*sal)/(sl*sl-nc*sll);
    if(a<0){
      for(unsigned int i=0;i<3;i++)result.dir[i]=-result.dir[i];
    }
  }
 
  result.chi2 = (chi2x+chi2y)/nc;
  result.rms = sqrt(rms/nc);

  return result;

}


fitResult ProtoCluster::FitAllHits(){

  std::vector<vec4>points;
  fitResult result;
  result.ok = false;
  if(_nHits<2)return result;
  
  for(int ilayer=_firstLayer;ilayer<=_lastLayer;++ilayer){
    for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size()>0;++ihit){
      vec4 point;
      point.x = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[0];
      point.y = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[1];
      point.z = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[2];
      point.ilayer = ilayer;
      point.padSize = _hitsByPseudoLayer[ilayer][ihit]->getHitSizeZ();
      points.push_back(point);
    }
  }

  result = this->FitPoints(points);

  return result;
}


fitResult ProtoCluster::FitLayers(int istart, int iend){

  std::vector<vec4>points;
  fitResult result;
  result.ok = false;
  if(_nHits<2)return result;
  
  for(int ilayer=istart;ilayer<=iend;++ilayer){
    for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size()>0;++ihit){
      vec4 point;
      point.x = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[0];
      point.y = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[1];
      point.z = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[2];
      point.ilayer = ilayer;
      point.padSize = _hitsByPseudoLayer[ilayer][ihit]->getHitSizeZ();
      points.push_back(point);
    }
  }

  result = this->FitPoints(points);

  return result;
}


fitResult ProtoCluster::FitStart(int nlayers){

  fitResult result;
  result.ok = false;
  if(_nHits<4)return result;

  int ilayer = _firstLayer;
  int layerCount = 0;
  std::vector<vec4>points;

  while(ilayer<=_lastLayer&&layerCount<nlayers){
    if(_hitsByPseudoLayer[ilayer].size()>0){
      layerCount++;
      for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size()>0;++ihit){
	vec4 point;
	point.x = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[0];
	point.y = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[1];
	point.z = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[2];
	point.ilayer = ilayer - _firstLayer;
	point.padSize = _hitsByPseudoLayer[ilayer][ihit]->getHitSizeZ();
	points.push_back(point);
      }
    }
    ilayer++;
  }

  if(layerCount<4)return result;

  result = this->FitPoints(points);

  return result;
}


fitResult ProtoCluster::FitEnd(int nlayers){

  fitResult result;
  result.ok = false;

  if(_nHits<4)return result;
  int ilayer = _lastLayer;
  int layerCount = 0;
  std::vector<vec4>points;

  while(ilayer>=0&&layerCount<nlayers){
    if(_hitsByPseudoLayer[ilayer].size()>0){
      layerCount++;
      for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size()>0;++ihit){
	vec4 point;
	point.x = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[0];
	point.y = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[1];
	point.z = _hitsByPseudoLayer[ilayer][ihit]->getPosition()[2];
	point.ilayer = ilayer;
	point.padSize = _hitsByPseudoLayer[ilayer][ihit]->getHitSizeZ();
	points.push_back(point);
      }
    }
    ilayer--;
  }

  if(layerCount<2||points.size()<4)return result;

  result = this->FitPoints(points);

  return result;
}


bool ProtoCluster::IsPrimaryPhoton(){

  if(_nHits==0)return false;
  if(_clusterClassification==PC_PRIMARY_PHOTON)return true;
  return false;

};

bool ProtoCluster::IsAttachable(){

  if(_defunct)return false;
  bool retval = !this->IsPrimaryPhoton();
  if(_mipFraction>0.7)return true;
  if(_clusterRms<5.0)return true;

  return retval;

};



bool ProtoCluster::IsPhoton(){

  if(_nHits==0)return false;
  if(_clusterClassification==PC_PRIMARY_PHOTON)return true;
  if(_clusterClassification==PC_SECONDARY_PHOTON)return true;
  return false;

};

void ProtoCluster::IsPrimaryPhoton(bool tf){

  if(tf)_clusterClassification=PC_PRIMARY_PHOTON;

};




void ProtoCluster::Assimilate(ProtoCluster* pC){
  
  for(int i=0;i<pC->Hits();i++){
    MyCaloHitExtended* hiti = pC->Hit(i);
    this->AddHit(hiti);
  }
  for(int i=0;i<pC->IsolatedHits();i++){
    MyCaloHitExtended* hiti = pC->IsolatedHit(i);
    this->AddIsolatedHit(hiti);
  }
  this->Update(_lastLayer);

};


float ProtoCluster::MipFraction(){

  if(_nHits==0)return 0;
  float mipFraction = static_cast<float>(_nMips)/static_cast<float>(_nHits);

  return mipFraction;

};



float ProtoCluster::MipFraction(int startLayer, int endLayer){

  if(_nHits==0)return 0;
  if(startLayer>=endLayer)return 0;

  float hitCount = 0.;
  float mipCount = 0.;
  if(startLayer<1)startLayer=1;
  if(endLayer>=MAX_NUMBER_OF_LAYERS)endLayer=MAX_NUMBER_OF_LAYERS-1;

  for(int ilayer=startLayer;ilayer<=endLayer;++ilayer){
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      hitCount+=1.;
      if(hiti->isCandidateMip())mipCount+=1.;
    }
  }

  float mipFraction = 0.;
  if(mipCount>0)mipFraction=mipCount/hitCount;

  return mipFraction;

};


float ProtoCluster::DistanceToClosestHit(fitResult fit, int iStartLayer, int iEndLayer){

  float retVal = 99999.;
  int istart = iStartLayer;
  int iend   = iEndLayer;
  if(istart<_firstLayer)istart=_firstLayer;
  if(iend>_lastLayer)iend=_lastLayer;
  float r[3];
  r[0] = fit.intercept[0];
  r[1] = fit.intercept[1];
  r[2] = fit.intercept[2];
  float dir[3];
  dir[0] = fit.dir[0];
  dir[1] = fit.dir[1];
  dir[2] = fit.dir[2];

  for(int ilayer=istart;ilayer<=iend;++ilayer){
    float dr[3];
    for(unsigned int ihit=0; ihit < _hitsByPseudoLayer[ilayer].size(); ++ihit){
      MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
      dr[0] = hiti->getPosition()[0]-r[0];
      dr[1] = hiti->getPosition()[1]-r[1];
      dr[2] = hiti->getPosition()[2]-r[2];
      float ddx = dr[1]*dir[2]-dr[2]*dir[1];
      float ddy = dr[2]*dir[0]-dr[0]*dir[2];
      float ddz = dr[0]*dir[1]-dr[1]*dir[0];
      float d = sqrt( ddx*ddx + ddy*ddy + ddz*ddz );
      if(d<retVal)retVal=d;
    }
  }

  return retVal;


}



float ProtoCluster::DistanceToClosestHit(ProtoCluster* pC){

  float retVal = 99999.;

  int nhits = pC->Hits();

  for(int ihit= 0; ihit<nhits;++ihit){
    MyCaloHitExtended* hiti = pC->Hit(ihit);
    float r[3];
    r[0] = hiti->getPosition()[0];
    r[1] = hiti->getPosition()[1];
    r[2] = hiti->getPosition()[2];
    for(unsigned int jhit=0; jhit < _hits.size(); ++jhit){
      MyCaloHitExtended* hitj = _hits[jhit];
      float dx = hitj->getPosition()[0]-r[0];
      float dy = hitj->getPosition()[1]-r[1];
      float dz = hitj->getPosition()[2]-r[2];
      float d = sqrt( dx*dx + dy*dy + dz*dz );
      if(d<retVal)retVal=d;
    }
  }

  return retVal;


}

float ProtoCluster::TruePhotonContribution(){ 

  float total = 0;

  for(unsigned int ihit= 0; ihit<_hits.size();++ihit){
    MyCaloHitExtended* hiti = _hits[ihit];
    if(hiti->getPhoton())total+=hiti->getEnergyEM();
  }
  
  return total;

}

float ProtoCluster::FractionOfCloseHits(ProtoCluster* pC, float distance){

  float retVal = 0.;

  int nhits = pC->Hits();
  int nclose = 0;

  for(int ihit= 0; ihit<nhits;++ihit){
    MyCaloHitExtended* hiti = pC->Hit(ihit);
    float r[3];
    r[0] = hiti->getPosition()[0];
    r[1] = hiti->getPosition()[1];
    r[2] = hiti->getPosition()[2];
    bool notFound = true;
    for(unsigned int jhit=0; notFound && jhit < _hits.size(); ++jhit){
      MyCaloHitExtended* hitj = _hits[jhit];
      float dx = hitj->getPosition()[0]-r[0];
      float dy = hitj->getPosition()[1]-r[1];
      float dz = hitj->getPosition()[2]-r[2];
      float d = sqrt( dx*dx + dy*dy + dz*dz );
      if(d<distance){
	nclose++;
	notFound = false;
      }
    }
  }
  
  retVal = static_cast<float>(nclose)/static_cast<float>(nhits);

  return retVal;

}





float ProtoCluster::DistanceToClosestCentroid(fitResult fit, int iStartLayer, int iEndLayer){

  float r[3];
  r[0] = fit.intercept[0];
  r[1] = fit.intercept[1];
  r[2] = fit.intercept[2];
  float dir[3];
  dir[0] = fit.dir[0];
  dir[1] = fit.dir[1];
  dir[2] = fit.dir[2];
  float retVal = 99999.;
  int istart = iStartLayer;
  int iend   = iEndLayer;
  if(istart<_firstLayer)istart=_firstLayer;
  if(iend>_lastLayer)iend=_lastLayer;
  for(int ilayer=istart;ilayer<=iend;++ilayer){
    if(_hitsByPseudoLayer[ilayer].size()>0){
      float dr[3];
      dr[0] = _centroidX[ilayer]-r[0];
      dr[1] = _centroidY[ilayer]-r[1];
      dr[2] = _centroidZ[ilayer]-r[2];
      float ddx = dr[1]*dir[2]-dr[2]*dir[1];
      float ddy = dr[2]*dir[0]-dr[0]*dir[2];
      float ddz = dr[0]*dir[1]-dr[1]*dir[0];
      float d = sqrt( ddx*ddx + ddy*ddy + ddz*ddz );
      if(d<retVal)retVal=d;
    }
  }

  return retVal;

}



void ProtoCluster::SoftPhotonID(){

  if(_fixedPhoton)_clusterClassification = PC_PRIMARY_PHOTON;
  if(_nHits==0)return;
  if(_nHits>=25)return;
  if(_energyEM>1.0)return;
  if(_firstLayer>15)return;
  bool softPhotonID = false;
  if(_dCosR>0.90)softPhotonID=true;
  if(_dCosR>0.80&&_energyEM<0.5)softPhotonID=true;
  if(softPhotonID)_clusterClassification = PC_PRIMARY_PHOTON;
  return;
}

void ProtoCluster::TagFixedPhoton(){

  _fixedPhoton = true;
  _clusterClassification = PC_PRIMARY_PHOTON;
  return;

}


void ProtoCluster::ClassifyYourself(){


  if(_nHits==0){
    _clusterClassification = PC_UNASSOCIATED_TRACK;
    return;
  }

  _mipFraction = this->MipFraction();
  fitResult fitAll = this->FitAllHits();
  if(fitAll.ok){
    _clusterRms = fitAll.rms;
  }else{
    _clusterRms = 0.;
  }

  float dx = _centroidX[_lastLayer]-_centroidX[_firstLayer];
  float dy = _centroidY[_lastLayer]-_centroidY[_firstLayer];
  float dz = _centroidZ[_lastLayer]-_centroidZ[_firstLayer];
  _showerLength = sqrt(dx*dx+dy*dy+dz*dz);
  
  float cosTheta = fabs(_centroidZ[_firstLayer])/
    sqrt(_centroidX[_firstLayer]*_centroidX[_firstLayer]+
	 _centroidY[_firstLayer]*_centroidY[_firstLayer]+
	 _centroidZ[_firstLayer]*_centroidZ[_firstLayer]);
  
  float sinTheta = sqrt(1.0-cosTheta*cosTheta);
  float rdotn = sinTheta;
  if(fabs(_centroidZ[_firstLayer])>_zOfEndcap-50.0)rdotn = cosTheta;
  
  _layer90 = 0;
  float emax = 0.;
  int   imax = 0;
  float e90=0.;
  bool  f90=false;
  float etot = 0.;
  for(int i=_firstLayer; i<=_lastLayer ;i++)etot+=_energy[i];
  
  _hitLayers = 0;
  for(int i=_firstLayer; i<=_lastLayer ;i++){
    if(_hitsByPseudoLayer[i].size()>0)_hitLayers++;
    e90+= _energy[i];
    if(_energy[i]>emax && _hitsByPseudoLayer[i].size()>0){
      emax = _energy[i];
      imax = i;
    }
    if(e90>0.9*etot && !f90){
      f90 = true;
     _layer90 = i;
    }
  }
  if(_layer90==0){
    streamlog_out( WARNING ) << " ProtoCluster::ClassifyYourself : Shower only in layer 0 - gear problem ? " << " " << _energyEM << " " << _firstLayer << " " << _lastLayer << std::endl;
  }

  _layerMax = imax;

  _showerMax= 0;
  if(imax>=_firstLayer){
    dx = _centroidX[imax]-_centroidX[_firstLayer]; 
    dy = _centroidY[imax]-_centroidY[_firstLayer]; 
    dz = _centroidZ[imax]-_centroidZ[_firstLayer]; 
    _showerMax = sqrt(dx*dx+dy*dy+dz*dz);
  }


  _showerLayer=_lastLayer;
  int currentShowerLayers = 0;
  bool showerStartNotFound = true;
  int showerStart=_lastLayer;
  for(int i=_firstLayer; i<=_lastLayer && showerStartNotFound ;i++){
    float localFrac=0.;
    if(_hitsByPseudoLayer[i].size()>0)localFrac = _mipHitsByPseudoLayer[i].size()/_hitsByPseudoLayer[i].size();
    if(localFrac>0.8){
      currentShowerLayers = 0;
    }else{
      currentShowerLayers++;
    }
    if(currentShowerLayers>=2){
      showerStartNotFound=false;
      if(_hitsByPseudoLayer[i].size()>0)showerStart = i;
    }
  }
  
  // found two consecutive shower layers 
  // now go back until find two consecutive MIP layers
  if(showerStartNotFound==false){
    _showerLayer = showerStart;
    int currentMipLayers = 0;
    bool mipEndNotFound = true;
    for(int i=showerStart; i>= _firstLayer && mipEndNotFound; --i){
      if(_hitsByPseudoLayer[i].size()>0){
	float localFrac = 0;
	localFrac = _mipHitsByPseudoLayer[i].size()/_hitsByPseudoLayer[i].size();
	if(localFrac>0.8){
	  currentMipLayers++;
	}else{
	  currentMipLayers=0;
	  _showerLayer = i;
	}
	if(currentMipLayers>=2){
	  mipEndNotFound=false;
	}
      }
    }  
  }



  _dCosR = 0.;
  if(fitAll.ok){
    float rx = _centroidX[_firstLayer];
    float ry = _centroidY[_firstLayer];
    float rz = _centroidZ[_firstLayer];
    float r = sqrt(rx*rx+ry*ry+rz*rz);
    _dCosR = rx/r*fitAll.dir[0]+ry/r*fitAll.dir[1]+rz/r*fitAll.dir[2];
  }



  float mipCut=0.9;
  float dCosCut=_photonID_dcoscut;
  float lCut1 =  0.;
  float lCut2 = 40.;
  float rmsCut = 40.;
  float layerCut1 = 5.;
  float layerCut2 = 40.;


  if(_energyEM>1.5){
    lCut1 =1.;
    mipCut = 0.7;
    dCosCut += 0.01;
  }
  if(_energyEM>3){
    lCut1  =3.0;
    mipCut = 0.6;
  }
  if(_energyEM>7.5){
    mipCut = 0.4;
  }
  if(_energyEM>15){
    mipCut = 0.3;
  }
  if(_energyEM>40){
    mipCut = 0.3;
    rmsCut = 50.;
    layerCut2 = 50.;
  }
  float firstLayerInX0= _calThicknessByPhysicalLayer[_firstLayer]/rdotn/3.5;

  if(!_associatedTrack&&_firstLayer<_nEcalLayers){
    if( firstLayerInX0<10                       &&
	_mipFraction<mipCut                     &&
	_clusterRms<rmsCut                      &&
	(_layerMax-_firstLayer)>lCut1*rdotn     &&
	(_layerMax-_firstLayer)<lCut2*rdotn     &&
	_dCosR>dCosCut                          &&
	(_layer90-_firstLayer)>layerCut1*rdotn  && 
	(_layer90-_firstLayer)<layerCut2*rdotn  &&
	_layer90<_nEcalLayers+10    ){


      _clusterClassification = PC_PRIMARY_PHOTON;
      if(firstLayerInX0>10)_clusterClassification = PC_SECONDARY_PHOTON;
    }
  }
  if(_fixedPhoton)_clusterClassification = PC_PRIMARY_PHOTON;


  if(_energyEM>7.0&&_energyEM<-7.2){
    std::cout << " FIRST : " << _firstLayer << std::endl;
    std::cout << " FIRSTX0 : " << _calThicknessByPhysicalLayer[_firstLayer]/rdotn/3.5 << std::endl;
    std::cout << " RDOTN : " << rdotn << std::endl;
    std::cout << " CLUSTER E : " << _energyEM << "(" << _energyHAD << ")"<<std::endl;
    std::cout << " MIP   : " << _mipFraction << " < " << mipCut << " ? " << std::endl;
    std::cout << " RMS   : " << _clusterRms  << " < " << rmsCut << " ? " << std::endl;
    std::cout << " DCOSR : " << _dCosR << " > " << dCosCut << " ? " << std::endl;
    std::cout << " LAY90 : " << layerCut1*rdotn << " < " << (_layer90-_firstLayer) << " < " << layerCut2*rdotn << std::endl;
    std::cout << " LAYMAX : " << lCut1*rdotn << " < " << (_layerMax-_firstLayer) << " < " << lCut2*rdotn << std::endl;
    if(_clusterClassification ==PC_PRIMARY_PHOTON)std::cout << " PRIMARY PHOTON CLUSTER " << std::endl;
    if(_clusterClassification ==PC_SECONDARY_PHOTON)std::cout << " SECONDARY PHOTON CLUSTER " << std::endl;
    std::cout << this->IsPhoton() << std::endl;
  }


  // fit part of cluster from start to the shower start
  fitResult mipFit;
  mipFit.ok = false;
  if(_showerLayer-_firstLayer>3)mipFit = this->FitLayers(_firstLayer,_showerLayer-1);
  _mipFitResult = mipFit;

}

void ProtoCluster::SetID(int id){
  _ID = id;
  _earliestRelationID = id;
}



void ProtoCluster::AddDaughter(ProtoCluster* cluster){

  _daughters.push_back(cluster);
  _connected = true;
  int relID = cluster->getEarliestRelationID();
  ProtoCluster* rel = cluster->getEarliestRelation();
  if(relID>_earliestRelationID){
    cluster->SetEarliestRelationID(_earliestRelationID);
    cluster->SetEarliestRelation(_earliestRelation);
  }
  if(relID<_earliestRelationID){
    this->SetEarliestRelationID(relID);
    this->SetEarliestRelation(rel);
  }

}

ProtoCluster* ProtoCluster::FragmentMip(){
  
  ProtoCluster* mipCluster = NULL;
  
  // start at end of cluster looking for protruding mip
  bool mipStartNotFound = true;
  int mipEnd   = _lastLayer;
  int mipStart = _lastLayer;
  int firstMip = _lastLayer;
  int currentShowerLayers = 0;
  for(int i=_lastLayer; i >= _firstLayer && mipStartNotFound ;--i){
    float localFrac=0.;
    if(_hitsByPseudoLayer[i].size()>0)localFrac = _mipHitsByPseudoLayer[i].size()/_hitsByPseudoLayer[i].size();
    if(localFrac>0.8 && _hitsByPseudoLayer[i].size()<=2){
      currentShowerLayers = 0;
      firstMip = i;
    }else{
      currentShowerLayers++;
    }
    if(currentShowerLayers>=2){
      mipStartNotFound=false;
      mipStart = firstMip;
    }
  }
  
  
  if(mipEnd-mipStart>7){
    fitResult mipFit = this->FitLayers(mipStart+2,mipEnd);
    if(mipFit.chi2<2.5){
      for(int ilayer=mipStart+2;ilayer<=mipEnd;++ilayer){
	for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size();++ihit){
	  MyCaloHitExtended * hiti = _hitsByPseudoLayer[ilayer][ihit];
	  if(mipCluster==NULL){
	    mipCluster = new ProtoCluster(hiti);
	  }else{
	    mipCluster->AddHit(hiti);
	  }
	}
      }
    }
  }
  return mipCluster;
  
}

ProtoCluster* ProtoCluster::FragmentPhoton(){
  
  ProtoCluster* photonCluster = NULL;
  
  float distance = 1000.;
  
  if(!_trackSeeded && !_associatedTrack)return photonCluster;
  int missed = 0;
  bool mipregion1   = true;
  bool mipregion2   = false;
  bool showerregion = false;
  int mipcount = 0;
  int showercount = 0;
  
  int mip1s = -999;
  int mip1e = -999;
  int mip2s = -999;
  int mip2e = -999;
  int showers = -999;
  int showere = -999;

  //  loop over layers until no hits found on track projection
  bool carryon = true;
  for(int ilayer=1; ilayer <= _lastLayer && carryon;ilayer++){
    bool trackhitfound = false;
    bool miptrackhitfound = false;
    bool showertrackhitfound = false;
    for(unsigned int ihit = 0; ihit < _hitsByPseudoLayer[ilayer].size();ihit++){
      MyCaloHitExtended *hiti = _hitsByPseudoLayer[ilayer][ihit];
      float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
      if(distToTrackSeed<1.0){
	trackhitfound=true;
	if(hiti->isCandidateMip())miptrackhitfound = true;
	if(!hiti->isCandidateMip())showertrackhitfound = true;
      }
    }
    if(trackhitfound)missed = 0;
    if(!trackhitfound)missed++;
    // if failed to find track hits in two subsequent layers terminate 
    if(miptrackhitfound&&mipregion1)mip1e=ilayer;
    if(miptrackhitfound&&mipregion2)mip2e=ilayer;
    if(showertrackhitfound&&showerregion)showere=ilayer;

    if(miptrackhitfound&&!showertrackhitfound){
      if(mipregion1&&mip1s<0)mip1s = ilayer;
      if(mipregion1||mipregion2)showercount=0;
      if(showerregion){
	mipcount++;
	if(mipcount==2){
	  mipregion2 = true;
	  showerregion = false;
	  showercount = 0;
	}else{
	  mip2s = ilayer;
	}
      }
    }

    if(!miptrackhitfound&&showertrackhitfound){
      if(showerregion)mipcount=0;
      if(mipregion1||mipregion2){
	showercount++;
	if(showercount==2){
	  if(mipregion1){
	    showerregion = true;
	    mipregion1=false;
	    showercount = 0;
	  }
	  if(mipregion2)carryon=false;
	}else{
	  if(mipregion1)showers = ilayer;
	}
      }
    }
    if(missed>1)carryon=false;
  }

  if(mip2e>0 && ( (mip2e-mip2s>4)  || (showers<5 && (showere-showers)>200))){
    // have found MIP stub of greater than 4 layers from shower
    if(showere-showers>4&showers<20){
      // also have a respectable shower
      bool debugPrint=false;
      if(debugPrint){
	std::cout << " CANDIDATE PHOTON " << std::endl;
	std::cout << " CLUSTER  : " << _centroidX[_firstLayer] << "," << _centroidY[_firstLayer] << "," << _centroidZ[_firstLayer] << std::endl;
	std::cout << " MIP 1  REGION " << mip1s << "-" << mip1e << std::endl;
	std::cout << " SHOWER REGION " << showers << "-" << showere << std::endl;
	std::cout << " MIP 2  REGION " << mip2s << "-" << mip2e << std::endl;
      }
      carryon = true;
      for(int ilayer=1; ilayer <= _lastLayer && carryon;ilayer++){
	bool trackhitfound = false;
	bool miptrackhitfound = false;
	bool showertrackhitfound = false;
	for(unsigned int ihit = 0; ihit < _hitsByPseudoLayer[ilayer].size();ihit++){
	  MyCaloHitExtended *hiti = _hitsByPseudoLayer[ilayer][ihit];
	  float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
	  if(distToTrackSeed<1.0){
	    trackhitfound=true;
	    if(hiti->isCandidateMip())miptrackhitfound = true;
	    if(!hiti->isCandidateMip())showertrackhitfound = true;
	  }
	}
	if(trackhitfound)missed = 0;
	if(!trackhitfound)missed++;
	// if failed to find track hits in two subsequent layers terminate 
	if(debugPrint){
	  
	  if(miptrackhitfound&&!showertrackhitfound)std::cout  << "M";
	  if(miptrackhitfound&&showertrackhitfound)std::cout   << "X";
	  if(!miptrackhitfound&&showertrackhitfound)std::cout  << "S";
	  if(!miptrackhitfound&&!showertrackhitfound)std::cout << " ";
	  if(missed>1)carryon=false;
	  std::cout << std::endl;
	}
      }

      int nhits = 0;
      ProtoCluster* photonCandidate = NULL;
      for(int ilayer=showers;ilayer<=showere;++ilayer){
	for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size();++ihit){
	  MyCaloHitExtended * hiti = _hitsByPseudoLayer[ilayer][ihit];
	  float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
	  if(distToTrackSeed>1.0){
	    nhits++;
	    if(photonCandidate==NULL){
	      photonCandidate= new ProtoCluster(hiti);
	    }else{
	      photonCandidate->AddHit(hiti);
	    }
	  }
	}
      }
      
      if(nhits>5){
	photonCandidate->ClassifyYourself();
	photonCluster= photonCandidate;
      }
    }
  }
  
  return photonCluster;
  
}




ProtoCluster* ProtoCluster::AgressiveFragmentPhoton(){
  
  ProtoCluster* photonCluster = NULL;
  
  float distance = 1000.;
  
  if(!_trackSeeded && !_associatedTrack)return photonCluster;
  int missed = 0;
  bool mipregion1   = true;
  bool mipregion2   = false;
  bool showerregion = false;
  int mipcount = 0;
  int showercount = 0;
  
  int mip1s = -999;
  int mip1e = -999;
  int mip2s = -999;
  int mip2e = -999;
  int showers = -999;
  int showere = -999;

  //  loop over layers until no hits found on track projection
  bool carryon = true;
  for(int ilayer=1; ilayer <= _lastLayer && carryon;ilayer++){
    bool trackhitfound = false;
    bool miptrackhitfound = false;
    bool showertrackhitfound = false;
    for(unsigned int ihit = 0; ihit < _hitsByPseudoLayer[ilayer].size();ihit++){
      MyCaloHitExtended *hiti = _hitsByPseudoLayer[ilayer][ihit];
      float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
      if(distToTrackSeed<1.0){
	trackhitfound=true;
	if(hiti->isCandidateMip())miptrackhitfound = true;
	if(!hiti->isCandidateMip())showertrackhitfound = true;
      }
    }
    if(trackhitfound)missed = 0;
    if(!trackhitfound)missed++;
    // if failed to find track hits in two subsequent layers terminate 
    if(miptrackhitfound&&mipregion1)mip1e=ilayer;
    if(miptrackhitfound&&mipregion2)mip2e=ilayer;
    if(showertrackhitfound&&showerregion)showere=ilayer;

    if(miptrackhitfound&&!showertrackhitfound){
      if(mipregion1&&mip1s<0)mip1s = ilayer;
      if(mipregion1||mipregion2)showercount=0;
      if(showerregion){
	mipcount++;
	if(mipcount==2){
	  mipregion2 = true;
	  showerregion = false;
	  showercount = 0;
	}else{
	  mip2s = ilayer;
	}
      }
    }

    if(!miptrackhitfound&&showertrackhitfound){
      if(showerregion)mipcount=0;
      if(mipregion1||mipregion2){
	showercount++;
	if(showercount==2){
	  if(mipregion1){
	    showerregion = true;
	    mipregion1=false;
	    showercount = 0;
	  }
	  if(mipregion2)carryon=false;
	}else{
	  if(mipregion1)showers = ilayer;
	}
      }
    }
    if(missed>1)carryon=false;
  }

  if(1==1){
    // mip2e>0 && ( (mip2e-mip2s>4)  || (showers<5 && (showere-showers)>200))){
    // have found MIP stub of greater than 4 layers from shower
    if(showere-showers>4&showers<20){
      // also have a respectable shower
      bool debugPrint=true;
      if(debugPrint){
	std::cout << " CANDIDATE PHOTON " << std::endl;
	std::cout << " CLUSTER  : " << _centroidX[_firstLayer] << "," << _centroidY[_firstLayer] << "," << _centroidZ[_firstLayer] << std::endl;
	std::cout << " MIP 1  REGION " << mip1s << "-" << mip1e << std::endl;
	std::cout << " SHOWER REGION " << showers << "-" << showere << std::endl;
	std::cout << " MIP 2  REGION " << mip2s << "-" << mip2e << std::endl;
      }
      carryon = true;
      for(int ilayer=1; ilayer <= _lastLayer && carryon;ilayer++){
	bool trackhitfound = false;
	bool miptrackhitfound = false;
	bool showertrackhitfound = false;
	for(unsigned int ihit = 0; ihit < _hitsByPseudoLayer[ilayer].size();ihit++){
	  MyCaloHitExtended *hiti = _hitsByPseudoLayer[ilayer][ihit];
	  float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
	  if(distToTrackSeed<1.0){
	    trackhitfound=true;
	    if(hiti->isCandidateMip())miptrackhitfound = true;
	    if(!hiti->isCandidateMip())showertrackhitfound = true;
	  }
	}
	if(trackhitfound)missed = 0;
	if(!trackhitfound)missed++;
	// if failed to find track hits in two subsequent layers terminate 
	if(debugPrint){
	  
	  if(miptrackhitfound&&!showertrackhitfound)std::cout  << "M";
	  if(miptrackhitfound&&showertrackhitfound)std::cout   << "X";
	  if(!miptrackhitfound&&showertrackhitfound)std::cout  << "S";
	  if(!miptrackhitfound&&!showertrackhitfound)std::cout << " ";
	  if(missed>1)carryon=false;
	  std::cout << std::endl;
	}
      }

      int nhits = 0;
      ProtoCluster* photonCandidate = NULL;
      for(int ilayer=showers;ilayer<=showere;++ilayer){
	for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size();++ihit){
	  MyCaloHitExtended * hiti = _hitsByPseudoLayer[ilayer][ihit];
	  float distToTrackSeed = genericDistanceToTrackSeed(hiti,distance);
	  if(distToTrackSeed>1.0){
	    nhits++;
	    if(photonCandidate==NULL){
	      photonCandidate= new ProtoCluster(hiti);
	    }else{
	      photonCandidate->AddHit(hiti);
	    }
	  }
	}
      }
      
      if(nhits>5){
	photonCandidate->ClassifyYourself();
	photonCluster= photonCandidate;
      }
    }
  }
  
  return photonCluster;
  
}


void ProtoCluster::AddParent(ProtoCluster* cluster){

  _parents.push_back(cluster);
  _connected = true;
  int relID = cluster->getEarliestRelationID();
  ProtoCluster* rel = cluster->getEarliestRelation();
  if(relID<_earliestRelationID){
    this->SetEarliestRelationID(relID);
    this->SetEarliestRelation(rel);
  }
  if(relID>_earliestRelationID){
    cluster->SetEarliestRelationID(_earliestRelationID);
    cluster->SetEarliestRelation(_earliestRelation);
  }
}


void ProtoCluster::SetEarliestRelationID(int id){

  
  if(_earliestRelationID==id)return;
  _earliestRelationID=id;
  for(unsigned int i=0;i<_daughters.size();++i)_daughters[i]->SetEarliestRelationID(id);
  for(unsigned int i=0;i<_parents.size();++i)_parents[i]->SetEarliestRelationID(id);
  
}


void ProtoCluster::SetEarliestRelation(ProtoCluster* pC){

  
  if(_earliestRelation==pC)return;
  _earliestRelation=pC;
  for(unsigned int i=0;i<_daughters.size();++i)_daughters[i]->SetEarliestRelation(pC);
  for(unsigned int i=0;i<_parents.size();++i)_parents[i]->SetEarliestRelation(pC);
  
}


float ProtoCluster::genericDistance(ProtoCluster* cluster){
  
  float retVal = 9999999.;

  for(int ilayer=_firstLayer;ilayer<=_lastLayer;++ilayer){
    float d = genericDistance(cluster, ilayer);
    if(d<retVal)retVal=d;
  }
  
  return retVal;

}



float ProtoCluster::genericDistance(ProtoCluster* cluster, int is, int ie){
  
  float retVal = 9999999.;

  int istart = is;
  if(istart<_firstLayer)istart=_firstLayer;
  int iend = ie;
  if(iend>_lastLayer)iend=_lastLayer;

  for(int ilayer=istart;ilayer<=iend;++ilayer){
    float d = genericDistance(cluster, ilayer);
    if(d<retVal)retVal=d;
  }
  
  return retVal;

}



float ProtoCluster::genericDistance(ProtoCluster* cluster, int ilayer){

  float retVal = 9999999.;
  
  for(unsigned int ihit=0;ihit<_hitsByPseudoLayer[ilayer].size();++ihit){
    MyCaloHitExtended* hiti = _hitsByPseudoLayer[ilayer][ihit];
    float xi = hiti->getPosition()[0];
    float yi = hiti->getPosition()[1];
    float zi = hiti->getPosition()[2];
    float ri = sqrt(xi*xi+yi*yi+zi*zi);
    float dir[3];
    dir[0] = xi/ri;
    dir[1] = yi/ri;
    dir[2] = zi/ri;
    int istart = ilayer-2;
    if(istart<0)istart=0;
    int iend   = ilayer+2;
    if(iend>MAX_NUMBER_OF_LAYERS)iend=MAX_NUMBER_OF_LAYERS;
    for(int jlayer=istart;jlayer<=iend;++jlayer){
      int nj = cluster->hitsInLayer(jlayer);
      for(int jhit=0;jhit<nj;++jhit){
	const MyCaloHitExtended* hitj = cluster->hitInLayer(jlayer,jhit);
	float xj = hitj->getPosition()[0];
	float yj = hitj->getPosition()[1];
	float zj = hitj->getPosition()[2];
	float dx = xi-xj;
	float dy = yi-yj;
	float dz = zi-zj;
	float ddx = dir[1]*dz - dir[2]*dy;
	float ddy = dir[2]*dx - dir[0]*dz;
	float ddz = dir[0]*dy - dir[1]*dx;      
	float distSq =  ddx*ddx + ddy*ddy+ ddz*ddz;
	float dPerp  = sqrt(distSq);  
	// Scalar product
	float dAlong = fabs(dir[0]*dx + dir[1]*dy + dir[2]*dz);
	float  dr = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
	dr = sqrt(dr);
	if(dAlong<1000 && dPerp<retVal)retVal=dPerp;
      }
    }
  }
  
  return retVal;
}



float ProtoCluster::FractionInCone(ProtoCluster* pCluster, float coneCosineHalfAngle){

  // look in cone starting from the point where cluster looks shower like
  //  29/05/08 Fix major bug in the calculation of the shower start point 
  // rsi = r0i + [(rs-r0).n]i   => rs = r0 + (rs-r0).n n

  float retVal =0.;
  bool useMip = true;

  if(_defunct)return retVal;
  if(_mipFitResult.ok==false){
    return retVal;
  }

  int nhits = pCluster->Hits();
  if(nhits<=0)return retVal;
  if(pCluster->FirstLayer()<_showerLayer)return retVal;
  ProtoCluster* pOrigin = pCluster->getEarliestRelation();
  if(pOrigin->FirstLayer()<_showerLayer)return retVal;
  //  if(_showerLayer==_lastLayer)return retVal;

  float dxi = _centroidX[_showerLayer] - _mipFitResult.intercept[0];
  float dyi = _centroidY[_showerLayer] - _mipFitResult.intercept[1];
  float dzi = _centroidZ[_showerLayer] - _mipFitResult.intercept[2];
  float dot = dxi*_mipFitResult.dir[0]+dyi*_mipFitResult.dir[1]+dzi*_mipFitResult.dir[2];


  float xi  = _mipFitResult.intercept[0] + dot*_mipFitResult.dir[0];
  float yi  = _mipFitResult.intercept[1] + dot*_mipFitResult.dir[1];
  float zi  = _mipFitResult.intercept[2] + dot*_mipFitResult.dir[2];

  float ddx = _centroidX[_showerLayer] - _centroidX[_firstLayer];
  float ddy = _centroidY[_showerLayer] - _centroidY[_firstLayer];
  float ddz = _centroidZ[_showerLayer] - _centroidZ[_firstLayer];
  float ddr = sqrt(ddx*ddx+ddy*ddy+ddz*ddz);
  ddx = ddx/ddr;
  ddy = ddy/ddr;
  ddz = ddz/ddr;
  float dcosCheck = ddx*_mipFitResult.dir[0]+ddy*_mipFitResult.dir[1]+ddz*_mipFitResult.dir[2];


  float ri = sqrt(xi*xi+yi*yi+zi*zi);

  float assoc = 0;
  float drmin= 99999.;
  for(int ilayer=pCluster->FirstLayer();ilayer<=pCluster->LastLayer();++ilayer){    for(int ihit = 0; ihit <pCluster->hitsInLayer(ilayer);++ihit){
      const MyCaloHitExtended* hiti = pCluster->hitInLayer(ilayer,ihit);
      float xj = hiti->getPosition()[0];
      float yj = hiti->getPosition()[1];
      float zj = hiti->getPosition()[2];
      float dx = xj-xi;
      float dy = yj-yi;
      float dz = zj-zi;
      float dr = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
      if(dr<drmin)drmin=dr;

      float cosTheta=0;
      if(useMip){
	cosTheta = (_mipFitResult.dir[0]*dx + _mipFitResult.dir[1]*dy + _mipFitResult.dir[2]*dz)/sqrt(dr);
      }else{
	cosTheta = (_initialDir[0]*dx + _initialDir[1]*dy + _initialDir[2]*dz)/sqrt(dr);
      }
      if(cosTheta>coneCosineHalfAngle)assoc++;
    }
  }
  float udotr=0;
  retVal = assoc/static_cast<float>(nhits);
  if(useMip){
    udotr = (xi*_mipFitResult.dir[0]+yi*_mipFitResult.dir[1]+zi*_mipFitResult.dir[2])/ri;
  }else{
    udotr = (xi*_initialDir[0]+yi*_initialDir[1]+zi*_initialDir[2])/ri;
  }

  // don't allow large distance associations at low angle
  if(udotr<0.5&&drmin>1000.0)retVal = 0.;
  if(udotr<0.75&&drmin>1500.0)retVal = 0.;
  if(udotr<0.25)retVal = 0.;

  //some debugging (turned off) 
  if(1==2){
    std::cout << "FL Centroid : " << _centroidX[_firstLayer] << "," << _centroidY[_firstLayer] << "," << _centroidZ[_firstLayer] << std::endl;
    std::cout << "SL Centroid : " << _centroidX[_showerLayer] << "," << _centroidY[_showerLayer] << "," << _centroidZ[_showerLayer] << std::endl;
    std::cout << "MIP Int : " << _mipFitResult.intercept[0] << "," << _mipFitResult.intercept[1]  << "," <<_mipFitResult.intercept[2]  << std::endl;    
    std::cout << "MIP Dir : " << _mipFitResult.dir[0] << "," << _mipFitResult.dir[1]  << "," <<_mipFitResult.dir[2]  << std::endl;    
    std::cout << "MIP FIT : " << xi << "," << yi << "," << zi << std::endl;    
    std::cout << "CHECK   : " << dcosCheck << std::endl;
  }


  return retVal;
}



float ProtoCluster::FractionInRadialCone(ProtoCluster* pCluster, float coneCosineHalfAngle){
 
  float retVal =0.;
  if(_defunct)return retVal;
  int nhits = pCluster->Hits();
  if(nhits<=0)return retVal;
  float xi = _centroidX[_showerLayer];
  float yi = _centroidY[_showerLayer];
  float zi = _centroidZ[_showerLayer];
  float ri = sqrt(xi*xi+yi*yi+zi*zi);
  float dir[3];
  dir[0] = xi/ri;
  dir[1] = yi/ri;
  dir[2] = zi/ri;

  if(_associatedTracks.size()>0){
    MyTrackExtended* track = _associatedTracks[0];
    xi =track->getSeedPosition()[0];
    yi =track->getSeedPosition()[1];
    zi =track->getSeedPosition()[2];
    dir[0]= track->getSeedDirection()[0];
    dir[1]= track->getSeedDirection()[1];
    dir[2]= track->getSeedDirection()[2];
  }
  

  float assoc = 0;
  float drmin= 99999.;
  for(int ilayer=pCluster->FirstLayer();ilayer<=pCluster->LastLayer();++ilayer){
    for(int ihit = 0; ihit <pCluster->hitsInLayer(ilayer);++ihit){
      const MyCaloHitExtended* hiti = pCluster->hitInLayer(ilayer,ihit);
      float xj = hiti->getPosition()[0];
      float yj = hiti->getPosition()[1];
      float zj = hiti->getPosition()[2];
      float dx = xj-xi;
      float dy = yj-yi;
      float dz = zj-zi;
      float dr = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
      if(dr<drmin)drmin=dr;
      float cosTheta = (dir[0]*dx + dir[1]*dy + dir[2]*dz)/sqrt(dr);
      if(cosTheta>coneCosineHalfAngle)assoc++;
    }
  }
  retVal = assoc/static_cast<float>(nhits);

  // don't allow large distance associations 
  //  if(drmin>1000.0)retVal = 0.;

  return retVal;
}



float ProtoCluster::FractionInTrackCone(MyTrackExtended* track, float coneCosineHalfAngle){
 
  float retVal =0.;
  if(_defunct)return retVal;
  int nhits = this->Hits();
  if(nhits<=0)return retVal;

  float dir[3];
  float xi =track->getSeedPosition()[0];
  float yi =track->getSeedPosition()[1];
  float zi =track->getSeedPosition()[2];
  dir[0]= track->getSeedDirection()[0];
  dir[1]= track->getSeedDirection()[1];
  dir[2]= track->getSeedDirection()[2];

 
  float assoc = 0;
  float drmin= 99999.;
  for(int ilayer=this->FirstLayer();ilayer<=this->LastLayer();++ilayer){
    for(int ihit = 0; ihit <this->hitsInLayer(ilayer);++ihit){
      const MyCaloHitExtended* hiti = this->hitInLayer(ilayer,ihit);
      float xj = hiti->getPosition()[0];
      float yj = hiti->getPosition()[1];
      float zj = hiti->getPosition()[2];
      float dx = xj-xi;
      float dy = yj-yi;
      float dz = zj-zi;
      float dr = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
      if(dr<drmin)drmin=dr;
      float cosTheta = (dir[0]*dx + dir[1]*dy + dir[2]*dz)/sqrt(dr);
      if(cosTheta>coneCosineHalfAngle)assoc++;
    }
  }
  retVal = assoc/static_cast<float>(nhits);

  return retVal;
}