Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59: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 #ifdef EDM_ML_DEBUG
0195   int ntr = 0;
0196 #endif
0197   std::vector<seedAtEC> VecSeedsatEC;
0198   //loop to select isolated tracks
0199   for (unsigned iS = 0; iS < pixelTrackRefs.size(); iS++) {
0200     bool vtxMatch = false;
0201     //associate to vertex (in Z)
0202     reco::VertexCollection::const_iterator vitSel;
0203     double minDZ = 1000;
0204     bool found(false);
0205     for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0206       if (std::abs(pixelTrackRefs[iS]->dz(vit->position())) < minDZ) {
0207         minDZ = std::abs(pixelTrackRefs[iS]->dz(vit->position()));
0208         vitSel = vit;
0209         found = true;
0210       }
0211     }
0212     //cut on dYX:
0213     if (found) {
0214       if (std::abs(pixelTrackRefs[iS]->dxy(vitSel->position())) < vtxCutSeed_)
0215         vtxMatch = true;
0216     } else {
0217       vtxMatch = true;
0218     }
0219 
0220     //check taujet matching
0221     bool tmatch = false;
0222     l1extra::L1JetParticleCollection::const_iterator selj;
0223     for (l1extra::L1JetParticleCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0224       if (reco::deltaR(pixelTrackRefs[iS]->momentum().eta(),
0225                        pixelTrackRefs[iS]->momentum().phi(),
0226                        tj->momentum().eta(),
0227                        tj->momentum().phi()) > drMaxL1Track_)
0228         continue;
0229       selj = tj;
0230       tmatch = true;
0231     }  //loop over L1 tau
0232 
0233     //propagate seed track to ECAL surface:
0234     std::pair<double, double> seedCooAtEC;
0235     // in case vertex is found:
0236     if (found)
0237       seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0238                                     pixelTrackRefs[iS]->phi(),
0239                                     pixelTrackRefs[iS]->pt(),
0240                                     pixelTrackRefs[iS]->charge(),
0241                                     vitSel->z());
0242     //in case vertex is not found:
0243     else
0244       seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0245                                     pixelTrackRefs[iS]->phi(),
0246                                     pixelTrackRefs[iS]->pt(),
0247                                     pixelTrackRefs[iS]->charge(),
0248                                     0);
0249     seedAtEC seed(iS, (tmatch || vtxMatch), seedCooAtEC.first, seedCooAtEC.second);
0250     VecSeedsatEC.push_back(seed);
0251   }
0252 #ifdef EDM_ML_DEBUG
0253   edm::LogVerbatim("HcalIsoTrack") << "IsolatedPixelTrakCandidate: " << VecSeedsatEC.size()
0254                                    << " seeds after propagation\n";
0255 #endif
0256 
0257   for (unsigned int i = 0; i < VecSeedsatEC.size(); i++) {
0258     unsigned int iSeed = VecSeedsatEC[i].index;
0259     if (!VecSeedsatEC[i].ok)
0260       continue;
0261     if (pixelTrackRefs[iSeed]->p() < minPTrackValue_)
0262       continue;
0263     l1extra::L1JetParticleCollection::const_iterator selj;
0264     for (l1extra::L1JetParticleCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0265       if (reco::deltaR(pixelTrackRefs[iSeed]->momentum().eta(),
0266                        pixelTrackRefs[iSeed]->momentum().phi(),
0267                        tj->momentum().eta(),
0268                        tj->momentum().phi()) > drMaxL1Track_)
0269         continue;
0270       selj = tj;
0271     }  //loop over L1 tau
0272     double maxP = 0;
0273     double sumP = 0;
0274     for (unsigned int j = 0; j < VecSeedsatEC.size(); j++) {
0275       if (i == j)
0276         continue;
0277       unsigned int iSurr = VecSeedsatEC[j].index;
0278       //define preliminary cone around seed track impact point from which tracks will be extrapolated:
0279       if (reco::deltaR(pixelTrackRefs[iSeed]->eta(),
0280                        pixelTrackRefs[iSeed]->phi(),
0281                        pixelTrackRefs[iSurr]->eta(),
0282                        pixelTrackRefs[iSurr]->phi()) > prelimCone_)
0283         continue;
0284       double minDZ2(1000);
0285       bool found(false);
0286       reco::VertexCollection::const_iterator vitSel2;
0287       for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0288         if (std::abs(pixelTrackRefs[iSurr]->dz(vit->position())) < minDZ2) {
0289           minDZ2 = std::abs(pixelTrackRefs[iSurr]->dz(vit->position()));
0290           vitSel2 = vit;
0291           found = true;
0292         }
0293       }
0294       //cut ot dXY:
0295       if (found && std::abs(pixelTrackRefs[iSurr]->dxy(vitSel2->position())) > vtxCutIsol_)
0296         continue;
0297       //calculate distance at ECAL surface and update isolation:
0298       if (getDistInCM(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi, VecSeedsatEC[j].eta, VecSeedsatEC[j].phi) <
0299           pixelIsolationConeSizeAtEC_) {
0300         sumP += pixelTrackRefs[iSurr]->p();
0301         if (pixelTrackRefs[iSurr]->p() > maxP)
0302           maxP = pixelTrackRefs[iSurr]->p();
0303       }
0304     }
0305     if (maxP < maxPForIsolationValue_) {
0306       reco::IsolatedPixelTrackCandidate newCandidate(
0307           pixelTrackRefs[iSeed], l1extra::L1JetParticleRef(l1eTauJets, selj - l1eTauJets->begin()), maxP, sumP);
0308       newCandidate.setEtaPhiEcal(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi);
0309       trackCollection->push_back(newCandidate);
0310 #ifdef EDM_ML_DEBUG
0311       ntr++;
0312 #endif
0313     }
0314   }
0315   // put the product in the event
0316   theEvent.put(std::move(trackCollection));
0317 #ifdef EDM_ML_DEBUG
0318   edm::LogVerbatim("HcalIsoTrack") << "IsolatedPixelTrackCandidate: Final # of candiates " << ntr << "\n";
0319 #endif
0320 }
0321 
0322 double IsolatedPixelTrackCandidateProducer::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0323   double Rec;
0324   double theta1 = 2 * atan(exp(-eta1));
0325   double theta2 = 2 * atan(exp(-eta2));
0326   if (std::abs(eta1) < 1.479)
0327     Rec = rEB_;  //radius of ECAL barrel
0328   else if (std::abs(eta1) > 1.479 && std::abs(eta1) < 7.0)
0329     Rec = tan(theta1) * zEE_;  //distance from IP to ECAL endcap
0330   else
0331     return 1000;
0332 
0333   //|vect| times tg of acos(scalar product)
0334   double angle =
0335       acos((sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) + cos(theta1) * cos(theta2)));
0336   if (angle < M_PI_2)
0337     return std::abs((Rec / sin(theta1)) * tan(angle));
0338   else
0339     return 1000;
0340 }
0341 
0342 std::pair<double, double> IsolatedPixelTrackCandidateProducer::GetEtaPhiAtEcal(
0343     double etaIP, double phiIP, double pT, int charge, double vtxZ) {
0344   double deltaPhi = 0;
0345   double etaEC = 100;
0346   double phiEC = 100;
0347 
0348   double Rcurv = 9999999;
0349   if (bfVal_ != 0)
0350     Rcurv = pT * 33.3 * 100 / (bfVal_ * 10);  //r(m)=pT(GeV)*33.3/B(kG)
0351 
0352   double ecDist = zEE_;  //distance to ECAL andcap from IP (cm), 317 - ecal (not preshower), preshower -300
0353   double ecRad = rEB_;   //radius of ECAL barrel (cm)
0354   double theta = 2 * atan(exp(-etaIP));
0355   double zNew = 0;
0356   if (theta > M_PI_2)
0357     theta = M_PI - theta;
0358   if (std::abs(etaIP) < ebEtaBoundary_) {
0359     if ((0.5 * ecRad / Rcurv) > 1) {
0360       etaEC = 10000;
0361       deltaPhi = 0;
0362     } else {
0363       deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0364       double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
0365       double z = ecRad / tan(theta);
0366       if (etaIP > 0)
0367         zNew = z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
0368       else
0369         zNew = -z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
0370       double zAbs = std::abs(zNew);
0371       if (zAbs < ecDist) {
0372         etaEC = -log(tan(0.5 * atan(ecRad / zAbs)));
0373         deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0374       }
0375       if (zAbs > ecDist) {
0376         zAbs = (std::abs(etaIP) / etaIP) * ecDist;
0377         double Zflight = std::abs(zAbs - vtxZ);
0378         double alpha = (Zflight * ecRad) / (z * Rcurv);
0379         double Rec = 2 * Rcurv * sin(alpha / 2);
0380         deltaPhi = -charge * alpha / 2;
0381         etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0382       }
0383     }
0384   } else {
0385     zNew = (std::abs(etaIP) / etaIP) * ecDist;
0386     double Zflight = std::abs(zNew - vtxZ);
0387     double Rvirt = std::abs(Zflight * tan(theta));
0388     double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
0389     deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
0390     etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0391   }
0392 
0393   if (zNew < 0)
0394     etaEC = -etaEC;
0395   phiEC = phiIP + deltaPhi;
0396 
0397   if (phiEC < -M_PI)
0398     phiEC += M_2_PI;
0399   if (phiEC > M_PI)
0400     phiEC -= M_2_PI;
0401 
0402   std::pair<double, double> retVal(etaEC, phiEC);
0403   return retVal;
0404 }
0405 
0406 #include "FWCore/PluginManager/interface/ModuleDef.h"
0407 #include "FWCore/Framework/interface/MakerMacros.h"
0408 
0409 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!