Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-22 02:24:08 

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