CMS 3D CMS Logo

Public Types | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | Friends

PFClusterAlgo Class Reference

Algorithm for particle flow clustering. More...

#include <PFClusterAlgo.h>

List of all members.

Public Types

enum  Color { NONE = 0, SEED, SPECIAL }
typedef std::multimap< double,
unsigned >::iterator 
EH
typedef std::map< unsigned,
unsigned >::const_iterator 
IDH
enum  Parameter {
  THRESH, SEED_THRESH, PT_THRESH, SEED_PT_THRESH,
  CLEAN_THRESH, CLEAN_S4S1, DOUBLESPIKE_THRESH, DOUBLESPIKE_S6S2
}
typedef edm::Handle
< reco::PFRecHitCollection
PFRecHitHandle
enum  SeedState { UNKNOWN = -1, NO = 0, YES = 1, CLEAN = 2 }

Public Member Functions

std::auto_ptr< std::vector
< reco::PFCluster > > & 
clusters ()
unsigned color (unsigned rhi) const
void doClustering (const reco::PFRecHitCollection &rechits)
 perform clustering
void doClustering (const reco::PFRecHitCollection &rechits, const std::vector< bool > &mask)
void doClustering (const PFRecHitHandle &rechitsHandle)
 perform clustering in full framework
void doClustering (const PFRecHitHandle &rechitsHandle, const std::vector< bool > &mask)
void enableDebugging (bool debug)
 set hits on which clustering will be performed
bool isSeed (unsigned rhi) const
bool masked (unsigned rhi) const
int nNeighbours () const
 get number of neighbours for
double parameter (Parameter paramtype, PFLayer::Layer layer, unsigned iCoeff=0, int iring0=0) const
 PFClusterAlgo ()
 constructor
int posCalcNCrystal () const
 get number of crystals for position calculation (-1 all,5, or 9)
double posCalcP1 () const
 get p1 for position calculation
const reco::PFRecHitrechit (unsigned i, const reco::PFRecHitCollection &rechits)
 ----------------------------------------------------------------
std::auto_ptr< std::vector
< reco::PFRecHit > > & 
rechitsCleaned ()
void setCleanRBXandHPDs (bool cleanRBXandHPDs)
 Activate cleaning of HCAL RBX's and HPD's.
void setHistos (TFile *file, TH2F *hB, TH2F *hE)
 set endcap clean threshold
void setNNeighbours (int n)
 set number of neighbours for
void setPosCalcNCrystal (int n)
 set number of crystals for position calculation (-1 all,5, or 9)
void setPosCalcP1 (double p1)
 set p1 for position calculation
void setS4S1CleanBarrel (const std::vector< double > &coeffs)
void setS4S1CleanEndcap (const std::vector< double > &coeffs)
void setS6S2DoubleSpikeBarrel (double cut)
void setS6S2DoubleSpikeEndcap (double cut)
void setShowerSigma (double sigma)
 set shower sigma for
void setThreshBarrel (double thresh)
 setters -------------------------------------------------------
void setThreshCleanBarrel (double thresh)
 set barrel clean threshold
void setThreshCleanEndcap (double thresh)
 set endcap clean threshold
void setThreshDoubleSpikeBarrel (double thresh)
 set endcap thresholds for double spike cleaning
void setThreshDoubleSpikeEndcap (double thresh)
 set endcap thresholds for double spike cleaning
void setThreshEndcap (double thresh)
 set endcap threshold
void setThreshPtBarrel (double thresh)
void setThreshPtEndcap (double thresh)
void setThreshPtSeedBarrel (double thresh)
void setThreshPtSeedEndcap (double thresh)
void setThreshSeedBarrel (double thresh)
 set barrel seed threshold
void setThreshSeedEndcap (double thresh)
 set endcap seed threshold
void setUseCornerCells (bool usecornercells)
 activate use of cells with a common corner to build topo-clusters
double showerSigma () const
 get shower sigma
double threshBarrel () const
 getters -------------------------------------------------------
double threshEndcap () const
 get endcap threshold
double threshSeedBarrel () const
 get barrel seed threshold
double threshSeedEndcap () const
 get endcap seed threshold
void write ()
 write histos
virtual ~PFClusterAlgo ()
 destructor

Private Member Functions

void buildPFClusters (const std::vector< unsigned > &cluster, const reco::PFRecHitCollection &rechits)
 build PFClusters from a topocluster
