DTRecHit/plugins/DTLinearDriftFromDBAlgo.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author S. Bolognesi - INFN Torino
0005  */
0006
0007 #include "RecoLocalMuon/DTRecHit/plugins/DTLinearDriftFromDBAlgo.h"
0008 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0009 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0010 #include "Geometry/DTGeometry/interface/DTLayer.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "CondFormats/DTObjects/interface/DTMtime.h"
0017 #include "CondFormats/DataRecord/interface/DTMtimeRcd.h"
0018 #include "CondFormats/DTObjects/interface/DTRecoConditions.h"
0019 #include "CondFormats/DataRecord/interface/DTRecoConditionsVdriftRcd.h"
0020 #include "CondFormats/DataRecord/interface/DTRecoConditionsUncertRcd.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "MagneticField/Engine/interface/MagneticField.h"
0023 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0024
0025 using namespace std;
0026 using namespace edm;
0027
0028 DTLinearDriftFromDBAlgo::DTLinearDriftFromDBAlgo(const ParameterSet& config, ConsumesCollector cc)
0029     : DTRecHitBaseAlgo(config, cc),
0030       mTimeMap(nullptr),
0031       vDriftMap(nullptr),
0032       field(nullptr),
0033       fieldToken_(cc.esConsumes()),
0034       nominalB(-1),
0035       minTime(config.getParameter<double>("minTime")),
0036       maxTime(config.getParameter<double>("maxTime")),
0037       doVdriftCorr(config.getParameter<bool>("doVdriftCorr")),
0038       // Option to force going back to digi time at Step 2
0039       stepTwoFromDigi(config.getParameter<bool>("stepTwoFromDigi")),
0040       useUncertDB(config.getParameter<bool>("useUncertDB")),
0041       readLegacyTTrigDB(config.getParameter<bool>("readLegacyTTrigDB")),
0042       readLegacyVDriftDB(config.getParameter<bool>("readLegacyVDriftDB")),
0043       // Set verbose output
0044       debug(config.getUntrackedParameter<bool>("debug")) {
0045   if (readLegacyVDriftDB) {
0046     mTimeMapToken_ = cc.esConsumes();
0047   } else {
0048     vDriftMapToken_ = cc.esConsumes();
0049   }
0050   if (useUncertDB) {
0051     uncertMapToken_ = cc.esConsumes();
0052   }
0053 }
0054
0055 DTLinearDriftFromDBAlgo::~DTLinearDriftFromDBAlgo() {}
0056
0057 void DTLinearDriftFromDBAlgo::setES(const EventSetup& setup) {
0058   if (debug)
0059     edm::LogVerbatim("DTLocalReco") << "[DTLinearDriftFromDBAlgo] setES called" << endl;
0060   theSync->setES(setup);
0061   // Get the map of ttrig from the Setup
0062   if (readLegacyVDriftDB) {
0063     mTimeMap = &setup.getData(mTimeMapToken_);
0064     vDriftMap = nullptr;
0065   } else {
0066     vDriftMap = &setup.getData(vDriftMapToken_);
0067     mTimeMap = nullptr;
0068
0069     // Consistency check: no parametrization is implemented for the time being
0070     int version = vDriftMap->version();
0071     if (version != 1) {
0072       throw cms::Exception("Configuration") << "only version 1 is presently supported for VDriftDB";
0073     }
0074   }
0075
0076   field = &setup.getData(fieldToken_);
0077   nominalB = field->nominalValue();
0078
0079   if (useUncertDB) {
0080     uncertMap = &setup.getData(uncertMapToken_);
0081     if (uncertMap->version() > 1)
0082       edm::LogError("NotImplemented") << "DT Uncertainty DB version unsupported: " << uncertMap->version();
0083   }
0084
0085   if (debug) {
0086     if (readLegacyVDriftDB) {
0087       edm::LogVerbatim("DTLocalReco") << "[DTLinearDriftFromDBAlgo] meanTimer version: " << mTimeMap->version() << endl;
0088     } else {
0089       edm::LogVerbatim("DTLocalReco") << "[DTLinearDriftFromDBAlgo] vDrift version: " << vDriftMap->version() << endl;
0090     }
0091
0092     if (useUncertDB)
0093       edm::LogVerbatim("DTLocalReco") << "                          uncertDB  version: " << uncertMap->version()
0094                                       << endl;
0095   }
0096 }
0097
0098 // First Step
0099 bool DTLinearDriftFromDBAlgo::compute(
0100     const DTLayer* layer, const DTDigi& digi, LocalPoint& leftPoint, LocalPoint& rightPoint, LocalError& error) const {
0101   // Get the wireId
0102   DTLayerId layerId = layer->id();
0103   const DTWireId wireId(layerId, digi.wire());
0104
0105   // Get Wire position
0106   if (!layer->specificTopology().isWireValid(digi.wire()))
0107     return false;
0108   LocalPoint locWirePos(layer->specificTopology().wirePosition(digi.wire()), 0, 0);
0109   const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
0110
0111   return compute(layer, wireId, digi.time(), globWirePos, leftPoint, rightPoint, error, 1);
0112 }
0113
0114 // Second step: the same as 1st step (optionally, redo 1st step starting from digi time)
0115 bool DTLinearDriftFromDBAlgo::compute(const DTLayer* layer,
0116                                       const DTRecHit1D& recHit1D,
0117                                       const float& angle,
0118                                       DTRecHit1D& newHit1D) const {
0119   if (!stepTwoFromDigi) {
0120     newHit1D.setPositionAndError(recHit1D.localPosition(), recHit1D.localPositionError());
0121     return true;
0122   }
0123
0124   const DTWireId wireId = recHit1D.wireId();
0125
0126   // Get Wire position
0127   if (!layer->specificTopology().isWireValid(wireId.wire()))
0128     return false;
0129   LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
0130   const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
0131
0132   return compute(layer, wireId, recHit1D.digiTime(), globWirePos, newHit1D, 2);
0133 }
0134
0135 // Third step.
0136 bool DTLinearDriftFromDBAlgo::compute(const DTLayer* layer,
0137                                       const DTRecHit1D& recHit1D,
0138                                       const float& angle,
0139                                       const GlobalPoint& globPos,
0140                                       DTRecHit1D& newHit1D) const {
0141   return compute(layer, recHit1D.wireId(), recHit1D.digiTime(), globPos, newHit1D, 3);
0142 }
0143
0144 // Do the actual work.
0145 bool DTLinearDriftFromDBAlgo::compute(const DTLayer* layer,
0146                                       const DTWireId& wireId,
0147                                       const float digiTime,
0148                                       const GlobalPoint& globPos,
0149                                       LocalPoint& leftPoint,
0150                                       LocalPoint& rightPoint,
0151                                       LocalError& error,
0152                                       int step) const {
0153   // Subtract the offset to the digi time accordingly to the DTTTrigBaseSync concrete instance
0154   float driftTime = digiTime - theSync->offset(layer, wireId, globPos);
0155
0156   // check for out-of-time
0157   if (driftTime < minTime || driftTime > maxTime) {
0158     if (debug)
0159       edm::LogWarning("DTLocalReco") << "[DTLinearDriftFromDBAlgo]*** Drift time out of window for in-time hits "
0160                                      << driftTime << endl;
0161
0162     if (step == 1) {  //FIXME: protection against failure at 2nd and 3rd steps, must be checked!!!
0163       // Hits are interpreted as coming from out-of-time pile-up and recHit
0164       // is ignored.
0165       return false;
0166     }
0167   }
0168
0169   // Small negative times interpreted as hits close to the wire.
0170   if (driftTime < 0.)
0171     driftTime = 0;
0172
0173   // Read the vDrift and reso for this wire
0174   float vDrift = 0;
0175   float hitResolution = 0;
0176
0177   if (readLegacyVDriftDB) {
0178     // vdrift is cm/ns , resolution is cm
0179     mTimeMap->get(wireId.superlayerId(),
0180                   vDrift,
0181                   hitResolution,  // Value from vdrift DB; replaced below if useUncertDB card is set
0182                   DTVelocityUnits::cm_per_ns);
0183   } else {
0184     // For v2, we will pass also: double args[1] = {(layer->toLocal(globPos)).y()};
0185     vDrift = vDriftMap->get(wireId);
0186   }
0187
0188   if (useUncertDB) {
0189     // Read the uncertainty from the DB for the given channel and step
0190     double args[1] = {double(step - 1)};
0191     hitResolution = uncertMap->get(wireId, args);
0192   }
0193
0194   //only in step 3
0195   if (doVdriftCorr && step == 3 && nominalB != 0) {
0196     if (abs(wireId.wheel()) == 2 && wireId.station() == 1 && wireId.superLayer() != 2) {
0197       // Variation of vdrift along Y due to B field,
0198       /// cf. http://arxiv.org/PS_cache/arxiv/pdf/0911/0911.4895v2.pdf
0199       // vdrift is lower  a negative Y (lower global |Z|)
0200       const float k_param = 1.2e-04;
0201       LocalPoint local_pos = layer->toLocal(globPos);
0202       vDrift = vDrift * (1. - k_param * local_pos.y());
0203     }
0204   }
0205
0206   // Compute the drift distance
0207   float drift = driftTime * vDrift;
0208
0209   // Get Wire position
0210   if (!layer->specificTopology().isWireValid(wireId.wire()))
0211     return false;
0212   LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
0213   //Build the two possible points and the error on the position
0214   leftPoint = LocalPoint(locWirePos.x() - drift, locWirePos.y(), locWirePos.z());
0215   rightPoint = LocalPoint(locWirePos.x() + drift, locWirePos.y(), locWirePos.z());
0216   error = LocalError(hitResolution * hitResolution, 0., 0.);
0217
0218   if (debug) {
0219     edm::LogWarning("DTLocalReco") << "[DTLinearDriftFromDBAlgo] Compute drift distance, for digi at wire: " << wireId
0220                                    << endl
0221                                    << "       Step:           " << step << endl
0222                                    << "       Digi time:      " << digiTime << endl
0223                                    << "       Drift time:     " << driftTime << endl
0224                                    << "       Drift distance: " << drift << endl
0225                                    << "       Hit Resolution: " << hitResolution << endl
0226                                    << "       Left point:     " << leftPoint << endl
0227                                    << "       Right point:    " << rightPoint << endl
0228                                    << "       Error:          " << error << endl;
0229   }
0230
0231   return true;
0232 }
0233
0234 // Interface to the method which does the actual work suited for 2nd and 3rd steps
0235 bool DTLinearDriftFromDBAlgo::compute(const DTLayer* layer,
0236                                       const DTWireId& wireId,
0237                                       const float digiTime,
0238                                       const GlobalPoint& globPos,
0239                                       DTRecHit1D& newHit1D,
0240                                       int step) const {
0241   LocalPoint leftPoint;
0242   LocalPoint rightPoint;
0243   LocalError error;
0244
0245   if (compute(layer, wireId, digiTime, globPos, leftPoint, rightPoint, error, step)) {
0246     // Set the position and the error of the rechit which is being updated
0247     switch (newHit1D.lrSide()) {
0248       case DTEnums::Left: {
0249         // Keep the original y position of newHit1D: for step==3, it's the
0250         // position along the wire. Needed for rotation alignment
0251         LocalPoint leftPoint3D(leftPoint.x(), newHit1D.localPosition().y(), leftPoint.z());
0252         newHit1D.setPositionAndError(leftPoint3D, error);
0253         break;
0254       }
0255
0256       case DTEnums::Right: {
0257         // as above: 3d position
0258         LocalPoint rightPoint3D(rightPoint.x(), newHit1D.localPosition().y(), rightPoint.z());
0259         newHit1D.setPositionAndError(rightPoint3D, error);
0260         break;
0261       }
0262
0263       default:
0264         throw cms::Exception("InvalidDTCellSide") << "[DTLinearDriftFromDBAlgo] Compute at Step " << step
0265                                                   << ", Hit side " << newHit1D.lrSide() << " is invalid!" << endl;
0266         return false;
0267     }
0268
0269     return true;
0270   } else {
0271     return false;
0272   }
0273 }