RecoAlgos/interface/CosmicTrackingParticleSelector.h

0001 #ifndef RecoSelectors_CosmicTrackingParticleSelector_h
0002 #define RecoSelectors_CosmicTrackingParticleSelector_h
0003 /* \class CosmicTrackingParticleSelector
0004  *
0005  * \author Yanyan Gao, FNAL
0006  *
0007  *  $Date: 2013/06/24 12:25:14 $
0008  *  $Revision: 1.5 $
0009  *
0010  */
0011 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0018 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0019 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0020 #include "MagneticField/Engine/interface/MagneticField.h"
0021 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0022
0023 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0025 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0026 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0027
0028 #include "DataFormats/GeometrySurface/interface/Surface.h"
0029 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
0030 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0031 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0032 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0033
0034 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036
0037 class TrajectoryStateClosestToBeamLineBuilder;
0038
0039 class CosmicTrackingParticleSelector {
0040 public:
0041   typedef TrackingParticleCollection collection;
0042   typedef std::vector<const TrackingParticle*> container;
0043   typedef container::const_iterator const_iterator;
0044
0045   CosmicTrackingParticleSelector() {}
0046
0047   CosmicTrackingParticleSelector(double ptMin,
0048                                  double minRapidity,
0049                                  double maxRapidity,
0050                                  double tip,
0051                                  double lip,
0052                                  int minHit,
0053                                  bool chargedOnly,
0054                                  const std::vector<int>& pdgId = std::vector<int>())
0055       : ptMin_(ptMin),
0056         minRapidity_(minRapidity),
0057         maxRapidity_(maxRapidity),
0058         tip_(tip),
0059         lip_(lip),
0060         minHit_(minHit),
0061         chargedOnly_(chargedOnly),
0062         pdgId_(pdgId) {}
0063
0064   CosmicTrackingParticleSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0065       : ptMin_(cfg.getParameter<double>("ptMin")),
0066         minRapidity_(cfg.getParameter<double>("minRapidity")),
0067         maxRapidity_(cfg.getParameter<double>("maxRapidity")),
0068         tip_(cfg.getParameter<double>("tip")),
0069         lip_(cfg.getParameter<double>("lip")),
0070         minHit_(cfg.getParameter<int>("minHit")),
0071         chargedOnly_(cfg.getParameter<bool>("chargedOnly")),
0072         pdgId_(cfg.getParameter<std::vector<int> >("pdgId")),
0073         beamSpotToken_(iC.consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
0074         globalTrackingGeomToken_(iC.esConsumes()),
0075         theMFToken_(iC.esConsumes()) {}
0076
0077   void select(const edm::Handle<collection>& c, const edm::Event& event, const edm::EventSetup& setup) {
0078     selected_.clear();
0079     edm::Handle<reco::BeamSpot> beamSpot;
0080     event.getByToken(beamSpotToken_, beamSpot);
0081     for (TrackingParticleCollection::const_iterator itp = c->begin(); itp != c->end(); ++itp)
0082       if (operator()(TrackingParticleRef(c, itp - c->begin()), beamSpot.product(), event, setup)) {
0083         selected_.push_back(&*itp);
0084       }
0085   }
0086
0087   const_iterator begin() const { return selected_.begin(); }
0088   const_iterator end() const { return selected_.end(); }
0089
0090   void initEvent(edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssocToSet) const {
0091     simHitsTPAssoc = simHitsTPAssocToSet;
0092   }
0093
0094   // Operator() performs the selection: e.g. if (tPSelector(tp, bs, event, evtsetup)) {...
0095   bool operator()(const TrackingParticleRef tpr,
0096                   const reco::BeamSpot* bs,
0097                   const edm::Event& iEvent,
0098                   const edm::EventSetup& iSetup) const {
0099     if (chargedOnly_ && tpr->charge() == 0)
0100       return false;  //select only if charge!=0
0101     //bool testId = false;
0102     //unsigned int idSize = pdgId_.size();
0103     //if (idSize==0) testId = true;
0104     //else for (unsigned int it=0;it!=idSize;++it){
0105     //if (tpr->pdgId()==pdgId_[it]) testId = true;
0106     //}
0107
0108     edm::ESHandle<GlobalTrackingGeometry> theGeometry = iSetup.getHandle(globalTrackingGeomToken_);
0109
0110     GlobalVector finalGV(0, 0, 0);
0111     GlobalPoint finalGP(0, 0, 0);
0112     GlobalVector momentum(0, 0, 0);  //At the PCA
0113     GlobalPoint vertex(0, 0, 0);     //At the PCA
0114     double radius(9999);
0115     bool found(false);
0116
0117     int ii = 0;
0118     DetId::Detector det;
0119     int subdet;
0120
0121     edm::LogVerbatim("CosmicTrackingParticleSelector")
0122         << "TOT Number of PSimHits = " << tpr->numberOfHits()
0123         << ", Number of Tracker PSimHits = " << tpr->numberOfTrackerHits() << "\n";
0124
0125     if (simHitsTPAssoc.isValid() == 0) {
0126       edm::LogError("CosmicTrackingParticleSelector") << "Invalid handle!";
0127       return false;
0128     }
0129     std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0130         tpr, TrackPSimHitRef());  //SimHit is dummy: for simHitTPAssociationListGreater
0131                                   // sorting only the cluster is needed
0132     auto range = std::equal_range(simHitsTPAssoc->begin(),
0133                                   simHitsTPAssoc->end(),
0134                                   clusterTPpairWithDummyTP,
0135                                   SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0136     for (auto ip = range.first; ip != range.second; ++ip) {
0137       TrackPSimHitRef it = ip->second;
0138       ++ii;
0139       const GeomDet* tmpDet = theGeometry->idToDet(DetId(it->detUnitId()));
0140       if (!tmpDet) {
0141         edm::LogVerbatim("CosmicTrackingParticleSelector")
0142             << "***WARNING:  PSimHit " << ii << ", no GeomDet for: " << it->detUnitId() << ". Skipping it.";
0143         continue;
0144       } else {
0145         det = DetId(it->detUnitId()).det();
0146         subdet = DetId(it->detUnitId()).subdetId();
0147       }
0148
0149       LocalVector lv = it->momentumAtEntry();
0150       Local3DPoint lp = it->localPosition();
0151       GlobalVector gv = tmpDet->surface().toGlobal(lv);
0152       GlobalPoint gp = tmpDet->surface().toGlobal(lp);
0153       edm::LogVerbatim("CosmicTrackingParticleSelector")
0154           << "PSimHit " << ii << ", Detector = " << det << ", subdet = " << subdet << "\t Radius = " << gp.perp()
0155           << ", z = " << gp.z() << "\t     pt = " << gv.perp() << ", pz = " << gv.z();
0156       edm::LogVerbatim("CosmicTrackingParticleSelector")
0157           << "\t trackId = " << it->trackId() << ", particleType = " << it->particleType()
0158           << ", processType = " << it->processType();
0159
0160       // discard hits related to low energy debris from the primary particle
0161       if (it->processType() != 0)
0162         continue;
0163
0164       if (gp.perp() < radius) {
0165         found = true;
0166         radius = gp.perp();
0167         finalGV = gv;
0168         finalGP = gp;
0169       }
0170     }
0171     edm::LogVerbatim("CosmicTrackingParticleSelector")
0172         << "\n"
0173         << "FINAL State at InnerMost Hit:        Radius = " << finalGP.perp() << ", z = " << finalGP.z()
0174         << ", pt = " << finalGV.perp() << ", pz = " << finalGV.z();
0175
0176     if (!found)
0177       return false;
0178     else {
0179       FreeTrajectoryState ftsAtProduction(finalGP, finalGV, TrackCharge(tpr->charge()), &iSetup.getData(theMFToken_));
0180       TSCBLBuilderNoMaterial tscblBuilder;
0181       //as in TrackProducerAlgorithm
0182       TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction, *bs);
0183       if (!tsAtClosestApproach.isValid()) {
0184         edm::LogVerbatim("CosmicTrackingParticleSelector")
0185             << "*** WARNING in CosmicTrackingParticleSelector: tsAtClosestApproach is not valid."
0186             << "\n";
0187         return false;
0188       } else {
0189         momentum = tsAtClosestApproach.trackStateAtPCA().momentum();
0190         vertex = tsAtClosestApproach.trackStateAtPCA().position();
0191
0192         edm::LogVerbatim("CosmicTrackingParticleSelector")
0193             << "FINAL State extrapolated at PCA: Radius = " << vertex.perp() << ", z = " << vertex.z()
0194             << ", pt = " << momentum.perp() << ", pz = " << momentum.z() << "\n";
0195
0196         return (tpr->numberOfTrackerLayers() >= minHit_ && sqrt(momentum.perp2()) >= ptMin_ &&
0197                 momentum.eta() >= minRapidity_ && momentum.eta() <= maxRapidity_ && sqrt(vertex.perp2()) <= tip_ &&
0198                 fabs(vertex.z()) <= lip_);
0199       }
0200     }
0201   }
0202
0203   size_t size() const { return selected_.size(); }
0204
0205 private:
0206   double ptMin_;
0207   double minRapidity_;
0208   double maxRapidity_;
0209   double tip_;
0210   double lip_;
0211   int minHit_;
0212   bool chargedOnly_;
0213   std::vector<int> pdgId_;
0214   container selected_;
0215   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0216   edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> globalTrackingGeomToken_;
0217   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken_;
0218
0219   mutable edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0220 };
0221
0222 #endif