void buildTopoCluster (std::vector< unsigned > &cluster, unsigned rhi, const reco::PFRecHitCollection &rechits)
 build a topocluster (recursive)
void buildTopoClusters (const reco::PFRecHitCollection &rechits)
 build topoclusters around seeds
void calculateClusterPosition (reco::PFCluster &cluster, reco::PFCluster &clusterwodepthcor, bool depcor=true, int posCalcNCrystal=0)
 calculate position of a cluster
void cleanRBXAndHPD (const reco::PFRecHitCollection &rechits)
 Clean HCAL readout box noise and HPD discharge.
reco::PFRecHitRef createRecHitRef (const reco::PFRecHitCollection &rechits, unsigned rhi)
std::pair< double, double > dCrack (double phi, double eta)
 distance to a crack in the ECAL barrel in eta and phi direction
void doClusteringWorker (const reco::PFRecHitCollection &rechits)
 perform clustering
void findSeeds (const reco::PFRecHitCollection &rechits)
 look for seeds
void paint (unsigned rhi, unsigned color=1)
 paint a rechit with a color.

Private Attributes

bool cleanRBXandHPDs_
 option to clean HCAL RBX's and HPD's
std::vector< unsigned > color_
 color, for all rechits
bool debug_
 debugging on/off
std::multimap< double,
unsigned, std::greater< double > > 
eRecHits_
 indices to rechits, sorted by decreasing E (not E_T)
TFile * file_
TH2F * hBNeighbour
TH2F * hENeighbour
std::set< unsigned > idUsedRecHits_
 ids of rechits used in seed search
std::vector< bool > mask_
std::vector< double > minS4S1Barrel_
std::vector< double > minS4S1Endcap_
double minS6S2DoubleSpikeBarrel_
double minS6S2DoubleSpikeEndcap_
int nNeighbours_
 number of neighbours
std::auto_ptr< std::vector
< reco::PFCluster > > 
pfClusters_
 all clusters
std::auto_ptr< std::vector
< reco::PFRecHit > > 
pfRecHitsCleaned_
 particle flow rechits cleaned
int posCalcNCrystal_
 number of crystals for position calculation
double posCalcP1_
 parameter for position calculation
PFRecHitHandle rechitsHandle_
std::vector< unsigned > seeds_
 vector of indices for seeds.
std::vector< SeedStateseedStates_
 seed state, for all rechits
double showerSigma_
 sigma of shower (cm)
double threshBarrel_
 barrel threshold
double threshCleanBarrel_
 Barrel cleaning threshold and S4/S1 smallest fractiom.
double threshCleanEndcap_
 Endcap cleaning threshold and S4/S1 smallest fractiom.
double threshDoubleSpikeBarrel_
 Barrel double-spike cleaning.
double threshDoubleSpikeEndcap_
 Endcap double-spike cleaning.
double threshEndcap_
 endcap threshold
double threshPtBarrel_
double threshPtEndcap_
double threshPtSeedBarrel_
double threshPtSeedEndcap_
double threshSeedBarrel_
 barrel seed threshold
double threshSeedEndcap_
 endcap seed threshold
std::vector< std::vector
< unsigned > > 
topoClusters_
 sets of cells having one common side, and energy over threshold
bool useCornerCells_
 option to use cells with a common corner to build topo-clusters
std::vector< bool > usedInTopo_
 used in topo cluster? for all rechits

Static Private Attributes

static unsigned prodNum_ = 1
 product number

Friends

std::ostream & operator<< (std::ostream &out, const PFClusterAlgo &algo)

Detailed Description

Algorithm for particle flow clustering.

This class takes as an input a map of pointers to PFRecHit's, and creates PFCluster's from these rechits. Clustering is implemented for ECAL, HCAL, and preshower.

Author:
Colin Bernet, Patrick Janot
Date:
July 2006

Definition at line 31 of file PFClusterAlgo.h.


Member Typedef Documentation

typedef std::multimap<double, unsigned >::iterator PFClusterAlgo::EH

Definition at line 204 of file PFClusterAlgo.h.

typedef std::map<unsigned, unsigned >::const_iterator PFClusterAlgo::IDH

Definition at line 203 of file PFClusterAlgo.h.

Definition at line 48 of file PFClusterAlgo.h.


