Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-28 01:32:08 

0001 /* \class IsolatedPixelTrackCandidateProducer
0002  *
0003  *  
0004  */
0005 
0006 #include <vector>
0007 #include <memory>
0008 #include <algorithm>
0009 #include <cmath>
0010 
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/Common/interface/Ref.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0017 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0018 #include "DataFormats/Common/interface/TriggerResults.h"
0019 // L1Extra
0020 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0021 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0022 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0023 
0024 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0025 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0026 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0027 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
0028 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0029 //vertices
0030 #include "DataFormats/VertexReco/interface/Vertex.h"
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 // Framework
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/stream/EDProducer.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "FWCore/Framework/interface/EventSetup.h"
0037 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0040 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0041 #include "FWCore/Utilities/interface/Exception.h"
0042 #include "FWCore/Utilities/interface/transform.h"
0043 
0044 // Math
0045 #include "Math/GenVector/VectorUtil.h"
0046 #include "Math/GenVector/PxPyPzE4D.h"
0047 #include "DataFormats/Math/interface/deltaR.h"
0048 
0049 //magF
0050 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0051 #include "MagneticField/Engine/interface/MagneticField.h"
0052 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0053 
0054 //for ECAL geometry
0055 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0056 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0057 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0058 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0059 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0060 
0061 //#define EDM_ML_DEBUG
0062 
0063 class IsolatedPixelTrackCandidateProducer : public edm::stream::EDProducer<> {
0064 public:
0065   IsolatedPixelTrackCandidateProducer(const edm::ParameterSet& ps);
0066   ~IsolatedPixelTrackCandidateProducer() override;
0067 
0068   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0069 
0070   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0071   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0072 
0073   double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0074   std::pair<double, double> GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ);
0075 
0076 private:
0077   struct seedAtEC {
0078     seedAtEC(unsigned int i, bool f, double et, double fi) : index(i), ok(f), eta(et), phi(fi) {}
0079     unsigned int index;
0080     bool ok;
0081     double eta, phi;
0082   };
0083 
0084   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0085   const edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_l1_;
0086   const edm::EDGetTokenT<reco::VertexCollection> tok_vert_;
0087   const std::vector<edm::EDGetTokenT<reco::TrackCollection> > toks_pix_;
0088   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0089   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0090 
0091   const std::string bfield_;
0092   const double prelimCone_;
0093   const double pixelIsolationConeSizeAtEC_;
0094   const double vtxCutSeed_;
0095   const double vtxCutIsol_;
0096   const double tauAssocCone_;
0097   const double tauUnbiasCone_;
0098   const double minPTrackValue_;
0099   const double maxPForIsolationValue_;
0100   const double ebEtaBoundary_;
0101 
0102   // these are read from the EventSetup, cannot be const
0103   double rEB_;
0104   double zEE_;
0105   double bfVal_;
0106 };
0107 
0108 IsolatedPixelTrackCandidateProducer::IsolatedPixelTrackCandidateProducer(const edm::ParameterSet& config)
0109     : tok_hlt_(consumes<trigger::TriggerFilterObjectWithRefs>(config.getParameter<edm::InputTag>("L1GTSeedLabel"))),
0110       tok_l1_(consumes<l1extra::L1JetParticleCollection>(config.getParameter<edm::InputTag>("L1eTauJetsSource"))),
0111       tok_vert_(consumes<reco::VertexCollection>(config.getParameter<edm::InputTag>("VertexLabel"))),
0112       toks_pix_(
0113           edm::vector_transform(config.getParameter<std::vector<edm::InputTag> >("PixelTracksSources"),
0114                                 [this](edm::InputTag const& tag) { return consumes<reco::TrackCollection>(tag); })),
0115       tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0116       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0117       bfield_(config.getParameter<std::string>("MagFieldRecordName")),
0118       prelimCone_(config.getParameter<double>("ExtrapolationConeSize")),
0119       pixelIsolationConeSizeAtEC_(config.getParameter<double>("PixelIsolationConeSizeAtEC")),
0120       vtxCutSeed_(config.getParameter<double>("MaxVtxDXYSeed")),
0121       vtxCutIsol_(config.getParameter<double>("MaxVtxDXYIsol")),
0122       tauAssocCone_(config.getParameter<double>("tauAssociationCone")),
0123       tauUnbiasCone_(config.getParameter<double>("tauUnbiasCone")),
0124       minPTrackValue_(config.getParameter<double>("minPTrack")),
0125       maxPForIsolationValue_(config.getParameter<double>("maxPTrackForIsolation")),
0126       ebEtaBoundary_(config.getParameter<double>("EBEtaBoundary")),
0127       rEB_(-1),
0128       zEE_(-1),
0129       bfVal_(0) {
0130   // Register the product
0131   produces<reco::IsolatedPixelTrackCandidateCollection>();
0132 }
0133 
0134 IsolatedPixelTrackCandidateProducer::~IsolatedPixelTrackCandidateProducer() {}
0135 
0136 void IsolatedPixelTrackCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0137   edm::ParameterSetDescription desc;
0138   std::vector<edm::InputTag> tracksrc = {edm::InputTag("hltPixelTracks")};
0139   desc.add<edm::InputTag>("L1eTauJetsSource", edm::InputTag("hltCaloStage2Digis", "Tau"));
0140   desc.add<double>("tauAssociationCone", 0.0);
0141   desc.add<double>("tauUnbiasCone", 1.2);
0142   desc.add<std::vector<edm::InputTag> >("PixelTracksSources", tracksrc);
0143   desc.add<double>("ExtrapolationConeSize", 1.0);
0144   desc.add<double>("PixelIsolationConeSizeAtEC", 40);
0145   desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("hltL1sIsoTrack"));
0146   desc.add<double>("MaxVtxDXYSeed", 101.0);
0147   desc.add<double>("MaxVtxDXYIsol", 101.0);
0148   desc.add<edm::InputTag>("VertexLabel", edm::InputTag("hltTrimmedPixelVertices"));
0149   desc.add<std::string>("MagFieldRecordName", "VolumeBasedMagneticField");
0150   desc.add<double>("minPTrack", 5.0);
0151   desc.add<double>("maxPTrackForIsolation", 3.0);
0152   desc.add<double>("EBEtaBoundary", 1.479);
0153   descriptions.add("isolPixelTrackProd", desc);
0154 }
0155 
0156 void IsolatedPixelTrackCandidateProducer::beginRun(const edm::Run& run, const edm::EventSetup& theEventSetup) {
0157   const CaloGeometry* pG = &theEventSetup.getData(tok_geom_);
0158 
0159   const double rad(dynamic_cast<const EcalBarrelGeometry*>(pG->getSubdetectorGeometry(DetId::Ecal, EcalBarrel))
0160                        ->avgRadiusXYFrontFaceCenter());
0161   const double zz(dynamic_cast<const EcalEndcapGeometry*>(pG->getSubdetectorGeometry(DetId::Ecal, EcalEndcap))
0162                       ->avgAbsZFrontFaceCenter());
0163 
0164   rEB_ = rad;
0165   zEE_ = zz;
0166 
0167   const MagneticField* vbfField = &theEventSetup.getData(tok_bFieldH_);
0168   const VolumeBasedMagneticField* vbfCPtr = dynamic_cast<const VolumeBasedMagneticField*>(&(*vbfField));
0169   GlobalVector BField = vbfCPtr->inTesla(GlobalPoint(0, 0, 0));
0170   bfVal_ = BField.mag();
0171 }
0172 
0173 void IsolatedPixelTrackCandidateProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0174   auto trackCollection = std::make_unique<reco::IsolatedPixelTrackCandidateCollection>();
0175 
0176   //create vector of refs from input collections
0177   std::vector<reco::TrackRef> pixelTrackRefs;
0178 #ifdef EDM_ML_DEBUG
0179   edm::LogVerbatim("HcalIsoTrack") << "IsolatedPixelTrakCandidate: with" << toks_pix_.size()
0180                                    << " candidates to start with\n";
0181 #endif
0182   for (unsigned int iPix = 0; iPix < toks_pix_.size(); iPix++) {
0183     const edm::Handle<reco::TrackCollection>& iPixCol = theEvent.getHandle(toks_pix_[iPix]);
0184     for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
0185       pixelTrackRefs.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0186     }
0187   }
0188 
0189   const edm::Handle<l1extra::L1JetParticleCollection>& l1eTauJets = theEvent.getHandle(tok_l1_);
0190 
0191   const edm::Handle<reco::VertexCollection>& pVert = theEvent.getHandle(tok_vert_);
0192 
0193   double drMaxL1Track_ = tauAssocCone_;
0194 
0195   int ntr = 0;
0196   std::vector<seedAtEC> VecSeedsatEC;
0197   //loop to select isolated tracks
0198   for (unsigned iS = 0; iS < pixelTrackRefs.size(); iS++) {
0199     bool vtxMatch = false;
0200     //associate to vertex (in Z)
0201     reco::VertexCollection::const_iterator vitSel;
0202     double minDZ = 1000;
0203     bool found(false);
0204     for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0205       if (std::abs(pixelTrackRefs[iS]->dz(vit->position())) < minDZ) {
0206         minDZ = std::abs(pixelTrackRefs[iS]->dz(vit->position()));
0207         vitSel = vit;
0208         found = true;
0209       }
0210     }
0211     //cut on dYX:
0212     if (found) {
0213       if (std::abs(pixelTrackRefs[iS]->dxy(vitSel->position())) < vtxCutSeed_)
0214         vtxMatch = true;
0215     } else {
0216       vtxMatch = true;
0217     }
0218 
0219     //check taujet matching
0220     bool tmatch = false;
0221     l1extra::L1JetParticleCollection::const_iterator selj;
0222     for (l1extra::L1JetParticleCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0223       if (reco::deltaR(pixelTrackRefs[iS]->momentum().eta(),
0224                        pixelTrackRefs[iS]->momentum().phi(),
0225                        tj->momentum().eta(),
0226                        tj->momentum().phi()) > drMaxL1Track_)
0227         continue;
0228       selj = tj;
0229       tmatch = true;
0230     }  //loop over L1 tau
0231 
0232     //propagate seed track to ECAL surface:
0233     std::pair<double, double> seedCooAtEC;
0234     // in case vertex is found:
0235     if (found)
0236       seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0237                                     pixelTrackRefs[iS]->phi(),
0238                                     pixelTrackRefs[iS]->pt(),
0239                                     pixelTrackRefs[iS]->charge(),
0240                                     vitSel->z());
0241     //in case vertex is not found:
0242     else
0243       seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0244                                     pixelTrackRefs[iS]->phi(),
0245                                     pixelTrackRefs[iS]->pt(),
0246                                     pixelTrackRefs[iS]->charge(),
0247                                     0);
0248     seedAtEC seed(iS, (tmatch || vtxMatch), seedCooAtEC.first, seedCooAtEC.second);
0249     VecSeedsatEC.push_back(seed);
0250   }
0251 #ifdef EDM_ML_DEBUG
0252   edm::LogVerbatim("HcalIsoTrack") << "IsolatedPixelTrakCandidate: " << VecSeedsatEC.size()
0253                                    << " seeds after propagation\n";
0254 #endif
0255 
0256   for (unsigned int i = 0; i < VecSeedsatEC.size(); i++) {
0257     unsigned int iSeed = VecSeedsatEC[i].index;
0258     if (!VecSeedsatEC[i].ok)
0259       continue;
0260     if (pixelTrackRefs[iSeed]->p() < minPTrackValue_)
0261       continue;
0262     l1extra::L1JetParticleCollection::const_iterator selj;
0263     for (l1extra::L1JetParticleCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0264       if (reco::deltaR(pixelTrackRefs[iSeed]->momentum().eta(),
0265                        pixelTrackRefs[iSeed]->momentum().phi(),
0266                        tj->momentum().eta(),
0267                        tj->momentum().phi()) > drMaxL1Track_)
0268         continue;
0269       selj = tj;
0270     }  //loop over L1 tau
0271     double maxP = 0;
0272     double sumP = 0;
0273     for (unsigned int j = 0; j < VecSeedsatEC.size(); j++) {
0274       if (i == j)
0275         continue;
0276       unsigned int iSurr = VecSeedsatEC[j].index;
0277       //define preliminary cone around seed track impact point from which tracks will be extrapolated:
0278       if (reco::deltaR(pixelTrackRefs[iSeed]->eta(),
0279                        pixelTrackRefs[iSeed]->phi(),
0280                        pixelTrackRefs[iSurr]->eta(),
0281                        pixelTrackRefs[iSurr]->phi()) > prelimCone_)
0282         continue;
0283       double minDZ2(1000);
0284       bool found(false);
0285       reco::VertexCollection::const_iterator vitSel2;
0286       for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0287         if (std::abs(pixelTrackRefs[iSurr]->dz(vit->position())) < minDZ2) {
0288           minDZ2 = std::abs(pixelTrackRefs[iSurr]->dz(vit->position()));
0289           vitSel2 = vit;
0290           found = true;
0291         }
0292       }
0293       //cut ot dXY:
0294       if (found && std::abs(pixelTrackRefs[iSurr]->dxy(vitSel2->position())) > vtxCutIsol_)
0295         continue;
0296       //calculate distance at ECAL surface and update isolation:
0297       if (getDistInCM(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi, VecSeedsatEC[j].eta, VecSeedsatEC[j].phi) <
0298           pixelIsolationConeSizeAtEC_) {
0299         sumP += pixelTrackRefs[iSurr]->p();
0300         if (pixelTrackRefs[iSurr]->p() > maxP)
0301           maxP = pixelTrackRefs[iSurr]->p();
0302       }
0303     }
0304     if (maxP < maxPForIsolationValue_) {
0305       reco::IsolatedPixelTrackCandidate newCandidate(
0306           pixelTrackRefs[iSeed], l1extra::L1JetParticleRef(l1eTauJets, selj - l1eTauJets->begin()), maxP, sumP);
0307       newCandidate.setEtaPhiEcal(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi);
0308       trackCollection->push_back(newCandidate);
0309       ntr++;
0310     }
0311   }
0312   // put the product in the event
0313   theEvent.put(std::move(trackCollection));
0314 #ifdef EDM_ML_DEBUG
0315   edm::LogVerbatim("HcalIsoTrack") << "IsolatedPixelTrackCandidate: Final # of candiates " << ntr << "\n";
0316 #endif
0317 }
0318 
0319 double IsolatedPixelTrackCandidateProducer::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0320   double Rec;
0321   double theta1 = 2 * atan(exp(-eta1));
0322   double theta2 = 2 * atan(exp(-eta2));
0323   if (std::abs(eta1) < 1.479)
0324     Rec = rEB_;  //radius of ECAL barrel
0325   else if (std::abs(eta1) > 1.479 && std::abs(eta1) < 7.0)
0326     Rec = tan(theta1) * zEE_;  //distance from IP to ECAL endcap
0327   else
0328     return 1000;
0329 
0330   //|vect| times tg of acos(scalar product)
0331   double angle =
0332       acos((sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) + cos(theta1) * cos(theta2)));
0333   if (angle < M_PI_2)
0334     return std::abs((Rec / sin(theta1)) * tan(angle));
0335   else
0336     return 1000;
0337 }
0338 
0339 std::pair<double, double> IsolatedPixelTrackCandidateProducer::GetEtaPhiAtEcal(
0340     double etaIP, double phiIP, double pT, int charge, double vtxZ) {
0341   double deltaPhi = 0;
0342   double etaEC = 100;
0343   double phiEC = 100;
0344 
0345   double Rcurv = 9999999;
0346   if (bfVal_ != 0)
0347     Rcurv = pT * 33.3 * 100 / (bfVal_ * 10);  //r(m)=pT(GeV)*33.3/B(kG)
0348 
0349   double ecDist = zEE_;  //distance to ECAL andcap from IP (cm), 317 - ecal (not preshower), preshower -300
0350   double ecRad = rEB_;   //radius of ECAL barrel (cm)
0351   double theta = 2 * atan(exp(-etaIP));
0352   double zNew = 0;
0353   if (theta > M_PI_2)
0354     theta = M_PI - theta;
0355   if (std::abs(etaIP) < ebEtaBoundary_) {
0356     if ((0.5 * ecRad / Rcurv) > 1) {
0357       etaEC = 10000;
0358       deltaPhi = 0;
0359     } else {
0360       deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0361       double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
0362       double z = ecRad / tan(theta);
0363       if (etaIP > 0)
0364         zNew = z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
0365       else
0366         zNew = -z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
0367       double zAbs = std::abs(zNew);
0368       if (zAbs < ecDist) {
0369         etaEC = -log(tan(0.5 * atan(ecRad / zAbs)));
0370         deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0371       }
0372       if (zAbs > ecDist) {
0373         zAbs = (std::abs(etaIP) / etaIP) * ecDist;
0374         double Zflight = std::abs(zAbs - vtxZ);
0375         double alpha = (Zflight * ecRad) / (z * Rcurv);
0376         double Rec = 2 * Rcurv * sin(alpha / 2);
0377         deltaPhi = -charge * alpha / 2;
0378         etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0379       }
0380     }
0381   } else {
0382     zNew = (std::abs(etaIP) / etaIP) * ecDist;
0383     double Zflight = std::abs(zNew - vtxZ);
0384     double Rvirt = std::abs(Zflight * tan(theta));
0385     double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
0386     deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
0387     etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0388   }
0389 
0390   if (zNew < 0)
0391     etaEC = -etaEC;
0392   phiEC = phiIP + deltaPhi;
0393 
0394   if (phiEC < -M_PI)
0395     phiEC += M_2_PI;
0396   if (phiEC > M_PI)
0397     phiEC -= M_2_PI;
0398 
0399   std::pair<double, double> retVal(etaEC, phiEC);
0400   return retVal;
0401 }
0402 
0403 #include "FWCore/PluginManager/interface/ModuleDef.h"
0404 #include "FWCore/Framework/interface/MakerMacros.h"
0405 
0406 DEFINE_FWK_MODULE(IsolatedPixelTrackCandidateProducer);
[ 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!