|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
#include "DataFormats/L1THGCal/interface/ClusterShapes.h"
#include <cmath>
using namespace l1t;
ClusterShapes ClusterShapes::operator+(const ClusterShapes& x) {
ClusterShapes cs(*this); // copy constructor
cs += x;
return cs;
}
void ClusterShapes::operator+=(const ClusterShapes& x) {
sum_e_ += x.sum_e_;
sum_e2_ += x.sum_e2_;
sum_logE_ += x.sum_logE_;
n_ += x.n_;
sum_w_ += x.sum_w_;
emax_ = (emax_ > x.emax_) ? emax_ : x.emax_;
// mid-point
sum_eta_ += x.sum_eta_;
sum_phi_0_ += x.sum_phi_0_; //
sum_phi_1_ += x.sum_phi_1_; //
sum_r_ += x.sum_r_;
// square
sum_eta2_ += x.sum_eta2_;
sum_phi2_0_ += x.sum_phi2_0_;
sum_phi2_1_ += x.sum_phi2_1_;
sum_r2_ += x.sum_r2_;
// off diagonal
sum_eta_r_ += x.sum_eta_r_;
sum_r_phi_0_ += x.sum_r_phi_0_;
sum_r_phi_1_ += x.sum_r_phi_1_;
sum_eta_phi_0_ += x.sum_eta_phi_0_;
sum_eta_phi_1_ += x.sum_eta_phi_1_;
}
// -------------- CLUSTER SHAPES ---------------
void ClusterShapes::Init(float e, float eta, float phi, float r) {
if (e <= 0)
return;
sum_e_ = e;
sum_e2_ = e * e;
sum_logE_ = std::log(e);
float w = e;
n_ = 1;
sum_w_ = w;
sum_phi_0_ = w * (phi);
sum_phi_1_ = w * (phi + M_PI);
sum_r_ = w * r;
sum_eta_ = w * eta;
//--
sum_r2_ += w * (r * r);
sum_eta2_ += w * (eta * eta);
sum_phi2_0_ += w * (phi * phi);
sum_phi2_1_ += w * (phi + M_PI) * (phi + M_PI);
// -- off diagonal
sum_eta_r_ += w * (r * eta);
sum_r_phi_0_ += w * (r * phi);
sum_r_phi_1_ += w * r * (phi + M_PI);
sum_eta_phi_0_ += w * (eta * phi);
sum_eta_phi_1_ += w * eta * (phi + M_PI);
}
// ------
float ClusterShapes::Eta() const { return sum_eta_ / sum_w_; }
float ClusterShapes::R() const { return sum_r_ / sum_w_; }
float ClusterShapes::SigmaEtaEta() const { return sum_eta2_ / sum_w_ - Eta() * Eta(); }
float ClusterShapes::SigmaRR() const { return sum_r2_ / sum_w_ - R() * R(); }
float ClusterShapes::SigmaPhiPhi() const {
float phi_0 = (sum_phi_0_ / sum_w_);
float phi_1 = (sum_phi_1_ / sum_w_);
float spp_0 = sum_phi2_0_ / sum_w_ - phi_0 * phi_0;
float spp_1 = sum_phi2_1_ / sum_w_ - phi_1 * phi_1;
if (spp_0 < spp_1) {
isPhi0_ = true;
return spp_0;
} else {
isPhi0_ = false;
return spp_1;
}
}
float ClusterShapes::Phi() const {
SigmaPhiPhi(); //update phi
if (isPhi0_)
return (sum_phi_0_ / sum_w_);
else
return (sum_phi_1_ / sum_w_);
}
// off - diagonal
float ClusterShapes::SigmaEtaR() const { return -(sum_eta_r_ / sum_w_ - Eta() * R()); }
float ClusterShapes::SigmaEtaPhi() const {
SigmaPhiPhi(); // decide which phi use, update phi
if (isPhi0_)
return -(sum_eta_phi_0_ / sum_w_ - Eta() * (sum_phi_0_ / sum_w_));
else
return -(sum_eta_phi_1_ / sum_w_ - Eta() * (sum_phi_1_ / sum_w_));
}
float ClusterShapes::SigmaRPhi() const {
SigmaPhiPhi(); // decide which phi use, update phi
if (isPhi0_)
return -(sum_r_phi_0_ / sum_w_ - R() * (sum_phi_0_ / sum_w_));
else
return -(sum_r_phi_1_ / sum_w_ - R() * (sum_phi_1_ / sum_w_));
}
|