Member Enumeration Documentation

Enumerator:
NONE 
SEED 
SPECIAL 

Definition at line 154 of file PFClusterAlgo.h.

             { 
    NONE=0,
    SEED, 
    SPECIAL 
  };
Returns:
threshold, seed threshold, (gaussian width, p1 ??) for a given zone (endcap, barrel, VFCAL ??)
Enumerator:
THRESH 
SEED_THRESH 
PT_THRESH 
SEED_PT_THRESH 
CLEAN_THRESH 
CLEAN_S4S1 
DOUBLESPIKE_THRESH 
DOUBLESPIKE_S6S2 

Definition at line 177 of file PFClusterAlgo.h.

Enumerator:
UNKNOWN 
NO 
YES 
CLEAN 

Definition at line 193 of file PFClusterAlgo.h.

                 {
    UNKNOWN=-1,
    NO=0,
    YES=1,
    CLEAN=2
  };

Constructor & Destructor Documentation

PFClusterAlgo::PFClusterAlgo ( )
virtual PFClusterAlgo::~PFClusterAlgo ( ) [inline, virtual]

destructor

Definition at line 39 of file PFClusterAlgo.h.

{;}

Member Function Documentation

void PFClusterAlgo::buildPFClusters ( const std::vector< unsigned > &  cluster,
const reco::PFRecHitCollection rechits 
) [private]

build PFClusters from a topocluster

Definition at line 1002 of file PFClusterAlgo.cc.

References reco::PFCluster::addRecHitFraction(), calculateClusterPosition(), dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, createRecHitRef(), debug_, delta, diffTreeTool::diff, relval_parameters_module::energy, create_public_lumi_plots::exp, cropTnPTrees::frac, i, isSeed(), reco::PFRecHit::layer(), max(), paint(), pfClusters_, posCalcNCrystal(), posCalcNCrystal_, position, reco::PFRecHit::position(), rechit(), seedStates_, showerSigma_, SPECIAL, tmp, v, and YES.

Referenced by doClusteringWorker().

{


  //  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::buildTopoCluster ( std::vector< unsigned > &  cluster,
unsigned  rhi,
const reco::PFRecHitCollection rechits 
) [private]

build a topocluster (recursive)

Definition at line 926 of file PFClusterAlgo.cc.

References abs, gather_cfg::cout, debug_, alignCSCRings::e, reco::PFRecHit::energy(), PFLayer::HCAL_BARREL2, i, reco::PFRecHit::layer(), masked(), reco::PFRecHit::neighbours4(), reco::PFRecHit::neighbours8(), parameter(), reco::PFRecHit::positionREP(), reco::PFRecHit::pt2(), PT_THRESH, rechit(), GOODCOLL_filter_cfg::thresh, THRESH, useCornerCells_, and usedInTopo_.

Referenced by buildTopoClusters().

                                                                        {


#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::buildTopoClusters ( const reco::PFRecHitCollection rechits) [private]

build topoclusters around seeds

Definition at line 882 of file PFClusterAlgo.cc.

References buildTopoCluster(), gather_cfg::cout, debug_, masked(), seeds_, topoClusters_, and usedInTopo_.

Referenced by doClusteringWorker().

                                                                            {

  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::calculateClusterPosition ( reco::PFCluster cluster,
reco::PFCluster clusterwodepthcor,
bool  depcor = true,
int  posCalcNCrystal = 0 
) [private]

calculate position of a cluster

Definition at line 1294 of file PFClusterAlgo.cc.

References abs, dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, reco::PFCluster::depthCorA_, reco::PFCluster::depthCorAp_, reco::PFCluster::depthCorB_, reco::PFCluster::depthCorBp_, reco::PFCluster::depthCorMode_, reco::PFRecHit::detId(), alignCSCRings::e, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, reco::PFRecHit::energy(), reco::CaloCluster::energy_, reco::PFRecHit::getAxisXYZ(), PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, reco::PFRecHit::isNeighbour4(), reco::PFRecHit::isNeighbour8(), isSeed(), reco::PFCluster::layer(), reco::PFRecHit::layer(), create_public_lumi_plots::log, max(), NONE, p1, posCalcNCrystal_, posCalcP1_, reco::PFRecHit::position(), reco::CaloCluster::position_, reco::PFCluster::posrep_, PFLayer::PS1, PFLayer::PS2, reco::PFCluster::rechits_, reco::PFCluster::setLayer(), mathSSE::sqrt(), threshBarrel_, threshEndcap_, x, detailsBasic3DVector::y, and z.

Referenced by buildPFClusters().

                                                             {

  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;
    }

  }
}
void PFClusterAlgo::cleanRBXAndHPD ( const reco::PFRecHitCollection rechits) [private]

