Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-02-05 23:51:09 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!