|
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
|
// Implementation of class EcalShowerContainmentCorrections
// Author: Stefano Argiro'
// $Id: EcalShowerContainmentCorrections.cc,v 1.1 2007/05/15 20:37:22 argiro Exp $
#include "CondFormats/EcalCorrections/interface/EcalShowerContainmentCorrections.h"
#include <DataFormats/EcalDetId/interface/EBDetId.h>
#include <FWCore/MessageLogger/interface/MessageLogger.h>
//#include <iostream>
#include "CondCore/CondDB/interface/Serialization.h"
#include "CondFormats/External/interface/EcalDetID.h"
const EcalShowerContainmentCorrections::Coefficients EcalShowerContainmentCorrections::correctionCoefficients(
const EBDetId& centerxtal) const {
GroupMap::const_iterator iter = groupmap_.find(centerxtal.rawId());
if (iter != groupmap_.end()) {
int group = iter->second;
return coefficients_[group - 1];
}
edm::LogError("ShowerContaiment Correction not found");
return Coefficients();
}
void EcalShowerContainmentCorrections::fillCorrectionCoefficients(const EBDetId& xtal,
int group,
const Coefficients& coefficients) {
// do not replace if we already have the xtal
if (groupmap_.find(xtal) != groupmap_.end())
return;
groupmap_[xtal] = group;
if (coefficients_.size() < (unsigned int)(group)) {
coefficients_.resize(group);
coefficients_[group - 1] = coefficients;
}
// we don't need to fill coefficients if the group has already been inserted
}
void EcalShowerContainmentCorrections::fillCorrectionCoefficients(const int supermodule,
const int module,
const Coefficients& coefficients) {
if (module > EBDetId::kModulesPerSM) {
edm::LogError("Invalid Module Number");
return;
}
// what is EBDetID::kModuleBoundaries ? we better redefine them here ...
const int kModuleLow[] = {1, 501, 901, 1301};
const int kModuleHigh[] = {500, 900, 1300, 1700};
for (int xtal = kModuleLow[module - 1]; xtal <= kModuleHigh[module - 1]; ++xtal) {
EBDetId detid(supermodule, xtal, EBDetId::SMCRYSTALMODE);
fillCorrectionCoefficients(detid, module, coefficients);
}
}
/** Calculate the correction for the given direction and type */
const double EcalShowerContainmentCorrections::correctionXY(const EBDetId& xtal,
double position,
EcalShowerContainmentCorrections::Direction dir,
EcalShowerContainmentCorrections::Type type) const {
GroupMap::const_iterator iter = groupmap_.find(xtal);
if (iter == groupmap_.end())
return -1;
int group = iter->second;
EcalShowerContainmentCorrections::Coefficients coeff = coefficients_[group - 1];
int offset = 0;
if (dir == eY)
offset += 2 * Coefficients::kPolynomialDegree;
if (position < 0)
offset += Coefficients::kPolynomialDegree;
if (type == e5x5)
offset += 4 * Coefficients::kPolynomialDegree;
double corr = 0;
for (int i = offset; i < offset + Coefficients::kPolynomialDegree; ++i) {
corr += coeff.data[i] * pow(position, i - offset);
}
return corr;
}
const double EcalShowerContainmentCorrections::correction3x3(const EBDetId& xtal, const math::XYZPoint& pos) const {
double x = pos.X() * 10; // correction functions use mm
double y = pos.Y() * 10;
double corrx = correctionXY(xtal, x, eX, e3x3);
double corry = correctionXY(xtal, y, eY, e3x3);
return corrx * corry;
}
const double EcalShowerContainmentCorrections::correction5x5(const EBDetId& xtal, const math::XYZPoint& pos) const {
double x = pos.X() * 10; // correction functions use mm
double y = pos.Y() * 10;
double corrx = correctionXY(xtal, x, eX, e5x5);
double corry = correctionXY(xtal, y, eY, e5x5);
return corrx * corry;
}
|