Clean HCAL readout box noise and HPD discharge.

Definition at line 290 of file PFClusterAlgo.cc.

References abs, reco::PFRecHit::detId(), reco::PFRecHit::energy(), eRecHits_, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, HcalDetId::ieta(), recoMuon::in, HcalDetId::iphi(), reco::PFRecHit::layer(), Association::map, mask_, masked(), reco::PFRecHit::neighbours4(), errorMatrix2Lands_multiChannel::nN, pfRecHitsCleaned_, reco::PFRecHit::position(), rechit(), mathSSE::sqrt(), and dtDQMClient_cfg::threshold.

Referenced by doClusteringWorker().

                                                                      {

  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;
          */
        }
      }
  }

}
std::auto_ptr< std::vector< reco::PFCluster > >& PFClusterAlgo::clusters ( ) [inline]
Returns:
particle flow clusters

Definition at line 167 of file PFClusterAlgo.h.

References pfClusters_.

Referenced by PFRootEventManager::clustering().

    {return pfClusters_;}
unsigned PFClusterAlgo::color ( unsigned  rhi) const
Returns:
color of the rechit

Definition at line 1687 of file PFClusterAlgo.cc.

References color_.

Referenced by DisplayManager::loadGRecHits(), and paint().

                                                {

  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];
}
reco::PFRecHitRef PFClusterAlgo::createRecHitRef ( const reco::PFRecHitCollection rechits,
unsigned  rhi 
) [private]

create a reference to a rechit. in case rechitsHandle_.isValid(), this reference is permanent.

Definition at line 1722 of file PFClusterAlgo.cc.

References edm::HandleBase::isValid(), and rechitsHandle_.

Referenced by buildPFClusters().

                                               {

  if( rechitsHandle_.isValid() ) {
    return reco::PFRecHitRef(  rechitsHandle_, rhi );
  } 
  else {
    return reco::PFRecHitRef(  &rechits, rhi );
  }
}
std::pair< double, double > PFClusterAlgo::dCrack ( double  phi,
double  eta 
) [private]

distance to a crack in the ECAL barrel in eta and phi direction

Definition at line 1772 of file PFClusterAlgo.cc.

References gather_cfg::cout, alignCSCRings::e, i, M_PI, min, and pi.

Referenced by findSeeds().

                                           {

  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);

}
void PFClusterAlgo::doClustering ( const PFRecHitHandle rechitsHandle)

perform clustering in full framework

Definition at line 66 of file PFClusterAlgo.cc.

References doClusteringWorker(), mask_, HI_PhotonSkim_cff::rechits, rechitsHandle_, and funct::true.

                                                                      {
  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 
)

Definition at line 80 of file PFClusterAlgo.cc.

References doClusteringWorker(), mask_, HI_PhotonSkim_cff::rechits, rechitsHandle_, and funct::true.

                                                                                                    {

  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)

perform clustering

Definition at line 102 of file PFClusterAlgo.cc.

References edm::HandleBase::clear(), doClusteringWorker(), mask_, rechitsHandle_, and funct::true.

Referenced by PFRootEventManager::clustering().

                                                                        {

  // 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 
)

Definition at line 116 of file PFClusterAlgo.cc.

References edm::HandleBase::clear(), doClusteringWorker(), mask_, rechitsHandle_, and funct::true.

                                                                                                      {
  // 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) [private]

perform clustering

Definition at line 137 of file PFClusterAlgo.cc.

References buildPFClusters(), buildTopoClusters(), cleanRBXAndHPD(), cleanRBXandHPDs_, color_, relval_parameters_module::energy, eRecHits_, funct::false, findSeeds(), i, pfClusters_, pfRecHitsCleaned_, rechit(), seedStates_, topoClusters_, UNKNOWN, and usedInTopo_.

Referenced by doClustering().

                                                                              {


  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 ); 

  }

}
void PFClusterAlgo::enableDebugging ( bool  debug) [inline]

set hits on which clustering will be performed

