Skip to content

Commit

Permalink
Merge pull request #139 from osahin/DBSCAN-3DCluster
Browse files Browse the repository at this point in the history 
DBSCAN 3D cluster algorithm
  • Loading branch information
@jbsauvan
jbsauvan committed Sep 1, 2017
2 parents f1d2669 + b5cb13b commit 7947a82
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 10 deletions.
14 changes: 13 additions & 1 deletion L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h
View file
Expand Up @@ -18,15 +18,27 @@ class HGCalMulticlusteringImpl{
const l1t::HGCalMulticluster & mclu,
double dR ) const;

void clusterize( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr,
void clusterizeDR( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr,
l1t::HGCalMulticlusterBxCollection & multiclusters);

void clusterizeDBSCAN( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr,
l1t::HGCalMulticlusterBxCollection & multiclusters);

private:

void findNeighbor( const edm::PtrVector<l1t::HGCalCluster> & clustersPtrs,
const l1t::HGCalCluster & cluster,
std::vector<int> & neighborList);

bool isNeighbor( const l1t::HGCalCluster & clu1,
const l1t::HGCalCluster & clu2) const;

double dr_;
double ptC3dThreshold_;
double calibSF_;
string multiclusterAlgoType_;
double distDbscan_ = 0.03;
unsigned minNDbscan_ = 3;

HGCalShowerShape shape_;

Expand Down
27 changes: 26 additions & 1 deletion L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc
View file
Expand Up @@ -28,6 +28,10 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
dRC2d,
NNC2d
};
enum MulticlusterType{
dRC3d,
DBSCANC3d
};

public:

Expand All @@ -52,6 +56,16 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
<< "'. Using nearest neighbor NNC2d instead.\n";
clusteringAlgorithmType_ = NNC2d;
}
std::string typeMulticluster(conf.getParameterSet("C3d_parameters").getParameter<std::string>("type_multicluster"));
if(typeMulticluster=="dRC3d"){
multiclusteringAlgoType_ = dRC3d;
}else if(typeMulticluster=="DBSCANC3d"){
multiclusteringAlgoType_ = DBSCANC3d;
}else {
edm::LogWarning("ParameterError") << "Unknown Multiclustering type '" << typeMulticluster
<< "'. Using Cone Algorithm instead.\n";
multiclusteringAlgoType_ = dRC3d;
}

}

Expand Down Expand Up @@ -97,6 +111,7 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
ClusterType clusteringAlgorithmType_;
double clustering_threshold_silicon_;
double clustering_threshold_scintillator_;
MulticlusterType multiclusteringAlgoType_;
};


Expand Down Expand Up @@ -170,7 +185,17 @@ void HGCClusterAlgo<FECODEC,DATA>::run(const l1t::HGCFETriggerDigiCollection & c
}

/* call to multiclustering and compute shower shape*/
multiclustering_.clusterize( clustersPtrs, *multicluster_product_ );
switch(multiclusteringAlgoType_){
case dRC3d :
multiclustering_.clusterizeDR( clustersPtrs, *multicluster_product_ );
break;
case DBSCANC3d:
multiclustering_.clusterizeDBSCAN( clustersPtrs, *multicluster_product_ );
break;
default:
// Should not happen, clustering type checked in constructor
break;
}

/* retrieve the orphan handle to the multiclusters collection and put the collection in the event */
multiclustersHandle = evt.put( std::move( multicluster_product_ ), "cluster3D");
Expand Down
6 changes: 5 additions & 1 deletion L1Trigger/L1THGCal/python/hgcalTriggerPrimitiveDigiProducer_cfi.py
View file
Expand Up @@ -65,7 +65,11 @@

