CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    MuonHIPAlignmentRefitter
00004 // Class:      MuonHIPAlignmentRefitter
00005 // 
00013 //
00014 // Original Author:  Jim Pivarski
00015 //         Created:  Wed Dec 12 13:31:55 CST 2007
00016 // $Id: MuonHIPAlignmentRefitter.cc,v 1.1 2008/02/14 15:39:58 pivarski Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 
00024 // user include files
00025 #include "FWCore/Framework/interface/Frameworkfwd.h"
00026 #include "FWCore/Framework/interface/EDProducer.h"
00027 
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/EventSetup.h"
00030 #include "FWCore/Framework/interface/MakerMacros.h"
00031 
00032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00033 
00034 #include "DataFormats/MuonReco/interface/Muon.h"
00035 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00036 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
00037 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00038 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00039 #include "DataFormats/TrackReco/interface/Track.h"
00040 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00041 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00042 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00043 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
00044 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00045 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00046 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00047 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00048 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00049 #include "MagneticField/Engine/interface/MagneticField.h"
00050 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00051 #include "CondFormats/Alignment/interface/Definitions.h"
00052 
00053 //
00054 // class decleration
00055 //
00056 
00057 class MuonHIPAlignmentRefitter : public edm::EDProducer {
00058    public:
00059       explicit MuonHIPAlignmentRefitter(const edm::ParameterSet&);
00060       ~MuonHIPAlignmentRefitter();
00061 
00062    private:
00063       virtual void beginJob(const edm::EventSetup&) ;
00064       virtual void produce(edm::Event&, const edm::EventSetup&);
00065       virtual void endJob() ;
00066       
00067       // ----------member data ---------------------------
00068 
00069       edm::InputTag m_muonSource;
00070       std::string m_propagatorSource;
00071 
00072       TrackTransformer *m_trackTransformer;
00073 };
00074 
00075 //
00076 // constants, enums and typedefs
00077 //
00078 
00079 //
00080 // static data member definitions
00081 //
00082 
00083 //
00084 // constructors and destructor
00085 //
00086 MuonHIPAlignmentRefitter::MuonHIPAlignmentRefitter(const edm::ParameterSet& iConfig)
00087 {
00088    m_muonSource = iConfig.getParameter<edm::InputTag>("MuonSource");
00089    m_propagatorSource = iConfig.getParameter<std::string>("MuonPropagator");
00090 
00091   m_trackTransformer = new TrackTransformer(iConfig.getParameter<edm::ParameterSet>("TrackerTrackTransformer"));
00092   
00093   produces<std::vector<Trajectory> >();
00094   produces<TrajTrackAssociationCollection>();
00095 }
00096 
00097 
00098 MuonHIPAlignmentRefitter::~MuonHIPAlignmentRefitter()
00099 {
00100  
00101    // do anything here that needs to be done at desctruction time
00102    // (e.g. close files, deallocate resources etc.)
00103 
00104 }
00105 
00106 
00107 //
00108 // member functions
00109 //
00110 
00111 // ------------ method called to produce the data  ------------
00112 void
00113 MuonHIPAlignmentRefitter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00114 {
00115    m_trackTransformer->setServices(iSetup);
00116 
00117    edm::Handle<reco::MuonCollection> muons;
00118    iEvent.getByLabel(m_muonSource, muons);
00119 
00120    edm::ESHandle<Propagator> propagator;
00121    iSetup.get<TrackingComponentsRecord>().get(m_propagatorSource, propagator);
00122 
00123    edm::ESHandle<DTGeometry> dtGeometry;
00124    iSetup.get<MuonGeometryRecord>().get(dtGeometry);
00125 
00126    edm::ESHandle<CSCGeometry> cscGeometry;
00127    iSetup.get<MuonGeometryRecord>().get(cscGeometry);
00128 
00129    edm::ESHandle<MagneticField> magneticField;
00130    iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
00131 
00132    for (reco::MuonCollection::const_iterator muon = muons->begin();  muon != muons->end();  ++muon) {
00133       std::vector<Trajectory> trackerTrajectories = m_trackTransformer->transform(*muon->track());
00134 
00135       if (trackerTrajectories.size() == 1) {
00136          const Trajectory trackerTrajectory = *(trackerTrajectories.begin());
00137          TrajectoryStateOnSurface tracker_tsos = trackerTrajectory.lastMeasurement().forwardPredictedState();
00138          TrajectoryStateOnSurface last_tsos = trackerTrajectory.lastMeasurement().forwardPredictedState();
00139    
00140          int last_chamber = 0;
00141          std::vector<double> last_x, last_xerr2, last_zforx;
00142          std::vector<double> last_y, last_yerr2, last_zfory;
00143          std::vector<const TrackingRecHit*> last_hits;
00144 
00145          for (trackingRecHit_iterator hit = muon->standAloneMuon()->recHitsBegin();  hit != muon->standAloneMuon()->recHitsEnd();  ++hit) {
00146             DetId id = (*hit)->geographicalId();
00147 
00148             if (id.det() == DetId::Muon  &&  id.subdetId() != MuonSubdetId::RPC) {
00149                int chamberId = -1;
00150                if (id.subdetId() == MuonSubdetId::DT) {
00151                   DTChamberId dtChamberId(id.rawId());
00152                   chamberId = dtChamberId.rawId();
00153                }
00154                else if (id.subdetId() == MuonSubdetId::CSC) {
00155                   CSCDetId cscChamberId(id.rawId());
00156                   cscChamberId = CSCDetId(cscChamberId.endcap(), cscChamberId.station(), cscChamberId.ring(), cscChamberId.chamber(), 0);
00157                   chamberId = cscChamberId.rawId();
00158                }
00159                else assert(false);
00160 
00161                // First chamber: set "last_chamber" variable
00162                if (last_chamber == 0) {
00163                   last_chamber = chamberId;
00164                }
00165                // New chamber: re-orient propagator based on linear fit to chamber's 1-12 dof
00166                else if (last_chamber != chamberId) {
00167                   GlobalPoint globalPosition = last_tsos.globalPosition();
00168                   GlobalVector globalDirection = last_tsos.globalDirection();
00169 
00170                   const Surface* surface;
00171                   if (DetId(last_chamber).subdetId() == MuonSubdetId::DT) {
00172                      surface = &(dtGeometry->idToDet(DetId(last_chamber))->surface());
00173                   }
00174                   else {
00175                      surface = &(cscGeometry->idToDet(DetId(last_chamber))->surface());
00176                   }
00177 
00178                   LocalPoint localPosition = surface->toLocal(globalPosition);
00179                   LocalVector localDirection = surface->toLocal(globalDirection);
00180                   // we want a normalization in which dz/dz == 1
00181                   localDirection /= localDirection.z();
00182 
00183                   if (last_zforx.size() > 3) {
00184                      // weighted linear fit to x(z)
00185                      double x_S = 0;
00186                      double x_Sz = 0;
00187                      double x_Sx = 0;
00188                      double x_Szz = 0;
00189                      double x_Szx = 0;
00190                      
00191                      for (unsigned int i = 0;  i < last_zforx.size();  i++) {
00192                         double x = last_x[i];
00193                         double z = last_zforx[i];
00194                         double x_weight = 1./last_xerr2[i];
00195 
00196                         x_S += x_weight;
00197                         x_Sz += z * x_weight;
00198                         x_Sx += x * x_weight;
00199                         x_Szz += z*z * x_weight;
00200                         x_Szx += z*x * x_weight;
00201                      }
00202 
00203                      double x_delta = x_S * x_Szz - x_Sz*x_Sz;
00204                      double x_intercept = (x_Szz * x_Sx - x_Sz * x_Szx) / x_delta;
00205 //                   double x_intercept_err2 = x_Szz / x_delta;
00206                      double x_slope = (x_S * x_Szx - x_Sz * x_Sx) / x_delta;
00207 //                   double x_slope_err2 = x_S / x_delta;
00208 //                   double x_covariance = -x_Sz / x_delta;
00209 
00210                      double x_on_surface = x_intercept + x_slope * localPosition.z();
00211 //                   double x_on_surface_err2 = x_intercept_err2 + localPosition.z()*localPosition.z() * x_slope_err2 + localPosition.z() * x_covariance;
00212 
00213                      localPosition = LocalPoint(x_on_surface, localPosition.y(), localPosition.z());
00214                      localDirection = LocalVector(x_slope, localDirection.y(), 1);
00215                   }
00216 
00217                   if (last_zfory.size() > 3) {
00218                      // weighted linear fit to y(z)
00219                      double y_S = 0;
00220                      double y_Sz = 0;
00221                      double y_Sy = 0;
00222                      double y_Szz = 0;
00223                      double y_Szy = 0;
00224                      for (unsigned int i = 0;  i < last_zfory.size();  i++) {
00225                         double y = last_y[i];
00226                         double z = last_zfory[i];
00227                         double y_weight = 1./last_yerr2[i];
00228 
00229                         y_S += y_weight;
00230                         y_Sz += z * y_weight;
00231                         y_Sy += y * y_weight;
00232                         y_Szz += z*z * y_weight;
00233                         y_Szy += z*y * y_weight;
00234                      }
00235 
00236                      double y_delta = y_S * y_Szz - y_Sz*y_Sz;
00237                      double y_intercept = (y_Szz * y_Sy - y_Sz * y_Szy) / y_delta;
00238 //                   double y_intercept_err2 = y_Szz / y_delta;
00239                      double y_slope = (y_S * y_Szy - y_Sz * y_Sy) / y_delta;
00240 //                   double y_slope_err2 = y_S / y_delta;
00241 //                   double y_covariance = -y_Sz / y_delta;
00242                   
00243                      double y_on_surface = y_intercept + y_slope * localPosition.z();
00244 //                   double y_on_surface_err2 = y_intercept_err2 + localPosition.z()*localPosition.z() * y_slope_err2 + localPosition.z() * y_covariance;
00245 
00246                      localPosition = LocalPoint(localPosition.x(), y_on_surface, localPosition.z());
00247                      localDirection = LocalVector(localDirection.x(), y_slope, 1);
00248                   }
00249 
00250                   if (last_zforx.size() > 3  ||  last_zfory.size() > 3) {
00251                      // make the direction normalized again
00252                      localDirection /= sqrt(localDirection.mag());
00253 
00254                      globalPosition = surface->toGlobal(localPosition);
00255                      globalDirection = surface->toGlobal(localDirection);
00256                      GlobalVector globalMomentum = globalDirection * fabs(1./last_tsos.signedInverseMomentum());
00257 
00258                      AlgebraicSymMatrix55 error(last_tsos.curvilinearError().matrix());
00259                      GlobalTrajectoryParameters globalTrajectoryParameters(globalPosition, globalMomentum, last_tsos.charge(), &*magneticField);
00260                      FreeTrajectoryState freeTrajectoryState(globalTrajectoryParameters, CurvilinearTrajectoryError(error));
00261                      last_tsos = TrajectoryStateOnSurface(freeTrajectoryState, last_tsos.surface());
00262                   }
00263 
00264                   last_x.clear();
00265                   last_y.clear();
00266                   last_zforx.clear();
00267                   last_zfory.clear();
00268                   last_xerr2.clear();
00269                   last_yerr2.clear();
00270                   last_chamber = chamberId;
00271 
00272                   last_hits.clear();
00273                }
00274 
00275                // Collect hits and propagate to surface
00276                LocalPoint chamberPoint;
00277                align::RotationType layerRot, chamberRot;
00278 
00279                last_hits.push_back(&**hit);
00280                GlobalPoint global_hit_position;
00281                TrajectoryStateOnSurface extrapolation, fromtracker;
00282                if (id.subdetId() == MuonSubdetId::DT) {
00283                   extrapolation = propagator->propagate(last_tsos, dtGeometry->idToDet(id)->surface());
00284                   fromtracker = propagator->propagate(tracker_tsos, dtGeometry->idToDet(id)->surface());
00285                   LocalPoint local = (*hit)->localPosition();
00286                   if (extrapolation.isValid()) {
00287                      // The extrapolated track is used to set a y position for the hit;
00288                      // y positions are always small corrections to chamber coordinates,
00289                      // applicable only when layers are rotated inside the chamber
00290                      local = LocalPoint(local.x(), extrapolation.localPosition().y(), 0);
00291                   }
00292 
00293                   GlobalPoint globalPoint = dtGeometry->idToDet(id)->surface().toGlobal(local);
00294                   chamberPoint = dtGeometry->idToDet(chamberId)->surface().toLocal(globalPoint);
00295                   layerRot = dtGeometry->idToDet(id)->surface().rotation();
00296                   chamberRot = dtGeometry->idToDet(chamberId)->surface().rotation();
00297                   global_hit_position = globalPoint;
00298                }
00299                else if (id.subdetId() == MuonSubdetId::CSC) {
00300                   extrapolation = propagator->propagate(last_tsos, cscGeometry->idToDet(id)->surface());
00301                   fromtracker = propagator->propagate(tracker_tsos, cscGeometry->idToDet(id)->surface());
00302                   GlobalPoint globalPoint = cscGeometry->idToDet(id)->surface().toGlobal((*hit)->localPosition());
00303                   chamberPoint = cscGeometry->idToDet(chamberId)->surface().toLocal(globalPoint);
00304                   layerRot = cscGeometry->idToDet(id)->surface().rotation();
00305                   chamberRot = cscGeometry->idToDet(chamberId)->surface().rotation();
00306                   global_hit_position = globalPoint;
00307                }
00308                else assert(false);
00309 
00310                double lRxx = layerRot.xx();
00311                double lRxy = layerRot.xy();
00312                double lRxz = layerRot.xz();
00313                double lRyx = layerRot.yx();
00314                double lRyy = layerRot.yy();
00315                double lRyz = layerRot.yz();
00316 //             double lRzx = layerRot.zx();
00317 //             double lRzy = layerRot.zy();
00318 //             double lRzz = layerRot.zz();
00319                double cRxx = chamberRot.xx();
00320                double cRxy = chamberRot.xy();
00321                double cRxz = chamberRot.xz();
00322                double cRyx = chamberRot.yx();
00323                double cRyy = chamberRot.yy();
00324                double cRyz = chamberRot.yz();
00325 //             double cRzx = chamberRot.zx();
00326 //             double cRzy = chamberRot.zy();
00327 //             double cRzz = chamberRot.zz();
00328                
00329                if (id.subdetId() == MuonSubdetId::DT) {
00330                   if (fabs(cRxx*lRxx + cRxy*lRxy + cRxz*lRxz) > fabs(cRxx*lRyx + cRxy*lRyy + cRxz*lRyz)) {
00331                      last_x.push_back(chamberPoint.x());
00332                      last_xerr2.push_back(pow((cRxx*lRxx + cRxy*lRxy + cRxz*lRxz), 2) * (*hit)->localPositionError().xx() +
00333                                           pow((cRxx*lRyx + cRxy*lRyy + cRxz*lRyz), 2) * (*hit)->localPositionError().yy());
00334                      last_zforx.push_back(chamberPoint.z());
00335                   }
00336                   else {
00337                      last_y.push_back(chamberPoint.y());
00338                      last_yerr2.push_back(pow((cRyx*lRxx + cRyy*lRxy + cRyz*lRxz), 2) * (*hit)->localPositionError().xx() +
00339                                           pow((cRyx*lRyx + cRyy*lRyy + cRyz*lRyz), 2) * (*hit)->localPositionError().yy());
00340                      last_zfory.push_back(chamberPoint.z());
00341                   }
00342                }
00343                else if (id.subdetId() == MuonSubdetId::CSC) {
00344                   last_x.push_back(chamberPoint.x());
00345                   last_xerr2.push_back(pow((cRxx*lRxx + cRxy*lRxy + cRxz*lRxz), 2) * (*hit)->localPositionError().xx() +
00346                                        pow((cRxx*lRyx + cRxy*lRyy + cRxz*lRyz), 2) * (*hit)->localPositionError().yy());
00347                   last_zforx.push_back(chamberPoint.z());
00348 
00349                   last_y.push_back(chamberPoint.y());
00350                   last_yerr2.push_back(pow((cRyx*lRxx + cRyy*lRxy + cRyz*lRxz), 2) * (*hit)->localPositionError().xx() +
00351                                        pow((cRyx*lRyx + cRyy*lRyy + cRyz*lRyz), 2) * (*hit)->localPositionError().yy());
00352                   last_zfory.push_back(chamberPoint.z());
00353                }
00354                else assert(false);
00355 
00356                if (extrapolation.isValid()) {
00357                   last_tsos = extrapolation;
00358                }
00359             } // end if Muon and not RPC
00360          } // end loop over standAlone hits
00361       } // end if successfully refit tracker track
00362    } // end loop over muons
00363 
00364    std::auto_ptr<std::vector<Trajectory> > trajectoryCollection(new std::vector<Trajectory>);
00365    std::auto_ptr<TrajTrackAssociationCollection> trajTrackMap(new TrajTrackAssociationCollection);
00366    iEvent.put(trajectoryCollection);
00367    iEvent.put(trajTrackMap);
00368 }
00369 
00370 // ------------ method called once each job just before starting event loop  ------------
00371 void 
00372 MuonHIPAlignmentRefitter::beginJob(const edm::EventSetup&)
00373 {
00374 }
00375 
00376 // ------------ method called once each job just after ending the event loop  ------------
00377 void 
00378 MuonHIPAlignmentRefitter::endJob() {
00379 }
00380 
00381 //define this as a plug-in
00382 DEFINE_FWK_MODULE(MuonHIPAlignmentRefitter);