source: trunk/KtJet/KtJetTable.cc@ 16

Last change on this file since 16 was 2, checked in by Pavel Demin, 16 years ago

first commit

File size: 6.6 KB
Line 
1#include "KtJet/KtJetTable.h"
2#include "KtJet/KtLorentzVector.h"
3#include "KtJet/KtDistanceInterface.h"
4#include "KtJet/KtRecomInterface.h"
5#include <vector>
6#include <utility>
7#include <algorithm>
8#include <iostream>
9
10namespace KtJet{
11
12KtJetTable::KtJetTable(const std::vector<KtLorentzVector> & p, KtDistance *ktdist, KtRecom *recom)
13 : m_fKtDist(ktdist), m_ktRecom(recom) {
14 m_nRows = p.size();
15 m_jets.reserve(m_nRows); // Reserve space for jet array
16 m_ddi.reserve(m_nRows); // Reserve space for di array
17 m_dPairs.resize(m_nRows); // Make space for dij table
18
19 std::vector<KtLorentzVector>::const_iterator pitr = p.begin();
20 for (; pitr != p.end(); ++pitr) { // Initialize jets with one particle each
21 KtLorentzVector j; // create new jet
22 j.add(*pitr); // add single particle to jet
23 m_jets.push_back(j);
24 }
25 for (int i = 0; i < m_nRows-1; ++i) { // Fill array of pair kt's
26 for (int j = i+1 ; j < m_nRows; ++j) {
27 KtFloat D = (*m_fKtDist)(m_jets[i],m_jets[j]);
28 m_dPairs(i,j) = D;
29 }
30 }
31 for (int i = 0; i < m_nRows; ++i){ // Fill vector of particle kt's
32 m_ddi.push_back((*m_fKtDist)(m_jets[i]));
33 }
34}
35
36KtJetTable::~KtJetTable() {}
37
38const KtLorentzVector & KtJetTable::getJet(int i) const {
39 return m_jets[i];
40}
41
42KtFloat KtJetTable::getD(int i, int j) const {
43 return m_dPairs(i,j);
44}
45
46KtFloat KtJetTable::getD(int i) const {
47 if (i<0 || i>=static_cast<int>(m_ddi.size())) {
48 std::cout << "ERROR in KtJetTable::getD(int)" << std::endl;
49 std::cout << " i, m_ddi.size() = " << i << ", " << m_ddi.size() << std::endl;
50 }
51 return m_ddi[i];
52}
53
54/***************************************************************
55 * Merge jets i and j, updating four-momentum and kt vectors *
56 ***************************************************************/
57void KtJetTable::mergeJets(int i, int j) {
58 int njet = getNJets();
59 if (i<0 || i>=njet || j<0 || j>=njet || i>=j) {
60 std::cout << "ERROR in KtJetTable::mergeJets" << std::endl;
61 std::cout << " Attempt to merge jets " << i << ", " << j << " in event with " << njet << " jets" << std::endl;
62 }
63 m_jets[i].add(m_jets[j],m_ktRecom); // Add constituents and merge 4-momenta using required scheme
64 for (int ii=0; ii<i; ++ii) { // Calculate pair kt's involving merged particles
65 KtFloat D = (*m_fKtDist)(m_jets[ii],m_jets[i]);
66 m_dPairs(ii,i) = D;
67 }
68 for (int jj=i+1; jj<njet; ++jj) {
69 KtFloat D = (*m_fKtDist)(m_jets[i],m_jets[jj]);
70 m_dPairs(i,jj) = D;
71 }
72 m_ddi[i] = (*m_fKtDist)(m_jets[i]); // Calculate kt of merged particles
73 killJet(j); // Now delete particle j
74}
75
76/***************************************************
77 * Delete jet i by moving last jet on top of it, *
78 * updating four momentum and kt vectors *
79 ***************************************************/
80void KtJetTable::killJet(int i) {
81 // std::cout << " KtJetTable::killJet, i = " << i << std::endl;
82 int njet = getNJets();
83 if (i<0 || i>=njet) {
84 std::cout << "ERROR in KtEvent::killJet" << std::endl;
85 std::cout << " Attempt to delete jet " << i << " in event with " << njet << " jets" << std::endl;
86 }
87 // std::cout << " njet = " << njet << std::endl;
88 m_jets[i] = m_jets[njet-1]; // move last jet into space left by jet i
89 for (int j=0; j<i; ++j) { // move pair kt's
90 m_dPairs(j,i) = m_dPairs(j,njet-1);
91 }
92 for (int j=i+1; j<njet-1; ++j) {
93 m_dPairs(i,j) = m_dPairs(j,njet-1);
94 }
95 m_ddi[i] = m_ddi[njet-1]; // move jet kt
96 m_dPairs.killJet(); // reduce size of kt array
97 m_jets.pop_back(); // delete last jet from vector
98 m_ddi.pop_back(); // delete last jet kt
99 // std::cout << " Finished in KtJetTable::killJet" << std::endl;
100}
101
102std::pair<int,int> KtJetTable::getMinDPair() const {
103 return m_dPairs.getMin();
104}
105
106int KtJetTable::getMinDJet() const {
107 return std::distance(m_ddi.begin(),std::min_element(m_ddi.begin(),m_ddi.end()));
108}
109
110/*************************************************************
111 * Now the functions for nested class KtJetTable::DijTable *
112 *************************************************************/
113
114KtJetTable::DijTable::DijTable(int nParticles) : m_nRows(nParticles), m_nJets(nParticles) {
115 m_table.resize(m_nRows*m_nRows);
116}
117
118KtJetTable::DijTable::~DijTable() {}
119
120void KtJetTable::DijTable::resize(int nParticles) {
121 /*************************************************************
122 * Reserve space for kt of pairs with nParticles particles *
123 *************************************************************/
124 m_nRows = nParticles;
125 m_nJets = nParticles;
126 m_table.resize(m_nRows*m_nRows);
127}
128
129std::pair<int,int> KtJetTable::DijTable::getMin() const {
130 /********************************************************
131 * Find position of smallest entry in table *
132 ********************************************************/
133 KtFloat d = m_table[1];
134 int k=0; int i=0; int j=1; // Initialize to first used element
135 for (int ii=0; ii<m_nJets-1; ++ii) {
136 for (int jj=ii+1; jj<m_nJets; ++jj) {
137 ++k;
138 if (m_table[k]<d) {
139 i = ii; j = jj; d = m_table[k];
140 }
141 }
142 k += 2 + ii + m_nRows - m_nJets;
143 }
144 return std::pair<int,int>(i,j);
145}
146
147void KtJetTable::DijTable::print() const {
148 /*****************************************
149 * Write out contents of table to cout *
150 *****************************************/
151 for (int i = 0; i<m_nRows-1; ++i){
152 for (int j = i+1; j<m_nRows; ++j){
153 KtFloat D = (*this)(i,j);
154 std::cout << i+1 << " " << j+1 << " " << D << '\n';
155 }
156 }
157 std::cout << std::endl;
158}
159
160KtFloat & KtJetTable::DijTable::operator() (int ii, int jj) {
161 int i = std::min(ii,jj);
162 int j = std::max(ii,jj);
163 if (i<0 || j<0 || i>=m_nJets || j>=m_nJets || i>=j) {
164 std::cout << "ERROR in KtJetTable::DijTable::operator()" << std::endl;
165 std::cout << " Attempt to access element (" << i << "," << j << ") in table with nJets, nRows = "
166 << m_nJets << ", " << m_nRows << std::endl;
167 }
168 return *(m_table.begin() + i*m_nRows + j);
169}
170
171
172KtFloat KtJetTable::DijTable::operator() (int i, int j) const {
173 if (i<0 || j<0 || i>=m_nJets || j>=m_nJets) {
174 std::cout << "ERROR in KtJetTable::DijTable::operator() const";
175 std::cout << " Attempt to access element (" << i << "," << j << ") in table with nJets, nRows = "
176 << m_nJets << ", " << m_nRows << std::endl;
177 }
178 return *(m_table.begin() + i*m_nRows + j);
179}
180
181void KtJetTable::DijTable::killJet() {
182 if (m_nJets<=0) {
183 std::cout << "ERROR in KtJetTable::DijTable::killJet()" << std::endl;
184 std::cout << " Called when m_nJets = " << m_nJets << std::endl;
185 }
186 --m_nJets;
187}
188}//end of namespace
Note: See TracBrowser for help on using the repository browser.