File indexing completed on 2023-10-25 09:50:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022 #include <vector>
0023 #include <fstream>
0024 #include <sstream>
0025 #include <string>
0026 #include <vector>
0027
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0037 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0038 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0039 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0040 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0041 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0043
0044 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0045 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
0046 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0047
0048
0049
0050
0051
0052 class TrackerMapTool : public edm::one::EDAnalyzer<> {
0053 public:
0054 explicit TrackerMapTool(const edm::ParameterSet&);
0055
0056 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0057
0058 private:
0059 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> pDDToken_;
0060 };
0061
0062
0063
0064
0065 TrackerMapTool::TrackerMapTool(const edm::ParameterSet& iConfig) : pDDToken_(esConsumes()) {}
0066
0067
0068
0069
0070
0071 int layerno(int subdet, int leftright, int layer) {
0072 if (subdet == 6 && leftright == 1)
0073 return (10 - layer);
0074 if (subdet == 6 && leftright == 2)
0075 return (layer + 21);
0076 if (subdet == 4 && leftright == 1)
0077 return (4 - layer + 9);
0078 if (subdet == 4 && leftright == 2)
0079 return (layer + 18);
0080 if (subdet == 2 && leftright == 1)
0081 return (4 - layer + 12);
0082 if (subdet == 2 && leftright == 2)
0083 return (layer + 15);
0084 if (subdet == 1)
0085 return (layer + 30);
0086 if (subdet == 3)
0087 return (layer + 33);
0088 if (subdet == 5)
0089 return (layer + 37);
0090
0091
0092
0093
0094 return 0;
0095 }
0096
0097 void TrackerMapTool::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0098 edm::LogInfo("TrackerMapTool") << "Here I am";
0099
0100
0101
0102
0103 auto const& pDD = iSetup.getData(pDDToken_);
0104 edm::LogInfo("TrackerMapTool") << " Geometry node for TrackerGeom is " << &pDD;
0105 edm::LogInfo("TrackerMapTool") << " I have " << pDD.dets().size() << " detectors";
0106 edm::LogInfo("TrackerMapTool") << " I have " << pDD.detTypes().size() << " types";
0107 int spicchif[] = {24, 24, 40, 56, 40, 56, 80};
0108 int spicchib[] = {20, 32, 44, 30, 38, 46, 56, 42, 48, 54, 60, 66, 74};
0109
0110 int nlay = 0;
0111 int layer, subdet, leftright = 0, ringno, petalno, moduleno, isStereo, pixel_strip, barrel_forward;
0112 std::string name0, name1, name2, name3, name4, name5;
0113 int ring = 0;
0114 int nmod = 0;
0115 int ntotmod = 0;
0116 float r;
0117 int bar_fow = 1;
0118 float phi, phi1;
0119 float rmedioS[] = {0.27665, 0.3671, 0.4474, 0.5617, 0.6768, 0.8189, 0.9907};
0120 float rmedioP[] = {0.0623081, 0.074111, 0.0870344, 0.103416, 0.115766, 0.132728, 0.140506};
0121 std::string nameDet;
0122 int last_layer = 0;
0123 float width, length, thickness, widthAtHalfLength;
0124 std::ofstream output("tracker.dat", std::ios::out);
0125
0126 auto begin = pDD.detUnits().begin();
0127 auto end = pDD.detUnits().end();
0128
0129 for (; begin != end; ++begin) {
0130 ntotmod++;
0131 subdet = (*begin)->geographicalId().subdetId();
0132 if (subdet == 1 || subdet == 3 || subdet == 5) {
0133 layer = ((*begin)->geographicalId().rawId() >> 16) & 0xF;
0134 leftright = 0;
0135 name0 = " ";
0136 } else {
0137 leftright = ((*begin)->geographicalId().rawId() >> 23) & 0x3;
0138 name0 = "+z";
0139 if (leftright == 1)
0140 name0 = "-z";
0141 layer = ((*begin)->geographicalId().rawId() >> 16) & 0xF;
0142 }
0143 isStereo = (*begin)->geographicalId().rawId() & 0x3;
0144 pixel_strip = 2;
0145 if (subdet <= 2)
0146 pixel_strip = 1;
0147 barrel_forward = 2;
0148 if (subdet == 2 || subdet == 4 || subdet == 6) {
0149 if (leftright == 1)
0150 barrel_forward = 1;
0151 else
0152 barrel_forward = 3;
0153 }
0154 nlay = layerno(subdet, leftright, layer);
0155 ringno = 0;
0156 petalno = 0;
0157 moduleno = 0;
0158 name1 = " ";
0159 if (subdet == 1) {
0160 nameDet = "BPIX";
0161 name1 = " ladder ";
0162 ringno = ((*begin)->geographicalId().rawId() >> 8) & 0xFF;
0163 }
0164 if (subdet == 2) {
0165 nameDet = "FPIX";
0166 name1 = " ring ";
0167 ringno = ((*begin)->geographicalId().rawId() >> 8) & 0xFF;
0168 }
0169 if (subdet == 3) {
0170 nameDet = "TIB";
0171 name1 = " string ";
0172 ringno = ((*begin)->geographicalId().rawId() >> 8) & 0x3F;
0173 }
0174 if (subdet == 5) {
0175 nameDet = "TOB";
0176 name1 = " rod ";
0177 ringno = ((*begin)->geographicalId().rawId() >> 8) & 0x7F;
0178 }
0179 if (subdet == 4) {
0180 nameDet = "TID";
0181 name1 = " ring ";
0182 ringno = ((*begin)->geographicalId().rawId() >> 8) & 0xFF;
0183 }
0184 if (subdet == 6) {
0185 nameDet = "TEC";
0186 name1 = " ring ";
0187 }
0188 name2 = " ";
0189 if (subdet == 6) {
0190 name2 = " petal ";
0191 petalno = ((*begin)->geographicalId().rawId() >> 8) & 0x7F;
0192 }
0193 if (subdet == 2) {
0194 name2 = " blade ";
0195 petalno = ((*begin)->geographicalId().rawId() >> 8) & 0x3F;
0196 }
0197 name3 = " ";
0198 if (subdet == 6) {
0199 name3 = " forward ";
0200 if ((((*begin)->geographicalId().rawId() >> 4) & 0x1) == 1)
0201 name3 = " backward ";
0202 }
0203 if (subdet == 2) {
0204 name3 = " forward ";
0205 if ((((*begin)->geographicalId().rawId() >> 4) & 0x1) == 1)
0206 name3 = " backward ";
0207 }
0208 name5 = " ";
0209 if (subdet == 6) {
0210 name5 = " forward ";
0211 if ((((*begin)->geographicalId().rawId() >> 15) & 0x1) == 1) {
0212 name5 = " backward ";
0213 }
0214 }
0215 if (subdet == 2) {
0216 name5 = " left ";
0217 if ((((*begin)->geographicalId().rawId() >> 14) & 0x1) == 1) {
0218 name5 = " right ";
0219 }
0220 }
0221 if (subdet == 3) {
0222 name2 = " neg ";
0223 if ((((*begin)->geographicalId().rawId() >> 15) & 0x1) == 1)
0224 name2 = " pos ";
0225 }
0226 if (subdet == 5) {
0227 name2 = " neg ";
0228 if ((((*begin)->geographicalId().rawId() >> 15) & 0x1) == 1)
0229 name2 = " pos ";
0230 }
0231 if (subdet == 3) {
0232 name3 = " internal ";
0233 if ((((*begin)->geographicalId().rawId() >> 14) & 0x1) == 1)
0234 name3 = " external ";
0235 }
0236 if (subdet == 4) {
0237 name3 = " forward ";
0238 if ((((*begin)->geographicalId().rawId() >> 7) & 0x1) == 1)
0239 name3 = " backward ";
0240 }
0241 if (subdet == 1) {
0242 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x3F;
0243 }
0244 if (subdet == 3) {
0245 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x3F;
0246 }
0247 if (subdet == 5) {
0248 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x3F;
0249 }
0250 if (subdet == 4) {
0251 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x1F;
0252 }
0253 if (subdet == 6) {
0254 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x3;
0255 }
0256 if (subdet == 2) {
0257 moduleno = ((*begin)->geographicalId().rawId() >> 2) & 0x7;
0258 }
0259 length = (*begin)->surface().bounds().length() / 100.0;
0260 width = (*begin)->surface().bounds().width() / 100.0;
0261 thickness = (*begin)->surface().bounds().thickness() / 100.0;
0262 widthAtHalfLength = (*begin)->surface().bounds().widthAtHalfLength() / 100.0;
0263
0264 if (nlay != last_layer) {
0265 ring = -1;
0266 last_layer = nlay;
0267 }
0268 bar_fow = 2;
0269 if (nlay < 16)
0270 bar_fow = 1;
0271 if (nlay > 15 && nlay < 31)
0272 bar_fow = 3;
0273 float posx = (*begin)->surface().position().x() / 100.0;
0274 float posy = (*begin)->surface().position().y() / 100.0;
0275 float posz = (*begin)->surface().position().z() / 100.0;
0276 r = pow(((double)(posx * posx) + posy * posy), 0.5);
0277 phi1 = atan(posy / posx);
0278 phi = phi1;
0279 if (posy < 0. && posx > 0)
0280 phi = phi1 + 2. * M_PI;
0281 if (posx < 0.)
0282 phi = phi1 + M_PI;
0283 if (fabs(posy) < 0.000001 && posx > 0)
0284 phi = 0;
0285 if (fabs(posy) < 0.000001 && posx < 0)
0286 phi = M_PI;
0287 if (fabs(posx) < 0.000001 && posy > 0)
0288 phi = M_PI / 2.;
0289 if (fabs(posx) < 0.000001 && posy < 0)
0290 phi = M_PI + M_PI / 2.;
0291
0292 if (bar_fow == 2) {
0293 if (subdet == 1)
0294 ring = moduleno;
0295 if (subdet == 5) {
0296 ring = moduleno;
0297 if ((((*begin)->geographicalId().rawId() >> 15) & 0x1) == 1)
0298 ring = ring + 6;
0299 }
0300 if (subdet == 3) {
0301 ring = moduleno;
0302 if (layer == 2 || layer == 4) {
0303 if ((((*begin)->geographicalId().rawId() >> 14) & 0x1) == 1) {
0304 ring = ring * 2;
0305 } else {
0306 ring = ring * 2 - 1;
0307 }
0308 } else {
0309 if ((((*begin)->geographicalId().rawId() >> 14) & 0x1) == 1) {
0310 ring = ring * 2 - 1;
0311 } else {
0312 ring = ring * 2;
0313 }
0314 }
0315 if ((((*begin)->geographicalId().rawId() >> 15) & 0x1) == 1)
0316 ring = ring + 6;
0317 }
0318
0319 nmod = (int)((phi / (2. * M_PI)) * spicchib[nlay - 31] + .1) + 1;
0320 if (nlay == 40)
0321 nmod = nmod - 1;
0322 if (subdet == 1)
0323 nmod = ringno;
0324 } else {
0325 if (subdet == 4 || subdet == 6) {
0326 for (int i = 0; i < 7; i++) {
0327 if (fabs(r - rmedioS[i]) < 0.015) {
0328 ring = i + 1;
0329 break;
0330 }
0331 }
0332 nmod = (int)((phi / (2. * M_PI)) * spicchif[ring - 1] + .1) + 1;
0333 } else {
0334 for (int i = 0; i < 7; i++) {
0335 if (fabs(r - rmedioP[i]) < 0.0001) {
0336 ring = i + 1;
0337 break;
0338 }
0339 }
0340 nmod = (int)((phi / (2. * M_PI)) * 24 + .1) + 1;
0341 }
0342 }
0343 if (isStereo == 1)
0344 nmod = nmod + 100;
0345 name4 = " ";
0346 if (isStereo == 1)
0347 name4 = " stereo";
0348 std::ostringstream outs;
0349 if (subdet == 6)
0350 outs << nameDet << " " << name0 << " layer " << layer << name1 << ringno << name2 << petalno << name5
0351 << " module " << moduleno << name3 << " " << name4;
0352 if (subdet == 2)
0353 outs << nameDet << " " << name0 << " disc " << layer << name1 << ring << name2 << petalno << " plaquette "
0354 << moduleno << name5;
0355 if (subdet == 1 || subdet == 3 || subdet == 4 || subdet == 5)
0356 outs << nameDet << " " << name0 << " layer " << layer << name1 << ringno << name2 << name3 << " module "
0357 << moduleno << " " << name4;
0358 char buffer[20];
0359 sprintf(buffer, "%X", (*begin)->geographicalId().rawId());
0360 int id = (*begin)->geographicalId().rawId();
0361 output << ntotmod << " " << pixel_strip << " " << barrel_forward << " " << nlay << " " << ring << " " << nmod << " "
0362 << posx << " " << posy << " " << posz << " " << length << " " << width << " " << thickness << " "
0363 << widthAtHalfLength << " "
0364 << " " << id << std::endl;
0365 output << outs.str() << std::endl;
0366 }
0367 }
0368
0369 DEFINE_FWK_MODULE(TrackerMapTool);