31 mDigitIndices.resize(geo.getNrows(), std::vector<int>(geo.getNcols(), -1));
32 mApplyCorrectionZ = applyCorrectionZ;
33 mApplyCorrectionE = applyCorrectionE;
35 mCrystalEnergyCorrectionPars.reserve(6);
36 mCrystalEnergyCorrectionPars[0] = 0.00444;
37 mCrystalEnergyCorrectionPars[1] = -1.322;
38 mCrystalEnergyCorrectionPars[2] = 1.021;
39 mCrystalEnergyCorrectionPars[3] = 0.0018;
40 mCrystalEnergyCorrectionPars[4] = 0.;
41 mCrystalEnergyCorrectionPars[5] = 0.;
43 mSamplingEnergyCorrectionPars.reserve(6);
44 mSamplingEnergyCorrectionPars[0] = 0.0033;
45 mSamplingEnergyCorrectionPars[1] = -2.09;
46 mSamplingEnergyCorrectionPars[2] = 1.007;
47 mSamplingEnergyCorrectionPars[3] = 0.0667;
48 mSamplingEnergyCorrectionPars[4] = -0.108;
49 mSamplingEnergyCorrectionPars[5] = 0.0566;
51 mCrystalZCorrectionPars.reserve(9);
52 mCrystalZCorrectionPars[0] = -0.005187;
53 mCrystalZCorrectionPars[1] = 0.7301;
54 mCrystalZCorrectionPars[2] = -0.7382;
55 mCrystalZCorrectionPars[3] = 0.;
56 mCrystalZCorrectionPars[4] = 0.;
57 mCrystalZCorrectionPars[5] = 0.;
58 mCrystalZCorrectionPars[6] = 0.;
59 mCrystalZCorrectionPars[7] = 0.;
60 mCrystalZCorrectionPars[8] = 0.;
62 mSamplingZCorrectionPars.reserve(9);
63 mSamplingZCorrectionPars[0] = -2.137;
64 mSamplingZCorrectionPars[1] = 6.400;
65 mSamplingZCorrectionPars[2] = -3.342;
66 mSamplingZCorrectionPars[3] = -0.1364;
67 mSamplingZCorrectionPars[4] = 0.4019;
68 mSamplingZCorrectionPars[5] = -0.1969;
69 mSamplingZCorrectionPars[6] = 0.008223;
70 mSamplingZCorrectionPars[7] = -0.02425;
71 mSamplingZCorrectionPars[8] = 0.01190;
73 fCrystalShowerShape =
new TF1(
"fCrystal",
"x<[1] ? [0]*exp([3]*x+[4]*x*x+[5]*x*x*x) : (x<[2] ? [0]*[6]*exp([7]*x+[8]*x*x) : [0]*[9]*exp([10]*x+[11]*x*x))", 0, 15);
88 fCrystalShowerShape->SetParameters(pc);
90 fSamplingShowerShape =
new TF1(
"fSampling",
"x<[1] ? [0]*exp([3]*x+[4]*x*x+[5]*x*x*x) : (x<[2] ? [0]*[6]*exp([7]*x+[8]*x*x) : [0]*[9]*exp([10]*x+[11]*x*x))", 0, 15);
104 fSamplingShowerShape->SetParameters(ps);
106 fCrystalRMS =
new TF1(
"fCrystalRMS",
"[0]*x*exp([1]*x+[2]*x*x+[3]*x*x*x)", 0, 2.2);
112 fCrystalRMS->SetParameters(p);
136 if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols()) {
139 int digitIndex = mDigitIndices[
row][
col];
140 LOGP(
debug,
" checking row={} and col={} digitIndex={}",
row,
col, digitIndex);
141 if (digitIndex < 0) {
148 auto [sector1, ch1] = geo.getSectorChamber(digit.
getTower());
149 auto [sector2, ch2] = geo.getSectorChamber(digit2.
getTower());
150 LOGP(
debug,
" checking if sector and chamber are the same: ({},{}) ({},{})", sector1, ch1, sector2, ch2);
151 if (sector1 != sector2 || ch1 != ch2) {
156 mDigitIndices[
row][
col] = -1;
158 LOGP(
debug,
" adding new digit at row={} and col={}",
row,
col);
250 std::vector<double>
x(mult);
251 std::vector<double>
y(mult);
252 std::vector<double>
z(mult);
253 std::vector<double> e(mult);
254 std::vector<std::vector<double>> eInClusters(mult, std::vector<double>(nMax));
259 for (
int idig = 0; idig < mult; idig++) {
265 std::vector<double> xMax(nMax);
266 std::vector<double>
yMax(nMax);
267 std::vector<double>
zMax(nMax);
268 std::vector<double> eMax(nMax);
270 for (
int iclu = 0; iclu < nMax; iclu++) {
271 xMax[iclu] =
x[digitId[iclu]];
272 yMax[iclu] =
y[digitId[iclu]];
273 zMax[iclu] =
z[digitId[iclu]];
274 eMax[iclu] = e[digitId[iclu]];
277 std::vector<double> prop(nMax);
281 bool insuficientAccuracy =
true;
283 while (insuficientAccuracy && nIterations < mNMaxIterations) {
286 for (
int idig = 0; idig < mult; idig++) {
287 double eEstimated = 0;
288 for (
int iclu = 0; iclu < nMax; iclu++) {
289 prop[iclu] = eMax[iclu] *
showerShape(std::sqrt((
x[idig] - xMax[iclu]) * (
x[idig] - xMax[iclu]) +
290 (
y[idig] -
yMax[iclu]) * (
y[idig] -
yMax[iclu])),
291 z[idig] -
zMax[iclu], isCrystal);
292 eEstimated += prop[iclu];
294 if (eEstimated == 0.) {
298 for (
int iclu = 0; iclu < nMax; iclu++) {
299 eInClusters[idig][iclu] = e[idig] * prop[iclu] / eEstimated;
303 insuficientAccuracy =
false;
304 for (
int iclu = 0; iclu < nMax; iclu++) {
305 double oldX = xMax[iclu];
306 double oldY =
yMax[iclu];
307 double oldZ =
zMax[iclu];
308 double oldE = eMax[iclu];
311 for (
int idig = 0; idig < mult; idig++) {
312 eMax[iclu] += eInClusters[idig][iclu];
318 for (
int idig = 0; idig < mult; idig++) {
319 double w = std::max(std::log(eInClusters[idig][iclu] / eMax[iclu]) + mLogWeight, 0.);
320 xMax[iclu] +=
x[idig] *
w;
321 yMax[iclu] +=
y[idig] *
w;
322 zMax[iclu] +=
z[idig] *
w;
331 insuficientAccuracy += (std::abs(eMax[iclu] - oldE) > mUnfogingEAccuracy);
332 insuficientAccuracy += (std::abs(xMax[iclu] - oldX) > mUnfogingXZAccuracy);
333 insuficientAccuracy += (std::abs(
yMax[iclu] - oldY) > mUnfogingXZAccuracy);
334 insuficientAccuracy += (std::abs(
zMax[iclu] - oldZ) > mUnfogingXZAccuracy);
340 for (
int iclu = 0; iclu < nMax; iclu++) {
341 auto&
clu = unfoldedClusters.emplace_back();
343 for (
int idig = 0; idig < mult; idig++) {
346 clu.
addDigit(jdigit, towerId, eInClusters[idig][iclu]);
360 double etot = cluster.getE();
361 for (
size_t i = 0;
i < cluster.getMultiplicity();
i++) {
362 float energy = cluster.getDigitEnergy(
i);
363 int towerId = cluster.getDigitTowerId(
i);
365 geo.detIdToGlobalPosition(towerId, xi, yi, zi);
366 double w = std::max(0., mLogWeight + std::log(energy / etot));
384 float ee = cluster.getE();
385 for (
size_t i = 0;
i < cluster.getMultiplicity();
i++) {
386 float energy = cluster.getDigitEnergy(
i);
387 int towerId = cluster.getDigitTowerId(
i);
389 geo.detIdToGlobalPosition(towerId, xi, yi, zi);
390 double r = std::sqrt((
x - xi) * (
x - xi) + (
y - yi) * (
y - yi) + (
z - zi) * (
z - zi));
394 double frac = fCrystalShowerShape->Eval(
r);
395 double rms = fCrystalRMS->Eval(
r);
396 chi2 += std::pow((energy / ee - frac) / rms, 2.);
399 cluster.setChi2(chi2 / ndf);
402 float eta = std::abs(cluster.getEta());
403 bool isCrystal = geo.isCrystal(cluster.getDigitTowerId(0));
404 if (mApplyCorrectionE) {
405 std::vector<double>& pe = isCrystal ? mCrystalEnergyCorrectionPars : mSamplingEnergyCorrectionPars;
406 ee *= pe[0] * std::pow(ee, pe[1]) + pe[2] + pe[3] * eta + pe[4] * eta * eta + pe[5] * eta * eta * eta;
409 if (mApplyCorrectionZ) {
410 std::vector<double>& pz = isCrystal ? mCrystalZCorrectionPars : mSamplingZCorrectionPars;
411 float zCor = (pz[0] + pz[1] * eta + pz[2] * eta * eta) + (pz[3] + pz[4] * eta + pz[5] * eta * eta) * ee + (pz[6] + pz[7] * eta + pz[8] * eta * eta) * ee * ee;
412 cluster.setZ(
z > 0 ?
z - zCor :
z + zCor);
417 for (
size_t i = 0;
i < cluster.getMultiplicity();
i++) {
418 int towerId = cluster.getDigitTowerId(
i);
419 if (geo.isAtTheEdge(towerId)) {
424 cluster.setEdgeFlag(isEdge);
426 LOGF(
debug,
"Cluster coordinates: (%6.2f,%6.2f,%6.2f)", cluster.getX(), cluster.getY(), cluster.getZ());