source: trunk/KtJet/KtDistance.cc@ 16

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

first commit

File size: 3.4 KB
Line 
1#include "KtJet/KtDistance.h"
2#include "KtJet/KtUtil.h"
3#include "KtJet/KtDistanceInterface.h"
4#include <cmath>
5#include <vector>
6#include <string>
7#include <iostream>
8namespace KtJet {
9
10
11KtDistance* getDistanceScheme(int angle, int collision_type) {
12 if (angle == 1) return new KtDistanceAngle(collision_type);
13 else if( angle == 2) return new KtDistanceDeltaR(collision_type);
14 else if (angle == 3) return new KtDistanceQCD(collision_type);
15 else{
16 std::cout << "WARNING, unreconised distance scheme specified!" << std::endl;
17 std::cout << "Distance Scheme set to KtDistanceAngle" << std::endl;
18 return new KtDistanceAngle(collision_type);
19 }
20}
21
22KtDistanceAngle::KtDistanceAngle(int collision_type) : m_type(collision_type), m_name("angle") {}
23 //KtDistanceAngle::~KtDistanceAngle() {}
24std::string KtDistanceAngle::name() const {return m_name;}
25
26KtFloat KtDistanceAngle::operator()(const KtLorentzVector & a) const {
27 KtFloat kt, r, costh;
28 const KtFloat small = 0.0001; // ??? Should be defined somewhere else?
29 switch (m_type) { // direction of beam depends on collision type
30 case 1:
31 return -1; // e+e- : no beam remnant, so result will be ignored anyway
32 break;
33 case 2: // ep (p beam -z direction)
34 costh = -(a.cosTheta());
35 break;
36 case 3: // pe (p beam +z direction)
37 costh = a.cosTheta();
38 break;
39 case 4: // pp (p beams in both directions)
40 costh = fabs(a.cosTheta());
41 break;
42 default: // type out of range - WARNING ???
43 costh = 0.;
44 break;
45 }
46 r = 2*(1-costh);
47 if (r<small) r = a.perp2()/a.vect().mag2(); // Use approx if close to beam
48 kt = a.e()*a.e() * r;
49 return kt;
50}
51
52KtFloat KtDistanceAngle::operator()(const KtLorentzVector & a, const KtLorentzVector & b) const {
53 KtFloat emin = std::min(a.e(),b.e());
54 KtFloat esq = emin*emin;
55 KtFloat costh = a.vect().cosTheta(b.vect());
56 return 2 * esq * (1 - costh);
57}
58
59
60KtDistanceDeltaR::KtDistanceDeltaR(int collision_type) : m_type(collision_type), m_name("DeltaR") {}
61 //KtDistanceDeltaR::~KtDistanceDeltaR() {}
62std::string KtDistanceDeltaR::name() const {return m_name;}
63
64KtFloat KtDistanceDeltaR::operator()(const KtLorentzVector & a) const {
65 return (m_type==1) ? -1 : a.perp2(); // If e+e-, no beam remnant, so result will be ignored anyway
66}
67
68KtFloat KtDistanceDeltaR::operator()(const KtLorentzVector & a, const KtLorentzVector & b) const {
69 KtFloat rsq,esq,kt,deltaEta,deltaPhi;
70 deltaEta = a.crapidity()-b.crapidity();
71 deltaPhi = phiAngle(a.phi()-b.phi());
72 rsq = deltaEta*deltaEta + deltaPhi*deltaPhi;
73 esq = std::min(a.perp2(),b.perp2());
74 kt = esq*rsq;
75 return kt;
76}
77
78
79KtDistanceQCD::KtDistanceQCD(int collision_type) : m_type(collision_type), m_name("QCD") {}
80 //KtDistanceQCD::~KtDistanceQCD() {}
81std::string KtDistanceQCD::name() const {return m_name;}
82
83KtFloat KtDistanceQCD::operator()(const KtLorentzVector & a) const {
84 return (m_type==1) ? -1 : a.perp2(); // If e+e-, no beam remnant, so result will be ignored anyway
85}
86
87KtFloat KtDistanceQCD::operator()(const KtLorentzVector & a, const KtLorentzVector & b) const {
88 KtFloat rsq,esq,kt,deltaEta,deltaPhi;
89 deltaEta = a.crapidity()-b.crapidity();
90 deltaPhi = phiAngle(a.phi()-b.phi());
91 rsq = 2 * (cosh(deltaEta)-cos(deltaPhi));
92 esq = std::min(a.perp2(),b.perp2());
93 kt = esq*rsq;
94 return kt;
95}
96
97}//end of namespace
Note: See TracBrowser for help on using the repository browser.