/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/RecoParticleFlow/PFClusterProducer/src/PFClusterAlgo.cc

Go to the documentation of this file.

#include "RecoParticleFlow/PFClusterProducer/interface/PFClusterAlgo.h"
#include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "Math/GenVector/VectorUtil.h"
#include "TFile.h"
#include "TH2F.h"
#include "TROOT.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <stdexcept>
#include <string>
#include <sstream>

using namespace std;

unsigned PFClusterAlgo::prodNum_ = 1;

//for debug only 
//#define PFLOW_DEBUG

PFClusterAlgo::PFClusterAlgo() :
  pfClusters_( new vector<reco::PFCluster> ),
  pfRecHitsCleaned_( new vector<reco::PFRecHit> ),
  threshBarrel_(0.),
  threshPtBarrel_(0.),
  threshSeedBarrel_(0.2),
  threshPtSeedBarrel_(0.),
  threshEndcap_(0.),
  threshPtEndcap_(0.),
  threshSeedEndcap_(0.6),
  threshPtSeedEndcap_(0.),
  threshCleanBarrel_(1E5),
  minS4S1Barrel_(0.),
  threshDoubleSpikeBarrel_(1E9),
  minS6S2DoubleSpikeBarrel_(-1.),
  threshCleanEndcap_(1E5),
  minS4S1Endcap_(0.),
  threshDoubleSpikeEndcap_(1E9),
  minS6S2DoubleSpikeEndcap_(-1.),
  nNeighbours_(4),
  posCalcNCrystal_(-1),
  posCalcP1_(-1),
  showerSigma_(5),
  useCornerCells_(false),
  cleanRBXandHPDs_(false),
  debug_(false) 
{
  file_ = 0;
  hBNeighbour = 0;
  hENeighbour = 0;
}

void
PFClusterAlgo::write() { 

  if ( file_ ) { 
    file_->Write();
    cout << "Benchmark output written to file " << file_->GetName() << endl;
    file_->Close();
  }

}


void PFClusterAlgo::doClustering( const PFRecHitHandle& rechitsHandle ) {
  const reco::PFRecHitCollection& rechits = * rechitsHandle;

  // cache the Handle to the rechits
  rechitsHandle_ = rechitsHandle;

  // clear rechits mask
  mask_.clear();
  mask_.resize( rechits.size(), true );

  // perform clustering
  doClusteringWorker( rechits );
}

void PFClusterAlgo::doClustering( const PFRecHitHandle& rechitsHandle, const std::vector<bool> & mask ) {

  const reco::PFRecHitCollection& rechits = * rechitsHandle;

  // cache the Handle to the rechits
  rechitsHandle_ = rechitsHandle;

  // use the specified mask, unless it doesn't match with the rechits
  mask_.clear();
  if (mask.size() == rechits.size()) {
      mask_.insert( mask_.end(), mask.begin(), mask.end() );
  } else {
      edm::LogError("PClusterAlgo::doClustering") << "map size should be " << rechits.size() << ". Will be reinitialized.";
      mask_.resize( rechits.size(), true );
  }

  // perform clustering

  doClusteringWorker( rechits );

}

void PFClusterAlgo::doClustering( const reco::PFRecHitCollection& rechits ) {

  // using rechits without a Handle, clear to avoid a stale member
  rechitsHandle_.clear();

  // clear rechits mask
  mask_.clear();
  mask_.resize( rechits.size(), true );

  // perform clustering
  doClusteringWorker( rechits );

}

void PFClusterAlgo::doClustering( const reco::PFRecHitCollection& rechits, const std::vector<bool> & mask ) {
  // using rechits without a Handle, clear to avoid a stale member

  rechitsHandle_.clear();

  // use the specified mask, unless it doesn't match with the rechits
  mask_.clear();

  if (mask.size() == rechits.size()) {
      mask_.insert( mask_.end(), mask.begin(), mask.end() );
  } else {
      edm::LogError("PClusterAlgo::doClustering") << "map size should be " << rechits.size() << ". Will be reinitialized.";
      mask_.resize( rechits.size(), true );
  }

  // perform clustering
  doClusteringWorker( rechits );

}


void PFClusterAlgo::doClusteringWorker( const reco::PFRecHitCollection& rechits ) {


  if ( pfClusters_.get() )
    pfClusters_->clear();
  else 
    pfClusters_.reset( new std::vector<reco::PFCluster> );


  if ( pfRecHitsCleaned_.get() )
    pfRecHitsCleaned_->clear();
  else 
    pfRecHitsCleaned_.reset( new std::vector<reco::PFRecHit> );


  eRecHits_.clear();

  for ( unsigned i = 0; i < rechits.size(); i++ ){

    eRecHits_.insert( make_pair( rechit(i, rechits).energy(), i) );
  }

  color_.clear(); 
  color_.resize( rechits.size(), 0 );

  seedStates_.clear();
  seedStates_.resize( rechits.size(), UNKNOWN );

  usedInTopo_.clear();
  usedInTopo_.resize( rechits.size(), false );

  if ( cleanRBXandHPDs_ ) cleanRBXAndHPD( rechits);

  // look for seeds.

  findSeeds( rechits );

  // build topological clusters around seeds
  buildTopoClusters( rechits );

  // look for PFClusters inside each topological cluster (one per seed)
  
  
  //  int ix=0;
  //  for (reco::PFRecHitCollection::const_iterator cand =rechits.begin(); cand<rechits.end(); cand++){
  //    cout <<ix++<<" "<< cand->layer()<<endl;
  //  }


  for(unsigned i=0; i<topoClusters_.size(); i++) {

    const std::vector< unsigned >& topocluster = topoClusters_[i];

    buildPFClusters( topocluster, rechits ); 

  }

}


double PFClusterAlgo::parameter( Parameter paramtype, 
                                 PFLayer::Layer layer,
                                 unsigned iCoeff, int iring )  const {
  

  double value = 0;

  switch( layer ) {

  case PFLayer::ECAL_BARREL:
  case PFLayer::HCAL_BARREL1:
  case PFLayer::HCAL_BARREL2: //  HO. 
                                 
    switch(paramtype) {
    case THRESH:
      value = (iring==0) ? threshBarrel_ : threshEndcap_; //For HO Ring0 and others
      break;
    case SEED_THRESH:
      value = (iring==0) ? threshSeedBarrel_ : threshSeedEndcap_;
      break;
    case PT_THRESH:
      value = (iring==0) ? threshPtBarrel_ : threshPtEndcap_;
      break;
    case SEED_PT_THRESH:
      value = (iring==0) ? threshPtSeedBarrel_ : threshPtSeedEndcap_;
      break;
    case CLEAN_THRESH:
      value = threshCleanBarrel_;
      break;
    case CLEAN_S4S1:
      value = minS4S1Barrel_[iCoeff];
      break;
    case DOUBLESPIKE_THRESH:
      value = threshDoubleSpikeBarrel_;
      break;
    case DOUBLESPIKE_S6S2:
      value = minS6S2DoubleSpikeBarrel_;
      break;
    default:
      cerr<<"PFClusterAlgo::parameter : unknown parameter type "
          <<paramtype<<endl;
      assert(0);
    }
    break;
  case PFLayer::ECAL_ENDCAP:
  case PFLayer::HCAL_ENDCAP:
  case PFLayer::PS1:
  case PFLayer::PS2:
  case PFLayer::HF_EM:
  case PFLayer::HF_HAD:
    // and no particle flow in VFCAL
    switch(paramtype) {
    case THRESH:
      value = threshEndcap_;
      break;
    case SEED_THRESH:
      value = threshSeedEndcap_;
      break;
    case PT_THRESH:
      value = threshPtEndcap_;
      break;
    case SEED_PT_THRESH:
      value = threshPtSeedEndcap_;
      break;
    case CLEAN_THRESH:
      value = threshCleanEndcap_;
      break;
    case CLEAN_S4S1:
      value = minS4S1Endcap_[iCoeff];
      break;
    case DOUBLESPIKE_THRESH:
      value = threshDoubleSpikeEndcap_;
      break;
    case DOUBLESPIKE_S6S2:
      value = minS6S2DoubleSpikeEndcap_;
      break;
    default:
      cerr<<"PFClusterAlgo::parameter : unknown parameter type "
          <<paramtype<<endl;
      assert(0);
    }
    break;
  default:
    cerr<<"PFClusterAlgo::parameter : unknown layer "<<layer<<endl;
    assert(0);
    break;
  }

  return value;
}


