MuonTools/plugins/MuDTSegmentExtTableProducer.cc

0001 /** \class MuDTSegmentExtTableProducer MuDTSegmentExtTableProducer.ccDPGAnalysis/MuonTools/src/MuDTSegmentExtTableProducer.cc
0002  *
0003  * Helper class : the segment TableProducer for Phase-1 / Phase2 segments (the DataFormat is the same)
0004  *
0005  * \author C. Battilana (INFN BO)
0006  *
0007 *
0008 */
0009
0010 #include "DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h"
0011
0012 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0013 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0014
0015 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
0016
0017 #include <vector>
0018 #include <array>
0019
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0022
0023 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0024 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0025
0026 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0027 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0028
0029 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0030 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0031
0032 class MuDTSegmentExtTableProducer : public MuBaseFlatTableProducer {
0033 public:
0034   /// Constructor
0035   MuDTSegmentExtTableProducer(const edm::ParameterSet&);
0036
0037   /// Fill descriptors
0038   static void fillDescriptions(edm::ConfigurationDescriptions&);
0039
0040 protected:
0041   /// Fill tree branches for a given event
0042   void fillTable(edm::Event&) final;
0043
0044   /// Get info from the ES by run
0045   void getFromES(const edm::Run&, const edm::EventSetup&) final;
0046
0047   /// Get info from the ES for a given event
0048   void getFromES(const edm::EventSetup&) final;
0049
0050 private:
0051   static const int FIRST_LAYER{1};
0052   static const int LAST_LAYER{4};
0053   static const int THETA_SL{2};
0054   /// The segment token
0055   nano_mu::EDTokenHandle<DTRecSegment4DCollection> m_token;
0056
0057   /// Fill rec-hit table
0058   bool m_fillHits;
0059
0060   /// Fill segment extrapolation  table
0061   bool m_fillExtr;
0062
0063   /// DT Geometry
0064   nano_mu::ESTokenHandle<DTGeometry, MuonGeometryRecord, edm::Transition::BeginRun> m_dtGeometry;
0065
0066   /// Tracking Geometry
0067   nano_mu::ESTokenHandle<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_trackingGeometry;
0068
0069   /// Handle DT trigger time pedestals
0070   std::unique_ptr<DTTTrigBaseSync> m_dtSync;
0071 };
0072
0073 MuDTSegmentExtTableProducer::MuDTSegmentExtTableProducer(const edm::ParameterSet& config)
0074     : MuBaseFlatTableProducer{config},
0075       m_token{config, consumesCollector(), "src"},
0076       m_fillHits{config.getParameter<bool>("fillHits")},
0077       m_fillExtr{config.getParameter<bool>("fillExtrapolation")},
0078       m_dtGeometry{consumesCollector()},
0079       m_trackingGeometry{consumesCollector()} {
0080   produces<nanoaod::FlatTable>();
0081   if (m_fillHits) {
0082     produces<nanoaod::FlatTable>("hits");
0083   }
0084   if (m_fillExtr) {
0085     produces<nanoaod::FlatTable>("extr");
0086   }
0087
0088   m_dtSync = DTTTrigSyncFactory::get()->create(config.getParameter<std::string>("tTrigMode"),
0089                                                config.getParameter<edm::ParameterSet>("tTrigModeConfig"),
0090                                                consumesCollector());
0091 }
0092
0093 void MuDTSegmentExtTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0094   edm::ParameterSetDescription desc;
0095
0096   desc.add<std::string>("name", "dtSegment");
0097   desc.add<edm::InputTag>("src", edm::InputTag{"dt4DSegments"});
0098   desc.add<bool>("fillExtrapolation", true);
0099   desc.add<bool>("fillHits", true);
0100
0101   desc.add<std::string>("tTrigMode", "DTTTrigSyncFromDB");
0102
0103   edm::ParameterSetDescription tTrigModeParams;
0104   tTrigModeParams.add<bool>("doTOFCorrection", true);
0105   tTrigModeParams.add<int>("tofCorrType", 2);
0106
0107   tTrigModeParams.add<double>("vPropWire", 24.4);
0108   tTrigModeParams.add<bool>("doWirePropCorrection", true);
0109   tTrigModeParams.add<int>("wirePropCorrType", 0);
0110
0111   tTrigModeParams.add<std::string>("tTrigLabel", "");
0112   tTrigModeParams.add<bool>("doT0Correction", true);
0113   tTrigModeParams.add<std::string>("t0Label", "");
0114   tTrigModeParams.addUntracked<bool>("debug", false);
0115
0116   desc.add<edm::ParameterSetDescription>("tTrigModeConfig", tTrigModeParams);
0117
0118   descriptions.addWithDefaultLabel(desc);
0119 }
0120
0121 void MuDTSegmentExtTableProducer::getFromES(const edm::Run& run, const edm::EventSetup& environment) {
0122   m_dtGeometry.getFromES(environment);
0123 }
0124
0125 void MuDTSegmentExtTableProducer::getFromES(const edm::EventSetup& environment) {
0126   m_trackingGeometry.getFromES(environment);
0127   m_dtSync->setES(environment);
0128 }
0129
0130 void MuDTSegmentExtTableProducer::fillTable(edm::Event& ev) {
0131   unsigned int nSegments{0};
0132
0133   std::vector<float> seg4D_posLoc_x_SL1;
0134   std::vector<float> seg4D_posLoc_x_SL3;
0135   std::vector<float> seg4D_posLoc_x_midPlane;
0136
0137   std::vector<uint32_t> seg4D_extr_begin;
0138   std::vector<uint32_t> seg4D_extr_end;
0139
0140   std::vector<uint32_t> seg2D_hits_begin;
0141   std::vector<uint32_t> seg2D_hits_end;
0142
0143   // segment extrapolation to DT layers filled, if m_fillExtr == true
0144   unsigned int nExtr{0};
0145
0146   std::vector<float> seg4D_hitsExpPos;
0147   std::vector<float> seg4D_hitsExpPosCh;
0148   std::vector<int16_t> seg4D_hitsExpWire;
0149
0150   // rec-hits vectors, filled if m_fillHits == true
0151   unsigned int nHits{0};
0152
0153   std::vector<float> seg2D_hits_pos;
0154   std::vector<float> seg2D_hits_posCh;
0155   std::vector<float> seg2D_hits_posErr;
0156   std::vector<int16_t> seg2D_hits_side;
0157   std::vector<int16_t> seg2D_hits_wire;
0158   std::vector<int16_t> seg2D_hits_wirePos;
0159   std::vector<int16_t> seg2D_hits_layer;
0160   std::vector<int16_t> seg2D_hits_superLayer;
0161   std::vector<float> seg2D_hits_time;
0162   std::vector<float> seg2D_hits_timeCali;
0163
0164   auto fillHits = [&](const DTRecSegment2D* seg, const GeomDet* chamb) {
0165     const auto& recHits = seg->specificRecHits();
0166
0167     for (const auto& recHit : recHits) {
0168       ++nHits;
0169
0170       const auto wireId = recHit.wireId();
0171       const auto layerId = wireId.layerId();
0172       const auto* layer = m_dtGeometry->layer(layerId);
0173
0174       seg2D_hits_pos.push_back(recHit.localPosition().x());
0175       seg2D_hits_posCh.push_back(chamb->toLocal(layer->toGlobal(recHit.localPosition())).x());
0176       seg2D_hits_posErr.push_back(recHit.localPositionError().xx());
0177
0178       seg2D_hits_side.push_back(recHit.lrSide());
0179       seg2D_hits_wire.push_back(wireId.wire());
0180       seg2D_hits_wirePos.push_back(layer->specificTopology().wirePosition(wireId.wire()));
0181       seg2D_hits_layer.push_back(layerId.layer());
0182       seg2D_hits_superLayer.push_back(layerId.superlayer());
0183
0184       auto digiTime = recHit.digiTime();
0185
0186       auto tTrig = m_dtSync->offset(wireId);
0187
0188       seg2D_hits_time.push_back(digiTime);
0189       seg2D_hits_timeCali.push_back(digiTime - tTrig);
0190     }
0191   };
0192
0193   auto segments4D = m_token.conditionalGet(ev);
0194
0195   if (segments4D.isValid()) {
0196     auto chambIt = segments4D->id_begin();
0197     const auto chambEnd = segments4D->id_end();
0198
0199     for (; chambIt != chambEnd; ++chambIt) {
0200       const auto& range = segments4D->get(*chambIt);
0201
0202       for (auto segment4D = range.first; segment4D != range.second; ++segment4D) {
0203         auto station = (*chambIt).station();
0204         auto wheel = (*chambIt).wheel();
0205         auto sector = (*chambIt).sector();
0206
0207         bool hasPhi = segment4D->hasPhi();
0208         bool hasZed = segment4D->hasZed();
0209
0210         auto pos = segment4D->localPosition();
0211         auto dir = segment4D->localDirection();
0212
0213         std::array<float, 2> xPosLocSL{{DEFAULT_DOUBLE_VAL, DEFAULT_DOUBLE_VAL}};
0214         std::array<bool, 2> hasPptSL{{false, false}};
0215         auto xPosLocMidPlane = DEFAULT_DOUBLE_VAL;
0216
0217         const auto* chamb = m_dtGeometry->chamber(*chambIt);
0218
0219         StraightLinePlaneCrossing segmentPlaneCrossing{
0220             chamb->toGlobal(pos).basicVector(), chamb->toGlobal(dir).basicVector(), anyDirection};
0221
0222         if (hasPhi) {
0223           for (int iSL = 0; iSL < 2; ++iSL) {
0224             const auto* sl = chamb->superLayer(1 + iSL * 2);
0225
0226             auto pptSL = segmentPlaneCrossing.position(sl->surface());
0227             hasPptSL[iSL] = pptSL.first;
0228
0229             if (hasPptSL[iSL]) {
0230               GlobalPoint segExrapolationToSL(pptSL.second);
0231               xPosLocSL[iSL] = chamb->toLocal(segExrapolationToSL).x();
0232             }
0233           }
0234         }
0235
0236         seg4D_posLoc_x_SL1.push_back(xPosLocSL[0]);
0237         seg4D_posLoc_x_SL3.push_back(xPosLocSL[1]);
0238
0239         if (hasPptSL[0] && hasPptSL[1])
0240           xPosLocMidPlane = (xPosLocSL[0] + xPosLocSL[1]) * 0.5;
0241
0242         seg4D_posLoc_x_midPlane.push_back(xPosLocMidPlane);
0243
0244         const auto begin = seg4D_hitsExpPos.size();
0245
0246         const auto size{station == 4 ? 8 : 12};
0247
0248         nExtr += size;
0249         seg4D_extr_begin.push_back(begin);
0250         seg4D_extr_end.push_back(begin + size);
0251
0252         const auto iSLs = station < 4 ? std::vector{1, 2, 3} : std::vector{1, 3};
0253
0254         for (int iL = FIRST_LAYER; iL <= LAST_LAYER; ++iL) {
0255           for (const auto iSL : iSLs) {
0256             auto* layer = m_dtGeometry->layer(DTWireId{wheel, station, sector, iSL, iL, 2});
0257             auto ppt = segmentPlaneCrossing.position(layer->surface());
0258
0259             bool success{ppt.first};  // check for failure
0260
0261             auto expPos{DEFAULT_DOUBLE_VAL};
0262             auto expPosCh{DEFAULT_DOUBLE_VAL};
0263             auto expWire{DEFAULT_INT_VAL_POS};
0264
0265             if (success) {
0266               GlobalPoint segExrapolationToLayer(ppt.second);
0267               LocalPoint segPosAtZWireLayer = layer->toLocal(segExrapolationToLayer);
0268               LocalPoint segPosAtZWireChamber = chamb->toLocal(segExrapolationToLayer);
0269
0270               if (hasPhi && iSL != THETA_SL) {
0271                 expPos = segPosAtZWireLayer.x();
0272                 expPosCh = segPosAtZWireChamber.x();
0273                 expWire = layer->specificTopology().channel(segPosAtZWireLayer);
0274               } else if (hasZed && iSL == THETA_SL) {
0275                 expPos = segPosAtZWireLayer.x();
0276                 expPosCh = segPosAtZWireChamber.y();
0277                 expWire = layer->specificTopology().channel(segPosAtZWireLayer);
0278               }
0279             }
0280
0281             seg4D_hitsExpPos.push_back(expPos);
0282             seg4D_hitsExpPosCh.push_back(expPosCh);
0283             seg4D_hitsExpWire.push_back(expWire);
0284           }
0285         }
0286
0287         seg2D_hits_begin.push_back(seg2D_hits_pos.size());
0288
0289         const GeomDet* geomDet = m_trackingGeometry->idToDet(segment4D->geographicalId());
0290         if (hasPhi) {
0291           fillHits(segment4D->phiSegment(), geomDet);
0292         }
0293
0294         if (hasZed) {
0295           fillHits(segment4D->zSegment(), geomDet);
0296         }
0297
0298         seg2D_hits_end.push_back(seg2D_hits_pos.size());
0299         ++nSegments;
0300       }
0301     }
0302   }
0303
0304   auto table = std::make_unique<nanoaod::FlatTable>(nSegments, m_name, false, true);
0305
0306   table->setDoc("DT segment information");
0307
0308   addColumn(table, "seg4D_posLoc_x_SL1", seg4D_posLoc_x_SL1, "position x at SL1 in local coordinates - cm");
0309   addColumn(table, "seg4D_posLoc_x_SL3", seg4D_posLoc_x_SL3, "position x at SL3 in local coordinates - cm");
0310   addColumn(table,
0311             "seg4D_posLoc_x_midPlane",
0312             seg4D_posLoc_x_midPlane,
0313             "position x at SL1 - SL3 mid plane in local coordinates - cm");
0314
0315   addColumn(table, "seg2D_hits_begin", seg2D_hits_begin, "begin() of range of quantities in the *_hits_* vectors");
0316   addColumn(table, "seg2D_hits_end", seg2D_hits_end, "end() of range of quantities in the *_hits_* vectors");
0317
0318   addColumn(table, "seg4D_extr_begin", seg4D_extr_begin, "begin() of range of quantities in the *_extr_* vectors");
0319   addColumn(table, "seg4D_extr_end", seg4D_extr_end, "end() of range of quantities in the *_extr_* vectors");
0320
0321   ev.put(std::move(table));
0322
0323   if (m_fillHits) {
0324     auto tabHits = std::make_unique<nanoaod::FlatTable>(nHits, m_name + "_hits", false, false);
0325
0326     tabHits->setDoc("Size of DT segment *_hits_* vectors");
0327
0328     addColumn(tabHits, "pos", seg2D_hits_pos, "local x position of a hit in layer local coordinates");
0329     addColumn(tabHits, "posCh", seg2D_hits_posCh, "local x position of a hit in chamber local coordinates");
0330     addColumn(tabHits,
0331               "posErr",
0332               seg2D_hits_posErr,
0333               "local position error of a hit in layer local coordinates - xx component of error matrix");
0334     addColumn(tabHits, "side", seg2D_hits_side, "is hit on L/R side of a cell wire - 1/2 is R/L");
0335     addColumn(tabHits, "wire", seg2D_hits_wire, "hit wire number - range depends on chamber size");
0336     addColumn(tabHits, "wirePos", seg2D_hits_wirePos, "hit wire x position in layer local coordinates");
0337     addColumn(tabHits, "layer", seg2D_hits_layer, "hit layer number - range [1:4]");
0338     addColumn(tabHits,
0339               "superLayer",
0340               seg2D_hits_superLayer,
0341               "hit superlayer - [1:3] range"
0342               "<br />SL 1 and 3 are phi SLs"
0343               "<br />SL 2 is theta SL");
0344     addColumn(tabHits, "time", seg2D_hits_time, "digi time - ns, pedestal not subtracted");
0345     addColumn(tabHits, "timeCali", seg2D_hits_timeCali, "digi time - ns, pedestal subtracted");
0346
0347     ev.put(std::move(tabHits), "hits");
0348   }
0349
0350   if (m_fillExtr) {
0351     auto tabExtr = std::make_unique<nanoaod::FlatTable>(nExtr, m_name + "_extr", false, false);
0352
0353     tabExtr->setDoc("Size of DT segment *_extr_* vectors");
0354     addColumn(tabExtr,
0355               "ExpPos",
0356               seg4D_hitsExpPos,
0357               "expected position of segment extrapolated"
0358               "<br />to a given layer in layer local coordinates - cm");
0359
0360     addColumn(tabExtr,
0361               "ExpPosCh",
0362               seg4D_hitsExpPosCh,
0363               "expected position of segment extrapolated"
0364               "<br />to a given layer in chhamber local coordinates - cm");
0365
0366     addColumn(tabExtr,
0367               "ExpWire",
0368               seg4D_hitsExpWire,
0369               "expected wire crossed by segment extrapolated"
0370               "<br />to a given layer - range depends on chamber size");
0371
0372     ev.put(std::move(tabExtr), "extr");
0373   }
0374 }
0375
0376 #include "FWCore/PluginManager/interface/ModuleDef.h"
0377 #include "FWCore/Framework/interface/MakerMacros.h"
0378
0379 DEFINE_FWK_MODULE(MuDTSegmentExtTableProducer);