enable/disable debugging

Definition at line 45 of file PFClusterAlgo.h.

References debug, and debug_.

Referenced by PFRootEventManager::readOptions().

{ debug_ = debug;}
void PFClusterAlgo::findSeeds ( const reco::PFRecHitCollection rechits) [private]

look for seeds

Definition at line 571 of file PFClusterAlgo.cc.

References abs, dtNoiseDBValidation_cfg::cerr, CLEAN, CLEAN_S4S1, CLEAN_THRESH, gather_cfg::cout, dCrack(), debug_, DOUBLESPIKE_S6S2, DOUBLESPIKE_THRESH, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, reco::PFRecHit::energy(), reco::PFRecHit::energyUp(), eRecHits_, eta(), file_, hBNeighbour, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, hENeighbour, PFLayer::HF_EM, PFLayer::HF_HAD, recoMuon::in, reco::PFRecHit::layer(), mask_, masked(), min, reco::PFRecHit::neighbours4(), reco::PFRecHit::neighbours8(), nNeighbours_, NO, paint(), parameter(), pfRecHitsCleaned_, phi, reco::PFRecHit::position(), reco::PFRecHit::positionREP(), PFLayer::PS1, PFLayer::PS2, reco::PFRecHit::pt2(), rechit(), SEED, SEED_PT_THRESH, SEED_THRESH, seeds_, seedStates_, and YES.

Referenced by doClusteringWorker().

                                                                     {

  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
}
bool PFClusterAlgo::isSeed ( unsigned  rhi) const
Returns:
seed flag (not seed state) for rechit with index rhi

Definition at line 1699 of file PFClusterAlgo.cc.

References seedStates_.

Referenced by buildPFClusters(), calculateClusterPosition(), and PFRootEventManager::printRecHits().

                                             {

  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;
}
bool PFClusterAlgo::masked ( unsigned  rhi) const
Returns:
mask flag for the rechit

Definition at line 1676 of file PFClusterAlgo.cc.

References mask_.

Referenced by buildTopoCluster(), buildTopoClusters(), cleanRBXAndHPD(), and findSeeds().

                                             {

  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];
}
int PFClusterAlgo::nNeighbours ( ) const [inline]

get number of neighbours for

Definition at line 130 of file PFClusterAlgo.h.

References nNeighbours_.

{ return nNeighbours_;}
void PFClusterAlgo::paint ( unsigned  rhi,
unsigned  color = 1 
) [private]

paint a rechit with a color.

Definition at line 1710 of file PFClusterAlgo.cc.

References color(), and color_.

Referenced by buildPFClusters(), and findSeeds().

                                                       {

  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;
}
double PFClusterAlgo::parameter ( Parameter  paramtype,
PFLayer::Layer  layer,
unsigned  iCoeff = 0,
int  iring0 = 0 
) const
Returns:
the value of a parameter of a given type, see enum Parameter, in a given layer.

Definition at line 153 of file JetResolution.cc.

References newFWLiteAna::bin, dtNoiseDBValidation_cfg::cerr, i, N, JetResolution::parameterFncs_, JetResolution::parameters_, query::result, and findQualityFiles::size.

Referenced by buildTopoCluster(), and findSeeds().

{
  TF1* result(0);
  for (unsigned i=0;i<parameterFncs_.size()&&result==0;i++) {
    string fncname = parameterFncs_[i]->GetName();
    if (fncname.find("f"+parameterName)==0) {
      stringstream ssname; ssname<<parameterFncs_[i]->GetName();
      for (unsigned ii=0;ii<x.size();ii++)
        ssname<<"_"<<parameters_[i]->definitions().binVar(ii)<<x[ii];
      result = (TF1*)parameterFncs_[i]->Clone();
      result->SetName(ssname.str().c_str());
      int N = parameters_[i]->definitions().nParVar();
      int bin = parameters_[i]->binIndex(x);
      assert(bin>=0);
      assert(bin<(int)parameters_[i]->size());
      const std::vector<float>& pars = parameters_[i]->record(bin).parameters();
      for (unsigned ii=2*N;ii<pars.size();ii++) result->SetParameter(ii-2*N,pars[ii]); 
    }
  }
  
  if (0==result) cerr<<"JetResolution::parameter() ERROR: no parameter "
                     <<parameterName<<" found."<<endl;
  
  return result;
}
int PFClusterAlgo::posCalcNCrystal ( ) const [inline]