void 
PFClusterAlgo::cleanRBXAndHPD(  const reco::PFRecHitCollection& rechits ) {

  std::map< int, std::vector<unsigned> > hpds;
  std::map< int, std::vector<unsigned> > rbxs;

  for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {

    unsigned  rhi      = ih->second; 

    if(! masked(rhi) ) continue;
    // rechit was asked to be processed
    const reco::PFRecHit& rhit = rechit( rhi, rechits);
    //double energy = rhit.energy();
    int layer = rhit.layer();
    if ( layer != PFLayer::HCAL_BARREL1 &&
         layer != PFLayer::HCAL_ENDCAP ) break; 
      // layer != PFLayer::HCAL_BARREL2) break; //BARREL2 for HO : need specific cleaning
    HcalDetId theHcalDetId = HcalDetId(rhit.detId());
    int ieta = theHcalDetId.ieta();
    int iphi = theHcalDetId.iphi();
    int ihpd = ieta < 0 ?  
      ( layer == PFLayer::HCAL_ENDCAP ?  -(iphi+1)/2-100 : -iphi ) : 
      ( layer == PFLayer::HCAL_ENDCAP ?   (iphi+1)/2+100 :  iphi ) ;      
    hpds[ihpd].push_back(rhi);
    int irbx = ieta < 0 ? 
      ( layer == PFLayer::HCAL_ENDCAP ?  -(iphi+5)/4 - 20 : -(iphi+5)/4 ) : 
      ( layer == PFLayer::HCAL_ENDCAP ?   (iphi+5)/4 + 20 :  (iphi+5)/4 ) ;      
    if ( irbx == 19 ) irbx = 1;
    else if ( irbx == -19 ) irbx = -1;
    else if ( irbx == 39 ) irbx = 21;
    else if ( irbx == -39 ) irbx = -21;
    rbxs[irbx].push_back(rhi);
  }

  // Loop on readout boxes
  for ( std::map<int, std::vector<unsigned> >::iterator itrbx = rbxs.begin();
        itrbx != rbxs.end(); ++itrbx ) { 
    if ( ( abs(itrbx->first)<20 && itrbx->second.size() > 30 ) || 
         ( abs(itrbx->first)>20 && itrbx->second.size() > 30 ) ) { 
      const std::vector<unsigned>& rhits = itrbx->second;
      double totalEta = 0.;
      double totalEtaW = 0.;
      double totalPhi = 0.;
      double totalPhiW = 0.;
      double totalEta2 = 1E-9;
      double totalEta2W = 1E-9;
      double totalPhi2 = 1E-9;
      double totalPhi2W = 1E-9;
      double totalEnergy = 0.;
      double totalEnergy2 = 1E-9;
      unsigned nSeeds = rhits.size();
      unsigned nSeeds0 = rhits.size();
      std::map< int,std::vector<unsigned> > theHPDs;
      std::multimap< double,unsigned > theEnergies;
      for ( unsigned jh=0; jh < rhits.size(); ++jh ) {
        const reco::PFRecHit& hit = rechit(rhits[jh], rechits);
        // Check if the hit is a seed
        unsigned nN = 0;
        bool isASeed = true;
        const vector<unsigned>& neighbours4 = *(& hit.neighbours4());
        for(unsigned in=0; in<neighbours4.size(); in++) {
          const reco::PFRecHit& neighbour = rechit( neighbours4[in], rechits ); 
          // one neighbour has a higher energy -> the tested rechit is not a seed
          if( neighbour.energy() > hit.energy() ) {
            --nSeeds;
            --nSeeds0;
            isASeed = false;
            break;
          } else {
            if ( neighbour.energy() > 0.4 ) ++nN;
          }
        }
        if ( isASeed && !nN ) --nSeeds0;

        HcalDetId theHcalDetId = HcalDetId(hit.detId());
        // int ieta = theHcalDetId.ieta();
        int iphi = theHcalDetId.iphi();
        // std::cout << "Hit : " << hit.energy() << " " << ieta << " " << iphi << std::endl;
        if ( hit.layer() == PFLayer::HCAL_BARREL1 )
          theHPDs[iphi].push_back(rhits[jh]);
        else
          theHPDs[(iphi-1)/2].push_back(rhits[jh]);
        theEnergies.insert(std::pair<double,unsigned>(hit.energy(),rhits[jh]));
        totalEnergy += hit.energy();
        totalPhi += fabs(hit.position().phi());
        totalPhiW += hit.energy()*fabs(hit.position().phi());
        totalEta += hit.position().eta();
        totalEtaW += hit.energy()*hit.position().eta();
        totalEnergy2 += hit.energy()*hit.energy();
        totalPhi2 += hit.position().phi()*hit.position().phi();
        totalPhi2W += hit.energy()*hit.position().phi()*hit.position().phi();
        totalEta2 += hit.position().eta()*hit.position().eta();
        totalEta2W += hit.energy()*hit.position().eta()*hit.position().eta();
      }
      // totalPhi /= totalEnergy;
      totalPhi /= rhits.size();
      totalEta /= rhits.size();
      totalPhiW /= totalEnergy;
      totalEtaW /= totalEnergy;
      totalPhi2 /= rhits.size();
      totalEta2 /= rhits.size();
      totalPhi2W /= totalEnergy;
      totalEta2W /= totalEnergy;
      totalPhi2 = std::sqrt(totalPhi2 - totalPhi*totalPhi);
      totalEta2 = std::sqrt(totalEta2 - totalEta*totalEta);
      totalPhi2W = std::sqrt(totalPhi2W - totalPhiW*totalPhiW);
      totalEta2W = std::sqrt(totalEta2W - totalEtaW*totalEtaW);
      totalEnergy /= rhits.size();
      totalEnergy2 /= rhits.size();
      totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
      //if ( totalPhi2W/totalEta2W < 0.18 ) { 
      if ( nSeeds0 > 6 ) {
        unsigned nHPD15 = 0;
        for ( std::map<int, std::vector<unsigned> >::iterator itHPD = theHPDs.begin();
              itHPD != theHPDs.end(); ++itHPD ) { 
          int hpdN = itHPD->first;
          const std::vector<unsigned>& hpdHits = itHPD->second;
          if ( ( abs(hpdN) < 100 && hpdHits.size() > 14 ) || 
               ( abs(hpdN) > 100 && hpdHits.size() > 14 ) ) ++nHPD15;
        }
        if ( nHPD15 > 1 ) {
          /*
          std::cout << "Read out box numero " << itrbx->first 
                    << " has " << itrbx->second.size() << " hits in it !"
                    << std::endl << "sigma Eta/Phi = " << totalEta2 << " " << totalPhi2 << " " << totalPhi2/totalEta2
                    << std::endl << "sigma EtaW/PhiW = " << totalEta2W << " " << totalPhi2W << " " << totalPhi2W/totalEta2W
                    << std::endl << "E = " << totalEnergy << " +/- " << totalEnergy2
                    << std::endl << "nSeeds = " << nSeeds << " " << nSeeds0
                    << std::endl;
          for ( std::map<int, std::vector<unsigned> >::iterator itHPD = theHPDs.begin();
                itHPD != theHPDs.end(); ++itHPD ) { 
            unsigned hpdN = itHPD->first;
            const std::vector<unsigned>& hpdHits = itHPD->second;
            std::cout << "HPD number " << hpdN << " contains " << hpdHits.size() << " hits" << std::endl;
          }
          */
          std::multimap<double, unsigned >::iterator ntEn = theEnergies.end();
          --ntEn;
          unsigned nn = 0;
          double threshold = 1.;
          for ( std::multimap<double, unsigned >::iterator itEn = theEnergies.begin();
                itEn != theEnergies.end(); ++itEn ) {
            ++nn;
            if ( nn < 5 ) { 
              mask_[itEn->second] = false;
            } else if ( nn == 5 ) { 
              threshold = itEn->first * 5.;
              mask_[itEn->second] = false;
            } else { 
              if ( itEn->first < threshold ) mask_[itEn->second] = false;
            }
            if ( !masked(itEn->second) ) { 
              reco::PFRecHit theCleanedHit(rechit(itEn->second, rechits));
              //theCleanedHit.setRescale(0.);
              pfRecHitsCleaned_->push_back(theCleanedHit);
            }
            /*
            if ( !masked(itEn->second) ) 
              std::cout << "Hit Energies = " << itEn->first 
                        << " " << ntEn->first/itEn->first << " masked " << std::endl;
            else 
              std::cout << "Hit Energies = " << itEn->first 
                        << " " << ntEn->first/itEn->first << " kept for clustering " << std::endl;
            */
          }
        }
      }
    }
  }

  // Loop on hpd's
  std::map<int, std::vector<unsigned> >::iterator neighbour1;
  std::map<int, std::vector<unsigned> >::iterator neighbour2;
  std::map<int, std::vector<unsigned> >::iterator neighbour0;
  std::map<int, std::vector<unsigned> >::iterator neighbour3;
  unsigned size1 = 0;
  unsigned size2 = 0;
  for ( std::map<int, std::vector<unsigned> >::iterator ithpd = hpds.begin();
        ithpd != hpds.end(); ++ithpd ) { 

    const std::vector<unsigned>& rhits = ithpd->second;
    std::multimap< double,unsigned > theEnergies;
    double totalEnergy = 0.;
    double totalEnergy2 = 1E-9;
    for ( unsigned jh=0; jh < rhits.size(); ++jh ) {
      const reco::PFRecHit& hit = rechit(rhits[jh], rechits);
      totalEnergy += hit.energy();
      totalEnergy2 += hit.energy()*hit.energy();
      theEnergies.insert(std::pair<double,unsigned>(hit.energy(),rhits[jh]));
    }
    totalEnergy /= rhits.size();
    totalEnergy2 /= rhits.size();
    totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);

    if ( ithpd->first == 1 ) neighbour1 = hpds.find(72);
    else if ( ithpd->first == -1 ) neighbour1 = hpds.find(-72);
    else if ( ithpd->first == 101 ) neighbour1 = hpds.find(136);
    else if ( ithpd->first == -101 ) neighbour1 = hpds.find(-136);
    else neighbour1 = ithpd->first > 0 ? hpds.find(ithpd->first-1) : hpds.find(ithpd->first+1) ;

    if ( ithpd->first == 72 ) neighbour2 = hpds.find(1);
    else if ( ithpd->first == -72 ) neighbour2 = hpds.find(-1);
    else if ( ithpd->first == 136 ) neighbour2 = hpds.find(101);
    else if ( ithpd->first == -136 ) neighbour2 = hpds.find(-101);
    else neighbour2 = ithpd->first > 0 ? hpds.find(ithpd->first+1) : hpds.find(ithpd->first-1) ;

    if ( neighbour1 != hpds.end() ) { 
      if ( neighbour1->first == 1 ) neighbour0 = hpds.find(72);
      else if ( neighbour1->first == -1 ) neighbour0 = hpds.find(-72);
      else if ( neighbour1->first == 101 ) neighbour0 = hpds.find(136);
      else if ( neighbour1->first == -101 ) neighbour0 = hpds.find(-136);
      else neighbour0 = neighbour1->first > 0 ? hpds.find(neighbour1->first-1) : hpds.find(neighbour1->first+1) ;
    } 

    if ( neighbour2 != hpds.end() ) { 
      if ( neighbour2->first == 72 ) neighbour3 = hpds.find(1);
      else if ( neighbour2->first == -72 ) neighbour3 = hpds.find(-1);
      else if ( neighbour2->first == 136 ) neighbour3 = hpds.find(101);
      else if ( neighbour2->first == -136 ) neighbour3 = hpds.find(-101);
      else neighbour3 = neighbour2->first > 0 ? hpds.find(neighbour2->first+1) : hpds.find(neighbour2->first-1) ;
    }

    size1 = neighbour1 != hpds.end() ? neighbour1->second.size() : 0;
    size2 = neighbour2 != hpds.end() ? neighbour2->second.size() : 0;

    // Also treat the case of two neighbouring HPD's not in the same RBX
    if ( size1 > 10 ) { 
      if ( ( abs(neighbour1->first) > 100 && neighbour1->second.size() > 15 ) || 
           ( abs(neighbour1->first) < 100 && neighbour1->second.size() > 12 ) ) 
        size1 = neighbour0 != hpds.end() ? neighbour0->second.size() : 0;
    }
    if ( size2 > 10 ) { 
      if ( ( abs(neighbour2->first) > 100 && neighbour2->second.size() > 15 ) || 
           ( abs(neighbour2->first) < 100 && neighbour2->second.size() > 12 ) ) 
        size2 = neighbour3 != hpds.end() ? neighbour3->second.size() : 0;
    }
    
    if ( ( abs(ithpd->first) > 100 && ithpd->second.size() > 15 ) || 
         ( abs(ithpd->first) < 100 && ithpd->second.size() > 12 ) )
      if ( (float)(size1 + size2)/(float)ithpd->second.size() < 1.0 ) {
        /* 
        std::cout << "HPD numero " << ithpd->first 
                  << " has " << ithpd->second.size() << " hits in it !" << std::endl
                  << "Neighbours : " << size1 << " " << size2
                  << std::endl;
        */
        std::multimap<double, unsigned >::iterator ntEn = theEnergies.end();
        --ntEn;
        unsigned nn = 0;
        double threshold = 1.;
        for ( std::multimap<double, unsigned >::iterator itEn = theEnergies.begin();
              itEn != theEnergies.end(); ++itEn ) {
          ++nn;
          if ( nn < 5 ) { 
            mask_[itEn->second] = false;
          } else if ( nn == 5 ) { 
            threshold = itEn->first * 2.5;
            mask_[itEn->second] = false;
          } else { 
            if ( itEn->first < threshold ) mask_[itEn->second] = false;
          }
          if ( !masked(itEn->second) ) { 
            reco::PFRecHit theCleanedHit(rechit(itEn->second, rechits));
            //theCleanedHit.setRescale(0.);
            pfRecHitsCleaned_->push_back(theCleanedHit);
          }
          /*
          if ( !masked(itEn->second) ) 
            std::cout << "Hit Energies = " << itEn->first 
                      << " " << ntEn->first/itEn->first << " masked " << std::endl;
          else 
            std::cout << "Hit Energies = " << itEn->first 
            << " " << ntEn->first/itEn->first << " kept for clustering " << std::endl;
          */
        }
      }
  }

}


