DTRecHits/src/DTHitQualityUtils.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author S. Bolognesi and G. Cerminara - INFN Torino
0006  */
0007
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0011 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0012 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0013
0014 #include <iostream>
0015
0016 using namespace std;
0017 using namespace edm;
0018
0019 std::atomic<bool> DTHitQualityUtils::debug{false};
0020
0021 // Return a map between simhits of a layer,superlayer or chamber and the wireId
0022 // of their cell
0023 map<DTWireId, PSimHitContainer> DTHitQualityUtils::mapSimHitsPerWire(const PSimHitContainer &simhits) {
0024   map<DTWireId, PSimHitContainer> hitWireMapResult;
0025
0026   for (PSimHitContainer::const_iterator simhit = simhits.begin(); simhit != simhits.end(); simhit++) {
0027     hitWireMapResult[DTWireId((*simhit).detUnitId())].push_back(*simhit);
0028   }
0029
0030   return hitWireMapResult;
0031 }
0032
0033 // Extract mu simhits from a map of simhits by wire and map them by wire
0034 map<DTWireId, const PSimHit *> DTHitQualityUtils::mapMuSimHitsPerWire(
0035     const map<DTWireId, PSimHitContainer> &simHitWireMap) {
0036   map<DTWireId, const PSimHit *> ret;
0037
0038   for (map<DTWireId, PSimHitContainer>::const_iterator wireAndSimHit = simHitWireMap.begin();
0039        wireAndSimHit != simHitWireMap.end();
0040        wireAndSimHit++) {
0041     const PSimHit *muHit = findMuSimHit((*wireAndSimHit).second);
0042     if (muHit != nullptr) {
0043       ret[(*wireAndSimHit).first] = (muHit);
0044     }
0045   }
0046   return ret;
0047 }
0048
0049 // Find the sim hit from muon track in a vector of simhits
0050 // If no mu sim hit is found then it returns a null pointer
0051 const PSimHit *DTHitQualityUtils::findMuSimHit(const PSimHitContainer &hits) {
0052   // PSimHitContainer muHits;
0053   vector<const PSimHit *> muHits;
0054
0055   // Loop over simhits
0056   for (PSimHitContainer::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0057     if (abs((*hit).particleType()) == 13)
0058       muHits.push_back(&(*hit));
0059   }
0060
0061   if (muHits.empty())
0062     return nullptr;  // FIXME: Throw of exception???
0063   else if (muHits.size() > 1)
0064     if (debug)
0065       cout << "[DTHitQualityUtils]***WARNING: # muSimHits in a wire = " << muHits.size() << endl;
0066
0067   return (muHits.front());
0068 }
0069
0070 // Find Innermost and outermost SimHit from Mu in a SL (they identify a
0071 // simulated segment)
0072 pair<const PSimHit *, const PSimHit *> DTHitQualityUtils::findMuSimSegment(
0073     const map<DTWireId, const PSimHit *> &mapWireAndMuSimHit) {
0074   int outSL = 0;
0075   int inSL = 4;
0076   int outLayer = 0;
0077   int inLayer = 5;
0078   const PSimHit *inSimHit = nullptr;
0079   const PSimHit *outSimHit = nullptr;
0080
0081   for (map<DTWireId, const PSimHit *>::const_iterator wireAndMuSimHit = mapWireAndMuSimHit.begin();
0082        wireAndMuSimHit != mapWireAndMuSimHit.end();
0083        wireAndMuSimHit++) {
0084     const DTWireId wireId = (*wireAndMuSimHit).first;
0085     const PSimHit *theMuHit = (*wireAndMuSimHit).second;
0086
0087     int sl = ((wireId.layerId()).superlayerId()).superLayer();
0088     int layer = (wireId.layerId()).layer();
0089
0090     if (sl == outSL) {
0091       if (layer > outLayer) {
0092         outLayer = layer;
0093         outSimHit = theMuHit;
0094       }
0095     }
0096     if (sl > outSL) {
0097       outSL = sl;
0098       outLayer = layer;
0099       outSimHit = theMuHit;
0100     }
0101     if (sl == inSL) {
0102       if (layer < inLayer) {
0103         inLayer = layer;
0104         inSimHit = theMuHit;
0105       }
0106     }
0107     if (sl < inSL) {
0108       inSL = sl;
0109       inLayer = layer;
0110       inSimHit = theMuHit;
0111     }
0112   }
0113
0114   if (inSimHit != nullptr) {
0115     if (debug)
0116       cout << "Innermost SimHit on SL: " << inSL << " layer: " << inLayer << endl;
0117   } else {
0118     cout << "[DTHitQualityUtils]***Error: No Innermost SimHit found!!!" << endl;
0119     abort();
0120   }
0121
0122   if (outSimHit != nullptr) {
0123     if (debug)
0124       cout << "Outermost SimHit on SL: " << outSL << " layer: " << outLayer << endl;
0125   } else {
0126     cout << "[DTHitQualityUtils]***Error: No Outermost SimHit found!!!" << endl;
0127     abort();
0128   }
0129
0130   // //Check that outermost and innermost SimHit are not the same
0131   // if(outSimHit == inSimHit) {
0132   //   cout << "[DTHitQualityUtils]***Warning: outermost and innermost SimHit
0133   //   are the same!" << endl; abort();
0134   //     }
0135   return make_pair(inSimHit, outSimHit);
0136 }
0137
0138 // Find direction and position of a segment (in local RF) from outer and inner
0139 // mu SimHit in the RF of object Det (Concrete implementation of Det are MuBarSL
0140 // and MuBarChamber)
0141 pair<LocalVector, LocalPoint> DTHitQualityUtils::findMuSimSegmentDirAndPos(
0142     const pair<const PSimHit *, const PSimHit *> &inAndOutSimHit, const DetId detId, const DTGeometry &muonGeom) {
0143   // FIXME: What should happen if outSimHit = inSimHit???? Now, this case is not
0144   // considered
0145   const PSimHit *innermostMuSimHit = inAndOutSimHit.first;
0146   const PSimHit *outermostMuSimHit = inAndOutSimHit.second;
0147
0148   // Find simulated segment direction from SimHits position
0149   const DTLayer *layerIn = muonGeom.layer((DTWireId(innermostMuSimHit->detUnitId())).layerId());
0150   const DTLayer *layerOut = muonGeom.layer((DTWireId(outermostMuSimHit->detUnitId())).layerId());
0151   GlobalPoint inGlobalPos = layerIn->toGlobal(innermostMuSimHit->localPosition());
0152   GlobalPoint outGlobalPos = layerOut->toGlobal(outermostMuSimHit->localPosition());
0153   LocalVector simHitDirection = (muonGeom.idToDet(detId))->toLocal(inGlobalPos - outGlobalPos);
0154   simHitDirection = -simHitDirection.unit();
0155
0156   // SimHit position extrapolated at z=0 in the Det RF
0157   LocalPoint outLocalPos = (muonGeom.idToDet(detId))->toLocal(outGlobalPos);
0158   LocalPoint simSegLocalPosition =
0159       outLocalPos + simHitDirection * (-outLocalPos.z() / (simHitDirection.mag() * cos(simHitDirection.theta())));
0160
0161   return make_pair(simHitDirection, simSegLocalPosition);
0162 }
0163
0164 // Find the angles from a segment direction:
0165 // atan(dx/dz) = "phi"   angle in the chamber RF
0166 // atan(dy/dz) = "theta" angle in the chamber RF (note: this has opposite sign
0167 // in the SLZ RF!)
0168 pair<double, double> DTHitQualityUtils::findSegmentAlphaAndBeta(const LocalVector &direction) {
0169   return make_pair(atan(direction.x() / direction.z()), atan(direction.y() / direction.z()));
0170 }
0171
0172 // Find error on angle (squared) from localDirectionError, which is the error on
0173 // tan(Angle)
0174 double DTHitQualityUtils::sigmaAngle(double Angle, double sigma2TanAngle) {
0175   double XdivZ = tan(Angle);
0176   double sigma2Angle = 1 / (1 + XdivZ * XdivZ);
0177   sigma2Angle *= sigma2Angle * sigma2TanAngle;
0178
0179   return sigma2Angle;
0180 }