CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/Validation/DTRecHits/src/DTHitQualityUtils.cc

Go to the documentation of this file.
00001 
00002 /*
00003  *  See header file for a description of this class.
00004  *
00005  *  $Date: 2008/12/03 12:17:41 $
00006  *  $Revision: 1.7 $
00007  *  \author S. Bolognesi and G. Cerminara - INFN Torino
00008  */
00009 
00010 
00011 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
00012 #include "DataFormats/DetId/interface/DetId.h"
00013 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00014 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00015 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
00016 
00017 #include <iostream>
00018 
00019 using namespace std;
00020 using namespace edm;
00021 
00022 bool DTHitQualityUtils::debug;
00023 
00024 // Constructor
00025 DTHitQualityUtils::DTHitQualityUtils(){
00026   //DTHitQualityUtils::setDebug(debug);
00027 }
00028 
00029 // Destructor
00030 DTHitQualityUtils::~DTHitQualityUtils(){
00031 }
00032 
00033 
00034 // Return a map between simhits of a layer,superlayer or chamber and the wireId of their cell
00035 map<DTWireId, PSimHitContainer >
00036 DTHitQualityUtils::mapSimHitsPerWire(const PSimHitContainer& simhits) {
00037   map<DTWireId, PSimHitContainer > hitWireMapResult;
00038    
00039   for(PSimHitContainer::const_iterator simhit = simhits.begin();
00040       simhit != simhits.end();
00041       simhit++) {
00042     hitWireMapResult[DTWireId((*simhit).detUnitId())].push_back(*simhit);
00043   }
00044    
00045   return hitWireMapResult;
00046 }
00047 
00048 // Extract mu simhits from a map of simhits by wire and map them by wire
00049 map<DTWireId, const PSimHit*>
00050 DTHitQualityUtils::mapMuSimHitsPerWire(const map<DTWireId, PSimHitContainer>& simHitWireMap) {
00051 
00052   map<DTWireId, const PSimHit*> ret;
00053   
00054   for(map<DTWireId, PSimHitContainer>::const_iterator wireAndSimHit = simHitWireMap.begin();
00055       wireAndSimHit != simHitWireMap.end();
00056       wireAndSimHit++) {
00057 
00058     const PSimHit* muHit = findMuSimHit((*wireAndSimHit).second);
00059     if(muHit != 0) {
00060       ret[(*wireAndSimHit).first]=(muHit);
00061     }
00062   }
00063   return ret;
00064 }
00065 
00066 
00067 // Find the sim hit from muon track in a vector of simhits
00068 // If no mu sim hit is found then it returns a null pointer
00069 const PSimHit* DTHitQualityUtils::findMuSimHit(const PSimHitContainer& hits) {
00070   //PSimHitContainer muHits;
00071   vector<const PSimHit*> muHits;
00072  
00073    // Loop over simhits
00074   for (PSimHitContainer::const_iterator hit=hits.begin();
00075        hit != hits.end(); hit++) {
00076     if (abs((*hit).particleType())==13) muHits.push_back(&(*hit));
00077   }
00078 
00079   if (muHits.size()==0)
00080     return 0; //FIXME: Throw of exception???
00081   else if (muHits.size()>1)
00082     if(debug)
00083       cout << "[DTHitQualityUtils]***WARNING: # muSimHits in a wire = " << muHits.size() << endl;
00084 
00085   return (muHits.front());
00086 }
00087 
00088 
00089 // Find Innermost and outermost SimHit from Mu in a SL (they identify a simulated segment)
00090 pair<const PSimHit*, const PSimHit*>
00091 DTHitQualityUtils::findMuSimSegment(const map<DTWireId, const PSimHit*>& mapWireAndMuSimHit) {
00092 
00093   int outSL = 0;
00094   int inSL = 4;
00095   int outLayer = 0;
00096   int inLayer = 5;
00097   const PSimHit *inSimHit = 0;
00098   const PSimHit *outSimHit = 0;
00099 
00100   for(map<DTWireId, const PSimHit*>::const_iterator wireAndMuSimHit = mapWireAndMuSimHit.begin();
00101       wireAndMuSimHit != mapWireAndMuSimHit.end();
00102       wireAndMuSimHit++) {
00103     
00104     const DTWireId wireId = (*wireAndMuSimHit).first;
00105     const PSimHit *theMuHit = (*wireAndMuSimHit).second;
00106 
00107     int sl = ((wireId.layerId()).superlayerId()).superLayer();
00108     int layer = (wireId.layerId()).layer();
00109 
00110     if(sl == outSL) {
00111       if(layer > outLayer) {
00112         outLayer = layer;
00113         outSimHit = theMuHit;
00114       }
00115     }
00116     if(sl > outSL) {
00117       outSL = sl;
00118       outLayer = layer;
00119       outSimHit = theMuHit;
00120     }
00121     if(sl == inSL) {
00122       if(layer < inLayer) {
00123         inLayer = layer;
00124         inSimHit = theMuHit;
00125       }
00126     }
00127     if(sl < inSL) {
00128       inSL = sl;
00129       inLayer = layer;
00130       inSimHit = theMuHit;
00131     }
00132   }
00133 
00134   if(inSimHit != 0) {
00135     if(debug)
00136       cout << "Innermost SimHit on SL: " << inSL << " layer: " << inLayer << endl;
00137   } else {
00138     cout << "[DTHitQualityUtils]***Error: No Innermost SimHit found!!!" << endl;
00139     abort();
00140   }
00141 
00142   if(outSimHit != 0) {
00143     if(debug)
00144       cout << "Outermost SimHit on SL: " << outSL << " layer: " << outLayer << endl;
00145   } else {
00146     cout << "[DTHitQualityUtils]***Error: No Outermost SimHit found!!!" << endl;
00147     abort();
00148   }
00149 
00150   // //Check that outermost and innermost SimHit are not the same
00151   // if(outSimHit == inSimHit) {
00152   //   cout << "[DTHitQualityUtils]***Warning: outermost and innermost SimHit are the same!" << endl;
00153   //   abort();
00154   //     }
00155   return make_pair(inSimHit, outSimHit);
00156 }
00157 
00158 
00159 
00160 // Find direction and position of a segment (in local RF) from outer and inner mu SimHit in the RF of object Det
00161 // (Concrete implementation of Det are MuBarSL and MuBarChamber)
00162 pair<LocalVector, LocalPoint>
00163 DTHitQualityUtils::findMuSimSegmentDirAndPos(const pair<const PSimHit*, const PSimHit*>& inAndOutSimHit,
00164                                              const DetId detId, const DTGeometry *muonGeom) {
00165 
00166   //FIXME: What should happen if outSimHit = inSimHit???? Now, this case is not considered
00167   const PSimHit* innermostMuSimHit = inAndOutSimHit.first;
00168   const PSimHit* outermostMuSimHit = inAndOutSimHit.second;
00169 
00170   //Find simulated segment direction from SimHits position
00171   const DTLayer* layerIn = muonGeom->layer((DTWireId(innermostMuSimHit->detUnitId())).layerId()); 
00172   const DTLayer* layerOut = muonGeom->layer((DTWireId(outermostMuSimHit->detUnitId())).layerId()); 
00173   GlobalPoint inGlobalPos = layerIn->toGlobal(innermostMuSimHit->localPosition());
00174   GlobalPoint outGlobalPos = layerOut->toGlobal(outermostMuSimHit->localPosition());
00175   LocalVector simHitDirection = (muonGeom->idToDet(detId))->toLocal(inGlobalPos - outGlobalPos);
00176   simHitDirection = -simHitDirection.unit();
00177 
00178   //SimHit position extrapolated at z=0 in the Det RF
00179   LocalPoint outLocalPos = (muonGeom->idToDet(detId))->toLocal(outGlobalPos);
00180   LocalPoint simSegLocalPosition =
00181     outLocalPos + simHitDirection*(-outLocalPos.z()/(simHitDirection.mag()*cos(simHitDirection.theta())));
00182 
00183   return make_pair(simHitDirection, simSegLocalPosition);
00184 }
00185 
00186 // Find the angles from a segment direction:
00187 // NB: For 4D RecHits: 
00188 //                        Alpha = angle measured by SL RPhi
00189 //                        Beta  = angle measured by SL RZ
00190 //     For 2D RecHits: only Alpha makes sense
00191 pair<double, double> DTHitQualityUtils::findSegmentAlphaAndBeta(const LocalVector& direction) {
00192   //return make_pair(atan(direction.x()/direction.z()), atan(direction.y()/direction.z()));
00193   return make_pair((direction.x()/direction.z()), (direction.y()/direction.z()));
00194 }
00195 
00196 
00197 
00198 
00199 
00200 
00201 
00202 
00203 
00204