void PFClusterAlgo::findSeeds( const reco::PFRecHitCollection& rechits ) {

  seeds_.clear();

  // should replace this by the message logger.
#ifdef PFLOW_DEBUG
  if(debug_) 
    cout<<"PFClusterAlgo::findSeeds : start"<<endl;
#endif


  // An empty list of neighbours
  const vector<unsigned> noNeighbours(0, static_cast<unsigned>(0));

  // loop on rechits (sorted by decreasing energy - not E_T)

  for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {

    unsigned  rhi      = ih->second; 

    if(! masked(rhi) ) continue;
    // rechit was asked to be processed

    double    rhenergy = ih->first;   
    const reco::PFRecHit& wannaBeSeed = rechit(rhi, rechits);
     
    if( seedStates_[rhi] == NO ) continue;
    // this hit was already tested, and is not a seed
 
    // determine seed energy threshold depending on the detector
    int layer = wannaBeSeed.layer();
    //for HO Ring0 and 1/2 boundary
    
    int iring = 0;
    if (layer==PFLayer::HCAL_BARREL2 && abs(wannaBeSeed.positionREP().Eta())>0.34) iring= 1;

    double seedThresh = parameter( SEED_THRESH, 
                                   static_cast<PFLayer::Layer>(layer), 0, iring );

    double seedPtThresh = parameter( SEED_PT_THRESH, 
                                     static_cast<PFLayer::Layer>(layer), 0., iring );

    double cleanThresh = parameter( CLEAN_THRESH, 
                                    static_cast<PFLayer::Layer>(layer), 0, iring );

    double minS4S1_a = parameter( CLEAN_S4S1, 
                                  static_cast<PFLayer::Layer>(layer), 0, iring );

    double minS4S1_b = parameter( CLEAN_S4S1, 
                                  static_cast<PFLayer::Layer>(layer), 1, iring );

    double doubleSpikeThresh = parameter( DOUBLESPIKE_THRESH, 
                                          static_cast<PFLayer::Layer>(layer), 0, iring );

    double doubleSpikeS6S2 = parameter( DOUBLESPIKE_S6S2, 
                                        static_cast<PFLayer::Layer>(layer), 0, iring );

#ifdef PFLOW_DEBUG
    if(debug_) 
      cout<<"layer:"<<layer<<" seedThresh:"<<seedThresh<<endl;
#endif


    if( rhenergy < seedThresh || (seedPtThresh>0. && wannaBeSeed.pt2() < seedPtThresh*seedPtThresh )) {
      seedStates_[rhi] = NO; 
      continue;
    } 

    // Find the cell unused neighbours
    const vector<unsigned>* nbp;
    double tighterE = 1.0;
    double tighterF = 1.0;

    switch ( layer ) { 
    case PFLayer::ECAL_BARREL:         
    case PFLayer::ECAL_ENDCAP:       
    case PFLayer::HCAL_BARREL1:
    case PFLayer::HCAL_BARREL2:
    case PFLayer::HCAL_ENDCAP:
      tighterE = 2.0;
      tighterF = 3.0;
    case PFLayer::HF_EM:
    case PFLayer::HF_HAD:
      if( nNeighbours_ == 4 ) {
        nbp = & wannaBeSeed.neighbours4();
      }
      else if( nNeighbours_ == 8 ) {
        nbp = & wannaBeSeed.neighbours8();
      }
      else if( nNeighbours_ == 0 ) {
        nbp = & noNeighbours;
        // Allows for no clustering at all: all rechits are clusters.
        // Useful for HF
      }
      else {
        cerr<<"you're not allowed to set n neighbours to "
            <<nNeighbours_<<endl;
        assert(0);
      }
      break;
    case PFLayer::PS1:       
    case PFLayer::PS2:     
      nbp = & wannaBeSeed.neighbours4();
      break;

    default:
      cerr<<"CellsEF::PhotonSeeds : unknown layer "<<layer<<endl;
      assert(0);
    }

    const vector<unsigned>& neighbours = *nbp;

      
    // Select as a seed if all neighbours have a smaller energy

    seedStates_[rhi] = YES;
    for(unsigned in=0; in<neighbours.size(); in++) {
        
      unsigned rhj =  neighbours[in];
      // Ignore neighbours already masked
      if ( !masked(rhj) ) continue;
      const reco::PFRecHit& neighbour = rechit( rhj, rechits ); 
        
      // one neighbour has a higher energy -> the tested rechit is not a seed
      if( neighbour.energy() > wannaBeSeed.energy() ) {
        seedStates_[rhi] = NO;
        break;
      }
    }
    
    // Cleaning : check energetic, isolated seeds, likely to come from erratic noise.
    if ( file_ || wannaBeSeed.energy() > cleanThresh ) { 
      
      const vector<unsigned>& neighbours4 = *(& wannaBeSeed.neighbours4());
      // Determine the fraction of surrounding energy
      double surroundingEnergy = wannaBeSeed.energyUp();
      double neighbourEnergy = 0.;
      double layerEnergy = 0.;
      for(unsigned in4=0; in4<neighbours4.size(); in4++) {
        unsigned rhj =  neighbours4[in4];
        // Ignore neighbours already masked
        if ( !masked(rhj) ) continue;
        const reco::PFRecHit& neighbour = rechit( rhj, rechits ); 
        surroundingEnergy += neighbour.energy() + neighbour.energyUp();
        neighbourEnergy += neighbour.energy() + neighbour.energyUp();
        layerEnergy += neighbour.energy();
      }
      // Fraction 0 is the balance between EM and HAD layer for this tower
      // double fraction0 = layer == PFLayer::HF_EM || layer == PFLayer::HF_HAD ? 
      //   wannaBeSeed.energyUp()/wannaBeSeed.energy() : 1.;
      // Fraction 1 is the balance between the hit and its neighbours from both layers
      double fraction1 = surroundingEnergy/wannaBeSeed.energy();
      // Fraction 2 is the balance between the tower and the tower neighbours
      // double fraction2 = neighbourEnergy/(wannaBeSeed.energy()+wannaBeSeed.energyUp());
      // Fraction 3 is the balance between the hits and the hits neighbours in the same layer.
      // double fraction3 = layerEnergy/(wannaBeSeed.energy());
      // Mask the seed and the hit if energetic/isolated rechit
      // if ( fraction0 < minS4S1 || fraction1 < minS4S1 || fraction2 < minS4S1 || fraction3 < minS4S1 ) {
      // if ( fraction1 < minS4S1 || ( wannaBeSeed.energy() > 1.5*cleanThresh && fraction0 + fraction3 < minS4S1 ) ) {
      
      if ( file_ ) { 
        if ( layer == PFLayer::ECAL_BARREL || layer == PFLayer::HCAL_BARREL1 || layer == PFLayer::HCAL_BARREL2) { //BARREL2 for HO 
          /*
          double eta = wannaBeSeed.position().eta();
          double phi = wannaBeSeed.position().phi();
          std::pair<double,double> dcr = dCrack(phi,eta);
          double dcrmin = std::min(dcr.first, dcr.second);
          if ( dcrmin > 1. ) 
          */
          hBNeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
        } else if ( fabs(wannaBeSeed.position().eta()) < 5.0  ) {
          if ( layer == PFLayer::ECAL_ENDCAP || layer == PFLayer::HCAL_ENDCAP ) { 
            if ( wannaBeSeed.energy() > 1000 ) { 
              if ( fabs(wannaBeSeed.position().eta()) < 2.8  ) 
                hENeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
            }
          } else { 
            hENeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
          }
        }
      }
      
      if ( wannaBeSeed.energy() > cleanThresh ) { 
        double f1Cut = minS4S1_a * log10(wannaBeSeed.energy()) + minS4S1_b;
        if ( fraction1 < f1Cut ) {
          // Double the energy cleaning threshold when close to the ECAL/HCAL - HF transition
          double eta = wannaBeSeed.position().eta();
          double phi = wannaBeSeed.position().phi();
          std::pair<double,double> dcr = dCrack(phi,eta);
          double dcrmin = layer == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second;
          eta = fabs(eta);
          if (   eta < 5.0 &&                         // No cleaning for the HF border 
                 ( ( eta < 2.85 && dcrmin > 1. ) || 
                   ( rhenergy > tighterE*cleanThresh && 
                     fraction1 < f1Cut/tighterF ) )  // Tighter cleaning for various cracks 
              ) { 
            seedStates_[rhi] = CLEAN;
            mask_[rhi] = false;
            reco::PFRecHit theCleanedHit(wannaBeSeed);
            //theCleanedHit.setRescale(0.);
            pfRecHitsCleaned_->push_back(theCleanedHit);
            /*
            std::cout << "A seed with E/pT/eta/phi = " << wannaBeSeed.energy() << " " << wannaBeSeed.energyUp() 
                      << " " << sqrt(wannaBeSeed.pt2()) << " " << wannaBeSeed.position().eta() << " " << phi 
                      << " and with surrounding fractions = " << fraction0 << " " << fraction1 
                      << " " << fraction2 << " " << fraction3
                      << " in layer " << layer 
                      << " had been cleaned " << std::endl
                      << " Distances to cracks : " << dcr.first << " " << dcr.second << std::endl
                      << "(Cuts were : " << cleanThresh << " and " << f1Cut << ")" << std::endl; 
            */
          }
        }
      }
    }

    // Clean double spikes
    if ( mask_[rhi] && wannaBeSeed.energy() > doubleSpikeThresh ) {
      // Determine energy surrounding the seed and the most energetic neighbour
      double surroundingEnergyi = 0.;
      double enmax = -999.;
      unsigned mostEnergeticNeighbour = 0;
      const vector<unsigned>& neighbours4i = *(& wannaBeSeed.neighbours4());
      for(unsigned in4=0; in4<neighbours4i.size(); in4++) {
        unsigned rhj =  neighbours4i[in4];
        if ( !masked(rhj) ) continue;
        const reco::PFRecHit& neighbouri = rechit( rhj, rechits );
        surroundingEnergyi += neighbouri.energy();
        if ( neighbouri.energy() > enmax ) { 
          enmax = neighbouri.energy();
          mostEnergeticNeighbour = rhj;
        }
      }
      // Is there an energetic neighbour ?
      if ( enmax > 0. ) { 
        unsigned rhj = mostEnergeticNeighbour;
        const reco::PFRecHit& neighbouri = rechit( rhj, rechits );
        double surroundingEnergyj = 0.;
        //if ( mask_[rhj] && neighbouri.energy() > doubleSpikeThresh ) {
        // Determine energy surrounding the energetic neighbour
        const vector<unsigned>& neighbours4j = *(& neighbouri.neighbours4());
        for(unsigned jn4=0; jn4<neighbours4j.size(); jn4++) {
          unsigned rhk =  neighbours4j[jn4];
          const reco::PFRecHit& neighbourj = rechit( rhk, rechits ); 
          surroundingEnergyj += neighbourj.energy();
        }
        // The energy surrounding the double spike candidate 
        double surroundingEnergyFraction = 
          (surroundingEnergyi+surroundingEnergyj) / (wannaBeSeed.energy()+neighbouri.energy()) - 1.;
        if ( surroundingEnergyFraction < doubleSpikeS6S2 ) { 
          double eta = wannaBeSeed.position().eta();
          double phi = wannaBeSeed.position().phi();
          std::pair<double,double> dcr = dCrack(phi,eta);
          double dcrmin = layer == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second;
          eta = fabs(eta);
          if (  ( eta < 5.0 && dcrmin > 1. ) ||
                ( wannaBeSeed.energy() > tighterE*doubleSpikeThresh &&
                  surroundingEnergyFraction < doubleSpikeS6S2/tighterF ) ) {
            /*
            std::cout << "Double spike cleaned : Energies = " << wannaBeSeed.energy()
                      << " " << neighbouri.energy() 
                      << " surrounded by only " << surroundingEnergyFraction*100.
                      << "% of the two spike energies at eta/phi/distance to closest crack = "  
                      <<        eta << " " << phi << " " << dcrmin
                      << std::endl;
            */
            // mask the seed
            seedStates_[rhi] = CLEAN;
            mask_[rhi] = false;
            reco::PFRecHit theCleanedSeed(wannaBeSeed);
            pfRecHitsCleaned_->push_back(theCleanedSeed);
            // mask the neighbour
            seedStates_[rhj] = CLEAN;
            mask_[rhj] = false;
            reco::PFRecHit theCleanedNeighbour(wannaBeSeed);
            pfRecHitsCleaned_->push_back(neighbouri);
          }
        }
      } else { 
        /*
        std::cout << "PFClusterAlgo : Double spike search : An isolated cell should have been killed " << std::endl 
                  << "but is going through the double spike search!" << std::endl;
        */
      }
    }

    if ( seedStates_[rhi] == YES ) {

      // seeds_ contains the indices of all seeds. 
      seeds_.push_back( rhi );
      
      // marking the rechit
      paint(rhi, SEED);
        
      // then all neighbours cannot be seeds and are flagged as such
      for(unsigned in=0; in<neighbours.size(); in++) {
        seedStates_[ neighbours[in] ] = NO;
      }
    }

  }  

#ifdef PFLOW_DEBUG
  if(debug_) 
    cout<<"PFClusterAlgo::findSeeds : done"<<endl;
#endif
}



  
void PFClusterAlgo::buildTopoClusters( const reco::PFRecHitCollection& rechits ){

  topoClusters_.clear(); 
  
#ifdef PFLOW_DEBUG
  if(debug_) 
    cout<<"PFClusterAlgo::buildTopoClusters start"<<endl;
#endif
  
  for(unsigned is = 0; is<seeds_.size(); is++) {
    
    unsigned rhi = seeds_[is];

    if( !masked(rhi) ) continue;
    // rechit was masked to be processed

    // already used in a topological cluster
    if( usedInTopo_[rhi] ) {
#ifdef PFLOW_DEBUG
      if(debug_) 
        cout<<rhi<<" used"<<endl; 
#endif
      continue;
    }
    
    vector< unsigned > topocluster;
    buildTopoCluster( topocluster, rhi, rechits );
   
    if(topocluster.empty() ) continue;
    
    topoClusters_.push_back( topocluster );

  }

#ifdef PFLOW_DEBUG
  if(debug_) 
    cout<<"PFClusterAlgo::buildTopoClusters done"<<endl;
#endif
  
  return;
}


