TrackerTreeGenerator.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    TrackerTreeGenerator
// Class:      TrackerTreeGenerator
//
//
// Original Author:  Johannes Hauk
//         Created:  Fri Jan 16 14:09:52 CET 2009
//         Modified by: Gregor Mittag (DESY)
//
//
//

// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CommonTools/Utils/interface/TFileDirectory.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
#include "Geometry/Records/interface/PTrackerParametersRcd.h"
#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
#include "Geometry/CommonTopologies/interface/StripTopology.h"

#include "Alignment/TrackerAlignment/interface/TrackerTreeVariables.h"
#include "Alignment/TrackerAlignment/interface/AlignableTracker.h"

#include "TTree.h"
//
// class decleration
//

class TrackerTreeGenerator : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit TrackerTreeGenerator(const edm::ParameterSet&);
  ~TrackerTreeGenerator() override = default;

private:
  void beginJob() override;
  void analyze(const edm::Event&, const edm::EventSetup&) override;
  void endJob() override;

  // ----------member data ---------------------------

  const bool createEntryForDoubleSidedModule_;
  std::vector<TrackerTreeVariables> vTkTreeVar_;
  edm::ParameterSet config_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
TrackerTreeGenerator::TrackerTreeGenerator(const edm::ParameterSet& config)
    : createEntryForDoubleSidedModule_(config.getParameter<bool>("createEntryForDoubleSidedModule")), config_(config) {
  usesResource(TFileService::kSharedResource);
}

//
// member functions
//

