TrackerAlignment/plugins/TrackerTreeGenerator.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackerTreeGenerator
0004 // Class:      TrackerTreeGenerator
0005 //
0006 /**\class TrackerTreeGenerator TrackerTreeGenerator.cc Alignment/TrackerAlignment/plugins/TrackerTreeGenerator.cc
0007  Description: <one line class summary>
0008  Implementation:
0009      <Notes on implementation>
0010 */
0011 //
0012 // Original Author:  Johannes Hauk
0013 //         Created:  Fri Jan 16 14:09:52 CET 2009
0014 //         Modified by: Gregor Mittag (DESY)
0015 //
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031
0032 #include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034 #include "CommonTools/Utils/interface/TFileDirectory.h"
0035
0036 #include "DataFormats/DetId/interface/DetId.h"
0037 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0038 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0039 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0040 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0041 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0042 #include "DataFormats/GeometrySurface/interface/Surface.h"
0043 #include "DataFormats/Math/interface/deltaPhi.h"
0044
0045 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0046 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0047 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0048 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0049 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
0050 #include "Geometry/Records/interface/PTrackerParametersRcd.h"
0051 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0052 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0053 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0054 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0055
0056 #include "Alignment/TrackerAlignment/interface/TrackerTreeVariables.h"
0057 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0058
0059 #include "TTree.h"
0060 //
0061 // class decleration
0062 //
0063
0064 class TrackerTreeGenerator : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0065 public:
0066   explicit TrackerTreeGenerator(const edm::ParameterSet&);
0067   ~TrackerTreeGenerator() override = default;
0068
0069 private:
0070   void beginJob() override;
0071   void analyze(const edm::Event&, const edm::EventSetup&) override;
0072   void endJob() override;
0073
0074   // ----------member data ---------------------------
0075   const edm::ESGetToken<GeometricDet, IdealGeometryRecord> geomDetToken_;
0076   const edm::ESGetToken<PTrackerParameters, PTrackerParametersRcd> ptpToken_;
0077   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0078
0079   const bool createEntryForDoubleSidedModule_;
0080   std::vector<TrackerTreeVariables> vTkTreeVar_;
0081   edm::ParameterSet config_;
0082 };
0083
0084 //
0085 // constants, enums and typedefs
0086 //
0087
0088 //
0089 // static data member definitions
0090 //
0091
0092 //
0093 // constructors and destructor
0094 //
0095 TrackerTreeGenerator::TrackerTreeGenerator(const edm::ParameterSet& config)
0096     : geomDetToken_(esConsumes()),
0097       ptpToken_(esConsumes()),
0098       topoToken_(esConsumes()),
0099       createEntryForDoubleSidedModule_(config.getParameter<bool>("createEntryForDoubleSidedModule")),
0100       config_(config) {
0101   usesResource(TFileService::kSharedResource);
0102 }
0103
0104 //
0105 // member functions
0106 //
0107
0108 // ------------ method called to for each event  ------------
0109 void TrackerTreeGenerator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110   // now try to take directly the ideal geometry independent of used geometry in Global Tag
0111   const GeometricDet* geometricDet = &iSetup.getData(geomDetToken_);
0112   const PTrackerParameters& ptp = iSetup.getData(ptpToken_);
0113   const TrackerTopology* tTopo = &iSetup.getData(topoToken_);
0114
0115   TrackerGeomBuilderFromGeometricDet trackerBuilder;
0116   const TrackerGeometry* tkGeom = trackerBuilder.build(geometricDet, ptp, tTopo);
0117   AlignableTracker alignableTracker{tkGeom, tTopo};
0118   const auto& ns = alignableTracker.trackerNameSpace();
0119
0120   edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
0121                                        << "There are " << tkGeom->detIds().size() << " dets and "
0122                                        << tkGeom->detUnitIds().size() << " detUnits in the Geometry Record";
0123
0124   if (createEntryForDoubleSidedModule_) {
0125     edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
0126                                          << "Create entry for each module AND one entry for virtual "
0127                                          << "double-sided module in addition";
0128   } else {
0129     edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
0130                                          << "Create one entry for each physical module, do NOT create additional "
0131                                          << "entry for virtual double-sided module";
0132   }
0133
0134   for (const auto& detId : tkGeom->detIds()) {
0135     const GeomDet& geomDet = *tkGeom->idToDet(detId);
0136     const Surface& surface = geomDet.surface();
0137
0138     TrackerTreeVariables tkTreeVar;
0139     const auto rawId = detId.rawId();
0140     tkTreeVar.rawId = rawId;
0141     tkTreeVar.subdetId = detId.subdetId();
0142
0143     switch (tkTreeVar.subdetId) {
0144       case PixelSubdetector::PixelBarrel:
0145         tkTreeVar.layer = tTopo->pxbLayer(detId);
0146         tkTreeVar.half = ns.tpb().halfBarrelNumber(rawId);
0147         tkTreeVar.rod = tTopo->pxbLadder(detId);  // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer!
0148         tkTreeVar.module = tTopo->pxbModule(detId);
0149         break;
0150       case PixelSubdetector::PixelEndcap:
0151         tkTreeVar.layer = tTopo->pxfDisk(detId);
0152         tkTreeVar.side = tTopo->pxfSide(detId);
0153         tkTreeVar.half = ns.tpe().halfCylinderNumber(rawId);
0154         tkTreeVar.blade = tTopo->pxfBlade(detId);
0155         tkTreeVar.panel = tTopo->pxfPanel(detId);
0156         tkTreeVar.module = tTopo->pxfModule(detId);
0157         break;
0158       case StripSubdetector::TIB:
0159         tkTreeVar.layer = tTopo->tibLayer(detId);
0160         tkTreeVar.side = tTopo->tibStringInfo(detId)[0];
0161         tkTreeVar.half = ns.tib().halfShellNumber(rawId);
0162         tkTreeVar.rod = tTopo->tibStringInfo(detId)[2];
0163         tkTreeVar.outerInner = tTopo->tibStringInfo(detId)[1];
0164         tkTreeVar.module = tTopo->tibModule(detId);
0165         tkTreeVar.isDoubleSide = tTopo->tibIsDoubleSide(detId);
0166         tkTreeVar.isRPhi = tTopo->tibIsRPhi(detId);
0167         tkTreeVar.isStereo = tTopo->tibIsStereo(detId);
0168         break;
0169       case StripSubdetector::TID:
0170         tkTreeVar.layer = tTopo->tidWheel(detId);
0171         tkTreeVar.side = tTopo->tidSide(detId);
0172         tkTreeVar.ring = tTopo->tidRing(detId);
0173         tkTreeVar.outerInner = tTopo->tidModuleInfo(detId)[0];
0174         tkTreeVar.module = tTopo->tidModuleInfo(detId)[1];
0175         tkTreeVar.isDoubleSide = tTopo->tidIsDoubleSide(detId);
0176         tkTreeVar.isRPhi = tTopo->tidIsRPhi(detId);
0177         tkTreeVar.isStereo = tTopo->tidIsStereo(detId);
0178         break;
0179       case StripSubdetector::TOB:
0180         tkTreeVar.layer = tTopo->tobLayer(detId);
0181         tkTreeVar.side = tTopo->tobRodInfo(detId)[0];
0182         tkTreeVar.rod = tTopo->tobRodInfo(detId)[1];
0183         tkTreeVar.module = tTopo->tobModule(detId);
0184         tkTreeVar.isDoubleSide = tTopo->tobIsDoubleSide(detId);
0185         tkTreeVar.isRPhi = tTopo->tobIsRPhi(detId);
0186         tkTreeVar.isStereo = tTopo->tobIsStereo(detId);
0187         break;
0188       case StripSubdetector::TEC:
0189         tkTreeVar.layer = tTopo->tecWheel(detId);
0190         tkTreeVar.side = tTopo->tecSide(detId);
0191         tkTreeVar.ring = tTopo->tecRing(detId);
0192         tkTreeVar.petal = tTopo->tecPetalInfo(detId)[1];
0193         tkTreeVar.outerInner = tTopo->tecPetalInfo(detId)[0];
0194         tkTreeVar.module = tTopo->tecModule(detId);
0195         tkTreeVar.isDoubleSide = tTopo->tecIsDoubleSide(detId);
0196         tkTreeVar.isRPhi = tTopo->tecIsRPhi(detId);
0197         tkTreeVar.isStereo = tTopo->tecIsStereo(detId);
0198         break;
0199     }
0200
0201     LocalPoint lPModule(0., 0., 0.), lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
0202     GlobalPoint gPModule = surface.toGlobal(lPModule), gUDirection = surface.toGlobal(lUDirection),
0203                 gVDirection = surface.toGlobal(lVDirection), gWDirection = surface.toGlobal(lWDirection);
0204     double dR(999.), dPhi(999.), dZ(999.);
0205     switch (tkTreeVar.subdetId) {
0206       case PixelSubdetector::PixelBarrel:
0207       case StripSubdetector::TIB:
0208       case StripSubdetector::TOB:
0209         dR = gWDirection.perp() - gPModule.perp();
0210         dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
0211         dZ = gVDirection.z() - gPModule.z();
0212         tkTreeVar.uDirection = dPhi > 0. ? 1 : -1;
0213         tkTreeVar.vDirection = dZ > 0. ? 1 : -1;
0214         tkTreeVar.wDirection = dR > 0. ? 1 : -1;
0215         break;
0216       case PixelSubdetector::PixelEndcap:
0217         dR = gUDirection.perp() - gPModule.perp();
0218         dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
0219         dZ = gWDirection.z() - gPModule.z();
0220         tkTreeVar.uDirection = dR > 0. ? 1 : -1;
0221         tkTreeVar.vDirection = dPhi > 0. ? 1 : -1;
0222         tkTreeVar.wDirection = dZ > 0. ? 1 : -1;
0223         break;
0224       case StripSubdetector::TID:
0225       case StripSubdetector::TEC:
0226         dR = gVDirection.perp() - gPModule.perp();
0227         dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
0228         dZ = gWDirection.z() - gPModule.z();
0229         tkTreeVar.uDirection = dPhi > 0. ? 1 : -1;
0230         tkTreeVar.vDirection = dR > 0. ? 1 : -1;
0231         tkTreeVar.wDirection = dZ > 0. ? 1 : -1;
0232         break;
0233     }
0234     tkTreeVar.posR = gPModule.perp();
0235     tkTreeVar.posPhi = gPModule.barePhi();  // = gPModule.barePhi().degrees();
0236     tkTreeVar.posEta = gPModule.eta();
0237     tkTreeVar.posX = gPModule.x();
0238     tkTreeVar.posY = gPModule.y();
0239     tkTreeVar.posZ = gPModule.z();
0240
0241     if (auto stripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(&geomDet)) {  //is it a single physical module?
0242       switch (tkTreeVar.subdetId) {
0243         case StripSubdetector::TIB:
0244         case StripSubdetector::TOB:
0245         case StripSubdetector::TID:
0246         case StripSubdetector::TEC:
0247           auto& topol = dynamic_cast<const StripTopology&>(stripGeomDetUnit->specificTopology());
0248           tkTreeVar.nStrips = topol.nstrips();
0249           break;
0250       }
0251     }
0252
0253     if (!createEntryForDoubleSidedModule_) {
0254       // do so only for individual modules and not also one entry for the combined doubleSided Module
0255       if (tkTreeVar.isDoubleSide)
0256         continue;
0257     }
0258     vTkTreeVar_.push_back(tkTreeVar);
0259   }
0260 }
0261
0262 // ------------ method called once each job just before starting event loop  ------------
0263 void TrackerTreeGenerator::beginJob() {}
0264
0265 // ------------ method called once each job just after ending the event loop  ------------
0266 void TrackerTreeGenerator::endJob() {
0267   UInt_t rawId(999), subdetId(999), layer(999), side(999), half(999), rod(999), ring(999), petal(999), blade(999),
0268       panel(999), outerInner(999), module(999), nStrips(999);
0269   Bool_t isDoubleSide(false), isRPhi(false), isStereo(false);
0270   Int_t uDirection(999), vDirection(999), wDirection(999);
0271   Float_t posR(999.F), posPhi(999.F), posEta(999.F), posX(999.F), posY(999.F), posZ(999.F);
0272
0273   edm::Service<TFileService> fileService;
0274   TFileDirectory treeDir = fileService->mkdir("TrackerTree");
0275   auto trackerTree{treeDir.make<TTree>("TrackerTree", "IDs of all modules (ideal geometry)")};
0276   trackerTree->Branch("RawId", &rawId, "RawId/i");
0277   trackerTree->Branch("SubdetId", &subdetId, "SubdetId/i");
0278   trackerTree->Branch("Layer", &layer, "Layer/i");                 // Barrel: Layer, Forward: Disk
0279   trackerTree->Branch("Side", &side, "Side/i");                    // Rod/Ring in +z or -z
0280   trackerTree->Branch("Half", &half, "Half/i");                    // PXB: HalfBarrel, PXF: HalfCylinder, TIB: HalfShell
0281   trackerTree->Branch("Rod", &rod, "Rod/i");                       // Barrel (Ladder or String or Rod)
0282   trackerTree->Branch("Ring", &ring, "Ring/i");                    // Forward
0283   trackerTree->Branch("Petal", &petal, "Petal/i");                 // TEC
0284   trackerTree->Branch("Blade", &blade, "Blade/i");                 // PXF
0285   trackerTree->Branch("Panel", &panel, "Panel/i");                 // PXF
0286   trackerTree->Branch("OuterInner", &outerInner, "OuterInner/i");  // front/back String,Ring,Petal
0287   trackerTree->Branch("Module", &module, "Module/i");              // Module ID
0288   trackerTree->Branch("NStrips", &nStrips, "NStrips/i");
0289   trackerTree->Branch("IsDoubleSide", &isDoubleSide, "IsDoubleSide/O");
0290   trackerTree->Branch("IsRPhi", &isRPhi, "IsRPhi/O");
0291   trackerTree->Branch("IsStereo", &isStereo, "IsStereo/O");
0292   trackerTree->Branch("UDirection", &uDirection, "UDirection/I");
0293   trackerTree->Branch("VDirection", &vDirection, "VDirection/I");
0294   trackerTree->Branch("WDirection", &wDirection, "WDirection/I");
0295   trackerTree->Branch("PosR", &posR, "PosR/F");
0296   trackerTree->Branch("PosPhi", &posPhi, "PosPhi/F");
0297   trackerTree->Branch("PosEta", &posEta, "PosEta/F");
0298   trackerTree->Branch("PosX", &posX, "PosX/F");
0299   trackerTree->Branch("PosY", &posY, "PosY/F");
0300   trackerTree->Branch("PosZ", &posZ, "PosZ/F");
0301
0302   for (const auto& iTree : vTkTreeVar_) {
0303     rawId = iTree.rawId;
0304     subdetId = iTree.subdetId;
0305     layer = iTree.layer;
0306     side = iTree.side;
0307     half = iTree.half;
0308     rod = iTree.rod;
0309     ring = iTree.ring;
0310     petal = iTree.petal;
0311     blade = iTree.blade;
0312     panel = iTree.panel;
0313     outerInner = iTree.outerInner;
0314     module = iTree.module;
0315     nStrips = iTree.nStrips;
0316     isDoubleSide = iTree.isDoubleSide;
0317     isRPhi = iTree.isRPhi;
0318     isStereo = iTree.isStereo;
0319     uDirection = iTree.uDirection;
0320     vDirection = iTree.vDirection;
0321     wDirection = iTree.wDirection;
0322     posR = iTree.posR;
0323     posPhi = iTree.posPhi;
0324     posEta = iTree.posEta;
0325     posX = iTree.posX;
0326     posY = iTree.posY;
0327     posZ = iTree.posZ;
0328
0329     trackerTree->Fill();
0330   }
0331   edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::endJob"
0332                                        << "TrackerTree contains " << vTkTreeVar_.size() << " entries overall";
0333 }
0334
0335 //define this as a plug-in
0336 DEFINE_FWK_MODULE(TrackerTreeGenerator);