From 439680eb7a71c02b605b9606c4a861e582e57596 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Sat, 20 Feb 2016 15:27:38 +0100 Subject: [PATCH 01/18] optimize Gsf --- .../TrackProducer/src/TrackProducerAlgorithm.cc | 1 + .../GsfTracking/src/GsfBetheHeitlerUpdator.cc | 16 +++++++++------- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc index 85348ce2b0506..5c1b8cebfd178 100644 --- a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc +++ b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc @@ -342,6 +342,7 @@ TrackProducerAlgorithm::buildTrack (const TrajectoryFitter * the math::XYZVector mom( p.x(), p.y(), p.z() ); LogDebug("GsfTrackProducer") << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag(); + std::cout << "GsfTrackProducer " << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag() << std::endl; auto theTrack = new reco::GsfTrack(theTraj->chiSquared(), int(ndof),//FIXME fix weight() in TrackingRecHit diff --git a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc index 1758715a5fcd2..8aa044ca42ea2 100644 --- a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc @@ -20,7 +20,7 @@ namespace { /// Second moment of the Bethe-Heitler distribution (in z=E/E0) inline float BetheHeitlerVariance (const float rl) { - const float l3ol2 = std::log(3.)/std::log(2.); + constexpr float l3ol2 = std::log(3.)/std::log(2.); float mean = BetheHeitlerMean(rl); return unsafe_expf<4>(-rl*l3ol2) - mean*mean; } @@ -143,10 +143,10 @@ GsfBetheHeitlerUpdator::compute (const TrajectoryStateOnSurface& TSoS, // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside), // then convert sig(p) to sig(1/p). // - effects[i].deltaP += p*(mixture[i].second-1); + effects[i].deltaP += p*(mixture[i].second-1.f); // float f = 1./p/mixture[i].second/mixture[i].second; // patch to ensure consistency between for- and backward propagation - float f = 1./p/mixture[i].second; + float f = 1.f/(p*mixture[i].second); varPinv = f*f*mixture[i].third; } else { @@ -154,8 +154,8 @@ GsfBetheHeitlerUpdator::compute (const TrajectoryStateOnSurface& TSoS, // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside // convert to obtain equivalent delta(p) // - effects[i].deltaP += p*(1/mixture[i].second-1); - varPinv = mixture[i].third/p/p; + effects[i].deltaP += p*(1.f/mixture[i].second-1.f); + varPinv = mixture[i].third/(p*p); } using namespace materialEffect; effects[i].deltaCov[elos] += varPinv; @@ -181,7 +181,8 @@ GsfBetheHeitlerUpdator::getMixtureParameters (const float rl, float vz = thePolyVars[i](rl); if ( theTransformationCode ) vz = unsafe_expf<4>(vz); - else vz = vz*vz; + else + vz = vz*vz; mixture[i]=Triplet(weight,z,vz); } @@ -202,8 +203,9 @@ GsfBetheHeitlerUpdator::correctWeights (GSContainer mixture[]) const // // rescale to obtain 1 // + wsum = 1.f/wsum; for ( int i=0; i Date: Sat, 20 Feb 2016 16:07:07 +0100 Subject: [PATCH 02/18] remove more divisions --- .../src/MultiGaussianStateCombiner1D.cc | 7 +-- .../src/MultiTrajectoryStateAssembler.cc | 6 +-- .../GsfTools/src/MultiTrajectoryStateMode.cc | 44 +++++++++---------- 3 files changed, 29 insertions(+), 28 deletions(-) diff --git a/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc b/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc index 5f07b8f76c977..6deaa5c892c2e 100644 --- a/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc +++ b/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc @@ -69,9 +69,10 @@ MultiGaussianStateCombiner1D::combine(const VSC& theComponents) const measCovar = 0.; weightSum = 0.; } else { - meanMean /= weightSum; - measCovar1 *= (1./weightSum); - measCovar2 *= (1./weightSum/weightSum); + weightSum = 1./weightSum; + meanMean *= weightSum; + measCovar1 *= weightSum; + measCovar2 *= weightSum*weightSum; measCovar = measCovar1 + measCovar2; // measCovar = measCovar1/weightSum + measCovar2/weightSum/weightSum; } diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index 071d6fa049582..538dfa85b41c8 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -136,7 +136,7 @@ MultiTrajectoryStateAssembler::prepareCombinedState () { // Check for minimum fraction of valid states // double allWeights(theValidWeightSum+theInvalidWeightSum); - if ( theInvalidWeightSum>0. && (theValidWeightSum/allWeights)0. && theValidWeightSum globalMomentum()); - double px = mom.x(); - double py = mom.y(); - double pz = mom.z(); - double p = mom.mag(); - double pt2 = mom.perp2(); - double phi = mom.phi(); - double eta = mom.eta(); + auto px = mom.x(); + auto py = mom.y(); + auto pz = mom.z(); + auto op = 1./mom.mag(); + auto opt2 = 1./mom.perp2(); + auto phi = mom.phi(); + auto eta = mom.eta(); // jacobian - jacobian(0,0) = px/p; - jacobian(0,1) = py/p; - jacobian(0,2) = pz/p; - jacobian(1,0) = py/pt2; - jacobian(1,1) = -px/pt2; + jacobian(0,0) = px*op; + jacobian(0,1) = py*op; + jacobian(0,2) = pz*op; + jacobian(1,0) = py*opt2; + jacobian(1,1) = -px*opt2; jacobian(1,2) = 0; - jacobian(2,0) = px*pz/(pt2*p); - jacobian(2,1) = py*pz/(pt2*p); - jacobian(2,2) = -1./p; + jacobian(2,0) = px*pz*opt2*op; + jacobian(2,1) = py*pz*opt2*op; + jacobian(2,2) = -op; // extraction of the momentum part from the 6x6 cartesian error matrix // and conversion to p-phi-eta covCart = ic->cartesianError().matrix().Sub(3,3); covPPhiEta = ROOT::Math::Similarity(jacobian,covCart); - pStates.push_back(SingleGaussianState1D(p,covPPhiEta(0,0),ic->weight())); + pStates.push_back(SingleGaussianState1D(1/op,covPPhiEta(0,0),ic->weight())); phiStates.push_back(SingleGaussianState1D(phi,covPPhiEta(1,1),ic->weight())); etaStates.push_back(SingleGaussianState1D(eta,covPPhiEta(2,2),ic->weight())); } @@ -315,13 +315,13 @@ MultiTrajectoryStateMode::momentumFromModePPhiEta (const TrajectoryStateOnSurfac // // parameters from mode (in case of failure: mean) // - double p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean(); - double phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean(); - double eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean(); + auto p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean(); + auto phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean(); + auto eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean(); // double theta = 2*atan(exp(-eta)); - double tanth2 = exp(-eta); - double pt = p*2*tanth2/(1+tanth2*tanth2); // p*sin(theta) - double pz = p*(1-tanth2*tanth2)/(1+tanth2*tanth2); // p*cos(theta) + auto tanth2 = std::exp(-eta); + auto pt = p*2*tanth2/(1+tanth2*tanth2); // p*sin(theta) + auto pz = p*(1-tanth2*tanth2)/(1+tanth2*tanth2); // p*cos(theta) // conversion to a cartesian momentum vector momentum = GlobalVector(pt*cos(phi),pt*sin(phi),pz); return true; From a571e5c52b0f564c75dcff14b0a836d46db9a2a0 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Sat, 20 Feb 2016 17:22:28 +0100 Subject: [PATCH 03/18] allow vectorization --- .../interface/GsfBetheHeitlerUpdator.h | 14 ++-- TrackingTools/GsfTracking/interface/Triplet.h | 47 ----------- .../GsfTracking/src/GsfBetheHeitlerUpdator.cc | 79 +++++++++++-------- 3 files changed, 50 insertions(+), 90 deletions(-) delete mode 100644 TrackingTools/GsfTracking/interface/Triplet.h diff --git a/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h b/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h index a7a9c604ed6d9..906314076b3c5 100644 --- a/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h +++ b/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h @@ -6,10 +6,6 @@ #include "DataFormats/Math/interface/approx_exp.h" #include "DataFormats/Math/interface/approx_log.h" - - -#include "TrackingTools/GsfTracking/interface/Triplet.h" - #include #include @@ -69,7 +65,7 @@ class GsfBetheHeitlerUpdator final: public GsfMaterialEffectsUpdator { GsfBetheHeitlerUpdator (const std::string fileName, const int correctionFlag); private: - typedef Triplet GSContainer; + struct GSContainer {float *first, *second, *third;}; /// Computation: generates vectors of weights, means and standard deviations virtual void compute (const TrajectoryStateOnSurface&, const PropagationDirection, Effect[]) const; @@ -83,13 +79,13 @@ class GsfBetheHeitlerUpdator final: public GsfMaterialEffectsUpdator { /// Filling of mixture (in terms of z=E/E0) - void getMixtureParameters (const float, GSContainer[]) const; + void getMixtureParameters (const float, GSContainer &) const; /// Correction for weight of component 1 - void correctWeights (GSContainer[]) const; + void correctWeights (GSContainer&) const; /// Correction for mean of component 1 - float correctedFirstMean (const float, const GSContainer[]) const; + float correctedFirstMean (const float, const GSContainer &) const; /// Correction for variance of component 1 - float correctedFirstVar (const float,const GSContainer[]) const; + float correctedFirstVar (const float,const GSContainer &) const; private: diff --git a/TrackingTools/GsfTracking/interface/Triplet.h b/TrackingTools/GsfTracking/interface/Triplet.h deleted file mode 100644 index ee4d889fe9bed..0000000000000 --- a/TrackingTools/GsfTracking/interface/Triplet.h +++ /dev/null @@ -1,47 +0,0 @@ -#ifndef Triplet_H -#define Triplet_H - -/** triplet is identical to stl::pair exceot that it has three members - * instead of two (first, second, third). - */ - -template -struct Triplet { - typedef T1 first_type; - typedef T2 second_type; - typedef T3 third_type; - - T1 first; - T2 second; - T3 third; - Triplet() : first(T1()), second(T2()), third(T3()) {} - Triplet(const T1& a, const T2& b, const T3& c) : - first(a), second(b), third(c) {} - - template - Triplet(const Triplet& p) : - first(p.first), second(p.second), third(p.third) {} -}; - -template -inline bool operator==(const Triplet& x, - const Triplet& y) { - return x.first == y.first && x.second == y.second && x.third==y.third; -} - -template -inline bool operator<(const Triplet& x, - const Triplet& y) { - bool pair_less = - x.first < y.first || (!(y.first < x.first) && x.second < y.second); - - return pair_less || - (!(y.first < x.first) && !(y.second < x.second) && x.third < y.third); -} - -template -inline Triplet make_Triplet(const T1& x, const T2& y, const T3& z) { - return Triplet(x, y, z); -} - -#endif diff --git a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc index 8aa044ca42ea2..1057b6a49be6b 100644 --- a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc @@ -4,12 +4,12 @@ #include "FWCore/ParameterSet/interface/FileInPath.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" + #include #include #include #include - namespace { /// Logistic function (needed for transformation of weight and mean) inline float logisticFunction (const float x) {return 1./(1.+unsafe_expf<4>(-x));} @@ -122,40 +122,37 @@ GsfBetheHeitlerUpdator::compute (const TrajectoryStateOnSurface& TSoS, if ( rl<0.01f ) rl = 0.01f; if ( rl>0.20f ) rl = 0.20f; - #if __clang__ - std::vector mixtureHolder(theNrComponents); - GSContainer *mixture = mixtureHolder.data(); - #else - GSContainer mixture[theNrComponents]; - #endif + float mixtureData[3][theNrComponents]; + GSContainer mixture{mixtureData[0],mixtureData[1],mixtureData[2]}; + getMixtureParameters(rl,mixture); correctWeights(mixture); if ( theCorrectionFlag>=1 ) - mixture[0].second = correctedFirstMean(rl,mixture); + mixture.second[0] = correctedFirstMean(rl,mixture); if ( theCorrectionFlag>=2 ) - mixture[0].third = correctedFirstVar(rl,mixture); + mixture.third[0] = correctedFirstVar(rl,mixture); for ( int i=0; i(vz); + + mixture.first[i]=weight; + mixture.second[i]=z; + mixture.third[i]=vz; + } + else // theTransformationCode for ( int i=0; i(vz); - else - vz = vz*vz; + vz = vz*vz; - mixture[i]=Triplet(weight,z,vz); + mixture.first[i]=weight; + mixture.second[i]=z; + mixture.third[i]=vz; } } @@ -192,57 +203,57 @@ GsfBetheHeitlerUpdator::getMixtureParameters (const float rl, // Correct weights // void -GsfBetheHeitlerUpdator::correctWeights (GSContainer mixture[]) const +GsfBetheHeitlerUpdator::correctWeights (GSContainer & mixture) const { // // get sum of weights // float wsum(0); for ( int i=0; i Date: Sat, 20 Feb 2016 18:23:26 +0100 Subject: [PATCH 04/18] vectorize more --- .../GsfTracking/src/GsfBetheHeitlerUpdator.cc | 39 +++++++------------ 1 file changed, 13 insertions(+), 26 deletions(-) diff --git a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc index 1057b6a49be6b..78bbe8ced1bd7 100644 --- a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc @@ -4,7 +4,6 @@ #include "FWCore/ParameterSet/interface/FileInPath.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" - #include #include #include @@ -166,36 +165,24 @@ void GsfBetheHeitlerUpdator::getMixtureParameters (const float rl, GSContainer & mixture) const { - + + float weight[theNrComponents], z[theNrComponents], vz[theNrComponents]; + for ( int i=0; i(vz); - - mixture.first[i]=weight; - mixture.second[i]=z; - mixture.third[i]=vz; + mixture.first[i]=logisticFunction(weight[i]); + mixture.second[i]=logisticFunction(z[i]); + mixture.third[i]=unsafe_expf<4>(vz[i]);; } else // theTransformationCode for ( int i=0; i Date: Sun, 21 Feb 2016 13:33:08 +0100 Subject: [PATCH 05/18] fix debug --- .../src/GsfMaterialEffectsUpdator.cc | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc b/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc index 46c66658644bd..1f36883af3024 100644 --- a/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc @@ -42,13 +42,15 @@ GsfMaterialEffectsUpdator::updateState (const TrajectoryStateOnSurface& TSoS, // // loop over components // - // edm::LogDebug("GsfMaterialEffectsUpdator") << "found " << Ws.size() << " components\n" - // << " input state has weight " << TSoS.weight(); + LogDebug("GsfMaterialEffectsUpdator") << "found " << size() << " components " + << " input state has weight " << TSoS.weight(); for ( auto const & effect : effects ) { - // edm::LogDebug("GsfMaterialEffectsUpdator") << "w, dp, sigp = " - // << Ws[ic] << ", " - // << dPs[ic] << ", " - // << sqrt((deltaErrors[ic])[0][0]); +#ifdef DEBUG_DETAIL + LogDebug("GsfMaterialEffectsUpdator") << "w, dp, sigp = " + << effect.weight << ", " + << effect.deltaP << ", " + << std::sqrt(effect.deltaCov[materialEffect::elos]); +#endif // // Update momentum. In case of failure: return invalid state. // Use deltaP method to ensure update of cache, if necessary! @@ -68,9 +70,10 @@ GsfMaterialEffectsUpdator::updateState (const TrajectoryStateOnSurface& TSoS, surface, &(TSoS.globalParameters().magneticField()), side)); - - // edm::LogDebug("GsfMaterialEffectsUpdator") - // << "adding state with weight " << weight*Ws[ic]; +#ifdef DEBUG_DETAIL + LogDebug("GsfMaterialEffectsUpdator") + << "adding state with weight " << weight*effect.weight; +#endif } else { result.addState(TrajectoryStateOnSurface(lp,surface, @@ -78,8 +81,8 @@ GsfMaterialEffectsUpdator::updateState (const TrajectoryStateOnSurface& TSoS, side)); } } - // edm::LogDebug("GsfMaterialEffectsUpdator") - // << " output state has weight " << result.combinedState().weight(); + LogDebug("GsfMaterialEffectsUpdator") + << " output state has weight " << result.combinedState().weight(); return result.combinedState(); } From df982864a7826ac390b35f3c8bd632807a2ce9fd Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Sun, 21 Feb 2016 14:24:54 +0100 Subject: [PATCH 06/18] add debug to gsffitter --- .../GsfTracking/src/GsfTrajectoryFitter.cc | 137 +++++++++++++----- 1 file changed, 99 insertions(+), 38 deletions(-) diff --git a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc index 7c8662789b69a..f38a44a445904 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc @@ -2,18 +2,97 @@ #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h" #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" -// #include "CommonDet/BasicDet/interface/Det.h" #include "DataFormats/GeometrySurface/interface/BoundPlane.h" #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h" -// #include "Utilities/Notification/interface/Verbose.h" -// #include "Utilities/Notification/interface/TimingReport.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" - #include "DataFormats/TrackerRecHit2D/interface/TkCloner.h" #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" +#ifdef EDM_ML_DEBUG + +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "DataFormats/MuonDetId/interface/MuonSubdetId.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" + +namespace { + void dump(TrackingRecHit const & hit, int hitcounter) { + if (hit.isValid()) { + LogTrace("GsfTrackFitters")<< " ----------------- HIT #" << hitcounter << " (VALID)-----------------------\n" + << " HIT IS AT R " << hit.globalPosition().perp() << "\n" + << " HIT IS AT Z " << hit.globalPosition().z() << "\n" + << " HIT IS AT Phi " << hit.globalPosition().phi() << "\n" + << " HIT IS AT Loc " << hit.localPosition() << "\n" + << " WITH LocError " << hit.localPositionError() << "\n" + << " HIT IS AT Glo " << hit.globalPosition() << "\n" + << "SURFACE POSITION" << "\n" + << hit.surface()->position()<<"\n" + << "SURFACE ROTATION" << "\n" + << hit.surface()->rotation() + << "dimension " << hit.dimension(); + + DetId hitId = hit.geographicalId(); + + LogDebug("GsfTrackFitters") << " hit det=" << hitId.rawId(); + + if(hitId.det() == DetId::Tracker) { + if (hitId.subdetId() == StripSubdetector::TIB ) + LogDebug("GsfTrackFitters") << " I am TIB " << TIBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TOB ) + LogDebug("GsfTrackFitters") << " I am TOB " << TOBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TEC ) + LogDebug("GsfTrackFitters") << " I am TEC " << TECDetId(hitId).wheel(); + else if (hitId.subdetId() == StripSubdetector::TID ) + LogDebug("GsfTrackFitters") << " I am TID " << TIDDetId(hitId).wheel(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) + LogDebug("GsfTrackFitters") << " I am PixBar " << PXBDetId(hitId).layer(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap ) + LogDebug("GsfTrackFitters") << " I am PixFwd " << PXFDetId(hitId).disk(); + else + LogDebug("GsfTrackFitters") << " UNKNOWN TRACKER HIT TYPE "; + } + else if(hitId.det() == DetId::Muon) { + if(hitId.subdetId() == MuonSubdetId::DT) + LogDebug("GsfTrackFitters") << " I am DT " << DTWireId(hitId); + else if (hitId.subdetId() == MuonSubdetId::CSC ) + LogDebug("GsfTrackFitters") << " I am CSC " << CSCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::RPC ) + LogDebug("GsfTrackFitters") << " I am RPC " << RPCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::GEM ) + LogDebug("GsfTrackFitters") << " I am GEM " << GEMDetId(hitId); + + else if (hitId.subdetId() == MuonSubdetId::ME0 ) + LogDebug("GsfTrackFitters") << " I am ME0 " << ME0DetId(hitId); + else + LogDebug("GsfTrackFitters") << " UNKNOWN MUON HIT TYPE "; + } + else + LogDebug("GsfTrackFitters") << " UNKNOWN HIT TYPE "; + + } else { + LogDebug("GsfTrackFitters") + << " ----------------- INVALID HIT #" << hitcounter << " -----------------------"; + } + } +} +#else +namespace { + inline void dump(TrackingRecHit const &, int) {} +} +#endif + + + GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator, const TrajectoryStateUpdator& aUpdator, const MeasurementEstimator& aEstimator, @@ -26,8 +105,6 @@ GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator, theGeometry(detLayerGeometry) { if(!theGeometry) theGeometry = &dummyGeometry; - // static SimpleConfigurable timeConf(false,"GsfTrajectoryFitter:activateTiming"); - // theTiming = timeConf.value(); } GsfTrajectoryFitter::~GsfTrajectoryFitter() { @@ -60,22 +137,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, const TrajectoryStateOnSurface& firstPredTsos, fitType) const { - // static TimingReport::Item* propTimer = - // &(*TimingReport::current())[string("GsfTrajectoryFitter:propagation")]; - // propTimer->switchCPU(false); - // if ( !theTiming ) propTimer->switchOn(false); - // static TimingReport::Item* updateTimer = - // &(*TimingReport::current())[string("GsfTrajectoryFitter:update")]; - // updateTimer->switchCPU(false); - // if ( !theTiming ) updateTimer->switchOn(false); - if(hits.empty()) return Trajectory(); Trajectory myTraj(aSeed, propagator()->propagationDirection()); TSOS predTsos(firstPredTsos); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitter") << "GsfTrajectoryFitter: predicted tsos of first measurement not valid!"; return Trajectory(); } @@ -87,8 +155,8 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, assert( (!(ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); assert( (!(ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast(ihit.get()))); auto preciseHit = theHitCloner->makeShared(ihit,predTsos); + dump(*preciseHit,1); { - // TimeMe t(*updateTimer,false); currTsos = updator()->update(predTsos, *preciseHit); } if (!predTsos.isValid() || !currTsos.isValid()){ @@ -106,8 +174,11 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, myTraj.push(TM(predTsos, *hits.begin(),0., theGeometry->idToLayer((*hits.begin())->geographicalId()) )); } + int hitcounter = 1; for(RecHitContainer::const_iterator ihit = hits.begin() + 1; - ihit != hits.end(); ihit++) { + ihit != hits.end(); ihit++) { + ++hitcounter; + // // temporary protection copied from KFTrajectoryFitter. // @@ -116,23 +187,6 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, continue; } - //!!! no invalid hits on cylinders anymore?? - // // - // // check type of surface in case of invalid hit - // // (in this version only propagations to planes are - // // supported for multi trajectory states) - // // - // if ( !(**ihit).isValid() ) { - // const BoundPlane* plane = - // dynamic_cast(&(**ihit).det().surface()); - // // - // // no plane: insert invalid measurement - // // - // if ( plane==0 ) { - // myTraj.push(TM(TrajectoryStateOnSurface(),&(**ihit))); - // continue; - // } - // } { // TimeMe t(*propTimer,false); predTsos = propagator()->propagate(currTsos, @@ -140,13 +194,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, } if(!predTsos.isValid()) { if ( myTraj.foundHits()>=3 ) { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitter") << "GsfTrajectoryFitter: predicted tsos not valid! \n" << "Returning trajectory with " << myTraj.foundHits() << " found hits."; return myTraj; } else { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitter") << "GsfTrajectoryFitter: predicted tsos not valid after " << myTraj.foundHits() << " hits, discarding candidate!"; return Trajectory(); @@ -159,8 +213,8 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast((*ihit).get()))); auto preciseHit = theHitCloner->makeShared(*ihit,predTsos); + dump(*preciseHit,hitcounter); { - // TimeMe t(*updateTimer,false); currTsos = updator()->update(predTsos, *preciseHit); } if (!predTsos.isValid() || !currTsos.isValid()){ @@ -178,6 +232,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, } myTraj.push(TM(predTsos, *ihit,0., theGeometry->idToLayer( (*ihit)->geographicalId()) )); } + LogTrace("GsfTrackFitters") + << "predTsos !" << "\n" + << predTsos + <<" with local position " << predTsos.localPosition() + <<"currTsos !" << "\n" + << currTsos + <<" with local position " << currTsos.localPosition(); } return myTraj; } From 0204f01780fb4d3a0b7b095425a7788483aa36f4 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Sun, 21 Feb 2016 16:26:08 +0100 Subject: [PATCH 07/18] more debug --- .../src/MultiTrajectoryStateAssembler.cc | 45 +++++++------------ .../GsfTracking/src/GsfMultiStateUpdator.cc | 12 +++-- 2 files changed, 20 insertions(+), 37 deletions(-) diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index 538dfa85b41c8..a6ff915131ee4 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -17,22 +17,9 @@ MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () : sortStates = false; minValidFraction = 0.01; minFractionalWeight = 1.e-6; - // // - // // Timers - // // - // if ( theTimerAdd==0 ) { - // theTimerAdd = - // &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::addState")]; - // theTimerAdd->switchCPU(false); - // theTimerComb = - // &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::combinedState")]; - // theTimerComb->switchCPU(false); - // } } void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) { - // // Timer - // TimeMe t(*theTimerAdd,false); // // refuse to add states after combination has been done // @@ -92,8 +79,6 @@ void MultiTrajectoryStateAssembler::addInvalidState (const double weight) { } TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () { - // // Timer - // TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // @@ -110,8 +95,6 @@ TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () { } TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState (const float newWeight) { - // // Timer - // TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // @@ -215,11 +198,11 @@ MultiTrajectoryStateAssembler::removeSmallWeights() void MultiTrajectoryStateAssembler::removeWrongPz () { - // edm::LogDebug("MultiTrajectoryStateAssembler") - // << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n" - // << " #state / weights before cleaning = " << theStates.size() - // << " / " << theValidWeightSum - // << " / " << theInvalidWeightSum; + LogDebug("GsfTrackFitters") + << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n" + << " #state / weights before cleaning = " << theStates.size() + << " / " << theValidWeightSum + << " / " << theInvalidWeightSum; // // Calculate average pz // @@ -227,7 +210,7 @@ MultiTrajectoryStateAssembler::removeWrongPz () { for ( MultiTSOS::const_iterator is=theStates.begin(); is!=theStates.end(); is++ ) { meanPz += is->weight()*is->localParameters().pzSign(); - // edm::LogDebug("MultiTrajectoryStateAssembler") + // LogDebug("GsfTrackFitters") // << " weight / pz / global position = " << is->weight() // << " " << is->localParameters().pzSign() // << " " << is->globalPosition(); @@ -248,13 +231,15 @@ MultiTrajectoryStateAssembler::removeWrongPz () { } else { theInvalidWeightSum += is->weight(); + LogDebug("GsfTrackFitters") + << "removing weight / pz / global position = " << is->weight() + << " " << is->localParameters().pzSign() + << " " << is->globalPosition(); + } } - // edm::LogDebug("MultiTrajectoryStateAssembler") - // << " #state / weights after cleaning = " << theStates.size() - // << " / " << theValidWeightSum - // << " / " << theInvalidWeightSum; + LogDebug("GsfTrackFitters") + << " #state / weights after cleaning = " << theStates.size() + << " / " << theValidWeightSum + << " / " << theInvalidWeightSum; } - -// TimingReport::Item * MultiTrajectoryStateAssembler::theTimerAdd(0); -// TimingReport::Item * MultiTrajectoryStateAssembler::theTimerComb(0); diff --git a/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc b/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc index b479d717e7b8d..97060654949b5 100644 --- a/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc @@ -28,22 +28,20 @@ TrajectoryStateOnSurface GsfMultiStateUpdator::update(const TrajectoryStateOnSur MultiTrajectoryStateAssembler result; int i = 0; - for (std::vector::const_iterator iter = predictedComponents.begin(); - iter != predictedComponents.end(); iter++) { - TrajectoryStateOnSurface updatedTSOS = KFUpdator().update(*iter, aRecHit); + for (auto const & tsosI : predictedComponents) { + TrajectoryStateOnSurface updatedTSOS = KFUpdator().update(tsosI, aRecHit); if (updatedTSOS.isValid()){ result.addState(TrajectoryStateOnSurface(weights[i], updatedTSOS.localParameters(), updatedTSOS.localError(), updatedTSOS.surface(), &(tsos.globalParameters().magneticField()), - (*iter).surfaceSide() + tsosI.surfaceSide() )); - i++; } else{ - edm::LogError("GsfMultiStateUpdator") << "one of the KF updated state is invalid. skipping."; + edm::LogError("GsfMultiStateUpdator") << "KF updated state " << i << " is invalid. skipping."; } - + ++i; } return result.combinedState(); From 33522e14b49f1e53012fafaaa0a361417013abf5 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Sun, 21 Feb 2016 16:37:40 +0100 Subject: [PATCH 08/18] add size to debug --- TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc | 2 +- TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index a6ff915131ee4..862ad9143664d 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -198,7 +198,7 @@ MultiTrajectoryStateAssembler::removeSmallWeights() void MultiTrajectoryStateAssembler::removeWrongPz () { - LogDebug("GsfTrackFitters") + LogDebug("GsfTrackFitters") << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n" << " #state / weights before cleaning = " << theStates.size() << " / " << theValidWeightSum diff --git a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc index f38a44a445904..8a908cf6be734 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc @@ -234,10 +234,10 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, } LogTrace("GsfTrackFitters") << "predTsos !" << "\n" - << predTsos + << predTsos << " size " << predTsos.components().size() <<" with local position " << predTsos.localPosition() - <<"currTsos !" << "\n" - << currTsos + <<"\ncurrTsos !" << "\n" + << currTsos << " size " << currTsos.components().size() <<" with local position " << currTsos.localPosition(); } return myTraj; From 7f7a80444914466c644721f0833d00462a97a8e4 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 26 Feb 2016 16:24:16 +0100 Subject: [PATCH 09/18] cleanup, make inv 1dim matrix trivial --- .../Math/interface/invertPosDefMatrix.h | 11 +++ .../interface/MultiGaussianStateAssembler.h | 3 - .../interface/MultiGaussianStateAssembler.icc | 98 +++---------------- .../GsfTools/src/BasicMultiTrajectoryState.cc | 10 +- .../GsfTools/src/MultiStatePropagation.icc | 13 +-- .../src/MultiTrajectoryStateAssembler.cc | 67 +++++-------- .../src/FullConvolutionWithMaterial.cc | 15 +-- .../GsfTracking/src/GsfMultiStateUpdator.cc | 4 +- .../src/GsfPropagatorWithMaterial.cc | 28 ------ .../GsfTracking/src/GsfTrajectoryFitter.cc | 50 ++++++---- .../src/MultiTrajectoryStateMerger.cc | 1 - .../src/PosteriorWeightsCalculator.cc | 31 +++--- .../src/TsosGaussianStateConversions.cc | 35 +++---- TrackingTools/KalmanUpdators/BuildFile.xml | 2 - 14 files changed, 122 insertions(+), 246 deletions(-) diff --git a/DataFormats/Math/interface/invertPosDefMatrix.h b/DataFormats/Math/interface/invertPosDefMatrix.h index 42d53e6d7f745..d6b1ad3cc5af3 100644 --- a/DataFormats/Math/interface/invertPosDefMatrix.h +++ b/DataFormats/Math/interface/invertPosDefMatrix.h @@ -17,6 +17,17 @@ inline bool invertPosDefMatrix(ROOT::Math::SMatrix +inline bool invertPosDefMatrix(ROOT::Math::SMatrix > & m) { + m(0,0) = 1./m(0,0); + return true; +} +template<> +inline bool invertPosDefMatrix(ROOT::Math::SMatrix > & m) { + m(0,0) = 1.f/m(0,0); + return true; +} + template inline bool invertPosDefMatrix(ROOT::Math::SMatrix > const & mIn, ROOT::Math::SMatrix > & mOut) { diff --git a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h index ec29ec4ee86eb..635d238efb513 100644 --- a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h +++ b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h @@ -65,15 +65,12 @@ class MultiGaussianStateAssembler { private: const MultiState theInitialState; -// bool sortStates; double minFractionalWeight; bool combinationDone; double theValidWeightSum; SingleStateContainer theStates; -// static TimingReport::Item * theTimerAdd; -// static TimingReport::Item * theTimerComb; }; diff --git a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc index 5e330a1a8f8aa..3db2235eee0a8 100644 --- a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc +++ b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc @@ -1,30 +1,12 @@ #include "TrackingTools/GsfTools/interface/GaussianStateLessWeight.h" #include "FWCore/Utilities/interface/Exception.h" + template MultiGaussianStateAssembler::MultiGaussianStateAssembler (const MultiState & state) : theInitialState(state), minFractionalWeight(1.e-16), combinationDone(false), theValidWeightSum(0.) -{ - // - // parameters (could be configurable) - // -// static SimpleConfigurable sortStatesConf(false,"MultiGaussianStateAssembler:sortStates"); -// sortStates = sortStatesConf.value(); -// static SimpleConfigurable minWeightConf(1.e-16,"MultiGaussianStateAssembler:minFractionalWeight"); -// minFractionalWeight = minWeightConf.value(); -// // -// // Timers -// // -// if ( theTimerAdd==0 ) { -// theTimerAdd = -// &(*TimingReport::current())[string("MultiGaussianStateAssembler::addState")]; -// theTimerAdd->switchCPU(false); -// theTimerComb = -// &(*TimingReport::current())[string("MultiGaussianStateAssembler::combinedState")]; -// theTimerComb->switchCPU(false); -// } -} +{} template void MultiGaussianStateAssembler::addState (const SingleStatePtr& sgs) @@ -44,8 +26,6 @@ void MultiGaussianStateAssembler::addState (const SingleStatePtr& sgs) template void MultiGaussianStateAssembler::addState (const MultiState& mgs) { -// // Timer -// TimeMe t(*theTimerAdd,false); // // refuse to add states after combination has been done // @@ -53,12 +33,6 @@ void MultiGaussianStateAssembler::addState (const MultiState& mgs) { throw cms::Exception("LogicError") << "MultiGaussianStateAssembler: trying to add states after combination"; // - // Verify validity of state to be added - //ThS: Is not possible for gaussian states. -// if ( !tsos.isValid() ) -// throw cms::Exception("LogicError") -// << "MultiGaussianStateAssembler: trying to add invalid state"; - // // Add components (i.e. state to be added can be single or multi state) // SingleStateContainer components(mgs.components()); @@ -77,15 +51,8 @@ void MultiGaussianStateAssembler::addStateVector (const SingleStateContainer& // // sum up weights (all components are supposed to be valid!!!) // - double sum(0.); - for ( typename SingleStateContainer::const_iterator i=states.begin(); - i!=states.end(); i++ ) { - //ThS: Is not possible for gaussian states: -// if ( !(i->isValid()) ) -// throw cms::Exception("LogicError") -// << "MultiGaussianStateAssembler: trying to add invalid state"; - sum += (**i).weight(); - } + double sum=0.; + for (auto const & is : states ) sum += (*is).weight(); theValidWeightSum += sum; // // add to vector of states @@ -97,31 +64,19 @@ void MultiGaussianStateAssembler::addStateVector (const SingleStateContainer& template MultiGaussianState MultiGaussianStateAssembler::combinedState () { -// // Timer -// TimeMe t(*theTimerComb,false); // // Prepare resulting state vector prepareCombinedState(); - // ThS: What to return here? I choose a RCMP with empty state vector... -// if ( !prepareCombinedState() ) return TSOS(); - // - // Return new multi state without reweighting - // -// return theInitialState.createNewState(theStates); return MultiState(theStates); } template MultiGaussianState MultiGaussianStateAssembler::combinedState (const float newWeight) { -// // Timer -// TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // prepareCombinedState(); - // ThS: What to return here? I choose a RCMP with empty state vector... -// if ( !prepareCombinedState() ) return TSOS(); // // return reweighted state // @@ -135,19 +90,12 @@ MultiGaussianStateAssembler::prepareCombinedState () { // remaining part to be done only once // if ( combinationDone ) return true; - else combinationDone = true; + combinationDone = true; // // Remove states with negligible weights // removeSmallWeights(); - if ( theStates.empty() ) return false; -// // -// // Sort output by weights? (currently deactivated - make it again configurable?) -// // -// if ( sortStates ) -// sort(theStates.begin(),theStates.end(),GaussianStateLessWeight()); - - return true; + return !theStates.empty(); } template @@ -162,16 +110,11 @@ MultiGaussianStateAssembler::reweightedCombinedState (const double newWeight) // SingleStateContainer reweightedStates; reweightedStates.reserve(theStates.size()); - for ( typename SingleStateContainer::const_iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - double oldWeight = (**i).weight(); -// reweightedStates.push_back( (**i).createNewState ((**i).mean(), -// (**i).covariance(), factor*oldWeight) ); - reweightedStates.push_back( SingleState((**i).mean(), - (**i).covariance(), factor*oldWeight) ); + for (auto const & is : theStates) { + auto oldWeight = (*is).weight(); + reweightedStates.emplace_back((*is).mean(), + (*is).covariance(), factor*oldWeight); } - -// return theInitialState.createNewState(reweightedStates); return MultiState(reweightedStates); } @@ -186,22 +129,7 @@ MultiGaussianStateAssembler::removeSmallWeights() theStates.clear(); return; } - // - // Loop until no more states are removed - // - bool redo; - do { - redo = false; - for ( typename SingleStateContainer::iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - if ( (**i).weight()/theValidWeightSum < minFractionalWeight ) { - theStates.erase(i); - redo = true; - break; - } - } - } while (redo); + theStates.erase(std::remove_if(theStates.begin(),theStates.end(), + [&](typename SingleStateContainer::value_type const & s){ return s->weight() < minFractionalWeight*theValidWeightSum;}),theStates.end() + ); } - -// TimingReport::Item * MultiGaussianStateAssembler::theTimerAdd(0); -// TimingReport::Item * MultiGaussianStateAssembler::theTimerComb(0); diff --git a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc index b049dfe66c13d..30e2c39337e6e 100644 --- a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc +++ b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc @@ -46,9 +46,7 @@ void BasicMultiTrajectoryState::rescaleError(double factor) { return; } - for (std::vector::iterator it = theStates.begin(); it != theStates.end(); it++) { - it->rescaleError(factor); - } + for (auto & is : theStates) is.rescaleError(factor); combine(); } @@ -88,13 +86,13 @@ BasicMultiTrajectoryState::combine() { AlgebraicVector5 mean; AlgebraicSymMatrix55 covarPart1, covarPart2, covtmp; for (auto it1 = tsos.begin(); it1 != tsos.end(); it1++) { - double weight = it1->weight(); - AlgebraicVector5 param = it1->localParameters().vector(); + auto weight = it1->weight(); + auto const & param = it1->localParameters().vector(); sumw += weight; mean += weight * param; covarPart1 += weight * it1->localError().matrix(); for (auto it2 = it1 + 1; it2 != tsos.end(); it2++) { - AlgebraicVector5 diff = param - it2->localParameters().vector(); + auto diff = param - it2->localParameters().vector(); ROOT::Math::AssignSym::Evaluate(covtmp,ROOT::Math::TensorProd(diff,diff)); covarPart2 += (weight * it2->weight()) * covtmp; } diff --git a/TrackingTools/GsfTools/src/MultiStatePropagation.icc b/TrackingTools/GsfTools/src/MultiStatePropagation.icc index e4e0c81a7add2..774a4a65b48cd 100644 --- a/TrackingTools/GsfTools/src/MultiStatePropagation.icc +++ b/TrackingTools/GsfTools/src/MultiStatePropagation.icc @@ -19,19 +19,19 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso // bool firstPropagation(true); SurfaceSideDefinition::SurfaceSide firstSide(SurfaceSideDefinition::atCenterOfSurface); - for ( MultiTSOS::const_iterator iTsos=input.begin(); - iTsos!=input.end(); iTsos++ ) { + for ( auto const & iTsos : input) { // // weight of component // - double weight(iTsos->weight()); + double weight(iTsos.weight()); // // geometrical propagation (assumption: only one output state!) // TsosWP newTsosWP = - thePropagator.propagateWithPath(*iTsos,surface); + thePropagator.propagateWithPath(iTsos,surface); // check validity if ( !(newTsosWP.first).isValid() ) { + LogDebug("MultiStatePropagation") << "adding invalid state"; result.addInvalidState(weight); continue; } @@ -77,14 +77,9 @@ template TrajectoryStateOnSurface MultiStatePropagation::setWeight (const TrajectoryStateOnSurface state, const double weight) const { - // if ( state.hasError() ) return TrajectoryStateOnSurface(weight, state.localParameters(), state.localError(), state.surface(), &(state.globalParameters().magneticField()), state.surfaceSide()); - // else - // return TrajectoryStateOnSurface(weight, state.localParameters(), - // state.surface()); - // return state; } diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index 862ad9143664d..c3021002da4b0 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -124,7 +124,7 @@ MultiTrajectoryStateAssembler::prepareCombinedState () { // remaining part to be done only once // if ( combinationDone ) return true; - else combinationDone = true; + combinationDone = true; // // Remove states with negligible weights // @@ -154,16 +154,15 @@ MultiTrajectoryStateAssembler::reweightedCombinedState (const double newWeight) // MultiTSOS reweightedStates; reweightedStates.reserve(theStates.size()); - for ( MultiTSOS::const_iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - double oldWeight = i->weight(); - reweightedStates.push_back(TrajectoryStateOnSurface(factor*oldWeight, - i->localParameters(), - i->localError(), - i->surface(), - &(i->globalParameters().magneticField()), - i->surfaceSide() - )); + for (auto const & is : theStates) { + auto oldWeight = is.weight(); + reweightedStates.emplace_back(factor*oldWeight, + is.localParameters(), + is.localError(), + is.surface(), + &(is.globalParameters().magneticField()), + is.surfaceSide() + ); } return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(reweightedStates))); } @@ -179,21 +178,9 @@ MultiTrajectoryStateAssembler::removeSmallWeights() theStates.clear(); return; } - // - // Loop until no more states are removed - // - bool redo; - do { - redo = false; - for ( MultiTSOS::iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - if ( (*i).weight() < minFractionalWeight*totalWeight ) { - theStates.erase(i); - redo = true; - break; - } - } - } while (redo); + theStates.erase(std::remove_if(theStates.begin(),theStates.end(), + [&](MultiTSOS::value_type const & s){ return s.weight() < minFractionalWeight*totalWeight;}), + theStates.end()); } void @@ -207,34 +194,26 @@ MultiTrajectoryStateAssembler::removeWrongPz () { // Calculate average pz // double meanPz(0.); - for ( MultiTSOS::const_iterator is=theStates.begin(); - is!=theStates.end(); is++ ) { - meanPz += is->weight()*is->localParameters().pzSign(); - // LogDebug("GsfTrackFitters") - // << " weight / pz / global position = " << is->weight() - // << " " << is->localParameters().pzSign() - // << " " << is->globalPosition(); - } + for (auto const & is : theStates) + meanPz += is.weight()*is.localParameters().pzSign(); meanPz /= theValidWeightSum; // // Now keep only states compatible with the average pz // - // double oldValidWeight(theValidWeightSum); theValidWeightSum = 0.; MultiTSOS oldStates(theStates); theStates.clear(); - for ( MultiTSOS::const_iterator is=oldStates.begin(); - is!=oldStates.end(); is++ ) { - if ( meanPz*is->localParameters().pzSign()>=0. ) { - theValidWeightSum += is->weight(); - theStates.push_back(*is); + for (auto const & is :oldStates) { + if ( meanPz*is.localParameters().pzSign()>=0. ) { + theValidWeightSum += is.weight(); + theStates.push_back(std::move(is)); } else { - theInvalidWeightSum += is->weight(); + theInvalidWeightSum += is.weight(); LogDebug("GsfTrackFitters") - << "removing weight / pz / global position = " << is->weight() - << " " << is->localParameters().pzSign() - << " " << is->globalPosition(); + << "removing weight / pz / global position = " << is.weight() + << " " << is.localParameters().pzSign() + << " " << is.globalPosition(); } } diff --git a/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc b/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc index f3a8d50d9e9fc..14bc37ab03914 100644 --- a/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc +++ b/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc @@ -1,5 +1,4 @@ #include "TrackingTools/GsfTracking/interface/FullConvolutionWithMaterial.h" - #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateAssembler.h" TrajectoryStateOnSurface @@ -8,8 +7,7 @@ FullConvolutionWithMaterial::operator() (const TrajectoryStateOnSurface& tsos, // // decomposition of input state // - typedef std::vector MultiTSOS; - MultiTSOS input(tsos.components()); + auto && input = tsos.components(); // // vector of result states // @@ -17,20 +15,15 @@ FullConvolutionWithMaterial::operator() (const TrajectoryStateOnSurface& tsos, // // now add material effects to each component // - for ( MultiTSOS::const_iterator iTsos=input.begin(); - iTsos!=input.end(); iTsos++ ) { - // + for (auto const & iTsos : input) { // add material // TrajectoryStateOnSurface updatedTSoS = - theMEUpdator->updateState(*iTsos,propDir); + theMEUpdator->updateState(iTsos,propDir); if ( updatedTSoS.isValid() ) result.addState(updatedTSoS); else - result.addInvalidState(iTsos->weight()); + result.addInvalidState(iTsos.weight()); } return result.combinedState(); } - -// TimingReport::Item* FullConvolutionWithMaterial::theTimer1(0); -// TimingReport::Item* FullConvolutionWithMaterial::theTimer2(0); diff --git a/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc b/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc index 97060654949b5..17236f3775d7c 100644 --- a/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfMultiStateUpdator.cc @@ -13,13 +13,13 @@ TrajectoryStateOnSurface GsfMultiStateUpdator::update(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& aRecHit) const { - std::vector && predictedComponents = tsos.components(); + auto && predictedComponents = tsos.components(); if (predictedComponents.empty()) { edm::LogError("GsfMultiStateUpdator") << "Trying to update trajectory state with zero components! " ; return TrajectoryStateOnSurface(); } - std::vector && weights = PosteriorWeightsCalculator(predictedComponents).weights(aRecHit); + auto && weights = PosteriorWeightsCalculator(predictedComponents).weights(aRecHit); if ( weights.empty() ) { edm::LogError("GsfMultiStateUpdator") << " no weights could be retreived. invalid updated state !."; return TrajectoryStateOnSurface(); diff --git a/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc b/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc index 9995ad824cd84..0189cf74e246d 100644 --- a/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc +++ b/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc @@ -15,7 +15,6 @@ GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const Propagator& aPropaga theConvolutor(new FullConvolutionWithMaterial(aMEUpdator)), theMaterialLocation(atDestination) { - // if ( propWithPathTimer1==0 ) defineTimer(); } GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const GsfPropagatorAdapter& aGsfPropagator, @@ -25,34 +24,11 @@ GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const GsfPropagatorAdapter theConvolutor(aConvolutor.clone()), theMaterialLocation(atDestination) { - // if ( propWithPathTimer1==0 ) defineTimer(); } -// void -// GsfPropagatorWithMaterial::defineTimer() -// { -// if ( propWithPathTimer1==0 ) { -// propWithPathTimer1 = -// &(*TimingReport::current())[string("GsfPropagatorWithMaterial:toPlane")]; -// propWithPathTimer2 = -// &(*TimingReport::current())[string("GsfPropagatorWithMaterial:toCylinder")]; -// static SimpleConfigurable timeConf(false,"GsfPropagatorWithMaterial:activateTiming"); -// if ( timeConf.value() ) { -// propWithPathTimer1->switchCPU(false); -// propWithPathTimer2->switchCPU(false); -// } -// else { -// propWithPathTimer1->switchOn(false); -// propWithPathTimer2->switchOn(false); -// } -// } -// } - std::pair GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& tsos, const Plane& plane) const { - // TimeMe t1(*propWithPathTimer1,false); - // // add material before propagation? // TrajectoryStateOnSurface stateAtSource; @@ -75,8 +51,6 @@ GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& ts std::pair GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& tsos, const Cylinder& cylinder) const { - // TimeMe t2(*propWithPathTimer2,false); - // // add material before propagation? // TrajectoryStateOnSurface stateAtSource; @@ -167,5 +141,3 @@ GsfPropagatorWithMaterial::materialAtSource() const { propagationDirection()==alongMomentum); } -// TimingReport::Item* GsfPropagatorWithMaterial::propWithPathTimer1(0); -// TimingReport::Item* GsfPropagatorWithMaterial::propWithPathTimer2(0); diff --git a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc index 8a908cf6be734..11d121da5a523 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc @@ -85,9 +85,22 @@ namespace { } } } +#include + inline void dump(TrajectoryStateOnSurface const & tsos, const char * header) { + std::ostringstream ss; ss<< " weights "; + for (auto const & c : tsos.components()) ss << c.weight() << '/'; + ss << "\nmomentums "; + for (auto const & c : tsos.components()) ss << c.globalMomentum().mag() << '/'; + LogTrace("GsfTrackFitters") + << header << "! size " << tsos.components().size() << ss.str() << "\n" + <<" with local position " << tsos.localPosition() << "\n" + << tsos; + } + #else namespace { inline void dump(TrackingRecHit const &, int) {} + inline void dump(TrajectoryStateOnSurface const &, const char *){} } #endif @@ -171,7 +184,7 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, edm::LogError("InvalidState")<<"first invalid hit"; return Trajectory(); } - myTraj.push(TM(predTsos, *hits.begin(),0., theGeometry->idToLayer((*hits.begin())->geographicalId()) )); + myTraj.push(TM(predTsos, hits.front(),0., theGeometry->idToLayer((hits.front())->geographicalId()) )); } int hitcounter = 1; @@ -183,7 +196,7 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, // temporary protection copied from KFTrajectoryFitter. // if ((**ihit).isValid() == false && (**ihit).det() == 0) { - LogDebug("GsfTrajectoryFitter") << " Error: invalid hit with no GeomDet attached .... skipping"; + LogDebug("GsfTrackFitters") << " Error: invalid hit with no GeomDet attached .... skipping"; continue; } @@ -194,13 +207,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, } if(!predTsos.isValid()) { if ( myTraj.foundHits()>=3 ) { - edm::LogInfo("GsfTrackFitter") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectoryFitter: predicted tsos not valid! \n" << "Returning trajectory with " << myTraj.foundHits() << " found hits."; return myTraj; } else { - edm::LogInfo("GsfTrackFitter") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectoryFitter: predicted tsos not valid after " << myTraj.foundHits() << " hits, discarding candidate!"; return Trajectory(); @@ -210,35 +223,30 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, if((**ihit).isValid()) { //update - assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); - assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast((*ihit).get()))); - auto preciseHit = theHitCloner->makeShared(*ihit,predTsos); - dump(*preciseHit,hitcounter); - { - currTsos = updator()->update(predTsos, *preciseHit); - } + assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); + assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast((*ihit).get()))); + auto preciseHit = theHitCloner->makeShared(*ihit,predTsos); + dump(*preciseHit,hitcounter); + currTsos = updator()->update(predTsos, *preciseHit); if (!predTsos.isValid() || !currTsos.isValid()){ edm::LogError("InvalidState")<<"inside hit"; return Trajectory(); } + auto chi2=estimator()->estimate(predTsos, *preciseHit).second; myTraj.push(TM(predTsos, currTsos, preciseHit, - estimator()->estimate(predTsos, *preciseHit).second, + chi2, theGeometry->idToLayer(preciseHit->geographicalId() ))); + LogDebug("GsfTrackFitters") << "added measurement with chi2 " << chi2; } else { currTsos = predTsos; if (!predTsos.isValid()){ - edm::LogError("InvalidState")<<"inside invalid hit"; - return Trajectory(); + edm::LogError("InvalidState")<<"inside invalid hit"; + return Trajectory(); } myTraj.push(TM(predTsos, *ihit,0., theGeometry->idToLayer( (*ihit)->geographicalId()) )); } - LogTrace("GsfTrackFitters") - << "predTsos !" << "\n" - << predTsos << " size " << predTsos.components().size() - <<" with local position " << predTsos.localPosition() - <<"\ncurrTsos !" << "\n" - << currTsos << " size " << currTsos.components().size() - <<" with local position " << currTsos.localPosition(); + dump(predTsos,"predTsos"); + dump(currTsos,"currTsos"); } return myTraj; } diff --git a/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc b/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc index 24652c09661ae..e79353a437763 100644 --- a/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc +++ b/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc @@ -1,5 +1,4 @@ #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h" - #include "TrackingTools/GsfTracking/interface/TsosGaussianStateConversions.h" TrajectoryStateOnSurface diff --git a/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc b/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc index d57f9d86374f8..59bcbacfce910 100644 --- a/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc +++ b/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc @@ -3,6 +3,8 @@ #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h" +#include "DataFormats/Math/interface/invertPosDefMatrix.h" +#include "DataFormats/Math/interface/ProjectMatrix.h" #include @@ -39,7 +41,6 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re VecD r, rMeas; SMatDD V(SMatrixNoInit{}), R(SMatrixNoInit{}); - AlgebraicVector5 x; ProjectMatrix p; // // calculate chi2 and determinant / component and find @@ -49,7 +50,7 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re for ( unsigned int i=0; i(&r, &V, &p, &rMeas, &R, x, predictedComponents[i].localError().matrix()); recHit.getKfComponents(holder); @@ -64,21 +65,23 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re } detRs.push_back(detR); - int ierr = ! R.Invert(); // if (ierr != 0) throw exception; - if ( ierr!=0 ) { + bool ok = invertPosDefMatrix(R); + if ( !ok ) { edm::LogError("PosteriorWeightsCalculator") - << "PosteriorWeightsCalculator: inversion failed, ierr = " << ierr; + << "PosteriorWeightsCalculator: inversion failed"; return std::vector(); } double chi2 = ROOT::Math::Similarity(r,R); chi2s.push_back(chi2); if ( chi2(); } + // // calculate weights (extracting a common factor // exp(-0.5*chi2Min) to avoid numerical problems @@ -96,16 +99,17 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re // Calculation of (non-normalised) weight. Common factors exp(-chi2Norm/2.) and // 1./sqrt(2*pi*recHit.dimension()) have been omitted // - tempWeight = priorWeight * sqrt(1./detRs[i]) * exp(-0.5 * chi2); + tempWeight = priorWeight * std::sqrt(1./detRs[i]) * std::exp(-0.5 * chi2); + } + else { + LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: detR < FLT_MIN !!"; } - // else { - // edm::LogInfo("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: detR < FLT_MIN !!"; - // } weights.push_back(tempWeight); sumWeights += tempWeight; } + if ( sumWeights(); } @@ -114,10 +118,7 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes (2)"; return std::vector(); } - for (std::vector::iterator iter = weights.begin(); - iter != weights.end(); iter++) { - (*iter) /= sumWeights; - } - + sumWeights = 1./sumWeights; + for (auto & w : weights) w *= sumWeights; return weights; } diff --git a/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc b/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc index d11adafc3185f..f1c5373f8131f 100644 --- a/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc +++ b/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc @@ -12,16 +12,15 @@ namespace GaussianStateConversions { { if ( !tsos.isValid() ) return MultiGaussianState<5>(); - typedef boost::shared_ptr< SingleGaussianState<5> > SingleStatePtr; - const std::vector& components = tsos.components(); + using SingleStatePtr = boost::shared_ptr>; + auto const & components = tsos.components(); MultiGaussianState<5>::SingleStateContainer singleStates; singleStates.reserve(components.size()); - for ( std::vector::const_iterator ic=components.begin(); - ic!=components.end(); ic ++ ) { - if ( ic->isValid() ) { - SingleStatePtr sgs(new SingleGaussianState<5>(ic->localParameters().vector(), - ic->localError().matrix(), - ic->weight())); + for (auto const & ic : components) { + if ( ic.isValid() ) { + SingleStatePtr sgs(new SingleGaussianState<5>(ic.localParameters().vector(), + ic.localError().matrix(), + ic.weight())); singleStates.push_back(sgs); } } @@ -32,23 +31,21 @@ namespace GaussianStateConversions { const TrajectoryStateOnSurface refTsos) { if ( multiState.components().empty() ) return TrajectoryStateOnSurface(); - const Surface& surface = refTsos.surface(); + const Surface & surface = refTsos.surface(); SurfaceSide side = refTsos.surfaceSide(); const MagneticField* field = refTsos.magneticField(); - TrajectoryStateOnSurface refTsos1 = refTsos.components().front(); - double pzSign = refTsos1.localParameters().pzSign(); + auto const & refTsos1 = refTsos.components().front(); + auto pzSign = refTsos1.localParameters().pzSign(); bool charged = refTsos1.charge()!=0; - const MultiGaussianState<5>::SingleStateContainer& singleStates = - multiState.components(); + auto const & singleStates = multiState.components(); std::vector components; components.reserve(singleStates.size()); - for ( MultiGaussianState<5>::SingleStateContainer::const_iterator ic=singleStates.begin(); - ic!=singleStates.end(); ic++ ) { - components.push_back(TrajectoryStateOnSurface((**ic).weight(), LocalTrajectoryParameters((**ic).mean(), - pzSign,charged), - LocalTrajectoryError((**ic).covariance()), - surface,field,side)); + for ( auto const & ic : singleStates ) { + components.emplace_back((*ic).weight(), + LocalTrajectoryParameters((*ic).mean(), pzSign,charged), + LocalTrajectoryError((*ic).covariance()), + surface,field,side); } return TrajectoryStateOnSurface((BasicTrajectoryState*)new BasicMultiTrajectoryState(components)); } diff --git a/TrackingTools/KalmanUpdators/BuildFile.xml b/TrackingTools/KalmanUpdators/BuildFile.xml index 3e9d8f8035384..4fed031aa2979 100644 --- a/TrackingTools/KalmanUpdators/BuildFile.xml +++ b/TrackingTools/KalmanUpdators/BuildFile.xml @@ -1,6 +1,4 @@ - - From 56a48e2b679a1c955196d24cc424d1b68e8d709a Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 26 Feb 2016 16:37:17 +0100 Subject: [PATCH 10/18] fix auto for samtrix --- TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc index 30e2c39337e6e..a6560787b4c69 100644 --- a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc +++ b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc @@ -92,7 +92,7 @@ BasicMultiTrajectoryState::combine() { mean += weight * param; covarPart1 += weight * it1->localError().matrix(); for (auto it2 = it1 + 1; it2 != tsos.end(); it2++) { - auto diff = param - it2->localParameters().vector(); + AlgebraicVector5 diff = param - it2->localParameters().vector(); ROOT::Math::AssignSym::Evaluate(covtmp,ROOT::Math::TensorProd(diff,diff)); covarPart2 += (weight * it2->weight()) * covtmp; } From 040e5c8a4e35a09d073f3a12e30fed28aab439c7 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 26 Feb 2016 17:25:30 +0100 Subject: [PATCH 11/18] stabilize gsf lowering minFractionalWeight --- TrackingTools/GsfTools/src/MultiStatePropagation.icc | 6 +++--- TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/TrackingTools/GsfTools/src/MultiStatePropagation.icc b/TrackingTools/GsfTools/src/MultiStatePropagation.icc index 774a4a65b48cd..eb5eabd3bae51 100644 --- a/TrackingTools/GsfTools/src/MultiStatePropagation.icc +++ b/TrackingTools/GsfTools/src/MultiStatePropagation.icc @@ -11,7 +11,7 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso double sWeight(0.); double sWeightPath(0.); // decomposition of input state - MultiTSOS input(tsos.components()); + MultiTSOS const & input = tsos.components(); // vector of result states MultiTrajectoryStateAssembler result; // @@ -31,7 +31,7 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso thePropagator.propagateWithPath(iTsos,surface); // check validity if ( !(newTsosWP.first).isValid() ) { - LogDebug("MultiStatePropagation") << "adding invalid state"; + LogDebug("GsfTrackFitters") << "adding invalid state"; result.addInvalidState(weight); continue; } @@ -44,7 +44,7 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso } else { if ( weightedTsos.surfaceSide()!=firstSide ) { - edm::LogInfo("MultiStatePropagation") + edm::LogInfo("GsfTrackFitters") << "MultiStatePropagation resulted in states " << "on different sides of the target surface"; return TsosWP(TrajectoryStateOnSurface(),0.); diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index c3021002da4b0..cf97338b1f1d0 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -16,7 +16,7 @@ MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () : // sortStates = false; minValidFraction = 0.01; - minFractionalWeight = 1.e-6; + minFractionalWeight = 1.e-4; // 6; } void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) { From 26963acda0ffada50515d5896d162812add8a722 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 09:59:37 +0100 Subject: [PATCH 12/18] fix few denug compilation isseus --- RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc | 9 ++++++--- RecoParticleFlow/PFTracking/plugins/PFV0Producer.cc | 2 +- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc index 8e9b221cd06eb..5c66a4dda7710 100644 --- a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc +++ b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc @@ -200,7 +200,7 @@ GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) auto const & Tj=*(tjCollection.product()); LogDebug("GoodSeedProducer")<<"Number of tracks in collection " - < V0coll; iEvent.getByToken(V0list_[il],V0coll); - LogDebug("PFV0Producer")<size()<<" V0 candidates "; + LogDebug("PFV0Producer")<< "V0list_[" << il <<"] contains "<size()<<" V0 candidates "; for (unsigned int iv=0;ivsize();iv++){ VertexCompositeCandidateRef V0(V0coll, iv); vector Tracks; From b128eba22ccc1c68177f89ea8e70236e502d469a Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 10:01:13 +0100 Subject: [PATCH 13/18] roll back tuning --- TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index cf97338b1f1d0..4ef4f1d1a5fc8 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -16,7 +16,7 @@ MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () : // sortStates = false; minValidFraction = 0.01; - minFractionalWeight = 1.e-4; // 6; + minFractionalWeight = 1.e-6; // 4; } void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) { From 83ac89fa4d2cb3997a8cb622952926cd1d723d48 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 11:36:43 +0100 Subject: [PATCH 14/18] add debug to smoother --- .../src/TrackProducerAlgorithm.cc | 58 +++++------ TrackingTools/GsfTracking/src/DebugHelpers.h | 96 +++++++++++++++++++ .../GsfTracking/src/GsfTrajectoryFitter.cc | 96 +------------------ .../GsfTracking/src/GsfTrajectorySmoother.cc | 87 +++++++---------- 4 files changed, 162 insertions(+), 175 deletions(-) create mode 100644 TrackingTools/GsfTracking/src/DebugHelpers.h diff --git a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc index 5c1b8cebfd178..1ba931e2a0454 100644 --- a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc +++ b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc @@ -30,6 +30,8 @@ #include "TrackingTools/TrackFitters/interface/RecHitSorter.h" #include "DataFormats/TrackReco/interface/TrackBase.h" +#include + // #define VI_DEBUG #define STAT_TSB @@ -257,33 +259,36 @@ TrackProducerAlgorithm::buildTrack (const TrajectoryFitter * the auto theTraj = new Trajectory( std::move(trajTmp) ); theTraj->setSeedRef(seedRef); - - // TrajectoryStateOnSurface innertsos; - // TrajectoryStateOnSurface outertsos; - // if (theTraj->direction() == alongMomentum) { - // innertsos = theTraj->firstMeasurement().updatedState(); - // outertsos = theTraj->lastMeasurement().updatedState(); - // } else { - // innertsos = theTraj->lastMeasurement().updatedState(); - // outertsos = theTraj->firstMeasurement().updatedState(); - // } - // std::cout - // << "Nr. of first / last states = " - // << innertsos.components().size() << " " - // << outertsos.components().size() << std::endl; - // std::vector components = - // innertsos.components(); - // double sinTheta = - // sin(innertsos.globalMomentum().theta()); - // for ( std::vector::const_iterator ic=components.begin(); - // ic!=components.end(); ic++ ) { - // std::cout << " comp " << ic-components.begin() << " " - // << (*ic).weight() << " " - // << (*ic).localParameters().vector()[0]/sinTheta << " " - // << sqrt((*ic).localError().matrix()[0][0])/sinTheta << std::endl; - // } - +#ifdef EDM_ML_DEBUG + TrajectoryStateOnSurface innertsos; + TrajectoryStateOnSurface outertsos; + + if (theTraj->direction() == alongMomentum) { + innertsos = theTraj->firstMeasurement().updatedState(); + outertsos = theTraj->lastMeasurement().updatedState(); + } else { + innertsos = theTraj->lastMeasurement().updatedState(); + outertsos = theTraj->firstMeasurement().updatedState(); + } + std::ostringstream ss; + auto dc = [&](TrajectoryStateOnSurface const & tsos){ + std::vector const & components = tsos.components(); + auto sinTheta = std::sin(tsos.globalMomentum().theta()); + for (auto const & ic : components) ss << ic.weight() << "/"; ss << "\n"; + for (auto const & ic : components) ss << ic.localParameters().vector()[0]/sinTheta << "/"; ss << "\n"; + for (auto const & ic : components) ss << std::sqrt(ic.localError().matrix()(0,0))/sinTheta << "/"; + }; + ss << "\ninner comps\n"; + dc(innertsos); + ss << "\nouter comps\n"; + dc(outertsos); + LogDebug("TrackProducer") + << "Nr. of first / last states = " + << innertsos.components().size() << " " + << outertsos.components().size() << ss.str(); +#endif + ndof = 0; for (auto const & tm : theTraj->measurements()) { auto const & h = tm.recHitR(); @@ -342,7 +347,6 @@ TrackProducerAlgorithm::buildTrack (const TrajectoryFitter * the math::XYZVector mom( p.x(), p.y(), p.z() ); LogDebug("GsfTrackProducer") << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag(); - std::cout << "GsfTrackProducer " << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag() << std::endl; auto theTrack = new reco::GsfTrack(theTraj->chiSquared(), int(ndof),//FIXME fix weight() in TrackingRecHit diff --git a/TrackingTools/GsfTracking/src/DebugHelpers.h b/TrackingTools/GsfTracking/src/DebugHelpers.h new file mode 100644 index 0000000000000..c88aee50f0076 --- /dev/null +++ b/TrackingTools/GsfTracking/src/DebugHelpers.h @@ -0,0 +1,96 @@ +#ifdef EDM_ML_DEBUG + +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "DataFormats/MuonDetId/interface/MuonSubdetId.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" + +namespace { + inline + void dump(TrackingRecHit const & hit, int hitcounter) { + if (hit.isValid()) { + LogTrace("GsfTrackFitters")<< " ----------------- HIT #" << hitcounter << " (VALID)-----------------------\n" + << " HIT IS AT R " << hit.globalPosition().perp() << "\n" + << " HIT IS AT Z " << hit.globalPosition().z() << "\n" + << " HIT IS AT Phi " << hit.globalPosition().phi() << "\n" + << " HIT IS AT Loc " << hit.localPosition() << "\n" + << " WITH LocError " << hit.localPositionError() << "\n" + << " HIT IS AT Glo " << hit.globalPosition() << "\n" + << "SURFACE POSITION" << "\n" + << hit.surface()->position()<<"\n" + << "SURFACE ROTATION" << "\n" + << hit.surface()->rotation() + << "dimension " << hit.dimension(); + + DetId hitId = hit.geographicalId(); + + LogDebug("GsfTrackFitters") << " hit det=" << hitId.rawId(); + + if(hitId.det() == DetId::Tracker) { + if (hitId.subdetId() == StripSubdetector::TIB ) + LogDebug("GsfTrackFitters") << " I am TIB " << TIBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TOB ) + LogDebug("GsfTrackFitters") << " I am TOB " << TOBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TEC ) + LogDebug("GsfTrackFitters") << " I am TEC " << TECDetId(hitId).wheel(); + else if (hitId.subdetId() == StripSubdetector::TID ) + LogDebug("GsfTrackFitters") << " I am TID " << TIDDetId(hitId).wheel(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) + LogDebug("GsfTrackFitters") << " I am PixBar " << PXBDetId(hitId).layer(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap ) + LogDebug("GsfTrackFitters") << " I am PixFwd " << PXFDetId(hitId).disk(); + else + LogDebug("GsfTrackFitters") << " UNKNOWN TRACKER HIT TYPE "; + } + else if(hitId.det() == DetId::Muon) { + if(hitId.subdetId() == MuonSubdetId::DT) + LogDebug("GsfTrackFitters") << " I am DT " << DTWireId(hitId); + else if (hitId.subdetId() == MuonSubdetId::CSC ) + LogDebug("GsfTrackFitters") << " I am CSC " << CSCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::RPC ) + LogDebug("GsfTrackFitters") << " I am RPC " << RPCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::GEM ) + LogDebug("GsfTrackFitters") << " I am GEM " << GEMDetId(hitId); + + else if (hitId.subdetId() == MuonSubdetId::ME0 ) + LogDebug("GsfTrackFitters") << " I am ME0 " << ME0DetId(hitId); + else + LogDebug("GsfTrackFitters") << " UNKNOWN MUON HIT TYPE "; + } + else + LogDebug("GsfTrackFitters") << " UNKNOWN HIT TYPE "; + + } else { + LogDebug("GsfTrackFitters") + << " ----------------- INVALID HIT #" << hitcounter << " -----------------------"; + } + } +} +#include + inline void dump(TrajectoryStateOnSurface const & tsos, const char * header) { + std::ostringstream ss; ss<< " weights "; + for (auto const & c : tsos.components()) ss << c.weight() << '/'; + ss << "\nmomentums "; + for (auto const & c : tsos.components()) ss << c.globalMomentum().mag() << '/'; + ss << "\ndeltap/p "; + for (auto const & c : tsos.components()) ss << std::sqrt(tsos.curvilinearError().matrix()(0,0))/c.globalMomentum().mag() << '/'; + LogTrace("GsfTrackFitters") + << header << "! size " << tsos.components().size() << ss.str() << "\n" + <<" with local position " << tsos.localPosition() << "\n" + << tsos; + } +#else +namespace { + inline void dump(TrackingRecHit const &, int) {} + inline void dump(TrajectoryStateOnSurface const &, const char *){} +} +#endif + diff --git a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc index 11d121da5a523..7631757f19251 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc @@ -9,101 +9,7 @@ #include "DataFormats/TrackerRecHit2D/interface/TkCloner.h" #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" - -#ifdef EDM_ML_DEBUG - -#include "DataFormats/MuonDetId/interface/CSCDetId.h" -#include "DataFormats/MuonDetId/interface/DTWireId.h" -#include "DataFormats/MuonDetId/interface/RPCDetId.h" -#include "DataFormats/MuonDetId/interface/GEMDetId.h" -#include "DataFormats/MuonDetId/interface/ME0DetId.h" -#include "DataFormats/MuonDetId/interface/MuonSubdetId.h" -#include "DataFormats/SiStripDetId/interface/TIBDetId.h" -#include "DataFormats/SiStripDetId/interface/TOBDetId.h" -#include "DataFormats/SiStripDetId/interface/TIDDetId.h" -#include "DataFormats/SiStripDetId/interface/TECDetId.h" -#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" -#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" - -namespace { - void dump(TrackingRecHit const & hit, int hitcounter) { - if (hit.isValid()) { - LogTrace("GsfTrackFitters")<< " ----------------- HIT #" << hitcounter << " (VALID)-----------------------\n" - << " HIT IS AT R " << hit.globalPosition().perp() << "\n" - << " HIT IS AT Z " << hit.globalPosition().z() << "\n" - << " HIT IS AT Phi " << hit.globalPosition().phi() << "\n" - << " HIT IS AT Loc " << hit.localPosition() << "\n" - << " WITH LocError " << hit.localPositionError() << "\n" - << " HIT IS AT Glo " << hit.globalPosition() << "\n" - << "SURFACE POSITION" << "\n" - << hit.surface()->position()<<"\n" - << "SURFACE ROTATION" << "\n" - << hit.surface()->rotation() - << "dimension " << hit.dimension(); - - DetId hitId = hit.geographicalId(); - - LogDebug("GsfTrackFitters") << " hit det=" << hitId.rawId(); - - if(hitId.det() == DetId::Tracker) { - if (hitId.subdetId() == StripSubdetector::TIB ) - LogDebug("GsfTrackFitters") << " I am TIB " << TIBDetId(hitId).layer(); - else if (hitId.subdetId() == StripSubdetector::TOB ) - LogDebug("GsfTrackFitters") << " I am TOB " << TOBDetId(hitId).layer(); - else if (hitId.subdetId() == StripSubdetector::TEC ) - LogDebug("GsfTrackFitters") << " I am TEC " << TECDetId(hitId).wheel(); - else if (hitId.subdetId() == StripSubdetector::TID ) - LogDebug("GsfTrackFitters") << " I am TID " << TIDDetId(hitId).wheel(); - else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) - LogDebug("GsfTrackFitters") << " I am PixBar " << PXBDetId(hitId).layer(); - else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap ) - LogDebug("GsfTrackFitters") << " I am PixFwd " << PXFDetId(hitId).disk(); - else - LogDebug("GsfTrackFitters") << " UNKNOWN TRACKER HIT TYPE "; - } - else if(hitId.det() == DetId::Muon) { - if(hitId.subdetId() == MuonSubdetId::DT) - LogDebug("GsfTrackFitters") << " I am DT " << DTWireId(hitId); - else if (hitId.subdetId() == MuonSubdetId::CSC ) - LogDebug("GsfTrackFitters") << " I am CSC " << CSCDetId(hitId); - else if (hitId.subdetId() == MuonSubdetId::RPC ) - LogDebug("GsfTrackFitters") << " I am RPC " << RPCDetId(hitId); - else if (hitId.subdetId() == MuonSubdetId::GEM ) - LogDebug("GsfTrackFitters") << " I am GEM " << GEMDetId(hitId); - - else if (hitId.subdetId() == MuonSubdetId::ME0 ) - LogDebug("GsfTrackFitters") << " I am ME0 " << ME0DetId(hitId); - else - LogDebug("GsfTrackFitters") << " UNKNOWN MUON HIT TYPE "; - } - else - LogDebug("GsfTrackFitters") << " UNKNOWN HIT TYPE "; - - } else { - LogDebug("GsfTrackFitters") - << " ----------------- INVALID HIT #" << hitcounter << " -----------------------"; - } - } -} -#include - inline void dump(TrajectoryStateOnSurface const & tsos, const char * header) { - std::ostringstream ss; ss<< " weights "; - for (auto const & c : tsos.components()) ss << c.weight() << '/'; - ss << "\nmomentums "; - for (auto const & c : tsos.components()) ss << c.globalMomentum().mag() << '/'; - LogTrace("GsfTrackFitters") - << header << "! size " << tsos.components().size() << ss.str() << "\n" - <<" with local position " << tsos.localPosition() << "\n" - << tsos; - } - -#else -namespace { - inline void dump(TrackingRecHit const &, int) {} - inline void dump(TrajectoryStateOnSurface const &, const char *){} -} -#endif - +#include "DebugHelpers.h" GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator, diff --git a/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc b/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc index 6bfca3a78c00e..07c5226873100 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc @@ -6,6 +6,9 @@ #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DebugHelpers.h" + + GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aPropagator, const TrajectoryStateUpdator& aUpdator, const MeasurementEstimator& aEstimator, @@ -37,8 +40,6 @@ GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aP if(!theGeometry) theGeometry = &dummyGeometry; - // static SimpleConfigurable timeConf(false,"GsfTrajectorySmoother:activateTiming"); - // theTiming = timeConf.value(); } GsfTrajectorySmoother::~GsfTrajectorySmoother() { @@ -54,15 +55,6 @@ GsfTrajectorySmoother::~GsfTrajectorySmoother() { Trajectory GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { - // static TimingReport::Item* propTimer = - // &(*TimingReport::current())[string("GsfTrajectorySmoother:propagation")]; - // propTimer->switchCPU(false); - // if ( !theTiming ) propTimer->switchOn(false); - // static TimingReport::Item* updateTimer = - // &(*TimingReport::current())[string("GsfTrajectorySmoother:update")]; - // updateTimer->switchCPU(false); - // if ( !theTiming ) updateTimer->switchOn(false); - if(aTraj.empty()) return Trajectory(); const Propagator* usePropagator = theAlongPropagator; @@ -81,7 +73,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { predTsos.rescaleError(theErrorRescaling); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos of last measurement not valid!"; return Trajectory(); } @@ -128,49 +120,30 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { TrajectoryStateCombiner combiner; + int hitcounter = avtm.size()-1; for(std::vector::const_reverse_iterator itm = avtm.rbegin() + 1; itm < avtm.rend() - 1; ++itm) { - { - // TimeMe t(*propTimer,false); - // // - // // check type of surface in case of invalid hit - // // (in this version only propagations to planes are - // // supported for multi trajectory states) - // // - // if ( !(*itm).recHit().isValid() ) { - // const BoundPlane* plane = - // dynamic_cast(&(*itm).recHit().det().surface()); - // // - // // no plane: insert invalid - // if ( plane==0 ) { - // myTraj.push(TM(TrajectoryStateOnSurface(), - // (*itm).recHit()); - // continue; - // } - // } - predTsos = usePropagator->propagate(currTsos, + + predTsos = usePropagator->propagate(currTsos, *(*itm).recHit()->surface()); - } if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate ) predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection()); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!"; + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!"; return Trajectory(); } if ( theMerger ) predTsos = theMerger->merge(predTsos); if((*itm).recHit()->isValid()) { //update - { - // TimeMe t(*updateTimer,false); currTsos = updator()->update(predTsos, *(*itm).recHit()); - } + if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate ) currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection()); if(!currTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!"; return Trajectory(); } @@ -180,7 +153,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { //3: combine bwd-prediction with fwd-filter TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState()); if(!combTsos.isValid()) { - LogDebug("GsfTrajectorySmoother") << + LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!\n"<< "pred Tsos pos: "<estimate(combTsos, *(*itm).recHit()).second; myTraj.push(TM((*itm).forwardPredictedState(), predTsos, smooTsos, (*itm).recHit(), - estimator()->estimate(combTsos, *(*itm).recHit()).second, + chi2, theGeometry->idToLayer((*itm).recHit()->geographicalId() ) ), (*itm).estimate()); + LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- <<" with chi2 " << chi2; + dump(predTsos,"predTsos"); + dump(smooTsos,"smooTsos"); + } else { currTsos = predTsos; TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState()); if(!combTsos.isValid()) { - LogDebug("GsfTrajectorySmoother") << + LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!"; return Trajectory(); } @@ -231,36 +208,35 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { (*itm).recHit(), 0., theGeometry->idToLayer((*itm).recHit()->geographicalId()) )); + LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--; + dump(predTsos,"predTsos"); + dump(combTsos,"smooTsos"); } if ( theMerger ) currTsos = theMerger->merge(currTsos); + + dump(currTsos,"currTsos"); } //last smoothed tm is last filtered - { - // TimeMe t(*propTimer,false); predTsos = usePropagator->propagate(currTsos, *avtm.front().recHit()->surface()); - } if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate ) predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection()); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!"; + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!"; return Trajectory(); } if ( theMerger ) predTsos = theMerger->merge(predTsos); if(avtm.front().recHit()->isValid()) { //update - { - // TimeMe t(*updateTimer,false); - currTsos = updator()->update(predTsos, *avtm.front().recHit()); - } + currTsos = updator()->update(predTsos, *avtm.front().recHit()); if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate ) currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection()); if(!currTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!"; return Trajectory(); } @@ -270,14 +246,18 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { return Trajectory(); } + auto chi2 = estimator()->estimate(predTsos, *avtm.front().recHit()).second; myTraj.push(TM(avtm.front().forwardPredictedState(), predTsos, currTsos, avtm.front().recHit(), - estimator()->estimate(predTsos, *avtm.front().recHit()).second, + chi2, theGeometry->idToLayer(avtm.front().recHit()->geographicalId() )), avtm.front().estimate()); - //estimator()->estimate(predTsos, avtm.front().recHit())); + LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- <<" with chi2 " << chi2; + dump(predTsos,"predTsos"); + dump(currTsos,"smooTsos"); + } else { if (!avtm.front().forwardPredictedState().isValid()){ @@ -288,6 +268,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { avtm.front().recHit(), 0., theGeometry->idToLayer(avtm.front().recHit()->geographicalId()) )); + LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--; } return myTraj; From 336115ca8b8eee0c36e64a6116fd003a58098ae7 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 11:46:23 +0100 Subject: [PATCH 15/18] add examples --- RecoTracker/TrackProducer/test/TrackRefit.py | 99 +++++---------- .../TrackProducer/test/TrackRefitOld.py | 80 ++++++++++++ RecoTracker/TrackProducer/test/recoDebug.py | 120 ++++++++++++++++++ 3 files changed, 230 insertions(+), 69 deletions(-) create mode 100644 RecoTracker/TrackProducer/test/TrackRefitOld.py create mode 100644 RecoTracker/TrackProducer/test/recoDebug.py diff --git a/RecoTracker/TrackProducer/test/TrackRefit.py b/RecoTracker/TrackProducer/test/TrackRefit.py index 64eab2d5c8f37..009ca4aac4e57 100644 --- a/RecoTracker/TrackProducer/test/TrackRefit.py +++ b/RecoTracker/TrackProducer/test/TrackRefit.py @@ -1,80 +1,41 @@ import FWCore.ParameterSet.Config as cms - -process = cms.Process("Refitting") -### standard MessageLoggerConfiguration -process.load("FWCore.MessageService.MessageLogger_cfi") +process = cms.Process("Refitting") ### Standard Configurations -process.load("Configuration.StandardSequences.Services_cff") -process.load('Configuration/StandardSequences/GeometryIdeal_cff') -process.load('Configuration/StandardSequences/Reconstruction_cff') -process.load('Configuration/StandardSequences/MagneticField_AutoFromDBCurrent_cff') - - -## Fitter-smoother: loosen outlier rejection as for first data-taking with LHC "collisions" -process.KFFittingSmootherWithOutliersRejectionAndRK.BreakTrajWith2ConsecutiveMissing = False -process.KFFittingSmootherWithOutliersRejectionAndRK.EstimateCut = 1000 +#process.load("RecoVertex.BeamSpotProducer.BeamSpot_cfi") +#process.load("Configuration.StandardSequences.MagneticField_cff") +#process.load('Configuration.Geometry.GeometryRecoDB_cff') +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") + +process.load('Configuration.StandardSequences.Services_cff') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') + +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +# choose! +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_data_GRun', '') +# process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc_GRun', '') - -### Conditions -process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") -#process.GlobalTag.globaltag = "IDEAL_V5::All" -process.GlobalTag.globaltag = 'GR09_P_V6::All' - ### Track refitter specific stuff -process.load("RecoTracker.TrackProducer.TrackRefitters_cff") #the correct one -#process.load("RecoTracker.TrackProducer.RefitterWithMaterial_cff") #the one for backward compatibility - - -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) -) - -process.source = cms.Source("PoolSource", -### tracks from collisions - fileNames = cms.untracked.vstring( -'rfio:/castor/cern.ch/user/c/chiochia/09_beam_commissioning/BSCskim_123151_Express.root') -#'/store/relval/CMSSW_2_1_10/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v1/0001/1E04FC31-F99A-DD11-94EE-0018F3D096DE.root') - -### tracks from cosmics -# fileNames = cms.untracked.vstring( -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/006F3A6A-0373-DD11-A8E7-00304876A0FF.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02CF5B1E-6476-DD11-A034-003048769E65.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02DF31C3-A775-DD11-91C2-001A92971BB8.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02F71F56-CE74-DD11-9DD0-001A92810AE4.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0446C89C-E072-DD11-A341-0018F3D0960C.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04750FC3-3E73-DD11-B054-00304876A147.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04DFD531-0473-DD11-964E-0018F3D096AE.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067111FB-3873-DD11-AD86-00304875A9C5.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067982F4-E175-DD11-99F7-001731AF6AC5.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0680EB9B-4F73-DD11-83F8-0018F3D0962E.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/06BF1AF3-E175-DD11-B467-00304876A147.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A3843F3-E175-DD11-8419-003048767EE7.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A5AAABA-3973-DD11-B949-003048767FA1.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A911B18-0273-DD11-A5A6-001731A283E1.root') - -### tracks from beam halo muons -) - -process.TRACKS = cms.OutputModule("PoolOutputModule", - outputCommands = cms.untracked.vstring('drop *_*_*_*', - 'keep recoTracks_*_*_*', - 'keep recoTrackExtras_*_*_*', - 'keep TrackingRecHitsOwned_*_*_*'), - - fileName = cms.untracked.string('refitting.root') - ) +process.load("RecoTracker.TrackProducer.TrackRefitter_cfi") +process.TrackRefitter.NavigationSchool = '' +process.TrackRefitter.Fitter = 'FlexibleKFFittingSmoother' -process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) +process.source = cms.Source ("PoolSource", + fileNames=cms.untracked.vstring('file:pickevents_1.root', + ), + skipEvents=cms.untracked.uint32(0) + ) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(-1)) -process.p1 = cms.Path(process.TrackRefitter - #process.TrackRefitterP5 - #process.TrackRefitterBHM -) -process.outpath = cms.EndPath(process.TRACKS) +process.Path = cms.Path(process.TrackRefitter) - diff --git a/RecoTracker/TrackProducer/test/TrackRefitOld.py b/RecoTracker/TrackProducer/test/TrackRefitOld.py new file mode 100644 index 0000000000000..64eab2d5c8f37 --- /dev/null +++ b/RecoTracker/TrackProducer/test/TrackRefitOld.py @@ -0,0 +1,80 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("Refitting") + +### standard MessageLoggerConfiguration +process.load("FWCore.MessageService.MessageLogger_cfi") + +### Standard Configurations +process.load("Configuration.StandardSequences.Services_cff") +process.load('Configuration/StandardSequences/GeometryIdeal_cff') +process.load('Configuration/StandardSequences/Reconstruction_cff') +process.load('Configuration/StandardSequences/MagneticField_AutoFromDBCurrent_cff') + + +## Fitter-smoother: loosen outlier rejection as for first data-taking with LHC "collisions" +process.KFFittingSmootherWithOutliersRejectionAndRK.BreakTrajWith2ConsecutiveMissing = False +process.KFFittingSmootherWithOutliersRejectionAndRK.EstimateCut = 1000 + + + +### Conditions +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#process.GlobalTag.globaltag = "IDEAL_V5::All" +process.GlobalTag.globaltag = 'GR09_P_V6::All' + +### Track refitter specific stuff +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") #the correct one +#process.load("RecoTracker.TrackProducer.RefitterWithMaterial_cff") #the one for backward compatibility + + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) +) + +process.source = cms.Source("PoolSource", +### tracks from collisions + fileNames = cms.untracked.vstring( +'rfio:/castor/cern.ch/user/c/chiochia/09_beam_commissioning/BSCskim_123151_Express.root') +#'/store/relval/CMSSW_2_1_10/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v1/0001/1E04FC31-F99A-DD11-94EE-0018F3D096DE.root') + +### tracks from cosmics +# fileNames = cms.untracked.vstring( +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/006F3A6A-0373-DD11-A8E7-00304876A0FF.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02CF5B1E-6476-DD11-A034-003048769E65.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02DF31C3-A775-DD11-91C2-001A92971BB8.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02F71F56-CE74-DD11-9DD0-001A92810AE4.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0446C89C-E072-DD11-A341-0018F3D0960C.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04750FC3-3E73-DD11-B054-00304876A147.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04DFD531-0473-DD11-964E-0018F3D096AE.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067111FB-3873-DD11-AD86-00304875A9C5.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067982F4-E175-DD11-99F7-001731AF6AC5.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0680EB9B-4F73-DD11-83F8-0018F3D0962E.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/06BF1AF3-E175-DD11-B467-00304876A147.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A3843F3-E175-DD11-8419-003048767EE7.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A5AAABA-3973-DD11-B949-003048767FA1.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A911B18-0273-DD11-A5A6-001731A283E1.root') + +### tracks from beam halo muons +) + +process.TRACKS = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *_*_*_*', + 'keep recoTracks_*_*_*', + 'keep recoTrackExtras_*_*_*', + 'keep TrackingRecHitsOwned_*_*_*'), + + fileName = cms.untracked.string('refitting.root') + ) + +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + + +process.p1 = cms.Path(process.TrackRefitter + #process.TrackRefitterP5 + #process.TrackRefitterBHM +) +process.outpath = cms.EndPath(process.TRACKS) + + diff --git a/RecoTracker/TrackProducer/test/recoDebug.py b/RecoTracker/TrackProducer/test/recoDebug.py new file mode 100644 index 0000000000000..035f13e282e3e --- /dev/null +++ b/RecoTracker/TrackProducer/test/recoDebug.py @@ -0,0 +1,120 @@ +# Auto generated configuration file +# using: +# Revision: 1.19 +# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v +# with command line options: RECO --conditions 76X_dataRun2_v1 -s RAW2DIGI,L1Reco,RECO,EI,PAT,DQM --runUnscheduled --data --eventcontent RECO,AOD,MINIAOD,DQM --scenario pp --datatier RECO,AOD,MINIAOD,DQMIO --customise Configuration/DataProcessing/RecoTLR.customiseDataRun2Common_25ns --filein /store/data/Run2015C/SingleElectron/RAW/v1/000/254/879/00000/8E51CA98-7349-E511-B9AE-02163E01427B.root -n 100 --no_exec --python_filename=RECO_Run2015C.py +import FWCore.ParameterSet.Config as cms + +process = cms.Process('RECO') + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') +process.load('Configuration.StandardSequences.RawToDigi_Data_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_Data_cff') +process.load('CommonTools.ParticleFlow.EITopPAG_cff') +process.load('DQMOffline.Configuration.DQMOffline_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +#Adding SimpleMemoryCheck service: +#process.SimpleMemoryCheck=cms.Service("SimpleMemoryCheck", +# ignoreTotal=cms.untracked.int32(1), +# oncePerEventMode=cms.untracked.bool(True)) + + +process.Timing = cms.Service("Timing" +# ,summaryOnly = cms.untracked.bool(True) +) + +# process.add_(cms.Service("Tracer")) + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(100) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('file:jetHT256630_RAW.root'), + secondaryFileNames = cms.untracked.vstring() +) + +process.options = cms.untracked.PSet( +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('RECO nevts:100'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('RECO_RAW2DIGI_L1Reco_RECO_EI_PAT_DQM.root'), + outputCommands = process.RECOEventContent.outputCommands, + splitLevel = cms.untracked.int32(0) +) + +# Additional output definition + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, + 'auto:run2_data_GRun', '') # if mc change in mc.... + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.reconstruction_step = cms.Path(process.reconstruction) +process.eventinterpretaion_step = cms.Path(process.EIsequence) +process.RECOoutput_step = cms.EndPath(process.RECOoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.reconstruction_step,process.eventinterpretaion_step,process.RECOoutput_step) + +# customisation of the process. + +# Automatic addition of the customisation function from Configuration.DataProcessing.RecoTLR +from Configuration.DataProcessing.RecoTLR import customiseDataRun2Common_25ns + +#call to customisation function customiseDataRun2Common_25ns imported from Configuration.DataProcessing.RecoTLR +process = customiseDataRun2Common_25ns(process) + +# End of customisation functions + +process.MessageLogger = cms.Service("MessageLogger", + destinations = cms.untracked.vstring("cout"), #1 + debugModules = cms.untracked.vstring("electronGsfTracks"),#initialStepTrackCandidates), #2 + categories = cms.untracked.vstring( + 'TrackProducer', + 'GsfTrackFitters', + # 'GsfMaterialEffectsUpdator', + 'AnalyticalPropagator','RKPropagatorInS', + # 'CkfPattern' + ), #3 + cout = cms.untracked.PSet(threshold = cms.untracked.string("DEBUG"), #4 + DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(0)), #5 + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), #6 + TrackProducer = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + GsfTrackFitters = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + # GsfMaterialEffectsUpdator = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + AnalyticalPropagator = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + RKPropagatorInS = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + # CkfPattern = cms.untracked.PSet(limit = cms.untracked.int32(-1)) #7 + ) + ) + + From 8cf84d412f000617e20688e941cdf5285c849f34 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 11:47:53 +0100 Subject: [PATCH 16/18] remove stat --- RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc index 1ba931e2a0454..29698615027a9 100644 --- a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc +++ b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc @@ -33,7 +33,7 @@ #include // #define VI_DEBUG -#define STAT_TSB +// #define STAT_TSB #ifdef VI_DEBUG #define DPRINT(x) std::cout << x << ": " From 8a41a66632dcd53bb77512a8fc2468c556a6a22b Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 12:24:35 +0100 Subject: [PATCH 17/18] remove ifdef --- .../GsfTracking/src/GsfMaterialEffectsUpdator.cc | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc b/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc index 1f36883af3024..9b4464520c2bf 100644 --- a/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfMaterialEffectsUpdator.cc @@ -45,12 +45,10 @@ GsfMaterialEffectsUpdator::updateState (const TrajectoryStateOnSurface& TSoS, LogDebug("GsfMaterialEffectsUpdator") << "found " << size() << " components " << " input state has weight " << TSoS.weight(); for ( auto const & effect : effects ) { -#ifdef DEBUG_DETAIL - LogDebug("GsfMaterialEffectsUpdator") << "w, dp, sigp = " + LogDebug("GsfMaterialEffectsUpdatorDETAIL") << "w, dp, sigp = " << effect.weight << ", " << effect.deltaP << ", " << std::sqrt(effect.deltaCov[materialEffect::elos]); -#endif // // Update momentum. In case of failure: return invalid state. // Use deltaP method to ensure update of cache, if necessary! @@ -70,10 +68,8 @@ GsfMaterialEffectsUpdator::updateState (const TrajectoryStateOnSurface& TSoS, surface, &(TSoS.globalParameters().magneticField()), side)); -#ifdef DEBUG_DETAIL - LogDebug("GsfMaterialEffectsUpdator") + LogDebug("GsfMaterialEffectsUpdatorDETAIL") << "adding state with weight " << weight*effect.weight; -#endif } else { result.addState(TrajectoryStateOnSurface(lp,surface, From d99ca5b25522b6797100db29be4b664c43502423 Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Fri, 4 Mar 2016 16:11:26 +0100 Subject: [PATCH 18/18] silence clang --- .../GsfTracking/src/GsfBetheHeitlerUpdator.cc | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc index 78bbe8ced1bd7..50bf93c418f48 100644 --- a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc @@ -7,6 +7,7 @@ #include #include #include + #include namespace { @@ -19,7 +20,12 @@ namespace { /// Second moment of the Bethe-Heitler distribution (in z=E/E0) inline float BetheHeitlerVariance (const float rl) { - constexpr float l3ol2 = std::log(3.)/std::log(2.); +#if __clang__ + const +#else + constexpr +#endif + float l3ol2 = std::log(3.)/std::log(2.); float mean = BetheHeitlerMean(rl); return unsafe_expf<4>(-rl*l3ol2) - mean*mean; } @@ -36,7 +42,12 @@ namespace { /// Second moment of the Bethe-Heitler distribution (in z=E/E0) inline float BetheHeitlerVariance (const float rl) { - constexpr float l3ol2 = std::log(3.)/std::log(2.); +#if __clang__ + const +#else + constexpr +#endif + float l3ol2 = std::log(3.)/std::log(2.); return std::exp(-rl*l3ol2) - std::exp(-2*rl); } }