Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // system include files
0002 #include <memory>
0003 #include <string>
0004 #include <vector>
0005 
0006 // Root objects
0007 #include "TH1D.h"
0008 
0009 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0010 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0011 
0012 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0013 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0014 
0015 //Tracks
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018 #include "DataFormats/TrackReco/interface/HitPattern.h"
0019 #include "DataFormats/TrackReco/interface/TrackBase.h"
0020 // Vertices
0021 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0024 //Generator information
0025 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0026 
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/ServiceRegistry/interface/Service.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0035 
0036 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0037 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0038 #include "Calibration/IsolatedParticles/interface/eCone.h"
0039 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0040 
0041 #include "MagneticField/Engine/interface/MagneticField.h"
0042 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0043 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0044 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0045 
0046 //#define EDM_ML_DEBUG
0047 
0048 class HcalIsoTrackAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0049 public:
0050   explicit HcalIsoTrackAnalysis(edm::ParameterSet const&);
0051   ~HcalIsoTrackAnalysis() override {}
0052 
0053   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0054 
0055 private:
0056   void analyze(edm::Event const&, edm::EventSetup const&) override;
0057   void beginJob() override;
0058   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0059   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0060 
0061   spr::trackSelectionParameters selectionParameter_;
0062   const std::string theTrackQuality_;
0063   const std::vector<double> maxDxyPV_, maxDzPV_, maxChi2_, maxDpOverP_;
0064   const std::vector<int> minOuterHit_, minLayerCrossed_;
0065   const std::vector<int> maxInMiss_, maxOutMiss_;
0066   const double a_coneR_, a_mipR_;
0067   const double pTrackLow_, pTrackHigh_;
0068   const int useRaw_, dataType_, etaMin_, etaMax_;
0069   const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0070   const double hitEthrEE2_, hitEthrEE3_;
0071   const double hitEthrEELo_, hitEthrEEHi_;
0072   const std::string labelGenTrack_, labelRecVtx_, labelEB_;
0073   const std::string labelEE_, labelHBHE_, labelBS_;
0074   const bool usePFThresh_;
0075   double a_charIsoR_;
0076 
0077   const edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0078   const edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0079   const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0080   const edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0081   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0082   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0083   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0084 
0085   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0086   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0087   const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0088 
0089   const EcalPFRecHitThresholds* eThresholds_;
0090 
0091   std::vector<TH1D*> h_eta_, h_eta0_, h_eta1_, h_rat0_, h_rat1_;
0092   TH1D *h_Dxy_, *h_Dz_, *h_Chi2_, *h_DpOverP_;
0093   TH1D *h_Layer_, *h_OutHit_, *h_InMiss_, *h_OutMiss_;
0094 };
0095 
0096 HcalIsoTrackAnalysis::HcalIsoTrackAnalysis(const edm::ParameterSet& iConfig)
0097     : theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0098       maxDxyPV_(iConfig.getParameter<std::vector<double>>("maxDxyPV")),
0099       maxDzPV_(iConfig.getParameter<std::vector<double>>("maxDzPV")),
0100       maxChi2_(iConfig.getParameter<std::vector<double>>("maxChi2")),
0101       maxDpOverP_(iConfig.getParameter<std::vector<double>>("maxDpOverP")),
0102       minOuterHit_(iConfig.getParameter<std::vector<int>>("minOuterHit")),
0103       minLayerCrossed_(iConfig.getParameter<std::vector<int>>("minLayerCrossed")),
0104       maxInMiss_(iConfig.getParameter<std::vector<int>>("maxInMiss")),
0105       maxOutMiss_(iConfig.getParameter<std::vector<int>>("maxOutMiss")),
0106       a_coneR_(iConfig.getParameter<double>("coneRadius")),
0107       a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0108       pTrackLow_(iConfig.getParameter<double>("momentumLow")),
0109       pTrackHigh_(iConfig.getParameter<double>("momentumHigh")),
0110       useRaw_(iConfig.getUntrackedParameter<int>("useRaw", 0)),
0111       dataType_(iConfig.getUntrackedParameter<int>("dataType", 0)),
0112       etaMin_(iConfig.getUntrackedParameter<int>("etaMin", -1)),
0113       etaMax_(iConfig.getUntrackedParameter<int>("etaMax", 10)),
0114       hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0115       hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0116       hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0117       hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0118       hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0119       hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0120       hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0121       labelGenTrack_(iConfig.getParameter<std::string>("labelTrack")),
0122       labelRecVtx_(iConfig.getParameter<std::string>("labelVertex")),
0123       labelEB_(iConfig.getParameter<std::string>("labelEBRecHit")),
0124       labelEE_(iConfig.getParameter<std::string>("labelEERecHit")),
0125       labelHBHE_(iConfig.getParameter<std::string>("labelHBHERecHit")),
0126       labelBS_(iConfig.getParameter<std::string>("labelBeamSpot")),
0127       usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")),
0128       tok_bs_(consumes<reco::BeamSpot>(labelBS_)),
0129       tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0130       tok_genTrack_(consumes<reco::TrackCollection>(labelGenTrack_)),
0131       tok_recVtx_(consumes<reco::VertexCollection>(labelRecVtx_)),
0132       tok_EB_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEB_))),
0133       tok_EE_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEE_))),
0134       tok_hbhe_(consumes<HBHERecHitCollection>(labelHBHE_)),
0135       tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0136       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0137       tok_ecalPFRecHitThresholds_(esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>()) {
0138   usesResource(TFileService::kSharedResource);
0139 
0140   //now do whatever initialization is needed
0141   const double isolationRadius(28.9);
0142   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0143   selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0144   selectionParameter_.minQuality = trackQuality_;
0145   a_charIsoR_ = a_coneR_ + isolationRadius;
0146   // Different isolation cuts are described in DN-2016/029
0147   // Tight cut uses 2 GeV; Loose cut uses 10 GeV
0148   // Eta dependent cut uses (maxRestrictionP_ * exp(|ieta|*log(2.5)/18))
0149   // with the factor for exponential slopeRestrictionP_ = log(2.5)/18
0150   // maxRestrictionP_ = 8 GeV as came from a study
0151 
0152   // tokens for access
0153   edm::LogVerbatim("HcalIsoTrack") << "Labels used " << labelBS_ << " " << labelRecVtx_ << " " << labelGenTrack_ << " "
0154                                    << edm::InputTag("ecalRecHit", labelEB_) << " "
0155                                    << edm::InputTag("ecalRecHit", labelEE_) << " " << labelHBHE_;
0156 
0157   edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0158                                    << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
0159                                    << theTrackQuality_ << "\t a_coneR " << a_coneR_ << "\t a_charIsoR " << a_charIsoR_
0160                                    << "\t a_mipR " << a_mipR_ << "\n\t momentumLow_ " << pTrackLow_
0161                                    << "\t momentumHigh_ " << pTrackHigh_ << "\t useRaw_ " << useRaw_
0162                                    << "\t dataType_      " << dataType_ << "\t etaLimit " << etaMin_ << ":" << etaMax_
0163                                    << "\nThreshold flag used " << usePFThresh_ << " value for EB " << hitEthrEB_
0164                                    << " EE " << hitEthrEE0_ << ":" << hitEthrEE1_ << ":" << hitEthrEE2_ << ":"
0165                                    << hitEthrEE3_ << ":" << hitEthrEELo_ << ":" << hitEthrEEHi_;
0166 }
0167 
0168 void HcalIsoTrackAnalysis::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0169 #ifdef EDM_ML_DEBUG
0170   edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event() << " type "
0171                                    << dataType_ << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0172                                    << iEvent.bunchCrossing();
0173 #endif
0174   //Get magnetic field
0175   const MagneticField* bField = &iSetup.getData(tok_bFieldH_);
0176 
0177   // get calogeometry
0178   const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0179 
0180   // get ECAL thresholds
0181   eThresholds_ = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0182 
0183   bool okC(true);
0184   //Get track collection
0185   edm::Handle<reco::TrackCollection> trkCollection = iEvent.getHandle(tok_genTrack_);
0186   if (!trkCollection.isValid()) {
0187     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0188     okC = false;
0189   }
0190 
0191   //event weight for FLAT sample
0192   const edm::Handle<GenEventInfoProduct> genEventInfo = iEvent.getHandle(tok_ew_);
0193   double wt = ((genEventInfo.isValid()) ? genEventInfo->weight() : 1.0);
0194 
0195   //Define the best vertex and the beamspot
0196   const edm::Handle<reco::VertexCollection> recVtxs = iEvent.getHandle(tok_recVtx_);
0197   const edm::Handle<reco::BeamSpot> beamSpotH = iEvent.getHandle(tok_bs_);
0198   math::XYZPoint leadPV(0, 0, 0);
0199   bool goodPV(false);
0200   if (recVtxs.isValid() && !(recVtxs->empty())) {
0201     for (unsigned int k = 0; k < recVtxs->size(); ++k) {
0202       if (!((*recVtxs)[k].isFake()) && ((*recVtxs)[k].ndof() > 4)) {
0203         leadPV = math::XYZPoint((*recVtxs)[k].x(), (*recVtxs)[k].y(), (*recVtxs)[k].z());
0204         goodPV = true;
0205         break;
0206       }
0207     }
0208   }
0209   if (!goodPV && beamSpotH.isValid()) {
0210     leadPV = beamSpotH->position();
0211   }
0212 #ifdef EDM_ML_DEBUG
0213   edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex (" << goodPV << ") " << leadPV;
0214   if (beamSpotH.isValid()) {
0215     edm::LogVerbatim("HcalIsoTrack") << " Beam Spot " << beamSpotH->position();
0216   }
0217 #endif
0218   // RecHits
0219   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0220   if (!barrelRecHitsHandle.isValid()) {
0221     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0222     okC = false;
0223   }
0224   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0225   if (!endcapRecHitsHandle.isValid()) {
0226     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0227     okC = false;
0228   }
0229   edm::Handle<HBHERecHitCollection> hbhe = iEvent.getHandle(tok_hbhe_);
0230   if (!hbhe.isValid()) {
0231     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0232     okC = false;
0233   }
0234 
0235   if (okC) {
0236     //Propagate tracks to calorimeter surface)
0237     std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0238     spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0239     std::vector<spr::propagatedTrackID> trkCaloDets;
0240     spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, false);
0241 
0242     //Loop over all tracks
0243     unsigned int nTracks(0);
0244     for (const auto& trkDetItr : trkCaloDirections) {
0245       const reco::Track* pTrack = &(*(trkDetItr.trkItr));
0246       double p = pTrack->p();
0247       if (p >= pTrackLow_ && p <= pTrackHigh_ && (trkDetItr.okHCAL)) {
0248         int ieta = (static_cast<HcalDetId>(trkDetItr.detIdHCAL)).ieta();
0249 
0250         ////////////////////////////////-Energy in ECAL-//////////////////////////
0251         std::vector<DetId> eIds;
0252         std::vector<double> eHit;
0253         double eMipDR = spr::eCone_ecal(geo,
0254                                         barrelRecHitsHandle,
0255                                         endcapRecHitsHandle,
0256                                         trkDetItr.pointHCAL,
0257                                         trkDetItr.pointECAL,
0258                                         a_mipR_,
0259                                         trkDetItr.directionECAL,
0260                                         eIds,
0261                                         eHit);
0262         double eEcal(0);
0263         for (unsigned int k = 0; k < eIds.size(); ++k) {
0264           double eThr(hitEthrEB_);
0265           if (usePFThresh_) {
0266             eThr = static_cast<double>((*eThresholds_)[eIds[k]]);
0267           } else {
0268             const GlobalPoint& pos = geo->getPosition(eIds[k]);
0269             double eta = std::abs(pos.eta());
0270             if (eIds[k].subdetId() != EcalBarrel) {
0271               eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
0272               if (eThr < hitEthrEELo_)
0273                 eThr = hitEthrEELo_;
0274               else if (eThr > hitEthrEEHi_)
0275                 eThr = hitEthrEEHi_;
0276             }
0277           }
0278           if (eHit[k] > eThr)
0279             eEcal += eHit[k];
0280         }
0281 #ifdef EDM_ML_DEBUG
0282         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR << ":" << eEcal;
0283 #endif
0284 
0285         ////////////////////////////////-Energy in HCAL-//////////////////////////
0286         int nRecHits(-999), nNearTRKs(0);
0287         std::vector<DetId> ids;
0288         std::vector<double> edet0;
0289         double eHcal = spr::eCone_hcal(geo,
0290                                        hbhe,
0291                                        trkDetItr.pointHCAL,
0292                                        trkDetItr.pointECAL,
0293                                        a_coneR_,
0294                                        trkDetItr.directionHCAL,
0295                                        nRecHits,
0296                                        ids,
0297                                        edet0,
0298                                        useRaw_);
0299         double ratio0 = eHcal / (p - eEcal);
0300         double ratio1 = eHcal / (p - eMipDR);
0301         double hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
0302         static const double tightCut(2.0), looseCut(2.0);
0303         bool tight = (hmaxNearP < tightCut);
0304         bool loose = (hmaxNearP < looseCut);
0305 #ifdef EDM_ML_DEBUG
0306         edm::LogVerbatim("HcalIsoTrack") << "eHcal and responses: " << eHcal << ":" << ratio0 << ":" << ratio1
0307                                          << " Isolation " << hmaxNearP << ":" << loose << ":" << tight;
0308 #endif
0309         //Different criteria for selection of good tracks
0310         if (std::abs(ieta) > etaMin_ && std::abs(ieta) < etaMax_) {
0311           unsigned id(0);
0312           h_eta_[id]->Fill(ieta, wt);
0313           h_rat0_[id]->Fill(ratio0, wt);
0314           h_rat1_[id]->Fill(ratio1, wt);
0315           if (loose)
0316             h_eta0_[id]->Fill(ieta, wt);
0317           if (tight)
0318             h_eta1_[id]->Fill(ieta, wt);
0319           for (unsigned int k1 = 0; k1 < maxDxyPV_.size(); ++k1) {
0320             for (unsigned int k2 = 0; k2 < maxDzPV_.size(); ++k2) {
0321               for (unsigned int k3 = 0; k3 < maxChi2_.size(); ++k3) {
0322                 for (unsigned int k4 = 0; k4 < maxDpOverP_.size(); ++k4) {
0323                   for (unsigned int k5 = 0; k5 < minOuterHit_.size(); ++k5) {
0324                     for (unsigned int k6 = 0; k6 < minLayerCrossed_.size(); ++k6) {
0325                       for (unsigned int k7 = 0; k7 < maxInMiss_.size(); ++k7) {
0326                         for (unsigned int k8 = 0; k8 < maxOutMiss_.size(); ++k8) {
0327                           ++id;
0328                           selectionParameter_.maxDxyPV = maxDxyPV_[k1];
0329                           selectionParameter_.maxDzPV = maxDzPV_[k2];
0330                           selectionParameter_.maxChi2 = maxChi2_[k3];
0331                           selectionParameter_.maxDpOverP = maxDpOverP_[k4];
0332                           selectionParameter_.minOuterHit = minOuterHit_[k5];
0333                           selectionParameter_.minLayerCrossed = minLayerCrossed_[k6];
0334                           selectionParameter_.maxInMiss = maxInMiss_[k7];
0335                           selectionParameter_.maxOutMiss = maxOutMiss_[k8];
0336                           if (spr::goodTrack(pTrack, leadPV, selectionParameter_, false)) {
0337                             h_eta_[id]->Fill(ieta, wt);
0338                             h_rat0_[id]->Fill(ratio0, wt);
0339                             h_rat1_[id]->Fill(ratio1, wt);
0340                             if (loose)
0341                               h_eta0_[id]->Fill(ieta, wt);
0342                             if (tight)
0343                               h_eta1_[id]->Fill(ieta, wt);
0344                             const reco::HitPattern& hitp = pTrack->hitPattern();
0345                             if ((k2 + k3 + k4 + k5 + k6 + k7 + k8 == 0) && (k1 + 1 == maxDxyPV_.size()))
0346                               h_Dxy_->Fill(pTrack->dxy(leadPV), wt);
0347                             if ((k1 + k3 + k4 + k5 + k6 + k7 + k8 == 0) && (k2 + 1 == maxDzPV_.size()))
0348                               h_Dz_->Fill(pTrack->dz(leadPV), wt);
0349                             if ((k1 + k2 + k4 + k5 + k6 + k7 + k8 == 0) && (k3 + 1 == maxChi2_.size()))
0350                               h_Chi2_->Fill(pTrack->normalizedChi2(), wt);
0351                             if ((k1 + k2 + k3 + k5 + k6 + k7 + k8 == 0) && (k4 + 1 == maxDpOverP_.size()))
0352                               h_DpOverP_->Fill(std::abs(pTrack->qoverpError() / pTrack->qoverp()), wt);
0353                             if ((k1 + k2 + k3 + k4 + k6 + k7 + k8 == 0) && (k5 + 1 == minOuterHit_.size()))
0354                               h_OutHit_->Fill(
0355                                   (hitp.stripTOBLayersWithMeasurement() + hitp.stripTECLayersWithMeasurement()), wt);
0356                             if ((k1 + k2 + k3 + k4 + k5 + k7 + k8 == 0) && (k6 + 1 == minLayerCrossed_.size()))
0357                               h_Layer_->Fill(hitp.trackerLayersWithMeasurement(), wt);
0358                             if ((k1 + k2 + k3 + k4 + k5 + k6 + k8 == 0) && (k7 + 1 == maxInMiss_.size()))
0359                               h_InMiss_->Fill(
0360                                   hitp.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS), wt);
0361                             if ((k1 + k2 + k3 + k4 + k5 + k6 + k7 == 0) && (k8 + 1 == maxOutMiss_.size()))
0362                               h_OutMiss_->Fill(
0363                                   hitp.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS), wt);
0364                           }
0365                         }
0366                       }
0367                     }
0368                   }
0369                 }
0370               }
0371             }
0372           }
0373         }
0374       }
0375       ++nTracks;
0376     }
0377   }
0378 }
0379 
0380 void HcalIsoTrackAnalysis::beginJob() {
0381   edm::Service<TFileService> fs;
0382   char name[100], title[200];
0383   h_eta_.emplace_back(fs->make<TH1D>("eta", "Track i#eta (All)", 60, -30, 30));
0384   h_eta0_.emplace_back(fs->make<TH1D>("eta", "Track i#eta (All Loose Isolation)", 60, -30, 30));
0385   h_eta1_.emplace_back(fs->make<TH1D>("eta", "Track i#eta (All Tight Isolation)", 60, -30, 30));
0386   h_rat0_.emplace_back(fs->make<TH1D>("rat0", "Response 0", 100, 0.0, 5.0));
0387   h_rat1_.emplace_back(fs->make<TH1D>("rat1", "Response 1", 100, 0.0, 5.0));
0388   for (unsigned int k1 = 0; k1 < maxDxyPV_.size(); ++k1) {
0389     for (unsigned int k2 = 0; k2 < maxDzPV_.size(); ++k2) {
0390       for (unsigned int k3 = 0; k3 < maxChi2_.size(); ++k3) {
0391         for (unsigned int k4 = 0; k4 < maxDpOverP_.size(); ++k4) {
0392           for (unsigned int k5 = 0; k5 < minOuterHit_.size(); ++k5) {
0393             for (unsigned int k6 = 0; k6 < minLayerCrossed_.size(); ++k6) {
0394               for (unsigned int k7 = 0; k7 < maxInMiss_.size(); ++k7) {
0395                 for (unsigned int k8 = 0; k8 < maxOutMiss_.size(); ++k8) {
0396                   sprintf(name, "eta%d%d%d%d%d%d%d%d", k1, k2, k3, k4, k5, k6, k7, k8);
0397                   sprintf(title,
0398                           "i#eta (d_{xy}=4.2%f, d_{z}=4.2%f, #chi^{2}=5.2%f, (#Delta p)/p=5.2%f, Hit_{out}=%d, "
0399                           "Layer=%d, Miss_{in}=%d, Miss_{out}=%d)",
0400                           maxDxyPV_[k1],
0401                           maxDzPV_[k2],
0402                           maxChi2_[k3],
0403                           maxDpOverP_[k4],
0404                           minOuterHit_[k5],
0405                           minLayerCrossed_[k6],
0406                           maxInMiss_[k7],
0407                           maxOutMiss_[k8]);
0408                   h_eta_.emplace_back(fs->make<TH1D>(name, title, 60, -30, 30));
0409                   sprintf(name, "eta0%d%d%d%d%d%d%d%d", k1, k2, k3, k4, k5, k6, k7, k8);
0410                   sprintf(title,
0411                           "i#eta (d_{xy}=4.2%f, d_{z}=4.2%f, #chi^{2}=5.2%f, (#Delta p)/p=5.2%f, Hit_{out}=%d, "
0412                           "Layer=%d, Miss_{in}=%d, Miss_{out}=%d, loose isolation)",
0413                           maxDxyPV_[k1],
0414                           maxDzPV_[k2],
0415                           maxChi2_[k3],
0416                           maxDpOverP_[k4],
0417                           minOuterHit_[k5],
0418                           minLayerCrossed_[k6],
0419                           maxInMiss_[k7],
0420                           maxOutMiss_[k8]);
0421                   h_eta0_.emplace_back(fs->make<TH1D>(name, title, 60, -30, 30));
0422                   sprintf(name, "eta1%d%d%d%d%d%d%d%d", k1, k2, k3, k4, k5, k6, k7, k8);
0423                   sprintf(title,
0424                           "i#eta (d_{xy}=4.2%f, d_{z}=4.2%f, #chi^{2}=5.2%f, (#Delta p)/p=5.2%f, Hit_{out}=%d, "
0425                           "Layer=%d, Miss_{in}=%d, Miss_{out}=%d, tight isolation)",
0426                           maxDxyPV_[k1],
0427                           maxDzPV_[k2],
0428                           maxChi2_[k3],
0429                           maxDpOverP_[k4],
0430                           minOuterHit_[k5],
0431                           minLayerCrossed_[k6],
0432                           maxInMiss_[k7],
0433                           maxOutMiss_[k8]);
0434                   h_eta1_.emplace_back(fs->make<TH1D>(name, title, 60, -30, 30));
0435                   sprintf(name, "rat0%d%d%d%d%d%d%d%d", k1, k2, k3, k4, k5, k6, k7, k8);
0436                   sprintf(title,
0437                           "Response 0 (d_{xy}=4.2%f, d_{z}=4.2%f, #chi^{2}=5.2%f, (#Delta p)/p=5.2%f, Hit_{out}=%d, "
0438                           "Layer=%d, Miss_{in}=%d, Miss_{out}=%d)",
0439                           maxDxyPV_[k1],
0440                           maxDzPV_[k2],
0441                           maxChi2_[k3],
0442                           maxDpOverP_[k4],
0443                           minOuterHit_[k5],
0444                           minLayerCrossed_[k6],
0445                           maxInMiss_[k7],
0446                           maxOutMiss_[k8]);
0447                   h_rat0_.emplace_back(fs->make<TH1D>(name, title, 100, 0.0, 5.0));
0448                   sprintf(name, "rat1%d%d%d%d%d%d%d%d", k1, k2, k3, k4, k5, k6, k7, k8);
0449                   sprintf(title,
0450                           "Response 1 (d_{xy}=4.2%f, d_{z}=4.2%f, #chi^{2}=5.2%f, (#Delta p)/p=5.2%f, Hit_{out}=%d, "
0451                           "Layer=%d, Miss_{in}=%d, Miss_{out}=%d)",
0452                           maxDxyPV_[k1],
0453                           maxDzPV_[k2],
0454                           maxChi2_[k3],
0455                           maxDpOverP_[k4],
0456                           minOuterHit_[k5],
0457                           minLayerCrossed_[k6],
0458                           maxInMiss_[k7],
0459                           maxOutMiss_[k8]);
0460                   h_rat1_.emplace_back(fs->make<TH1D>(name, title, 100, 0.0, 5.0));
0461                 }
0462               }
0463             }
0464           }
0465         }
0466       }
0467     }
0468   }
0469   h_Dxy_ = fs->make<TH1D>("Dxy", "d_{xy}", 100, 0.0, 1.0);
0470   h_Dz_ = fs->make<TH1D>("Dz", "d_{z}", 100, 0.0, 1.0);
0471   h_Chi2_ = fs->make<TH1D>("Chi2", "#chi^{2}", 100, 0.0, 20.0);
0472   h_DpOverP_ = fs->make<TH1D>("DpOverP", "#frac{#Delta p}{p}", 100, 0.0, 1.0);
0473   h_Layer_ = fs->make<TH1D>("Layer", "Layers Crossed", 50, 0.0, 50.0);
0474   h_OutHit_ = fs->make<TH1D>("OutHit", "Outer Layers Hit", 20, 0.0, 20.0);
0475   h_InMiss_ = fs->make<TH1D>("InMiss", "Missed Inner Hits", 20, 0.0, 20.0);
0476   h_OutMiss_ = fs->make<TH1D>("OutMiss", "Missed Outer Hits", 20, 0.0, 20.0);
0477 }
0478 
0479 void HcalIsoTrackAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0480   edm::ParameterSetDescription desc;
0481   // following 10 parameters are parameters to select good tracks
0482   desc.add<std::string>("trackQuality", "highPurity");
0483   desc.add<double>("minTrackPt", 1.0);
0484   std::vector<double> maxdxy = {0.02, 0.01, 0.05, 0.10};
0485   std::vector<double> maxdz = {0.02, 0.01, 0.04, 0.50};
0486   std::vector<double> maxchi2 = {5.0, 2.0, 10.0, 20.0};
0487   std::vector<double> maxdpoverp = {0.1, 0.02, 0.05, 0.4};
0488   std::vector<int> minouterhit = {4, 2, 1, 0};
0489   std::vector<int> minlayercrossed = {8, 4, 2, 0};
0490   std::vector<int> maxinmiss = {0, 1, 2, 4};
0491   std::vector<int> maxoutmiss = {0, 1, 2, 4};
0492   desc.add<std::vector<double>>("maxDxyPV", maxdxy);
0493   desc.add<std::vector<double>>("maxDzPV", maxdz);
0494   desc.add<std::vector<double>>("maxChi2", maxchi2);
0495   desc.add<std::vector<double>>("maxDpOverP", maxdpoverp);
0496   desc.add<std::vector<int>>("minOuterHit", minouterhit);
0497   desc.add<std::vector<int>>("minLayerCrossed", minlayercrossed);
0498   desc.add<std::vector<int>>("maxInMiss", maxinmiss);
0499   desc.add<std::vector<int>>("maxOutMiss", maxoutmiss);
0500   // Signal zone in HCAL and ECAL
0501   desc.add<double>("coneRadius", 34.98);
0502   desc.add<double>("coneRadiusMIP", 14.0);
0503   // energy thershold for ECAL (from Egamma group)
0504   desc.add<double>("EBHitEnergyThreshold", 0.08);
0505   desc.add<double>("EEHitEnergyThreshold0", 0.30);
0506   desc.add<double>("EEHitEnergyThreshold1", 0.00);
0507   desc.add<double>("EEHitEnergyThreshold2", 0.00);
0508   desc.add<double>("EEHitEnergyThreshold3", 0.00);
0509   desc.add<double>("EEHitEnergyThresholdLow", 0.30);
0510   desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0511   // prescale factors
0512   desc.add<double>("momentumLow", 40.0);
0513   desc.add<double>("momentumHigh", 60.0);
0514   // various labels for collections used in the code
0515   desc.add<std::string>("labelTrack", "generalTracks");
0516   desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0517   desc.add<std::string>("labelEBRecHit", "EcalRecHitsEB");
0518   desc.add<std::string>("labelEERecHit", "EcalRecHitsEE");
0519   desc.add<std::string>("labelHBHERecHit", "hbhereco");
0520   desc.add<std::string>("labelBeamSpot", "offlineBeamSpot");
0521   //  Various flags used for selecting tracks, choice of energy Method2/0
0522   //  Data type 0/1 for single jet trigger or others
0523   desc.addUntracked<int>("useRaw", 0);
0524   desc.addUntracked<int>("dataType", 0);
0525   desc.addUntracked<int>("etaMin", -1);
0526   desc.addUntracked<int>("etaMax", 10);
0527   desc.add<bool>("usePFThreshold", true);
0528   descriptions.add("hcalIsoTrackAnalysis", desc);
0529 }
0530 
0531 //define this as a plug-in
0532 DEFINE_FWK_MODULE(HcalIsoTrackAnalysis);
[ 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!