Tp sim hit assoc fastsim compliant

CommonTools/RecoAlgos/interface/CosmicTrackingParticleSelector.h

@@ -4,8 +4,8 @@
  *
  * \author Yanyan Gao, FNAL
  *
- *  $Date: 2013/05/14 14:26:32 $
- *  $Revision: 1.4.4.1 $
+ *  $Date: 2013/06/24 12:25:14 $
+ *  $Revision: 1.5 $
  *
  */
 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
@@ -27,6 +27,9 @@
 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
 
+#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
 class TrajectoryStateClosestToBeamLineBuilder;
 
 class CosmicTrackingParticleSelector {
@@ -62,22 +65,26 @@ class CosmicTrackingParticleSelector {
 	event.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot); 
 	for( TrackingParticleCollection::const_iterator itp = c->begin(); 
 	     itp != c->end(); ++ itp )
-	  if ( operator()(*itp,beamSpot.product(),event,setup) ) {
+	  if ( operator()(TrackingParticleRef(c,itp-c->begin()),beamSpot.product(),event,setup) ) {
 	    selected_.push_back( & * itp );
 	  }
       }
 
     const_iterator begin() const { return selected_.begin(); }
     const_iterator end() const { return selected_.end(); }
-
+
+    void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) const {
+      simHitsTPAssoc = simHitsTPAssocToSet;
+    }
+
     // Operator() performs the selection: e.g. if (tPSelector(tp, bs, event, evtsetup)) {...