void 
PFClusterAlgo::buildTopoCluster( vector< unsigned >& cluster,
                                 unsigned rhi, 
                                 const reco::PFRecHitCollection& rechits ){


#ifdef PFLOW_DEBUG
  if(debug_)
    cout<<"PFClusterAlgo::buildTopoCluster in"<<endl;
#endif

  const reco::PFRecHit& rh = rechit( rhi, rechits); 

  double e = rh.energy();
  int layer = rh.layer();
  

  int iring = 0;
  if (layer==PFLayer::HCAL_BARREL2 && abs(rh.positionREP().Eta())>0.34) iring= 1;

  double thresh = parameter( THRESH, 
                             static_cast<PFLayer::Layer>(layer), 0, iring );
  double ptThresh = parameter( PT_THRESH, 
                               static_cast<PFLayer::Layer>(layer), 0, iring );


  if( e < thresh ||  (ptThresh > 0. && rh.pt2() < ptThresh*ptThresh) ) {
#ifdef PFLOW_DEBUG
    if(debug_)
      cout<<"return : "<<e<<"<"<<thresh<<endl; 
#endif
    return;
  }

  // add hit to cluster

  cluster.push_back( rhi );
  // idUsedRecHits_.insert( rh.detId() );

  usedInTopo_[ rhi ] = true;

  //   cout<<" hit ptr "<<hit<<endl;

  // get neighbours, either with one side in common, 
  // or with one corner in common (if useCornerCells_)
  const std::vector< unsigned >& nbs 
    = useCornerCells_ ? rh.neighbours8() : rh.neighbours4();
  
  for(unsigned i=0; i<nbs.size(); i++) {

//     const reco::PFRecHit& neighbour = rechit( nbs[i], rechits );

//     set<unsigned>::iterator used 
//       = idUsedRecHits_.find( neighbour.detId() );
//     if(used != idUsedRecHits_.end() ) continue;
    
    // already used
    if( usedInTopo_[ nbs[i] ] ) {
#ifdef PFLOW_DEBUG
      if(debug_) 
        cout<<rhi<<" used"<<endl; 
#endif
      continue;
    }
                             
    if( !masked(nbs[i]) ) continue;
    buildTopoCluster( cluster, nbs[i], rechits );
  }
#ifdef PFLOW_DEBUG
  if(debug_)
    cout<<"PFClusterAlgo::buildTopoCluster out"<<endl;
#endif

}