get number of crystals for position calculation (-1 all,5, or 9)

Definition at line 136 of file PFClusterAlgo.h.

References posCalcNCrystal_.

Referenced by buildPFClusters().

{return  posCalcNCrystal_;}
double PFClusterAlgo::posCalcP1 ( ) const [inline]

get p1 for position calculation

Definition at line 133 of file PFClusterAlgo.h.

References posCalcP1_.

{ return posCalcP1_; }
const reco::PFRecHit & PFClusterAlgo::rechit ( unsigned  i,
const reco::PFRecHitCollection rechits 
)

----------------------------------------------------------------

Returns:
hit with index i. performs a test on the vector range

Definition at line 1663 of file PFClusterAlgo.cc.

References i.

Referenced by buildPFClusters(), buildTopoCluster(), cleanRBXAndHPD(), doClusteringWorker(), and findSeeds().

                                                              {

  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];
}
std::auto_ptr< std::vector< reco::PFRecHit > >& PFClusterAlgo::rechitsCleaned ( ) [inline]
Returns:
cleaned rechits

Definition at line 171 of file PFClusterAlgo.h.

References pfRecHitsCleaned_.

    {return pfRecHitsCleaned_;}
void PFClusterAlgo::setCleanRBXandHPDs ( bool  cleanRBXandHPDs) [inline]

Activate cleaning of HCAL RBX's and HPD's.

Definition at line 111 of file PFClusterAlgo.h.

References cleanRBXandHPDs_.

Referenced by PFRootEventManager::readOptions().

{ cleanRBXandHPDs_ = cleanRBXandHPDs; }
void PFClusterAlgo::setHistos ( TFile *  file,
TH2F *  hB,
TH2F *  hE 
) [inline]

set endcap clean threshold

Definition at line 93 of file PFClusterAlgo.h.

References mergeVDriftHistosByStation::file, file_, hBNeighbour, and hENeighbour.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setNNeighbours ( int  n) [inline]

set number of neighbours for

Definition at line 96 of file PFClusterAlgo.h.

References n, and nNeighbours_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setPosCalcNCrystal ( int  n) [inline]

set number of crystals for position calculation (-1 all,5, or 9)

Definition at line 102 of file PFClusterAlgo.h.

References n, and posCalcNCrystal_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setPosCalcP1 ( double  p1) [inline]

set p1 for position calculation

Definition at line 99 of file PFClusterAlgo.h.

References p1, and posCalcP1_.

Referenced by PFRootEventManager::readOptions().

{ posCalcP1_ = p1; }
void PFClusterAlgo::setS4S1CleanBarrel ( const std::vector< double > &  coeffs) [inline]

Definition at line 70 of file PFClusterAlgo.h.

References minS4S1Barrel_.

Referenced by PFRootEventManager::readOptions().

{minS4S1Barrel_ = coeffs;}
void PFClusterAlgo::setS4S1CleanEndcap ( const std::vector< double > &  coeffs) [inline]

Definition at line 86 of file PFClusterAlgo.h.

References minS4S1Endcap_.

Referenced by PFRootEventManager::readOptions().

{minS4S1Endcap_ = coeffs;}
void PFClusterAlgo::setS6S2DoubleSpikeBarrel ( double  cut) [inline]
void PFClusterAlgo::setS6S2DoubleSpikeEndcap ( double  cut) [inline]
void PFClusterAlgo::setShowerSigma ( double  sigma) [inline]

set shower sigma for

Definition at line 105 of file PFClusterAlgo.h.

References showerSigma_.

Referenced by PFRootEventManager::readOptions().

{ showerSigma_ = sigma;}
void PFClusterAlgo::setThreshBarrel ( double  thresh) [inline]

setters -------------------------------------------------------

set barrel threshold

Definition at line 61 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshBarrel_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshCleanBarrel ( double  thresh) [inline]

set barrel clean threshold

Definition at line 69 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshCleanBarrel_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshCleanEndcap ( double  thresh) [inline]

set endcap clean threshold

Definition at line 85 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshCleanEndcap_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshDoubleSpikeBarrel ( double  thresh) [inline]

set endcap thresholds for double spike cleaning

Definition at line 73 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshDoubleSpikeBarrel_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshDoubleSpikeEndcap ( double  thresh) [inline]