-    bool operator()( const TrackingParticle & tp, const reco::BeamSpot* bs, const edm::Event &iEvent, const edm::EventSetup& iSetup ) const { 
-      if (chargedOnly_ && tp.charge()==0) return false;//select only if charge!=0
+    bool operator()( const TrackingParticleRef tpr, const reco::BeamSpot* bs, const edm::Event &iEvent, const edm::EventSetup& iSetup ) const { 
+      if (chargedOnly_ && tpr->charge()==0) return false;//select only if charge!=0
       //bool testId = false;
       //unsigned int idSize = pdgId_.size();
       //if (idSize==0) testId = true;
       //else for (unsigned int it=0;it!=idSize;++it){
-	//if (tp.pdgId()==pdgId_[it]) testId = true;
+	//if (tpr->pdgId()==pdgId_[it]) testId = true;
       //}
 
       edm::ESHandle<TrackerGeometry> tracker;
@@ -90,16 +97,20 @@ class CosmicTrackingParticleSelector {
       GlobalPoint finalGP(0,0,0);      
       GlobalVector momentum(0,0,0);//At the PCA
       GlobalPoint vertex(0,0,0);//At the PCA
-#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
-#ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
       double radius(9999);
-#endif
       bool found(0);
-
-#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
-#ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
-      const std::vector<PSimHit> & simHits = tp.trackPSimHit(DetId::Tracker);
-      for(std::vector<PSimHit>::const_iterator it=simHits.begin(); it!=simHits.end(); ++it){
+
+
+      if (simHitsTPAssoc.isValid()==0) {
+	edm::LogError("TrackAssociation") << "Invalid handle!";
+	return false;
+      }
+      std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(tpr,TrackPSimHitRef());//SimHit is dummy: for simHitTPAssociationListGreater 
+                                                                                                      // sorting only the cluster is needed
+      auto range = std::equal_range(simHitsTPAssoc->begin(), simHitsTPAssoc->end(), 
+				    clusterTPpairWithDummyTP, SimHitTPAssociationProducer::simHitTPAssociationListGreater);
+      for(auto ip = range.first; ip != range.second; ++ip) {
+	TrackPSimHitRef it = ip->second;
 	const GeomDet* tmpDet  = tracker->idToDet( DetId(it->detUnitId()) ) ;
 	LocalVector  lv = it->momentumAtEntry();
 	Local3DPoint lp = it->localPosition ();
@@ -112,11 +123,10 @@ class CosmicTrackingParticleSelector {
 	  finalGP = gp;
 	}
       }
-#endif
       if(!found) return 0;
       else
 	{
-	  FreeTrajectoryState ftsAtProduction(finalGP,finalGV,TrackCharge(tp.charge()),theMF.product());
+	  FreeTrajectoryState ftsAtProduction(finalGP,finalGV,TrackCharge(tpr->charge()),theMF.product());
 	  TSCBLBuilderNoMaterial tscblBuilder;
 	  TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction,*bs);//as in TrackProducerAlgorithm
 	  if(!tsAtClosestApproach.isValid()){
@@ -128,7 +138,7 @@ class CosmicTrackingParticleSelector {
 	      momentum = tsAtClosestApproach.trackStateAtPCA().momentum();
 	      vertex = tsAtClosestApproach.trackStateAtPCA().position();
 	      return (
-		      tp.numberOfTrackerLayers() >= minHit_ &&
+		      tpr->numberOfTrackerLayers() >= minHit_ &&
 		      sqrt(momentum.perp2()) >= ptMin_ && 
 		      momentum.eta() >= minRapidity_ && momentum.eta() <= maxRapidity_ && 
 		      sqrt(vertex.perp2()) <= tip_ &&
@@ -151,7 +161,8 @@ class CosmicTrackingParticleSelector {
     bool chargedOnly_;
       std::vector<int> pdgId_;
       container selected_;
-
+
+    mutable edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;      
 
 };
 

FastSimulation/Validation/python/globalValidation_cff.py

@@ -3,6 +3,12 @@
 # Tracking particle module
 #from FastSimulation.Validation.trackingParticlesFastSim_cfi import * # now deprecated
 
+# TrackingParticle-SimHit associator
+from SimGeneral.TrackingAnalysis.simHitTPAssociation_cfi import * 
+simHitTPAssocProducer.simHitSrc = cms.VInputTag(cms.InputTag('famosSimHits','TrackerHits'),
+                                                cms.InputTag("MuonSimHits","MuonCSCHits"),
+                                                cms.InputTag("MuonSimHits","MuonDTHits"),
+                                                cms.InputTag("MuonSimHits","MuonRPCHits"))
 
 from Validation.RecoMET.METRelValForDQM_cff import *
 
@@ -21,7 +27,10 @@
 
 
 #globalAssociation = cms.Sequence(trackingParticles + recoMuonAssociationFastSim + tracksValidationSelectors + prebTagSequence)
-globalAssociation = cms.Sequence(recoMuonAssociationFastSim + tracksValidationSelectors + prebTagSequence)
+globalAssociation = cms.Sequence(recoMuonAssociationFastSim
+                                 + simHitTPAssocProducer
+                                 + tracksValidationSelectors
+                                 + prebTagSequence)
 
 globalValidation = cms.Sequence(trackingTruthValid
                                 +tracksValidationFS

SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h

@@ -0,0 +1,34 @@
+#ifndef SimGeneral_TrackingAnalysis_SimHitTPAssociationProducer_h
+#define SimGeneral_TrackingAnalysis_SimHitTPAssociationProducer_h
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
+
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
+#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
+
+class SimHitTPAssociationProducer : public edm::EDProducer 
+{
+public:
+
+  typedef std::pair<TrackingParticleRef, TrackPSimHitRef> SimHitTPPair;
+  typedef std::vector<SimHitTPPair> SimHitTPAssociationList;
+
+  explicit SimHitTPAssociationProducer(const edm::ParameterSet&);
+  ~SimHitTPAssociationProducer();
+
+  static bool simHitTPAssociationListGreater(SimHitTPPair i,SimHitTPPair j) { return (i.first.key()>j.first.key()); }
+
+private:
+  virtual void beginJob() {}
+  virtual void produce(edm::Event&, const edm::EventSetup&);
+  virtual void endJob() {}
+
+  std::vector<edm::InputTag> _simHitSrc;
+  edm::InputTag _trackingParticleSrc;
+};
+#endif

SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc

@@ -0,0 +1,68 @@
+#include <memory>
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <utility>
+
+#include "FWCore/Utilities/interface/InputTag.h"
+
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/Common/interface/Handle.h"
+
+#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
+
+SimHitTPAssociationProducer::SimHitTPAssociationProducer(const edm::ParameterSet & cfg) 
+  : _simHitSrc(cfg.getParameter<std::vector<edm::InputTag> >("simHitSrc")),
+    _trackingParticleSrc(cfg.getParameter<edm::InputTag>("trackingParticleSrc"))
+{
+  produces<SimHitTPAssociationList>();
+}
+
+SimHitTPAssociationProducer::~SimHitTPAssociationProducer() {
+}
+
+void SimHitTPAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& es) {
+  std::auto_ptr<SimHitTPAssociationList> simHitTPList(new SimHitTPAssociationList);
+
+  // TrackingParticle
+  edm::Handle<TrackingParticleCollection>  TPCollectionH;
+  iEvent.getByLabel(_trackingParticleSrc,  TPCollectionH);
+
+  // prepare temporary map between SimTrackId and TrackingParticle index
+  std::map<std::pair<size_t, EncodedEventId>, TrackingParticleRef> mapping;
+  for (TrackingParticleCollection::size_type itp = 0; itp < TPCollectionH.product()->size(); ++itp) {
+    TrackingParticleRef trackingParticle(TPCollectionH, itp);
+    // SimTracks inside TrackingParticle
+    EncodedEventId eid(trackingParticle->eventId());
+    for (auto itrk  = trackingParticle->g4Track_begin(); itrk != trackingParticle->g4Track_end(); ++itrk) {
+      std::pair<uint32_t, EncodedEventId> trkid(itrk->trackId(), eid);
+      mapping.insert(std::make_pair(trkid, trackingParticle));
+    }
+  }
+
+  // PSimHits
+  for (auto psit=_simHitSrc.begin();psit<_simHitSrc.end();++psit) {
+    edm::Handle<edm::PSimHitContainer>  PSimHitCollectionH;
+    iEvent.getByLabel(*psit,  PSimHitCollectionH);
+    for (unsigned int psimHit = 0;psimHit != PSimHitCollectionH->size();++psimHit) {
+      TrackPSimHitRef pSimHitRef(PSimHitCollectionH,psimHit);
+      std::pair<uint32_t, EncodedEventId> simTkIds(pSimHitRef->trackId(),pSimHitRef->eventId()); 
+      auto ipos = mapping.find(simTkIds);
+      if (ipos != mapping.end()) {
+	simHitTPList->push_back(std::make_pair(ipos->second,pSimHitRef));
+      }
+    }
+  } 
+
+  std::sort(simHitTPList->begin(),simHitTPList->end(),simHitTPAssociationListGreater);
+  iEvent.put(simHitTPList);
+
+}
+
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+DEFINE_FWK_MODULE(SimHitTPAssociationProducer);

SimGeneral/TrackingAnalysis/python/simHitTPAssociation_cfi.py

@@ -0,0 +1,20 @@
+import FWCore.ParameterSet.Config as cms
+
+simHitTPAssocProducer = cms.EDProducer("SimHitTPAssociationProducer",
+  simHitSrc = cms.VInputTag(cms.InputTag('g4SimHits','MuonDTHits'),
+                            cms.InputTag('g4SimHits','MuonCSCHits'),
+                            cms.InputTag('g4SimHits','MuonRPCHits'),
+                            cms.InputTag('g4SimHits','TrackerHitsTIBLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTIBHighTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTIDLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTIDHighTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTOBLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTOBHighTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTECLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsTECHighTof'),
+                            cms.InputTag( 'g4SimHits','TrackerHitsPixelBarrelLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsPixelBarrelHighTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsPixelEndcapLowTof'),
+                            cms.InputTag('g4SimHits','TrackerHitsPixelEndcapHighTof') ),
+  trackingParticleSrc = cms.InputTag('mix', 'MergedTrackTruth')
+)

SimGeneral/TrackingAnalysis/src/classes.h

@@ -0,0 +1,13 @@
+#include "DataFormats/Common/interface/Wrapper.h"
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
+#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
+
+namespace {
+  struct dictionary {
+    std::pair<TrackingParticleRef, TrackPSimHitRef> dummy13;
+    edm::Wrapper<std::pair<TrackingParticleRef, TrackPSimHitRef> > dummy14;
+    std::vector<std::pair<TrackingParticleRef, TrackPSimHitRef> > dummy07;
+    edm::Wrapper<std::vector<std::pair<TrackingParticleRef, TrackPSimHitRef> > > dummy08;
+  };
+}

SimGeneral/TrackingAnalysis/src/classes_def.xml

@@ -0,0 +1,5 @@
+<lcgdict>
+  <class name="std::pair<TrackingParticleRef, TrackPSimHitRef>" persistent="true" /> 
+  <class name="std::vector<std::pair<TrackingParticleRef, TrackPSimHitRef> >" persistent="true" /> 
+  <class name="edm::Wrapper<std::vector<std::pair<TrackingParticleRef, TrackPSimHitRef> > >" />
+</lcgdict>

SimTracker/TrackAssociation/interface/CosmicParametersDefinerForTP.h

@@ -16,8 +16,8 @@ class CosmicParametersDefinerForTP : public ParametersDefinerForTP {
   CosmicParametersDefinerForTP(){};
   virtual ~CosmicParametersDefinerForTP() {};
 
-  virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticle& tp) const;
-  virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticle& tp) const;
+  virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticleRef tpr) const;
+  virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticleRef tpr) const;
 
   virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
 	const Charge ch, const Point & vertex, const LorentzVector& lv) const {
@@ -29,7 +29,12 @@ class CosmicParametersDefinerForTP : public ParametersDefinerForTP {
     return TrackingParticle::Point();
   }
 
+  void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) const {
+    simHitsTPAssoc = simHitsTPAssocToSet;
+  }
 
+ private:
+  mutable edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
 };
 
 

SimTracker/TrackAssociation/interface/ParametersDefinerForTP.h

@@ -8,6 +8,7 @@
  */
 
 #include <SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h>
+#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
 #include "DataFormats/Candidate/interface/Candidate.h"
 #include "FWCore/Framework/interface/Event.h"
 #include "FWCore/Framework/interface/ESHandle.h"      
@@ -23,12 +24,11 @@ class ParametersDefinerForTP {
     typedef math::XYZPointD Point; ///< point in the space
     typedef math::XYZTLorentzVectorD LorentzVector; ///< Lorentz vector
 
-
   virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
 	const Charge ch, const Point & vtx, const LorentzVector& lv) const;
 
-  virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticle& tp) const{
-    return momentum(iEvent, iSetup, tp.charge(),tp.vertex(),tp.p4());
+  virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticleRef tpr) const{
+    return momentum(iEvent, iSetup, tpr->charge(),tpr->vertex(),tpr->p4());
   }
 
   virtual TrackingParticle::Vector momentum(const edm::Event& iEvent, const edm::EventSetup& iSetup, const reco::Candidate& tp) const {
@@ -37,15 +37,17 @@ class ParametersDefinerForTP {
 
   virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup,
 	const Charge ch, const Point & vtx, const LorentzVector& lv) const;
-
-  virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticle& tp) const{
-    return vertex(iEvent, iSetup, tp.charge(),tp.vertex(),tp.p4());
+ 
+  virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup, const TrackingParticleRef& tpr) const{
+    return vertex(iEvent, iSetup, tpr->charge(),tpr->vertex(),tpr->p4());
   }
 
   virtual TrackingParticle::Point vertex(const edm::Event& iEvent, const edm::EventSetup& iSetup, const reco::Candidate& tp) const {
     return vertex(iEvent, iSetup, tp.charge(),tp.vertex(),tp.p4());
   }
 
+  virtual void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) const { }
+
 };
 
 

SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h

@@ -96,6 +96,8 @@ class TrackAssociatorByHits : public TrackAssociatorBase {
 
   const TrackingRecHit* getHitPtr(edm::OwnVector<TrackingRecHit>::const_iterator iter) const {return &*iter;}
   const TrackingRecHit* getHitPtr(trackingRecHit_iterator iter) const {return &**iter;}
+
+  edm::InputTag _simHitTpMapTag;
 };
 
 #endif

SimTracker/TrackAssociation/interface/TrackAssociatorByPosition.h

@@ -22,6 +22,8 @@
 
 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h>
 
+#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
+
 #include<map>
 
 //Note that the Association Map is filled with -ch2 and not chi2 because it is ordered using std::greater:
@@ -50,6 +52,7 @@ class TrackAssociatorByPosition : public TrackAssociatorBase {
        edm::LogError("TrackAssociatorByPosition")<<meth<<" mothed not recognized. Use dr or chi2.";     }
 
      theConsiderAllSimHits = iConfig.getParameter<bool>("ConsiderAllSimHits");
+     _simHitTpMapTag = iConfig.getParameter<edm::InputTag>("simHitTpMapTag");
    };
 
 
@@ -84,8 +87,9 @@ class TrackAssociatorByPosition : public TrackAssociatorBase {
   bool theConsiderAllSimHits;
 
   FreeTrajectoryState getState(const reco::Track &) const;
-  TrajectoryStateOnSurface getState(const TrackingParticle &)const;
-
+  TrajectoryStateOnSurface getState(const TrackingParticleRef)const;
+  mutable edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
+  edm::InputTag _simHitTpMapTag;
 };
 
 #endif