File indexing completed on 2024-04-06 12:09:06
0001
0002
0003
0004
0005
0006 #include "DQM/TrackingMonitor/interface/dEdxHitAnalyzer.h"
0007
0008 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0009 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0011 #include "MagneticField/Engine/interface/MagneticField.h"
0012
0013 #include "CommonTools/TriggerUtils/interface/GenericTriggerEventFlag.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "DQMServices/Core/interface/DQMStore.h"
0018
0019 #include <string>
0020 #include "TMath.h"
0021
0022 dEdxHitAnalyzer::dEdxHitAnalyzer(const edm::ParameterSet& iConfig)
0023 : fullconf_(iConfig),
0024 conf_(fullconf_.getParameter<edm::ParameterSet>("dEdxParameters")),
0025 doAllPlots_(conf_.getParameter<bool>("doAllPlots")),
0026 doDeDxPlots_(conf_.getParameter<bool>("doDeDxPlots")),
0027 genTriggerEventFlag_(new GenericTriggerEventFlag(
0028 conf_.getParameter<edm::ParameterSet>("genericTriggerEventPSet"), consumesCollector(), *this)) {
0029 trackInputTag_ = edm::InputTag(conf_.getParameter<std::string>("TracksForDeDx"));
0030 trackToken_ = consumes<reco::TrackCollection>(trackInputTag_);
0031
0032 dEdxInputList_ = conf_.getParameter<std::vector<std::string> >("deDxHitProducers");
0033 for (auto const& tag : dEdxInputList_) {
0034 dEdxTokenList_.push_back(consumes<reco::DeDxHitInfoAss>(edm::InputTag(tag)));
0035 }
0036
0037
0038 MEFolderName = conf_.getParameter<std::string>("FolderName");
0039
0040 dEdxNHitBin = conf_.getParameter<int>("dEdxNHitBin");
0041 dEdxNHitMin = conf_.getParameter<double>("dEdxNHitMin");
0042 dEdxNHitMax = conf_.getParameter<double>("dEdxNHitMax");
0043
0044 dEdxStripBin = conf_.getParameter<int>("dEdxStripBin");
0045 dEdxStripMin = conf_.getParameter<double>("dEdxStripMin");
0046 dEdxStripMax = conf_.getParameter<double>("dEdxStripMax");
0047
0048 dEdxPixelBin = conf_.getParameter<int>("dEdxPixelBin");
0049 dEdxPixelMin = conf_.getParameter<double>("dEdxPixelMin");
0050 dEdxPixelMax = conf_.getParameter<double>("dEdxPixelMax");
0051
0052 dEdxHarm2Bin = conf_.getParameter<int>("dEdxHarm2Bin");
0053 dEdxHarm2Min = conf_.getParameter<double>("dEdxHarm2Min");
0054 dEdxHarm2Max = conf_.getParameter<double>("dEdxHarm2Max");
0055 }
0056
0057 dEdxHitAnalyzer::~dEdxHitAnalyzer() {
0058 if (genTriggerEventFlag_)
0059 delete genTriggerEventFlag_;
0060 }
0061
0062
0063
0064 void dEdxHitAnalyzer::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0065
0066 if (genTriggerEventFlag_->on())
0067 genTriggerEventFlag_->initRun(iRun, iSetup);
0068 }
0069
0070 void dEdxHitAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0071 ibooker.setCurrentFolder(MEFolderName);
0072
0073
0074
0075
0076 if (doDeDxPlots_ || doAllPlots_) {
0077 for (unsigned int i = 0; i < dEdxInputList_.size(); i++) {
0078 ibooker.setCurrentFolder(MEFolderName + "/" + dEdxInputList_[i]);
0079 dEdxMEsVector.push_back(dEdxMEs());
0080
0081 histname = "Strip_dEdxPerCluster_";
0082 dEdxMEsVector[i].ME_StripHitDeDx = ibooker.book1D(histname, histname, dEdxStripBin, dEdxStripMin, dEdxStripMax);
0083 dEdxMEsVector[i].ME_StripHitDeDx->setAxisTitle("dEdx of on-track strip cluster (ADC)");
0084 dEdxMEsVector[i].ME_StripHitDeDx->setAxisTitle("Number of Strip clusters", 2);
0085
0086 histname = "Pixel_dEdxPerCluster_";
0087 dEdxMEsVector[i].ME_PixelHitDeDx = ibooker.book1D(histname, histname, dEdxPixelBin, dEdxPixelMin, dEdxPixelMax);
0088 dEdxMEsVector[i].ME_PixelHitDeDx->setAxisTitle("dEdx of on-track pixel cluster (ADC)");
0089 dEdxMEsVector[i].ME_PixelHitDeDx->setAxisTitle("Number of Pixel clusters", 2);
0090
0091 histname = "NumberOfdEdxHitsPerTrack_";
0092 dEdxMEsVector[i].ME_NHitDeDx = ibooker.book1D(histname, histname, dEdxNHitBin, dEdxNHitMin, dEdxNHitMax);
0093 dEdxMEsVector[i].ME_NHitDeDx->setAxisTitle("Number of dEdxHits per Track");
0094 dEdxMEsVector[i].ME_NHitDeDx->setAxisTitle("Number of Tracks", 2);
0095
0096 histname = "Harm2_dEdxPerTrack_";
0097 dEdxMEsVector[i].ME_Harm2DeDx = ibooker.book1D(histname, histname, dEdxHarm2Bin, dEdxHarm2Min, dEdxHarm2Max);
0098 dEdxMEsVector[i].ME_Harm2DeDx->setAxisTitle("Harmonic2 dEdx estimator for each Track");
0099 dEdxMEsVector[i].ME_Harm2DeDx->setAxisTitle("Number of Tracks", 2);
0100 }
0101 }
0102 }
0103
0104 double dEdxHitAnalyzer::harmonic2(const reco::DeDxHitInfo* dedxHits) {
0105 if (!dedxHits)
0106 return -1;
0107 std::vector<double> vect;
0108 for (unsigned int h = 0; h < dedxHits->size(); h++) {
0109 DetId detid(dedxHits->detId(h));
0110 double Norm = (detid.subdetId() < 3) ? 3.61e-06 : 3.61e-06 * 265;
0111 double ChargeOverPathlength = Norm * dedxHits->charge(h) / dedxHits->pathlength(h);
0112 vect.push_back(ChargeOverPathlength);
0113 }
0114
0115 int size = vect.size();
0116 if (size <= 0)
0117 return -1;
0118 double result = 0;
0119 double expo = -2;
0120 for (int i = 0; i < size; i++) {
0121 result += pow(vect[i], expo);
0122 }
0123 return pow(result / size, 1. / expo);
0124 }
0125
0126
0127
0128 void dEdxHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0129
0130 if (genTriggerEventFlag_->on() && !genTriggerEventFlag_->accept(iEvent, iSetup))
0131 return;
0132
0133 if (doDeDxPlots_ || doAllPlots_) {
0134 edm::Handle<reco::TrackCollection> trackCollectionHandle = iEvent.getHandle(trackToken_);
0135 if (!trackCollectionHandle.isValid())
0136 return;
0137
0138 for (unsigned int i = 0; i < dEdxInputList_.size(); i++) {
0139 edm::Handle<reco::DeDxHitInfoAss> dEdxObjectHandle = iEvent.getHandle(dEdxTokenList_[i]);
0140 if (!dEdxObjectHandle.isValid())
0141 continue;
0142
0143 for (unsigned int t = 0; t < trackCollectionHandle->size(); t++) {
0144 reco::TrackRef track = reco::TrackRef(trackCollectionHandle, t);
0145
0146 if (track->quality(reco::TrackBase::highPurity)) {
0147 const reco::DeDxHitInfo* dedxHits = nullptr;
0148 if (!track.isNull()) {
0149 reco::DeDxHitInfoRef dedxHitsRef = (*dEdxObjectHandle)[track];
0150 if (!dedxHitsRef.isNull())
0151 dedxHits = &(*dedxHitsRef);
0152 }
0153 if (!dedxHits)
0154 continue;
0155
0156 for (unsigned int h = 0; h < dedxHits->size(); h++) {
0157 DetId detid(dedxHits->detId(h));
0158 if (detid.subdetId() >= 3)
0159 dEdxMEsVector[i].ME_StripHitDeDx->Fill(dedxHits->charge(h));
0160 if (detid.subdetId() < 3)
0161 dEdxMEsVector[i].ME_PixelHitDeDx->Fill(dedxHits->charge(h));
0162 }
0163 dEdxMEsVector[i].ME_NHitDeDx->Fill(dedxHits->size());
0164 dEdxMEsVector[i].ME_Harm2DeDx->Fill(harmonic2(dedxHits));
0165 }
0166 }
0167 }
0168 }
0169 }
0170
0171
0172 void dEdxHitAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0173
0174
0175 edm::ParameterSetDescription desc;
0176 desc.setUnknown();
0177 descriptions.addDefault(desc);
0178 }
0179
0180 DEFINE_FWK_MODULE(dEdxHitAnalyzer);