OfflineValidation/plugins/EopElecTreeWriter.cc

0001 // -*- C++ -*- //
0002 // Package:    EopTreeWriter
0003 // Class:      EopTreeWriter
0004 //
0005 /**\class EopTreeWriter EopTreeWriter.cc Alignment/OfflineValidation/plugins/EopTreeWriter.cc
0006  //
0007  Description: <one line class summary>
0008
0009  Implementation:
0010  <Notes on implementation>
0011 */
0012 //
0013 // Original Author:  Holger Enderle
0014 //         Created:  Thu Dec  4 11:22:48 CET 2008
0015 // $Id: EopElecTreeWriter.cc,v 1.8 2012/08/30 08:40:51 cgoetzma Exp $
0016 //
0017 //
0018
0019 // framework include files
0020 #include "Alignment/OfflineValidation/interface/EopElecVariables.h"
0021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0022 #include "CommonTools/Utils/interface/TFileDirectory.h"
0023 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0024 #include "DataFormats/Common/interface/Handle.h"
0025 #include "DataFormats/Common/interface/TriggerResults.h"
0026 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0027 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0028 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0029 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0030 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0031 #include "DataFormats/EgammaCandidates/interface/GsfElectronCore.h"
0032 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0033 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0034 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0035 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0036 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0037 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0038 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0039 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0040 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0041 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0042 #include "DataFormats/JetReco/interface/CaloJet.h"
0043 #include "DataFormats/Math/interface/deltaR.h"
0044 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0045 #include "DataFormats/RecoCandidate/interface/RecoCaloTowerCandidate.h"
0046 #include "DataFormats/TrackReco/interface/Track.h"
0047 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0048 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0049 #include "DataFormats/VertexReco/interface/Vertex.h"
0050 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0051 #include "FWCore/Framework/interface/ESHandle.h"
0052 #include "FWCore/Framework/interface/Event.h"
0053 #include "FWCore/Framework/interface/EventSetup.h"
0054 #include "FWCore/Framework/interface/Frameworkfwd.h"
0055 #include "FWCore/Framework/interface/MakerMacros.h"
0056 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0057 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0058 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0059 #include "FWCore/ServiceRegistry/interface/Service.h"
0060 #include "FWCore/Utilities/interface/EDGetToken.h"
0061 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0062 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0063 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0064 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0065 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0066 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0067 #include "MagneticField/Engine/interface/MagneticField.h"
0068 #include "RecoEcal/EgammaCoreTools/interface/SuperClusterShapeAlgo.h"
0069 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0070 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0071 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0072 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
0073 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
0074 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0075 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0076
0077 #include "TCanvas.h"
0078 #include "TH1.h"
0079 #include "TH2D.h"
0080 #include "TLorentzVector.h"
0081 #include "TMath.h"
0082 #include "TTree.h"
0083
0084 struct EopTriggerType {
0085   bool fired;
0086   double prescale;
0087   int index;
0088
0089   EopTriggerType() {
0090     fired = false;
0091     prescale = 0;
0092     index = 0;
0093   }
0094 };
0095
0096 //
0097 // class decleration
0098 //
0099 using namespace std;
0100
0101 class EopElecTreeWriter : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0102 public:
0103   explicit EopElecTreeWriter(const edm::ParameterSet&);
0104   ~EopElecTreeWriter() override;
0105
0106   static void fillDescriptions(edm::ConfigurationDescriptions&);
0107
0108 private:
0109   // methods
0110   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0111   void endRun(const edm::Run&, const edm::EventSetup&) override {}
0112   void analyze(const edm::Event&, const edm::EventSetup&) override;
0113   void endJob() override;
0114
0115   // ----------member data ---------------------------
0116   // ES tokens
0117   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0118   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0119   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0120
0121   // EDM tokens
0122   const edm::EDGetTokenT<reco::VertexCollection> theVertexCollectionToken_;
0123   const edm::EDGetTokenT<HBHERecHitCollection> theHBHERecHitCollectionToken_;
0124   const edm::EDGetTokenT<EcalRecHitCollection> theEcalRecHitCollectionToken_;
0125   const edm::EDGetTokenT<reco::SuperClusterCollection> theBarrelSupClusCollectionToken_;
0126   const edm::EDGetTokenT<reco::SuperClusterCollection> theEndCapSupClusCollectionToken_;
0127   const edm::EDGetTokenT<edm::TriggerResults> theTriggerResultsToken_;
0128   const edm::EDGetTokenT<trigger::TriggerEvent> theTriggerEventToken_;
0129   const edm::EDGetTokenT<reco::GsfTrackCollection> theGsfTrackCollectionToken_;
0130   const edm::EDGetTokenT<reco::GsfElectronCoreCollection> theGsfElectronCoreCollectionToken_;
0131
0132   // some static constants
0133   static constexpr float k_etaBarrel = 1.55;
0134   static constexpr float k_etaEndcap = 1.44;
0135
0136   // other member data
0137   std::string theTrigger_;
0138   std::string theFilter_;
0139   bool debugTriggerSelection_;
0140   edm::Service<TFileService> fs_;
0141   TTree* tree_;
0142   EopElecVariables* treeMemPtr_;
0143   HLTConfigProvider hltConfig_;
0144   std::vector<std::string> triggerNames_;
0145   MultiTrajectoryStateTransform* mtsTransform_;
0146
0147   // histograms
0148
0149   // Cut flow (events number)
0150   TH1D* h_nEvents;
0151   TH1D* h_nEventsWithVertex;
0152   TH1D* h_nEventsTriggered;
0153   TH1D* h_nEventsHLTFilter;
0154   TH1D* h_nEventsHLTelectron;
0155   TH1D* h_nEventsHLTrejected;
0156   TH1D* h_nEvents2Elec;
0157   TH1D* h_nHLTelectrons;
0158   TH1D* h_nTrkRejectedPerEvt;
0159   TH1D* h_nTrkSelectedPerEvt;
0160
0161   // Cut flow (tracks number)
0162   TH1D* h_nTracks;
0163   TH1D* h_nTracksFiltered;
0164   TH1D* h_cut_Ptmin;
0165   TH1D* h_cut_OneSCmatch;
0166
0167   TH1D* h_counter1;
0168   TH1D* h_counter2;
0169
0170   TH1D* h_distToClosestSCgsf;
0171   TH1D* h_distToClosestSC;
0172   TH1D* h_EcalEnergy;
0173   TH1D* h_Momentum;
0174   TH1D* HcalEnergy;
0175   TH1D* h_fBREM;
0176   TH1D* h_Eop_InnerNegative;
0177   TH1D* h_Eop_InnerPositive;
0178
0179   TH2D* HcalVSEcal;
0180 };
0181
0182 namespace eopUtils {
0183
0184   static constexpr float R_ECAL = 136.5;
0185   static constexpr float Z_Endcap = 328.0;
0186   static constexpr float etaBarrelEndcap = 1.479;
0187
0188   // Function to convert the eta of the track to a detector eta (for matching with SC)
0189   float ecalEta(float EtaParticle, float Zvertex, float RhoVertex) {
0190     if (EtaParticle != 0.) {
0191       float Theta = 0.0;
0192       float ZEcal = (R_ECAL - RhoVertex) * sinh(EtaParticle) + Zvertex;
0193
0194       if (ZEcal != 0.0)
0195         Theta = atan(R_ECAL / ZEcal);
0196       if (Theta < 0.0)
0197         Theta = Theta + M_PI;
0198
0199       float ETA = -log(tan(0.5 * Theta));
0200
0201       if (fabs(ETA) > etaBarrelEndcap) {
0202         float Zend = Z_Endcap;
0203         if (EtaParticle < 0.0)
0204           Zend = -Zend;
0205         float Zlen = Zend - Zvertex;
0206         float RR = Zlen / sinh(EtaParticle);
0207         Theta = atan((RR + RhoVertex) / Zend);
0208         if (Theta < 0.0)
0209           Theta = Theta + M_PI;
0210         ETA = -log(tan(0.5 * Theta));
0211       }
0212       return ETA;
0213     } else {
0214       edm::LogWarning("") << "[EcalPositionFromTrack::etaTransformation] Warning: Eta equals to zero, not correcting";
0215       return EtaParticle;
0216     }
0217   }
0218 }  // namespace eopUtils
0219
0220 // constructors and destructor
0221
0222 EopElecTreeWriter::EopElecTreeWriter(const edm::ParameterSet& iConfig)
0223     : magFieldToken_(esConsumes()),
0224       tkGeomToken_(esConsumes()),
0225       caloGeomToken_(esConsumes()),
0226       theVertexCollectionToken_(consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"))),
0227       theHBHERecHitCollectionToken_(consumes<HBHERecHitCollection>(edm::InputTag("hbhereco", ""))),
0228       theEcalRecHitCollectionToken_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"))),
0229       theBarrelSupClusCollectionToken_(
0230           consumes<reco::SuperClusterCollection>(edm::InputTag("hybridSuperClusters", ""))),
0231       theEndCapSupClusCollectionToken_(consumes<reco::SuperClusterCollection>(
0232           edm::InputTag("multi5x5SuperClusters", "multi5x5EndcapSuperClusters"))),
0233       theTriggerResultsToken_(consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", "HLT"))),
0234       theTriggerEventToken_(consumes<trigger::TriggerEvent>(edm::InputTag("hltTriggerSummaryAOD"))),
0235       theGsfTrackCollectionToken_(consumes<reco::GsfTrackCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0236       theGsfElectronCoreCollectionToken_(
0237           consumes<reco::GsfElectronCoreCollection>(edm::InputTag("gedGsfElectronCores"))),
0238       theTrigger_(iConfig.getParameter<std::string>("triggerPath")),
0239       theFilter_(iConfig.getParameter<std::string>("hltFilter")),
0240       debugTriggerSelection_(iConfig.getParameter<bool>("debugTriggerSelection")) {
0241   usesResource(TFileService::kSharedResource);
0242
0243   // TTree creation
0244   tree_ = fs_->make<TTree>("EopTree", "EopTree");
0245   treeMemPtr_ = new EopElecVariables;
0246   tree_->Branch("EopElecVariables", &treeMemPtr_);  // address of pointer!
0247
0248   // Control histograms declaration
0249   h_distToClosestSC = fs_->make<TH1D>("distToClosestSC", "distToClosestSC", 100, 0, 0.1);
0250   h_distToClosestSCgsf = fs_->make<TH1D>("distToClosestSCgsf", "distToClosestSCgsf", 100, 0, 0.1);
0251   h_EcalEnergy = fs_->make<TH1D>("EcalEnergy", "EcalEnergy", 100, 0, 200);
0252   h_Momentum = fs_->make<TH1D>("Momentum", "Momentum", 100, 0, 200);
0253   HcalEnergy = fs_->make<TH1D>("HcalEnergy", "HcalEnergy", 100, 0, 40);
0254   h_fBREM = fs_->make<TH1D>("fBREM", "fBREM", 100, -0.2, 1);
0255   h_Eop_InnerNegative = fs_->make<TH1D>("Eop_InnerNegative", "Eop_InnerNegative", 100, 0, 3);
0256   h_Eop_InnerPositive = fs_->make<TH1D>("Eop_InnerPositive", "Eop_InnerPositive", 100, 0, 3);
0257   HcalVSEcal = fs_->make<TH2D>("HcalVSEcal", "HcalVSEcal", 100, 0, 160, 100, 0, 10);
0258
0259   h_nEvents = fs_->make<TH1D>("nEvents", "nEvents", 1, 0, 1);
0260   h_nEventsWithVertex = fs_->make<TH1D>("nEventsWithVertex", "nEventsWithVertex", 1, 0, 1);
0261   h_nEventsTriggered = fs_->make<TH1D>("nEventsTriggered", "nEventsTriggered", 1, 0, 1);
0262   h_nEventsHLTFilter = fs_->make<TH1D>("nEventsHLTFilter", "nEventsHLTFilter", 1, 0, 1);
0263   h_nEventsHLTelectron = fs_->make<TH1D>("nEventsHLTelectron", "nEventsHLTelectron", 1, 0, 1);
0264   h_nEventsHLTrejected = fs_->make<TH1D>("nEventsHLTrejected", "nEventsHLTrejected", 1, 0, 1);
0265   h_nEvents2Elec = fs_->make<TH1D>("nEvents2Elec", "nEvents2Elec", 1, 0, 1);
0266
0267   h_nHLTelectrons = fs_->make<TH1D>("nHLTelectrons", "nHLTelectrons", 20, 0, 20);
0268   h_nTrkRejectedPerEvt = fs_->make<TH1D>("nTrkRejectedPerEvt", "nTrkRejectedPerEvt", 20, 0, 20);
0269   h_nTrkSelectedPerEvt = fs_->make<TH1D>("nTrkSelectedPerEvt", "nTrkSelectedPerEvt", 20, 0, 20);
0270
0271   h_nTracks = fs_->make<TH1D>("nTracks", "nTracks", 1, 0, 1);
0272   h_nTracksFiltered = fs_->make<TH1D>("nTracksFiltered", "nTracksFiltered", 1, 0, 1);
0273   h_cut_Ptmin = fs_->make<TH1D>("cut_Ptmin", "cut_Ptmin", 1, 0, 1);
0274   h_cut_OneSCmatch = fs_->make<TH1D>("cut_OneSCmatch", "cut_OneSCmatch", 1, 0, 1);
0275
0276   h_counter1 = fs_->make<TH1D>("counter1", "counter1", 1, 0, 1);
0277   h_counter2 = fs_->make<TH1D>("counter2", "counter2", 1, 0, 1);
0278 }
0279
0280 EopElecTreeWriter::~EopElecTreeWriter() {
0281   // control histograms
0282   h_distToClosestSC->SetXTitle("distance from track to closest SuperCluster in eta-phi plan (weighted matching)");
0283   h_distToClosestSC->SetYTitle("# Tracks");
0284
0285   h_distToClosestSCgsf->SetXTitle("distance from track to closest SuperCluster in eta-phi plan (gsfElectronCore)");
0286   h_distToClosestSCgsf->SetYTitle("# Tracks");
0287
0288   h_EcalEnergy->SetXTitle("Ecal energy deposit (GeV)");
0289   h_EcalEnergy->SetYTitle("# tracks");
0290
0291   HcalEnergy->SetXTitle("Hcal energy deposit (GeV)");
0292   HcalEnergy->SetYTitle("# tracks");
0293
0294   h_Momentum->SetXTitle("Momentum magnitude (GeV)");
0295   h_Momentum->SetYTitle("# tracks");
0296
0297   h_Eop_InnerNegative->SetXTitle("E/p");
0298   h_Eop_InnerNegative->SetYTitle("# tracks");
0299
0300   h_Eop_InnerPositive->SetXTitle("E/p");
0301   h_Eop_InnerPositive->SetYTitle("# tracks");
0302
0303   HcalVSEcal->SetXTitle("Ecal energy (GeV)");
0304   HcalVSEcal->SetYTitle("Hcal energy (GeV)");
0305 }
0306
0307 //###########################################
0308 //#     method called to for each event     #
0309 //###########################################
0310
0311 void EopElecTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0312   using namespace edm;
0313   h_nEvents->Fill(0.5);
0314
0315   Double_t EnergyHcalIn01;
0316   Double_t EnergyHcalIn02;
0317   Double_t EnergyHcalIn03;
0318   Double_t EnergyHcalIn04;
0319   Double_t EnergyHcalIn05;
0320   Double_t etaWidth;
0321   Double_t phiWidth;
0322   Int_t algo_ID;
0323   Double_t EnergyEcal;
0324   Double_t fbrem;
0325   Double_t pin;
0326   Double_t etaIn;
0327   Double_t phiIn;
0328   Double_t pout;
0329   Double_t etaOut;
0330   Double_t phiOut;
0331   Double_t MaxPtIn01;
0332   Double_t SumPtIn01;
0333   Bool_t NoTrackIn0015;
0334   Double_t MaxPtIn02;
0335   Double_t SumPtIn02;
0336   Bool_t NoTrackIn0020;
0337   Double_t MaxPtIn03;
0338   Double_t SumPtIn03;
0339   Bool_t NoTrackIn0025;
0340   Double_t MaxPtIn04;
0341   Double_t SumPtIn04;
0342   Bool_t NoTrackIn0030;
0343   Double_t MaxPtIn05;
0344   Double_t SumPtIn05;
0345   Bool_t NoTrackIn0035;
0346   Bool_t NoTrackIn0040;
0347   Double_t dRSC_first;
0348   Double_t dRSC_second;
0349   Double_t etaSC;
0350   Double_t phiSC;
0351   Int_t nbSC;        //to count the nb of SuperCluster matching with a given track
0352   Int_t nBasicClus;  //to count the nb of basic cluster in a given superCluster
0353   Bool_t isEcalDriven;
0354   Bool_t isTrackerDriven;
0355   Bool_t isBarrel;
0356   Bool_t isEndcap;
0357   Bool_t TrigTag;
0358
0359   //---------------   for GsfTrack propagation through tracker  ---------------
0360
0361   const MagneticField* magField_ = &iSetup.getData(magFieldToken_);
0362   const TrackerGeometry* trackerGeom_ = &iSetup.getData(tkGeomToken_);
0363   MultiTrajectoryStateTransform mtsTransform(trackerGeom_, magField_);
0364
0365   //---------------    Super Cluster    -----------------
0366
0367   // getting primary vertex (necessary to convert eta track to eta detector
0368   const reco::VertexCollection& vertex = iEvent.get(theVertexCollectionToken_);
0369
0370   if (vertex.empty()) {
0371     edm::LogError("EopElecTreeWriter") << "Error: no primary vertex found!";
0372     return;
0373   }
0374   const reco::Vertex& vert = vertex.front();
0375   h_nEventsWithVertex->Fill(0.5);
0376
0377   // getting calorimeter geometry
0378   const CaloGeometry* geo = &iSetup.getData(caloGeomToken_);
0379   const CaloSubdetectorGeometry* subGeo = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0380   if (subGeo == nullptr)
0381     edm::LogError("EopElecTreeWriter") << "ERROR: unable to find SubDetector geometry!!!";
0382
0383   // getting Hcal rechits
0384   const HBHERecHitCollection* HcalHits = &iEvent.get(theHBHERecHitCollectionToken_);
0385
0386   // getting Ecal rechits
0387   const EcalRecHitCollection* rhc = &iEvent.get(theEcalRecHitCollectionToken_);
0388   if (rhc == nullptr)
0389     edm::LogError("EopElecTreeWriter") << "ERROR: unable to find the EcalRecHit collection !!!";
0390
0391   // getting SuperCluster
0392   const reco::SuperClusterCollection* BarrelSupClusCollection = &iEvent.get(theBarrelSupClusCollectionToken_);
0393   const reco::SuperClusterCollection* EndcapSupClusCollection = &iEvent.get(theEndCapSupClusCollectionToken_);
0394
0395   // necessary to re-calculate phi and eta width of SuperClusters
0396   SuperClusterShapeAlgo SCShape(rhc, subGeo);
0397
0398   //---------------    Trigger   -----------------
0399   TrigTag = false;
0400   const edm::TriggerResults* trigRes = &iEvent.get(theTriggerResultsToken_);
0401
0402   // trigger event
0403   const trigger::TriggerEvent* triggerEvent = &iEvent.get(theTriggerEventToken_);
0404
0405   // our trigger table
0406   std::map<std::string, EopTriggerType> HLTpaths;
0407   for (const auto& triggerName : triggerNames_) {
0408     if (triggerName.find(theTrigger_) != 0)
0409       continue;
0410     EopTriggerType myTrigger;
0411
0412     const unsigned int prescaleSize = hltConfig_.prescaleSize();
0413     for (unsigned int ps = 0; ps < prescaleSize; ps++) {
0414       auto const prescaleValue = hltConfig_.prescaleValue<double>(ps, triggerName);
0415       if (prescaleValue != 1) {
0416         myTrigger.prescale = prescaleValue;
0417       }
0418     }
0419
0420     myTrigger.index = hltConfig_.triggerIndex(triggerName);
0421     if (myTrigger.index == -1)
0422       continue;
0423     myTrigger.fired =
0424         trigRes->wasrun(myTrigger.index) && trigRes->accept(myTrigger.index) && !trigRes->error(myTrigger.index);
0425     HLTpaths[triggerName] = myTrigger;
0426   }
0427
0428   // First cut : trigger cut
0429   std::string firstFiredPath = "";
0430   for (const auto& it : HLTpaths) {
0431     if (it.second.fired) {
0432       TrigTag = true;
0433       firstFiredPath = it.first;
0434       break;
0435     }
0436   }
0437   if (!TrigTag)
0438     return;
0439   h_nEventsTriggered->Fill(0.5);
0440
0441   // Displaying filters label from the first fired trigger
0442   // Useful for finding the good filter label
0443   std::vector<std::string> filters = hltConfig_.moduleLabels(firstFiredPath);
0444
0445   if (debugTriggerSelection_) {
0446     edm::LogInfo("EopElecTreeWriter") << "filters : ";
0447     for (unsigned int i = 0; i < filters.size(); i++) {
0448       edm::LogInfo("EopElecTreeWriter") << filters[i] << " ";
0449     }
0450   }
0451
0452   // Getting HLT electrons
0453   edm::InputTag testTag(theFilter_, "", "HLT");
0454   int testindex = triggerEvent->filterIndex(testTag);
0455
0456   if (testindex >= triggerEvent->sizeFilters())
0457     return;
0458   h_nEventsHLTFilter->Fill(0.5);
0459
0460   const trigger::Keys& KEYS_el(triggerEvent->filterKeys(testindex));
0461
0462   std::vector<const trigger::TriggerObject*> HLTelectrons;
0463   for (unsigned int i = 0; i != KEYS_el.size(); ++i) {
0464     const trigger::TriggerObject* triggerObject_el = &(triggerEvent->getObjects().at(KEYS_el[i]));
0465     HLTelectrons.push_back(triggerObject_el);
0466   }
0467
0468   h_nHLTelectrons->Fill(HLTelectrons.size());
0469
0470   if (HLTelectrons.empty())
0471     return;
0472   h_nEventsHLTelectron->Fill(0.5);
0473
0474   // finding the HLT electron with highest pt and saving the corresponding index
0475   unsigned int HighPtIndex = 0;
0476   double maxPtHLT = -5.;
0477   for (unsigned int j = 0; j < HLTelectrons.size(); j++) {
0478     if (HLTelectrons[j]->pt() > maxPtHLT) {
0479       maxPtHLT = HLTelectrons[j]->pt();
0480       HighPtIndex = j;
0481     }
0482   }
0483
0484   //-----------------   Tracks   -------------------
0485
0486   // getting GsfTrack
0487   const reco::GsfTrackCollection& tracks = iEvent.get(theGsfTrackCollectionToken_);
0488
0489   // filtering track
0490   int nRejected = 0;
0491   int nSelected = 0;
0492   std::vector<const reco::GsfTrack*> filterTracks;
0493   for (const auto& track : tracks) {
0494     h_nTracks->Fill(0.5);
0495     double deltar =
0496         reco::deltaR(track.eta(), track.phi(), HLTelectrons[HighPtIndex]->eta(), HLTelectrons[HighPtIndex]->phi());
0497     // remove the triggered electron with highest pt
0498     if (deltar < 0.025) {
0499       treeMemPtr_->px_rejected_track = track.px();
0500       treeMemPtr_->py_rejected_track = track.py();
0501       treeMemPtr_->pz_rejected_track = track.pz();
0502       treeMemPtr_->p_rejected_track = track.p();
0503       nRejected++;
0504       continue;
0505     }
0506     filterTracks.push_back(&track);  // we use all the others
0507     nSelected++;
0508     h_nTracksFiltered->Fill(0.5);
0509   }
0510   h_nTrkRejectedPerEvt->Fill(nRejected);
0511   h_nTrkSelectedPerEvt->Fill(nSelected);
0512
0513   if (nRejected == 0)
0514     return;
0515   h_nEventsHLTrejected->Fill(0.5);
0516
0517   if (filterTracks.empty())
0518     return;
0519   h_nEvents2Elec->Fill(0.5);
0520
0521   //-------- test:Matching SC/track using gsfElectonCore collection --------
0522
0523   const reco::GsfElectronCoreCollection* electrons = &iEvent.get(theGsfElectronCoreCollectionToken_);
0524   for (const auto& elec : *electrons) {
0525     double etaGSF = eopUtils::ecalEta((elec.gsfTrack())->eta(), vert.z(), (vert.position()).rho());
0526     if ((elec.gsfTrack())->pt() < 10.)
0527       continue;
0528
0529     double DELTAR = 0;
0530     DELTAR = reco::deltaR((elec.superCluster())->eta(), (elec.superCluster())->phi(), etaGSF, (elec.gsfTrack())->phi());
0531
0532     if (DELTAR < 0.1)
0533       h_distToClosestSCgsf->Fill(DELTAR);
0534   }
0535
0536   //------------------------------------------------------------
0537   //--------------    Loop on tracks   -------------------------
0538
0539   for (const auto& track : filterTracks) {
0540     // initializing variables
0541     isEcalDriven = false;
0542     isTrackerDriven = false;
0543     isBarrel = false;
0544     isEndcap = false;
0545     etaWidth = 0;
0546     phiWidth = 0;
0547     etaSC = 0;
0548     phiSC = 0;
0549     fbrem = 0;
0550     pin = 0;
0551     etaIn = 0;
0552     phiIn = 0;
0553     pout = 0;
0554     etaOut = 0;
0555     phiOut = 0;
0556     algo_ID = 0;
0557     EnergyEcal = 0;
0558     dRSC_first = 999;
0559     dRSC_second = 9999;
0560     nbSC = 0;
0561     nBasicClus = 0;
0562
0563     // First cut on momentum magnitude
0564     h_Momentum->Fill(track->p());
0565     if (track->pt() < 10.)
0566       continue;
0567     h_cut_Ptmin->Fill(0.5);
0568
0569     // calculating track parameters at innermost and outermost for Gsf tracks
0570     TrajectoryStateOnSurface inTSOS = mtsTransform.innerStateOnSurface(*track);
0571     TrajectoryStateOnSurface outTSOS = mtsTransform.outerStateOnSurface(*track);
0572     if (inTSOS.isValid() && outTSOS.isValid()) {
0573       GlobalVector inMom;
0574       //multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS, inMom_);
0575       multiTrajectoryStateMode::momentumFromModePPhiEta(inTSOS, inMom);
0576       pin = inMom.mag();
0577       etaIn = inMom.eta();
0578       phiIn = inMom.phi();
0579       GlobalVector outMom;
0580       //multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, outMom_);
0581       multiTrajectoryStateMode::momentumFromModePPhiEta(outTSOS, outMom);
0582       pout = outMom.mag();
0583       etaOut = outMom.eta();
0584       phiOut = outMom.phi();
0585       fbrem = (pin - pout) / pin;
0586       h_fBREM->Fill(fbrem);
0587     }
0588
0589     // Matching track with Hcal rec hits
0590     EnergyHcalIn01 = 0;
0591     EnergyHcalIn02 = 0;
0592     EnergyHcalIn03 = 0;
0593     EnergyHcalIn04 = 0;
0594     EnergyHcalIn05 = 0;
0595
0596     //for (std::vector<HBHERecHit>::const_iterator hcal = (*HcalHits).begin(); hcal != (*HcalHits).end(); hcal++) {
0597     for (const auto& hcal : *HcalHits) {
0598       GlobalPoint posH = geo->getPosition(hcal.detid());
0599       double phihit = posH.phi();
0600       double etahit = posH.eta();
0601       double dR = reco::deltaR(etahit, phihit, etaOut, phiOut);
0602
0603       // saving Hcal energy deposit measured for different eta-phi radius
0604       if (dR < 0.1)
0605         EnergyHcalIn01 += hcal.energy();
0606       if (dR < 0.2)
0607         EnergyHcalIn02 += hcal.energy();
0608       if (dR < 0.3)
0609         EnergyHcalIn03 += hcal.energy();
0610       if (dR < 0.4)
0611         EnergyHcalIn04 += hcal.energy();
0612       if (dR < 0.5)
0613         EnergyHcalIn05 += hcal.energy();
0614     }
0615
0616     HcalEnergy->Fill(EnergyHcalIn02);
0617
0618     //Isolation against charged particles
0619     MaxPtIn01 = 0.;
0620     SumPtIn01 = 0.;
0621     NoTrackIn0015 = true;
0622     MaxPtIn02 = 0.;
0623     SumPtIn02 = 0.;
0624     NoTrackIn0020 = true;
0625     MaxPtIn03 = 0.;
0626     SumPtIn03 = 0.;
0627     NoTrackIn0025 = true;
0628     MaxPtIn04 = 0.;
0629     SumPtIn04 = 0.;
0630     NoTrackIn0030 = true;
0631     MaxPtIn05 = 0.;
0632     SumPtIn05 = 0.;
0633     NoTrackIn0035 = true;
0634     NoTrackIn0040 = true;
0635
0636     for (const auto& track1 : filterTracks) {
0637       if (track == track1)
0638         continue;
0639
0640       double etaIn1 = 0.;
0641       double phiIn1 = 0.;
0642
0643       TrajectoryStateOnSurface inTSOS1 = mtsTransform.innerStateOnSurface(*track1);
0644       if (inTSOS1.isValid()) {
0645         GlobalVector inMom1;
0646         multiTrajectoryStateMode::momentumFromModePPhiEta(inTSOS1, inMom1);
0647         etaIn1 = inMom1.eta();
0648         phiIn1 = inMom1.phi();
0649       }
0650
0651       if (etaIn1 == 0 && phiIn1 == 0)
0652         continue;
0653
0654       double dR = reco::deltaR(etaIn1, phiIn1, etaIn, phiIn);
0655
0656       // different radius of inner isolation cone
0657       if (dR < 0.015)
0658         NoTrackIn0015 = false;
0659       if (dR < 0.020)
0660         NoTrackIn0020 = false;
0661       if (dR < 0.025)
0662         NoTrackIn0025 = false;
0663       if (dR < 0.030)
0664         NoTrackIn0030 = false;
0665       if (dR < 0.035)
0666         NoTrackIn0035 = false;
0667       if (dR < 0.040)
0668         NoTrackIn0040 = false;
0669
0670       //calculate maximum Pt and sum Pt inside cones of different radius
0671       if (dR < 0.1) {
0672         SumPtIn01 += track1->pt();
0673         if (track1->pt() > MaxPtIn01) {
0674           MaxPtIn01 = track1->pt();
0675         }
0676       }
0677
0678       if (dR < 0.2) {
0679         SumPtIn02 += track1->pt();
0680         if (track1->pt() > MaxPtIn02) {
0681           MaxPtIn02 = track1->pt();
0682         }
0683       }
0684       if (dR < 0.3) {
0685         SumPtIn03 += track1->pt();
0686         if (track1->pt() > MaxPtIn03) {
0687           MaxPtIn03 = track1->pt();
0688         }
0689       }
0690       if (dR < 0.4) {
0691         SumPtIn04 += track1->pt();
0692         if (track1->pt() > MaxPtIn04) {
0693           MaxPtIn04 = track1->pt();
0694         }
0695       }
0696       if (dR < 0.5) {
0697         SumPtIn05 += track1->pt();
0698         if (track1->pt() > MaxPtIn05) {
0699           MaxPtIn05 = track1->pt();
0700         }
0701       }
0702     }
0703
0704     // Track-SuperCluster matching
0705
0706     double dRSC;
0707     double etaAtEcal = 0;  // to convert eta from track to detector frame
0708     etaAtEcal = eopUtils::ecalEta(etaIn, vert.z(), (vert.position()).rho());
0709
0710     //Barrel
0711     if (std::abs(track->eta()) < k_etaBarrel) {
0712       for (const auto& SC : *BarrelSupClusCollection) {
0713         dRSC = reco::deltaR(SC.eta(), SC.phi(), etaAtEcal, phiIn);
0714
0715         if (dRSC < dRSC_first) {
0716           dRSC_first = dRSC;
0717         }  // distance in eta-phi plan to closest SC
0718         else if (dRSC < dRSC_second) {
0719           dRSC_second = dRSC;
0720         }  // to next closest SC
0721
0722         if (dRSC < 0.09) {
0723           //Calculate phiWidth & etaWidth for associated SuperClusters
0724           SCShape.Calculate_Covariances(SC);
0725           phiWidth = SCShape.phiWidth();
0726           etaWidth = SCShape.etaWidth();
0727
0728           etaSC = SC.eta();
0729           phiSC = SC.phi();
0730
0731           algo_ID = SC.algoID();
0732           EnergyEcal = SC.energy();
0733           nBasicClus = SC.clustersSize();
0734           nbSC++;
0735           isBarrel = true;
0736         }
0737       }
0738     }
0739
0740     // Endcap
0741     if (std::abs(track->eta()) > k_etaEndcap) {
0742       for (const auto& SC : *EndcapSupClusCollection) {
0743         dRSC = reco::deltaR(SC.eta(), SC.phi(), etaAtEcal, phiIn);
0744
0745         if (dRSC < dRSC_first) {
0746           dRSC_first = dRSC;
0747         }  // distance in eta-phi plan to closest SC
0748         else if (dRSC < dRSC_second) {
0749           dRSC_second = dRSC;
0750         }  // to next closest SC
0751
0752         if (dRSC < 0.09) {
0753           //Calculate phiWidth & etaWidth for associated SuperClusters
0754           SCShape.Calculate_Covariances(SC);
0755           phiWidth = SCShape.phiWidth();
0756           etaWidth = SCShape.etaWidth();
0757
0758           etaSC = SC.eta();
0759           phiSC = SC.phi();
0760
0761           algo_ID = SC.algoID();
0762           EnergyEcal = SC.energy();
0763           nBasicClus = SC.clustersSize();
0764           nbSC++;
0765           isEndcap = true;
0766         }
0767       }
0768     }
0769
0770     if (dRSC_first < 0.1)
0771       h_distToClosestSC->Fill(dRSC_first);
0772     if (nbSC == 1)
0773       h_counter1->Fill(0.5);
0774     if (nbSC == 0)
0775       h_counter2->Fill(0.5);
0776
0777     if (nbSC > 1 || nbSC == 0)
0778       continue;
0779     h_cut_OneSCmatch->Fill(0.5);
0780
0781     if (isBarrel && isEndcap) {
0782       edm::LogError("EopElecTreeWriter") << "Error: Super Cluster double matching!";
0783       return;
0784     }
0785
0786     h_EcalEnergy->Fill(EnergyEcal);
0787     HcalVSEcal->Fill(EnergyEcal, EnergyHcalIn02);
0788
0789     // E over p plots
0790     if (track->charge() < 0)
0791       h_Eop_InnerNegative->Fill(EnergyEcal / pin);
0792     if (track->charge() > 0)
0793       h_Eop_InnerPositive->Fill(EnergyEcal / pin);
0794
0795     //Check if track-SuperCluster matching is Ecal or Tracker driven
0796     edm::RefToBase<TrajectorySeed> seed = track->extra()->seedRef();
0797     if (seed.isNull()) {
0798       edm::LogError("GsfElectronCore") << "The GsfTrack has no seed ?!";
0799     } else {
0800       reco::ElectronSeedRef elseed = seed.castTo<reco::ElectronSeedRef>();
0801       if (elseed.isNull()) {
0802         edm::LogError("GsfElectronCore") << "The GsfTrack seed is not an ElectronSeed ?!";
0803       } else {
0804         isEcalDriven = elseed->isEcalDriven();
0805         isTrackerDriven = elseed->isTrackerDriven();
0806       }
0807     }
0808
0809     treeMemPtr_->charge = track->charge();
0810     treeMemPtr_->nHits = track->numberOfValidHits();
0811     treeMemPtr_->nLostHits = track->numberOfLostHits();
0812     treeMemPtr_->innerOk = track->innerOk();
0813     treeMemPtr_->outerRadius = track->outerRadius();
0814     treeMemPtr_->chi2 = track->chi2();
0815     treeMemPtr_->normalizedChi2 = track->normalizedChi2();
0816     treeMemPtr_->px = track->px();
0817     treeMemPtr_->py = track->py();
0818     treeMemPtr_->pz = track->pz();
0819     treeMemPtr_->p = track->p();
0820     treeMemPtr_->pIn = pin;
0821     treeMemPtr_->etaIn = etaIn;
0822     treeMemPtr_->phiIn = phiIn;
0823     treeMemPtr_->pOut = pout;
0824     treeMemPtr_->etaOut = etaOut;
0825     treeMemPtr_->phiOut = phiOut;
0826     treeMemPtr_->pt = track->pt();
0827     treeMemPtr_->ptError = track->ptError();
0828     treeMemPtr_->theta = track->theta();
0829     treeMemPtr_->eta = track->eta();
0830     treeMemPtr_->phi = track->phi();
0831     treeMemPtr_->fbrem = fbrem;
0832     treeMemPtr_->MaxPtIn01 = MaxPtIn01;
0833     treeMemPtr_->SumPtIn01 = SumPtIn01;
0834     treeMemPtr_->NoTrackIn0015 = NoTrackIn0015;
0835     treeMemPtr_->MaxPtIn02 = MaxPtIn02;
0836     treeMemPtr_->SumPtIn02 = SumPtIn02;
0837     treeMemPtr_->NoTrackIn0020 = NoTrackIn0020;
0838     treeMemPtr_->MaxPtIn03 = MaxPtIn03;
0839     treeMemPtr_->SumPtIn03 = SumPtIn03;
0840     treeMemPtr_->NoTrackIn0025 = NoTrackIn0025;
0841     treeMemPtr_->MaxPtIn04 = MaxPtIn04;
0842     treeMemPtr_->SumPtIn04 = SumPtIn04;
0843     treeMemPtr_->NoTrackIn0030 = NoTrackIn0030;
0844     treeMemPtr_->MaxPtIn05 = MaxPtIn05;
0845     treeMemPtr_->SumPtIn05 = SumPtIn05;
0846     treeMemPtr_->NoTrackIn0035 = NoTrackIn0035;
0847     treeMemPtr_->NoTrackIn0040 = NoTrackIn0040;
0848     treeMemPtr_->SC_algoID = algo_ID;
0849     treeMemPtr_->SC_energy = EnergyEcal;
0850     treeMemPtr_->SC_nBasicClus = nBasicClus;
0851     treeMemPtr_->SC_etaWidth = etaWidth;
0852     treeMemPtr_->SC_phiWidth = phiWidth;
0853     treeMemPtr_->SC_eta = etaSC;
0854     treeMemPtr_->SC_phi = phiSC;
0855     treeMemPtr_->SC_isBarrel = isBarrel;
0856     treeMemPtr_->SC_isEndcap = isEndcap;
0857     treeMemPtr_->dRto1stSC = dRSC_first;
0858     treeMemPtr_->dRto2ndSC = dRSC_second;
0859     treeMemPtr_->HcalEnergyIn01 = EnergyHcalIn01;
0860     treeMemPtr_->HcalEnergyIn02 = EnergyHcalIn02;
0861     treeMemPtr_->HcalEnergyIn03 = EnergyHcalIn03;
0862     treeMemPtr_->HcalEnergyIn04 = EnergyHcalIn04;
0863     treeMemPtr_->HcalEnergyIn05 = EnergyHcalIn05;
0864     treeMemPtr_->isEcalDriven = isEcalDriven;
0865     treeMemPtr_->isTrackerDriven = isTrackerDriven;
0866     treeMemPtr_->RunNumber = iEvent.id().run();
0867     treeMemPtr_->EvtNumber = iEvent.id().event();
0868
0869     tree_->Fill();
0870
0871   }  // loop on tracks
0872 }  // analyze
0873
0874 // ------------ method called once each job just after ending the event loop  ------------
0875 void EopElecTreeWriter::endJob() {
0876   delete treeMemPtr_;
0877   treeMemPtr_ = nullptr;
0878 }
0879
0880 // ------------ method called once each job just before starting event loop  ------------
0881 void EopElecTreeWriter::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0882   bool changed = true;
0883
0884   // Load trigger configuration at each begin of run
0885   // Fill the trigger names
0886   if (hltConfig_.init(iRun, iSetup, "HLT", changed)) {
0887     if (changed) {
0888       triggerNames_ = hltConfig_.triggerNames();
0889     }
0890   }
0891
0892   // Displaying the trigger names
0893   if (debugTriggerSelection_) {
0894     unsigned int i = 0;
0895     for (const auto& it : triggerNames_) {
0896       edm::LogInfo("EopElecTreeWriter") << "HLTpath: " << (i++) << " = " << it;
0897     }
0898   }
0899 }
0900
0901 //*************************************************************
0902 void EopElecTreeWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0903 //*************************************************************
0904 {
0905   edm::ParameterSetDescription desc;
0906   desc.setComment("Generate tree for Tracker Alignment E/p validation");
0907   desc.add<edm::InputTag>("src", edm::InputTag("electronGsfTracks"));
0908   desc.add<std::string>("triggerPath", "HLT_Ele");
0909   desc.add<std::string>("hltFilter", "hltDiEle27L1DoubleEGWPTightHcalIsoFilter");
0910   desc.add<bool>("debugTriggerSelection", false);
0911   descriptions.addWithDefaultLabel(desc);
0912 }
0913
0914 //define this as a plug-in
0915 DEFINE_FWK_MODULE(EopElecTreeWriter);