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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
|
#include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
#include "DataFormats/DetId/interface/DetId.h"
//Includes need to read from geometry
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
//by default is not initialized, gets initialized at first call
std::atomic<bool> EcalRingCalibrationTools::isInitializedFromGeometry_(false);
short EcalRingCalibrationTools::endcapRingIndex_[EEDetId::IX_MAX][EEDetId::IY_MAX];
std::once_flag EcalRingCalibrationTools::once_;
constexpr short EcalRingCalibrationTools::N_RING_TOTAL;
constexpr short EcalRingCalibrationTools::N_RING_BARREL;
constexpr short EcalRingCalibrationTools::N_RING_ENDCAP;
constexpr short EcalRingCalibrationTools::N_MODULES_BARREL;
short EcalRingCalibrationTools::getRingIndex(DetId id) {
if (id.det() != DetId::Ecal)
return -1;
if (id.subdetId() == EcalBarrel) {
if (EBDetId(id).ieta() < 0)
return EBDetId(id).ieta() + 85;
else
return EBDetId(id).ieta() + 84;
}
if (id.subdetId() == EcalEndcap) {
//needed only for the EE, it can be replaced at some point with something smarter
if (not isInitializedFromGeometry_)
throw std::logic_error(
"EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
EEDetId eid(id);
short endcapRingIndex = endcapRingIndex_[eid.ix() - 1][eid.iy() - 1] + N_RING_BARREL;
if (eid.zside() == 1)
endcapRingIndex += N_RING_ENDCAP / 2;
return endcapRingIndex;
}
return -1;
}
short EcalRingCalibrationTools::getModuleIndex(DetId id) {
if (id.det() != DetId::Ecal)
return -1;
if (id.subdetId() == EcalBarrel) {
short module = 4 * (EBDetId(id).ism() - 1) + EBDetId(id).im() - 1;
// std::cout<<"SM construction # : "<<EBDetId(id).ism() <<" SM installation # : "<< installationSMNumber[ EBDetId(id).ism() -1 ]<< "Xtal is Module :"<<module<< std::endl;
return module;
}
if (id.subdetId() == EcalEndcap) {
return -1;
}
return -1;
}
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInRing(short etaIndex) {
std::vector<DetId> ringIds;
if (etaIndex < 0)
return ringIds;
if (etaIndex < N_RING_BARREL) {
int k = 0;
if (etaIndex < 85)
k = -85 + etaIndex;
else
k = etaIndex - 84;
for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
if (EBDetId::validDetId(k, iphi))
ringIds.push_back(EBDetId(k, iphi));
}
else if (etaIndex < N_RING_TOTAL) {
//needed only for the EE, it can be replaced at some point maybe with something smarter
if (not isInitializedFromGeometry_)
throw std::logic_error(
"EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
int zside = (etaIndex < N_RING_BARREL + (N_RING_ENDCAP / 2)) ? -1 : 1;
short eeEtaIndex = (etaIndex - N_RING_BARREL) % (N_RING_ENDCAP / 2);
for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
if (endcapRingIndex_[ix][iy] == eeEtaIndex)
ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
}
return ringIds;
}
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInECAL() {
std::vector<DetId> ringIds;
for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta)
for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
if (EBDetId::validDetId(ieta, iphi))
ringIds.push_back(EBDetId(ieta, iphi));
//needed only for the EE, it can be replaced at some point maybe with something smarter
if (not isInitializedFromGeometry_)
throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
for (int zside = -1; zside < 2; zside += 2)
if (EEDetId::validDetId(ix + 1, iy + 1, zside))
ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
// std::cout<<" [EcalRingCalibrationTools::getDetIdsInECAL()] DetId.size() is "<<ringIds.size()<<std::endl;
return ringIds;
}
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInModule(short moduleIndex) {
std::vector<DetId> ringIds;
if (moduleIndex < 0)
return ringIds;
short moduleBound[5] = {1, 26, 46, 66, 86};
if (moduleIndex < EcalRingCalibrationTools::N_MODULES_BARREL) {
////////////////////////////nuovo ciclo
short sm, moduleInSm, zsm;
short minModuleiphi, maxModuleiphi, minModuleieta = 360, maxModuleieta = 0;
// if(moduleIndex%4 != 0 )
sm = moduleIndex / 4 + 1;
// else
//sm = moduleIndex/4;//i.e. module 8 belongs to sm=3, not sm=3
//if(moduleIndex%4 != 0 )
moduleInSm = moduleIndex % 4;
//else
//moduleInSm = 4;//moduleInSm is [1,2,3,4]
if (moduleIndex > 71)
zsm = -1;
else
zsm = 1;
minModuleiphi = ((sm - 1) % 18 + 1) * 20 - 19;
maxModuleiphi = ((sm - 1) % 18 + 1) * 20;
if (zsm == 1) {
minModuleieta = moduleBound[moduleInSm];
maxModuleieta = moduleBound[moduleInSm + 1] - 1;
} else if (zsm == -1) {
minModuleieta = -moduleBound[moduleInSm + 1] + 1;
maxModuleieta = -moduleBound[moduleInSm];
}
////////////////////////////nuovo ciclo
std::cout << "Called moduleIndex " << moduleIndex << std::endl;
std::cout << "minModuleieta " << minModuleieta << " maxModuleieta " << maxModuleieta << " minModuleiphi "
<< minModuleiphi << " maxModuleiphi " << maxModuleiphi << std::endl;
for (int ieta = minModuleieta; ieta <= maxModuleieta; ++ieta) {
for (int iphi = minModuleiphi; iphi <= maxModuleiphi; ++iphi) {
ringIds.push_back(EBDetId(ieta, iphi));
// std::cout<<"Putting Xtal with ieta: "<<ieta<<" iphi "<<iphi<<" of SM "<<sm<<" into Module "<<moduleIndex<<std::endl;
} //close loop on phi
} //close loop on eta
} //close if ( moduleInstallationNumber < 144)
return ringIds;
}
void EcalRingCalibrationTools::setCaloGeometry(const CaloGeometry* geometry) {
std::call_once(once_, EcalRingCalibrationTools::initializeFromGeometry, geometry);
}
void EcalRingCalibrationTools::initializeFromGeometry(CaloGeometry const* geometry) {
if (not geometry)
throw std::invalid_argument("EcalRingCalibrationTools::initializeFromGeometry called with a nullptr argument");
float cellPosEta[EEDetId::IX_MAX][EEDetId::IY_MAX];
for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
for (int iy = 0; iy < EEDetId::IY_MAX; ++iy) {
cellPosEta[ix][iy] = -1.;
endcapRingIndex_[ix][iy] = -9;
}
CaloSubdetectorGeometry const* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
if (not endcapGeometry)
throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry not found");
std::vector<DetId> const& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
for (std::vector<DetId>::const_iterator endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); ++endcapIt) {
EEDetId ee(*endcapIt);
if (ee.zside() == -1)
continue; //Just using +side to fill absEta x,y map
auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
int ics = ee.ix() - 1;
int ips = ee.iy() - 1;
cellPosEta[ics][ips] = fabs(cellGeometry->getPosition().eta());
//std::cout<<"EE Xtal, |eta| is "<<fabs(cellGeometry->getPosition().eta())<<std::endl;
}
float eta_ring[N_RING_ENDCAP / 2];
for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring)
eta_ring[ring] = cellPosEta[ring][50];
double etaBoundary[N_RING_ENDCAP / 2 + 1];
etaBoundary[0] = 1.47;
etaBoundary[N_RING_ENDCAP / 2] = 4.0;
for (int ring = 1; ring < N_RING_ENDCAP / 2; ++ring)
etaBoundary[ring] = (eta_ring[ring] + eta_ring[ring - 1]) / 2.;
for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring) {
// std::cout<<"***********************EE ring: "<<ring<<" eta "<<(etaBoundary[ring] + etaBoundary[ring+1])/2.<<std::endl;
for (int ix = 0; ix < EEDetId::IX_MAX; ix++)
for (int iy = 0; iy < EEDetId::IY_MAX; iy++)
if (cellPosEta[ix][iy] > etaBoundary[ring] && cellPosEta[ix][iy] < etaBoundary[ring + 1]) {
endcapRingIndex_[ix][iy] = ring;
//std::cout<<"endcapRing_["<<ix+1<<"]["<<iy+1<<"] = "<<ring<<";"<<std::endl;
}
}
std::vector<DetId> const& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
for (std::vector<DetId>::const_iterator barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); ++barrelIt) {
EBDetId eb(*barrelIt);
}
//EB
isInitializedFromGeometry_ = true;
}
|