Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:15 

0001 // System include files
0002 #include <cmath>
0003 #include <memory>
0004 
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008 
0009 // root objects
0010 #include "TROOT.h"
0011 #include "TSystem.h"
0012 #include "TFile.h"
0013 #include "TH1F.h"
0014 #include "TH2F.h"
0015 #include "TProfile.h"
0016 #include "TDirectory.h"
0017 #include "TTree.h"
0018 
0019 #include <Math/GenVector/VectorUtil.h>
0020 
0021 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0022 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0023 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0024 #include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
0025 #include "Calibration/IsolatedParticles/interface/CaloSimInfo.h"
0026 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0027 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0028 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0029 #include "Calibration/IsolatedParticles/interface/eCone.h"
0030 #include "Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h"
0031 
0032 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0033 
0034 #include "DataFormats/Common/interface/Ref.h"
0035 #include "DataFormats/Common/interface/Handle.h"
0036 #include "DataFormats/Math/interface/Point3D.h"
0037 #include "DataFormats/Candidate/interface/Candidate.h"
0038 
0039 // muons and tracks
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0042 #include "DataFormats/TrackReco/interface/HitPattern.h"
0043 #include "DataFormats/TrackReco/interface/TrackBase.h"
0044 // Vertices
0045 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0046 #include "DataFormats/VertexReco/interface/Vertex.h"
0047 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0048 // Calorimeters
0049 #include "DataFormats/DetId/interface/DetId.h"
0050 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0051 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0052 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0053 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0054 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0055 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0056 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0057 
0058 // Jets in the event
0059 #include "DataFormats/JetReco/interface/CaloJet.h"
0060 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0061 #include "DataFormats/JetReco/interface/JetExtendedAssociation.h"
0062 
0063 #include "FWCore/Framework/interface/Frameworkfwd.h"
0064 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0065 #include "FWCore/Framework/interface/Event.h"
0066 #include "FWCore/Framework/interface/EventSetup.h"
0067 #include "FWCore/Framework/interface/MakerMacros.h"
0068 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0069 
0070 // TFile Service
0071 #include "FWCore/ServiceRegistry/interface/Service.h"
0072 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0073 
0074 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0075 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0076 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0077 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0078 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0079 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0080 
0081 #include "MagneticField/Engine/interface/MagneticField.h"
0082 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0083 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0084 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0085 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0086 
0087 // SimHit
0088 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0089 #include "SimDataFormats/Track/interface/SimTrack.h"
0090 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0091 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0092 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0093 
0094 // track associator
0095 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0096 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0097 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0098 
0099 class IsolatedTracksHcalScale : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0100 public:
0101   explicit IsolatedTracksHcalScale(const edm::ParameterSet &);
0102   ~IsolatedTracksHcalScale() override {}
0103 
0104   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0105 
0106 private:
0107   void beginJob() override;
0108   void analyze(const edm::Event &, const edm::EventSetup &) override;
0109   void endJob() override;
0110 
0111   void clearTreeVectors();
0112 
0113 private:
0114   bool doMC_;
0115   int myverbose_;
0116   std::string theTrackQuality_, minQuality_;
0117   spr::trackSelectionParameters selectionParameters_;
0118   double a_mipR_, a_coneR_, a_charIsoR_, a_neutIsoR_;
0119   double tMinE_, tMaxE_;
0120 
0121   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0122 
0123   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0124   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0125   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0126   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0127   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0128   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0129 
0130   edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0131   edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0132   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEB_;
0133   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEE_;
0134   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloHH_;
0135 
0136   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0137   edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0138   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0139   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0140   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0141 
0142   int nEventProc_;
0143 
0144   TTree *tree_;
0145 
0146   int t_nTracks, t_RunNo, t_EvtNo, t_Lumi, t_Bunch;
0147   std::vector<double> *t_trackP, *t_trackPt, *t_trackEta, *t_trackPhi;
0148   std::vector<double> *t_trackHcalEta, *t_trackHcalPhi, *t_eHCALDR;
0149   std::vector<double> *t_hCone, *t_conehmaxNearP, *t_eMipDR, *t_eECALDR;
0150   std::vector<double> *t_e11x11_20Sig, *t_e15x15_20Sig;
0151   std::vector<double> *t_eMipDR_1, *t_eECALDR_1, *t_eMipDR_2, *t_eECALDR_2;
0152   std::vector<double> *t_hConeHB, *t_eHCALDRHB;
0153   std::vector<double> *t_hsimInfoMatched, *t_hsimInfoRest, *t_hsimInfoPhoton;
0154   std::vector<double> *t_hsimInfoNeutHad, *t_hsimInfoCharHad, *t_hsimInfoPdgMatched;
0155   std::vector<double> *t_hsimInfoTotal, *t_hsim;
0156   std::vector<int> *t_hsimInfoNMatched, *t_hsimInfoNTotal, *t_hsimInfoNNeutHad;
0157   std::vector<int> *t_hsimInfoNCharHad, *t_hsimInfoNPhoton, *t_hsimInfoNRest;
0158   std::vector<int> *t_nSimHits;
0159 };
0160 
0161 IsolatedTracksHcalScale::IsolatedTracksHcalScale(const edm::ParameterSet &iConfig)
0162     : doMC_(iConfig.getUntrackedParameter<bool>("DoMC", false)),
0163       myverbose_(iConfig.getUntrackedParameter<int>("Verbosity", 5)),
0164       theTrackQuality_(iConfig.getUntrackedParameter<std::string>("TrackQuality", "highPurity")),
0165       a_mipR_(iConfig.getUntrackedParameter<double>("ConeRadiusMIP", 14.0)),
0166       a_coneR_(iConfig.getUntrackedParameter<double>("ConeRadius", 34.98)),
0167       tMinE_(iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.)),
0168       tMaxE_(iConfig.getUntrackedParameter<double>("TimeMaxCutECAL", 500.)),
0169       trackerHitAssociatorConfig_(consumesCollector()) {
0170   usesResource(TFileService::kSharedResource);
0171 
0172   //now do what ever initialization is needed
0173   reco::TrackBase::TrackQuality trackQuality = reco::TrackBase::qualityByName(theTrackQuality_);
0174   selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
0175   selectionParameters_.minQuality = trackQuality;
0176   selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
0177   selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("MaxDzPV", 5.0);
0178   selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("MaxChi2", 5.0);
0179   selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("MaxDpOverP", 0.1);
0180   selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
0181   selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
0182   selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
0183   selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
0184   a_charIsoR_ = a_coneR_ + 28.9;
0185   a_neutIsoR_ = a_charIsoR_ * 0.726;
0186 
0187   tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0188   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0189   tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0190   tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0191   tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0192   tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0193   tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0194   tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
0195   tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
0196   tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"));
0197   tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
0198 
0199   if (myverbose_ >= 0) {
0200     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0201                                  << " doMC " << doMC_ << "\t myverbose " << myverbose_ << "\t minPt "
0202                                  << selectionParameters_.minPt << "\t theTrackQuality " << theTrackQuality_
0203                                  << "\t minQuality " << selectionParameters_.minQuality << "\t maxDxyPV "
0204                                  << selectionParameters_.maxDxyPV << "\t maxDzPV " << selectionParameters_.maxDzPV
0205                                  << "\t maxChi2 " << selectionParameters_.maxChi2 << "\t maxDpOverP "
0206                                  << selectionParameters_.maxDpOverP << "\t minOuterHit "
0207                                  << selectionParameters_.minOuterHit << "\t minLayerCrossed "
0208                                  << selectionParameters_.minLayerCrossed << "\t maxInMiss "
0209                                  << selectionParameters_.maxInMiss << "\t maxOutMiss "
0210                                  << selectionParameters_.maxOutMiss << "\t a_coneR " << a_coneR_ << "\t a_charIsoR "
0211                                  << a_charIsoR_ << "\t a_neutIsoR " << a_neutIsoR_ << "\t a_mipR " << a_mipR_
0212                                  << "\t time Range (" << tMinE_ << ":" << tMaxE_ << ")";
0213   }
0214 
0215   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0216   tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0217   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0218   tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0219   tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0220 }
0221 
0222 void IsolatedTracksHcalScale::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0223   edm::ParameterSetDescription desc;
0224   desc.addUntracked<bool>("doMC", false);
0225   desc.addUntracked<int>("Verbosity", 0);
0226   desc.addUntracked<std::string>("TrackQuality", "highPurity");
0227   desc.addUntracked<double>("MinTrackPt", 10.0);
0228   desc.addUntracked<double>("MaxDxyPV", 0.02);
0229   desc.addUntracked<double>("MaxDzPV", 0.02);
0230   desc.addUntracked<double>("MaxChi2", 5.0);
0231   desc.addUntracked<double>("MaxDpOverP", 0.1);
0232   desc.addUntracked<int>("MinOuterHit", 4);
0233   desc.addUntracked<int>("MinLayerCrossed", 8);
0234   desc.addUntracked<int>("MaxInMiss", 0);
0235   desc.addUntracked<int>("MaxOutMiss", 0);
0236   desc.addUntracked<double>("ConeRadius", 34.98);
0237   desc.addUntracked<double>("ConeRadiusMIP", 14.0);
0238   desc.addUntracked<double>("TimeMinCutECAL", -500.0);
0239   desc.addUntracked<double>("TimeMaxCutECAL", 500.0);
0240   descriptions.add("isolatedTracksHcalScale", desc);
0241 }
0242 
0243 void IsolatedTracksHcalScale::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0244   const CaloGeometry *geo = &iSetup.getData(tok_geom_);
0245   const MagneticField *bField = &iSetup.getData(tok_magField_);
0246   const EcalChannelStatus *theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0247   const EcalSeverityLevelAlgo *sevlv = &iSetup.getData(tok_sevlv_);
0248   const CaloTopology *caloTopology = &iSetup.getData(tok_caloTopology_);
0249 
0250   clearTreeVectors();
0251 
0252   ++nEventProc_;
0253 
0254   t_RunNo = iEvent.id().run();
0255   t_EvtNo = iEvent.id().event();
0256   t_Lumi = iEvent.luminosityBlock();
0257   t_Bunch = iEvent.bunchCrossing();
0258   if (myverbose_ > 0)
0259     edm::LogVerbatim("IsoTrack") << nEventProc_ << " Run " << t_RunNo << " Event " << t_EvtNo << " Lumi " << t_Lumi
0260                                  << " Bunch " << t_Bunch;
0261 
0262   edm::Handle<reco::TrackCollection> trkCollection;
0263   iEvent.getByToken(tok_genTrack_, trkCollection);
0264 
0265   edm::Handle<reco::VertexCollection> recVtxs;
0266   iEvent.getByToken(tok_recVtx_, recVtxs);
0267 
0268   // Get the beamspot
0269   edm::Handle<reco::BeamSpot> beamSpotH;
0270   iEvent.getByToken(tok_bs_, beamSpotH);
0271 
0272   math::XYZPoint leadPV(0, 0, 0);
0273   if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
0274     leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
0275   } else if (beamSpotH.isValid()) {
0276     leadPV = beamSpotH->position();
0277   }
0278 
0279   if (myverbose_ > 0) {
0280     edm::LogVerbatim("IsoTrack") << "Primary Vertex " << leadPV;
0281     if (beamSpotH.isValid())
0282       edm::LogVerbatim("IsoTrack") << "Beam Spot " << beamSpotH->position();
0283   }
0284 
0285   std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0286   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, (myverbose_ > 2));
0287   std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0288 
0289   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0290   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0291   iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0292   iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0293 
0294   edm::Handle<HBHERecHitCollection> hbhe;
0295   iEvent.getByToken(tok_hbhe_, hbhe);
0296   const HBHERecHitCollection Hithbhe = *(hbhe.product());
0297 
0298   //get Handles to SimTracks and SimHits
0299   edm::Handle<edm::SimTrackContainer> SimTk;
0300   edm::Handle<edm::SimVertexContainer> SimVtx;
0301 
0302   //get Handles to PCaloHitContainers of eb/ee/hbhe
0303   edm::Handle<edm::PCaloHitContainer> pcaloeb;
0304   edm::Handle<edm::PCaloHitContainer> pcaloee;
0305   edm::Handle<edm::PCaloHitContainer> pcalohh;
0306 
0307   //associates tracker rechits/simhits to a track
0308   std::unique_ptr<TrackerHitAssociator> associate;
0309 
0310   if (doMC_) {
0311     iEvent.getByToken(tok_simTk_, SimTk);
0312     iEvent.getByToken(tok_simVtx_, SimVtx);
0313     iEvent.getByToken(tok_caloEB_, pcaloeb);
0314     iEvent.getByToken(tok_caloEE_, pcaloee);
0315     iEvent.getByToken(tok_caloHH_, pcalohh);
0316     associate = std::make_unique<TrackerHitAssociator>(iEvent, trackerHitAssociatorConfig_);
0317   }
0318 
0319   unsigned int nTracks = 0;
0320   for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0321        trkDetItr++, nTracks++) {
0322     const reco::Track *pTrack = &(*(trkDetItr->trkItr));
0323     if (spr::goodTrack(pTrack, leadPV, selectionParameters_, (myverbose_ > 2)) && trkDetItr->okECAL &&
0324         trkDetItr->okHCAL) {
0325       int nRH_eMipDR = 0, nRH_eDR = 0, nNearTRKs = 0, nRecHitsCone = -99;
0326       double distFromHotCell = -99.0, distFromHotCell2 = -99.0;
0327       int ietaHotCell = -99, iphiHotCell = -99;
0328       int ietaHotCell2 = -99, iphiHotCell2 = -99;
0329       GlobalPoint gposHotCell(0., 0., 0.), gposHotCell2(0., 0., 0.);
0330       std::vector<DetId> coneRecHitDetIds, coneRecHitDetIds2;
0331       std::pair<double, bool> e11x11_20SigP, e15x15_20SigP;
0332       double hCone = spr::eCone_hcal(geo,
0333                                      hbhe,
0334                                      trkDetItr->pointHCAL,
0335                                      trkDetItr->pointECAL,
0336                                      a_coneR_,
0337                                      trkDetItr->directionHCAL,
0338                                      nRecHitsCone,
0339                                      coneRecHitDetIds,
0340                                      distFromHotCell,
0341                                      ietaHotCell,
0342                                      iphiHotCell,
0343                                      gposHotCell,
0344                                      -1);
0345       double hConeHB = spr::eCone_hcal(geo,
0346                                        hbhe,
0347                                        trkDetItr->pointHCAL,
0348                                        trkDetItr->pointECAL,
0349                                        a_coneR_,
0350                                        trkDetItr->directionHCAL,
0351                                        nRecHitsCone,
0352                                        coneRecHitDetIds,
0353                                        distFromHotCell,
0354                                        ietaHotCell,
0355                                        iphiHotCell,
0356                                        gposHotCell,
0357                                        (int)(HcalBarrel));
0358       double eHCALDR = spr::eCone_hcal(geo,
0359                                        hbhe,
0360                                        trkDetItr->pointHCAL,
0361                                        trkDetItr->pointECAL,
0362                                        a_charIsoR_,
0363                                        trkDetItr->directionHCAL,
0364                                        nRecHitsCone,
0365                                        coneRecHitDetIds2,
0366                                        distFromHotCell2,
0367                                        ietaHotCell2,
0368                                        iphiHotCell2,
0369                                        gposHotCell2,
0370                                        -1);
0371       double eHCALDRHB = spr::eCone_hcal(geo,
0372                                          hbhe,
0373                                          trkDetItr->pointHCAL,
0374                                          trkDetItr->pointECAL,
0375                                          a_charIsoR_,
0376                                          trkDetItr->directionHCAL,
0377                                          nRecHitsCone,
0378                                          coneRecHitDetIds2,
0379                                          distFromHotCell2,
0380                                          ietaHotCell2,
0381                                          iphiHotCell2,
0382                                          gposHotCell2,
0383                                          (int)(HcalBarrel));
0384 
0385       double conehmaxNearP =
0386           spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, (myverbose_ > 3));
0387 
0388       double eMipDR = spr::eCone_ecal(geo,
0389                                       barrelRecHitsHandle,
0390                                       endcapRecHitsHandle,
0391                                       trkDetItr->pointHCAL,
0392                                       trkDetItr->pointECAL,
0393                                       a_mipR_,
0394                                       trkDetItr->directionECAL,
0395                                       nRH_eMipDR);
0396       double eECALDR = spr::eCone_ecal(geo,
0397                                        barrelRecHitsHandle,
0398                                        endcapRecHitsHandle,
0399                                        trkDetItr->pointHCAL,
0400                                        trkDetItr->pointECAL,
0401                                        a_neutIsoR_,
0402                                        trkDetItr->directionECAL,
0403                                        nRH_eDR);
0404       double eMipDR_1 = spr::eCone_ecal(geo,
0405                                         barrelRecHitsHandle,
0406                                         endcapRecHitsHandle,
0407                                         trkDetItr->pointHCAL,
0408                                         trkDetItr->pointECAL,
0409                                         a_mipR_,
0410                                         trkDetItr->directionECAL,
0411                                         nRH_eMipDR,
0412                                         0.030,
0413                                         0.150);
0414       double eECALDR_1 = spr::eCone_ecal(geo,
0415                                          barrelRecHitsHandle,
0416                                          endcapRecHitsHandle,
0417                                          trkDetItr->pointHCAL,
0418                                          trkDetItr->pointECAL,
0419                                          a_neutIsoR_,
0420                                          trkDetItr->directionECAL,
0421                                          nRH_eDR,
0422                                          0.030,
0423                                          0.150);
0424       double eMipDR_2 = spr::eCone_ecal(geo,
0425                                         barrelRecHitsHandle,
0426                                         endcapRecHitsHandle,
0427                                         trkDetItr->pointHCAL,
0428                                         trkDetItr->pointECAL,
0429                                         a_mipR_,
0430                                         trkDetItr->directionECAL,
0431                                         nRH_eMipDR,
0432                                         0.060,
0433                                         0.300);
0434       double eECALDR_2 = spr::eCone_ecal(geo,
0435                                          barrelRecHitsHandle,
0436                                          endcapRecHitsHandle,
0437                                          trkDetItr->pointHCAL,
0438                                          trkDetItr->pointECAL,
0439                                          a_neutIsoR_,
0440                                          trkDetItr->directionECAL,
0441                                          nRH_eDR,
0442                                          0.060,
0443                                          0.300);
0444 
0445       HcalDetId closestCell = (HcalDetId)(trkDetItr->detIdHCAL);
0446 
0447       e11x11_20SigP = spr::eECALmatrix(trkDetItr->detIdECAL,
0448                                        barrelRecHitsHandle,
0449                                        endcapRecHitsHandle,
0450                                        *theEcalChStatus,
0451                                        geo,
0452                                        caloTopology,
0453                                        sevlv,
0454                                        5,
0455                                        5,
0456                                        0.060,
0457                                        0.300,
0458                                        tMinE_,
0459                                        tMaxE_);
0460       e15x15_20SigP = spr::eECALmatrix(trkDetItr->detIdECAL,
0461                                        barrelRecHitsHandle,
0462                                        endcapRecHitsHandle,
0463                                        *theEcalChStatus,
0464                                        geo,
0465                                        caloTopology,
0466                                        sevlv,
0467                                        7,
0468                                        7,
0469                                        0.060,
0470                                        0.300,
0471                                        tMinE_,
0472                                        tMaxE_);
0473 
0474       // Fill the tree Branches here
0475       t_trackP->push_back(pTrack->p());
0476       t_trackPt->push_back(pTrack->pt());
0477       t_trackEta->push_back(pTrack->momentum().eta());
0478       t_trackPhi->push_back(pTrack->momentum().phi());
0479       t_trackHcalEta->push_back(closestCell.ieta());
0480       t_trackHcalPhi->push_back(closestCell.iphi());
0481       t_hCone->push_back(hCone);
0482       t_conehmaxNearP->push_back(conehmaxNearP);
0483       t_eMipDR->push_back(eMipDR);
0484       t_eECALDR->push_back(eECALDR);
0485       t_eHCALDR->push_back(eHCALDR);
0486       t_e11x11_20Sig->push_back(e11x11_20SigP.first);
0487       t_e15x15_20Sig->push_back(e15x15_20SigP.first);
0488       t_eMipDR_1->push_back(eMipDR_1);
0489       t_eECALDR_1->push_back(eECALDR_1);
0490       t_eMipDR_2->push_back(eMipDR_2);
0491       t_eECALDR_2->push_back(eECALDR_2);
0492       t_hConeHB->push_back(hConeHB);
0493       t_eHCALDRHB->push_back(eHCALDRHB);
0494 
0495       if (myverbose_ > 0) {
0496         edm::LogVerbatim("IsoTrack") << "Track p " << pTrack->p() << " pt " << pTrack->pt() << " eta "
0497                                      << pTrack->momentum().eta() << " phi " << pTrack->momentum().phi()
0498                                      << " ieta/iphi (" << closestCell.ieta() << ", " << closestCell.iphi()
0499                                      << ") Energy in cone " << hCone << " Charge Isolation " << conehmaxNearP
0500                                      << " eMIP (" << eMipDR << ", " << eMipDR_1 << ", " << eMipDR_2 << ")"
0501                                      << " Neutral isolation (ECAL) (" << eECALDR - eMipDR << ", "
0502                                      << eECALDR_1 - eMipDR_1 << ", " << eECALDR_2 - eMipDR_2 << ") (ECAL NxN) "
0503                                      << e15x15_20SigP.first - e11x11_20SigP.first << " (HCAL) " << eHCALDR - hCone;
0504       }
0505 
0506       if (doMC_) {
0507         int nSimHits = -999;
0508         double hsim;
0509         std::map<std::string, double> hsimInfo;
0510         std::vector<int> multiplicity;
0511         hsim = spr::eCone_hcal(
0512             geo, pcalohh, trkDetItr->pointHCAL, trkDetItr->pointECAL, a_coneR_, trkDetItr->directionHCAL, nSimHits);
0513         hsimInfo = spr::eHCALSimInfoCone(iEvent,
0514                                          pcalohh,
0515                                          SimTk,
0516                                          SimVtx,
0517                                          pTrack,
0518                                          *associate,
0519                                          geo,
0520                                          trkDetItr->pointHCAL,
0521                                          trkDetItr->pointECAL,
0522                                          a_coneR_,
0523                                          trkDetItr->directionHCAL,
0524                                          multiplicity);
0525 
0526         t_hsimInfoMatched->push_back(hsimInfo["eMatched"]);
0527         t_hsimInfoRest->push_back(hsimInfo["eRest"]);
0528         t_hsimInfoPhoton->push_back(hsimInfo["eGamma"]);
0529         t_hsimInfoNeutHad->push_back(hsimInfo["eNeutralHad"]);
0530         t_hsimInfoCharHad->push_back(hsimInfo["eChargedHad"]);
0531         t_hsimInfoPdgMatched->push_back(hsimInfo["pdgMatched"]);
0532         t_hsimInfoTotal->push_back(hsimInfo["eTotal"]);
0533 
0534         t_hsimInfoNMatched->push_back(multiplicity.at(0));
0535         t_hsimInfoNTotal->push_back(multiplicity.at(1));
0536         t_hsimInfoNNeutHad->push_back(multiplicity.at(2));
0537         t_hsimInfoNCharHad->push_back(multiplicity.at(3));
0538         t_hsimInfoNPhoton->push_back(multiplicity.at(4));
0539         t_hsimInfoNRest->push_back(multiplicity.at(5));
0540 
0541         t_hsim->push_back(hsim);
0542         t_nSimHits->push_back(nSimHits);
0543 
0544         if (myverbose_ > 0) {
0545           edm::LogVerbatim("IsoTrack") << "Matched (E) " << hsimInfo["eMatched"] << " (N) " << multiplicity.at(0)
0546                                        << " Rest (E) " << hsimInfo["eRest"] << " (N) " << multiplicity.at(1)
0547                                        << " Gamma (E) " << hsimInfo["eGamma"] << " (N) " << multiplicity.at(2)
0548                                        << " Neutral Had (E) " << hsimInfo["eNeutralHad"] << " (N) "
0549                                        << multiplicity.at(3) << " Charged Had (E) " << hsimInfo["eChargedHad"]
0550                                        << " (N) " << multiplicity.at(4) << " Total (E) " << hsimInfo["eTotal"]
0551                                        << " (N) " << multiplicity.at(5) << " PDG " << hsimInfo["pdgMatched"]
0552                                        << " Total E " << hsim << " NHit " << nSimHits;
0553         }
0554       }
0555     }
0556   }
0557 
0558   tree_->Fill();
0559 }
0560 
0561 void IsolatedTracksHcalScale::beginJob() {
0562   nEventProc_ = 0;
0563   edm::Service<TFileService> fs;
0564 
0565   //////Book Tree
0566   tree_ = fs->make<TTree>("tree", "tree");
0567   tree_->SetAutoSave(10000);
0568 
0569   tree_->Branch("t_RunNo", &t_RunNo, "t_RunNo/I");
0570   tree_->Branch("t_Lumi", &t_Lumi, "t_Lumi/I");
0571   tree_->Branch("t_Bunch", &t_Bunch, "t_Bunch/I");
0572 
0573   t_trackP = new std::vector<double>();
0574   t_trackPt = new std::vector<double>();
0575   t_trackEta = new std::vector<double>();
0576   t_trackPhi = new std::vector<double>();
0577   t_trackHcalEta = new std::vector<double>();
0578   t_trackHcalPhi = new std::vector<double>();
0579   t_hCone = new std::vector<double>();
0580   t_conehmaxNearP = new std::vector<double>();
0581   t_eMipDR = new std::vector<double>();
0582   t_eECALDR = new std::vector<double>();
0583   t_eHCALDR = new std::vector<double>();
0584   t_e11x11_20Sig = new std::vector<double>();
0585   t_e15x15_20Sig = new std::vector<double>();
0586   t_eMipDR_1 = new std::vector<double>();
0587   t_eECALDR_1 = new std::vector<double>();
0588   t_eMipDR_2 = new std::vector<double>();
0589   t_eECALDR_2 = new std::vector<double>();
0590   t_hConeHB = new std::vector<double>();
0591   t_eHCALDRHB = new std::vector<double>();
0592 
0593   tree_->Branch("t_trackP", "std::vector<double>", &t_trackP);
0594   tree_->Branch("t_trackPt", "std::vector<double>", &t_trackPt);
0595   tree_->Branch("t_trackEta", "std::vector<double>", &t_trackEta);
0596   tree_->Branch("t_trackPhi", "std::vector<double>", &t_trackPhi);
0597   tree_->Branch("t_trackHcalEta", "std::vector<double>", &t_trackHcalEta);
0598   tree_->Branch("t_trackHcalPhi", "std::vector<double>", &t_trackHcalPhi);
0599   tree_->Branch("t_hCone", "std::vector<double>", &t_hCone);
0600   tree_->Branch("t_conehmaxNearP", "std::vector<double>", &t_conehmaxNearP);
0601   tree_->Branch("t_eMipDR", "std::vector<double>", &t_eMipDR);
0602   tree_->Branch("t_eECALDR", "std::vector<double>", &t_eECALDR);
0603   tree_->Branch("t_eHCALDR", "std::vector<double>", &t_eHCALDR);
0604   tree_->Branch("t_e11x11_20Sig", "std::vector<double>", &t_e11x11_20Sig);
0605   tree_->Branch("t_e15x15_20Sig", "std::vector<double>", &t_e15x15_20Sig);
0606   tree_->Branch("t_eMipDR_1", "std::vector<double>", &t_eMipDR_1);
0607   tree_->Branch("t_eECALDR_1", "std::vector<double>", &t_eECALDR_1);
0608   tree_->Branch("t_eMipDR_2", "std::vector<double>", &t_eMipDR_2);
0609   tree_->Branch("t_eECALDR_2", "std::vector<double>", &t_eECALDR_2);
0610   tree_->Branch("t_hConeHB", "std::vector<double>", &t_hConeHB);
0611   tree_->Branch("t_eHCALDRHB", "std::vector<double>", &t_eHCALDRHB);
0612 
0613   if (doMC_) {
0614     t_hsimInfoMatched = new std::vector<double>();
0615     t_hsimInfoRest = new std::vector<double>();
0616     t_hsimInfoPhoton = new std::vector<double>();
0617     t_hsimInfoNeutHad = new std::vector<double>();
0618     t_hsimInfoCharHad = new std::vector<double>();
0619     t_hsimInfoPdgMatched = new std::vector<double>();
0620     t_hsimInfoTotal = new std::vector<double>();
0621     t_hsimInfoNMatched = new std::vector<int>();
0622     t_hsimInfoNTotal = new std::vector<int>();
0623     t_hsimInfoNNeutHad = new std::vector<int>();
0624     t_hsimInfoNCharHad = new std::vector<int>();
0625     t_hsimInfoNPhoton = new std::vector<int>();
0626     t_hsimInfoNRest = new std::vector<int>();
0627     t_hsim = new std::vector<double>();
0628     t_nSimHits = new std::vector<int>();
0629 
0630     tree_->Branch("t_hsimInfoMatched", "std::vector<double>", &t_hsimInfoMatched);
0631     tree_->Branch("t_hsimInfoRest", "std::vector<double>", &t_hsimInfoRest);
0632     tree_->Branch("t_hsimInfoPhoton", "std::vector<double>", &t_hsimInfoPhoton);
0633     tree_->Branch("t_hsimInfoNeutHad", "std::vector<double>", &t_hsimInfoNeutHad);
0634     tree_->Branch("t_hsimInfoCharHad", "std::vector<double>", &t_hsimInfoCharHad);
0635     tree_->Branch("t_hsimInfoPdgMatched", "std::vector<double>", &t_hsimInfoPdgMatched);
0636     tree_->Branch("t_hsimInfoTotal", "std::vector<double>", &t_hsimInfoTotal);
0637     tree_->Branch("t_hsimInfoNMatched", "std::vector<int>", &t_hsimInfoNMatched);
0638     tree_->Branch("t_hsimInfoNTotal", "std::vector<int>", &t_hsimInfoNTotal);
0639     tree_->Branch("t_hsimInfoNNeutHad", "std::vector<int>", &t_hsimInfoNNeutHad);
0640     tree_->Branch("t_hsimInfoNCharHad", "std::vector<int>", &t_hsimInfoNCharHad);
0641     tree_->Branch("t_hsimInfoNPhoton", "std::vector<int>", &t_hsimInfoNPhoton);
0642     tree_->Branch("t_hsimInfoNRest", "std::vector<int>", &t_hsimInfoNRest);
0643     tree_->Branch("t_hsim", "std::vector<double>", &t_hsim);
0644     tree_->Branch("t_nSimHits", "std::vector<int>", &t_nSimHits);
0645   }
0646 }
0647 
0648 void IsolatedTracksHcalScale::endJob() { edm::LogVerbatim("IsoTrack") << "Number of Events Processed " << nEventProc_; }
0649 
0650 void IsolatedTracksHcalScale::clearTreeVectors() {
0651   t_trackP->clear();
0652   t_trackPt->clear();
0653   t_trackEta->clear();
0654   t_trackPhi->clear();
0655   t_trackHcalEta->clear();
0656   t_trackHcalPhi->clear();
0657   t_hCone->clear();
0658   t_conehmaxNearP->clear();
0659   t_eMipDR->clear();
0660   t_eECALDR->clear();
0661   t_eHCALDR->clear();
0662   t_e11x11_20Sig->clear();
0663   t_e15x15_20Sig->clear();
0664   t_eMipDR_1->clear();
0665   t_eECALDR_1->clear();
0666   t_eMipDR_2->clear();
0667   t_eECALDR_2->clear();
0668   t_hConeHB->clear();
0669   t_eHCALDRHB->clear();
0670 
0671   if (doMC_) {
0672     t_hsimInfoMatched->clear();
0673     t_hsimInfoRest->clear();
0674     t_hsimInfoPhoton->clear();
0675     t_hsimInfoNeutHad->clear();
0676     t_hsimInfoCharHad->clear();
0677     t_hsimInfoPdgMatched->clear();
0678     t_hsimInfoTotal->clear();
0679     t_hsimInfoNMatched->clear();
0680     t_hsimInfoNTotal->clear();
0681     t_hsimInfoNNeutHad->clear();
0682     t_hsimInfoNCharHad->clear();
0683     t_hsimInfoNPhoton->clear();
0684     t_hsimInfoNRest->clear();
0685     t_hsim->clear();
0686     t_nSimHits->clear();
0687   }
0688 }
0689 
0690 //define this as a plug-in
0691 DEFINE_FWK_MODULE(IsolatedTracksHcalScale);
[ 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!