File indexing completed on 2024-09-07 04:34:37
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
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
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
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
0116
0117 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0118 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0119 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0120
0121
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
0133 static constexpr float k_etaBarrel = 1.55;
0134 static constexpr float k_etaEndcap = 1.44;
0135
0136
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
0148
0149
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
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
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 }
0219
0220
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
0244 tree_ = fs_->make<TTree>("EopTree", "EopTree");
0245 treeMemPtr_ = new EopElecVariables;
0246 tree_->Branch("EopElecVariables", &treeMemPtr_);
0247
0248
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
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
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;
0352 Int_t nBasicClus;
0353 Bool_t isEcalDriven;
0354 Bool_t isTrackerDriven;
0355 Bool_t isBarrel;
0356 Bool_t isEndcap;
0357 Bool_t TrigTag;
0358
0359
0360
0361 const MagneticField* magField_ = &iSetup.getData(magFieldToken_);
0362 const TrackerGeometry* trackerGeom_ = &iSetup.getData(tkGeomToken_);
0363 MultiTrajectoryStateTransform mtsTransform(trackerGeom_, magField_);
0364
0365
0366
0367
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
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
0384 const HBHERecHitCollection* HcalHits = &iEvent.get(theHBHERecHitCollectionToken_);
0385
0386
0387 const EcalRecHitCollection* rhc = &iEvent.get(theEcalRecHitCollectionToken_);
0388 if (rhc == nullptr)
0389 edm::LogError("EopElecTreeWriter") << "ERROR: unable to find the EcalRecHit collection !!!";
0390
0391
0392 const reco::SuperClusterCollection* BarrelSupClusCollection = &iEvent.get(theBarrelSupClusCollectionToken_);
0393 const reco::SuperClusterCollection* EndcapSupClusCollection = &iEvent.get(theEndCapSupClusCollectionToken_);
0394
0395
0396 SuperClusterShapeAlgo SCShape(rhc, subGeo);
0397
0398
0399 TrigTag = false;
0400 const edm::TriggerResults* trigRes = &iEvent.get(theTriggerResultsToken_);
0401
0402
0403 const trigger::TriggerEvent* triggerEvent = &iEvent.get(theTriggerEventToken_);
0404
0405
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
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
0442
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
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
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
0485
0486
0487 const reco::GsfTrackCollection& tracks = iEvent.get(theGsfTrackCollectionToken_);
0488
0489
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
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);
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
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
0538
0539 for (const auto& track : filterTracks) {
0540
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
0564 h_Momentum->Fill(track->p());
0565 if (track->pt() < 10.)
0566 continue;
0567 h_cut_Ptmin->Fill(0.5);
0568
0569
0570 TrajectoryStateOnSurface inTSOS = mtsTransform.innerStateOnSurface(*track);
0571 TrajectoryStateOnSurface outTSOS = mtsTransform.outerStateOnSurface(*track);
0572 if (inTSOS.isValid() && outTSOS.isValid()) {
0573 GlobalVector inMom;
0574
0575 multiTrajectoryStateMode::momentumFromModePPhiEta(inTSOS, inMom);
0576 pin = inMom.mag();
0577 etaIn = inMom.eta();
0578 phiIn = inMom.phi();
0579 GlobalVector outMom;
0580
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
0590 EnergyHcalIn01 = 0;
0591 EnergyHcalIn02 = 0;
0592 EnergyHcalIn03 = 0;
0593 EnergyHcalIn04 = 0;
0594 EnergyHcalIn05 = 0;
0595
0596
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
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
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
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
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
0705
0706 double dRSC;
0707 double etaAtEcal = 0;
0708 etaAtEcal = eopUtils::ecalEta(etaIn, vert.z(), (vert.position()).rho());
0709
0710
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 }
0718 else if (dRSC < dRSC_second) {
0719 dRSC_second = dRSC;
0720 }
0721
0722 if (dRSC < 0.09) {
0723
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
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 }
0748 else if (dRSC < dRSC_second) {
0749 dRSC_second = dRSC;
0750 }
0751
0752 if (dRSC < 0.09) {
0753
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
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
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 }
0872 }
0873
0874
0875 void EopElecTreeWriter::endJob() {
0876 delete treeMemPtr_;
0877 treeMemPtr_ = nullptr;
0878 }
0879
0880
0881 void EopElecTreeWriter::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0882 bool changed = true;
0883
0884
0885
0886 if (hltConfig_.init(iRun, iSetup, "HLT", changed)) {
0887 if (changed) {
0888 triggerNames_ = hltConfig_.triggerNames();
0889 }
0890 }
0891
0892
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
0915 DEFINE_FWK_MODULE(EopElecTreeWriter);