source: trunk/KtJet/KtEvent.cc

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

first commit
File size: 13.8 KB
RevLine 
[2]1#include "KtJet/KtEvent.h"
2#include "KtJet/KtLorentzVector.h"
3#include "KtJet/KtJetTable.h"
4#include "KtJet/KtDistance.h"
5#include "KtJet/KtDistanceInterface.h"
6#include "KtJet/KtRecom.h"
7#include "KtJet/KtRecomInterface.h"
8#include <vector>
9#include <utility>
10#include <iostream>
11#include <string>
12#include <algorithm>
13#include <cmath>
14
15namespace KtJet {
16 //using CLHEP::HepLorentzVector;
17 using namespace CLHEP;
18
19 /** Constructor for inclusive method */
20 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
21 int type, int angle, int recom,
22 KtFloat rParameter)
23 : m_type(type),
24 m_rParameterSq(rParameter*rParameter), m_inclusive(true),
25 m_jets() {
26 init(constituents,0,angle,0,recom);
27 }
28
29 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
30 int type, KtDistance *angle, KtRecom *recom,
31 KtFloat rParameter)
32 : m_type(type),
33 m_rParameterSq(rParameter*rParameter), m_inclusive(true),
34 m_jets() {
35 init(constituents,angle,0,recom,0);
36 }
37
38 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
39 int type, KtDistance *angle, int recom,
40 KtFloat rParameter)
41 : m_type(type),
42 m_rParameterSq(rParameter*rParameter), m_inclusive(true),
43 m_jets() {
44 init(constituents,angle,0,0,recom);
45 }
46
47 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
48 int type, int angle, KtRecom *recom,
49 KtFloat rParameter)
50 : m_type(type),
51 m_rParameterSq(rParameter*rParameter), m_inclusive(true),
52 m_jets() {
53 init(constituents,0,angle,recom,0);
54 }
55
56 KtEvent::KtEvent(const std::vector<HepLorentzVector> & constituents,
57 int type, int angle, int recom,
58 KtFloat rParameter)
59 : m_type(type),
60 m_rParameterSq(rParameter*rParameter), m_inclusive(true),
61 m_jets() {
62 std::vector<KtLorentzVector> konstituents;
63 makeKtFromHepLV(konstituents,constituents);
64 init(konstituents,0,angle,0,recom);
65 }
66
67 /** Constructor for exclusive method */
68 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
69 int type, int angle, int recom)
70 : m_type(type),
71 m_rParameterSq(1), m_inclusive(false),
72 m_jets() {
73 init(constituents,0,angle,0,recom);
74 }
75
76 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
77 int type, KtDistance *angle, KtRecom *recom)
78 : m_type(type),
79 m_rParameterSq(1), m_inclusive(false),
80 m_jets() {
81 init(constituents,angle,0,recom,0);
82 }
83
84 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
85 int type, KtDistance *angle, int recom)
86 : m_type(type),
87 m_rParameterSq(1), m_inclusive(false),
88 m_jets() {
89 init(constituents,angle,0,0,recom);
90 }
91
92 KtEvent::KtEvent(const std::vector<KtLorentzVector> & constituents,
93 int type, int angle, KtRecom *recom)
94 : m_type(type),
95 m_rParameterSq(1), m_inclusive(false),
96 m_jets() {
97 init(constituents,0,angle,recom,0);
98 }
99
100 KtEvent::KtEvent(const std::vector<HepLorentzVector> & constituents,
101 int type, int angle, int recom)
102 : m_type(type),
103 m_rParameterSq(1), m_inclusive(false),
104 m_jets() {
105 std::vector<KtLorentzVector> konstituents;
106 makeKtFromHepLV(konstituents,constituents);
107 init(konstituents,0,angle,0,recom);
108 }
109
110 /** Constructor for (exclusive) subjet method */
111 KtEvent::KtEvent(const KtLorentzVector & jet, int angle, int recom)
112 : m_type(1),
113 m_rParameterSq(1), m_inclusive(false),
114 m_jets() {
115 init(jet,0,angle,0,recom);
116 }
117
118 KtEvent::KtEvent(const KtLorentzVector & jet, KtDistance *angle, KtRecom *recom)
119 : m_type(1),
120 m_rParameterSq(1), m_inclusive(false),
121 m_jets() {
122 init(jet,angle,0,recom,0);
123 }
124
125 KtEvent::KtEvent(const KtLorentzVector & jet, KtDistance *angle, int recom)
126 : m_type(1),
127 m_rParameterSq(1), m_inclusive(false),
128 m_jets() {
129 init(jet,angle,0,0,recom);
130 }
131
132 KtEvent::KtEvent(const KtLorentzVector & jet, int angle, KtRecom *recom)
133 : m_type(1),
134 m_rParameterSq(1), m_inclusive(false),
135 m_jets() {
136 init(jet,0,angle,recom,0);
137 }
138
139 KtEvent::~KtEvent() {
140 delete m_ktDist;
141 delete m_ktRecom;
142 }
143
144 void KtEvent::findJetsN(int nJets) {
145 // Uses merging history created by makeJets
146 if (m_inclusive) {
147 std::cout << "WARNING in KtEvent::findJetsN : Trying to find exclusive jets in inclusive event\n";
148 }
149 m_jets.clear();
150 KtJetTable jt(m_constituents,m_ktDist,m_ktRecom);
151 int nParticles = m_constituents.size();
152 int njet = m_constituents.size();
153 while (njet>nJets) {
154 int hist = m_hist[njet-1];
155 if (hist>=0) {
156 jt.killJet(hist);
157 } else {
158 // Merge two jets
159 hist = -hist;
160 int iPair = hist % nParticles;
161 int jPair = hist / nParticles;
162 jt.mergeJets(iPair,jPair);
163 }
164 --njet;
165 }
166 // Now jets remaining in jt should be required jets: copy to m_jets
167 for (int i=0; i<jt.getNJets(); ++i) { m_jets.push_back(jt.getJet(i)); }
168 }
169
170 /** Do exclusive jet-finding to scale dCut */
171 void KtEvent::findJetsD(KtFloat dCut) {
172 // Find number of jets at scale dCut, then call findJetsN to do merging
173 int njets = 1;
174 for (int iloop = m_constituents.size(); iloop>1; --iloop) {
175 if (m_dMerge[iloop-1]>=dCut) {
176 njets = iloop;
177 break;
178 }
179 }
180 findJetsN(njets);
181 }
182
183 /** Do exclusive jet-finding to scale yCut */
184 void KtEvent::findJetsY(KtFloat yCut) {
185 // Find jets at scale yCut by rescaling to give dCut, then calling findJetsD
186 KtFloat dCut = yCut * m_eCut * m_eCut;
187 findJetsD(dCut);
188 }
189
190 /** Get jets */
191 std::vector<KtLorentzVector> KtEvent::getJets() {
192 return m_jets;
193 }
194
195 std::vector<KtLorentzVector> KtEvent::getJetsE() {
196 std::vector<KtLorentzVector> a(m_jets);
197 std::sort(a.begin(),a.end(),greaterE);
198 return a;
199 }
200
201 std::vector<KtLorentzVector> KtEvent::getJetsEt() {
202 std::vector<KtLorentzVector> a(m_jets);
203 std::sort(a.begin(),a.end(),greaterEt);
204 return a;
205 }
206
207 std::vector<KtLorentzVector> KtEvent::getJetsPt() {
208 std::vector<KtLorentzVector> a(m_jets);
209 std::sort(a.begin(),a.end(),greaterPt);
210 return a;
211 }
212
213 std::vector<KtLorentzVector> KtEvent::getJetsRapidity() {
214 std::vector<KtLorentzVector> a(m_jets);
215 std::sort(a.begin(),a.end(),greaterRapidity);
216 return a;
217 }
218
219 std::vector<KtLorentzVector> KtEvent::getJetsEta() {
220 std::vector<KtLorentzVector> a(m_jets);
221 std::sort(a.begin(),a.end(),greaterEta);
222 return a;
223 }
224
225 std::vector<const KtLorentzVector *> KtEvent::getConstituents() const {
226 std::vector<const KtLorentzVector *> a;
227 std::vector<KtLorentzVector>::const_iterator itr = m_constituents.begin();
228 for (; itr != m_constituents.end(); ++itr) {
229 a.push_back(&*itr); // ??? Converted to pointer?
230 }
231 return a;
232 }
233
234 std::vector<KtLorentzVector> KtEvent::copyConstituents() const {
235 std::vector<KtLorentzVector> a(m_constituents);
236 return a;
237 }
238
239 KtLorentzVector KtEvent::getJet(const KtLorentzVector & a) const {
240 std::vector<KtLorentzVector>::const_iterator itr = m_jets.begin();
241 for (; itr != m_jets.end(); ++itr) {
242 if (a == *itr || itr->contains(a)) return *itr;
243 }
244 std::cout << "ERROR in KtEvent::getJet : particle not in event" << std::endl;
245 return a; // ??? Do something more sensible here?
246 }
247
248
249 /*********************
250 * Private methods *
251 *********************/
252
253 /** Initialization for event jet analysis */
254 void KtEvent::init(const std::vector<KtLorentzVector> & constituents,
255 KtDistance* dist, int idist, KtRecom* recom, int irecom) {
256 static bool first_inclusive = true, first_exclusive = true;
257 listAuthors();
258 setSchemes(dist,idist,recom,irecom);
259 m_constituents.clear();
260 std::vector<KtLorentzVector>::const_iterator itr;
261 for (itr=constituents.begin(); itr!=constituents.end(); ++itr) {
262 KtLorentzVector v = (*m_ktRecom)(*itr);
263 m_constituents.push_back(v);
264 }
265 addEnergy();
266 setECut(getETot());
267 if (first_inclusive && m_inclusive) {
268 first_inclusive = false;
269 printSteering("inclusive");
270 }
271 if (first_exclusive && !m_inclusive) {
272 first_exclusive = false;
273 printSteering("exclusive");
274 }
275 makeJets();
276 if (!isInclusive()) m_jets.clear(); // Only need merging history in exclusive case.
277 // Jets reconstructed later with fixed d cut or njets.
278 }
279
280 /** Initialization for subjet analysis */
281 void KtEvent::init(const KtLorentzVector & jet, KtDistance* dist, int idist, KtRecom* recom, int irecom) {
282 static bool first_subjet = true;
283 listAuthors();
284 setSchemes(dist,idist,recom,irecom);
285 std::vector<const KtLorentzVector*> constituents = jet.getConstituents();
286 m_constituents.clear();
287 std::vector<const KtLorentzVector *>::const_iterator itr;
288 for (itr=constituents.begin(); itr!=constituents.end(); ++itr) {
289 KtLorentzVector v = (*m_ktRecom)(**itr);
290 m_constituents.push_back(v);
291 }
292 addEnergy();
293 setECut(jet.perp());
294 if (first_subjet) {
295 first_subjet = false;
296 printSteering("subjet");
297 }
298 makeJets();
299 m_jets.clear(); // Only need merging history. Jets reconstructed later with fixed d cut or njets.
300 }
301
302 void KtEvent::addEnergy() {
303 m_eTot = 0;
304 std::vector<KtLorentzVector>::const_iterator itr;
305 for (itr = m_constituents.begin(); itr != m_constituents.end(); ++itr) {
306 m_eTot += itr->e(); // Add up energy in event
307 }
308 }
309
310 void KtEvent::makeJets() {
311 int nParticles = m_constituents.size();
312 if (nParticles==0) return; // Do nothing if no input particles
313
314 m_dMerge.resize(nParticles); // Reserve space for D-cut vector
315 m_hist.resize(nParticles); // Reserve space for merging history vector
316 KtJetTable jt(m_constituents,m_ktDist,m_ktRecom);
317
318 int njet;
319 // KtFloat dMax = 0;
320 while ((njet=jt.getNJets())>1) { // Keep merging until only one jet left
321 int iPairMin, jPairMin, iJetMin;
322 KtFloat dPairMin, dJetMin, dMin;
323 std::pair<int,int> pPairMin = jt.getMinDPair(); // Find jet pair with minimum D
324 iPairMin = pPairMin.first;
325 jPairMin = pPairMin.second;
326 dPairMin = jt.getD(iPairMin,jPairMin);
327 iJetMin = jt.getMinDJet(); // Find jet with minimum D to beam
328 dJetMin = jt.getD(iJetMin);
329 dJetMin *= m_rParameterSq; // Scale D to beam by rParameter
330 bool mergeWithBeam = ((m_type != 1 && (dJetMin<=dPairMin)) || njet==1); // Merge with beam remnant?
331 dMin = mergeWithBeam ? dJetMin : dPairMin;
332
333 if (mergeWithBeam) { // Merge jet with beam remnant
334 m_jets.push_back(jt.getJet(iJetMin)); // Add jet to inclusive list
335 m_hist[njet-1] = iJetMin; // Record which jet killed (+ve to indicate merge with beam)
336 jt.killJet(iJetMin); // Chuck out jet
337 } else { // Merge pair of jets
338 m_hist[njet-1] = -(jPairMin*nParticles+iPairMin); // Record which jet pair merged (-ve to indicate pair)
339 jt.mergeJets(iPairMin,jPairMin); // Merge jets
340 }
341 m_dMerge[njet-1] = dMin; // Store D where jets merged
342 }
343 // End of loop: now should have njet = 1
344 m_jets.push_back(jt.getJet(0)); // Add last jet to list
345 m_dMerge[0] = jt.getD(0); // Put last jet kt in vectors ???
346 }
347
348 void KtEvent::makeKtFromHepLV(std::vector<KtLorentzVector> & kt,
349 const std::vector<HepLorentzVector> & hep) {
350 std::vector<HepLorentzVector>::const_iterator itr = hep.begin();
351 for (; itr != hep.end(); ++itr) {
352 KtLorentzVector ktvec(*itr);
353 kt.push_back(ktvec);
354 }
355 }
356
357 void KtEvent::setSchemes(KtDistance* dist, int idist, KtRecom* recom, int irecom) {
358 m_angle = idist;
359 m_ktDist = (dist) ? dist : getDistanceScheme(idist,m_type);
360 m_recom = irecom;
361 m_ktRecom = (recom) ? recom : getRecomScheme(irecom);
362 }
363
364 void KtEvent::printSteering(std::string mode) const {
365 static std::string collisionType[4] = {"e+e-","ep","pe","pp"};
366 std::string type = (m_type>=1 && m_type <= 4) ? collisionType[m_type-1] : "unknown";
367 std::string angle = m_ktDist->name();
368 std::string recom = m_ktRecom->name();
369 mode.resize(10,' ');
370 type.resize(16,' ');
371 angle.resize(16,' ');
372 recom.resize(16,' ');
373 std::cout << "******************************************\n";
374 std::cout << "* KtEvent constructor called: " << mode << " *\n";
375 std::cout << "* Collision type: " << type << " *\n";
376 std::cout << "* Kt scheme: " << angle << " *\n";
377 std::cout << "* Recombination scheme: " << recom << " *\n";
378#ifdef KTDOUBLEPRECISION
379 std::cout << "* Compiled to use double precision. *\n";
380#else
381 std::cout << "* Compiled to use single precision. *\n";
382#endif
383 std::cout << "******************************************\n";
384 }
385
386 void KtEvent::listAuthors() const {
387 static bool first = true;
388 if (first) {
389 first = false;
390 std::cout << "***********************************************\n";
391 std::cout << "* Package KtJet written by: *\n";
392 std::cout << "* Jon Butterworth *\n";
393 std::cout << "* Jon Couchman *\n";
394 std::cout << "* Brian Cox *\n";
395 std::cout << "* Ben Waugh *\n";
396 std::cout << "* See documentation at <http://www.ktjet.org> *\n";
397 std::cout << "***********************************************\n";
398 }
399 }
400
401}//end of namespace
Note: See TracBrowser for help on using the repository browser.