DTRecHit/plugins/DTParametrizedDriftAlgo.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - INFN Torino
0005  */
0006
0007 #include "RecoLocalMuon/DTRecHit/plugins/DTParametrizedDriftAlgo.h"
0008 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0009 #include "RecoLocalMuon/DTRecHit/plugins/DTTime2DriftParametrization.h"
0010
0011 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0012 #include "Geometry/DTGeometry/interface/DTLayer.h"
0013
0014 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0015 #include "MagneticField/Engine/interface/MagneticField.h"
0016
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/ConsumesCollector.h"
0021
0022 #include <iostream>
0023 #include <iomanip>
0024
0025 using namespace std;
0026 using namespace edm;
0027
0028 DTParametrizedDriftAlgo::DTParametrizedDriftAlgo(const ParameterSet& config, ConsumesCollector cc)
0029     : DTRecHitBaseAlgo(config, cc),
0030       interpolate(config.getParameter<bool>("interpolate")),
0031       minTime(config.getParameter<double>("minTime")),  // FIXME: Default was -3 ns
0032       maxTime(config.getParameter<double>("maxTime")),  // FIXME: Default was 415 ns
0033       // Set verbose output
0034       debug(config.getUntrackedParameter<bool>("debug", "false")),
0035       magField(nullptr),
0036       magFieldToken_(cc.esConsumes()) {}
0037
0038 DTParametrizedDriftAlgo::~DTParametrizedDriftAlgo() {}
0039
0040 void DTParametrizedDriftAlgo::setES(const EventSetup& setup) {
0041   theSync->setES(setup);
0042   // Access the magnetic field
0043   magField = &setup.getData(magFieldToken_);
0044 }
0045
0046 // First Step
0047 bool DTParametrizedDriftAlgo::compute(
0048     const DTLayer* layer, const DTDigi& digi, LocalPoint& leftPoint, LocalPoint& rightPoint, LocalError& error) const {
0049   // Get the layerId
0050   DTLayerId layerId = layer->id();  //FIXME: pass it instead of get it from layer
0051   const DTWireId wireId(layerId, digi.wire());
0052
0053   // Get Wire position
0054   if (!layer->specificTopology().isWireValid(wireId.wire()))
0055     return false;
0056   LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
0057   const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
0058
0059   // impact angle on DT chamber
0060   // compute the impact angle using the centre of the wire,
0061   // only for RZ superlayer (eta). In the other cases theta is very close to 0.
0062   float angle = 0.0;
0063
0064   if (layerId.superlayer() == 2) {
0065     //GlobalPoint lPos=layer->position();
0066     GlobalVector lDir = (GlobalPoint() - globWirePos).unit();
0067     LocalVector lDirLoc = layer->toLocal(lDir);
0068
0069     angle = atan(lDirLoc.x() / -lDirLoc.z());
0070   }
0071
0072   return compute(layer, wireId, digi.time(), angle, globWirePos, leftPoint, rightPoint, error, 1);
0073 }
0074
0075 // Second step
0076 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
0077                                       const DTRecHit1D& recHit1D,
0078                                       const float& angle,
0079                                       DTRecHit1D& newHit1D) const {
0080   const DTWireId wireId = recHit1D.wireId();
0081
0082   // Get Wire position
0083   if (!layer->specificTopology().isWireValid(wireId.wire()))
0084     return false;
0085   LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
0086   const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
0087
0088   return compute(layer, wireId, recHit1D.digiTime(), angle, globWirePos, newHit1D, 2);
0089 }
0090
0091 // Third step
0092 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
0093                                       const DTRecHit1D& recHit1D,
0094                                       const float& angle,
0095                                       const GlobalPoint& globPos,
0096                                       DTRecHit1D& newHit1D) const {
0097   return compute(layer, recHit1D.wireId(), recHit1D.digiTime(), angle, globPos, newHit1D, 3);
0098 }
0099
0100 // Do the actual work.
0101 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
0102                                       const DTWireId& wireId,
0103                                       const float digiTime,
0104                                       const float& angle,
0105                                       const GlobalPoint& globPos,
0106                                       LocalPoint& leftPoint,
0107                                       LocalPoint& rightPoint,
0108                                       LocalError& error,
0109                                       int step) const {
0110   // Subtract Offset according to DTTTrigBaseSync concrete instance
0111   // chosen with the 'tZeroMode' parameter
0112   float driftTime = digiTime - theSync->offset(layer, wireId, globPos);
0113
0114   // check for out-of-time only at step 1
0115   if (step == 1 && (driftTime < minTime || driftTime > maxTime)) {
0116     if (debug)
0117       cout << "*** Drift time out of window for in-time hits " << driftTime << endl;
0118
0119     if (step == 1) {  //FIXME: protection against failure at 2nd and 3rd steps, must be checked!!!
0120       // Hits are interpreted as coming from out-of-time pile-up and recHit
0121       // is ignored.
0122       return false;
0123     }
0124   }
0125
0126   // Small negative times interpreted as hits close to the wire.
0127   if (driftTime < 0.)
0128     driftTime = 0;
0129
0130   //----------------------------------------------------------------------
0131   // Magnetic Field in layer frame
0132   const LocalVector BLoc = layer->toLocal(magField->inTesla(globPos));
0133
0134   float By = BLoc.y();
0135   float Bz = BLoc.z();
0136
0137   //--------------------------------------------------------------------
0138   // Calculate the drift distance and the resolution from the parametrization
0139
0140   DTTime2DriftParametrization::drift_distance DX;
0141   static const DTTime2DriftParametrization par;
0142
0143   bool parStatus = par.computeDriftDistance_mean(driftTime, angle, By, Bz, interpolate, &DX);
0144
0145   if (!parStatus) {
0146     if (debug)
0147       cout << "[DTParametrizedDriftAlgo]*** WARNING: call to parametrization failed" << endl;
0148     return false;
0149   }
0150
0151   float sigma_l = DX.x_width_m;
0152   float sigma_r = DX.x_width_p;
0153   float drift = DX.x_drift;
0154
0155   float reso = 0;
0156
0157   bool variableSigma = false;
0158   // By defualt the errors are obtained from a fit of the residuals in the various
0159   // stations/SL.
0160   // The error returned by DTTime2DriftParametrization can not be used easily
0161   // to determine the hit error due to the way the parametrization of the error
0162   // is done (please contcat the authors for details).
0163   // Anyhow changing variableSigma==true, an attempt is done to set a variable error
0164   // according to the sigma calculated by DTTime2DriftParametrization.
0165   // Additionally, contributions from uncertaionties in TOF and signal propagation
0166   // corrections are added.
0167   // Uncertainty in the determination of incident angle and hit position
0168   // (ie B value) are NOT accounted.
0169   // This is not the default since it does not give good results...
0170
0171   int wheel = abs(wireId.wheel());
0172   if (variableSigma) {
0173     float sTDelays = 0;
0174     if (step == 1) {                    // 1. step
0175       reso = (sigma_l + sigma_r) / 2.;  // FIXME: theta/B are not yet known...
0176       if (wireId.superlayer() == 2) {   // RZ
0177         sTDelays = 2.92;
0178       } else {  // RPhi
0179         if (wheel == 0) {
0180           sTDelays = 2.74;
0181         } else if (wheel == 1) {
0182           sTDelays = 1.83;
0183         } else if (wheel == 2) {
0184           sTDelays = 1.25;
0185         }
0186       }
0187     } else if (step == 2) {             // 2. step
0188       reso = (sigma_l + sigma_r) / 2.;  // FIXME: B is not yet known...
0189       if (wireId.superlayer() == 2) {   // RZ
0190         sTDelays = 0.096;
0191       } else {  // RPhi
0192         if (wheel == 0) {
0193           sTDelays = 0.022;
0194         } else if (wheel == 1) {
0195           sTDelays = 0.079;
0196         } else if (wheel == 2) {
0197           sTDelays = 0.124;
0198         }
0199       }
0200     } else if (step == 3) {  // 3. step
0201       reso = (sigma_l + sigma_r) / 2.;
0202       if (wireId.superlayer() == 2) {  // RZ
0203         sTDelays = 0.096;
0204       } else {  // RPhi
0205         if (wheel == 0) {
0206           sTDelays = 0.022;
0207         } else if (wheel == 1) {
0208           sTDelays = 0.079;
0209         } else if (wheel == 2) {
0210           sTDelays = 0.124;
0211         }
0212       }
0213     }
0214     float sXDelays = sTDelays * DX.v_drift / 10.;
0215     reso = sqrt(reso * reso + sXDelays * sXDelays);
0216   } else {            // Use a fixed sigma, from fit of residuals.
0217     if (step == 1) {  // 1. step
0218       if (wireId.superlayer() == 2) {
0219         if (wheel == 0) {
0220           reso = 0.0250;
0221         } else if (wheel == 1) {
0222           reso = 0.0271;
0223         } else if (wheel == 2) {
0224           reso = 0.0308;
0225         }
0226       } else {
0227         reso = 0.0237;
0228       }
0229     } else if (step == 2) {  // 2. step //FIXME
0230       if (wireId.superlayer() == 2) {
0231         if (wheel == 0) {
0232           reso = 0.0250;
0233         } else if (wheel == 1) {
0234           reso = 0.0271;
0235         } else if (wheel == 2) {
0236           reso = 0.0305;
0237         }
0238       } else {
0239         reso = 0.0231;
0240       }
0241     } else if (step == 3) {  // 3. step
0242       if (wireId.superlayer() == 2) {
0243         if (wheel == 0) {
0244           reso = 0.0196;
0245         } else if (wheel == 1) {
0246           reso = 0.0210;
0247         } else if (wheel == 2) {
0248           reso = 0.0228;
0249         }
0250       } else {
0251         reso = 0.0207;
0252       }
0253     }
0254   }
0255   //--------------------------------------------------------------------
0256
0257   error = LocalError(reso * reso, 0., 0.);
0258
0259   // Get Wire position
0260   if (!layer->specificTopology().isWireValid(wireId.wire()))
0261     return false;
0262   LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
0263
0264   //Build the two possible points and the error on the position
0265   leftPoint = LocalPoint(locWirePos.x() - drift, locWirePos.y(), locWirePos.z());
0266   rightPoint = LocalPoint(locWirePos.x() + drift, locWirePos.y(), locWirePos.z());
0267
0268   if (debug) {
0269     int prevW = cout.width();
0270     cout << setiosflags(ios::showpoint | ios::fixed) << setw(3) << "[DTParametrizedDriftAlgo]: step " << step << endl
0271          << "  Global Position  " << globPos << endl
0272          << "  Local Position   " << layer->toLocal(globPos)
0273          << endl
0274          //          << "  y along Wire     " << wireCoord << endl
0275          << "  Digi time        " << digiTime
0276          << endl
0277          //          << "  dpropDelay       " << propDelay << endl
0278          //          << "  deltaTOF         " << deltaTOF << endl
0279          << " >Drif time        " << driftTime << endl
0280          << " >Angle            " << angle * 180. / M_PI << endl
0281          << " <Drift distance   " << drift << endl
0282          << " <sigma_l          " << sigma_l << endl
0283          << " <sigma_r          " << sigma_r << endl
0284          << "  reso             " << reso << endl
0285          << "  leftPoint        " << leftPoint << endl
0286          << "  rightPoint       " << rightPoint << endl
0287          << "  error            " << error << resetiosflags(ios::showpoint | ios::fixed) << setw(prevW) << endl
0288          << endl;
0289   }
0290
0291   return true;
0292 }
0293
0294 // Interface to the method which does the actual work suited for 2nd and 3rd steps
0295 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
0296                                       const DTWireId& wireId,
0297                                       const float digiTime,
0298                                       const float& angle,
0299                                       const GlobalPoint& globPos,
0300                                       DTRecHit1D& newHit1D,
0301                                       int step) const {
0302   LocalPoint leftPoint;
0303   LocalPoint rightPoint;
0304   LocalError error;
0305
0306   if (compute(layer, wireId, digiTime, angle, globPos, leftPoint, rightPoint, error, step)) {
0307     // Set the position and the error of the rechit which is being updated
0308     switch (newHit1D.lrSide()) {
0309       case DTEnums::Left: {
0310         // Keep the original y position of newHit1D: for step==3, it's the
0311         // position along the wire. Needed for rotation alignment
0312         LocalPoint leftPoint3D(leftPoint.x(), newHit1D.localPosition().y(), leftPoint.z());
0313         newHit1D.setPositionAndError(leftPoint3D, error);
0314         break;
0315       }
0316
0317       case DTEnums::Right: {
0318         // as above: 3d position
0319         LocalPoint rightPoint3D(rightPoint.x(), newHit1D.localPosition().y(), rightPoint.z());
0320         newHit1D.setPositionAndError(rightPoint3D, error);
0321         break;
0322       }
0323
0324       default:
0325         throw cms::Exception("InvalidDTCellSide") << "[DTParametrizedDriftAlgo] Compute at Step " << step
0326                                                   << ", Hit side " << newHit1D.lrSide() << " is invalid!" << endl;
0327         return false;
0328     }
0329
0330     return true;
0331   } else {
0332     return false;
0333   }
0334 }