C3d_parValues = cms.PSet( dR_multicluster = cms.double(0.01), # dR in normalized plane used to clusterize C2d
minPt_multicluster = cms.double(0.5), # minimum pt of the multicluster (GeV)
calibSF_multicluster = cms.double(1.084)
calibSF_multicluster = cms.double(1.084),
type_multicluster = cms.string('dRC3d'), #'DBSCANC3d' for the DBSCAN algorithm
dist_dbscan_multicluster = cms.double(0.01),
minN_dbscan_multicluster = cms.uint32(2)

)
cluster_algo = cms.PSet( AlgorithmName = cms.string('HGCClusterAlgoThreshold'),
FECodec = fe_codec.clone(),
Expand Down
125 changes: 118 additions & 7 deletions L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc
View file
Expand Up @@ -6,12 +6,18 @@
HGCalMulticlusteringImpl::HGCalMulticlusteringImpl( const edm::ParameterSet& conf ) :
dr_(conf.getParameter<double>("dR_multicluster")),
ptC3dThreshold_(conf.getParameter<double>("minPt_multicluster")),
calibSF_(conf.getParameter<double>("calibSF_multicluster"))
calibSF_(conf.getParameter<double>("calibSF_multicluster")),
multiclusterAlgoType_(conf.getParameter<string>("type_multicluster")),
distDbscan_(conf.getParameter<double>("dist_dbscan_multicluster")),
minNDbscan_(conf.getParameter<unsigned>("minN_dbscan_multicluster"))
{
edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster dR for Near Neighbour search: " << dr_;
edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster minimum transverse-momentum: " << ptC3dThreshold_;
edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster global calibration factor: " << calibSF_;

edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster DBSCAN Clustering distance: " << distDbscan_;
edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster clustering min number of subclusters: " << minNDbscan_;
edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster type of multiclustering algortihm: " << multiclusterAlgoType_;

}


Expand All @@ -33,10 +39,43 @@ bool HGCalMulticlusteringImpl::isPertinent( const l1t::HGCalCluster & clu,
}


void HGCalMulticlusteringImpl::clusterize( const edm::PtrVector<l1t::HGCalCluster> & clustersPtrs,
bool HGCalMulticlusteringImpl::isNeighbor( const l1t::HGCalCluster & clu1,
const l1t::HGCalCluster & clu2) const
{
HGCalDetId cluDetId( clu2.detId() );
HGCalDetId firstClusterDetId( clu1.detId() );

if( cluDetId.zside() != firstClusterDetId.zside() ){
return false;
}
double dist = (clu1.centreProj() - clu2.centreProj() ).mag();

if(dist < distDbscan_ ) {
return true;
}
return false;
}


void HGCalMulticlusteringImpl::findNeighbor( const edm::PtrVector<l1t::HGCalCluster> & clustersPtrs,
const l1t::HGCalCluster & cluster,
std::vector<int> & neighbors
)
{
int iclu = 0;

for(edm::PtrVector<l1t::HGCalCluster>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end(); ++clu, ++iclu){

if(isNeighbor(cluster, **clu)){
neighbors.push_back(iclu);
}
}
}

void HGCalMulticlusteringImpl::clusterizeDR( const edm::PtrVector<l1t::HGCalCluster> & clustersPtrs,
l1t::HGCalMulticlusterBxCollection & multiclusters)
{

std::vector<l1t::HGCalMulticluster> multiclustersTmp;

int iclu = 0;
Expand Down Expand Up @@ -95,6 +134,78 @@ void HGCalMulticlusteringImpl::clusterize( const edm::PtrVector<l1t::HGCalCluste
}

}



void HGCalMulticlusteringImpl::clusterizeDBSCAN( const edm::PtrVector<l1t::HGCalCluster> & clustersPtrs,
l1t::HGCalMulticlusterBxCollection & multiclusters)
{

std::vector<l1t::HGCalMulticluster> multiclustersTmp;
l1t::HGCalMulticluster mcluTmp;
std::vector<bool> visited(clustersPtrs.size(),false);
std::vector<bool> merged (clustersPtrs.size(),false);
std::vector<std::vector<int>> neighborList;
int iclu = 0, imclu = 0, neighNo = 0;

for(edm::PtrVector<l1t::HGCalCluster>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end(); ++clu, ++iclu){
std::vector<int> neighbors;

if(!visited.at(iclu)){
visited.at(iclu)=true;
findNeighbor(clustersPtrs, **clu, neighbors);
neighborList.push_back(std::move(neighbors));

if(neighborList.at(iclu).size() > minNDbscan_) {
multiclustersTmp.emplace_back( *clu );
merged.at(iclu) = true;
/* dynamic range loop: range-based loop syntax cannot be employed */
for(unsigned int neighInd = 0; neighInd < neighborList.at(iclu).size(); neighInd++){

neighNo = neighborList.at(iclu).at(neighInd);

if(!visited.at(neighNo)){
visited.at(neighNo) = true;
std::vector<int> secNeighbors;
findNeighbor(clustersPtrs,*(clustersPtrs[neighNo]), secNeighbors);
multiclustersTmp.at(imclu).addConstituent( clustersPtrs[neighNo]);
merged.at(neighNo) = true;

if(secNeighbors.size() > minNDbscan_){
neighborList.at(iclu).insert(neighborList.at(iclu).end(), secNeighbors.begin(), secNeighbors.end());
}

} else if(!merged.at(neighNo) ){
merged.at(neighNo) = true;
multiclustersTmp.at(imclu).addConstituent( clustersPtrs[neighNo] );
}
}
imclu++;
}
}

else neighborList.push_back(std::move(neighbors));
}

/* making the collection of multiclusters */
for( unsigned i(0); i<multiclustersTmp.size(); ++i ){
math::PtEtaPhiMLorentzVector calibP4( multiclustersTmp.at(i).pt() * calibSF_,
multiclustersTmp.at(i).eta(),
multiclustersTmp.at(i).phi(),
0. );
// overwriting the 4p with the calibrated 4p
multiclustersTmp.at(i).setP4( calibP4 );

if( multiclustersTmp.at(i).pt() > ptC3dThreshold_ ){

//compute shower shape
multiclustersTmp.at(i).set_showerLength(shape_.showerLength(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_firstLayer(shape_.firstLayer(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_sigmaEtaEtaTot(shape_.sigmaEtaEtaTot(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_sigmaEtaEtaMax(shape_.sigmaEtaEtaMax(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_sigmaZZ(shape_.sigmaZZ(multiclustersTmp.at(i)));
multiclustersTmp.at(i).set_eMax(shape_.eMax(multiclustersTmp.at(i)));

multiclusters.push_back( 0, multiclustersTmp.at(i));
}
}
}

0 comments on commit 7947a82

Please sign in to comment.