diff --git a/[refs] b/[refs] index 18202300c3107..25c354fdc70e2 100644 --- a/[refs] +++ b/[refs] @@ -1,2 +1,2 @@ --- -refs/heads/CMSSW_7_2_X: 764abbbd874100a0f1af780bc986038e3c5f7af9 +refs/heads/CMSSW_7_2_X: c8c4ce1babb141cad7c04e96911c5f56b805dd6d diff --git a/trunk/RecoMuon/MuonSeedGenerator/BuildFile b/trunk/RecoMuon/MuonSeedGenerator/BuildFile index e21fd9fc498f4..3bda615d177f8 100644 --- a/trunk/RecoMuon/MuonSeedGenerator/BuildFile +++ b/trunk/RecoMuon/MuonSeedGenerator/BuildFile @@ -8,6 +8,7 @@ + @@ -26,6 +27,7 @@ + diff --git a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.cc b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.cc index 89f2bfc3de632..91fdd0440461f 100644 --- a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.cc +++ b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.cc @@ -1,8 +1,8 @@ /** * See header file for a description of this class. * - * $Date: 2006/06/21 17:54:37 $ - * $Revision: 1.8 $ + * $Date: 2006/06/27 13:42:41 $ + * $Revision: 1.9 $ * \author A. Vitelli - INFN Torino, V.Palichik * */ @@ -148,7 +148,85 @@ vector MuonSeedFinder::seeds(const edm::EventSetup& eSetup) cons } bool -MuonSeedFinder::createEndcapSeed(MuonTransientTrackingRecHit *me, +MuonSeedFinder::createEndcapSeed(MuonTransientTrackingRecHit *last, + vector& theSeeds, + const edm::EventSetup& eSetup) const { + + std::string metname = "Muon|RecoMuon|MuonSeedFinder"; + + edm::ESHandle field; + eSetup.get().get(field); + + AlgebraicSymMatrix mat(5,0) ; + + // this perform H.T() * parErr * H, which is the projection of the + // the measurement error (rechit rf) to the state error (TSOS rf) + // Legenda: + // H => is the 4x5 projection matrix + // parError the 4x4 parameter error matrix of the RecHit + + mat = last->parametersError().similarityT( last->projectionMatrix() ); + + // We want pT but it's not in RecHit interface, so we've put it within this class + float momentum = computePt(last,&*field); + // FIXME + float smomentum = 0.25; // FIXME!!!! + + MuonSeedFromRecHits seedCreator; + TrajectorySeed cscSeed = seedCreator.createSeed(momentum,smomentum,last,eSetup); + + theSeeds.push_back(cscSeed); + + // FIXME + return true; +} + + +float MuonSeedFinder::computePt(const MuonTransientTrackingRecHit *muon, const MagneticField *field) const { +// assume dZ = dPhi*R*C, here C = pZ/pT +// ======================================================================= +// ptc: I suspect the following comment should really be +// dZ/dPsi = 0.5*dz/dPhi +// which I can derive if I assume the particle has travelled in a circle +// projected onto the global xy plane, starting at the origin on the z-axis. +// Here Psi is the angle traced out in the xy plane by the projection of the +// helical path of the charged particle. The axis of the helix is assumed +// parallel to the main B field of the solenoid. +// ======================================================================= +// dZ/dPhi = 0.5*dZ/dPsi, here phi = atan2(y,x), psi = rho*s + +// ptc: If the local direction is effectively (0,0,1) or (0,0,-1) +// then it's ridiculous to follow this algorithm... just set some +// arbitrary 'high' value and note the sign is undetermined + +//@@ DO SOMETHING SANE WITH THESE TRAP VALUES + static float small = 1.e-06; + static float big = 1.e+10; + + LocalVector lod = muon->localDirection(); + if ( fabs(lod.x())globalPosition(); + GlobalVector gv = muon->globalDirection(); + float getx0 = gp.x(); + float getay = gv.y()/gv.z(); + float gety0 = gp.y(); + float getax = gv.x()/gv.z(); + float getz0 = gp.z(); + + float dZdPhi = 0.5f*gp.perp2()/(getx0*getay - gety0*getax); + float dZdT = getz0/gp.perp(); + float rho = dZdT/dZdPhi; + + // convert to pT (watch the sign !) + GlobalVector fld = field->inInverseGeV( gp ); + return -fld.z()/rho; +} + +bool +MuonSeedFinder::createEndcapSeed_OLD(MuonTransientTrackingRecHit *me, vector& theSeeds, const edm::EventSetup& eSetup) const { @@ -276,47 +354,3 @@ MuonSeedFinder::createEndcapSeed(MuonTransientTrackingRecHit *me, delete propagator; return result; } - - -float MuonSeedFinder::computePt(const MuonTransientTrackingRecHit *muon, const MagneticField *field) const { -// assume dZ = dPhi*R*C, here C = pZ/pT -// ======================================================================= -// ptc: I suspect the following comment should really be -// dZ/dPsi = 0.5*dz/dPhi -// which I can derive if I assume the particle has travelled in a circle -// projected onto the global xy plane, starting at the origin on the z-axis. -// Here Psi is the angle traced out in the xy plane by the projection of the -// helical path of the charged particle. The axis of the helix is assumed -// parallel to the main B field of the solenoid. -// ======================================================================= -// dZ/dPhi = 0.5*dZ/dPsi, here phi = atan2(y,x), psi = rho*s - -// ptc: If the local direction is effectively (0,0,1) or (0,0,-1) -// then it's ridiculous to follow this algorithm... just set some -// arbitrary 'high' value and note the sign is undetermined - -//@@ DO SOMETHING SANE WITH THESE TRAP VALUES - static float small = 1.e-06; - static float big = 1.e+10; - - LocalVector lod = muon->localDirection(); - if ( fabs(lod.x())globalPosition(); - GlobalVector gv = muon->globalDirection(); - float getx0 = gp.x(); - float getay = gv.y()/gv.z(); - float gety0 = gp.y(); - float getax = gv.x()/gv.z(); - float getz0 = gp.z(); - - float dZdPhi = 0.5f*gp.perp2()/(getx0*getay - gety0*getax); - float dZdT = getz0/gp.perp(); - float rho = dZdT/dZdPhi; - - // convert to pT (watch the sign !) - GlobalVector fld = field->inInverseGeV( gp ); - return -fld.z()/rho; -} diff --git a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.h b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.h index b47d53570f8ed..40d68ec1fb185 100644 --- a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.h +++ b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFinder.h @@ -8,8 +8,8 @@ * \author A. Vitelli - INFN Torino * \author porting R. Bellan - INFN Torino * - * $Date: 2006/05/30 13:50:23 $ - * $Revision: 1.5 $ + * $Date: 2006/06/27 13:42:41 $ + * $Revision: 1.6 $ * */ @@ -53,9 +53,13 @@ class MuonSeedFinder { bool createEndcapSeed(MuonTransientTrackingRecHit *me, std::vector& theSeeds, const edm::EventSetup& eSetup) const; + + bool createEndcapSeed_OLD(MuonTransientTrackingRecHit *me, + std::vector& theSeeds, + const edm::EventSetup& eSetup) const; float computePt(const MuonTransientTrackingRecHit *muon, const MagneticField *field) const; - + RecHitContainer theRhits; // put a parameterSet instead of diff --git a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.cc b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.cc index 2354217275156..0f639a74f1757 100644 --- a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.cc +++ b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.cc @@ -2,35 +2,44 @@ * See header file for a description of this class. * * - * $Date: 2006/06/21 17:48:52 $ - * $Revision: 1.6 $ + * $Date: 2006/06/27 13:42:41 $ + * $Revision: 1.7 $ * \author A. Vitelli - INFN Torino, V.Palichik * */ #include "RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.h" #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h" +#include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h" +#include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h" +#include "RecoMuon/Navigation/interface/MuonNavigationSchool.h" +#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h" #include "Geometry/Surface/interface/BoundPlane.h" #include "Geometry/Surface/interface/SimpleCylinderBounds.h" #include "Geometry/Surface/interface/BoundCylinder.h" #include "Geometry/CommonDetUnit/interface/GeomDet.h" -// FIXME -#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" - -#include "MagneticField/Engine/interface/MagneticField.h" -#include "DataFormats/Common/interface/OwnVector.h" #include "MagneticField/Engine/interface/MagneticField.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/DetLayers/interface/NavigationSetter.h" + +#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h" +#include "DataFormats/Common/interface/OwnVector.h" +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" + #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h" -#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" - #include "gsl/gsl_statistics.h" using namespace std; @@ -50,14 +59,14 @@ TrajectorySeed MuonSeedFromRecHits::seed(const edm::EventSetup& eSetup) const { computePtWithoutVtx(pt, spt); // some dump... - LogDebug(metname) << " Pt MB1 @vtx: " << pt[0] << " w: " << spt[0] << endl + LogDebug(metname) << " Pt MB1 @vtx: " << pt[0] << " w: " << spt[0] << "\n" << " Pt MB2 @vtx: " << pt[1] << " w: " << spt[1]<< endl ; - LogDebug(metname) << " Pt MB2-MB1 " << pt[2] << " w: " << spt[2]<< endl - << " Pt MB3-MB1 " << pt[3] << " w: " << spt[3]<< endl - << " Pt MB3-MB2 " << pt[4] << " w: " << spt[4]<< endl - << " Pt MB4-MB1 " << pt[5] << " w: " << spt[5]<< endl - << " Pt MB4-MB2 " << pt[6] << " w: " << spt[6]<< endl + LogDebug(metname) << " Pt MB2-MB1 " << pt[2] << " w: " << spt[2]<< "\n" + << " Pt MB3-MB1 " << pt[3] << " w: " << spt[3]<< "\n" + << " Pt MB3-MB2 " << pt[4] << " w: " << spt[4]<< "\n" + << " Pt MB4-MB1 " << pt[5] << " w: " << spt[5]<< "\n" + << " Pt MB4-MB2 " << pt[6] << " w: " << spt[6]<< "\n" << " Pt MB4-MB3 " << pt[7] << " w: " << spt[7]<< endl ; /// now combine all pt estimate @@ -67,7 +76,9 @@ TrajectorySeed MuonSeedFromRecHits::seed(const edm::EventSetup& eSetup) const { LogDebug(metname) << " Seed Pt :" << ptmean << "+/-" << sptmean << endl; - return createSeed(ptmean, sptmean,eSetup); + // take the best candidate + MuonTransientTrackingRecHit* last = best_cand(); + return createSeed(ptmean, sptmean,last,eSetup); } MuonTransientTrackingRecHit* MuonSeedFromRecHits::best_cand() const { @@ -461,32 +472,47 @@ void MuonSeedFromRecHits::computeBestPt(double* pt, TrajectorySeed MuonSeedFromRecHits::createSeed(float ptmean, float sptmean, + MuonTransientTrackingRecHit* last, const edm::EventSetup& eSetup) const{ - + std::string metname = "Muon|RecoMuon|MuonSeedFromRecHits"; + + MuonPatternRecoDumper debug; edm::ESHandle field; eSetup.get().get(field); - TrajectorySeed result; - if ( fabs(ptmean) < 3. ) ptmean = 100. * ptmean/fabs(ptmean) ; + // FIXME: put it into a parameter set + edm::ESHandle estimator; + eSetup.get().get("Chi2MeasurementEstimator",estimator); + + // FIXME: put it into a parameter set + edm::ESHandle propagator; + eSetup.get().get("SteppingHelixPropagatorOpposite",propagator); - AlgebraicVector t(4); - AlgebraicSymMatrix mat(5,0) ; + // FIXME: put it into a parameter set! + double theMinMomentum = 3.0; - MuonTransientTrackingRecHit* last = best_cand(); + // Take the whole tracking geometry + edm::ESHandle muonLayers; + eSetup.get().get(muonLayers); + + // set the proper navigation school + MuonNavigationSchool school(&*muonLayers); + NavigationSetter setter(school); + + TrajectorySeed result; - // float p_x = last.localDirection().x(); - // float p_y = last.localDirection().y(); - // float p_z = last.localDirection().z(); + // Minimal pt + if ( fabs(ptmean) < theMinMomentum ) ptmean = theMinMomentum * ptmean/fabs(ptmean) ; - // float p_norm = ptmean/sqrt(p_x*p_x + p_z*p_z); - // t[0] = p_x/p_z; - // t[1] = p_y/p_z; - // t[2] = last.localPosition().x(); - // t[3] = last.localPosition().y(); - // LocalTrajectoryParameters param1(1./p_norm,t[0],t[1],t[2],t[3],-1); + // FIXME!! was: + // if ( fabs(ptmean) < 3. ) ptmean = 100. * ptmean/fabs(ptmean) ; + AlgebraicVector t(4); + AlgebraicSymMatrix mat(5,0) ; + + // Fill the LocalTrajectoryParameters LocalPoint segPos=last->localPosition(); GlobalVector mom=last->globalPosition()-GlobalPoint(); GlobalVector polar(GlobalVector::Spherical(mom.theta(), @@ -495,81 +521,112 @@ TrajectorySeed MuonSeedFromRecHits::createSeed(float ptmean, polar *=fabs(ptmean)/polar.perp(); LocalVector segDirFromPos=last->det()->toLocal(polar); int charge=(int)(ptmean/fabs(ptmean)); + LocalTrajectoryParameters param(segPos,segDirFromPos, charge); - // mat[1][1] = last->parametersError()[1][1]; - // mat[3][3] = last->parametersError()[0][0]; - // mat[1][3] = last->parametersError()[1][0]; - // mat[2][2] = last->parametersError()[1][1]; - // mat[4][4] = last->parametersError()[0][0]; - // mat[2][4] = last->parametersError()[1][0]; - // this perform H.T() * parErr * H, which is the projection of the // the measurement error (rechit rf) to the state error (TSOS rf) // Legenda: // H => is the 4x5 projection matrix // parError the 4x4 parameter error matrix of the RecHit - - // FIXME Use this!!!!!!!! mat = last->parametersError().similarityT( last->projectionMatrix() ); - float p_err = sqr(sptmean/(ptmean*ptmean)); mat[0][0]= p_err; LocalTrajectoryError error(mat); + // Create the TrajectoryStateOnSurface TrajectoryStateOnSurface tsos(param, error, last->det()->surface(),&*field); const FreeTrajectoryState state = *(tsos.freeState()); - - LogDebug(metname) << " Before extr.: pos. :" << state.position() - <<" eta "<geographicalId(); + const DetLayer *initialLayer = muonLayers->idToLayer(id); + + if(last->isDT()){ + DTChamberId chamberId(id.rawId()); + LogDebug(metname)<<"Starting from: "<isCSC()){ + CSCDetId chamberId(id.rawId()); + LogDebug(metname)<<"Starting from: "<isRPC()){ + RPCDetId chamberId(id.rawId()); + LogDebug(metname)<<"Starting from: "< cyl= - new BoundCylinder( pos, rot, SimpleCylinderBounds( 399., 401., -1200., 1200.)); - - // FIXME - SteppingHelixPropagator prop(&*field,oppositeToMomentum); + // ask for compatible layers + vector detLayers = initialLayer->compatibleLayers(state,oppositeToMomentum); - const TrajectoryStateOnSurface trj = prop.propagate( state, *cyl ); - if ( trj.isValid() ) { - const FreeTrajectoryState e_state = *trj.freeTrajectoryState(); + std::vector detsWithStates; + + if(detLayers.size()){ + // ask for compatible dets + debug.dumpLayer(detLayers.back(),metname); + detsWithStates = detLayers.back()->compatibleDets(tsos, *propagator, *estimator); + LogDebug(metname)<<"Number of compatible dets: "<isDT()){ +// DTChamberId internalChamberId( newTSOSDet->geographicalId().rawId() ); +// LogDebug(metname)<<"Most compatible det: "<isCSC()){ +// CSCDetId internalChamberId( newTSOSDet->geographicalId().rawId() ); +// LogDebug(metname)<<"Most compatible det: "<isRPC()){ +// RPCDetId internalChamberId( newTSOSDet->geographicalId().rawId() ); +// LogDebug(metname)<<"Most compatible det: "<geographicalId().rawId()); + + edm::OwnVector container; + TrajectorySeed theSeed(*seedTSOS,container,oppositeToMomentum); + + result = theSeed; + } + else LogDebug(metname) << " Extrapolation failed" << endl; + } + else{ // Transform it in a TrajectoryStateOnSurface TrajectoryStateTransform tsTransform; - - // FIXME FIXME TEST -// PTrajectoryStateOnDet *seedTSOS = -// tsTransform.persistentState( trj ,last->geographicalId().rawId()); - // FIXME the tsos is defined on the "me" surface, this must be changed!!! PTrajectoryStateOnDet *seedTSOS = tsTransform.persistentState( tsos ,last->geographicalId().rawId()); - //<< FIXME would be: - - // TrajectorySeed theSeed(e_state, rechitcontainer,oppositeToMomentum); - // But is: edm::OwnVector container; - // container.push_back(last->hit()->clone()); - TrajectorySeed theSeed(*seedTSOS,container,oppositeToMomentum); - //>> is it right?? - + result = theSeed; - } else { - LogDebug(metname) << " Extr. failed" << endl; } - + return result; } diff --git a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.h b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.h index 81da556f3f162..0a5ba7b0a0e3d 100644 --- a/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.h +++ b/trunk/RecoMuon/MuonSeedGenerator/src/MuonSeedFromRecHits.h @@ -22,12 +22,14 @@ class MuonTransientTrackingRecHit; class RecHit; class BoundPlane; +class GeomDet; namespace edm {class EventSetup;} class MuonSeedFromRecHits { typedef std::vector RecHitContainer; typedef RecHitContainer::const_iterator RecHitIterator; + typedef std::pair DetWithState; public: MuonSeedFromRecHits(){ @@ -39,6 +41,8 @@ class MuonSeedFromRecHits { unsigned int nrhit() const { return theRhits.size(); } private: + friend class MuonSeedFinder; + MuonTransientTrackingRecHit *best_cand() const; // was // TrackingRecHit best_cand() const; @@ -47,7 +51,9 @@ class MuonSeedFromRecHits { void computePtWithoutVtx(double* pt, double* spt) const; void computeBestPt(double* pt, double* spt, float& ptmean, float& sptmean) const; - TrajectorySeed createSeed(float ptmean, float sptmean,const edm::EventSetup& eSetup) const; + TrajectorySeed createSeed(float ptmean, float sptmean, + MuonTransientTrackingRecHit *last, + const edm::EventSetup& eSetup) const; private: RecHitContainer theRhits;