void 
PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,
                                const reco::PFRecHitCollection& rechits ) 
{


  //  bool debug = false;


  // several rechits may be seeds. initialize PFClusters on these seeds. 
  
  vector<reco::PFCluster> curpfclusters;
  vector<reco::PFCluster> curpfclusterswodepthcor;
  vector< unsigned > seedsintopocluster;


  for(unsigned i=0; i<topocluster.size(); i++ ) {

    unsigned rhi = topocluster[i];

    if( seedStates_[rhi] == YES ) {

      reco::PFCluster cluster;
      reco::PFCluster clusterwodepthcor;

      double fraction = 1.0; 
      
      reco::PFRecHitRef  recHitRef = createRecHitRef( rechits, rhi ); 

      cluster.addRecHitFraction( reco::PFRecHitFraction( recHitRef, 
                                                         fraction ) );

    // cluster.addRecHit( rhi, fraction );

      calculateClusterPosition( cluster,
                                clusterwodepthcor, 
                                true );    


//       cout<<"PFClusterAlgo: 2"<<endl;
      curpfclusters.push_back( cluster );
      curpfclusterswodepthcor.push_back( clusterwodepthcor );
#ifdef PFLOW_DEBUG
      if(debug_) {
        cout << "PFClusterAlgo::buildPFClusters: seed "
             << rechit( rhi, rechits) <<endl;
        cout << "PFClusterAlgo::buildPFClusters: pfcluster initialized : "
             << cluster <<endl;
      }
#endif

      // keep track of the seed of each topocluster
      seedsintopocluster.push_back( rhi );
    }
  }

  // if only one seed in the topocluster, use all crystals
  // in the position calculation (posCalcNCrystal = -1)
  // otherwise, use the user specified value
  int posCalcNCrystal = seedsintopocluster.size()>1 ? posCalcNCrystal_:-1;
  double ns2 = std::max(1.,(double)(seedsintopocluster.size())-1.);
  ns2 *= ns2;
    
  // Find iteratively the energy and position
  // of each pfcluster in the topological cluster
  unsigned iter = 0;
  unsigned niter = 50;
  double diff = ns2;

  // if(debug_) niter=2;
  vector<double> ener;
  vector<double> dist;
  vector<double> frac;
  vector<math::XYZVector> tmp;

  while ( iter++ < niter && diff > 1E-8*ns2 ) {

    // Store previous iteration's result and reset pfclusters     
    ener.clear();
    tmp.clear();

    for ( unsigned ic=0; ic<curpfclusters.size(); ic++ ) {
      ener.push_back( curpfclusters[ic].energy() );
      
      math::XYZVector v;
      v = curpfclusters[ic].position();

      tmp.push_back( v );

#ifdef PFLOW_DEBUG
      if(debug_)  {
        cout<<"saving photon pos "<<ic<<" "<<curpfclusters[ic]<<endl;
        cout<<tmp[ic].X()<<" "<<tmp[ic].Y()<<" "<<tmp[ic].Z()<<endl;
      }
#endif

      curpfclusters[ic].reset();
    }

    // Loop over topocluster cells
    for( unsigned irh=0; irh<topocluster.size(); irh++ ) {
      unsigned rhindex = topocluster[irh];
      
      const reco::PFRecHit& rh = rechit( rhindex, rechits);
      
      // int layer = rh.layer();
             
      dist.clear();
      frac.clear();
      double fractot = 0.;

      bool isaseed = isSeed(rhindex);

      math::XYZVector cposxyzcell;
      cposxyzcell = rh.position();

#ifdef PFLOW_DEBUG
      if(debug_) { 
        cout<<rh<<endl;
        cout<<"start loop on curpfclusters"<<endl;
      }
#endif

      // Loop over pfclusters
      for ( unsigned ic=0; ic<tmp.size(); ic++) {
        
#ifdef PFLOW_DEBUG
        if(debug_) cout<<"pfcluster "<<ic<<endl;
#endif
        
        double frc=0.;
        bool seedexclusion=true;

        // convert cluster coordinates to xyz
        //math::XYZVector cposxyzclust( tmp[ic].X(), tmp[ic].Y(), tmp[ic].Z() );
        // cluster position used to compute distance with cell
        math::XYZVector cposxyzclust;
        cposxyzclust = curpfclusterswodepthcor[ic].position();

#ifdef PFLOW_DEBUG
        if(debug_) {
          
          cout<<"CLUSTER "<<cposxyzclust.X()<<","
              <<cposxyzclust.Y()<<","
              <<cposxyzclust.Z()<<"\t\t"
              <<"CELL "<<cposxyzcell.X()<<","
              <<cposxyzcell.Y()<<","
              <<cposxyzcell.Z()<<endl;
        }  
#endif
        
        // Compute the distance between the current cell 
        // and the current PF cluster, normalized to a 
        // number of "sigma"
        math::XYZVector deltav = cposxyzclust;
        deltav -= cposxyzcell;
        double d = deltav.R() / showerSigma_;
        
        // if distance cell-cluster is too large, it means that 
        // we're right on the junction between 2 subdetectors (HCAL/VFCAL...)
        // in this case, distance is calculated in the xy plane
        // could also be a large supercluster... 
#ifdef PFLOW_DEBUG
        if( d > 10. && debug_ ) { 
          paint(rhindex, SPECIAL);
          cout<<"PFClusterAlgo Warning: distance too large"<<d<<endl;
        }
#endif
        dist.push_back( d );

        // the current cell is the seed from the current photon.
        if( rhindex == seedsintopocluster[ic] && seedexclusion ) {
          frc = 1.;
#ifdef PFLOW_DEBUG
          if(debug_) cout<<"this cell is a seed for the current photon"<<endl;
#endif
        }
        else if( isaseed && seedexclusion ) {
          frc = 0.;
#ifdef PFLOW_DEBUG
          if(debug_) cout<<"this cell is a seed for another photon"<<endl;
#endif
        }
        else {
          // Compute the fractions of the cell energy to be assigned to 
          // each curpfclusters in the cluster.
          frc = ener[ic] * exp ( - dist[ic]*dist[ic] / 2. );

#ifdef PFLOW_DEBUG
          if(debug_) {
            cout<<"dist["<<ic<<"] "<<dist[ic]
              //                <<", sigma="<<sigma
                <<", frc="<<frc<<endl;
          }  
#endif
        
        }
        fractot += frc;
        frac.push_back(frc);
      }      

      // Add the relevant fraction of the cell to the curpfclusters
#ifdef PFLOW_DEBUG
      if(debug_) cout<<"start add cell"<<endl;
#endif
      for ( unsigned ic=0; ic<tmp.size(); ++ic ) {
#ifdef PFLOW_DEBUG
        if(debug_) 
          cout<<" frac["<<ic<<"] "<<frac[ic]<<" "<<fractot<<" "<<rh<<endl;
#endif

        if( fractot ) 
          frac[ic] /= fractot;
        else { 
#ifdef PFLOW_DEBUG
          if( debug_ ) {
            int layer = rh.layer();
            cerr<<"fractot = 0 ! "<<layer<<endl;
            
            for( unsigned trh=0; trh<topocluster.size(); trh++ ) {
              unsigned tindex = topocluster[trh];
              const reco::PFRecHit& rh = rechit( tindex, rechits);
              cout<<rh<<endl;
            }

            // assert(0)
          }
#endif

          continue;
        }

        // if the fraction has been set to 0, the cell 
        // is now added to the cluster - careful ! (PJ, 19/07/08)
        // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
        // (PJ, 15/09/08 <- similar to what existed before the 
        // previous bug fix, but keeps the close seeds inside, 
        // even if their fraction was set to zero.)
        // Also add a protection to keep the seed in the cluster 
        // when the latter gets far from the former. These cases
        // (about 1% of the clusters) need to be studied, as 
        // they create fake photons, in general.
        // (PJ, 16/09/08) 
        if ( dist[ic] < 10. || frac[ic] > 0.99999 ) { 
          // if ( dist[ic] > 6. ) cout << "Warning : PCluster is getting very far from its seeding cell" << endl;
          reco::PFRecHitRef  recHitRef = createRecHitRef( rechits, rhindex ); 
          reco::PFRecHitFraction rhf( recHitRef,frac[ic] );
          curpfclusters[ic].addRecHitFraction( rhf );
        }
      }
      // if(debug_) cout<<" end add cell"<<endl;
    }

    // Determine the new cluster position and check 
    // the distance with the previous iteration
    diff = 0.;
    for (  unsigned ic=0; ic<tmp.size(); ++ic ) {

      calculateClusterPosition( curpfclusters[ic], curpfclusterswodepthcor[ic], 
                                true, posCalcNCrystal );
#ifdef PFLOW_DEBUG
      if(debug_) cout<<"new iter "<<ic<<endl;
      if(debug_) cout<<curpfclusters[ic]<<endl;
#endif

      double delta = ROOT::Math::VectorUtil::DeltaR(curpfclusters[ic].position(),tmp[ic]);
      if ( delta > diff ) diff = delta;
    }
  }
  
  // Issue a warning message if the number of iterations 
  // exceeds 50
#ifdef PFLOW_DEBUG
  if ( iter >= 50 && debug_ ) 
    cout << "PFClusterAlgo Warning: "
         << "more than "<<niter<<" iterations in pfcluster finding: " 
         <<  setprecision(10) << diff << endl;
#endif
  
  // There we go
  // add all clusters to the list of pfClusters.
  for(unsigned ic=0; ic<curpfclusters.size(); ic++) {

    calculateClusterPosition(curpfclusters[ic], curpfclusterswodepthcor[ic], 
                             true, posCalcNCrystal);

    pfClusters_->push_back(curpfclusters[ic]); 
  }
}