// ------------ method called to for each event  ------------
void TrackerTreeGenerator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  //iSetup.get<TrackerDigiGeometryRecord>().get(tkGeom);
  // now try to take directly the ideal geometry independent of used geometry in Global Tag
  edm::ESHandle<GeometricDet> geometricDet;
  iSetup.get<IdealGeometryRecord>().get(geometricDet);

  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology* const tTopo = tTopoHandle.product();

  edm::ESHandle<PTrackerParameters> ptp;
  iSetup.get<PTrackerParametersRcd>().get(ptp);
  TrackerGeomBuilderFromGeometricDet trackerBuilder;
  const TrackerGeometry* tkGeom = trackerBuilder.build(&(*geometricDet), *ptp, tTopo);
  AlignableTracker alignableTracker{tkGeom, tTopo};
  const auto& ns = alignableTracker.trackerNameSpace();

  edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
                                       << "There are " << tkGeom->detIds().size() << " dets and "
                                       << tkGeom->detUnitIds().size() << " detUnits in the Geometry Record";

  if (createEntryForDoubleSidedModule_) {
    edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
                                         << "Create entry for each module AND one entry for virtual "
                                         << "double-sided module in addition";
  } else {
    edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::analyze"
                                         << "Create one entry for each physical module, do NOT create additional "
                                         << "entry for virtual double-sided module";
  }

  for (const auto& detId : tkGeom->detIds()) {
    const GeomDet& geomDet = *tkGeom->idToDet(detId);
    const Surface& surface = geomDet.surface();

    TrackerTreeVariables tkTreeVar;
    const auto rawId = detId.rawId();
    tkTreeVar.rawId = rawId;
    tkTreeVar.subdetId = detId.subdetId();

    switch (tkTreeVar.subdetId) {
      case PixelSubdetector::PixelBarrel:
        tkTreeVar.layer = tTopo->pxbLayer(detId);
        tkTreeVar.half = ns.tpb().halfBarrelNumber(rawId);
        tkTreeVar.rod = tTopo->pxbLadder(detId);  // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer!
        tkTreeVar.module = tTopo->pxbModule(detId);
        break;
      case PixelSubdetector::PixelEndcap:
        tkTreeVar.layer = tTopo->pxfDisk(detId);
        tkTreeVar.side = tTopo->pxfSide(detId);
        tkTreeVar.half = ns.tpe().halfCylinderNumber(rawId);
        tkTreeVar.blade = tTopo->pxfBlade(detId);
        tkTreeVar.panel = tTopo->pxfPanel(detId);
        tkTreeVar.module = tTopo->pxfModule(detId);
        break;
      case StripSubdetector::TIB:
        tkTreeVar.layer = tTopo->tibLayer(detId);
        tkTreeVar.side = tTopo->tibStringInfo(detId)[0];
        tkTreeVar.half = ns.tib().halfShellNumber(rawId);
        tkTreeVar.rod = tTopo->tibStringInfo(detId)[2];
        tkTreeVar.outerInner = tTopo->tibStringInfo(detId)[1];
        tkTreeVar.module = tTopo->tibModule(detId);
        tkTreeVar.isDoubleSide = tTopo->tibIsDoubleSide(detId);
        tkTreeVar.isRPhi = tTopo->tibIsRPhi(detId);
        tkTreeVar.isStereo = tTopo->tibIsStereo(detId);
        break;
      case StripSubdetector::TID:
        tkTreeVar.layer = tTopo->tidWheel(detId);
        tkTreeVar.side = tTopo->tidSide(detId);
        tkTreeVar.ring = tTopo->tidRing(detId);
        tkTreeVar.outerInner = tTopo->tidModuleInfo(detId)[0];
        tkTreeVar.module = tTopo->tidModuleInfo(detId)[1];
        tkTreeVar.isDoubleSide = tTopo->tidIsDoubleSide(detId);
        tkTreeVar.isRPhi = tTopo->tidIsRPhi(detId);
        tkTreeVar.isStereo = tTopo->tidIsStereo(detId);
        break;
      case StripSubdetector::TOB:
        tkTreeVar.layer = tTopo->tobLayer(detId);
        tkTreeVar.side = tTopo->tobRodInfo(detId)[0];
        tkTreeVar.rod = tTopo->tobRodInfo(detId)[1];
        tkTreeVar.module = tTopo->tobModule(detId);
        tkTreeVar.isDoubleSide = tTopo->tobIsDoubleSide(detId);
        tkTreeVar.isRPhi = tTopo->tobIsRPhi(detId);
        tkTreeVar.isStereo = tTopo->tobIsStereo(detId);
        break;
      case StripSubdetector::TEC:
        tkTreeVar.layer = tTopo->tecWheel(detId);
        tkTreeVar.side = tTopo->tecSide(detId);
        tkTreeVar.ring = tTopo->tecRing(detId);
        tkTreeVar.petal = tTopo->tecPetalInfo(detId)[1];
        tkTreeVar.outerInner = tTopo->tecPetalInfo(detId)[0];
        tkTreeVar.module = tTopo->tecModule(detId);
        tkTreeVar.isDoubleSide = tTopo->tecIsDoubleSide(detId);
        tkTreeVar.isRPhi = tTopo->tecIsRPhi(detId);
        tkTreeVar.isStereo = tTopo->tecIsStereo(detId);
        break;
    }

    LocalPoint lPModule(0., 0., 0.), lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
    GlobalPoint gPModule = surface.toGlobal(lPModule), gUDirection = surface.toGlobal(lUDirection),
                gVDirection = surface.toGlobal(lVDirection), gWDirection = surface.toGlobal(lWDirection);
    double dR(999.), dPhi(999.), dZ(999.);
    switch (tkTreeVar.subdetId) {
      case PixelSubdetector::PixelBarrel:
      case StripSubdetector::TIB:
      case StripSubdetector::TOB:
        dR = gWDirection.perp() - gPModule.perp();
        dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
        dZ = gVDirection.z() - gPModule.z();
        tkTreeVar.uDirection = dPhi > 0. ? 1 : -1;
        tkTreeVar.vDirection = dZ > 0. ? 1 : -1;
        tkTreeVar.wDirection = dR > 0. ? 1 : -1;
        break;
      case PixelSubdetector::PixelEndcap:
        dR = gUDirection.perp() - gPModule.perp();
        dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
        dZ = gWDirection.z() - gPModule.z();
        tkTreeVar.uDirection = dR > 0. ? 1 : -1;
        tkTreeVar.vDirection = dPhi > 0. ? 1 : -1;
        tkTreeVar.wDirection = dZ > 0. ? 1 : -1;
        break;
      case StripSubdetector::TID:
      case StripSubdetector::TEC:
        dR = gVDirection.perp() - gPModule.perp();
        dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
        dZ = gWDirection.z() - gPModule.z();
        tkTreeVar.uDirection = dPhi > 0. ? 1 : -1;
        tkTreeVar.vDirection = dR > 0. ? 1 : -1;
        tkTreeVar.wDirection = dZ > 0. ? 1 : -1;
        break;
    }
    tkTreeVar.posR = gPModule.perp();
    tkTreeVar.posPhi = gPModule.barePhi();  // = gPModule.barePhi().degrees();
    tkTreeVar.posEta = gPModule.eta();
    tkTreeVar.posX = gPModule.x();
    tkTreeVar.posY = gPModule.y();
    tkTreeVar.posZ = gPModule.z();

    if (auto stripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(&geomDet)) {  //is it a single physical module?
      switch (tkTreeVar.subdetId) {
        case StripSubdetector::TIB:
        case StripSubdetector::TOB:
        case StripSubdetector::TID:
        case StripSubdetector::TEC:
          auto& topol = dynamic_cast<const StripTopology&>(stripGeomDetUnit->specificTopology());
          tkTreeVar.nStrips = topol.nstrips();
          break;
      }
    }

    if (!createEntryForDoubleSidedModule_) {
      // do so only for individual modules and not also one entry for the combined doubleSided Module
      if (tkTreeVar.isDoubleSide)
        continue;
    }
    vTkTreeVar_.push_back(tkTreeVar);
  }
}