set endcap thresholds for double spike cleaning

Definition at line 89 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshDoubleSpikeEndcap_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshEndcap ( double  thresh) [inline]

set endcap threshold

Definition at line 77 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshEndcap_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshPtBarrel ( double  thresh) [inline]
void PFClusterAlgo::setThreshPtEndcap ( double  thresh) [inline]
void PFClusterAlgo::setThreshPtSeedBarrel ( double  thresh) [inline]
void PFClusterAlgo::setThreshPtSeedEndcap ( double  thresh) [inline]
void PFClusterAlgo::setThreshSeedBarrel ( double  thresh) [inline]

set barrel seed threshold

Definition at line 65 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshSeedBarrel_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setThreshSeedEndcap ( double  thresh) [inline]

set endcap seed threshold

Definition at line 81 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshSeedEndcap_.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setUseCornerCells ( bool  usecornercells) [inline]

activate use of cells with a common corner to build topo-clusters

Definition at line 108 of file PFClusterAlgo.h.

References useCornerCells_.

Referenced by PFRootEventManager::readOptions().

{ useCornerCells_ = usecornercells;}
double PFClusterAlgo::showerSigma ( ) const [inline]

get shower sigma

Definition at line 139 of file PFClusterAlgo.h.

References showerSigma_.

{ return showerSigma_ ;}
double PFClusterAlgo::threshBarrel ( ) const [inline]

getters -------------------------------------------------------

get barrel threshold

Definition at line 116 of file PFClusterAlgo.h.

References threshBarrel_.

Referenced by DisplayManager::createGRecHit().

{return threshBarrel_;}
double PFClusterAlgo::threshEndcap ( ) const [inline]

get endcap threshold

Definition at line 123 of file PFClusterAlgo.h.

References threshEndcap_.

Referenced by DisplayManager::createGRecHit().

{return threshEndcap_;}
double PFClusterAlgo::threshSeedBarrel ( ) const [inline]

get barrel seed threshold

Definition at line 119 of file PFClusterAlgo.h.

References threshSeedBarrel_.

double PFClusterAlgo::threshSeedEndcap ( ) const [inline]

get endcap seed threshold

Definition at line 126 of file PFClusterAlgo.h.

References threshSeedEndcap_.

void PFClusterAlgo::write ( void  )

write histos

Definition at line 55 of file PFClusterAlgo.cc.

References gather_cfg::cout, and file_.

Referenced by PFRootEventManager::write().

                     { 

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

}

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  out,
const PFClusterAlgo algo 
) [friend]

Member Data Documentation

option to clean HCAL RBX's and HPD's

Definition at line 330 of file PFClusterAlgo.h.

Referenced by doClusteringWorker(), and setCleanRBXandHPDs().

std::vector< unsigned > PFClusterAlgo::color_ [private]

color, for all rechits

Definition at line 259 of file PFClusterAlgo.h.

Referenced by color(), doClusteringWorker(), and paint().

bool PFClusterAlgo::debug_ [private]

debugging on/off

Definition at line 333 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), buildTopoCluster(), buildTopoClusters(), enableDebugging(), and findSeeds().

std::multimap<double, unsigned, std::greater<double> > PFClusterAlgo::eRecHits_ [private]

indices to rechits, sorted by decreasing E (not E_T)

Definition at line 252 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClusteringWorker(), and findSeeds().

TFile* PFClusterAlgo::file_ [private]

Definition at line 342 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), setHistos(), and write().

TH2F* PFClusterAlgo::hBNeighbour [private]

Definition at line 340 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), and setHistos().

TH2F* PFClusterAlgo::hENeighbour [private]

Definition at line 341 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), and setHistos().

std::set<unsigned> PFClusterAlgo::idUsedRecHits_ [private]

ids of rechits used in seed search

Definition at line 249 of file PFClusterAlgo.h.

std::vector< bool > PFClusterAlgo::mask_ [private]

mask, for all rechits. only masked rechits will be clustered by default, all rechits are masked. see setMask function.

Definition at line 256 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClustering(), findSeeds(), and masked().

std::vector<double> PFClusterAlgo::minS4S1Barrel_ [private]

Definition at line 300 of file PFClusterAlgo.h.

Referenced by operator<<(), and setS4S1CleanBarrel().

