Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-12 00:01:22

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackerMapTool
0004 // Class:      TrackerMapTool
0005 //
0006 /**\class TrackerMapTool TrackerMapTool.cc 
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Filippo Ambroglini
0015 //         Created:  Tue Jul 26 08:47:57 CEST 2005
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 #include <vector>
0023 #include <fstream>
0024 #include <sstream>
0025 #include <string>
0026 #include <vector>
0027 // user include files
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 // class decleration
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 // constructors and destructor
0064 //
0065 TrackerMapTool::TrackerMapTool(const edm::ParameterSet& iConfig) : pDDToken_(esConsumes()) {}
0066 
0067 //
0068 // member functions
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   // 2009-08-26 Michael Case: to get rid of a compiler warning about control reaching
0091   // the end of a non-void function I put return -1 here.  This changed the output
0092   // of the test so I changed it to return 0 to match the "before my changes" run
0093   // of trackerMap_cfg.py.
0094   return 0;  // this was added.  No checks have been mad where layerno is used.
0095 }
0096 // ------------ method called to produce the data  ------------
0097 void TrackerMapTool::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0098   edm::LogInfo("TrackerMapTool") << "Here I am";
0099 
0100   //
0101   // get the TrackerGeom
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) {  //barrel
0133       layer = ((*begin)->geographicalId().rawId() >> 16) & 0xF;
0134       leftright = 0;
0135       name0 = " ";
0136     } else {  //endcap
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;        // cm -> m
0260     width = (*begin)->surface().bounds().width() / 100.0;          // cm -> mo
0261     thickness = (*begin)->surface().bounds().thickness() / 100.0;  // cm -> m
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;  // cm -> m
0274     float posy = (*begin)->surface().position().y() / 100.0;  // cm -> m
0275     float posz = (*begin)->surface().position().z() / 100.0;  // cm -> m
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) {  //barrel
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 {  // endcap
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     }  //end of endcap part
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 //define this as a plug-in
0369 DEFINE_FWK_MODULE(TrackerMapTool);