void 
PFClusterAlgo::calculateClusterPosition(reco::PFCluster& cluster,
                                        reco::PFCluster& clusterwodepthcor,
                                        bool depcor, 
                                        int posCalcNCrystal) {

  if( posCalcNCrystal_ != -1 && 
      posCalcNCrystal_ != 5 && 
      posCalcNCrystal_ != 9 ) {
    throw "PFCluster::calculatePosition : posCalcNCrystal_ must be -1, 5, or 9.";
  }  


  if(!posCalcNCrystal) posCalcNCrystal = posCalcNCrystal_; 

  cluster.position_.SetXYZ(0,0,0);

  cluster.energy_ = 0;
  
  double normalize = 0;

  // calculate total energy, average layer, and look for seed  ---------- //


  // double layer = 0;
  map <PFLayer::Layer, double> layers; 

  unsigned seedIndex = 0;
  bool     seedIndexFound = false;

  //Colin: the following can be simplified!

  // loop on rechit fractions
  for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {

    unsigned rhi = cluster.rechits_[ic].recHitRef().index();
    // const reco::PFRecHit& rh = rechit( rhi, rechits );

    const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
    double fraction =  cluster.rechits_[ic].fraction();

    // Find the seed of this sub-cluster (excluding other seeds found in the topological
    // cluster, the energy fraction of which were set to 0 fpr the position determination.
    if( isSeed(rhi) && fraction > 1e-9 ) {
      seedIndex = rhi;
      seedIndexFound = true;
    }

    double recHitEnergy = rh.energy() * fraction;

    // is nan ? 
    if( recHitEnergy!=recHitEnergy ) {
      ostringstream ostr;
      edm::LogError("PFClusterAlgo")<<"rechit "<<rh.detId()<<" has a NaN energy... The input of the particle flow clustering seems to be corrupted.";
    }

    cluster.energy_ += recHitEnergy;

    // sum energy in each layer
    PFLayer::Layer layer = rh.layer();  

    map <PFLayer::Layer, double>:: iterator it = layers.find(layer);

    if (it != layers.end()) {
      it->second += recHitEnergy;
    } else {

      layers.insert(make_pair(layer, recHitEnergy));

    }

  }  

  assert(seedIndexFound);

  // loop over pairs to find layer with max energy          
  double Emax = 0.;
  PFLayer::Layer layer = PFLayer::NONE;
  for (map<PFLayer::Layer, double>::iterator it = layers.begin();
       it != layers.end(); ++it) {
    double e = it->second;
    if(e > Emax){ 
      Emax = e; 
      layer = it->first;
    }
  }
  

  //setlayer here
  cluster.setLayer( layer ); // take layer with max energy

  // layer /= cluster.energy_;
  // cluster.layer_ = lrintf(layer); // nearest integer



  double p1 =  posCalcP1_;
  if( p1 < 0 ) { 
    // automatic (and hopefully best !) determination of the parameter
    // for position determination.
    
    // Remove the ad-hoc determination of p1, and set it to the 
    // seed threshold.
    switch(cluster.layer() ) {
    case PFLayer::ECAL_BARREL:
    case PFLayer::HCAL_BARREL1:
    case PFLayer::HCAL_BARREL2:
      //    case PFLayer::HCAL_HO:
      p1 = threshBarrel_;
      break;
    case PFLayer::ECAL_ENDCAP:
    case PFLayer::HCAL_ENDCAP:
    case PFLayer::HF_EM:
    case PFLayer::HF_HAD:
    case PFLayer::PS1:
    case PFLayer::PS2:
      p1 = threshEndcap_;
      break;

    /*
    switch(cluster.layer() ) {
    case PFLayer::ECAL_BARREL:
      p1 = 0.004 + 0.022*cluster.energy_; // 27 feb 2006 
      break;
    case PFLayer::ECAL_ENDCAP:
      p1 = 0.014 + 0.009*cluster.energy_; // 27 feb 2006 
      break;
    case PFLayer::HCAL_BARREL1:
    case PFLayer::HCAL_BARREL2:
    case PFLayer::HCAL_ENDCAP:
    case PFLayer::HCAL_HF:
      p1 = 5.41215e-01 * log( cluster.energy_ / 1.29803e+01 );
      if(p1<0.01) p1 = 0.01;
      break;
    */

    default:
      cerr<<"Clusters weight_p1 -1 not yet allowed for layer "<<layer
          <<". Chose a better value in the opt file"<<endl;
      assert(0);
      break;
    }
  } 
  else if( p1< 1e-9 ) { // will divide by p1 later on
    p1 = 1e-9;
  }
  // calculate uncorrected cluster position --------------------------------

  reco::PFCluster::REPPoint clusterpos;   // uncorrected cluster pos 
  math::XYZPoint clusterposxyz;           // idem, xyz coord 
  math::XYZPoint firstrechitposxyz;       // pos of the rechit with highest E

  double maxe = -9999;
  double x = 0;
  double y = 0;
  double z = 0;
  
  for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
    
    unsigned rhi = cluster.rechits_[ic].recHitRef().index();
//     const reco::PFRecHit& rh = rechit( rhi, rechits );

    const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());

    if(rhi != seedIndex) { // not the seed
      if( posCalcNCrystal == 5 ) { // pos calculated from the 5 neighbours only
        if(!rh.isNeighbour4(seedIndex) ) {
          continue;
        }
      }
      if( posCalcNCrystal == 9 ) { // pos calculated from the 9 neighbours only
        if(!rh.isNeighbour8(seedIndex) ) {
          continue;
        }
      }
    }
    double fraction =  cluster.rechits_[ic].fraction();
    double recHitEnergy = rh.energy() * fraction;

    double norm = fraction < 1E-9 ? 0. : max(0., log(recHitEnergy/p1 ));
    
    const math::XYZPoint& rechitposxyz = rh.position();
    
    if( recHitEnergy > maxe ) {
      firstrechitposxyz = rechitposxyz;
      maxe = recHitEnergy;
    }

    x += rechitposxyz.X() * norm;
    y += rechitposxyz.Y() * norm;
    z += rechitposxyz.Z() * norm;
    
    // clusterposxyz += rechitposxyz * norm;
    normalize += norm;
  }

  // normalize uncorrected position
  // assert(normalize);
  if( normalize < 1e-9 ) {
    //    cerr<<"--------------------"<<endl;
    //    cerr<<(*this)<<endl;
    cout << "Watch out : cluster too far from its seeding cell, set to 0,0,0" << endl;
    clusterposxyz.SetXYZ(0,0,0);
    clusterpos.SetCoordinates(0,0,0);
    return;
  }
  else {
    x /= normalize;
    y /= normalize; 
    z /= normalize; 

    clusterposxyz.SetCoordinates( x, y, z);

    clusterpos.SetCoordinates( clusterposxyz.Rho(), clusterposxyz.Eta(), clusterposxyz.Phi() );

  }  

  cluster.posrep_ = clusterpos;

  cluster.position_ = clusterposxyz;

  clusterwodepthcor = cluster;


  // correctio of the rechit position, 
  // according to the depth, only for ECAL 


  if( depcor &&   // correction requested and ECAL
      ( cluster.layer() == PFLayer::ECAL_BARREL ||       
        cluster.layer() == PFLayer::ECAL_ENDCAP ) ) {


    double corra = reco::PFCluster::depthCorA_;
    double corrb = reco::PFCluster::depthCorB_;
    if( abs(clusterpos.Eta() ) < 2.6 && 
        abs(clusterpos.Eta() ) > 1.65   ) { 
      // if crystals under preshower, correction is not the same  
      // (shower depth smaller)
      corra = reco::PFCluster::depthCorAp_;
      corrb = reco::PFCluster::depthCorBp_;
    }

    double depth = 0;

    switch( reco::PFCluster::depthCorMode_ ) {
    case 1: // for e/gamma 
      depth = corra * ( corrb + log(cluster.energy_) ); 
      break;
    case 2: // for hadrons
      depth = corra;
      break;
    default:
      cerr<<"PFClusterAlgo::calculateClusterPosition : unknown function for depth correction! "<<endl;
      assert(0);
    }


    // calculate depth vector:
    // its mag is depth
    // its direction is the cluster direction (uncorrected)

//     double xdepthv = clusterposxyz.X();
//     double ydepthv = clusterposxyz.Y();
//     double zdepthv = clusterposxyz.Z();
//     double mag = sqrt( xdepthv*xdepthv + 
//                     ydepthv*ydepthv + 
//                     zdepthv*zdepthv );
    

//     math::XYZPoint depthv(clusterposxyz); 
//     depthv.SetMag(depth);
    
    
    math::XYZVector depthv( clusterposxyz.X(), 
                            clusterposxyz.Y(),
                            clusterposxyz.Z() );
    depthv /= sqrt(depthv.Mag2() );
    depthv *= depth;


    // now calculate corrected cluster position:    
    math::XYZPoint clusterposxyzcor;

    maxe = -9999;
    x = 0;
    y = 0;
    z = 0;
    cluster.posrep_.SetXYZ(0,0,0);
    normalize = 0;
    for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {

      unsigned rhi = cluster.rechits_[ic].recHitRef().index();
//       const reco::PFRecHit& rh = rechit( rhi, rechits );
      
      const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());

      if(rhi != seedIndex) {
        if( posCalcNCrystal == 5 ) {
          if(!rh.isNeighbour4(seedIndex) ) {
            continue;
          }
        }
        if( posCalcNCrystal == 9 ) {
          if(!rh.isNeighbour8(seedIndex) ) {
            continue;
          }
        }
      }

      double fraction =  cluster.rechits_[ic].fraction();
      double recHitEnergy = rh.energy() * fraction;
      
      const math::XYZPoint&  rechitposxyz = rh.position();

      // rechit axis not correct ! 
      math::XYZVector rechitaxis = rh.getAxisXYZ();
      // rechitaxis -= math::XYZVector( rechitposxyz.X(), rechitposxyz.Y(), rechitposxyz.Z() );
      
      math::XYZVector rechitaxisu( rechitaxis );
      rechitaxisu /= sqrt( rechitaxis.Mag2() );

      math::XYZVector displacement( rechitaxisu );
      // displacement /= sqrt( displacement.Mag2() );    
      displacement *= rechitaxisu.Dot( depthv );
      
      math::XYZPoint rechitposxyzcor( rechitposxyz );
      rechitposxyzcor += displacement;

      if( recHitEnergy > maxe ) {
        firstrechitposxyz = rechitposxyzcor;
        maxe = recHitEnergy;
      }

      double norm = fraction < 1E-9 ? 0. : max(0., log(recHitEnergy/p1 ));
      
      x += rechitposxyzcor.X() * norm;
      y += rechitposxyzcor.Y() * norm;
      z += rechitposxyzcor.Z() * norm;
      
      // clusterposxyzcor += rechitposxyzcor * norm;
      normalize += norm;
    }

    // normalize
    if(normalize < 1e-9) {
      cerr<<"--------------------"<<endl;
      cerr<< cluster <<endl;
      assert(0);
    }
    else {
      x /= normalize;
      y /= normalize;
      z /= normalize;
      

      clusterposxyzcor.SetCoordinates(x,y,z);
      cluster.posrep_.SetCoordinates( clusterposxyzcor.Rho(), 
                                      clusterposxyzcor.Eta(), 
                                      clusterposxyzcor.Phi() );
      cluster.position_  = clusterposxyzcor;
      clusterposxyz = clusterposxyzcor;
    }

  }
}



