Page MenuHomeHEPForge
Files F10664295

SISConeFinder.cc
No OneTemporary
Actions

SISConeFinder.cc
View Options

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ActiveAreaSpec.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/ClusterSequenceActiveArea.hh"
#include "fastjet/ClusterSequencePassiveArea.hh"
//#include "SISCone/SISConePlugin.hh"
#include "fastjet/SISConePlugin.hh"
#include "JetCore/JetMomentMap.hh"
#include "JetCore/CommonUtils.hh"
#include <vector>
#include "SISConeFinder.hh"
// Declare new constructor for PseudoJet... implementation is in
// FastJetFinder ... not very clean
namespace fastjet {
template<>
PseudoJet::PseudoJet(const SpartyJet::Jet &v);
};
namespace SpartyJet {
namespace fastjet {
SISConeFinder::SISConeFinder(char* name,bool a) {
area = a;
areaChoice = 1;
inclusive = true;
exclusive = false;
coneRadius = 0.7;
overlapThreshold = 0.75;
nPassMax = 0;
protojetptmin = 0.0;
caching = 0;
ghost_etamax = 6.0;
repeat = 5;
ghost_area = 0.1;
grid_scatter = 1E-4;
kt_scatter = 0.01;
mean_ghost_kt = 1E-100;
saveSplitMerged = false;
m_name = (std::string)name;
init();
}
SISConeFinder::SISConeFinder(bool a) {
area = a;
areaChoice = 1;
inclusive = true;
exclusive = false;
coneRadius = 0.7;
overlapThreshold = 0.75;
nPassMax = 0;
protojetptmin = 0.0;
caching = 0;
ghost_etamax = 6.0;
repeat = 5;
ghost_area = 0.1;
grid_scatter = 1E-4;
kt_scatter = 0.01;
mean_ghost_kt = 1E-100;
saveSplitMerged = false;
m_name = "mySISCone";
init();
}
SISConeFinder::SISConeFinder(std::string name,bool a) {
area = a;
areaChoice = 1;
inclusive = true;
exclusive = false;
coneRadius = 0.7;
overlapThreshold = 0.75;
nPassMax = 0;
protojetptmin = 0.0;
caching = 0;
ghost_etamax = 6.0;
repeat = 5;
ghost_area = 0.1;
grid_scatter = 1E-4;
kt_scatter = 0.01;
mean_ghost_kt = 1E-100;
saveSplitMerged = false;
m_name = name;
init();
}
void SISConeFinder::init(JetMomentMap *mmap) {
if(mmap != NULL) {
// if(saveSplitMerged) mmap->schedule_jet_moment("split_merged");
if(area) {
mmap->schedule_jet_moment("area");
mmap->schedule_jet_moment("area_error");
}
}
}
void SISConeFinder::configure(bool a,bool in, bool ex,double cr,double ot,int n,double pt,int c) {
area = a;
inclusive = in;
exclusive = ex;
coneRadius = cr;
overlapThreshold = ot;
nPassMax = n;
protojetptmin = pt;
caching = c;
init();
}
void SISConeFinder::configure_area(double e,int r,double a,
double s,double ks,double kt) {
ghost_etamax = e;
repeat = r;
ghost_area = a;
grid_scatter = s;
kt_scatter = ks;
mean_ghost_kt = kt;
}
} // namespace fastjet
} // namespace SpartyJet
// To avoid conflict with original fastjet namespace
namespace fastjet_notSJ = fastjet ;
void SpartyJet::fastjet::SISConeFinder::execute(SpartyJet::JetCollection &inputJets) {
JetCollection::iterator iter = inputJets.begin();
// Prevent problems with empty events :
if (inputJets.size() == 0) return;
// convert jetlist to clusters and clear inputJets
std::vector<fastjet_notSJ::PseudoJet> towers;
towers.reserve(inputJets.size()); // we know what will be the size
Jet *tjet = NULL;
while(iter != inputJets.end()){
tjet = *iter;
towers.push_back(fastjet_notSJ::PseudoJet(*(*iter)));
++iter;
}
// run the algorithm
fastjet_notSJ::JetDefinition::Plugin * plugin;
plugin = new fastjet_notSJ::SISConePlugin(coneRadius,overlapThreshold,nPassMax,protojetptmin,caching);
fastjet_notSJ::JetDefinition jet_def(plugin);
// retrieve the map of the collection :
JetMomentMap * themap = inputJets.get_JetMomentMap();
// if (saveSplitMerged && (themap->num_jet_moment() ==0)) return;
if (area && (themap->num_jet_moment() ==0)) return;
int areapos = 0;
int area_errorpos = 0;
//int split_merged_pos = 0;
//if(saveSplitMerged) split_merged_pos = themap->get_jet_momentPos("split_merged");
fastjet_notSJ::ClusterSequence * clust_seq = 0;
fastjet_notSJ::ClusterSequenceArea * clust_seq_area = 0;
if(!area) {
clust_seq = new fastjet_notSJ::ClusterSequence(towers,jet_def);
}
else {
fastjet_notSJ::AreaType aType;
fastjet_notSJ::GhostedAreaSpec Garea_spec;
fastjet_notSJ::VoronoiAreaSpec Varea_spec;
fastjet_notSJ::AreaDefinition area_def;
switch (areaChoice)
{
case 1:
Garea_spec = fastjet_notSJ::GhostedAreaSpec(ghost_etamax,repeat,ghost_area,grid_scatter,kt_scatter,mean_ghost_kt);
aType = fastjet_notSJ::active_area;
area_def = fastjet_notSJ::AreaDefinition(Garea_spec,aType);
break;
case 2:
Garea_spec = fastjet_notSJ::GhostedAreaSpec(ghost_etamax,repeat,ghost_area,grid_scatter,kt_scatter,mean_ghost_kt);
aType = fastjet_notSJ::active_area_explicit_ghosts;
area_def = fastjet_notSJ::AreaDefinition(Garea_spec,aType);
break;
case 3:
Garea_spec = fastjet_notSJ::GhostedAreaSpec(ghost_etamax,repeat,ghost_area,grid_scatter,kt_scatter,mean_ghost_kt);
aType = fastjet_notSJ::passive_area;
area_def = fastjet_notSJ::AreaDefinition(Garea_spec,aType);
break;
case 4:
Garea_spec = fastjet_notSJ::GhostedAreaSpec(ghost_etamax,repeat,ghost_area,grid_scatter,kt_scatter,mean_ghost_kt);
aType = fastjet_notSJ::one_ghost_passive_area;
area_def = fastjet_notSJ::AreaDefinition(Garea_spec,aType);
break;
case 5:
area_def = fastjet_notSJ::AreaDefinition(Varea_spec);
break;
}
clust_seq_area = new fastjet_notSJ::ClusterSequenceArea(towers, jet_def,area_def);
areapos = themap->get_jet_momentPos("area");
area_errorpos = themap->get_jet_momentPos("area_error");
}
std::vector<fastjet_notSJ::PseudoJet> output_jets;
std::vector<fastjet_notSJ::PseudoJet> temp_jets;
if(inclusive && !exclusive) { // only inclusive
if (!area) output_jets = sorted_by_pt(clust_seq->inclusive_jets(protojetptmin));
else output_jets = sorted_by_pt(clust_seq_area->inclusive_jets(protojetptmin));
}
if(exclusive && !inclusive){ // only exclusive
if (!area) output_jets = sorted_by_pt(clust_seq->exclusive_jets(protojetptmin));
else output_jets = sorted_by_pt(clust_seq_area->exclusive_jets(protojetptmin));
}
else if(exclusive && inclusive){ // both inclusive and exclusive
if (!area) temp_jets = sorted_by_pt(clust_seq->inclusive_jets(protojetptmin));
else temp_jets = sorted_by_pt(clust_seq_area->inclusive_jets(protojetptmin));
for(unsigned i = 0; i < temp_jets.size(); i++)
output_jets.push_back(temp_jets[i]);
temp_jets.clear();
if (!area) temp_jets = sorted_by_pt(clust_seq->exclusive_jets(protojetptmin));
else temp_jets = sorted_by_pt(clust_seq_area->exclusive_jets(protojetptmin));
for(unsigned i = 0; i < temp_jets.size(); i++)
output_jets.push_back(temp_jets[i]);
}
JetCollection tmp_list;
for(unsigned i = 0; i < output_jets.size(); i++) {
fastjet_notSJ::PseudoJet & jet = output_jets[i];
// convert back to jet class
Jet * t2jet = new Jet(jet.px(),jet.py(),jet.pz(),jet.E());
// if(saveSplitMerged) { // doesn't work anymore
// // if(t2jet->split_merged)
// themap->set_jet_moment(split_merged_pos,t2jet, (float)t2jet->split_merged);// split/merged
// // else
// // themap->set_jet_moment(split_merged_pos,t2jet,0); // not split/merged
// }
if(area){
themap->set_jet_moment(areapos,t2jet,(float)clust_seq_area->area(jet)); // setting the area for the momentmap
themap->set_jet_moment(area_errorpos,t2jet,(float)clust_seq_area->area_error(jet));
}
// Retrieve jets' constituents
std::vector<fastjet_notSJ::PseudoJet> constituents;
if (!area) constituents = clust_seq->constituents(jet);
else {
constituents = clust_seq_area->constituents(jet);
}
for(unsigned c=0;c<constituents.size();c++){
int index = constituents[c].user_index();
// get original constituents of this jet (!!! there maybe more
// than 1, following is not 100% safe !!!)
Jet* constit = *(inputJets[index]->firstConstituent());
t2jet->addConstituent_notMoment(constit);
}
tmp_list.push_back(t2jet);
}
// we can now delete the input collection :
clear_list(inputJets);
// and update it :
inputJets = tmp_list;
// clear --------------------------------
if(area){
delete clust_seq_area;
}else{
delete clust_seq;
}
delete plugin;
}

File Metadata

Mime Type
text/x-c
Expires
Thu, Apr 24, 6:35 AM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4869244
Default Alt Text
SISConeFinder.cc (8 KB)

Event Timeline