// ------------ method called once each job just before starting event loop  ------------
void TrackerTreeGenerator::beginJob() {}

// ------------ method called once each job just after ending the event loop  ------------
void TrackerTreeGenerator::endJob() {
  UInt_t rawId(999), subdetId(999), layer(999), side(999), half(999), rod(999), ring(999), petal(999), blade(999),
      panel(999), outerInner(999), module(999), nStrips(999);
  Bool_t isDoubleSide(false), isRPhi(false), isStereo(false);
  Int_t uDirection(999), vDirection(999), wDirection(999);
  Float_t posR(999.F), posPhi(999.F), posEta(999.F), posX(999.F), posY(999.F), posZ(999.F);

  edm::Service<TFileService> fileService;
  TFileDirectory treeDir = fileService->mkdir("TrackerTree");
  auto trackerTree{treeDir.make<TTree>("TrackerTree", "IDs of all modules (ideal geometry)")};
  trackerTree->Branch("RawId", &rawId, "RawId/i");
  trackerTree->Branch("SubdetId", &subdetId, "SubdetId/i");
  trackerTree->Branch("Layer", &layer, "Layer/i");                 // Barrel: Layer, Forward: Disk
  trackerTree->Branch("Side", &side, "Side/i");                    // Rod/Ring in +z or -z
  trackerTree->Branch("Half", &half, "Half/i");                    // PXB: HalfBarrel, PXF: HalfCylinder, TIB: HalfShell
  trackerTree->Branch("Rod", &rod, "Rod/i");                       // Barrel (Ladder or String or Rod)
  trackerTree->Branch("Ring", &ring, "Ring/i");                    // Forward
  trackerTree->Branch("Petal", &petal, "Petal/i");                 // TEC
  trackerTree->Branch("Blade", &blade, "Blade/i");                 // PXF
  trackerTree->Branch("Panel", &panel, "Panel/i");                 // PXF
  trackerTree->Branch("OuterInner", &outerInner, "OuterInner/i");  // front/back String,Ring,Petal
  trackerTree->Branch("Module", &module, "Module/i");              // Module ID
  trackerTree->Branch("NStrips", &nStrips, "NStrips/i");
  trackerTree->Branch("IsDoubleSide", &isDoubleSide, "IsDoubleSide/O");
  trackerTree->Branch("IsRPhi", &isRPhi, "IsRPhi/O");
  trackerTree->Branch("IsStereo", &isStereo, "IsStereo/O");
  trackerTree->Branch("UDirection", &uDirection, "UDirection/I");
  trackerTree->Branch("VDirection", &vDirection, "VDirection/I");
  trackerTree->Branch("WDirection", &wDirection, "WDirection/I");
  trackerTree->Branch("PosR", &posR, "PosR/F");
  trackerTree->Branch("PosPhi", &posPhi, "PosPhi/F");
  trackerTree->Branch("PosEta", &posEta, "PosEta/F");
  trackerTree->Branch("PosX", &posX, "PosX/F");
  trackerTree->Branch("PosY", &posY, "PosY/F");
  trackerTree->Branch("PosZ", &posZ, "PosZ/F");

  for (const auto& iTree : vTkTreeVar_) {
    rawId = iTree.rawId;
    subdetId = iTree.subdetId;
    layer = iTree.layer;
    side = iTree.side;
    half = iTree.half;
    rod = iTree.rod;
    ring = iTree.ring;
    petal = iTree.petal;
    blade = iTree.blade;
    panel = iTree.panel;
    outerInner = iTree.outerInner;
    module = iTree.module;
    nStrips = iTree.nStrips;
    isDoubleSide = iTree.isDoubleSide;
    isRPhi = iTree.isRPhi;
    isStereo = iTree.isStereo;
    uDirection = iTree.uDirection;
    vDirection = iTree.vDirection;
    wDirection = iTree.wDirection;
    posR = iTree.posR;
    posPhi = iTree.posPhi;
    posEta = iTree.posEta;
    posX = iTree.posX;
    posY = iTree.posY;
    posZ = iTree.posZ;

    trackerTree->Fill();
  }
  edm::LogInfo("TrackerTreeGenerator") << "@SUB=TrackerTreeGenerator::endJob"
                                       << "TrackerTree contains " << vTkTreeVar_.size() << " entries overall";
}

//define this as a plug-in
DEFINE_FWK_MODULE(TrackerTreeGenerator);