CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/Alignment/OfflineValidation/src/TrackerValidationVariables.cc

Go to the documentation of this file.
00001 
00002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00003 #include "FWCore/Framework/interface/Frameworkfwd.h"
00004 #include "FWCore/Framework/interface/EDAnalyzer.h"
00005 #include "FWCore/Framework/interface/MakerMacros.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 
00010 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00011 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00012 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00013 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00014 #include "Geometry/CommonTopologies/interface/RadialStripTopology.h"
00015 
00016 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00017 #include "MagneticField/Engine/interface/MagneticField.h"
00018 
00019 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00020 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00021 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00022 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00023 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00025 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00026 
00027 #include "CondFormats/Alignment/interface/Definitions.h"
00028 
00029 #include "DataFormats/Math/interface/deltaPhi.h"
00030 #include "DataFormats/TrackReco/interface/Track.h"
00031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00033 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00034 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
00035 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
00036 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
00037 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
00038 
00039 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
00040 
00041 #include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
00042 
00043 #include "TMath.h"
00044 
00045 TrackerValidationVariables::TrackerValidationVariables()
00046 {
00047 
00048 }
00049 
00050 TrackerValidationVariables::TrackerValidationVariables(const edm::EventSetup& es, const edm::ParameterSet& iSetup) 
00051   : conf_(iSetup)
00052 {
00053   es.get<TrackerDigiGeometryRecord>().get(tkGeom_);
00054   es.get<IdealMagneticFieldRecord>().get(magneticField_);
00055 }
00056 
00057 TrackerValidationVariables::~TrackerValidationVariables()
00058 {
00059 
00060 }
00061 
00062 void
00063 TrackerValidationVariables::fillHitQuantities(const Trajectory* trajectory, std::vector<AVHitStruct> & v_avhitout)
00064 {
00065   TrajectoryStateCombiner tsoscomb;
00066 
00067   const std::vector<TrajectoryMeasurement> &tmColl = trajectory->measurements();
00068   for(std::vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin();
00069       itTraj != tmColl.end();
00070       ++itTraj) {
00071     
00072     if (!itTraj->updatedState().isValid()) continue;
00073     
00074     TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() );
00075     TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
00076     
00077     if(!hit->isValid() || hit->geographicalId().det() != DetId::Tracker) continue;
00078     
00079     AVHitStruct hitStruct;
00080     const DetId& hit_detId = hit->geographicalId();
00081     unsigned int IntRawDetID = (hit_detId.rawId());     
00082     unsigned int IntSubDetID = (hit_detId.subdetId());
00083     
00084     if(IntSubDetID == 0) continue;
00085     
00086     //first calculate residuals in cartesian coordinates in the local module coordinate system
00087     
00088     LocalPoint lPHit = hit->localPosition();
00089     LocalPoint lPTrk = tsos.localPosition();
00090     LocalVector lVTrk = tsos.localDirection();
00091 
00092     hitStruct.localAlpha = atan2(lVTrk.x(), lVTrk.z()); // wrt. normal tg(alpha)=x/z
00093     hitStruct.localBeta  = atan2(lVTrk.y(), lVTrk.z()); // wrt. normal tg(beta)= y/z
00094 
00095     //LocalError errHit = hit->localPositionError();
00096     // adding APE to hitError
00097     AlgebraicROOTObject<2>::SymMatrix mat = asSMatrix<2>(hit->parametersError());
00098     LocalError errHit = LocalError( mat(0,0),mat(0,1),mat(1,1) );
00099     LocalError errTrk = tsos.localError().positionError();
00100     
00101     //check for negative error values: track error can have negative value, if matrix inversion fails (very rare case)
00102     //hit error should always give positive values
00103     if(errHit.xx()<0. || errHit.yy()<0. || errTrk.xx()<0. || errTrk.yy()<0.){
00104       edm::LogError("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
00105                                                   << "One of the squared error methods gives negative result"
00106                                                   << "\n\terrHit.xx()\terrHit.yy()\terrTrk.xx()\terrTrk.yy()"
00107                                                   << "\n\t" << errHit.xx()
00108                                                   << "\t" << errHit.yy()
00109                                                   << "\t" << errTrk.xx()
00110                                                   << "\t" << errTrk.yy();
00111       continue;
00112     }
00113     
00114     align::LocalVector res = lPTrk - lPHit;
00115     
00116     float resXErr = std::sqrt( errHit.xx() + errTrk.xx() );
00117     float resYErr = std::sqrt( errHit.yy() + errTrk.yy() );
00118     
00119     hitStruct.resX = res.x();
00120     hitStruct.resY = res.y();
00121     hitStruct.resErrX = resXErr;
00122     hitStruct.resErrY = resYErr;
00123 
00124     // hitStruct.localX = lPhit.x();
00125     // hitStruct.localY = lPhit.y();
00126     // EM: use predictions for local coordinates
00127     hitStruct.localX = lPTrk.x();
00128     hitStruct.localY = lPTrk.y();
00129 
00130     // now calculate residuals taking global orientation of modules and radial topology in TID/TEC into account
00131     float resXprime(999.F), resYprime(999.F);
00132     float resXprimeErr(999.F), resYprimeErr(999.F);
00133     
00134     if(hit->detUnit()){ // is it a single physical module?
00135       const GeomDetUnit& detUnit = *(hit->detUnit());
00136       float uOrientation(-999.F), vOrientation(-999.F);
00137       float resXTopol(999.F), resYTopol(999.F);
00138       
00139       const Surface& surface = hit->detUnit()->surface();
00140       const BoundPlane& boundplane = hit->detUnit()->surface();
00141       const Bounds& bound = boundplane.bounds();
00142       
00143       float length = 0;
00144       float width = 0;
00145       float widthAtHalfLength = 0;
00146 
00147       LocalPoint lPModule(0.,0.,0.), lUDirection(1.,0.,0.), lVDirection(0.,1.,0.);
00148       GlobalPoint gPModule    = surface.toGlobal(lPModule),
00149                   gUDirection = surface.toGlobal(lUDirection),
00150                   gVDirection = surface.toGlobal(lVDirection);
00151       
00152       if (IntSubDetID == PixelSubdetector::PixelBarrel ||
00153           IntSubDetID == StripSubdetector::TIB || 
00154           IntSubDetID == StripSubdetector::TOB) {
00155         
00156         uOrientation = deltaPhi(gUDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
00157         vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;
00158         resXTopol = res.x();
00159         resYTopol = res.y();
00160         resXprimeErr = resXErr;
00161         resYprimeErr = resYErr;
00162 
00163         const RectangularPlaneBounds *rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
00164         hitStruct.inside = rectangularBound->inside(lPTrk);
00165         length = rectangularBound->length();
00166         width = rectangularBound->width();
00167         hitStruct.localXnorm = 2*hitStruct.localX/width;
00168         hitStruct.localYnorm = 2*hitStruct.localY/length;
00169         
00170       } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {
00171         
00172         uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
00173         vOrientation = deltaPhi(gVDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
00174         resXTopol = res.x();
00175         resYTopol = res.y();
00176         resXprimeErr = resXErr;
00177         resYprimeErr = resYErr;
00178         
00179         const RectangularPlaneBounds *rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
00180         hitStruct.inside = rectangularBound->inside(lPTrk);
00181         length = rectangularBound->length();
00182         width = rectangularBound->width();
00183         hitStruct.localXnorm = 2*hitStruct.localX/width;
00184         hitStruct.localYnorm = 2*hitStruct.localY/length;
00185         
00186       } else if (IntSubDetID == StripSubdetector::TID ||
00187                  IntSubDetID == StripSubdetector::TEC) {
00188         
00189         uOrientation = deltaPhi(gUDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
00190         vOrientation = gVDirection.perp() - gPModule.perp() >= 0. ? +1.F : -1.F;
00191         
00192         if (!dynamic_cast<const RadialStripTopology*>(&detUnit.type().topology()))continue;
00193         const RadialStripTopology& topol = dynamic_cast<const RadialStripTopology&>(detUnit.type().topology());
00194         
00195         MeasurementPoint measHitPos = topol.measurementPosition(lPHit);
00196         MeasurementPoint measTrkPos = topol.measurementPosition(lPTrk);
00197         
00198         MeasurementError measHitErr = topol.measurementError(lPHit,errHit);
00199         MeasurementError measTrkErr = topol.measurementError(lPTrk,errTrk);
00200         
00201         if (measHitErr.uu()<0. ||
00202             measHitErr.vv()<0. ||
00203             measTrkErr.uu()<0. ||
00204             measTrkErr.vv()<0.){
00205           edm::LogError("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
00206                                                       << "One of the squared error methods gives negative result"
00207                                                       << "\n\tmeasHitErr.uu()\tmeasHitErr.vv()\tmeasTrkErr.uu()\tmeasTrkErr.vv()"
00208                                                       << "\n\t" << measHitErr.uu()
00209                                                       << "\t" << measHitErr.vv()
00210                                                       << "\t" << measTrkErr.uu()
00211                                                       << "\t" << measTrkErr.vv();
00212           continue;
00213         }
00214           
00215         float localStripLengthHit = topol.localStripLength(lPHit);
00216         float localStripLengthTrk = topol.localStripLength(lPTrk);
00217         float phiHit = topol.stripAngle(measHitPos.x());
00218         float phiTrk = topol.stripAngle(measTrkPos.x());
00219         float r_0 = topol.originToIntersection();
00220         
00221         resXTopol = (phiTrk-phiHit)*r_0;
00222         //resYTopol = measTrkPos.y()*localStripLengthTrk - measHitPos.y()*localStripLengthHit;
00223         float cosPhiHit(cos(phiHit)), cosPhiTrk(cos(phiTrk)),
00224               sinPhiHit(sin(phiHit)), sinPhiTrk(sin(phiTrk));
00225         float l_0 = r_0 - topol.detHeight()/2;
00226         resYTopol = measTrkPos.y()*localStripLengthTrk - measHitPos.y()*localStripLengthHit + l_0*(1/cosPhiTrk - 1/cosPhiHit);
00227         
00228         resXprimeErr = std::sqrt(measHitErr.uu()+measTrkErr.uu())*topol.angularWidth()*r_0;
00229         //resYprimeErr = std::sqrt(measHitErr.vv()*localStripLengthHit*localStripLengthHit + measTrkErr.vv()*localStripLengthTrk*localStripLengthTrk);
00230         float helpSummand = l_0*l_0*topol.angularWidth()*topol.angularWidth()*(sinPhiHit*sinPhiHit/pow(cosPhiHit,4)*measHitErr.uu()
00231                                                                                + sinPhiTrk*sinPhiTrk/pow(cosPhiTrk,4)*measTrkErr.uu() );
00232         resYprimeErr = std::sqrt(measHitErr.vv()*localStripLengthHit*localStripLengthHit
00233                                  + measTrkErr.vv()*localStripLengthTrk*localStripLengthTrk + helpSummand );  
00234 
00235 
00236         const TrapezoidalPlaneBounds *trapezoidalBound = dynamic_cast < const TrapezoidalPlaneBounds * >(& bound);
00237         hitStruct.inside = trapezoidalBound->inside(lPTrk);
00238         length = trapezoidalBound->length();
00239         width  = trapezoidalBound->width();
00240         widthAtHalfLength = trapezoidalBound->widthAtHalfLength();
00241 
00242         int yAxisOrientation=trapezoidalBound->yAxisOrientation(); 
00243 // for trapezoidal shape modules, scale with as function of local y coordinate 
00244         float widthAtlocalY=width-(1-yAxisOrientation*2*lPTrk.y()/length)*(width-widthAtHalfLength); 
00245         hitStruct.localXnorm = 2*hitStruct.localX/widthAtlocalY;  
00246         hitStruct.localYnorm = 2*hitStruct.localY/length;
00247 
00248       } else {
00249         edm::LogWarning("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities" 
00250                                                       << "No valid tracker subdetector " << IntSubDetID;
00251         continue;
00252       }  
00253       
00254       resXprime = resXTopol*uOrientation;
00255       resYprime = resYTopol*vOrientation;
00256       
00257     } else { // not a detUnit, so must be a virtual 2D-Module
00258       // FIXME: at present only for det units residuals are calculated and filled in the hitStruct
00259       // But in principle this method should also be useable for the gluedDets (2D modules in TIB, TID, TOB, TEC)
00260       // In this case, only orientation should be taken into account for primeResiduals, but not the radial topology
00261       // At present, default values (999.F) are given out
00262     }
00263     
00264     hitStruct.resXprime = resXprime;
00265     hitStruct.resYprime = resYprime;
00266     hitStruct.resXprimeErr = resXprimeErr;
00267     hitStruct.resYprimeErr = resYprimeErr;
00268     
00269     hitStruct.rawDetId = IntRawDetID;
00270     hitStruct.phi = tsos.globalDirection().phi();
00271     hitStruct.eta = tsos.globalDirection().eta();
00272     
00273     // first try for overlap residuals
00274     // based on Code from Keith and Wolfgang
00275     if (itTraj+1 != tmColl.end()) {
00276       TransientTrackingRecHit::ConstRecHitPointer hit2 = (itTraj+1)->recHit();
00277       TrackerAlignableId ali1, ali2;
00278       if (hit2->isValid() && 
00279           ali1.typeAndLayerFromDetId(hit->geographicalId()) == ali2.typeAndLayerFromDetId(hit2->geographicalId())  &&
00280           hit2->geographicalId().rawId() !=  SiStripDetId(IntRawDetID).partnerDetId()) {            
00281           
00282         float overlapPath_;
00283         TrajectoryStateCombiner combiner_;
00284         AnalyticalPropagator propagator(&(*magneticField_));
00285         // forward and backward predicted states at module 1
00286         TrajectoryStateOnSurface fwdPred1 = (itTraj)->forwardPredictedState();
00287         TrajectoryStateOnSurface bwdPred1 = (itTraj)->backwardPredictedState();
00288         if ( !fwdPred1.isValid() || !bwdPred1.isValid() ) continue;
00289         // backward predicted state at module 2
00290         TrajectoryStateOnSurface bwdPred2 = (itTraj+1)->backwardPredictedState();
00291         TrajectoryStateOnSurface fwdPred2 = (itTraj+1)->forwardPredictedState();
00292         if ( !bwdPred2.isValid() ) continue;
00293         // extrapolation bwdPred2 to module 1
00294         TrajectoryStateOnSurface bwdPred2At1 = propagator.propagate(bwdPred2,fwdPred1.surface());
00295         if ( !bwdPred2At1.isValid() ) continue;
00296         // combination with fwdPred1 (ref. state, best estimate without hits 1 and 2)
00297         TrajectoryStateOnSurface comb1 = combiner_.combine(fwdPred1,bwdPred2At1);
00298         if ( !comb1.isValid() ) continue;
00299           
00300         //
00301         // propagation of reference parameters to module 2
00302         //
00303         std::pair<TrajectoryStateOnSurface,double> tsosWithS =
00304           propagator.propagateWithPath(comb1,bwdPred2.surface());
00305         TrajectoryStateOnSurface comb1At2 = tsosWithS.first;
00306         
00307         // Alternative possibility, not used at present
00308         // TrajectoryStateOnSurface comb1At2 = propagator.propagate(comb1,bwdPred2.surface());
00309           
00310         if ( !comb1At2.isValid() ) continue;
00311         overlapPath_ = tsosWithS.second;
00312           
00313         std::vector<GlobalPoint> predictedPositions;
00314         predictedPositions.push_back(comb1.globalPosition());
00315         predictedPositions.push_back(comb1At2.globalPosition());
00316         
00317         GlobalVector diff_pred = predictedPositions[0] - predictedPositions[1];
00318         
00319         TrajectoryStateOnSurface tsos2 = tsoscomb( (itTraj+1)->forwardPredictedState(), (itTraj+1)->backwardPredictedState() );
00320         align::LocalVector res2 = tsos2.localPosition() - hit2->localPosition();
00321         //float overlapresidual = res2.x() - res.x();
00322         float overlapresidual = diff_pred.x();
00323         
00324         hitStruct.overlapres = std::make_pair(hit2->geographicalId().rawId(),overlapresidual);
00325       }
00326     }
00327    
00328     v_avhitout.push_back(hitStruct);
00329   } 
00330 }
00331 
00332 void
00333 TrackerValidationVariables::fillTrackQuantities(const edm::Event& event, std::vector<AVTrackStruct> & v_avtrackout)
00334 {
00335   edm::InputTag TrjTrackTag = conf_.getParameter<edm::InputTag>("Tracks");
00336   edm::Handle<TrajTrackAssociationCollection> TrajTracksMap;
00337   event.getByLabel(TrjTrackTag, TrajTracksMap);
00338   LogDebug("TrackerValidationVariables") << "TrajTrack collection size " << TrajTracksMap->size();
00339   
00340   const Trajectory* trajectory;
00341   const reco::Track* track;
00342   
00343   for ( TrajTrackAssociationCollection::const_iterator iPair = TrajTracksMap->begin();
00344         iPair != TrajTracksMap->end();
00345         ++iPair) {
00346     
00347     trajectory = &(*(*iPair).key);
00348     track = &(*(*iPair).val);
00349     
00350     AVTrackStruct trackStruct;
00351     
00352     trackStruct.p = track->p();
00353     trackStruct.pt = track->pt();
00354     trackStruct.ptError = track->ptError();
00355     trackStruct.px = track->px();
00356     trackStruct.py = track->py();
00357     trackStruct.pz = track->pz();
00358     trackStruct.eta = track->eta();
00359     trackStruct.phi = track->phi();
00360     trackStruct.chi2 = track->chi2();
00361     trackStruct.chi2Prob= TMath::Prob(track->chi2(),track->ndof());
00362     trackStruct.normchi2 = track->normalizedChi2();
00363     GlobalPoint gPoint(track->vx(), track->vy(), track->vz());
00364     double theLocalMagFieldInInverseGeV = magneticField_->inInverseGeV(gPoint).z();
00365     trackStruct.kappa = -track->charge()*theLocalMagFieldInInverseGeV/track->pt();
00366     trackStruct.charge = track->charge();
00367     trackStruct.d0 = track->d0();
00368     trackStruct.dz = track->dz();
00369     trackStruct.numberOfValidHits = track->numberOfValidHits();
00370     trackStruct.numberOfLostHits = track->numberOfLostHits();
00371     
00372     fillHitQuantities(trajectory, trackStruct.hits);
00373     
00374     v_avtrackout.push_back(trackStruct);
00375   }
00376 }
00377 
00378 void
00379 TrackerValidationVariables::fillHitQuantities(const edm::Event& event, std::vector<AVHitStruct> & v_avhitout)
00380 {
00381   edm::Handle<std::vector<Trajectory> > trajCollectionHandle;
00382   event.getByLabel(conf_.getParameter<std::string>("trajectoryInput"), trajCollectionHandle);
00383   
00384   LogDebug("TrackerValidationVariables") << "trajColl->size(): " << trajCollectionHandle->size() ;
00385 
00386   for (std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(), itEnd = trajCollectionHandle->end(); 
00387        it!=itEnd;
00388        ++it) {
00389     
00390     fillHitQuantities(&(*it), v_avhitout);
00391   }
00392 }