std::vector<double> PFClusterAlgo::minS4S1Endcap_ [private]

Definition at line 308 of file PFClusterAlgo.h.

Referenced by operator<<(), and setS4S1CleanEndcap().

Definition at line 304 of file PFClusterAlgo.h.

Referenced by setS6S2DoubleSpikeBarrel().

Definition at line 312 of file PFClusterAlgo.h.

Referenced by setS6S2DoubleSpikeEndcap().

number of neighbours

Definition at line 315 of file PFClusterAlgo.h.

Referenced by findSeeds(), nNeighbours(), operator<<(), and setNNeighbours().

std::auto_ptr< std::vector<reco::PFCluster> > PFClusterAlgo::pfClusters_ [private]

all clusters

particle flow clusters

Definition at line 277 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), clusters(), doClusteringWorker(), and operator<<().

std::auto_ptr< std::vector<reco::PFRecHit> > PFClusterAlgo::pfRecHitsCleaned_ [private]

particle flow rechits cleaned

Definition at line 280 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClusteringWorker(), findSeeds(), and rechitsCleaned().

number of crystals for position calculation

Definition at line 318 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), calculateClusterPosition(), operator<<(), posCalcNCrystal(), and setPosCalcNCrystal().

double PFClusterAlgo::posCalcP1_ [private]

parameter for position calculation

Definition at line 321 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), posCalcP1(), and setPosCalcP1().

unsigned PFClusterAlgo::prodNum_ = 1 [static, private]

product number

Definition at line 337 of file PFClusterAlgo.h.

Definition at line 246 of file PFClusterAlgo.h.

Referenced by createRecHitRef(), and doClustering().

std::vector< unsigned > PFClusterAlgo::seeds_ [private]

vector of indices for seeds.

Definition at line 268 of file PFClusterAlgo.h.

Referenced by buildTopoClusters(), and findSeeds().

std::vector< SeedState > PFClusterAlgo::seedStates_ [private]

seed state, for all rechits

Definition at line 262 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), doClusteringWorker(), findSeeds(), and isSeed().

double PFClusterAlgo::showerSigma_ [private]

sigma of shower (cm)

Definition at line 324 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), operator<<(), setShowerSigma(), and showerSigma().

double PFClusterAlgo::threshBarrel_ [private]

barrel threshold

Definition at line 283 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), setThreshBarrel(), and threshBarrel().

Barrel cleaning threshold and S4/S1 smallest fractiom.

Definition at line 299 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshCleanBarrel().

Endcap cleaning threshold and S4/S1 smallest fractiom.

Definition at line 307 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshCleanEndcap().

Barrel double-spike cleaning.

Definition at line 303 of file PFClusterAlgo.h.

Referenced by setThreshDoubleSpikeBarrel().

Endcap double-spike cleaning.

Definition at line 311 of file PFClusterAlgo.h.

Referenced by setThreshDoubleSpikeEndcap().

double PFClusterAlgo::threshEndcap_ [private]

endcap threshold

Definition at line 291 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), setThreshEndcap(), and threshEndcap().

Definition at line 284 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshPtBarrel().

Definition at line 292 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshPtEndcap().

Definition at line 288 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshPtSeedBarrel().

Definition at line 296 of file PFClusterAlgo.h.

Referenced by operator<<(), and setThreshPtSeedEndcap().

barrel seed threshold

Definition at line 287 of file PFClusterAlgo.h.

Referenced by operator<<(), setThreshSeedBarrel(), and threshSeedBarrel().

endcap seed threshold

Definition at line 295 of file PFClusterAlgo.h.

Referenced by operator<<(), setThreshSeedEndcap(), and threshSeedEndcap().

std::vector< std::vector< unsigned > > PFClusterAlgo::topoClusters_ [private]

sets of cells having one common side, and energy over threshold

Definition at line 271 of file PFClusterAlgo.h.

Referenced by buildTopoClusters(), and doClusteringWorker().

option to use cells with a common corner to build topo-clusters

Definition at line 327 of file PFClusterAlgo.h.

Referenced by buildTopoCluster(), operator<<(), and setUseCornerCells().

std::vector< bool > PFClusterAlgo::usedInTopo_ [private]

used in topo cluster? for all rechits

Definition at line 265 of file PFClusterAlgo.h.

Referenced by buildTopoCluster(), buildTopoClusters(), and doClusteringWorker().