const reco::PFRecHit& 
PFClusterAlgo::rechit(unsigned i, 
                      const reco::PFRecHitCollection& rechits ) {

  if(i >= rechits.size() ) { // i >= 0, since i is unsigned
    string err = "PFClusterAlgo::rechit : out of range";
    throw std::out_of_range(err);
  }
  
  return rechits[i];
}



bool PFClusterAlgo::masked(unsigned rhi) const {

  if(rhi>=mask_.size() ) { // rhi >= 0, since rhi is unsigned
    string err = "PFClusterAlgo::masked : out of range";
    throw std::out_of_range(err);
  }
  
  return mask_[rhi];
}


unsigned PFClusterAlgo::color(unsigned rhi) const {

  if(rhi>=color_.size() ) { // rhi >= 0, since rhi is unsigned
    string err = "PFClusterAlgo::color : out of range";
    throw std::out_of_range(err);
  }
  
  return color_[rhi];
}



bool PFClusterAlgo::isSeed(unsigned rhi) const {

  if(rhi>=seedStates_.size() ) { // rhi >= 0, since rhi is unsigned
    string err = "PFClusterAlgo::isSeed : out of range";
    throw std::out_of_range(err);
  }
  
  return seedStates_[rhi]>0 ? true : false;
}


void PFClusterAlgo::paint(unsigned rhi, unsigned color ) {

  if(rhi>=color_.size() ) { // rhi >= 0, since rhi is unsigned
    string err = "PFClusterAlgo::color : out of range";
    throw std::out_of_range(err);
  }
  
  color_[rhi] = color;
}


