Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-11-05 00:08:12 

0001 // -*- C++ -*-
0002 //
0003 // Package:    PFTracking
0004 // Class:      GoodSeedProducer
0005 //
0006 // Original Author:  Michele Pioppi
0007 // March 2010. F. Beaudette. Produce PreId information
0008 
0009 /// \brief Abstract
0010 /*!
0011 \author Michele Pioppi
0012 \date January 2007
0013 
0014  GoodSeedProducer is the class
0015  for electron preidentification in PFLow FW.
0016  It reads refitted tracks and PFCluster collection, 
0017  and following some criteria divides electrons from hadrons.
0018  Then it saves the seed of the tracks preidentified as electrons.
0019  It also transform  all the tracks in the first PFRecTrack collection.
0020 */
0021 
0022 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0023 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0024 #include "CondFormats/GBRForest/interface/GBRForest.h"
0025 #include "DataFormats/Common/interface/ValueMap.h"
0026 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0027 #include "DataFormats/Math/interface/deltaR.h"
0028 #include "DataFormats/Math/interface/LorentzVector.h"
0029 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PreId.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PreIdFwd.h"
0032 #include "DataFormats/TrackReco/interface/Track.h"
0033 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0034 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0035 #include "FWCore/Framework/interface/ESHandle.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/stream/EDProducer.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040 #include "FWCore/ParameterSet/interface/FileInPath.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"
0045 #include "RecoParticleFlow/PFTracking/interface/PFGeometry.h"
0046 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
0047 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0048 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0049 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0050 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0051 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0052 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0053 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0054 
0055 #include "TMath.h"
0056 #include "Math/VectorUtil.h"
0057 
0058 #include <fstream>
0059 
0060 namespace goodseedhelpers {
0061   struct HeavyObjectCache {
0062     constexpr static unsigned int kMaxWeights = 9;
0063     HeavyObjectCache(const edm::ParameterSet& conf);
0064     std::array<std::unique_ptr<const GBRForest>, kMaxWeights> gbr;
0065   };
0066 }  // namespace goodseedhelpers
0067 
0068 class GoodSeedProducer final : public edm::stream::EDProducer<edm::GlobalCache<goodseedhelpers::HeavyObjectCache>> {
0069   typedef TrajectoryStateOnSurface TSOS;
0070 
0071 public:
0072   explicit GoodSeedProducer(const edm::ParameterSet&, const goodseedhelpers::HeavyObjectCache*);
0073 
0074   static std::unique_ptr<goodseedhelpers::HeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& conf) {
0075     return std::make_unique<goodseedhelpers::HeavyObjectCache>(conf);
0076   }
0077 
0078   static void globalEndJob(goodseedhelpers::HeavyObjectCache const*) {}
0079 
0080 private:
0081   void beginRun(const edm::Run& run, const edm::EventSetup&) override;
0082   void produce(edm::Event&, const edm::EventSetup&) override;
0083 
0084   ///Find the bin in pt and eta
0085   int getBin(float, float);
0086 
0087   void fillPreIdRefValueMap(edm::Handle<reco::TrackCollection> tkhandle,
0088                             const edm::OrphanHandle<reco::PreIdCollection>&,
0089                             edm::ValueMap<reco::PreIdRef>::Filler& filler);
0090   // ----------member data ---------------------------
0091 
0092   ///Name of the Seed(Ckf) Collection
0093   std::string preidckf_;
0094 
0095   ///Name of the Seed(Gsf) Collection
0096   std::string preidgsf_;
0097 
0098   ///Name of the preid Collection (FB)
0099   std::string preidname_;
0100 
0101   // needed by the above
0102   TkClonerImpl hitCloner;
0103 
0104   ///PFTrackTransformer
0105   std::unique_ptr<PFTrackTransformer> pfTransformer_;
0106 
0107   ///Number of hits in the seed;
0108   int nHitsInSeed_;
0109 
0110   ///Minimum transverse momentum and maximum pseudorapidity
0111   double minPt_;
0112   double maxPt_;
0113   double maxEta_;
0114 
0115   double HcalIsolWindow_;
0116   double EcalStripSumE_minClusEnergy_;
0117   double EcalStripSumE_deltaEta_;
0118   double EcalStripSumE_deltaPhiOverQ_minValue_;
0119   double EcalStripSumE_deltaPhiOverQ_maxValue_;
0120   double minEoverP_;
0121   double maxHoverP_;
0122 
0123   ///Cut on the energy of the clusters
0124   double clusThreshold_;
0125 
0126   ///Min and MAx allowed values forEoverP
0127   double minEp_;
0128   double maxEp_;
0129 
0130   ///Produce the Seed for Ckf tracks?
0131   bool produceCkfseed_;
0132 
0133   ///  switch to disable the pre-id
0134   bool disablePreId_;
0135 
0136   ///Produce the pre-id debugging collection
0137   bool producePreId_;
0138 
0139   /// Threshold to save Pre Idinfo
0140   double PtThresholdSavePredId_;
0141 
0142   ///vector of thresholds for different bins of eta and pt
0143   float thr[150];
0144 
0145   // ----------access to event data
0146   edm::ParameterSet conf_;
0147   edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagPSLabel_;
0148   edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagECLabel_;
0149   edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagHCLabel_;
0150   std::vector<edm::EDGetTokenT<std::vector<Trajectory>>> trajContainers_;
0151   std::vector<edm::EDGetTokenT<reco::TrackCollection>> tracksContainers_;
0152 
0153   std::string fitterName_;
0154   std::string smootherName_;
0155   std::string propagatorName_;
0156   std::string trackerRecHitBuilderName_;
0157 
0158   std::unique_ptr<PFResolutionMap> resMapEtaECAL_;
0159   std::unique_ptr<PFResolutionMap> resMapPhiECAL_;
0160 
0161   const edm::ESGetToken<TrajectoryFitter, TrajectoryFitter::Record> fitterToken_;
0162   const edm::ESGetToken<TrajectorySmoother, TrajectoryFitter::Record> smootherToken_;
0163   const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> trackerRecHitBuilderToken_;
0164   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0165   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldTokenBeginRun_;
0166 
0167   ///TRACK QUALITY
0168   bool useQuality_;
0169   reco::TrackBase::TrackQuality trackQuality_;
0170 
0171   ///VARIABLES NEEDED FOR TMVA
0172   float nhit, dpt, chired, chiRatio;
0173   float chikfred, trk_ecalDeta, trk_ecalDphi;
0174   double Min_dr_;
0175 
0176   ///USE OF TMVA
0177   bool useTmva_;
0178 
0179   ///B field
0180   math::XYZVector B_;
0181 
0182   /// Map used to create the TrackRef, PreIdRef value map
0183   std::map<reco::TrackRef, unsigned> refMap_;
0184 };
0185 
0186 using namespace edm;
0187 using namespace std;
0188 using namespace reco;
0189 
0190 GoodSeedProducer::GoodSeedProducer(const ParameterSet& iConfig, const goodseedhelpers::HeavyObjectCache*)
0191     : pfTransformer_(nullptr),
0192       conf_(iConfig),
0193       resMapEtaECAL_(nullptr),
0194       resMapPhiECAL_(nullptr),
0195       fitterToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Fitter")))),
0196       smootherToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Smoother")))),
0197       trackerRecHitBuilderToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("TTRHBuilder")))),
0198       magneticFieldToken_(esConsumes()),
0199       magneticFieldTokenBeginRun_(esConsumes<edm::Transition::BeginRun>()) {
0200   LogInfo("GoodSeedProducer") << "Electron PreIdentification started  ";
0201 
0202   //now do what ever initialization is needed
0203   std::vector<edm::InputTag> tags = iConfig.getParameter<vector<InputTag>>("TkColList");
0204   for (unsigned int i = 0; i < tags.size(); ++i) {
0205     trajContainers_.push_back(consumes<vector<Trajectory>>(tags[i]));
0206     tracksContainers_.push_back(consumes<reco::TrackCollection>(tags[i]));
0207   }
0208 
0209   minPt_ = iConfig.getParameter<double>("MinPt");
0210   maxPt_ = iConfig.getParameter<double>("MaxPt");
0211   maxEta_ = iConfig.getParameter<double>("MaxEta");
0212 
0213   HcalIsolWindow_ = iConfig.getParameter<double>("HcalWindow");
0214   EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
0215   EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
0216   EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
0217   EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
0218   minEoverP_ = iConfig.getParameter<double>("EOverPLead_minValue");
0219   maxHoverP_ = iConfig.getParameter<double>("HOverPLead_maxValue");
0220 
0221   pfCLusTagECLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFEcalClusterLabel"));
0222 
0223   pfCLusTagHCLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFHcalClusterLabel"));
0224 
0225   pfCLusTagPSLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFPSClusterLabel"));
0226 
0227   preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
0228   preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
0229   preidname_ = iConfig.getParameter<string>("PreIdLabel");
0230 
0231   fitterName_ = iConfig.getParameter<string>("Fitter");
0232   smootherName_ = iConfig.getParameter<string>("Smoother");
0233 
0234   nHitsInSeed_ = iConfig.getParameter<int>("NHitsInSeed");
0235 
0236   clusThreshold_ = iConfig.getParameter<double>("ClusterThreshold");
0237 
0238   minEp_ = iConfig.getParameter<double>("MinEOverP");
0239   maxEp_ = iConfig.getParameter<double>("MaxEOverP");
0240 
0241   //collection to produce
0242   produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed", false);
0243 
0244   // to disable the electron part (for HI collisions for examples)
0245   disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId", false);
0246 
0247   producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId", true);
0248   //  if no electron, cannot produce the preid
0249   if (disablePreId_)
0250     producePreId_ = false;
0251   PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId", 1.);
0252 
0253   LogDebug("GoodSeedProducer") << "Seeds for GSF will be produced ";
0254 
0255   // no disablePreId_ switch here. The collection will be empty if it is true
0256   produces<ElectronSeedCollection>(preidgsf_);
0257 
0258   if (produceCkfseed_) {
0259     LogDebug("GoodSeedProducer") << "Seeds for CKF will be produced ";
0260     produces<TrajectorySeedCollection>(preidckf_);
0261   }
0262 
0263   if (producePreId_) {
0264     LogDebug("GoodSeedProducer") << "PreId debugging information will be produced ";
0265 
0266     produces<PreIdCollection>(preidname_);
0267     if (tracksContainers_.size() == 1)  // do not make a value map if more than one input track collection
0268       produces<edm::ValueMap<reco::PreIdRef>>(preidname_);
0269   }
0270 
0271   useQuality_ = iConfig.getParameter<bool>("UseQuality");
0272   trackQuality_ = TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
0273 
0274   useTmva_ = iConfig.getUntrackedParameter<bool>("UseTMVA", false);
0275 
0276   Min_dr_ = iConfig.getParameter<double>("Min_dr");
0277 
0278   trackerRecHitBuilderName_ = iConfig.getParameter<std::string>("TTRHBuilder");
0279 }
0280 
0281 //
0282 // member functions
0283 //
0284 
0285 // ------------ method called to produce the data  ------------
0286 void GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0287   LogDebug("GoodSeedProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run();
0288   //Create empty output collections
0289   auto output_preid = std::make_unique<ElectronSeedCollection>();
0290   auto output_nopre = std::make_unique<TrajectorySeedCollection>();
0291   auto output_preidinfo = std::make_unique<PreIdCollection>();
0292   auto preIdMap_p = std::make_unique<edm::ValueMap<reco::PreIdRef>>();
0293   edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);
0294 
0295   std::unique_ptr<TrajectoryFitter> fitter;
0296   std::unique_ptr<TrajectorySmoother> smoother;
0297 
0298   //Tracking Tools
0299   if (!disablePreId_) {
0300     auto const& aFitter = &iSetup.getData(fitterToken_);
0301     auto const& aSmoother = &iSetup.getData(smootherToken_);
0302 
0303     smoother.reset(aSmoother->clone());
0304     fitter = aFitter->clone();
0305     auto const& theTrackerRecHitBuilder = &iSetup.getData(trackerRecHitBuilderToken_);
0306     hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder)->cloner();
0307     fitter->setHitCloner(&hitCloner);
0308     smoother->setHitCloner(&hitCloner);
0309   }
0310 
0311   // clear temporary maps
0312   refMap_.clear();
0313 
0314   //Magnetic Field
0315   auto const& magneticField = &iSetup.getData(magneticFieldToken_);
0316 
0317   //Handle input collections
0318   //ECAL clusters
0319   Handle<PFClusterCollection> theECPfClustCollection;
0320   iEvent.getByToken(pfCLusTagECLabel_, theECPfClustCollection);
0321 
0322   vector<PFCluster const*> basClus;
0323   for (auto const& klus : *theECPfClustCollection.product()) {
0324     if (klus.correctedEnergy() > clusThreshold_)
0325       basClus.push_back(&klus);
0326   }
0327 
0328   //HCAL clusters
0329   Handle<PFClusterCollection> theHCPfClustCollection;
0330   iEvent.getByToken(pfCLusTagHCLabel_, theHCPfClustCollection);
0331 
0332   //PS clusters
0333   Handle<PFClusterCollection> thePSPfClustCollection;
0334   iEvent.getByToken(pfCLusTagPSLabel_, thePSPfClustCollection);
0335 
0336   //Vector of track collections
0337   for (unsigned int istr = 0; istr < tracksContainers_.size(); ++istr) {
0338     //Track collection
0339     Handle<TrackCollection> tkRefCollection;
0340     iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
0341     const TrackCollection& Tk = *(tkRefCollection.product());
0342 
0343     LogDebug("GoodSeedProducer") << "Number of tracks in collection "
0344                                  << "tracksContainers_[" << istr << "] to be analyzed " << Tk.size();
0345 
0346     //loop over the track collection
0347     for (unsigned int i = 0; i < Tk.size(); ++i) {
0348       if (useQuality_ && (!(Tk[i].quality(trackQuality_))))
0349         continue;
0350 
0351       reco::PreId myPreId;
0352       bool GoodPreId = false;
0353 
0354       TrackRef trackRef(tkRefCollection, i);
0355       math::XYZVectorF tkmom(Tk[i].momentum());
0356       auto tketa = tkmom.eta();
0357       auto tkpt = std::sqrt(tkmom.perp2());
0358       auto const& Seed = (*trackRef->seedRef());
0359 
0360       if (!disablePreId_) {
0361         int ipteta = getBin(Tk[i].eta(), Tk[i].pt());
0362         int ibin = ipteta * 9;
0363 
0364         // FIXME the original code was buggy should be outerMomentum...
0365         float oPTOB = 1.f / std::sqrt(Tk[i].innerMomentum().mag2());
0366         //  float chikfred=Tk[i].normalizedChi2();
0367         float nchi = Tk[i].normalizedChi2();
0368 
0369         int nhitpi = Tk[i].found();
0370         float EP = 0;
0371 
0372         // set track info
0373         myPreId.setTrack(trackRef);
0374         //CLUSTERS - TRACK matching
0375 
0376         auto pfmass = 0.0005;
0377         auto pfoutenergy = sqrt((pfmass * pfmass) + Tk[i].outerMomentum().Mag2());
0378 
0379         XYZTLorentzVector mom = XYZTLorentzVector(
0380             Tk[i].outerMomentum().x(), Tk[i].outerMomentum().y(), Tk[i].outerMomentum().z(), pfoutenergy);
0381         XYZTLorentzVector pos =
0382             XYZTLorentzVector(Tk[i].outerPosition().x(), Tk[i].outerPosition().y(), Tk[i].outerPosition().z(), 0.);
0383 
0384         BaseParticlePropagator theOutParticle(RawParticle(mom, pos, Tk[i].charge()), 0, 0, B_.z());
0385 
0386         theOutParticle.propagateToEcalEntrance(false);
0387 
0388         float toteta = 1000.f;
0389         float totphi = 1000.f;
0390         float dr = 1000.f;
0391         float EE = 0.f;
0392         float feta = 0.f;
0393         GlobalPoint ElecTrkEcalPos(0, 0, 0);
0394 
0395         PFClusterRef clusterRef;
0396         math::XYZPoint meanShowerSaved;
0397         if (theOutParticle.getSuccess() != 0) {
0398           ElecTrkEcalPos = GlobalPoint(theOutParticle.particle().vertex().x(),
0399                                        theOutParticle.particle().vertex().y(),
0400                                        theOutParticle.particle().vertex().z());
0401 
0402           constexpr float psLim = 2.50746495928f;  // std::sinh(1.65f);
0403           bool isBelowPS = (ElecTrkEcalPos.z() * ElecTrkEcalPos.z()) > (psLim * psLim) * ElecTrkEcalPos.perp2();
0404           // bool isBelowPS=(std::abs(ElecTrkEcalPos.eta())>1.65f);
0405 
0406           unsigned clusCounter = 0;
0407           float max_ee = 0;
0408           for (auto aClus : basClus) {
0409             float tmp_ep = float(aClus->correctedEnergy()) * oPTOB;
0410             if ((tmp_ep < minEp_) | (tmp_ep > maxEp_)) {
0411               ++clusCounter;
0412               continue;
0413             }
0414 
0415             double ecalShowerDepth = PFCluster::getDepthCorrection(aClus->correctedEnergy(), isBelowPS, false);
0416             auto mom = theOutParticle.particle().momentum().Vect();
0417             auto meanShower = ElecTrkEcalPos + GlobalVector(mom.x(), mom.y(), mom.z()).unit() * ecalShowerDepth;
0418 
0419             float etarec = meanShower.eta();
0420             float phirec = meanShower.phi();
0421 
0422             float tmp_phi = std::abs(aClus->positionREP().phi() - phirec);
0423             if (tmp_phi > float(TMath::Pi()))
0424               tmp_phi -= float(TMath::TwoPi());
0425 
0426             float tmp_dr = std::sqrt(std::pow(tmp_phi, 2.f) + std::pow(aClus->positionREP().eta() - etarec, 2.f));
0427 
0428             if (tmp_dr < dr) {
0429               dr = tmp_dr;
0430               if (dr < Min_dr_) {  // find the most closest and energetic ECAL cluster
0431                 if (aClus->correctedEnergy() > max_ee) {
0432                   toteta = aClus->positionREP().eta() - etarec;
0433                   totphi = tmp_phi;
0434                   EP = tmp_ep;
0435                   EE = aClus->correctedEnergy();
0436                   feta = aClus->positionREP().eta();
0437                   clusterRef = PFClusterRef(theECPfClustCollection, clusCounter);
0438                   meanShowerSaved = meanShower;
0439                 }
0440               }
0441             }
0442             ++clusCounter;
0443           }
0444         }
0445         float trk_ecalDeta_ = fabs(toteta);
0446         float trk_ecalDphi_ = fabs(totphi);
0447 
0448         //Resolution maps
0449         auto ecaletares = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta, EE));
0450         auto ecalphires = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta, EE));
0451 
0452         //geomatrical compatibility
0453         float chieta = (toteta != 1000.f) ? toteta / ecaletares : toteta;
0454         float chiphi = (totphi != 1000.f) ? totphi / ecalphires : totphi;
0455         float chichi = sqrt(chieta * chieta + chiphi * chiphi);
0456 
0457         //Matching criteria
0458         float eta_cut = thr[ibin + 0];
0459         float phi_cut = thr[ibin + 1];
0460         float ep_cutmin = thr[ibin + 2];
0461         bool GoodMatching =
0462             ((trk_ecalDeta_ < eta_cut) && (trk_ecalDphi_ < phi_cut) && (EP > ep_cutmin) && (nhitpi > 10));
0463 
0464         bool EcalMatching = GoodMatching;
0465 
0466         if (tkpt > maxPt_)
0467           GoodMatching = true;
0468         if (tkpt < minPt_)
0469           GoodMatching = false;
0470 
0471         math::XYZPoint myPoint(ElecTrkEcalPos.x(), ElecTrkEcalPos.y(), ElecTrkEcalPos.z());
0472         myPreId.setECALMatchingProperties(
0473             clusterRef, myPoint, meanShowerSaved, std::abs(toteta), std::abs(totphi), chieta, chiphi, chichi, EP);
0474         myPreId.setECALMatching(EcalMatching);
0475 
0476         bool GoodRange = ((std::abs(tketa) < maxEta_) & (tkpt > minPt_));
0477         //KF FILTERING FOR UNMATCHED EVENTS
0478         int hit1max = int(thr[ibin + 3]);
0479         float chiredmin = thr[ibin + 4];
0480         bool GoodKFFiltering = ((nchi > chiredmin) | (nhitpi < hit1max));
0481 
0482         myPreId.setTrackFiltering(GoodKFFiltering);
0483 
0484         bool GoodTkId = false;
0485 
0486         if ((!GoodMatching) && (GoodKFFiltering) && (GoodRange)) {
0487           chired = 1000;
0488           chiRatio = 1000;
0489           dpt = 0;
0490           nhit = nhitpi;
0491           chikfred = nchi;
0492           trk_ecalDeta = trk_ecalDeta_;
0493           trk_ecalDphi = trk_ecalDphi_;
0494 
0495           Trajectory::ConstRecHitContainer tmp;
0496           for (auto const& hit : Tk[i].recHits())
0497             tmp.push_back(hit->cloneSH());
0498           auto const& theTrack = Tk[i];
0499           GlobalVector gv(theTrack.innerMomentum().x(), theTrack.innerMomentum().y(), theTrack.innerMomentum().z());
0500           GlobalPoint gp(theTrack.innerPosition().x(), theTrack.innerPosition().y(), theTrack.innerPosition().z());
0501           GlobalTrajectoryParameters gtps(gp, gv, theTrack.charge(), &*magneticField);
0502           TrajectoryStateOnSurface tsos(gtps, theTrack.innerStateCovariance(), *tmp[0]->surface());
0503           Trajectory&& FitTjs = fitter->fitOne(Seed, tmp, tsos);
0504 
0505           if (FitTjs.isValid()) {
0506             Trajectory&& SmooTjs = smoother->trajectory(FitTjs);
0507             if (SmooTjs.isValid()) {
0508               //Track refitted with electron hypothesis
0509 
0510               float pt_out = SmooTjs.firstMeasurement().updatedState().globalMomentum().perp();
0511               float pt_in = SmooTjs.lastMeasurement().updatedState().globalMomentum().perp();
0512               dpt = (pt_in > 0) ? fabs(pt_out - pt_in) / pt_in : 0.;
0513               // the following is simply the number of degrees of freedom
0514               chiRatio = SmooTjs.chiSquared() / Tk[i].chi2();
0515               chired = chiRatio * chikfred;
0516             }
0517           }
0518 
0519           //TMVA Analysis
0520           if (useTmva_) {
0521             float vars[10] = {nhit, chikfred, dpt, EP, chiRatio, chired, trk_ecalDeta, trk_ecalDphi, tkpt, tketa};
0522 
0523             float Ytmva = globalCache()->gbr[ipteta]->GetClassifier(vars);
0524 
0525             float BDTcut = thr[ibin + 5];
0526             if (Ytmva > BDTcut)
0527               GoodTkId = true;
0528             myPreId.setMVA(GoodTkId, Ytmva);
0529             myPreId.setTrackProperties(chired, chiRatio, dpt);
0530           } else {
0531             float chiratiocut = thr[ibin + 6];
0532             float gschicut = thr[ibin + 7];
0533             float gsptmin = thr[ibin + 8];
0534 
0535             GoodTkId = ((dpt > gsptmin) & (chired < gschicut) & (chiRatio < chiratiocut));
0536           }
0537         }
0538 
0539         GoodPreId = GoodTkId | GoodMatching;
0540 
0541         myPreId.setFinalDecision(GoodPreId);
0542 
0543 #ifdef EDM_ML_DEBUG
0544         if (GoodPreId)
0545           LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
0546                                        << ") preidentified for agreement between  track and ECAL cluster";
0547         if (GoodPreId && (!GoodMatching))
0548           LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
0549                                        << ") preidentified only for track properties";
0550 #endif
0551 
0552       }  // end of !disablePreId_
0553 
0554       if (GoodPreId) {
0555         //NEW SEED with n hits
0556         ElectronSeed NewSeed(Seed);
0557         NewSeed.setCtfTrack(trackRef);
0558         output_preid->push_back(NewSeed);
0559       } else {
0560         if (produceCkfseed_) {
0561           output_nopre->push_back(Seed);
0562         }
0563       }
0564       if (producePreId_ && myPreId.pt() > PtThresholdSavePredId_) {
0565         // save the index of the PreId object as to be able to create a Ref later
0566         refMap_[trackRef] = output_preidinfo->size();
0567         output_preidinfo->push_back(myPreId);
0568       }
0569     }  //end loop on track collection
0570   }    //end loop on the vector of track collections
0571 
0572   // no disablePreId_ switch, it is simpler to have an empty collection rather than no collection
0573   iEvent.put(std::move(output_preid), preidgsf_);
0574   if (produceCkfseed_)
0575     iEvent.put(std::move(output_nopre), preidckf_);
0576   if (producePreId_) {
0577     const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(std::move(output_preidinfo), preidname_);
0578     // now make the Value Map, but only if one input collection
0579     if (tracksContainers_.size() == 1) {
0580       Handle<TrackCollection> tkRefCollection;
0581       iEvent.getByToken(tracksContainers_[0], tkRefCollection);
0582       fillPreIdRefValueMap(tkRefCollection, preIdRefProd, mapFiller);
0583       mapFiller.fill();
0584       iEvent.put(std::move(preIdMap_p), preidname_);
0585     }
0586   }
0587   // clear temporary maps
0588   refMap_.clear();
0589 }
0590 
0591 // intialize the cross-thread cache to hold gbr trees and resolution maps
0592 namespace goodseedhelpers {
0593   HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& conf) {
0594     // mvas
0595     const bool useTmva = conf.getUntrackedParameter<bool>("UseTMVA", false);
0596 
0597     if (useTmva) {
0598       std::array<edm::FileInPath, kMaxWeights> weights = {{edm::FileInPath(conf.getParameter<string>("Weights1")),
0599                                                            edm::FileInPath(conf.getParameter<string>("Weights2")),
0600                                                            edm::FileInPath(conf.getParameter<string>("Weights3")),
0601                                                            edm::FileInPath(conf.getParameter<string>("Weights4")),
0602                                                            edm::FileInPath(conf.getParameter<string>("Weights5")),
0603                                                            edm::FileInPath(conf.getParameter<string>("Weights6")),
0604                                                            edm::FileInPath(conf.getParameter<string>("Weights7")),
0605                                                            edm::FileInPath(conf.getParameter<string>("Weights8")),
0606                                                            edm::FileInPath(conf.getParameter<string>("Weights9"))}};
0607 
0608       for (UInt_t j = 0; j < gbr.size(); ++j) {
0609         gbr[j] = createGBRForest(weights[j]);
0610       }
0611     }
0612   }
0613 }  // namespace goodseedhelpers
0614 
0615 // ------------ method called once each job just before starting event loop  ------------
0616 void GoodSeedProducer::beginRun(const edm::Run& run, const EventSetup& es) {
0617   //Magnetic Field
0618   auto const& magneticField = &es.getData(magneticFieldTokenBeginRun_);
0619   B_ = magneticField->inTesla(GlobalPoint(0, 0, 0));
0620 
0621   pfTransformer_ = std::make_unique<PFTrackTransformer>(B_);
0622   pfTransformer_->OnlyProp();
0623 
0624   //Resolution maps
0625   FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
0626   FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
0627   resMapEtaECAL_ = std::make_unique<PFResolutionMap>("ECAL_eta", ecalEtaMap.fullPath().c_str());
0628   resMapPhiECAL_ = std::make_unique<PFResolutionMap>("ECAL_phi", ecalPhiMap.fullPath().c_str());
0629 
0630   //read threshold
0631   FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
0632   ifstream ifs(parFile.fullPath().c_str());
0633   for (int iy = 0; iy < 81; ++iy)
0634     ifs >> thr[iy];
0635 }
0636 
0637 int GoodSeedProducer::getBin(float eta, float pt) {
0638   int ie = 0;
0639   int ip = 0;
0640   if (fabs(eta) < 0.8)
0641     ie = 0;
0642   else {
0643     if (fabs(eta) < 1.479)
0644       ie = 1;
0645     else
0646       ie = 2;
0647   }
0648   if (pt < 6)
0649     ip = 0;
0650   else {
0651     if (pt < 12)
0652       ip = 1;
0653     else
0654       ip = 2;
0655   }
0656   int iep = ie * 3 + ip;
0657   LogDebug("GoodSeedProducer") << "Track pt =" << pt << " eta=" << eta << " bin=" << iep;
0658   return iep;
0659 }
0660 
0661 void GoodSeedProducer::fillPreIdRefValueMap(Handle<TrackCollection> tracks,
0662                                             const edm::OrphanHandle<reco::PreIdCollection>& preidhandle,
0663                                             edm::ValueMap<reco::PreIdRef>::Filler& filler) {
0664   std::vector<reco::PreIdRef> values;
0665 
0666   unsigned ntracks = tracks->size();
0667   for (unsigned itrack = 0; itrack < ntracks; ++itrack) {
0668     reco::TrackRef theTrackRef(tracks, itrack);
0669     std::map<reco::TrackRef, unsigned>::const_iterator itcheck = refMap_.find(theTrackRef);
0670     if (itcheck == refMap_.end()) {
0671       // the track has been early discarded
0672       values.push_back(reco::PreIdRef());
0673     } else {
0674       edm::Ref<reco::PreIdCollection> preIdRef(preidhandle, itcheck->second);
0675       values.push_back(preIdRef);
0676       //     std::cout << " Checking Refs " << (theTrackRef==preIdRef->trackRef()) << std::endl;
0677     }
0678   }
0679   filler.insert(tracks, values.begin(), values.end());
0680 }
0681 
0682 DEFINE_FWK_MODULE(GoodSeedProducer);
[ 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!