Page MenuHomeHEPForge
Files F10664391

CellJet.cc
No OneTemporary
Actions

CellJet.cc
View Options

#include "CellJet.hh"
namespace SpartyJet {
namespace Pythia8 {
//**************************************************************************
// CellJet class.
// This class performs a cone jet search in (eta, phi, E_T) space.
//*********
// Analyze event.
bool CellJet::analyze(CellJetEvent& event, double eTjetMinIn,
double coneRadiusIn, double eTseedIn) {
bool quiet = true;
if(!quiet)cout << "eTjetMinIn: " << eTjetMinIn << endl;
if(!quiet)cout << "coneRadiusIn: " << coneRadiusIn << endl;
if(!quiet)cout << "eTseedIn: " << eTseedIn << endl;
if(!quiet)cout << "smear: " << smear << endl;
// Input values. Initial values zero.
eTjetMin = eTjetMinIn;
coneRadius = coneRadiusIn;
eTseed = eTseedIn;
jets.resize(0);
vector<SingleCell> cells;
// Loop over desired particles in the event.
for (int i = 0; i < event.size(); ++i) {
// if (event[i].isFinal()) {
// if (select > 2 && event[i].isNeutral() ) continue;
// if (select == 2 && !event[i].isVisible() ) continue;
// Find particle position in (eta, phi, pT) space.
double etaNow = event[i].eta();
if (abs(etaNow) > etaMax) continue;
double phiNow = event[i].phi();
double pTnow = event[i].pT();
int iEtaNow = max(1, min( nEta, 1 + int(nEta * 0.5
* (1. + etaNow / etaMax) ) ) );
int iPhiNow = max(1, min( nPhi, 1 + int(nPhi * 0.5
* (1. + phiNow / M_PI) ) ) );
int iCell = nPhi * iEtaNow + iPhiNow;
// Add pT to cell already hit or book a new cell.
bool found = false;
for (int j = 0; j < int(cells.size()); ++j) {
if (iCell == cells[j].iCell) {
found = true;
++cells[j].multiplicity;
cells[j].eTcell += pTnow;
if(!quiet)cout << "found " << j << " multiplicity: " << cells[j].multiplicity << " eTcell: "<< cells[j].eTcell << endl;
continue;
}
}
if (!found) {
double etaCell = (etaMax / nEta) * (2 * iEtaNow - 1 - nEta);
double phiCell = (M_PI / nPhi) * (2 * iPhiNow - 1 - nPhi);
cells.push_back( SingleCell( iCell, etaCell, phiCell, pTnow, 1) );
if(!quiet)cout << "not found, cells.size: " << cells.size() << " new cell: " << iCell << " " << etaCell << " " << phiCell << pTnow << endl;
}
}
// Smear true bin content by calorimeter resolution.
if (smear > 0)
for (int j = 0; j < int(cells.size()); ++j) {
double eTeConv = (smear < 2) ? 1. : cosh( cells[j].etaCell );
double eBef = cells[j].eTcell * eTeConv;
double eAft = 0.;
do eAft = eBef + resolution * sqrt(eBef) * Rndm::gauss();
while (eAft < 0 || eAft > upperCut * eBef);
cells[j].eTcell = eAft / eTeConv;
}
// Remove cells below threshold for seed or for use at all.
for (int j = 0; j < int(cells.size()); ++j) {
if (cells[j].eTcell < eTseed) cells[j].canBeSeed = false;
if (cells[j].eTcell < threshold) cells[j].isUsed = true;
}
if(!quiet)cout << "removed cells below threshold, now cells.size = " << cells.size() << endl;
// Find seed cell: the one with highest pT of not yet probed ones.
for ( ; ; ) {
int jMax = 0;
double eTmax = 0.;
for (int j = 0; j < int(cells.size()); ++j)
if (cells[j].canBeSeed && cells[j].eTcell > eTmax) {
jMax = j;
eTmax = cells[j].eTcell;
}
if(!quiet)cout << "max found: " << jMax << " " << eTmax << endl;
if(!quiet)cout << "comparing: " << eTmax << " " << eTseed << endl;
// If too small cell eT then done, else start new trial jet.
if (eTmax < eTseed) break;
double etaCenter = cells[jMax].etaCell;
double phiCenter = cells[jMax].phiCell;
double eTjet = 0.;
// Sum up unused cells within required distance of seed.
for (int j = 0; j < int(cells.size()); ++j) {
if (cells[j].isUsed) continue;
double dEta = abs( cells[j].etaCell - etaCenter );
if (dEta > coneRadius) continue;
double dPhi = abs( cells[j].phiCell - phiCenter );
if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
if (dPhi > coneRadius) continue;
if (pow2(dEta) + pow2(dPhi) > pow2(coneRadius)) continue;
cells[j].isAssigned = true;
eTjet += cells[j].eTcell;
if(!quiet)cout << "assigning cell " << j << " to jet, eTjet now " << eTjet << endl;
}
// Reject cluster below minimum ET.
if(!quiet)cout << "comparing eTjet and eTjetMin " << eTjet << " " << eTjetMin << endl;
if (eTjet < eTjetMin) {
if(!quiet)cout << "jet eT too small, not adding jet " << endl;
cells[jMax].canBeSeed = false;
for (int j = 0; j < int(cells.size()); ++j)
cells[j].isAssigned = false;
// Else find new jet properties.
} else {
if(!quiet) cout << "adding new jet" << endl;
double etaWeighted = 0.;
double phiWeighted = 0.;
int multiplicity = 0;
Vec4 pMassive;
for (int j = 0; j < int(cells.size()); ++j)
if (cells[j].isAssigned) {
cells[j].canBeSeed = false;
cells[j].isUsed = true;
cells[j].isAssigned = false;
etaWeighted += cells[j].eTcell * cells[j].etaCell;
double phiCell = cells[j].phiCell;
if (abs(phiCell - phiCenter) > M_PI)
phiCell += (phiCenter > 0.) ? 2. * M_PI : -2. * M_PI;
phiWeighted += cells[j].eTcell * phiCell;
multiplicity += cells[j].multiplicity;
pMassive += cells[j].eTcell * Vec4( cos(cells[j].phiCell),
sin(cells[j].phiCell), sinh(cells[j].etaCell),
cosh(cells[j].etaCell) );
}
etaWeighted /= eTjet;
phiWeighted /= eTjet;
// Bookkeep new jet, in decreasing ET order.
jets.push_back( SingleCellJet( eTjet, etaCenter, phiCenter,
etaWeighted, phiWeighted, multiplicity, pMassive) );
for (int i = int(jets.size()) - 1; i > 0; --i) {
if (jets[i-1].eTjet > jets[i].eTjet) break;
swap( jets[i-1], jets[i]);
}
}
}
// Done.
return true;
//}
}
//*********
// Provide a listing of the info.
void CellJet::list(ostream& os) {
// Header.
os << "\n -------- PYTHIA CellJet Listing, eTjetMin = "
<< fixed << setprecision(3) << setw(8) << eTjetMin
<< ", coneRadius = " << setw(5) << coneRadius
<< " ------------------------------ \n \n no "
<< " eTjet etaCtr phiCtr etaWt phiWt mult p_x"
<< " p_y p_z e m \n";
// The jets.
for (int i = 0; i < int(jets.size()); ++i) {
os << setw(4) << i << setw(10) << jets[i].eTjet << setw(8)
<< jets[i].etaCenter << setw(8) << jets[i].phiCenter << setw(8)
<< jets[i].etaWeighted << setw(8) << jets[i].phiWeighted
<< setw(5) << jets[i].multiplicity << setw(11)
<< jets[i].pMassive.px() << setw(11) << jets[i].pMassive.py()
<< setw(11) << jets[i].pMassive.pz() << setw(11)
<< jets[i].pMassive.e() << setw(11)
<< jets[i].pMassive.mCalc() << "\n";
}
// Listing finished.
os << "\n -------- End PYTHIA CellJet Listing ------------------"
<< "-------------------------------------------------"
<< endl;
}
//**************************************************************************
} // namespace SpartyJet
} // end namespace Pythia8

File Metadata

Mime Type
text/x-c
Expires
Thu, Apr 24, 6:38 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4802592
Default Alt Text
CellJet.cc (7 KB)

Event Timeline