Page MenuHomeHEPForge
Files F10664419

FastKtHull.cxx
No OneTemporary
Actions

FastKtHull.cxx
View Options

#include "FastKtHull.h"
//#include <cstdlib>
namespace SpartyJet {
namespace FastKtJet {
KtListHull::KtListHull(jetcollection_t *p, KtDistance * ktdist, KtRecom *recom, bool reversed_mode){
m_fKtDist = ktdist;
m_ktRecom = recom;
KtNNOperation *tt = new KtNNOperation();
m_NNoperator = (KtNNOperation*) tt;
m_reversedMode = reversed_mode;
buildKtJetInfo(p);
}
void KtListHull::buildKtJetInfo(jetcollection_t *jetColl){
jetcollection_t::iterator itktv = jetColl->begin();
jetcollection_t::iterator itktv_e = jetColl->end();
//m_ktlist.reserve(jetColl.size());
KtJetInfo_iterator it_list ;
bool first = true;
m_nRemaining = 0;
m_ktlist.clear();
float eTot= 0.0;
// Fill the KtJetInfo list
for( ; itktv != itktv_e; ++itktv){
//Jet *j = (*itktv);
if ( std::isnan( (*itktv)->e()) ) {
continue;
}
KtHullInfo *ki = 0;
try{ // take care of bad input vectors
//std::cout << m_nRemaining << " jet e = "<< (*itktv)->e() << std::endl;
ki = new KtHullInfo( (*itktv)); // object created here will be deleted in mergJet or passed to other objects (throught getMinJet) which must delete them. (ex. KtEvent)
ki->m_hull.push_back( point_t(ki->eta,ki->phi));
}
catch(...){
//MsgStream log( Athena::getMessageSvc(), "FastKt::FastKtUtils" );
std::cout << " Exception thrown for input vector "<< m_nRemaining <<" (E="<< (*itktv)->e() << ", pz="<< (*itktv)->pz() << ") this vector won't be used by jet finder"<< std::endl;
if(ki) delete ki;
continue;
}
m_ktlist.push_back(ki);
// Make sure it_list corresponds to this ki
if(first){
it_list = m_ktlist.begin();
first = false;
} else {
++it_list;
}
// compute individual dist for this jet and store it
float d_beam = (*m_fKtDist)(ki);
ki->pos_in_diList = insertValueAndIter(m_diList, d_beam, it_list);
ki->index = 0 ;
(*it_list)->index = m_nRemaining ;
eTot += (*itktv)->e();
++m_nRemaining;
}
m_ktlist.sort( KtLists::KtJetInfoEtaPhiComp );
// The sorted ktlist is done, we can init the Neighbourhood builder :
m_NNoperator->init(&m_ktlist , m_fKtDist);
// Now build the neigbours
it_list = m_ktlist.begin();
KtJetInfo_iterator it_listE = m_ktlist.end() ;
for( ; it_list != it_listE; ++it_list){
KtJetInfo_iterator nb = m_NNoperator->computeNN(it_list);
(*it_list)->closest_neighb = nb; // assign closest_neigb;
// compute distance associated with neighbour and sort it in the dij map
float dGi = (*m_fKtDist)((*it_list), (*nb) );
(*it_list)->pos_in_dijList = insertValueAndIter(m_dijList, dGi, it_list);
}
//std::cout << " Total E in input = " << eTot << std::endl;
}
void KtListHull::buildKtJetInfoPass2(KtJetInfo_list * jetlist){
m_ktlist.clear();
m_diList.clear();
m_dijList.clear();
m_ktlist = *jetlist;
std::cout << " building for pass 2 "<< jetlist->size() << " , "<< m_ktlist.size() << std::endl;
m_ktlist.sort( KtLists::KtJetInfoEtaPhiComp );
}
// ************************************************************************
// ************************************************************************
KtAlgoHull::KtAlgoHull(double rp): KtAlgoStandard(sqrt(rp)){
m_Rmin=1.3;
m_do2steps=true;
m_ConeMin = (0.1)*(0.1);
m_ConeMax = 2*(0.7)*M_PI;
};
void KtAlgoHull::init(jetcollection_t *constituents) {
if(m_ktList) delete m_ktList;
m_ktList = new KtListHull(constituents, m_ktDist, m_ktRecom, false);
m_N= constituents->size();
m_secondpass = false;
}
bool KtAlgoHull::continueClustering(){
bool docontinue = (m_ktList->getNJets() > 1);
if (!docontinue) {
if(m_secondpass|| !m_do2steps)return false;
if (m_ktList->getNJets() == 1) {
m_extractedJet = m_ktList->getMinJet();
m_finalJets->push_back(m_extractedJet);
m_ktList->killMinJet();
}
//std::cout <<
secondPass();
return false;
}
processStep();
if(m_doExtractJet) {
m_finalJets->push_back(m_extractedJet);
}
return true;
}
void KtAlgoHull::setProperty(std::string name,float val){
if (name =="rmin") m_Rmin =val;
if (name =="conemin") m_ConeMin =val*val;
if (name =="conemax") m_ConeMax =val;
if (name =="2steps"){ m_do2steps = (val>1) ;}
}
void KtAlgoHull::secondPass(){
std::cout << " second pass. finaljets = "<< m_finalJets->size()<< " m_ktlist: "<< m_ktList->getNJets() << std::endl;
// sort by eta
m_finalJets->sort(KtLists::KtJetInfoEtaPhiComp );
// Fill intermediate list
std::list<KtHullInfo*> m_firstpassList;
KtJetInfo_iterator fit = m_finalJets->begin();
KtJetInfo_iterator fitE = m_finalJets->end();
for(; fit!= fitE; ++fit) m_firstpassList.push_back(dynamic_cast<KtHullInfo*>(*fit));
// clear.
m_finalJets->clear();
// Loop on intermediate list
std::list<KtHullInfo*>::iterator it = m_firstpassList.begin();
std::list<KtHullInfo*>::iterator itE = m_firstpassList.end();
float distMax = m_ConeMax/M_PI;
while(it!=itE){
std::list<KtHullInfo*>::iterator next = it;
next++;
if (next== itE) {m_finalJets->push_back(*it); break;}
// else check if 'it' is overalpping :
bool merge=false, mergeback = false;
point_list_t newhull; float phicenter;
float ri = (*it)->getRadius();
do {
float dist = m_ktDist->geodist(*it, *next);
if ( fabs((*it)->eta-(*next)->eta) > distMax) break;
float rn = (*next)->getRadius();
if( (rn>dist) || (ri>dist)){
// jets may overlap :
merge = doMergeJet(*it,*next,newhull,phicenter);
if(merge) break;
}
next++;
} while(next!=m_firstpassList.end());
if((it!=m_firstpassList.begin()) && !merge){
newhull.clear(); phicenter=-10;
// redo the same on left side
next =it;
do {
next--;
float dist = m_ktDist->geodist(*it, *next);
//if (dist > distMax) break;
if ( fabs((*it)->eta-(*next)->eta) > distMax) break;
float rn = (*next)->getRadius();
if( (rn>dist) || (ri>dist)){
// jets may overlap :
merge = doMergeJet(*it,*next,newhull,phicenter);
if(merge) {mergeback=true;break;}
}
} while(next!=m_firstpassList.begin());
}
// Choose between merging or final insertion
if(merge){
(*it)->m_hull.clear();
recenter_set(newhull,(*it)->m_hull,-phicenter);
m_ktRecom->combine(*it,*next);
if(mergeback){ // pointer has been inserted. remove it
m_finalJets->remove(*next);
}
delete *next;
m_firstpassList.erase(next);
//it++;
}else{
m_finalJets->push_back(*it);
it++;
}
}
std::cout << " End second pass. finaljets = "<< m_finalJets->size()<< " m_ktlist: "<< m_ktList->getNJets() << std::endl;
}
bool KtAlgoHull::doMergeJet(KtHullInfo* j1, KtHullInfo* j2, point_list_t &newhull,float &phicenter ){
point_list_t inset(j1->m_hull), inset_recenter;
inset.insert(inset.end(),j2->m_hull.begin(), j2->m_hull.end());
phicenter = getMeanPhi(inset);
bool domerge;
if(phicenter>-9){
recenter_set(inset, inset_recenter,phicenter);
findConvexHull(inset_recenter,newhull);
float l1 = polygon_lenght(j1->m_hull);
float l2 = polygon_lenght(j2->m_hull);
float ln = polygon_lenght(newhull);
//float l1 = 2, l2 =3,ln =3;
float lp = l1 > l2 ? l1 : l2;
//int n=5;
int n = j1->constit_list.size()+j2->constit_list.size();
domerge = lp==0 ? true : ((ln/lp) <= testR(n)) ;
domerge = ((ln < m_ConeMax) && domerge) || (m_lastDPair <= m_ConeMin ) ;
// if(m_secondpass) std::cout <<m_lastIndex1 <<","<<m_lastIndex2 << " || dpair ="<<m_lastDPair << " n1="<< j1->constit_list.size() << " n2="<< j2->constit_list.size() << " | "
// << ln/lp << " testfr=" << testR(n) << " "<< std::endl;
}else{ // can't recenter set
domerge = false;
}
return domerge;
}
void KtAlgoHull::processStep(){
m_lastDPair = m_ktList->getMinDPair();
m_lastDJet = m_ktList->getMinDJet() ;
KtHullInfo *j1 = dynamic_cast<KtHullInfo*>(m_ktList->getMinPairJet());
KtHullInfo *j2 = dynamic_cast<KtHullInfo*>(*j1->closest_neighb);
KtJetInfo_iterator j2_it = j1->closest_neighb;
point_list_t newhull;
float phicenter;
m_doExtractJet = !doMergeJet(j1,j2,newhull,phicenter);
std::pair<int,int> p = m_ktList->getMinPairIndex();
m_lastIndex1 = p.first;
m_lastIndex2 = p.second;
if(m_doExtractJet){
m_extractedJet = j1;
m_ktList->killMinPairFirst();
// extract 2 jets
m_ktList->killJet(j2_it);
m_finalJets->push_back(j2);
}else{
j1->m_hull.clear();
recenter_set(newhull,j1->m_hull,-phicenter);
m_ktList->mergeMinJets();
//std::cout << m_ktList->getNJets() << std::endl;
}
}
void KtAlgoHull::endClustering(){
//std::cout<< " endClustering : "<<m_finalJets->size()+m_ktList->getNJets() << std::endl;
if(m_ktList->getNJets()==0) return;
m_lastDJet = m_ktList->getMinDPair() ;
m_extractedJet = m_ktList->getMinJet();
m_finalJets->push_back(m_extractedJet);
}
} // namespace SpartyJet
};

File Metadata

Mime Type
text/x-c
Expires
Thu, Apr 24, 6:39 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4837523
Default Alt Text
FastKtHull.cxx (9 KB)

Event Timeline