15 #include <ecl/modules/eclCovarianceMatrix/ECLCovarianceMatrixModule.h>
21 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
24 #include <ecl/dataobjects/ECLShower.h>
39 m_eclShowers(eclShowerArrayName()),
40 m_eventLevelClusteringInfo(eventLevelClusteringInfoName())
43 setDescription(
"ECLCovarianceMatrix: Sets the ECL photon shower covariance matrix.");
44 setPropertyFlags(c_ParallelProcessingCertified);
70 double background = 0.0;
81 const double energy = eclShower.getEnergy();
85 if (energy > 0.) invEnergy = 1. / energy;
88 double invRoot2Energy = 0;
89 if (energy > 0.) invRoot2Energy = 1. / TMath::Power(energy, 0.5);
92 double invRoot4Energy = 0;
93 if (energy > 0.) invRoot4Energy = 1. / TMath::Power(energy, 0.25);
95 int detregion = eclShower.getDetectorRegion();
96 if (detregion == 11 or detregion == 13) detregion = 2;
98 double sigmaEnergy = 0.;
99 double sigmaTheta = 0.;
100 double sigmaPhi = 0.;
103 if (detregion == 1 and background < 0.05) {
104 if (energy <= 0.022) {sigmaEnergy = energy * (0.0449188); }
105 else {sigmaEnergy = energy * (-0.0912379 * invEnergy + 1.91849 * invRoot2Energy + -2.82169 * invRoot4Energy + 3.03119) / 100.;}
106 }
else if (detregion == 1 and background <= 0.3) {
107 if (energy <= 0.022) {sigmaEnergy = energy * (0.100769); }
108 else {sigmaEnergy = energy * (-0.0725928 * invEnergy + 3.4685 * invRoot2Energy + -5.44127 * invRoot4Energy + 4.12045) / 100.;}
109 }
else if (detregion == 1 and background > 0.3) {
110 if (energy <= 0.022) {sigmaEnergy = energy * (0.148462); }
111 else {sigmaEnergy = energy * (-0.239931 * invEnergy + 6.94958 * invRoot2Energy + -10.4085 * invRoot4Energy + 5.92412) / 100.;}
114 if (detregion == 1 and background < 0.05) {
115 if (energy <= 0.0289609) { sigmaTheta = 1e-3 * (6.4008);}
116 else {sigmaTheta = 1e-3 * (-0.110397 * invEnergy + 0.753603 * invRoot2Energy + 2.63652 * invRoot4Energy + -0.606703);}
117 }
else if (detregion == 1 and background <= 0.3) {
118 if (energy <= 0.0274562) { sigmaTheta = 1e-3 * (7.39868);}
119 else {sigmaTheta = 1e-3 * (-0.207278 * invEnergy + 2.68616 * invRoot2Energy + -0.905487 * invRoot4Energy + 0.961485);}
120 }
else if (detregion == 1 and background > 0.3) {
121 if (energy <= 0.022) { sigmaTheta = 1e-3 * (8.83505);}
122 else {sigmaTheta = 1e-3 * (-0.160921 * invEnergy + 2.35311 * invRoot2Energy + -0.107975 * invRoot4Energy + 0.565367);}
125 if (detregion == 1 and background < 0.05) {
126 if (energy <= 0.0391279) {sigmaPhi = 1e-3 * (15.4147);}
127 else {sigmaPhi = 1e-3 * (-0.35934 * invEnergy + 2.01807 * invRoot2Energy + 7.26313 * invRoot4Energy + -1.93438);}
128 }
else if (detregion == 1 and background <= 0.3) {
129 if (energy <= 0.022) {sigmaPhi = 1e-3 * (18.1096);}
130 else {sigmaPhi = 1e-3 * (-0.178258 * invEnergy + 1.39155 * invRoot2Energy + 6.97462 * invRoot4Energy + -1.2795);}
131 }
else if (detregion == 1 and background > 0.3) {
132 if (energy <= 0.0224565) {sigmaPhi = 1e-3 * (20.1605);}
133 else {sigmaPhi = 1e-3 * (-0.344617 * invEnergy + 3.80823 * invRoot2Energy + 4.08729 * invRoot4Energy + -0.464714);}
136 if (detregion == 2 and background < 0.05) {
137 if (energy <= 0.022) {sigmaEnergy = energy * (0.0466552); }
138 else {sigmaEnergy = energy * (-0.0473491 * invEnergy + 1.02015 * invRoot2Energy + -0.705164 * invRoot4Energy + 1.77086) / 100.;}
139 }
else if (detregion == 2 and background <= 0.3) {
140 if (energy <= 0.022) {sigmaEnergy = energy * (0.112299); }
141 else {sigmaEnergy = energy * (0.0469588 * invEnergy + 1.81959 * invRoot2Energy + -2.13666 * invRoot4Energy + 2.37568) / 100.;}
142 }
else if (detregion == 2 and background > 0.3) {
143 if (energy <= 0.022) {sigmaEnergy = energy * (0.164881); }
144 else {sigmaEnergy = energy * (-0.267069 * invEnergy + 7.6176 * invRoot2Energy + -11.0341 * invRoot4Energy + 5.9203) / 100.;}
147 if (detregion == 2 and background < 0.05) {
148 if (energy <= 0.0309077) { sigmaTheta = 1e-3 * (9.31098);}
149 else {sigmaTheta = 1e-3 * (-0.221073 * invEnergy + 2.19876 * invRoot2Energy + 1.50833 * invRoot4Energy + 0.359609);}
150 }
else if (detregion == 2 and background <= 0.3) {
151 if (energy <= 0.022) { sigmaTheta = 1e-3 * (11.6166);}
152 else {sigmaTheta = 1e-3 * (-0.126603 * invEnergy + 1.81898 * invRoot2Energy + 1.85132 * invRoot4Energy + 0.300728);}
153 }
else if (detregion == 2 and background > 0.3) {
154 if (energy <= 0.022) { sigmaTheta = 1e-3 * (13.159);}
155 else {sigmaTheta = 1e-3 * (-0.316561 * invEnergy + 4.94686 * invRoot2Energy + -3.12199 * invRoot4Energy + 2.30275);}
158 if (detregion == 2 and background < 0.05) {
159 if (energy <= 0.0381178) {sigmaPhi = 1e-3 * (10.7661);}
160 else {sigmaPhi = 1e-3 * (-0.21842 * invEnergy + 0.648976 * invRoot2Energy + 7.1901 * invRoot4Energy + -3.10025);}
161 }
else if (detregion == 2 and background <= 0.3) {
162 if (energy <= 0.022237) {sigmaPhi = 1e-3 * (12.6804);}
163 else {sigmaPhi = 1e-3 * (-0.141507 * invEnergy + 0.488293 * invRoot2Energy + 7.30055 * invRoot4Energy + -3.13591);}
164 }
else if (detregion == 2 and background > 0.3) {
165 if (energy <= 0.024905) {sigmaPhi = 1e-3 * (14.8726);}
166 else {sigmaPhi = 1e-3 * (-0.323752 * invEnergy + 3.60933 * invRoot2Energy + 2.48526 * invRoot4Energy + -1.25493);}
169 if (detregion == 3 and background < 0.05) {
170 if (energy <= 0.022) {sigmaEnergy = energy * (0.0466662); }
171 else {sigmaEnergy = energy * (-0.106729 * invEnergy + 2.38163 * invRoot2Energy + -4.11299 * invRoot4Energy + 4.14056) / 100.;}
172 }
else if (detregion == 3 and background <= 0.3) {
173 if (energy <= 0.022) {sigmaEnergy = energy * (0.128851); }
174 else {sigmaEnergy = energy * (-0.204362 * invEnergy + 6.43 * invRoot2Energy + -10.8673 * invRoot4Energy + 7.04059) / 100.;}
175 }
else if (detregion == 3 and background > 0.3) {
176 if (energy <= 0.022) {sigmaEnergy = energy * (0.255153); }
177 else {sigmaEnergy = energy * (-0.316245 * invEnergy + 11.5211 * invRoot2Energy + -18.3626 * invRoot4Energy + 9.89422) / 100.;}
180 if (detregion == 3 and background < 0.05) {
181 if (energy <= 0.0318079) { sigmaTheta = 1e-3 * (10.4446);}
182 else {sigmaTheta = 1e-3 * (-0.295827 * invEnergy + 3.37748 * invRoot2Energy + -0.284446 * invRoot4Energy + 1.48098);}
183 }
else if (detregion == 3 and background <= 0.3) {
184 if (energy <= 0.022) { sigmaTheta = 1e-3 * (13.828);}
185 else {sigmaTheta = 1e-3 * (-0.252479 * invEnergy + 4.18833 * invRoot2Energy + -2.04822 * invRoot4Energy + 2.38486);}
186 }
else if (detregion == 3 and background > 0.3) {
187 if (energy <= 0.022) { sigmaTheta = 1e-3 * (17.0354);}
188 else {sigmaTheta = 1e-3 * (-0.420215 * invEnergy + 7.08904 * invRoot2Energy + -5.88269 * invRoot4Energy + 3.61643);}
191 if (detregion == 3 and background < 0.05) {
192 if (energy <= 0.022) {sigmaPhi = 1e-3 * (15.129);}
193 else {sigmaPhi = 1e-3 * (-0.0534441 * invEnergy + -1.72466 * invRoot2Energy + 12.8625 * invRoot4Energy + -4.21203);}
194 }
else if (detregion == 3 and background <= 0.3) {
195 if (energy <= 0.022) {sigmaPhi = 1e-3 * (20.1408);}
196 else {sigmaPhi = 1e-3 * (0.0869293 * invEnergy + -1.93352 * invRoot2Energy + 12.8105 * invRoot4Energy + -4.03756);}
197 }
else if (detregion == 3 and background > 0.3) {
198 if (energy <= 0.022) {sigmaPhi = 1e-3 * (23.2771);}
199 else {sigmaPhi = 1e-3 * (-0.439611 * invEnergy + 5.74196 * invRoot2Energy + 1.76693 * invRoot4Energy + -0.0407866);}
202 B2DEBUG(175,
"energy=" << energy <<
", detector region=" << detregion);
203 B2DEBUG(175,
"sigmaEnergy=" << sigmaEnergy <<
", sigmaPhi=" << sigmaPhi <<
", sigmaTheta=" << sigmaTheta);
205 double covMatrix[6] = {sigmaEnergy * sigmaEnergy, 0.0, sigmaPhi * sigmaPhi, 0.0, 0.0, sigmaTheta * sigmaTheta};
206 eclShower.setCovarianceMatrix(covMatrix);