reco::PFRecHitRef 
PFClusterAlgo::createRecHitRef( const reco::PFRecHitCollection& rechits, 
                                unsigned rhi ) {

  if( rechitsHandle_.isValid() ) {
    return reco::PFRecHitRef(  rechitsHandle_, rhi );
  } 
  else {
    return reco::PFRecHitRef(  &rechits, rhi );
  }
}


ostream& operator<<(ostream& out,const PFClusterAlgo& algo) {
  if(!out) return out;
  out<<"PFClusterAlgo parameters : "<<endl;
  out<<"-----------------------------------------------------"<<endl;
  out<<"threshBarrel       : "<<algo.threshBarrel_       <<endl;
  out<<"threshSeedBarrel   : "<<algo.threshSeedBarrel_   <<endl;
  out<<"threshPtBarrel     : "<<algo.threshPtBarrel_     <<endl;
  out<<"threshPtSeedBarrel : "<<algo.threshPtSeedBarrel_ <<endl;
  out<<"threshCleanBarrel  : "<<algo.threshCleanBarrel_  <<endl;
  out<<"minS4S1Barrel      : "<<algo.minS4S1Barrel_[0]<<" x log10(E) + "<<algo.minS4S1Barrel_[1]<<endl;
  out<<"threshEndcap       : "<<algo.threshEndcap_       <<endl;
  out<<"threshSeedEndcap   : "<<algo.threshSeedEndcap_   <<endl;
  out<<"threshPtEndcap     : "<<algo.threshPtEndcap_     <<endl;
  out<<"threshPtSeedEndcap : "<<algo.threshPtSeedEndcap_ <<endl;
  out<<"threshEndcap       : "<<algo.threshEndcap_       <<endl;
  out<<"threshCleanEndcap  : "<<algo.threshCleanEndcap_  <<endl;
  out<<"minS4S1Endcap      : "<<algo.minS4S1Endcap_[0]<<" x log10(E) + "<<algo.minS4S1Endcap_[1]<<endl;
  out<<"nNeighbours        : "<<algo.nNeighbours_        <<endl;
  out<<"posCalcNCrystal    : "<<algo.posCalcNCrystal_    <<endl;
  out<<"posCalcP1          : "<<algo.posCalcP1_          <<endl;
  out<<"showerSigma        : "<<algo.showerSigma_        <<endl;
  out<<"useCornerCells     : "<<algo.useCornerCells_     <<endl;

  out<<endl;
  out<<algo.pfClusters_->size()<<" clusters:"<<endl;

  for(unsigned i=0; i<algo.pfClusters_->size(); i++) {
    out<<(*algo.pfClusters_)[i]<<endl;
    
    if(!out) return out;
  }
  
  return out;
}


//compute the unsigned distance to the closest phi-crack in the barrel
std::pair<double,double>
PFClusterAlgo::dCrack(double phi, double eta){

  static double pi= M_PI;// 3.14159265358979323846;
  
  //Location of the 18 phi-cracks
  static std::vector<double> cPhi;
  if(cPhi.size()==0)
    {
      cPhi.resize(18,0);
      cPhi[0]=2.97025;
      for(unsigned i=1;i<=17;++i) cPhi[i]=cPhi[0]-2*i*pi/18;
    }

  //Shift of this location if eta<0
  static double delta_cPhi=0.00638;

  double defi; //the result

  //the location is shifted
  if(eta<0) phi +=delta_cPhi;

  if (phi>=-pi && phi<=pi){

    //the problem of the extrema
    if (phi<cPhi[17] || phi>=cPhi[0]){
      if (phi<0) phi+= 2*pi;
      defi = std::min(fabs(phi -cPhi[0]),fabs(phi-cPhi[17]-2*pi));              
    }

    //between these extrema...
    else{
      bool OK = false;
      unsigned i=16;
      while(!OK){
        if (phi<cPhi[i]){
          defi=std::min(fabs(phi-cPhi[i+1]),fabs(phi-cPhi[i]));
          OK=true;
        }
        else i-=1;
      }
    }
  }
  else{
    defi=0.;        //if there is a problem, we assum that we are in a crack
    std::cout<<"Problem in dminphi"<<std::endl;
  }
  //if(eta<0) defi=-defi;   //because of the disymetry

  static std::vector<double> cEta;
  if ( cEta.size() == 0 ) { 
    cEta.push_back(0.0);
    cEta.push_back(4.44747e-01)  ;   cEta.push_back(-4.44747e-01)  ;
    cEta.push_back(7.92824e-01)  ;   cEta.push_back(-7.92824e-01) ;
    cEta.push_back(1.14090e+00)  ;   cEta.push_back(-1.14090e+00)  ;
    cEta.push_back(1.47464e+00)  ;   cEta.push_back(-1.47464e+00)  ;
  }
  double deta = 999.; // the other result

  for ( unsigned ieta=0; ieta<cEta.size(); ++ieta ) { 
    deta = std::min(deta,fabs(eta-cEta[ieta]));
  }
  
  defi /= 0.0175;
  deta /= 0.0175;
  return std::pair<double,double>(defi,deta);

}